Hyperuniformity and anti-hyperuniformity in one-dimensional substitution tilings

This work examines the long-wavelength scaling properties of self-similar substitution tilings, placing them in their hyperuniformity classes. Quasiperiodic, non-PV (Pisot–Vijayaraghavan number) and limit-periodic examples are analyzed. Novel behavior is demonstrated for certain limit-periodic cases.


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). Crystals and quasicrystals are hyperuniform, as are a variety of disordered systems, including certain equilibrium structures, products of nonequilibrium selfassembly protocols and fabricated metamaterials. [For examples, see Man et al. (2013), Haberko & Scheffold (2013), Dreyfus et al. (2015), Torquato et al. (2015), Hexner & Levine (2015), Castro-Lopez et al. (2017), Torquato (2018).] 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). 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 onedimensional (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 AEk 0 n=p m , where n, m and p are positive integers (Godrè che, 1989;Baake & Grimm, 2011. 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;Godrè che & Luck, 1992;Baake et al., 2017). (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). Finally, there are cases for which SðkÞ is absolutely continuous (Baake & Grimm, 2012) 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 reviews the definition of the scaling exponent associated with both ZðkÞ and the variance 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 reviews the substitution method for creating tilings, using the well known Fibonacci tiling as an illustrative example. The substitution matrix M is defined and straightforward results for tile densities are derived. Section 4 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): strongly hyperuniform (class I), weakly hyperuniform (class III) and anti-hyperuniform. Section 5 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 considers substitutions that produce limit-periodic tilings. Examples are presented of four distinct classes: logarithmic hyperuniform (class II), weakly hyperuniform (class III), antihyperuniform, and an anomalous class in which ZðkÞ approaches zero faster than any power law. Finally, Section 7 provides a summary of the key results, including a table showing which types of tilings can exhibit the various classes of (anti-)hyperuniformity.

