Gaussian orthogonal ensemble for quasiperiodic tilings without unfolding: r -value statistics

We study the level-spacing statistics for non-interacting Hamiltonians defined on the two-dimensional quasiperiodic Ammann--Beenker (AB) tiling. When applying the numerical procedure of"unfolding", these spectral properties in each irreducible sector are known to be well-described by the universal Gaussian orthogonal random matrix ensemble. However, the validity and numerical stability of the unfolding procedure has occasionally been questioned due to the fractal self-similarity in the density of states for such quasiperiodic systems. Here, using the so-called $r$-value statistics for random matrices, $P(r)$, for which no unfolding is needed, we show that the Gaussian orthogonal ensemble again emerges as the most convincing level statistics for each irreducible sector. The results are extended to random-AB tilings where random flips of vertex connections lead to the irreducibility.

The statistical description of energy levels in complex systems has a long and distinguished history [1,2]. For interacting systems such as, e.g., highly excited heavy nuclei, already Bethe recognized the intrinsic difficulty in obtaining such statistics [1] where the spacing of levels is often complicated and the Hilbert space exponentially large. Much progress has been made, nevertheless, when ignoring details of the underlying Hamiltonian and instead assuming a random matrix structure [3][4][5][6][7]. In such a situation, it is the invariance of the matrix under specific symmetry operations that determines the functional form of the level spacing distribution P(s), with Gaussian orthogonal, unitary, and symplectic ensembles (GOE, GUE, GSE, respectively) being the most famous examples [7].
Difficulties in using level statistics also arise for quasiperiodic (QP) systems when electronic degrees of freedom are described by noninteracting tight-binding Hamiltonians defined on QP tilings. It is well known that for such systems, the density of states (DOS) typically has self-similar (fractal) characteristics while a straightforward comparison with the Gaussian ensembles requires a flat DOS. A procedure known as "unfolding" is often used to convert the fractal DOS into the required flat behavior but this of course represents a severe change in the energy spectrum. Indeed, it has been shown that P(s) can change when only partial unfolding is being done [8]. Using the integrated DOS instead as the basis for an unfolding procedure yields more stable results which demonstrate that P(s) for a time-invariant and spin-independent QP Hamiltonian is very well described by P GOE (s) [9]. In fact, the P(s) was shown to follow P GOE (s) better than the celebrated Wigner expression P Wigner (s) = π/2 exp (−π s 2 /4) of GOE which is based on a (2×2)-matrix surmise [2].
Removing the dependence on the unfolding procedure, Oganesyan and Huse [10] introduced the so-called r-value distribution in the context of disordered many-body systems. With s n = E n − E n−1 denoting the spacing of two consecutive energy levels E n , E n−1 , they define r n as 0 r n = min{s n , s n−1 }/max{s n , s n−1 } 1. (1) For an uncorrelated Poisson spectrum, one has P Poisson = 2/(1 + r) 2 with mean r Poisson = 0.386 while for the GOE, with r GOE ∼ 0.5307(1) [11]. Expression (2) has the same status as the Wigner surmise quoted above, i.e., has been derived for the smallest possible (3×3) matrix in the GOE while r GOE is based on high-precision numerics. Here, we shall also use a surmise based on a (5×5) matrix (see Appendix) which improves upon (2) such that the deviation to r GOE is reduced from 1% to < 0.4%. In the context of many-body localization (MBL), r-value statistics has proven its worth by allowing the numerical determination of the transition from MBL to the so-called ergodic phase at weak disorders without the need to unfold spectra [10,12,13].
In this Letter, we employ the r-value statistics to QP Hamiltonians defined on the Ammann-Beenker (AB) octagonal tiling [14] and also to randomized versions of the tiling  Table I. In both panels, the dashed (green) lines gives P Poisson (r) while the solid (blue) lines represent the improved P GOE (r) and the dotted (black) line is (2). when certain connections have been allowed to flip. As we show in Fig. 1, in both cases, we find that the computed P(r) agrees very well with the predictions of the GOE when compared to P GOE (r). Hence, within the numerical accuracy available to us, we can conclude that the GOE ensemble is indeed the correct descriptor of level statistics for noninteracting tight-binding hopping Hamiltonians in QP tilings, whether that result has been computed (i) by unfolding DOS [8] or integrated density of state (IDOS) [9], (ii) by restricting the analysis to regions in the spectrum that have a flat DOS and hence do not need unfolding [15], or (iii) by circumventing unfolding with the r-value statistics.
As in Ref. [9], we shall consider the octagonal (or Ammann-Beenker) [14,16] tiling consisting of squares and rhombi with equal edge lengths as in Fig. 2(a) (see [17] for more on this tiling and its properties). Besides this perfect quasicrystal, we shall also study a randomized version in which triangular connections in the rhombi are allowed to flip randomly. Such structures are often used to model imperfect quasicrystals (see, e.g., [18]). Both the AB tilings and the random-AB tilings increase in size exponentially via infla- where I ∈ N denotes the generation of the inflation with I = 0 corresponding to the initial patch (the seed). We include data for patches of the AB tiling corresponding to I = 2, 3, 4, and 5 inflation steps with 833, 4713, 27 137, and 157 369 vertices, respectively. For the random-AB tiling, we use I = 1, 2, 3, and 4 inflation steps with 82, 478, 2786, and 16 238 vertices, respectively. On these tilings, we define H = i, j |i j| as a Hamiltonian with free boundary conditions for the AB tilings and periodic boundary conditions for the random-AB tilings. Here, |i is indicating the Wannier state at vertex i while pairs   of r-value estimates and RMSD for (top) AB and (bottom) random-AB tilings for different inflation levels I. The column "sector/flips" gives the phase (and parity) of the subspectrum for a particular AB tiling, while for the random-AB tilings, it indicates the chosen average value of flips per vertex. Labels indicate whether the statistics of different sectors were "combined" for r values from different sectors or whether "mult"(iple) energy spectra were analyzed together. N (E = 0) indicates the number of nonzero energy levels, while N (r) counts the number of r values used to construct P(r). The column "parts/samples" labels the number of subspectra for each AB tiling and the total number of spectra for the random-AB tilings. The average is given by r = 1 0 r p(r)dr and RMSD GOE/Poisson show the RMSD with respect to either GOE or Poisson ensemble. The error estimates for AB tilings are computed as error of mean; for the random-AB tilings, these are less than 10 −6 and not indicated here. The columns headed by % show how the computed RMSD GOE/Poisson deviate from RMSD max . The averages for AB tilings with combined statistics for I = 5 and for the random-AB tiling with I = 4 are highlighted in bold.  Fig. 2(a) has the symmetry of the regular octagon, corresponding to the dihedral group D 8 . Hence the Hamiltonian matrix splits into ten blocks according to the irreducible representations of D 8 : Using the rotational symmetry, one obtains eight blocks, two of which split further under reflection, while the remaining six form three pairs with identical spectra. This gives a total of seven independent subspectra. The starting point in the generation of the random-AB tilings is the perfect periodic approximant shown in Fig. 2(b). Note that this is the periodic approximant of highest exact symmetry for this tiling, since eightfold symmetry is forbidden in a periodic structure. As such, it has D 4 symmetry and five independent spectra. We can introduce randomness by flipping the arrangement of hexagonal patches consisting of a square on two rhombi that meet in a three-valent vertex. These "simpleton flips" are ergodic in the sense that repeated applications of this flip explore the entire ensemble of random tilings of square and rhombi with the given ratio of the two tiles. We shall study cases with an increasing number of flips per vertex. Note that the flips generally will remove any exact symmetries, such that the whole matrix becomes an irreducible block, while the statistical eightfold symmetry (in the sense that local configurations are equally likely to appear in any of the eight directions) is maintained (see [17] for details). An example is given in Fig. 2(c). For both AB tilings and random-AB tilings, the Hamiltonian is diagonalized independently for each of the irreducible spectra; each spectrum is symmetric about E = 0, because of the bipartiteness of the AB and random-AB tilings. Furthermore, a finite fraction of the states is degenerate at E = 0, corresponding to compactly localized states [19][20][21]; they do not contribute to the universal statistics, and we neglect them.

