aperiodic 2018\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Journal logoFOUNDATIONS
ADVANCES
ISSN: 2053-2733

Hyperuniformity and anti-hyperuniformity in one-dimensional substitution tilings

CROSSMARK_Color_square_no_text.svg

aSchool of Mechanical Engineering and The Sackler Center for Computational Molecular and Materials Science, Tel Aviv University, Tel Aviv, Israel, bPhysics Department, Duke University, Durham, NC 27708, USA, cDepartment of Physics and Princeton Center for Theoretical Science, Princeton University, Princeton, NJ 08544, USA, dCenter for Cosmology and Particle Physics, Department of Physics, New York University, New York, NY, 10003, USA, and eDepartment of Chemistry, Department of Physics, Princeton Institute for the Science and Technology of Materials and Program in Applied and Computational Mathematics, Princeton University, Princeton, New Jersey 08540, USA
*Correspondence e-mail: socolar@phy.duke.edu

Edited by P. A. Thiel, Iowa State University, USA (Received 21 August 2018; accepted 2 November 2018)

This work considers the scaling properties characterizing the hyperuniformity (or anti-hyperuniformity) of long-wavelength fluctuations in a broad class of one-dimensional substitution tilings. A simple argument is presented which predicts the exponent α governing the scaling of Fourier intensities at small wavenumbers, tilings with α > 0 being hyperuniform, and numerical computations confirm that the predictions are accurate for quasiperiodic tilings, tilings with singular continuous spectra and limit-periodic tilings. Quasiperiodic or singular continuous cases can be constructed with α arbitrarily close to any given value between −1 and 3. Limit-periodic tilings can be constructed with α between −1 and 1 or with Fourier intensities that approach zero faster than any power law.

1. Introduction

Recent work has shown that spatial structures with density fluctuations weaker at long wavelengths than those of a typical random point set may have desirable physical properties, and such structures are said to be hyperuniform (Torquato & Stillinger, 2003[Torquato, S. & Stillinger, F. H. (2003). Phys. Rev. E, 68, 041113.]). Crystals and quasicrystals are hyperuniform, as are a variety of disordered systems, including certain equilibrium structures, products of nonequilibrium self-assembly protocols and fabricated metamaterials. [For examples, see Man et al. (2013[Man, W., Florescu, M., Williamson, E. P., He, Y., Hashemizad, S. R., Leung, B. Y. C., Liner, D. R., Torquato, S., Chaikin, P. M. & Steinhardt, P. J. (2013). Proc. Natl Acad. Sci. USA, 110, 15886-15891.]), Haberko & Scheffold (2013[Haberko, J. & Scheffold, F. (2013). Opt. Express, 21, 1057-1065.]), Dreyfus et al. (2015[Dreyfus, R., Xu, Y., Still, T., Hough, L. A., Yodh, A. G. & Torquato, S. (2015). Phys. Rev. E, 91, 012302.]), Torquato et al. (2015[Torquato, S., Zhang, G. & Stillinger, F. H. (2015). Phys. Rev. X, 5, 021020.]), Hexner & Levine (2015[Hexner, D. & Levine, D. (2015). Phys. Rev. Lett. 114, 110602.]), Castro-Lopez et al. (2017[Castro-Lopez, M., Gaio, M., Sellers, S., Gkantzounis, G., Florescu, M. & Sapienza, R. (2017). APL Photonics, 2, 061302.]), Torquato (2018[Torquato, S. (2018). Phys. Rep. 745, 1-95.]).] One approach to generating point sets with nontrivial spatial fluctuations is to use substitution tilings as templates. Our aim in this article is to characterize the degree of hyperuniformity in such systems and thereby provide design principles for creating hyperuniform (or anti-hyperuniform) point sets with desired scaling properties.

Substitution tilings are self-similar, space-filling tilings generated by repeated application of a rule that replaces each of a finite set of tile types with scaled copies of some or all of the tiles in the set (Frank, 2008[Frank, N. P. (2008). Expo. Math. 26, 295-326.]). We are interested in the properties of point sets formed by decorating each tile of the same type in the same way. We consider here only one-dimensional (1D) tilings. Although generalization to higher dimensions would be of great interest, the 1D case already reveals important conceptual features.