Classes of hyperuniformity
For systems having a structure factor SðkÞ that is a smooth function of the wavenumber k, SðkÞ tends to zero as k tends to zero (Torquato & Stillinger, 2003), typically scaling as a power law: 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: The factor of 2 is inserted for consistency with higherdimensional generalizations where q is treated as a radial coordinate. In both cases, may be defined by the relation (Og uz et al., 2017) ZðkÞ $ k 1þ as k ! 0: Systems with > 0 have long-wavelength spatial fluctuations that are suppressed compared with Poisson point sets and are said to be hyperuniform (Torquato & Stillinger, 2003). Prototypical strongly hyperuniform systems (with > 1) include crystals and quasicrystals. We refer to systems with < 0 as anti-hyperuniform (Torquato, 2018). 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 2 ðRÞ. If 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 2 ðRÞ is closely related to the behavior of ZðkÞ for small k (Torquato & Stillinger, 2003;Og uz et al., 2017). For a general point configuration in 1D with a well-defined average number density , 2 ðRÞ can be expressed in terms of SðkÞ and the Fourier transform ðk; RÞ of a uniform density interval of length 2R: For any 1D system with a smooth or quasicrystalline structure factor, the scaling of 2 ðRÞ for large R is determined by as follows (Torquato & Stillinger, 2003;Zachary & Torquato, 2009;Torquato, 2018 : For hyperuniform systems, we have > 0, and the distinct behaviors of 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 < 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 (Og uz et al., 2017). Other quasiperiodic point sets (not obtainable by substitution) are known to belong to class II (Kesten, 1966;Aubry et al., 1987;Og uz et al., 2017).

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 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: which acts on the column vector ðN S ; N L Þ 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 , 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. Given a finite sequence with N S tiles of length S and N L 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 ðN S ; N L Þ.
More generally, substitution rules can be defined for systems with more than two tile types, leading to substitution matrices with dimension D greater than 2. We present explicit reasoning here only for the D = 2 case. A substitution rule for two tile types is characterized by a substitution matrix: The associated rule may be the following: 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 produces the periodic tiling . . . SLSLSL . . ., while the rule produces the more complicated sequence discussed below in Section 6. 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 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 n S short intervals and n L long ones. We are interested in point sets formed by decorating each L tile with ' points and each S tile with s points. The total number of points at the mth iteration is and the length of the tiling at the same step is Let 1 and 2 be the eigenvalues of M, with 1 being the largest, and let v 1 and v 2 be the associated eigenvectors. We have The unit vectors ð1; 0Þ and ð0; 1Þ may be expressed as follows: 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.
The density of tile vertices after m iterations, m ¼ N m =X m , is thus

Scaling properties of 1D substitution tilings
As long as the coefficient of ð 2 = 1 Þ m in equation (22) does not vanish, the deviations of from for portions of the tiling that are mapped into each other by substitution are related by 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).
Conjecture 1. We take equation (23) 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 in equation (23) can be interpreted as that characteristic amplitude.
Conjecture 2. For the present purposes, we define the realspace density gðxÞ ¼ P n ðx À x n Þ to consist of unit-strength -functions placed at every point x n where two tiles meet, and consider the Fourier amplitudes 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: This implies the form AðkÞ $ k ðÀ ln j 2 = 1 j= ln j 1 jÞ ¼ k 1Àðln j 2 j= ln j 1 jÞ : Squaring to get SðkÞ, we have This conjecture may not hold when interference effects are important, as in the case discussed in Section 6 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Þ (Og uz et al., 2017). 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;Godreche & Luck, 1990), we conjecture that it remains sufficiently singular for the relation to hold. Thus we immediately obtain Note that this calculation of the scaling exponent makes no reference to the distinction between substitutions with j 2 j < 1 and those with j 2 j > 1. In the former case, 1 is a Pisot-Vijayaraghavan (PV) number, SðkÞ consists of Bragg peaks and 2 ðRÞ remains bounded for all R. In the latter case, the form of SðkÞ is more complex (Bombieri & Taylor, 1986), and quantities closely related to 2 ðRÞ, including the 'wandering exponent' associated with lifts of the sequence onto a higherdimensional hypercubic lattice, are known to show nontrivial scaling exponents (Godreche & Luck, 1990). From equation (28), we see that the hyperuniformity condition > 0 requires j 2 j < ð 1 Þ 1=2 . Though the result was obtained for substitutions with only D = 2 tile types, it holds for D > 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 j 2 j= 1 . This distinction between hyperuniform and antihyperuniform 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), 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 > 2, does not detect any qualitative difference between the cases n ¼ 3 and n ¼ 5. The former case is hyperuniform, with ¼ ð1=2Þ½1 AE ð13Þ 1=2 and ' 0:37, while the latter is anti-hyperuniform, with ¼ ð1=2Þ½1 AE ð21Þ 1=2 and ' À0:14.
For the Fibonacci case, we have 1 ¼ and 2 ¼ À1=, yielding ¼ 3, which agrees with the explicit calculation in Og uz et al. (2017). 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 1 ¼ n and 2 ¼ Àn=, yielding For n ! 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 n 4, our calculation predicts 0 < < 1 and hence 2 ðRÞ $ R 1À . 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 ¼ ' ¼ 1). Fig. 2 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 2 ðRÞ $ R 0:36 , in close agreement with the predicted exponent from equation (27): 1 À ¼ 2ðln 2 À ln Þ=ðln 2 þ ln Þ ' 0.36094. For n ! 5, the calculated value of is negative, approaching À1 as n approaches infinity. The point set is therefore antihyperuniform; it contains density fluctuations at long wavelengths that are stronger than those of a Poisson point set. For n ¼ 5, we have ¼ À0:0793 . . .. Fig. 3 shows a log-log plot of the computed number variance along with the line corresponding to 2 ðRÞ $ R 1À . 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 > À 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) and Godreche & Luck (1990) 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 below.
For rules that yield rational values of the length ratio , the coefficient of ð 2 = 1 Þ m in equation (22) can vanish for appropriate choices of n S and n L , 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 ! L and L ! SLS, which generates the periodic sequence . . . SLSLSL . . ., with ¼ 2, 1 ¼ 2 and 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 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 ¼ 1 and therefore a logarithmic scaling of 2 ðRÞ. This case is presented in more detail in Section 6 below, and the logarithmic scaling is confirmed.

Achievable values of a
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 M is full rank, always lies between À1 and 3.
First, note that the maximum value of j 2 = 1 j is 1, by definition, which sets the lower bound on via equation (28). The upper bound on is obtained when j 2 j is as small as possible, but there is a limit on how small this can be. The product of the eigenvalues of M is equal to det M, so j 2 j cannot be smaller than ðj det Mj= 1 Þ 1=ðDÀ1Þ . But j det Mj is an integer, and the smallest nonzero value it can take is 1. (The case 2 ¼ 0 is discussed in Section 6 below. For D ! 3, one can have det M ¼ 0 with nonzero 2 . The analysis of such cases is beyond our present scope.) Hence we have 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 À ' 0:36.

