The GOE 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 non-interacting 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 behaviour 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 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 For an uncorrelated Poisson spectrum, one has P Poisson = 2/(1 + r) 2 with mean r Poisson = 0.386 while for the GOE, (2) 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 paper, 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 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 non-interacting tight-binding hopping Hamiltonians in QP tilings, whether that result has been computed (i) by unfolding DOS [8] or IDOS [9], (ii) by restricting the analysis to regions in the spectrum that  I .In both panels, the dashed (green) lines gives PPoisson(r) while the solid (blue) lines represent the improved PGOE(r) and the dotted (black) line is (2).
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 inflation as for I → ∞, 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, 27137, and 157369 vertices, respectively.For the random-AB tiling, we use I = 1, 2, 3 and 4 inflation steps with 82, 478, 2786 and 16238 vertices, respectively.On these tilings, we define H = i,j |i j| as 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 neighboring vertices connected by an edge of unit length are denoted as i, j .
The AB tiling in 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.As 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.
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 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 behaviour of P GOE (r).Similarly, the bulk behaviour 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.526654 (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) [28].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 behaviour 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) = 699579, 4132302 and 2421369 for I = 2, 3 and 4, respectively, are already considerable larger than the N (r) = 45049 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.530347 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, 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 non-zero energy levels, while N (r) counts to 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 rp(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 RMSDmax.The averages for AB tilings with combined statistics for I = 5 and for the random-AB tiling with I = 4 are highlighted in bold.
we find r = 0.41971(4) for I = 5, when combining parities ±1 for sector 0 (2 irreducible blocks) as well as r = 0.391728 (7) when using all seven sectors of the AA tiling (cf.Table I).The results from Ref. [23] give r = 0.423415 and 0.391048, 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 posses non-standard level statistics [24,25].However, our results allow us to conclude that both P (s) [9] and P (r) statistics for twodimensional, QP tight-binding 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.

FIG. 1 .
FIG.1.Plot of P (r) for (a) the combined P (r) (red dots) for all seven sectors of the AB tiling with inflation step I = 5 corresponding to N = 157369 vertices and (b) the combined P (r) (red dots) for the 300 (combined) realizations of the largest random-AB tiling at I = 4 as detailed in TableI.In both panels, the dashed (green) lines gives PPoisson(r) while the solid (blue) lines represent the improved PGOE(r) and the dotted (black) line is (2).

FIG. 2 .
FIG. 2. Schematics of (a) the octagonal AB tiling with N = 833 vertices (red dots) and exact D8 symmetry around the central vertex.Background shades/colors indicate successive inflation steps (I = 0, 1, 2) of the central (orange/dark) octagon.Details of the inflation procedure are given by the two small figures which show how to inflate each rhombus and square.Arrows along the diagonals of the squares provide directions that are required to define the inflation rule, since the dissection of the square breaks its symmetry.The results tiling, when starting from a symmetric seed like in part (a), keeps the D8 symmetry for each I.(b) A perfect periodic approximant of the AB tiling with N = 478 and (c) a random-AB tiling, also of size N = 478.The small figure on the right shows an individual "singelton flip" which induces the randomness.All tilings (a)-(c) correspond to I = 2.The fundamental regions for (b) and (c) are enclosed within the blue lines, forming a square.The green dots on the blue line denote the periodically replicated vertices for periodic boundary conditions.The thin solid lines connecting vertices i, j indicate the neighbor connection i, j as used in the Hamiltonian.

TABLE I .
Table of r-value estimates and RMSD for (top) AB and (bottom) random-AB tilings for different inflation levels I.