Substitution rules are known to produce a variety of structures with qualitatively different types of structure factors S(k). Some rules generate periodic or quasiperiodic tilings, in which case S(k) consists of Bragg peaks on a reciprocal-space lattice supported at sums and differences of a (small) set of basis wavevectors, which in the quasiperiodic case form a dense set. Others produce limit-periodic structures consisting of Bragg peaks located on a different type of dense set consisting of wavenumbers of the form [\pm k_0 n/p^m], where n, m and p are positive integers (Godrèche, 1989[Godrèche, C. (1989). J. Phys. A Math. Gen. 22, L1163-L1166.]; Baake & Grimm, 2011[Baake, M. & Grimm, U. (2011). Philos. Mag. 91, 2661-2670.], 2013[Baake, M. & Grimm, U. (2013). Aperiodic Order, Vol. 1 of Encyclopedia of Mathematics and its Applications. Cambridge University Press.]). Still others produce structures for which S(k) is singular at a dense set of points but does not consist of Bragg peaks (Bombieri & Taylor, 1986[Bombieri, E. & Taylor, J. E. (1986). J. Phys. Coll. 47, C3-C19.]; Godrèche & Luck, 1992[Godrèche, C. & Luck, J. M. (1992). Phys. Rev. B, 45, 176-185.]; Baake et al., 2017[Baake, M., Frank, N. P., Grimm, U. & Robinson Jr, E. A. (2017). https://arxiv.org/abs/1706.03976.]). (A singular continuous spectrum has support on some set of zero Lebesgue measure, but has no finite weight at any single point.) We note in particular that a detailed analysis of the spectrum of substitution tilings with non-PV properties (defined below) reveals multifractal scaling laws (Godrèche & Luck, 1992[Godrèche, C. & Luck, J. M. (1992). Phys. Rev. B, 45, 176-185.]). Finally, there are cases for which S(k) is absolutely continuous (Baake & Grimm, 2012[Baake, M. & Grimm, U. (2012). Chem. Soc. Rev. 41, 6821-6843.]) or the nature of the spectrum has not been clearly described.

In this article, we present a simple ansatz that predicts the scaling properties relevant for assessing the hyperuniformity (or anti-hyperuniformity) of 1D substitution tilings. We illustrate the validity of the ansatz via numerical computations for a variety of example tilings that fall in different classes with respect to hyperuniformity measures. We also delineate the full range of behaviors that can be obtained using the substitution construction method, including a novel class in which the integrated Fourier intensity Z(k) decays faster than any power as k approaches zero.

Section 2[link] reviews the definition of the scaling exponent α associated with both Z(k) and the variance [\sigma^2(R)] in the number of points covered by a randomly placed interval of length 2R. We then review the classification of tilings based on the value of α. Section 3[link] reviews the substitution method for creating tilings, using the well known Fibonacci tiling as an illustrative example. The substitution matrix [{\bf M}] is defined and straightforward results for tile densities are derived. Section 4[link] presents a heuristic discussion of the link between density fluctuations in the tilings and the behaviors of S(k) and Z(k), which leads to a prediction for α. The prediction is shown to be accurate for example tilings of three qualitatively distinct types (Torquato, 2018[Torquato, S. (2018). Phys. Rep. 745, 1-95.]): strongly hyperuniform (class I), weakly hyperuniform (class III) and anti-hyperuniform. Section 5[link] shows, based on the heuristic theory, that the range of possible values of α produced by 1D substitution rules is [-1,3] and that this interval is densely filled. Section 6[link] considers substitutions that produce limit-periodic tilings. Examples are presented of four distinct classes: logarithmic hyperuniform (class II), weakly hyperuniform (class III), anti-hyperuniform, and an anomalous class in which Z(k) approaches zero faster than any power law. Finally, Section 7[link] provides a summary of the key results, including a table showing which types of tilings can exhibit the various classes of (anti-)hyperuniformity.

2. Classes of hyperuniformity

For systems having a structure factor [S({\bf k})] that is a smooth function of the wavenumber k, [S({\bf k})] tends to zero as k tends to zero (Torquato & Stillinger, 2003[Torquato, S. & Stillinger, F. H. (2003). Phys. Rev. E, 68, 041113.]), typically scaling as a power law:

[S({\bf k}) \sim k^{\alpha}.\eqno(1)]

In 1D, a unified treatment of standard cases with smooth S(k) and quasicrystals with dense but discontinuous S(k) is obtained by defining α in terms of the scaling of the integrated Fourier intensity:

[Z(k) = 2 \textstyle\int\limits_0^k S(q)\,{\rm d}q.\eqno(2)]

The factor of 2 is inserted for consistency with higher-dimensional generalizations where q is treated as a radial coordinate. In both cases, α may be defined by the relation (Oğuz et al., 2017[Oğuz, E. C., Socolar, J. E. S., Steinhardt, P. J. & Torquato, S. (2017). Phys. Rev. B, 95, 054119.])

[Z(k) \sim k^{1+\alpha} \,\,{\rm as} \,\, k \rightarrow 0.\eqno(3)]

Systems with [\alpha\,\gt\, 0] have long-wavelength spatial fluctuations that are suppressed compared with Poisson point sets and are said to be hyperuniform (Torquato & Stillinger, 2003[Torquato, S. & Stillinger, F. H. (2003). Phys. Rev. E, 68, 041113.]). Prototypical strongly hyperuniform systems (with [\alpha\,\gt\, 1]) include crystals and quasicrystals. We refer to systems with [\alpha\,\lt\, 0] as anti-hyperuniform (Torquato, 2018[Torquato, S. (2018). Phys. Rep. 745, 1-95.]). Prototypical examples of anti-hyperuniformity include systems at thermal critical points.

An alternate measure of hyperuniformity is based on the local number variance of particles within a spherical observation window of radius R (an interval of length 2R in the 1D case), denoted by [\sigma^2(R)]. If [\sigma^2(R)] grows more slowly than the window volume (proportional to R in 1D) in the large-R limit, the system is hyperuniform. The scaling behavior of [\sigma^2(R)] is closely related to the behavior of Z(k) for small k (Torquato & Stillinger, 2003[Torquato, S. & Stillinger, F. H. (2003). Phys. Rev. E, 68, 041113.]; Oğuz et al., 2017[Oğuz, E. C., Socolar, J. E. S., Steinhardt, P. J. & Torquato, S. (2017). Phys. Rev. B, 95, 054119.]). For a general point configuration in 1D with a well-defined average number density ρ, [\sigma^2(R)] can be expressed in terms of S(k) and the Fourier transform [{\tilde \mu}(k\semi R)] of a uniform density interval of length 2R:

[\sigma^2(R) = 2R\rho \left[{{1} \over {2\pi}} \int\limits_{-\infty}^{\infty} S({k}) {\tilde \mu}(k\semi R) \,{\rm d}{k}\right]\eqno(4)]

with

[{\tilde \mu}(k\semi R) = 2{{\sin^2(k R)} \over {k^2 R}},\eqno(5)]

where ρ is the density. [See Torquato & Stillinger (2003[Torquato, S. & Stillinger, F. H. (2003). Phys. Rev. E, 68, 041113.]) for the generalization to higher Euclidean space dimensions.] One can express the number variance alternatively in terms of the integrated intensity (Oğuz et al., 2017[Oğuz, E. C., Socolar, J. E. S., Steinhardt, P. J. & Torquato, S. (2017). Phys. Rev. B, 95, 054119.]):

[\sigma^2(R) = - 2R\rho\left[{{1}\over{2\pi}} \int\limits_0^\infty Z(k) {{\partial{\tilde\mu}(k\semi R)}\over{\partial k}} \,{\rm d}k \right].\eqno(6)]

For any 1D system with a smooth or quasicrystalline structure factor, the scaling of [\sigma^2(R)] for large R is determined by α as follows (Torquato & Stillinger, 2003[Torquato, S. & Stillinger, F. H. (2003). Phys. Rev. E, 68, 041113.]; Zachary & Torquato, 2009[Zachary, C. E. & Torquato, S. (2009). J. Stat. Mech. 2009, P12015.]; Torquato, 2018[Torquato, S. (2018). Phys. Rep. 745, 1-95.]):

[\sigma^2(R) \sim \left\{\matrix{ R^0,\hfill & \alpha \,\gt\, 1\,\, {\rm (class\ I)} \hfill\cr \ln R, \hfill& \alpha = 1 \,\,{\rm (class\ II)} \hfill\cr R^{1-\alpha}, & \alpha \,\lt\, 1 \,\,{\rm (class\ III)}}\right..\eqno(7)]

For hyperuniform systems, we have [\alpha\,\gt\, 0], and the distinct behaviors of [\sigma^2(R)] define the three classes, which we refer to as strongly hyperuniform (class I), logarithmic hyperuniform (class II) and weakly hyperuniform (class III). As mentioned above, systems with [\alpha \,\lt\, 0] are classified as anti-hyperuniform.

The bounded number fluctuations of class I occur trivially for 1D periodic point sets (crystals) and are also known to occur for certain quasicrystals, including the canonical Fibonacci tiling described below (Oğuz et al., 2017[Oğuz, E. C., Socolar, J. E. S., Steinhardt, P. J. & Torquato, S. (2017). Phys. Rev. B, 95, 054119.]). Other quasiperiodic point sets (not obtainable by substitution) are known to belong to class II (Kesten, 1966[Kesten, H. (1966). Acta Arith. 12, 193-212.]; Aubry et al., 1987[Aubry, S., Godrèche, C. & Vallet, F. (1987). J. Phys. Fr. 48, 327-334.]; Oğuz et al., 2017[Oğuz, E. C., Socolar, J. E. S., Steinhardt, P. J. & Torquato, S. (2017). Phys. Rev. B, 95, 054119.]).

3. Substitution tilings and the substitution matrix

A classic example of a substitution tiling is the 1D Fibonacci tiling composed of two intervals (tiles) of length L and S. The tiling is generated by the rule

[L\rightarrow LS, \quad S\rightarrow L,\eqno(8)]

which leads to a quasiperiodic sequence of L and S intervals. An important construct for characterizing the properties of the tiling is the substitution matrix:

[{\bf M} = \left(\matrix{0 & 1 \cr 1 & 1 }\right),\eqno(9)]

which acts on the column vector (NS,NL) to give the numbers of S and L tiles resulting from the substitution operation.

If the lengths L and S are chosen such that the ratio L/S remains fixed, which in the present case requires [L/S = (1+5^{1/2})/2\equiv\tau], the substitution operation can be viewed as an affine stretching of the original tiling by a factor of τ followed by the division of each stretched L tile into an LS pair, as illustrated in Fig. 1[link]. Given a finite sequence with NS tiles of length S and NL tiles of length L, the numbers of L's and S's in the system after one iteration of the substitution rule are given by the action of the substitution matrix on the column vector (NS,NL).

[Figure 1]
Figure 1
The Fibonacci substitution rule. The tiling on the upper line is uniformly stretched, then additional points are added to form tiles congruent to the originals.

More generally, substitution rules can be defined for systems with more than two tile types, leading to substitution matrices with dimension [{\cal D}] greater than 2. We present explicit reasoning here only for the [{\cal D}] = 2 case. A substitution rule for two tile types is characterized by a substitution matrix:

[{\bf M} = \left(\matrix{ a & b \cr c & d }\right).\eqno(10)]

The associated rule may be the following:

[S\rightarrow \underbrace{SS\ldots S}_{a}\underbrace{LL\ldots L}_{c},\,\, L\rightarrow \underbrace{SS\ldots S}_{b}\underbrace{LL\ldots L}_{d},\eqno(11)]

but different orderings of the tiles in the substituted strings are possible, and the choice can have dramatic effects. Note, for example, that the rule

[S\rightarrow SL, \quad L\rightarrow SLSL\eqno(12)]

produces the periodic tiling [\ldots SLSLSL \ldots], while the rule

[S\rightarrow SL, \quad L\rightarrow SLLS\eqno(13)]

produces the more complicated sequence discussed below in Section 6[link].

Defining the substitution tiling requires assigning finite lengths to S and L. We let ξ denote the length ratio L/S, and we consider only cases where the substitution rule preserves this ratio [i.e. (bS + dL)/(aS + cL) = L/S] so that the rule can be realized by affine stretching followed by subdivision. This requires

[\xi = {{d-a+[(a-d)^2 + 4bc]^{1/2}} \over {2c}}.\eqno(14)]

For all discussions and plots below, we measure lengths in units of the short tile length, S.

The SL sequence generated by a substitution rule is obtained by repeated application of that operation to some seed, which we will take to be a string containing nS short intervals and nL long ones. We are interested in point sets formed by decorating each L tile with l points and each S tile with s points. The total number of points at the mth iteration is

[{\cal N}_m = (s,\ell)\cdot {\bf M}^m \cdot (n_S,n_L),\eqno(15)]

and the length of the tiling at the same step is

[{\cal X}_m = (1,\xi) \cdot {\bf M}^m \cdot (n_S,n_L).\eqno(16)]

Let [\lambda_1] and [\lambda_2] be the eigenvalues of [{\bf M}], with [\lambda_1] being the largest, and let [{{\bf v}}_1] and [{{\bf v}}_2] be the associated eigenvectors. We have

[\lambda_1 = a + c\xi\semi \,\, \lambda_2 = d - c\xi\semi \eqno(17)]

[{{\bf v}}_1 = (b/c, \xi)\semi \,\,{{\bf v}}_2 = (-\xi, 1).\eqno(18)]

The unit vectors (1,0) and (0,1) may be expressed as follows:

[(1,0) = u(c{{\bf v}}_1 - c\xi {{\bf v}}_2),\eqno(19)]

[(0,1) = u(c\xi {{\bf v}}_1 + b {{\bf v}}_2),\eqno(20)]

where [u = 1/(b+c\xi^2)]. We then have

[\eqalignno{{\bf M}^m \cdot(n_S,n_L) & = {\bf M}^m \cdot\left[n_S(1,0)+n_L(0,1)\right] &\cr & = u \bigg[\lambda_1^m (cn_S+ c\xi n_L){{\bf v}}_1 &\cr &\quad + \lambda_2^m (-c\xi n_S+ bn_L){{\bf v}}_2 \bigg]. &(21)}]

The density of tile vertices after m iterations, [\rho_m = {\cal N}_m/{\cal X}_m], is thus

[\rho_m = \overline{\rho} + \left[{{\xi(s\xi-\ell)} \over {b+c\xi^2}}\right]\left({{c\xi n_S - b n_L} \over {n_S+\xi n_L }}\right) \left({{\lambda_2}\over{\lambda_1}} \right)^m,\eqno(22)]

with [\overline{\rho} = (b s+c\ell\xi)/(b+c\xi^2)], where we have used the fact that [(1,\xi)\cdot{{\bf v}}_2 = 0].

4. Scaling properties of 1D substitution tilings

As long as the coefficient of [(\lambda_2/\lambda_1)^m] in equation (22)[link] does not vanish, the deviations of ρ from [\overline{\rho}] for portions of the tiling that are mapped into each other by substitution are related by

[\delta\rho_{m+1} = {{\lambda_2}\over{\lambda_1}}\delta\rho_m.\eqno(23) ]

If the coefficient does vanish, which requires that ξ be rational, the tiling may be periodic, but the ordering of the intervals in the seed becomes important. We will revisit this point below. For now we assume that the tiling is not periodic.

We make three conjectures regarding nonperiodic substitution tilings, supported, as we shall see, by numerical experiments. The results are closely related to recently derived rigorous results (Baake, Gaehler et al., 2018[Baake, M., Gaehler, F. & Manibo, N. (2018). See Section 3.4, Equation (18). https://arxiv.org/abs/1805.09650.]).

Conjecture 1. We take equation (23)[link] to be the dominant behavior of density fluctuations throughout the system, not just for the special intervals that are directly related by substitution. That is, we assume that there exists a characteristic amplitude of the density fluctuations at a given length scale after averaging over all intervals of that length, and that the [\delta\rho] in equation (23)[link] can be interpreted as that characteristic amplitude.

Conjecture 2. For the present purposes, we define the real-space density [g(x) = \sum_n \delta(x-x_n)] to consist of unit-strength δ-functions placed at every point xn where two tiles meet, and consider the Fourier amplitudes

[A(k) = \lim_{h\rightarrow\infty}\left|{{1} \over {2h}}\int\limits_{-h}^{h} g(x) \exp(i k x)\,{\rm d}x \, \right|.\eqno(24)]

With this definition, the average density over a large domain is equal to the number density of tiles, ρ, in that domain. We assume that A(k) scales the same way as the density fluctuations at the corresponding length scale:

[A(k/\lambda_1) = {{\lambda_2} \over {\lambda_1}}A(k).\eqno(25)]

This implies the form

[A(k) \sim k^{(-\ln |\lambda_2/\lambda_1 |/ \ln |\lambda_1 |)} = k^{1-(\ln |\lambda_2 |/ \ln |\lambda_1 |)}. \eqno(26)]

Squaring to get S(k), we have

[S(k) \sim k^{(2 - 2 \ln |\lambda_2 |/ \ln |\lambda_1 |)}.\eqno(27)]

This conjecture may not hold when interference effects are important, as in the case discussed in Section 6[link] below.

Conjecture 3. While Z(k) is an integral of S(k), the exponent must be calculated carefully when S(k) consists of singular peaks. In the Fibonacci projection cases, the scaling of peak positions and intensities conspires to make Z(k) scale with the same exponent as the envelope of S(k) (Oğuz et al., 2017[Oğuz, E. C., Socolar, J. E. S., Steinhardt, P. J. & Torquato, S. (2017). Phys. Rev. B, 95, 054119.]). We assume that this property carries over to substitution tilings with more than one eigenvalue greater than unity. Though the diffraction pattern is not made up of Bragg peaks (Bombieri & Taylor, 1986[Bombieri, E. & Taylor, J. E. (1986). J. Phys. Coll. 47, C3-C19.]; Godreche & Luck, 1990[Godreche, C. & Luck, J. M. (1990). J. Phys. A Math. Gen. 23, 3769-3797.]), we conjecture that it remains sufficiently singular for the relation to hold. Thus we immediately obtain

[\alpha = 1 - 2\left({{\ln |\lambda_2 |} \over {\ln \lambda_1}}\right).\eqno(28)]

Note that this calculation of the scaling exponent makes no reference to the distinction between substitutions with [|\lambda_2|\,\lt\,1] and those with [|\lambda_2|\,\gt\,1]. In the former case, [\lambda_1] is a Pisot–Vijayaraghavan (PV) number, S(k) consists of Bragg peaks and [\sigma^2(R)] remains bounded for all R. In the latter case, the form of S(k) is more complex (Bombieri & Taylor, 1986[Bombieri, E. & Taylor, J. E. (1986). J. Phys. Coll. 47, C3-C19.]), and quantities closely related to [\sigma^2(R)], including the `wandering exponent' associated with lifts of the sequence onto a higher-dimensional hypercubic lattice, are known to show nontrivial scaling exponents (Godreche & Luck, 1990[Godreche, C. & Luck, J. M. (1990). J. Phys. A Math. Gen. 23, 3769-3797.]).

From equation (28)[link], we see that the hyperuniformity condition [\alpha \,\gt\, 0] requires [|\lambda_2| \,\lt\, ({\lambda_1})^{1/2}]. Though the result was obtained for substitutions with only [{\cal D}] = 2 tile types, it holds for [{\cal D} \,\gt\, 2] as well, so long as all ratios of tile lengths are preserved by the substitution rules; i.e. the dominant contribution to the long-wavelength fluctuations still scales like [|\lambda_2| / \lambda_1]. This distinction between hyperuniform and anti-hyperuniform substitution tilings thus divides the non-PV numbers into two classes which, to our knowledge, have not previously been identified as significantly different. We note, for example, that the analysis presented in Baake, Grimm et al. (2018[Baake, M., Grimm, U. & Manibo, N. (2018). https://arxiv.org/abs/1709.09083.]), which treats substitution matrices of the form (0,n,1,1) and shows that they have singular continuous spectra (having no Bragg component or absolutely continuous component) for [n\,\gt\, 2], does not detect any qualitative difference between the cases n = 3 and n = 5. The former case is hyperuniform, with [\lambda = (1/2)[1\pm(13)^{1/2}]] and [\alpha \simeq 0.37], while the latter is anti-hyperuniform, with [\lambda = (1/2)[1\pm(21)^{1/2}]] and [\alpha \simeq -0.14].

For the Fibonacci case, we have [\lambda_1 = \tau] and [\lambda_2 = -1/\tau], yielding [\alpha = 3], which agrees with the explicit calculation in Oğuz et al. (2017[Oğuz, E. C., Socolar, J. E. S., Steinhardt, P. J. & Torquato, S. (2017). Phys. Rev. B, 95, 054119.]). Considering (a,b,c,d) of the form (0,n,n,n) for arbitrary n, we find cases that allow explicit checks of our predictions for α for both hyperuniform and anti-hyperuniform systems. We have [\lambda_1 = n\tau] and [\lambda_2 = -n/\tau], yielding

[\alpha=1-2\left({{\ln n-\ln \tau}\over{\ln n+\ln\tau}}\right).\eqno(29)]

For [n\ge 2], the presence of more than one eigenvalue with magnitude greater than unity gives rise to more complex spectral features, possibly including a singular continuous component. For [2 \le n \le 4], our calculation predicts [0 \,\lt\, \alpha \,\lt\, 1] and hence [\sigma^2(R) \sim R^{1-\alpha}]. We numerically verify the latter result for n = 2 using a set of 954 369 points generated by 12 iterations of the substitution tiling, where the decoration consists of placing one point at the rightmost edge of each tile (with s = l = 1). Fig. 2[link] shows the computed number variance. For each point, a window of length 2R is moved continuously along the sequence and averages are computed by weighting the number of points in the window by the interval length over which that number does not change. A regression analysis yields [\sigma^2 (R) \sim R^{0.36}], in close agreement with the predicted exponent from equation (27)[link]: [1-\alpha =] [2 (\ln 2 - \ln \tau) / (\ln 2 + \ln \tau)] ≃ 0.36094.

[Figure 2]
Figure 2
Log–log plot of the number variance (black dots) for a non-PV substitution tiling corresponding to (a,b,c,d) = (0,2,2,2) decorated with points of equal weight at each tile boundary. The variance was computed numerically for the tiling created by 12 iterations of the substitution on the initial seed SL. The red dashed line has the predicted slope [1-\alpha \simeq 0.36].

For [n \ge 5], the calculated value of α is negative, approaching -1 as n approaches infinity. The point set is therefore anti-hyperuniform; it contains density fluctuations at long wavelengths that are stronger than those of a Poisson point set. For n = 5, we have [\alpha = -0.0793\ldots]. Fig. 3[link] shows a log–log plot of the computed number variance along with the line corresponding to [\sigma^2 (R) \sim R^{1-\alpha}]. Again, the agreement between the numerical result and the predicted value is quite good. Intuition derived from theories based on nonsingular forms of S(k) suggests that a negative value of α should be associated with a divergence in S(k) for small k, though it remains true that Z(k) converges to zero for [\alpha \,\gt\, -1]. For singular spectra, the envelope of S(k) scales like Z(k), and we do not expect any dramatic change in the behavior of S(k) as α crosses from positive (hyperuniform) to negative (anti-hyperuniform). The theories presented in Baake et al. (2017[Baake, M., Frank, N. P., Grimm, U. & Robinson Jr, E. A. (2017). https://arxiv.org/abs/1706.03976.]) and Godreche & Luck (1990[Godreche, C. & Luck, J. M. (1990). J. Phys. A Math. Gen. 23, 3769-3797.]) may provide a path to the computation of scaling properties of S(k) in these cases. It is worth noting, however, that the various classes of behavior can be realized by substitutions that produce limit-periodic tilings with S(k) consisting entirely of Bragg peaks with no singular-continuous component, as shown in Section 6[link] below.

[Figure 3]
Figure 3
Log–log plot of the number variance (black dots) for an anti-hyperuniform substitution tiling corresponding to (a,b,c,d) = (0, 5, 5, 5) decorated with points of equal weight at each tile boundary. The variance was computed numerically for the tiling created by six iterations of the substitution on the initial seed SL. The red dashed line has the predicted slope [1-\alpha \simeq 1.08].

For rules that yield rational values of the length ratio ξ, the coefficient of [(\lambda_2/\lambda_1)^m] in equation (22)[link] can vanish for appropriate choices of nS and nL, suggesting that there are no fluctuations about the average density that scale with wavelength. This reflects the fact that the sequence of intervals associated with the substitutions can be chosen to generate a periodic pattern. (A simple example is [S\rightarrow L] and [L\rightarrow SLS], which generates the periodic sequence [\ldots SLSLSL\ldots], with [\xi = 2], [\lambda_1 = 2] and [\lambda_2 = -1].) For such cases, S(k) is identically 0 for all k smaller than the reciprocal-lattice basis vector. For other interval sequence choices corresponding to the same [{\bf M}], the tiling can be limit-periodic, and we would expect the scaling to be given by applying the above considerations with generic choices of the ordering, which would yield [\alpha = 1] and therefore a logarithmic scaling of [\sigma^2(R)]. This case is presented in more detail in Section 6[link] below, and the logarithmic scaling is confirmed.

5. Achievable values of α

Beyond establishing that substitution tilings exist for each hyperuniformity class, it is natural to ask whether any desired value of α can be realized by this construction method. Here we show that if [{\bf M}] is full rank, α always lies between -1 and 3.

First, note that the maximum value of [|\lambda_2/\lambda_1|] is 1, by definition, which sets the lower bound on α via equation (28)[link]. The upper bound on α is obtained when [|\lambda_2|] is as small as possible, but there is a limit on how small this can be. The product of the eigenvalues of [{\bf M}] is equal to [\det {\bf M}], so [|\lambda_2|] cannot be smaller than [(|\det{\bf M}| / \lambda_1)^{{{1} / {({\cal D}-1)}}}]. But [|\det{\bf M}|] is an integer, and the smallest nonzero value it can take is 1. (The case [\lambda_2 = 0] is discussed in Section 6[link] below. For [{\cal D} \ge 3], one can have [\det{\bf M} = 0] with nonzero [\lambda_2]. The analysis of such cases is beyond our present scope.) Hence we have

[|\lambda_2| \ge \lambda_1^{-1/({\cal D}-1)}, \eqno(30)]

implying

[\alpha = 1 - 2 {{\ln|\lambda_2|} \over {\ln\lambda_1}} \le {{{\cal D}+1} \over {{\cal D}-1}}.\eqno(31)]

Thus the maximum value of α obtainable by this construction method is 3, which can occur for [{\cal D}] = 2, as in the Fibonacci case.

The family of substitutions considered in Section 4[link] above produces a discrete set of values of α ranging from -1 to 3. By considering two additional families, we can show that the possible values of α densely fill this interval. For

[{\bf M} = \left(\matrix{ a & 0 \cr c & d}\right)\eqno(32)]

with [d\,\gt\, a+1] and [2c\,\lt\, (d-a)], we have [\lambda_1 = a] and [\lambda_2 = d]. Note that [\xi = (d-a)/c] is rational here; we assume that the substitution sequences for the two tiles are chosen so as to avoid periodicity. We have

[\alpha = 1 - 2 {{\ln a} \over {\ln d}}.\eqno(33)]

For fixed a, d can range from a+2 to [\infty]. As d approaches infinity, α approaches 1. For d = a+2, as a approaches infinity, α approaches -1. For sufficiently large d, the values of a between 1 and d-2 yield an arbitrarily dense set of α's between -1 and 1.

Another class of [{\bf M}]'s produces α's between 1 and 3. For

[{\bf M} = \left(\matrix{ 0 & b \cr b & nb }\right)\eqno(34)]

with [n\,\gt\, b], we have

[\xi = {{1}\over{2}}[n + (n^2+4)^{1/2}],\eqno(35)]

[\lambda_{1,2} & = {{b}\over{2}}[n \pm (n^2+4)^{1/2}].\eqno(36)]

We thus obtain

[\alpha = 1 - 2 {{\ln b - \ln 2 + \ln [(n^2+4)^{1/2}-n]} \over {\ln b - \ln 2 + \ln [(n^2+4)^{1/2}+n]}}.\eqno(37)]

For large n, we have

[\alpha \simeq 1 - 2 {{\ln b - 2\ln 2 - \ln n} \over {\ln b + \ln n}},\eqno(38)]

which approaches 3 for [b \ll n] and approaches 1 for b = n. By making n as large as desired, the values of b between 1 and n give α's that fill the interval between 1 and 3 with arbitrarily high density.

6. Limit-periodic tilings

For a limit-periodic tiling, the set of tiles is a union of periodic patterns with ever-increasing lattice constants of the form a pn, where p is an integer and n runs over all positive definite integers (Godrèche, 1989[Godrèche, C. (1989). J. Phys. A Math. Gen. 22, L1163-L1166.]; Baake & Grimm, 2011[Baake, M. & Grimm, U. (2011). Philos. Mag. 91, 2661-2670.], 2013[Baake, M. & Grimm, U. (2013). Aperiodic Order, Vol. 1 of Encyclopedia of Mathematics and its Applications. Cambridge University Press.]; Socolar & Taylor, 2011[Socolar, J. E. & Taylor, J. M. (2011). J. Combin. Theory Ser. A, 118, 2207-2231.]). We show here that there exist limit-periodic tilings of four hyperuniformity classes: logarithmic (class II), weakly hyperuniform (class III), anti-hyperuniform, and an anomalous case in which Z(k) decays to zero faster than any power law as k goes to zero. The latter corresponds to a rule for which [\det{\bf M} = 0] (and [\lambda_2 = 0]), in which case α is not well defined. The existence of anti-hyperuniform limit-periodic tilings shows that anti-hyperuniformity does not require exotic singularities in S(k) for small k. Generally, it requires only that Z(k) grows sub-linearly with k.

6.1. The logarithmic case (α = 1)

The rule [L\rightarrow LSS], [S\rightarrow L] with S = 1 and L = 2 yields the well-known `period doubling' limit-periodic tiling. The eigenvalues of the substitution matrix are [\lambda_1 = 2] and [\lambda_2 = -1], leading to the prediction [\alpha = 1] and therefore quadratic scaling of Z(k) and logarithmic scaling of [\sigma^2(r)]. Numerical results for [\sigma^2(R)] are in good agreement with this prediction (Torquato et al., 2018[Torquato, S., Zhang, G. & De Courcy-Ireland, M. (2018). https://arxiv.org/abs/1804.06279]). In fact, one can show explicitly via direct calculation of [\sigma^2(R)] that the scaling is logarithmic. The calculation outlined in Appendix A[link] shows that

[\sigma^2(R) = {{1} \over {3}}\sum_{n = 0}^\infty \left\{{{w} \over {2^n}}\right\}\left(1-\left\{{{w} \over {2^n}}\right\}\right),\eqno(39)]

where w = 2R and {x} denotes the fractional part of x. From this it follows that for R = 2n-1/3 with [n\geq 1] we have

[\sigma^2(R) = {{2} \over {27}} \left({{13} \over {3}} + n \right),\eqno(40)]

demonstrating clear logarithmic growth for this special sequence of R values. One can also derive an upper bound over the interval [2^{n-1}\,\lt\,R\leq 2^n] by assuming that the summand in equation (39)[link] takes its maximum possible value on the intervals (0,1/2], (1/2,1], (2m-1,2m], for [m\leq n], and maximizing the possible sum of the exponentially decaying remaining contributions. The result is

[\sigma^2(R) \,\lt\, {{1} \over {4}} \left(1 + {{n} \over {3}}\right)\,\, {\rm for\ } 2^{n-1}\,\lt\,R\leq 2^n.\eqno(41)]

This upper bound also grows logarithmically and is shown as a series of dashed lines in Fig. 4[link].

[Figure 4]
Figure 4
The [\alpha = 1] (period doubling) limit-periodic tiling. Top: the tile boundaries with each point plotted at a height corresponding to the value of n for the sublattice to which it belongs. Bottom: plot of the number variance. The horizontal dotted line marks [\sigma^2 = 2/9], which is obtained for every R of the form 2n with integer [n\geq-1]. The dashed lines indicate upper bounds, and the open circles are analytically calculated values for R = 2n/3. See text for details.

It is instructive to carry out a more detailed analysis of S(k) for this particularly simple case as well. [See also Torquato et al. (2018[Torquato, S., Zhang, G. & De Courcy-Ireland, M. (2018). https://arxiv.org/abs/1804.06279]).] The tiling generated by applying the substitution rule repeatedly to a single L with its left edge at x = 1 consists of points located at positions 4l (2j+1), where l and j range over all positive integers (including zero). The structure factor therefore consists of peaks at [k_{mn} = 2\pi m/(a p^n)], with a = 2 and p = 4, for arbitrarily large n and all integer m. For m not a multiple of 4n-1, the peak at kmn gets nonzero contributions only from the lattices with [\ell \geq n]. These can be summed as follows:

[\eqalignno{S(k_{mn}) & = \lim_ {N\rightarrow\infty} \bigg| \sum_{\ell = n}^{\infty} \left( {{1}\over{2\times 4^{\ell}}}\right) &\cr & \quad\times {{1}\over{N}} \sum_{j=0}^{N-1} \exp\left[{{2\pi i m 4^{\ell}(2j+1)}\over{2\times4^n}}\right] \bigg|^2 &(42)}]

[= \left({{1}\over{9 \times 4^{2n}}}\right) 4^{\rm{mod}_2(m+1)},\eqno(43)]

where the factor of 1/(2×4l) in the first line is the density of the sublattice with that lattice constant. Applying this reasoning to each value of n gives a result that can be compactly expressed as

[S(k_{m\nu}) = \left[{{{\rm GCD}(2^{\nu},m)} \over {3 \times 4^{\nu}}} \right]^2,\eqno(44)]

where ν is an arbitrarily large integer, GCD() is the greatest common denominator function, and m can now take any positive integer value. Fig. 5[link] shows plots of S(k) and Z(k) for this tiling. [See also Torquato et al. (2018[Torquato, S., Zhang, G. & De Courcy-Ireland, M. (2018). https://arxiv.org/abs/1804.06279]) for an explicit expression for Z(k) and proof of the quadratic scaling.] Note that the apparent repeating unit in the plot of Z(k) spans only a factor of 2, even though the scaling factor for the lattice constants is 4. A similar effect occurs in the Poisson and anti-hyperuniform cases below. In the present case, the construction in Appendix A[link] showing that the density can be expressed using lattice constants 1/2n explains the origin of the effect.

[Figure 5]
Figure 5
The [\alpha = 1] limit-periodic tiling. Top: a logarithmic plot of the analytically computed S(k) (arbitrarily scaled) including kmn with [n\leq 3]. Bottom: a log–log plot of Z(k) computed numerically from S(k). The dashed red line shows the expected quadratic scaling law.

6.2. A Poisson scaling example (α = 0) and weak hyper­uniformity (0 < α < 1)

The substitution rule

[S\rightarrow LL, \quad L\rightarrow LLSSSS \eqno(45)]

with S = 1 and L = 2 produces a limit-periodic tiling with a = 2 and p = 16. Equation (28)[link] yields [\alpha = 0], which is the value corresponding to a Poisson system. Fig. 6[link] shows the result of direct computations of Z(k) including all of the Bragg peaks at [k = 2\pi n/ (a p^3)] and of [\sigma^2(R)]. Values of [\sigma^2] were computed from a sequence of 21 889 points obtained by seven iterations of the substitution rule on an initial L tile. For each point, a window of length 2R is moved continuously along the sequence for the computation of the averages.

[Figure 6]
Figure 6
Comparison of direct computation of Z(k) and [\sigma^2(R)] with the predicted scaling laws for a limit-periodic tiling with [\alpha = 0]. The dashed red lines show the expected linear scaling laws. The inset shows the piecewise parabolic behavior of [\sigma^2(R)] over a small span of R values.

Limit-periodic examples of weak hyperuniformity (class III) are afforded by substitutions of the form

[{\bf M} = \left[\matrix{ 0 & 2n \cr 2 & 2(n-1) }\right],\eqno(46)]

with [n\geq 3] with L/S = n, which yields

[\alpha = {{\ln n - \ln 2} \over {\ln n + \ln 2}} \eqno(47)]

[=\{0.226294, 1/3, 0.39794, 0.442114, \ldots\}.\eqno(48)]

6.3. Anti-hyperuniformity (α < 0)

The substitution rule

[S\rightarrow LLL,\, L\rightarrow LLLSSSSSS \eqno(49)]

with S = 1 and L = 2 produces a limit-periodic tiling with a = 2 and p = 36. Equation (28)[link] yields

[\alpha = 1-2{{\ln 3} \over {\ln 6}} = -0.226294\ldots,\eqno(50)]

which indicates anti-hyperuniform fluctuations. Fig. 7[link] shows the result of a direct computation of Z(k) including all of the Bragg peaks at [k = 2\pi n/ (a p^3)].

[Figure 7]
Figure 7
Comparison of direct computation of Z(k) with the predicted scaling law for a limit-periodic tiling with [\alpha = -0.226294\ldots]. The dashed red line shows the expected scaling law with slope [1+\alpha].

More generally, substitution matrices of the form

[{\bf M} = \left(\matrix{ 0 & 2 n \cr n & n }\right)\eqno(51)]

with [n\geq 3] and L/S = 2 yield limit-periodic anti-hyperuniform tilings with

[\alpha = 1 - 2{{\ln n} \over {\ln 2n}} = {{\ln 2 - \ln n} \over {\ln 2 + \ln n}}\eqno(52)]

[= \{-0.226294, -1/3, -0.39794, \ldots\}.\eqno(53)]

6.4. A λ2 = 0 case (α undefined)

A special class of tilings is derived from substitution matrices of dimension [{\cal D}] = 2 that have [\lambda_2 = 0] (and hence [\det{\bf M} = 0]). Such rules can produce periodic tilings, limit-periodic ones or more complex structures. The criteria for limit-periodicity can be obtained by analyzing constant-length substitution rules in which each L is considered to be made up of two tiles of unit length: L = AB. If the induced substitution rule on S, A and B exhibits appropriate coincidences, the tiling is limit-periodic (Dekking, 1978[Dekking, F. (1978). Z. Wahrscheinlichkeitstheorie verwandte Geb. 41, 221.]; Queffelec, 1995[Queffelec, M. (1995). In Beyond Quasicrystals, edited by F. Axel & D. Gratias, pp. 201-213. Berlin: Springer-Verlag.]). For the substitution matrix (1,1,2,2), the rule [(S\rightarrow SL]; [L\rightarrow SLSL)] produces a periodic tiling, and [(S\rightarrow SL]; [L\rightarrow SLLS)], for example, produces a limit-periodic tiling.

For the limit-periodic cases, the analysis above would suggest [\alpha \rightarrow \infty], or, more properly, α is not well defined. We present here an analysis of a particular case for which the convergence of Z(k) to zero is indeed observed to be faster than any power law.

The substitution rule

[S\rightarrow SL, L\rightarrow SLLS \eqno(54)]

with S = 1 and L = 2 produces a limit-periodic tiling with a = 1 and p = 3. Inspection of the point set (displayed in Fig. 8[link]) reveals that the number of points in the basis of each periodic subset for [n\ge 2] is 2n-2. The density of points in subset [n\ge 2] is (1/4)(2/3)n. The substitution matrix [{\bf M} = (1,1,2,2)] has eigenvalues [\lambda_1 = 3] and [\lambda_2 = 0].

[Figure 8]
Figure 8
Top: periodic sublattices of the limit-periodic point set generated by equation (54)[link]. Each point is plotted at a height n corresponding to the subset that contains it. Points of the same color form a periodic pattern with period 3n. Second: deviation of [|A(k_n)|] from 1/3n. Third: the integrated structure factor for the limit-periodic tiling with [\lambda_2 = 0], computed from subsets with [n\le 8]. The straight red (dashed) line of slope 5 is a guide to the eye for observing the concavity of the curve. Bottom: plot of the number variance for the limit-periodic tiling with [\lambda_2 = 0].

The unusual scaling in this case arises from interference effects associated with the form factors of the different periodic subsets. Let [k_n = 2\pi/3^n], the fundamental wavenumber for the nth subset, and let Xn denote the set of points in a single unit cell of the nth subset. S(kn) has contributions coming from all subsets of order n and higher. (Subsets of lower order do not contribute, as their fundamental wavenumber is larger than kn.) After some algebra, we find

[\eqalignno{ A(k_n) &= {{1}\over{3^n}} \sum_{x\in X_n} \exp(2\pi i x/3^n)&\cr &\quad + {{1}\over{3^{n+1}}}\left[\sum_{x\in X_{n+1}}\!\!\! \exp(2\pi i x/3^n) + \sum_{x\in X_{n+2}}\!\!\! \exp(2\pi i x/3^n) \right].&\cr &&(55)}]

Numerical evaluation of the sums over the unit-cell bases reveals that A(kn) is suppressed by the interference from subsets of higher order. Fig. 8[link] shows the behavior of the quantity [F_n = 3^n |A(k_n)|], revealing a rapid decay for small kn. The red (dashed) line shows the curve Fx = (1/3) (3x)-ln(9x)/2, which appears to fit the points well. An analytic calculation of Z(kn) is beyond our present reach. The middle panel of Fig. 8[link] shows the results of a numerical computation that includes all peaks [k = 2\pi m/p^6], with p = 3. It is clear that Z(k) is concave downwards on the log–log plot, consistent with the expectation that Z(k) goes to zero faster than any power of k. Note that the curve is not reliable for the smallest values of k due to the cutoff on the resolution of k values that are included. The deviation from power-law scaling is most easily seen in the increasing with n of the step sizes of the large jumps at [k = 2\pi/3^n]. (Compare with the constant step sizes in Figs. 5[link], 6[link] and 7[link].)

For completeness, the bottom panel of Fig. 8[link] also shows a plot of the number variance for this tiling. As expected, [\sigma^2(R)] is bounded from above. We note that the curve appears to be piecewise parabolic, which is also the case for the standard Fibonacci quasicrystal (Oğuz et al., 2017[Oğuz, E. C., Socolar, J. E. S., Steinhardt, P. J. & Torquato, S. (2017). Phys. Rev. B, 95, 054119.]), though the technique for calculating [\sigma^2(R)] based on projecting the tiling vertices from a 2D lattice is not applicable here.

7. Discussion

We have presented a heuristic method for calculating the hyperuniformity exponent α characterizing point sets generated by substitution rules that preserve the length ratios of the intervals between points. The calculation relies only on the relevant substitution matrix and an assumption that the tile order under substitution does not lead to a periodic tiling. The method performs well in that it yields a value of α consistent with direct measurements of the scaling of [\sigma^2(R)] in several representative cases. This allows for a straightforward construction of point sets with any value of α between -1 and 3.

It is well known that substitution rules can be divided into distinct classes corresponding to substitution matrices with eigenvalues that are not PV numbers leading to structure factors S(k) that are singular continuous (Bombieri & Taylor, 1986[Bombieri, E. & Taylor, J. E. (1986). J. Phys. Coll. 47, C3-C19.]; Baake et al., 2017[Baake, M., Frank, N. P., Grimm, U. & Robinson Jr, E. A. (2017). https://arxiv.org/abs/1706.03976.]), while substitution rules for which [|\lambda_2| \,\lt\, 1] yield Bragg peaks. Our analysis shows that this distinction corresponds to α greater than or less than unity, respectively. From the perspective of hyperuniformity, on the other hand, the critical value of α is zero, which corresponds to [|\lambda_2| = (\lambda_1)^{1/2}]. To achieve [\alpha\,\lt\, 0], a naïve comparison to scaling theories for systems with continuous spectra would suggest that S(k) must diverge for small k. We find, however, that anti-hyperuniformity, which does require sub-linear scaling of Z(k), can occur without any divergence both in cases where the spectrum is singular continuous, as for non-PV substitutions, and in cases where the spectrum consists of a dense set of Bragg peaks, as in some limit-periodic systems.

Finally, our investigations led us to consider the results of applying substitution rules for which [\lambda_2 = 0], which turned up a novel case of a limit-periodic tiling for which S(k) approaches zero faster than any power law. The physical implications of this type of scaling have yet to be explored.

The different tiling types and their hyperuniformity properties are summarized in Table 1[link]. Examples of quasiperiodic tilings in classes I and II are presented in Oğuz et al. (2017[Oğuz, E. C., Socolar, J. E. S., Steinhardt, P. J. & Torquato, S. (2017). Phys. Rev. B, 95, 054119.]). Note, however, that the class II case is not a substitution tiling. We do not know whether some other construction methods might yield quasiperiodic tilings that are in class III, anti-hyperuniform, or even anomalous. For non-PV tilings (which are substitution tilings by definition), at least two eigenvalues of the substitution matrix must be greater than unity, which rules out class II and class I. We conjecture that there are no limit-periodic tilings in class I. We can prove this for [{\cal D}] = 2 substitutions based on the fact that limit-periodicity requires the two eigenvalues to be rational and the fact that [{\bf M}] has only integer elements requires their sum and product to be integers, but we do not have a proof for [{\cal D} \,\gt\, 2].

Table 1
Types of 1D tilings and their possible hyperuniformity classes

A tick indicates that tilings of the given type exist, a dash that there are no such tilings, and a question mark that we are not sure whether such tilings exist.

  Anti-hyperuniform Weakly hyperuniform (class III) Logarithmically hyperuniform (class II) Strongly hyperuniform
  Class I Anomalous Gapped
  [-1\le \alpha \le 0] [0\,\lt\,\alpha \,\lt\, 1] [\alpha =1] [1\,\lt\,\alpha \le 3] [\alpha \rightarrow \infty] α irrelevant
Periodic [\surd]
Quasiperiodic ? ? [\surd] [\surd] ?
Non-PV [\surd] [\surd]
Limit-periodic [\surd] [\surd] [\surd] ? [\surd]

APPENDIX A

Calculation of σ2(R) for the period doubling limit-periodic tiling

The substitution rule [L\rightarrow LSS], [S\rightarrow L] (with S = 1 and L = 2) applied to an initial L with its left boundary at x = 1 produces a tiling with tile boundaries at all positions of the form xm,j = (2i + 1)4m, with j and m both running over all positive integers (including zero). We are interested in computing [\sigma^2(R)] for the density

[\rho(x) = \textstyle\sum\limits_{m,j = 0}^{\infty} \delta\left[x - (2j + 1)4^m\right].\eqno(56)]

Recall that [\sigma^2(R)] is the variance in the number of points covered by a window of length 2R placed with its left edge at x with uniform probability over all positive real values of x.

In Torquato et al. (2018[Torquato, S., Zhang, G. & De Courcy-Ireland, M. (2018). https://arxiv.org/abs/1804.06279]), an expression for [\sigma^2(R)] is derived using equation (4)[link] above. Here we show how [\sigma^2(R)] can be computed directly, thereby confirming the validity of equation (4)[link] for this limit-periodic system and arriving at a particularly simple expression that can be analyzed in detail.

We first note that we can rewrite [\rho(x)] as follows:

[\rho(x) = \textstyle\sum\limits_{n,i = 0}^{\infty} (-1)^n \delta (x - i 2^n).\eqno(57)]

To see this, first note that if x is an odd integer, then the only term that contributes is n = 0, i = x, which gives a +1. This is the m = 0 lattice of equation (56)[link]. More generally, if x is an odd multiple of 2p, there are contributions only from all [n \leq p], and these have alternating signs. If p is odd, the number of such contributions is even, yielding a density of zero. If p is even, the sum of the contributions is +1. The even values of p correspond to all integer values of m in equation (56)[link].

Let w = 2R be the length of the window and let Nn(x) be the number of points in the nth lattice covered by the window. It is convenient to take the window to be open at its left edge and closed at its right edge. Define {w}n as the fractional part of w/2n. It is convenient to think of w as being expressed in base 2. {w}n is then given by the first n digits to the left of the decimal point, plus all of the digits to the right. Note that Nn(x) depends on x only through {w}n; the integer part of w/2n adds the same number of points independent of the value of x. Furthermore, there are only two possible values of Nn(x), which differ by unity. For the purpose of computing the variance, we take these to be 0 and 1, and we work with the densities of these values rather than the full values of Nn.

If the window is placed with its left edge at x = j 2n, the contribution to the density from the nth lattice is 0. In order for the window to cover an additional point, the left edge must be placed such that [\{x\}_n \,\gt\, 1 - \{w\}_n]. The average density covered by the window of length w is thus

[\langle\mu\rangle = \textstyle\sum\limits_{n = 0}^{\infty} (-1)^n \{w\}_n.\eqno(58)]

The density squared is

[\rho^2(x) = \textstyle\sum\limits_{n,i = 0}^{\infty}\sum\limits_{\ell,j = 0}^{\infty} (-1)^{n+\ell} \delta (x - i 2^n) \delta (x - j 2^{\ell}).\eqno(59)]

A nonzero contribution to [\langle\mu^2\rangle] arises from an individual term if and only if both the n and l lattices contribute. For [\ell \,\gt\, n], the fraction of x's for which this is true is

[\{w\}_n\left[\{w\}_{\ell} + 2^{n-\ell} (1-\{w\}_{n})\right].\eqno(60)]

The first term accounts for window placements that give a contribution from the nth lattice. The light-gray bars in Fig. 9[link] show the values of x where the left edge of the window can be placed to produce a nonzero contribution. The term in parentheses counts the number of such intervals that occur within a region that contribute from the lth lattice (indicated by dark-gray bars in the figure) divided by the number of bars in one lattice spacing of the lth lattice. We thus have

[\eqalignno{ \langle\mu^2\rangle &= \textstyle\sum\limits_{n = 0}^{\infty} \{w\}_n (-1)^{2n} &\cr &\quad + 2 \textstyle\sum\limits_{n = 0}^{\infty} \sum\limits_{\ell \,\gt\, n}^{\infty} \{w\}_n\left[\{w\}_{\ell} + 2^{n-\ell} (1-\{w\}_{n})\right]. &(61) }]

Using (-1)2n = 1 and writing out [\langle\mu\rangle^2] as a sum over n plus a double sum with [\ell \,\gt\, n], straightforward algebra with convenient cancelations of some of the double sums yields

[\sigma^2 = \langle\mu^2\rangle - \langle\mu\rangle^2 \eqno(62)]

[\eqalignno{&=\sum_{n = 0}^{\infty} (\{w\}_n - \{w\}_n^2) & \cr &\quad + 2 \sum_{n = 0}^{\infty}\sum_{\ell = n+1}^{\infty} (\{w\}_n - \{w\}_n^2)\left({{-1} \over {2}}\right)^{\ell - n}&(63)}]

[={{1} \over {3}} \sum_{n = 0}^{\infty} (\{w\}_n - \{w\}_n^2).\eqno(64)]

This result has been confirmed to be in perfect agreement with direct computations.

[Figure 9]
Figure 9
Illustration for explaining the computation of a term in the double sum expression for [\langle\mu^2\rangle].

Equation (64)[link] describes a piecewise quadratic function of w (see Fig. 5[link]). One immediately sees that all values of w of the form 2l give the same result; they give {w}n = 0 for [n\leq \ell] and the same infinite series for [n \,\gt\, \ell]. Recalling that R = w/2, the shared value is

[\sigma^2(2^{n-1}) = {{1} \over {3}}\sum_{m = 1}^{\infty} \left({{1} \over {2}}\right)^m - \left({{1} \over {4}}\right)^m = {{2} \over {9}}.\eqno(65)]

To show that the upper envelope of [\sigma^2] grows logarithmically, we first prove an upper bound that grows only logarithmically, then identify a special sequence of window length values for which the growth is logarithmic. The upper bound is obtained by replacing all terms {w}n - {w}n2 in the sum with the maximum value 1/4 for all [n \,\lt\, 1+ \log_2 d], then replacing the remaining infinite series with its maximal value, obtained by maximizing [\sum_n (x/2^n - x^2/4^n)]. The result is

[\sigma^2(2^{n-2}\,\lt\,R\,\lt\,2^{n-1}) \,\lt\, {{1} \over {12}}(3+n).\eqno(66)]

To show that there is a sequence of R values for which [\sigma^2] grows logarithmically, consider w of the form 2n/3. Note that the binary representation of w is [101\ldots 01.0101 \ldots] for n odd and [\ldots 101\ldots 0.1010 \ldots] for n even. Again we consider the contributions from [\ell \leq n], then sum the remaining series. The value of {w}l oscillates between 1/3 and 2/3 for [\ell\,\lt\,n]. Straightforward algebra yields

[\sigma^2(R = 2^{n}/3) = {{2n} \over {27}} + {{26} \over {81}},\eqno(67)]

which clearly grows logarithmically with R. Note that the coefficient of n here is 2/27, reasonably close to the coefficient of 1/12 derived for the upper bound, and that, as shown in Fig. 5[link], these points are quite close to the true maxima for [w\,\lt\,2^n].

Acknowledgements

JESS thanks Michael Baake, Franz Gähler, Uwe Grimm and Lorenzo Sadun for helpful conversations at a workshop sponsored by the International Centre for Mathematical Sciences in Edinburgh. PJS thanks the Simons Foundation for its support and New York University for their generous hospitality during his leave at the Center for Cosmology and Particle Physics where this work was completed.

Funding information

ST was supported by the National Science Foundation under award No. DMR-1714722.

References

First citationAubry, S., Godrèche, C. & Vallet, F. (1987). J. Phys. Fr. 48, 327–334.  CrossRef Google Scholar
First citationBaake, M., Frank, N. P., Grimm, U. & Robinson Jr, E. A. (2017). https://arxiv.org/abs/1706.03976Google Scholar
First citationBaake, M., Gaehler, F. & Manibo, N. (2018). See Section 3.4, Equation (18). https://arxiv.org/abs/1805.09650Google Scholar
First citationBaake, M. & Grimm, U. (2011). Philos. Mag. 91, 2661–2670.  CrossRef Google Scholar
First citationBaake, M. & Grimm, U. (2012). Chem. Soc. Rev. 41, 6821–6843.  Web of Science CrossRef CAS PubMed Google Scholar
First citationBaake, M. & Grimm, U. (2013). Aperiodic Order, Vol. 1 of Encyclopedia of Mathematics and its Applications. Cambridge University Press.  Google Scholar
First citationBaake, M., Grimm, U. & Manibo, N. (2018). https://arxiv.org/abs/1709.09083Google Scholar
First citationBombieri, E. & Taylor, J. E. (1986). J. Phys. Coll. 47, C3–C19.  Google Scholar
First citationCastro-Lopez, M., Gaio, M., Sellers, S., Gkantzounis, G., Florescu, M. & Sapienza, R. (2017). APL Photonics, 2, 061302.  Google Scholar
First citationDekking, F. (1978). Z. Wahrscheinlichkeitstheorie verwandte Geb. 41, 221.  CrossRef Google Scholar
First citationDreyfus, R., Xu, Y., Still, T., Hough, L. A., Yodh, A. G. & Torquato, S. (2015). Phys. Rev. E, 91, 012302.  CrossRef Google Scholar
First citationFrank, N. P. (2008). Expo. Math. 26, 295–326.  CrossRef Google Scholar
First citationGodrèche, C. (1989). J. Phys. A Math. Gen. 22, L1163–L1166.  Google Scholar
First citationGodreche, C. & Luck, J. M. (1990). J. Phys. A Math. Gen. 23, 3769–3797.  CrossRef Google Scholar
First citationGodrèche, C. & Luck, J. M. (1992). Phys. Rev. B, 45, 176–185.  Google Scholar
First citationHaberko, J. & Scheffold, F. (2013). Opt. Express, 21, 1057–1065.  CrossRef Google Scholar
First citationHexner, D. & Levine, D. (2015). Phys. Rev. Lett. 114, 110602.  CrossRef Google Scholar
First citationKesten, H. (1966). Acta Arith. 12, 193–212.  CrossRef Google Scholar
First citationMan, W., Florescu, M., Williamson, E. P., He, Y., Hashemizad, S. R., Leung, B. Y. C., Liner, D. R., Torquato, S., Chaikin, P. M. & Steinhardt, P. J. (2013). Proc. Natl Acad. Sci. USA, 110, 15886–15891.  CrossRef Google Scholar
First citationOğuz, E. C., Socolar, J. E. S., Steinhardt, P. J. & Torquato, S. (2017). Phys. Rev. B, 95, 054119.  Google Scholar
First citationQueffelec, M. (1995). In Beyond Quasicrystals, edited by F. Axel & D. Gratias, pp. 201–213. Berlin: Springer-Verlag.  Google Scholar
First citationSocolar, J. E. & Taylor, J. M. (2011). J. Combin. Theory Ser. A, 118, 2207–2231.  CrossRef Google Scholar
First citationTorquato, S. (2018). Phys. Rep. 745, 1–95.  CrossRef Google Scholar
First citationTorquato, S. & Stillinger, F. H. (2003). Phys. Rev. E, 68, 041113.  Web of Science CrossRef Google Scholar
First citationTorquato, S., Zhang, G. & De Courcy-Ireland, M. (2018). https://arxiv.org/abs/1804.06279  Google Scholar
First citationTorquato, S., Zhang, G. & Stillinger, F. H. (2015). Phys. Rev. X, 5, 021020.  Google Scholar
First citationZachary, C. E. & Torquato, S. (2009). J. Stat. Mech. 2009, P12015.  Google Scholar

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.

Journal logoFOUNDATIONS
ADVANCES
ISSN: 2053-2733
Follow Acta Cryst. A
Sign up for e-alerts
Follow Acta Cryst. on Twitter
Follow us on facebook
Sign up for RSS feeds