Figure 3
Log-log plot of the number variance (black dots) for an antihyperuniform 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 À ' 1:08.
Thus the maximum value of obtainable by this construction method is 3, which can occur for D = 2, as in the Fibonacci case.
The family of substitutions considered in Section 4 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 with d > a þ 1 and 2c < ðd À aÞ, we have 1 ¼ a and 2 ¼ d.
Note that ¼ ð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 For fixed a, d can range from a þ 2 to 1. 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 M's produces 's between 1 and 3. For with n > b, we have We thus obtain ¼ 1 À 2 ln b À ln 2 þ ln½ðn 2 þ 4Þ 1=2 À n ln b À ln 2 þ ln½ðn 2 þ 4Þ 1=2 þ n : For large n, we have which approaches 3 for b ( 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.

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 ap n , where p is an integer and n runs over all positive definite integers (Godrè che, 1989; Baake & Grimm, 2011Socolar & Taylor, 2011). We show here that there exist limitperiodic 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 M ¼ 0 (and 2 ¼ 0), in which case is not well defined. The existence of anti-hyperuniform limitperiodic 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 (a = 1) The rule L ! LSS, S ! L with S ¼ 1 and L ¼ 2 yields the well-known 'period doubling' limit-periodic tiling. The eigenvalues of the substitution matrix are 1 ¼ 2 and 2 ¼ À1, leading to the prediction ¼ 1 and therefore quadratic scaling of ZðkÞ and logarithmic scaling of 2 ðrÞ. Numerical results for 2 ðRÞ are in good agreement with this prediction . In fact, one can show explicitly via direct calculation of 2 ðRÞ that the scaling is logarithmic. The calculation outlined in Appendix A shows that where w ¼ 2R and fxg denotes the fractional part of x. From this it follows that for R ¼ 2 nÀ1 =3 with n ! 1 we have demonstrating clear logarithmic growth for this special sequence of R values. One can also derive an upper bound over the interval 2 nÀ1 < R 2 n by assuming that the summand The ¼ 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 2 ¼ 2=9, which is obtained for every R of the form 2 n with integer n ! À1. in equation (39) takes its maximum possible value on the intervals ð0; 1=2; ð1=2; 1; ð2 mÀ1 ; 2 m , for m n, and maximizing the possible sum of the exponentially decaying remaining contributions. The result is This upper bound also grows logarithmically and is shown as a series of dashed lines in Fig. 4. 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).] 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 4 ' ð2j þ 1Þ, where ' and j range over all positive integers (including zero). The structure factor therefore consists of peaks at k mn ¼ 2m=ðap n Þ, with a ¼ 2 and p ¼ 4, for arbitrarily large n and all integer m. For m not a multiple of 4 nÀ1 , the peak at k mn gets nonzero contributions only from the lattices with ' ! n. These can be summed as follows: where the factor of 1=ð2 Â 4 ' Þ 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 where is an arbitrarily large integer, GCDðÞ is the greatest common denominator function, and m can now take any positive integer value. Fig. 5 shows plots of SðkÞ and ZðkÞ for this tiling. [See also Torquato et al. (2018) 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 antihyperuniform cases below. In the present case, the construction in Appendix A showing that the density can be expressed using lattice constants 1=2 n explains the origin of the effect.
6.2. A Poisson scaling example (a = 0) and weak hyperuniformity (0 < a < 1) The substitution rule with S ¼ 1 and L ¼ 2 produces a limit-periodic tiling with a ¼ 2 and p ¼ 16. Equation (28) yields ¼ 0, which is the value corresponding to a Poisson system. The ¼ 1 limit-periodic tiling. Top: a logarithmic plot of the analytically computed SðkÞ (arbitrarily scaled) including k mn with n 3. Bottom: a log-log plot of ZðkÞ computed numerically from SðkÞ. The dashed red line shows the expected quadratic scaling law.

Figure 6
Comparison of direct computation of ZðkÞ and 2 ðRÞ with the predicted scaling laws for a limit-periodic tiling with ¼ 0. The dashed red lines show the expected linear scaling laws. The inset shows the piecewise parabolic behavior of 2 ðRÞ over a small span of R values.
peaks at k ¼ 2n=ðap 3 Þ and of 2 ðRÞ. Values of 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. Limit-periodic examples of weak hyperuniformity (class III) are afforded by substitutions of the form with n ! 3 with L=S ¼ n, which yields ¼ ln n À ln 2 ln n þ ln 2 ð47Þ ¼ f0:226294; 1=3; 0:39794; 0:442114; . . .g: 6.3. Anti-hyperuniformity (a < 0) The substitution rule with S ¼ 1 and L ¼ 2 produces a limit-periodic tiling with a ¼ 2 and p ¼ 36. Equation (28) yields which indicates anti-hyperuniform fluctuations. Fig. 7 shows the result of a direct computation of ZðkÞ including all of the Bragg peaks at k ¼ 2n=ðap 3 Þ. More generally, substitution matrices of the form with n ! 3 and L=S ¼ 2 yield limit-periodic antihyperuniform tilings with ¼ 1 À 2 ln n ln 2n ¼ ln 2 À ln n ln 2 þ ln n ð52Þ ¼ fÀ0:226294; À1=3; À0:39794; . . .g: 6.4. A k 2 = 0 case (a undefined) A special class of tilings is derived from substitution matrices of dimension D = 2 that have 2 ¼ 0 (and hence det M ¼ 0). Such rules can produce periodic tilings, limitperiodic 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 Comparison of direct computation of ZðkÞ with the predicted scaling law for a limit-periodic tiling with ¼ À0:226294 . . .. The dashed red line shows the expected scaling law with slope 1 þ .

Figure 8
Top: periodic sublattices of the limit-periodic point set generated by equation (54). 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 3 n . Second: deviation of jAðk n Þj from 1=3 n . Third: the integrated structure factor for the limit-periodic tiling with 2 ¼ 0, computed from subsets with n 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 2 ¼ 0.
For the limit-periodic cases, the analysis above would suggest ! 1, 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 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) reveals that the number of points in the basis of each periodic subset for n ! 2 is 2 nÀ2 . The density of points in subset n ! 2 is ð1=4Þð2=3Þ n . The substitution matrix M ¼ ð1; 1; 2; 2Þ has eigenvalues 1 ¼ 3 and 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=3 n , the fundamental wavenumber for the nth subset, and let X n denote the set of points in a single unit cell of the nth subset. Sðk n Þ 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 k n .) After some algebra, we find Numerical evaluation of the sums over the unit-cell bases reveals that Aðk n Þ is suppressed by the interference from subsets of higher order. Fig. 8 shows the behavior of the quantity F n ¼ 3 n jAðk n Þj, revealing a rapid decay for small k n . The red (dashed) line shows the curve F x ¼ ð1=3Þð3xÞ À lnð9xÞ=2 , which appears to fit the points well. An analytic calculation of Zðk n Þ is beyond our present reach. The middle panel of Fig. 8 shows the results of a numerical computation that includes all peaks k ¼ 2m=p 6 , with p ¼ 3. It is clear that ZðkÞ is concave downwards on the log-log plot, consistent with the expecta-tion 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=3 n . (Compare with the constant step sizes in Figs. 5, 6 and 7.) For completeness, the bottom panel of Fig. 8 also shows a plot of the number variance for this tiling. As expected, 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 (Og uz et al., 2017), though the technique for calculating 2 ðRÞ based on projecting the tiling vertices from a 2D lattice is not applicable here.

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 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;Baake et al., 2017), while substitution rules for which j 2 j < 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 j 2 j ¼ ð 1 Þ 1=2 . To achieve < 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 antihyperuniformity, 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,  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.
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 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. Examples of quasiperiodic tilings in classes I and II are presented in Og uz et al. (2017). 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, antihyperuniform, 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 D = 2 substitutions based on the fact that limit-periodicity requires the two eigenvalues to be rational and the fact that M has only integer elements requires their sum and product to be integers, but we do not have a proof for D > 2.
APPENDIX A Calculation of r 2 (R) for the period doubling limitperiodic tiling The substitution rule L ! LSS, S ! 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 x m;j ¼ ð2i þ 1Þ4 m , with j and m both running over all positive integers (including zero). We are interested in computing 2 ðRÞ for the density Recall that 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), an expression for 2 ðRÞ is derived using equation (4) above. Here we show how 2 ðRÞ can be computed directly, thereby confirming the validity of equation (4) 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 ðxÞ as follows: ðxÞ ¼ P 1 n;i¼0 ðÀ1Þ n ðx À i2 n Þ: 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). More generally, if x is an odd multiple of 2 p , there are contributions only from all n 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).
Let w ¼ 2R be the length of the window and let N n ð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 fwg n as the fractional part of w=2 n . It is convenient to think of w as being expressed in base 2. fwg 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 N n ðxÞ depends on x only through fwg n ; the integer part of w=2 n adds the same number of points independent of the value of x. Furthermore, there are only two possible values of N n ð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 N n .
If the window is placed with its left edge at x ¼ j2 n , 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 fxg n > 1 À fwg n . The average density covered by the window of length w is thus The density squared is 2 ðxÞ ¼ P 1 n;i¼0 P 1 ';j¼0 ðÀ1Þ nþ' ðx À i2 n Þðx À j2 ' Þ: ð59Þ A nonzero contribution to h 2 i arises from an individual term if and only if both the n and ' lattices contribute. For ' > n, the fraction of x's for which this is true is fwg n fwg ' þ 2 nÀ' ð1 À fwg n Þ Â Ã : ð60Þ The first term accounts for window placements that give a contribution from the nth lattice. The light-gray bars in Fig. 9 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 'th lattice (indicated by dark-gray bars in the figure) divided by the number of bars in one lattice spacing of the 'th lattice. We thus have Using ðÀ1Þ 2n ¼ 1 and writing out hi 2 as a sum over n plus a double sum with ' > n, straightforward algebra with convenient cancelations of some of the double sums yields