I
As is well known, the DOS for each irreducible spectrum is rather spiky [8], while the IDOS is already rather smooth [9] (results not shown here). We proceed without any unfolding. Only eigenvalues |E n | > 10 −10 are included in the further analysis, and in Table I we give the number of these as N (E = 0). A further restriction to positive E n 's also removes the double degeneracy, for the s n , resulting from the bipartiteness of the tilings. This leads to the available number of r n values given in Table I as N (r). For each spectrum, we independently compute P(r) and r . For L060201-3 the AB tilings, we find that the P(r) distributions for all I show level repulsion such that P(r) ∝ r for r → 0 [11]. With increasing I, the slope of the level repulsion decreases slightly, level repulsion increases and rapidly approaches the small-r behavior of P GOE (r). Similarly, the bulk behavior of P(r) follows P GOE (r) ever more closely for increasing I. In order to ascertain quantitatively how well the estimated P(r) follows either P Poisson (r) or P GOE (r), we establish the root-mean-squared deviation (RMSD) defined as { ∞ 0 [P(r) − Q(r)] 2 dr} 1/2 with Q(r) either P Poisson (r) or P GOE (r). In Table I we also express the RMSD values in percentage by comparing with the RMSD max = { ∞ 0 [P Poisson (r) − P GOE (r)] 2 dr} 1/2 = 0.630508 between P Poisson (r) and P GOE (r). We find for P(r) that while the RMSD to GOE values, RMSD GOE , are roughly in the ∼10%-20% range, the RMSD Poisson values are at ∼90%-100% of RMSD max . Hence for each subspectrum, we see that P(r) already very nicely follows P GOE (r) while P Poisson (r) is certainly ruled out. This conclusion is also corroborated when studying r for the AB tilings with all estimates within 2% of r GOE ∼ 0.5307(1). The best agreement is found when we combine all r values for the largest system I = 5. The resulting P(r) is the one shown in Fig. 1 with r = 0.526 654(1).
Due to the exponential growth of the size of the Hamiltonian matrix with I, a further increase of I is computationally very challenging (for I = 6, we have N = 1, 657, 756, 990). We therefore now turn to the random-AB tilings introduced above where we can increase N (r) simply by computing many random realizations. In these systems, the DOS is also somewhat less spiky, but still retains considerable variation across the energy spectrum that would still require significant unfolding when studying P(s). We summarize the results in Table I and show the behavior of P(r) in Fig. 1. For I = 3, we give results for all 1000 samples when flipping each triangular connection [cf. Fig. 2(c)] on average 1, 10, 100, or 1000 times for a thoroughly randomized tiling. We find that the differences in P(r) between these cases, even with just a single flip per vertex (on average), are very small, and no major influence of the underlying exact D 4 symmetry of the approximant can be seen anymore. We therefore also present "combined" statistics where r values for 10, 100, and 1000 flips have been analyzed together. For I = 2 we only show these combined results while full details are given also for I = 4. In this case, the computational effort is already considerable for each sample so that only 100 samples have been calculated for all flips. The overall combined N (r) = 699 579, 4 132 302, and 2 421 369 for I = 2, 3, and 4, respectively, are already considerably larger than the N (r) = 45 049 available for I = 5 in the case of the AB tilings. With this increased statistical sample, we find that the agreement with P GOE (r) is now even better, particularly for I = 4. The final r = 0.530 347 value is indeed within ∼0.07% of r GOE .
In presenting the results shown in this work, we have been careful to only show level statistics computed for spectra consisting of irreducible blocks of the Hamiltonians. If we were to not separate these irreducible sectors (according to phase and parity) we would of course get a P(r) that becomes progressively closer to P Poisson (r) just as is the case for P(s) statistics [9]. Surmises for such spectra are only known for P(s) [22] and not yet for P(r) [11], but reliable estimates for r exist [23]. In good agreement with these latter results, we find r = 0.419 71(4) for I = 5, when combining parities ±1 for sector 0 (two irreducible blocks) as well as r = 0.391 728(7) when using all seven sectors of the AA tiling (cf. Table I). The results from Ref. [23] give r = 0.423 415 and 0.391 048, respectively, for these two cases when weighting according to N (E 0).
Quasicrystals represent a material class between periodic crystals and aperiodic solids. As such, it has earlier been speculated that they might possess nonstandard level statistics [24,25]. However, our results allow us to conclude that both P(s) [9] and P(r) statistics for two-dimensional, QP tightbinding models are, within the numerical accuracy currently achievable, very well described by the GOE ensemble. For P(s), this holds after unfolding [8] such that even the small difference between P GOE (s) and P Wigner (s) is resolved. For P(r), as we show here, also the unfolding procedure becomes superfluous to reach the same conclusion.
The data accompanying this publication are available in [26].
We are grateful to Fabien Alet for pointing us to Ref. [23] and providing r values for combinations of irreducible blocks. We thank Warwick's Scientific Computing Research Technology Platform for computing time and support. U.G. gratefully acknowledges support from EPSRC through Grant No. EP/S010335/1.

APPENDIX
To compute an improved approximation to P GOE (r), we follow Ref. [11] and perform the analogous calculation for the (5×5)-matrix case. The joint probability distribution for the GOE ensemble for the (5×5)-matrix case is [7] (e 1 , . . . , e 5 ) = C 5 1 i< j 5 where C 5 is the normalization constant. The distribution P (5) GOE (r) for the (5×5)-matrix case can then be computed as where we considered the eigenvalues to be ordered, with e 1 e 2 e 3 e 4 e 5 , and concentrated on the spacing around the central eigenvalue e 3 , which we believe to provide a better approximation than taking the average over the three terms arising from the spacings around e 2 , e 3 , and e 4 . When computing r = 1 0 rP (5) GOE (r)dr = 0.532 592, we see that the result is even closer to the high-precision numerical estimate r GOE ∼ 0.5307(1) than the numerical corrections using δP(r) (0.524 912) as proposed in Ref. [11]. However, evaluating the five integrals results in a lengthy expression for P (5) GOE (r); for details of the computation and the result, we refer to a Mathematica notebook [26]. We note that a systematic study for increasing (N×N ) matrices has been done previously in Ref. [27].