Quantum phase transition between hyperuniform density distributions

We study an electron distribution under a quasiperiodic potential in light of hyperuniformity, aiming to establish a classification and analysis method for aperiodic but orderly density distributions realized in, e.g., quasicrystals. Using the Aubry-Andre-Harper model, we first reveal that the electron-charge distribution changes its character as the increased quasiperiodic potential alters the eigenstates from extended to localized ones. While these changes of the charge distribution are characterized by neither multifractality nor translational-symmetry breaking, they are characterized by hyperuniformity class and its order metric. We find a nontrivial relationship between the density of states at the Fermi level, a charge-distribution histogram, and the hyperuniformity class. The change to a different hyperuniformity class occurs as a first-order phase transition except for an electron-hole symmetric point, where the transition is of the third order. Moreover, we generalize the hyperuniformity order metric to a function, to capture more detailed features of the density distribution, in some analogy with a generalization of the fractal dimension to a multifractal one.

In Ref. 13, we showed that the electron-charge distribution on the Penrose tiling, as well as of the AAH model, is characterized by hyperuniformity. Hyperuniformity, coined by Torquato and his collaborators [29,30], is a framework to quantify the regularity of the spatial distribution of a point set and has been generalized to a random scalar field [30][31][32]. It measures a density fluctuation of a given point set or scalar field distributed in a d-dimensional space and distinguishes different distributions according to the strength of the density fluctuation at a large length scale. Periodic and quasiperiodic point sets (i.e., lattice) are known to be hyperuniform. Namely, they possess significantly small density fluctuations thanks to the regularity of the lattices.
Various quasiperiodic lattices (as point sets) have then been classified in terms of hyperuniformity classes and its order metric [30,33,34], which quantify the degree of regularity of a hyperuniform distribution. The relevance of the order metric to a band-gap size of photonic quasicrystals has also been suggested [35]. In contrast, the nature of hyperuniform electron states or distributions (as scalar fields) realized on quasiperiodic structures remains largely unexplored. In particular, unlike periodic systems, where the change of the charge distribution occurs as a phase transition accompanied by the translational-symmetry breaking, it is unclear if such a change on quasiperiodic lattices occurs as a phase transition since the translational symmetry is absent in the first place.
In this paper, we scrutinize the AAH model [10,11,36], which is a prototypical quasiperiodic model in one dimension and has been realized experimentally in ultracold atoms [37] and photonic quasicrystals [38], in light of the hyperuniformity. Because the AAH model exhibits extended, critical, and localized eigenstates according to the strength of the quasiperiodic potential [39], it has long been studied in the context of quantum localization. Here, we utilize this property to study the relationship between the electron-localization strength and hyperuniformity, focusing on the charge distributions rather than the eigenstates.
We find that the charge distribution in the AAH model is always hyperuniform but its class and the order metric change according to the quasiperiodic potential and the Fermi level. For a weak potential, where the eigenstates are extended, the charge distribution has no jump in its histogram, exhibiting Class-I hyperuniform behavior. At and above the self-dual point [10,36], where the eigenstates are critical and localized, respectively, the charge distribution has no jump and Class-I hyperuniform only when the DOS at the Fermi level vanishes; otherwise, it has a vanishing point or a jump in the histogram and belongs to Class II. We thus reveal a nontrivial relationship between the DOS, charge distribution, and hyperuniformity class. We then clarify that the change of the hyperuniformity class is the first-order phase transition except for the electron-hole symmetric point where it is of the third order. These results, in turn, uncover a significant difference of the AAH model from random systems, where the localized states do not constitute a hyperuniform charge distribution.
Furthermore, we generalize the order metric to a function for a hyperuniform scalar-field distribution, in some analogy with the generalization of the fractal dimension to the multifractal one [9]. This generalization allows us to quantify more detailed features of density distributions. This "multi-hyperuniformity" [40] would be useful to characterize various density distributions, which are not multifractal but hyperuniform.
The rest of the paper is organized as follows. In Sec. II, we introduce the AAH Hamiltonian and the method to calculate its electron distribution and hyperuniformity. In Secs. III A and III B, we show the results of the DOS and charge distribution for various strengths of the quasiperiodic potential, revealing a relation between them. In Sec. III C, we discuss the results of hyperuniformity for the charge distributions and find that its class changes with the DOS at the Fermi energy as well as the continuity of the charge-distribution histogram. In Sec. III D, we reveal that the abrupt change of the hyperuniformty is indeed a phase transition. In Sec. III E, we introduce the "multi-hyperuniformity" to characterize more details of the density distribution. Section IV summarizes the paper. In Appendix A, we compare the results with those obtained for the Fibonacci models, to find a similarity to the critical case of the AAH model. Appendix B is devoted to demonstrate that the charge distribution under a random potential is not hyperuniform. Appendix C shows the results for an integrated intensity function, which gives an alternative way to calculate the hyperuniformity class. In Appendix D, we demonstrate that Class-II hyperuniform distributions remain Class II in our definition of the "multi-hyperuniformity". Appendix E presents the results of a local variance.

A. Aubry-André-Harper model
The AAH Hamiltonian [10,11] reads whereĉ i (ĉ † i ) annihilates (creates) a spinless fermion (which we call electron in this paper) at a site i on a one-dimensional chain with the lattice constant a = 1. t represents the hopping integral to the neighboring sites and λ does the strength of the quasiperiodic potential, where τ = √ 5+1 2 is the golden ratio. We set t = 1 and use it as the unit of energy. We have added the chemical potential (µ) term to the original Hamiltonian, to discuss the relation between the spectrum and the charge distribution. Considering zero temperature, we definen ≡ 1 N i n i with the number of sites N and where ψ α and E α are the eigenstates and eigenenergies of the Hamiltonian (1), respectively. This model is known to be self-dual at λ = 2t. Namely, the form of the Hamiltonian does not change after the Fourier transformation to momentum space with exchanging λ and 2t. As a consequence, the eigenfunctions are extended (localized) in real space for λ < 2t (λ > 2t) and critical at λ = 2t [see We numerically diagonalize a one-dimensional chain of N = F n sites, where F n is the n-th Fibonacci number. We use N = F 24 = 75025 with a periodic boundary condition, where τ in Eq. (1) is approximated by Fn Fn−1 . In this case, the phase shift φ does not play a significant role in the eigenvalues, unlike the topological surface states observed for the open-boundary condition [41]. We therefore set φ = 0 in the following. By comparing the results with those obtained with other sizes, we have confirmed that N = F 24 is sufficiently large to infer the infinite-size limit.

B. Hyperuniformity
Hyperuniformity is a framework to distinguish and quantify various spatial distributions. It was invented by Torquato and Stillinger [29] originally for point patterns distributed in space and has been generalized to various types of distribution including a random scalar field [30][31][32].
In one dimension, we consider a window of a range [−R, R] and count the number of points (or the total value of the scalar field) contained in the window. Namely, denoting the center position of the window as r c , we calculate the quantity, with the Heaviside step function Θ(r). Then, its variance is given by where Q represents the average of Q with respect to the center position r c over the system. While σ 2 (R) is proportional to R d (with d = 1 in the present case) for a random distribution of n i , the distribution with σ 2 (R) < O(R d ) is called hyperuniform, which means that a bulk contribution to the variance vanishes. Hyperuniform distributions are further classified into several classes: In one dimension, a distribution is called Class-I and Class-II hyperuniform, respectively, when the large-R behavior of σ 2 (R) is constant and proportional to log R [30]. Point distributions (i.e., n i ≡ 1) on periodic and quasiperiodic lattices are known to be hyperuniform [29,30].
To judge if a one-dimensional distribution is hyperuniform from a finite-size calculation, we consider the following function, If A(R) goes to zero as R increases, the distribution is hyperuniform. In particular, when it is Class I, i.e., σ 2 (R) = const. for a large R, we define [30] B(R) ≡ 1 Namely, we average over [0, R] to infer the order metric B(∞) because σ 2 (R) typically oscillates with R around its mean value. The factor 1/n 2 is just to eliminate a trivial contribution fromn to σ 2 (R).

A. Density of states and eigenfunctions
We first review known results for the density of states (DOS) and the eigenfunctions of the AAH model, showing calculated results. As shown in Figs. 1(a) and 1(d), when the quasiperiodic potential is weak (λ < 2), the electron state is extended in real space and the DOS has a continuous spectrum (though it is separated by gaps). At λ = 2, eigenfunctions exhibit self-similar distributions and the DOS is singular continuous [Figs. 1(b) and 1(e)]. For λ > 2, eigenfunctions are localized and the DOS is a dense set of δ functions. Note that the maximum (minimum) eigenvalue of the Hamiltonian (1) at µ = 0 is E max (−E max ) with E max = 2.1441, 2.5975, and 3.3862 for λ = 1, 2, and 3, respectively. We hence vary µ only within [−E max , E max ] in the following.

B. Charge distribution
We find that the charge distribution {n i }, which is a sum over eigenstates below the Fermi energy [i.e., Eq. (2)], also changes its character with λ. Figures 1(g-i) show the results for µ = 0. At λ = 1, n i continuously distributes from its minimum to the maximum [ Fig. 1(g)]. At λ = 2, however, the population of n i decreases around the center of the distribution [ Fig. 1(h)]. At λ = 3, the distribution bifurcates into roughly two values and shows a gap between them [ Fig. 1 Note that these (and following) results at λ = 1 (λ = 3) are representative of the results for λ < 2 (λ > 2) as we have obtained essentially the same results for various values of λ < 2 (λ > 2) though not explicitly shown.
As µ increases, the average fillingn increases unless it is in the gap of the DOS, wheren does not change [ Figs decreases and vanishes at a single point [ Fig. 2(e) middle and right panels], so that ∆ max vanishes. However, for λ = 3, when µ is located at a support of the DOS, the distribution has a jump [ Fig. 2(f) middle and right panels], so that ∆ max is finite as plotted in Fig. 3(b). Remarkably, the support of ∆ max completely agrees with that of the DOS shown in Fig. 1(c). Namely, we find a nontrivial relationship between the spectrum and the charge-distribution histogram.

C. Hyperuniformity
We analyze these charge distributions in terms of hyperuniformity. We first plot in Fig. 4 A(R) of Eq. (5) for various values of λ and µ. Irrespective of the potential strength λ and whether µ is in a gap of the DOS, A(R) always decreases in a power law and goes to zero in the large-R limit. Therefore, the charge distribution of the AAH model is always hyperuni- form. We point out here that this fact discriminates the AAH model from random systems, where the charge distribution is not hyperuniform (see Appendix B).
We then calculateB(R) of Eq. (6) for the same parameters and plot them in Fig. 5. In panels (a-d)B(R) converges to a constant value at a large R, which means that these distributions are Class-I hyperuniform [30]. As we have seen in Fig. 2, all these distributions have a histogram without a jump.
On the other hand, in Figs. 5(e) and 5(f),B(R) increases logarithmically with R, which means that these distributions are Class-II hyperuniform [30]. We have confirmed this point with another calculation in momentum space, too (see Appendix C). As we have seen in Figs. 2 and 3, these distributions have a point or a region where the fraction in the histogram becomes zero: For λ = 2 and µ = 0, the fraction vanishes at n c = 0.5 while for λ = 3 and µ = 0 the distribution has a finite jump in the histogram.
To quantify the above argument, we calculate the difference ofB(R)'s calculated at R = 1000 and 2000. For λ < 2, this quantity is virtually zero while for λ ≥ 2 it can be finite depending on µ. Figures 6(a) and 6(b) show the results at λ = 2 and 3, respectively. For λ = 3, the µ values giving a finite difference ofB(R) completely agree with the µ values of a finite ∆ max in Fig. 3(b), as well as with the support of the DOS in Fig. 1(c). For λ = 2, corresponding to the singular continuous spectrum in Fig. 1(b), the difference ofB(R) shows peaks of measure zero at the support of the DOS.
Thus, we obtain i) Class-I hyperuniformity for λ < 2, ii) Class-I or II hyperuniformity for λ ≥ 2, depending on the location of µ in the DOS. For the Class-I hyperuniformity, the order metricB =B(∞) represents the degree of regularity. Generally speaking,B is smaller for a simpler distribution [30].
We evaluateB at R = 2000 and plot it for λ = 1, 2, and 3 in Fig. 7. For λ = 1,B is always defined and relatively small. Since n i ≡ 1 for µ ≥ E max ,B at µ = E max agrees with that of the point distribution of the integer lattice, 1/6 [30]. On the other hand, as µ approaches −E max ,B goes to ∼ 0.31, which is consistent with the value obtained for the lowestenergy eigenstate in Ref. [13].B changes significantly when µ moves within the support of the DOS, while it is constant for µ moving within a gap. As we see in Fig. 1(a), the DOS has a sharp (δ-functional) peak at the gap edges. A general trend is that the inclusion of the states around the upper edge of a gap reducesB significantly while the states around the lower edge of a gap increaseB relatively less significantly.
For λ = 2, the distribution is Class-II hyperuniform when µ is located on the support of the singular continuous DOS. However, as µ crosses it,B changes significantly, whileB is constant for µ inside a gap. For µ = E max ,B = 1/6 for the same reason as above. On the other hand,B becomes extremely large for µ → −E max . This is because the lowestenergy eigenstate at λ = 2 is multifractal and is not hyperuniform. Although other eigenstates are also multifractal, the charge distribution, which is a sum over many eigenstates below the Fermi level, is hyperuniform. As µ increases, more states contribute to n i , makingB tend to decrease.
For λ = 3, the region of Class-II hyperuniformity expands, corresponding to the DOS in Fig. 1(c). In the Class-II region, the histogram of n i has a jump, as we have seen in Figs. 2(f) and 3(b). In other regions, the histogram has no jump and B is well defined. It tends to decrease as µ increases across the support of the DOS, except for the region around µ = 0. WhileB = 1/6 at µ = E max , the region slightly above µ = −E max is Class-II hyperuniform. However, in the limit of µ → −E max , it is not hyperuniform [13] since the lowestenergy eigenstate is localized.
Performing similar calculations for various values of µ and λ, we summarize the results ofB and hyperuniformity class in Fig. 8. The black region shows Class II while the color in other regions represents the order metricB of Class I. A general trend is thatB is larger for a smaller µ and larger λ.
In Appendix A, we show that the Fibonacci models, where the eigenstates are always critical, show a behavior similar to the λ = 2 case of the AAH model.

D. Phase transitions and criticality
As we have found in Fig. 8, the hyperuniformity class and order metric change with µ and λ. In particular, abrupt changes occur at the border of the Class-I and II regions. In this section, we examine whether these changes manifest themselves as a phase transition. We numerically calculate the total energy, and its derivative (evaluated by a difference between neighboring two data points) with respect to λ for fixed µ's at zero temperature. Since the distribution of {E α } is electron-hole symmetric, we concentrate only on the µ ≥ 0 side. First, for λ < 2, the order metric in Fig. 7(a) shows abrupt changes when µ crosses the gap edge. Since the eigenstates are extended, this is a metal-insulator transition. As a function of λ, too, E tot shows a kink and its first derivative shows a jump, as shown in Figs. 9(a) and 9(b). Here, we have chosen several µ values which show a singularity around λ = 1.5. The first derivative shows a rapid increase around the critical point presumably because of the large DOS at the gap edges in one dimension. Aside from these singularities, E tot and dE tot /dλ curves are smooth, showing no phase transition, even thoughB changes.
For λ > 2, on the other hand, all the eigenstates are localized, so that no metal-insulator transition occurs. Nevertheless, E tot plotted against λ still shows a kink, as shown in Fig. 9(e), where we have chosen several µ values crossing the border of Class-I and -II regions in Fig. 8. Notice that the relatively flat side corresponds to Class I. The presence of the kink is evidenced in the plots of dE tot /dλ in Fig. 9(f). This transition may be viewed as a transition from a band insulator (in the sense that the DOS vanishes though a 'band' is not well defined) to an Anderson insulator (where the DOS is finite but the mobility vanishes though the potential is not random but quasiperiodic).
Around λ = 2, we need a more careful analysis because the DOS is singular continuous. We have fine-tuned the µ values to several eigenenergies at λ = 2 and plotted E tot and dE tot /dλ in Figs. 9(c) and 9(d), respectively. We find kinks in E tot and jumps in dE tot /dλ at λ ∼ 2 for all the µ values except µ = 0. We see several additional kinks for λ 2, which are due to the crossing of the eigenenergies with a very small measure in this region.
All the above results except for µ = 0 show the first-order transition between the gapped and ungapped regions. On the other hand, at µ = 0, where the hyperuniformity class changes at λ = 2, no jump is observed in dE tot /dλ [ Fig. 10(a)]. In fact, µ = 0 is special because the DOS does never vanish for any λ due to the electron-hole symmetry and the self-duality. We then calculate d 2 E tot /dλ 2 and d 3 E tot /dλ 3 (by a difference between neighboring two data points), plotting them in Fig. 10(b). We find that d 2 E tot /dλ 2 is still continuous but has a kink whereas d 3 E tot /dλ 3 shows a jump. This weak singularity may be attributed to the singular continuous DOS at λ = 2. We have thus revealed a third-order criticality at λ = 2 and µ = 0.
These results clarify whether and where a phase transition occurs between electronic states with different inhomogeneous but orderly charge patterns. For λ < 2, while the observed phase transition is attributed to the metal-insulator one and is not so surprising, an important observation here is the absence of the phase transition in other regions whereB (and hence the charge distribution) smoothly changes. For λ > 2, the phase transition occurs between two different insulating phases characterized by different hyperuniformity classes; no phase transition occurs within the same hyperuniformity class. These results in turn prove an essential role of hyperuniformity analysis, which allows us to detect the phase transition in aperiodic systems independently of the total-energy calculation, like the role played by the order parameter in periodic systems. E. Multihyperuniformity

Straightforward extension
So far the Class-I hyperuniform distributions have been characterized by just one scalarB. Here, with a simple extension of the definition (3) of N (R), we generalize the order metric to a function that should capture more detailed information on the density distribution. Namely, we define and then σ 2 q (R) in the same way as Eq. (4). In analogy with the multifractal dimension [9], the exponent q works as a filter to emphasize the contribution from a large (small) n i for q >  to zero as R increases, i.e., {n q i } is also hyperuniform. Then, the latter shows thatB q (R) converges to finite values for all q's in the large-R limit. Namely, {n q i } is Class-I hyperuniform for all q's. We have obtained the same conclusion for other values of µ as far as {n i } belongs to Class I, as one may infer from the moderate values ofB q in Fig. 13 below. Note that when {n i } is Class-II hyperuniform, {n q i }(q = 0) remains Class II for all the parameters we studied (Appendix D).

0(< 0). Corresponding to Eqs. (5) and (6), we define
In Figs. 13(a-c), we plotB q (measured at R = 1000) against q for various λ and µ. We find thatB q takes the minimum at q = 0, whereB q agrees with the value (1/6) for a point distribution, and is convex downward around q = 0. As |q| increases,B q monotonically increases on each side of q > 0 and q < 0. This reflects the larger spatial fluctuation for a larger |q|.
At λ = 1,B q is larger on the q < 0 (q > 0) side for µ > 0 (µ < 0). This is reasonable because for µ > 0 (µ < 0) small (large) values of n q i can be further away fromn q (and hence more irregular) and q < 0 (q > 0) emphasizes these contributions. For q > 0,B q tends to decrease as µ increases, as is expected from the behavior at q = 1 displayed in Fig. 7(a); at µ = 2, all the sites are almost completely filled, so thatB q is nearly flat for q > 0. For q < 0, on the other hand,B q shows a complicated dependence on µ though it should approach 1/6 eventually for µ → E max . In particular, the largē B q for µ = 1 is interesting because this means that the charge distribution is significantly inhomogeneous even for this relatively large value of µ. In fact, as we see in Fig. 22(a) in Appendix E, the fluctuation measured by the local variance is maximized around µ = 1.
As λ increases,B q tends to increase, reflecting the larger fluctuation and consequent irregularity, in particular on the q < 0 side. On the q > 0 side,B q tends to decrease with µ in accord with Figs. 7(b) and 7(c) for q = 1.B q shows a more complicated dependence on µ on the q < 0 side. It is interesting thatB q at λ = 3 is always larger for q < 0 than for q > 0. For µ < 0, this is opposite to what we have seen at λ = 1. This is presumably because {n i } for µ < 0 reflects more directly the structure of localized eigenfunctions, which have vanishingly small amplitudes at most sites.

Symmetric definition
In Fig. 13(a), we see thatB q at µ = 0 (black curve) is asymmetric with respect to q = 0. However, as the charge distribution at µ = 0 is symmetric with respect to n c = 0.5 [see Figs. 1(g) and 2(d)], it may be preferable to define an order metric to reflect this symmetry. The asymmetry ofB q defined by Eq. (10) comes from the fact that (0.5 + δ) q does not agree with (0.5 − δ) −q , where δ represents a deviation from the average value 0.5. Hence, to remedy this asymmetry, we define s i ≡ n i /(1 − n i ) and Notice that s i at n i = 0.5 + δ equals s −1 i at n i = 0.5 − δ. Then, we define σ sym q 2 (R) in the same way as Eq. (4) and agrees with the order metric of the point distribution andB sym q (R) is symmetric with respect to the transformation (µ, q) ↔ (−µ, −q) as far as the DOS for µ = 0 is symmetric with respect to ω = 0.
As was done above, we first check the large-R behavior of A sym q (R) in Figs. 11(d-f). The results show that {s q i } is hyperuniform for all the q values studied. We then plot in Figs. 12(d-f) the correspondingB sym q (R) against R. We find that {s q i } belongs to Class I for all the parameters for which {n i } belongs to Class I. We have obtained the same conclusion for all other choices of µ that we study though not shown. Note that, when {n i } belongs to Class II, {s q i }(q = 0) also shows Class-II behavior (Appendix D).
We plotB sym q measured at R = 1000 in Figs. 13(d-f). First, for λ = 1 and µ = 0 (black curve), we see that the curve is symmetric with respect to q ↔ −q, as expected.B sym q takes the minimum of 1/6 at q = 0. Second, all the curves are symmetric against the simultaneous sign reversal of µ and q, i.e., (µ, q) ↔ (−µ, −q). Therefore, the asymmetry of thē B sym q curves for µ = 0 correctly represents the asymmetric distribution of {n i } aroundn.
Another notable difference from theB q curves is that theB sym q curves do not approach a flat curve for µ → E max (see blue curves). This is due to the denominator of n i /(1 − n i ), which amplifies more the sites closer to n i = 1. Namely,B sym q for µ → E max reflects the structure of the highest-energy eigenfunction, just asB q for µ → −E max does for the lowest-energy eigenfunction. Note thatB sym q for µ → −E max still reflects the structure of the lowest-energy eigenfunction though its contribution toB sym q differs from that toB q due to the difference between n i [in Eq. At λ = 1,B sym q is larger for q < 0 (q > 0) for µ > 0 (µ < 0) for the same reason described above forB q . The same occurs for λ = 2 and even for λ = 3 and µ = ±1. For λ = 3 and µ = ±2.8, while the same occurs for |q| 3, it is reversed for |q| 3. This is likely because the structure of the highest-or lowest-energy eigenstates (rather than the filling controlled by µ) becomes more relevant for µ close to ±E max as mentioned above.
In Appendix A, we calculateB sym q for the Fibonacci models. The convex-down behavior around q = 0 and a monotonic increase with |q|, as well as a large enhancement at µ's close to ±E max , are common to the Fibonacci models.

Application to the critical regions
Here, we apply the multihyperuniformity analysis to a critical behavior around the phase transition discussed in Sec. III D. Our aim is to clarify how the inhomogeneous charge distribution changes around the critical point, by quantifying it through the generalized order metric.
In Fig. 14(a), we focus on the continuous transition point at λ = 2 and µ = 0. Since the order metric is defined only in the Class-I hyperuniform region, we calculateB sym q only for λ < 2. We find a rapid increase ofB sym q for large |q|'s as λ approaches the critical point. This behavior means an increasing irregularity of the sites with a particularly large or small electron density. Notice thatB [Eq. (6)] alone cannot distinguish such a behavior from an overall increase of irregularity. In the inset, we plotB sym q against 2 − λ in a logarithmic scale for several values of q. For each q,B sym q increases in a power law as λ approaches the critical point, 2. The power seems to weakly depend on q, e.g., −0.166 at q = 1 and −0.248 at q = 4 for 2 − λ < 0.005.
Bifurcated by a jump II By contrast, Fig. 14(b) shows thatB sym q does not change on the Class-I side of the first-order phase transition at λ ≃ 2.866 and µ = 2.5. The three curves are almost completely overlapping here. This of course means no significant change in the charge distribution up to the transition point and a jump there.
Our generalization thus offers a useful tool to analyze inhomogeneous density distributions, which are not multifractal but hyperuniform, and their changes by quantifying the irregularity of each contribution.

IV. SUMMARY AND PERSPECTIVES
We have studied the charge distribution in the Aubry-André-Harper model in light of hyperuniformity. According to the strength λ of the quasiperiodic potential, the model is known to exhibit extended, critical, and localized electron states. In this paper, we have revealed that the inhomogeneous distribution of electron charge n i , which is neither periodic nor multifractal but still orderly, also changes its character with λ. The character is quantified in the framework of hy- peruniformity generalized to density distributions.
First, we have found a nontrivial relationship between λ, the DOS at the Fermi level, {n i }, and the hyperuniformity class, as summarized in Table I. For λ < 2, where eigenstates are extended, the charge distribution has no jump and is Class-I hyperuniform. There is no phase transition as far as the order metric changes smoothly while a first-order metal-insulator transition occurs in concomitance with an abrupt change of the order metric when the Fermi level µ crosses the gap edge. For λ > 2, where eigenstates are localized, the charge distribution has no jump and Class-I hyperuniform only when µ resides in the gap of the DOS; otherwise, the charge distribution has a jump in its histogram and belongs to Class II. While all the electron states are insulating in this region, the change from Class I to II manifests itself as a first-order phase transition. At λ = 2, where eigenstates are critical, the charge distribution has no jump and Class-I hyperuniform only when µ is in the gap of the DOS; otherwise, it shows a behavior vanishing at a single point in the histogram and belongs to Class II. The transition is of the third order at µ = 0 and the first order otherwise. For the Class-I hyperuniform distributions, we have also revealed the dependence of the order metric on λ and µ.
The hyperuniform charge distributions for λ > 2 discriminate the AAH model from random systems, where the eigenstates are localized but the charge distribution is not hyperuniform (Appendix B). In addition to this, the eigenstates for λ < 2 are also hyperuniform in the AAH model [13]. These facts may make a significant difference between the localization-delocalization transition at λ = 2 in the AAH model and the Anderson transition discussed in random systems in higher dimensions.
Since various extensions [42][43][44][45][46][47][48][49][50][51] have been proposed for the AAH model, it is intriguing to explore these models in light of hyperuniformity. Of particular interest is the coexistence of localized and extended states at the same quasiperiodic potential observed in several models preserving a selfduality. The hyperuniformity analysis of the charge distribution in these models constitutes an important future issue.
Although the order metric seems to represent well a regularity of the aperiodic density distributions, it is obvious that much information about the distribution is lost in this quantification. We therefore extend the order metric to a function, in analogy with the extension of the fractal dimension to the multifractal one [9]. In both the straightforward extension and a symmetric definition, we first confirm that the order-metric function is well defined, i.e., {n q i } and {s q i } belong to Class I when {n i } belongs to Class I. Thanks to the filtering effect of the power q, the order-metric function,B q orB sym q , represents the regularity of differently weighted subsets of {n i }. In particular,B sym q can correctly capture the asymmetry of the distribution.
This generalization applies to any density distribution ranging from 0 to 1 (i.e., probability distribution). As mentioned in the introduction, there are various density distributions, which are known to be neither random nor multifractal, on quasicrystalline structures. Some of them may be hyperuniform. For instance, when an electron property on a quasiperiodic lattice is determined by short-range physics, it is likely hyperuniform. To analyze such distributions, the generalized ordermetric function will be a useful tool. when the chemical potential resides in a gap of the DOS, and Class-II hyperuniform otherwise.
We examine the above expectation for the following two types of the Fibonacci model.
where V i = +V or −V according to the Fibonacci sequence. Off-diagonal model: where t i = t L or t S according to the Fibonacci sequence. We numerically diagonalize the Hamiltonian for N = F 24 = 75025 sites under periodic boundary conditions. For µ = 0, these models show the DOS of Figs. 15(a) and 15(b), respectively. We see that the DOS at the Fermi level (ω = 0) is zero for µ = 0 in the diagonal model with V = t = 1 and for µ = −0.5 in the off-diagonal model with t S = 2t L = 1. On the other hand, µ = −1.5 in the diagonal model and µ = −1 in the off-diagonal model are very close to the support of the DOS, whose measure is zero.
After confirming that A(R) of Eq. (5) goes to zero in the large-R limit, we plot in Fig. 16B of Eq. (6) against R. We show Class-II hyperuniformity. Note that a possible deviation from the expected log R behavior at large R is attributed to the slight deviation of µ from the support of the DOS. These results are fully consistent with those obtained for the AAH model at λ = 2.
A recent study [52] of the Fibonacci model revealed that the charge-density oscillation in the perpendicular space is related to the topological property when µ resides in a gap. Its relation with the Class-I hyperuniformity in the physical space is an interesting subject of future research.
In Fig. 17, we plotB sym q of Eq. (13) for the (a) diagonal and (b) off-diagonal models, where we select µ values residing eight major gaps in the DOS of Fig. 15. All the curves take the minimum at q = 0 and are convex downward around it, similarly to the results for the AAH model [ Figs. 13(d-f)]. We also see thatB sym q tends to be large for µ close to ±E max . For the diagonal model,B sym q shows a complicated dependence on µ. This would be at least partly due to the asymmetry in the DOS.
For the off-diagonal model, on the other hand,B sym q shows a symmetry with respect to the exchange of (µ, q) and (−µ, −q), due to the electron-hole symmetry of the DOS. Interestingly, for µ > 0(< 0),B sym q is larger on the q > 0(< 0) side, on the contrary to the behavior in the AAH model at λ = 2 [ Fig. 13(e)], suggesting a large irregularity of the eigenstates. In addition, the µ = 0.1 and 0.5 curves show a similar behavior to each other. This may be related to a self-similarity of the model, because the gaps around these two µ points are related by a self-similar transformation. Here, we demonstrate that the charge distribution in the localized phase in a random system is not hyperuniform. We consider the following one-dimensional Anderson model [53], where W i is a random potential independently and uniformly distributed in the range [− W 2 , W 2 ] (W > 0). All the states are localized for W = 0 [54,55]. We numerically diagonalize the above Hamiltonian for 50000 sites and calculate the charge density at each site based on Eq. (2). We then calculate A(R) of Eq. (5) for the charge distribution.
The results for W = 1 and 2 are plotted in Fig. 18. We see that A(R) remains finite at a large R. This means that the charge distribution of the model (B1) is not hyperuniform, unlike that of the AAH model.
The above results show that even in the localized (λ > 2) region of the AAH model, there is a significant difference from the random system in light of the hyperuniformity of the charge distribution: In the AAH model, it is either Class-I or II hyperuniform while it is not hyperuniform in a random system. This difference may be used to distinguish a localization in quasiperiodic systems from that in random systems experimentally.

Appendix C: Integrated intensity function
Here, we study the behavior of the structure factor, at the long-wavelength limit (k → 0). The asymptotic behavior, S(k) ∼ k α for k ∼ 0, is characterized by α > 1 for a Class-I and α = 1 for a Class-II hyperuniformity [33]. Because this classification based on α does not rely on a window used in Sec. II B, it gives an independent check for the determination of the hyperuniformity classes. For quasiperiodic systems, where S(k) consists of a dense set of Bragg peaks, an integrated intensity function, is smoother and hence more useful than S(k) [33]. Because Z(k) behaves as k α+1 for k ∼ 0, we plot it for (a) λ = 2 and µ = 0 and (b) λ = 3 and µ = 0 in a logarithmic scale in Fig. 19. We see that the results are consistent with α = 1 in both cases, supporting that the charge distributions for these parameters are Class-II hyperuniform. hyperuniformity, which characterizes the long-range density fluctuation.
Here, we study how this local variance changes with µ and λ. Figure 22 shows the results for λ = 1, 2 and 3. An overall trend is that the local variance is maximized around µ = 0 and decreases as µ approaches ±E max , as anticipated. However, for λ = 1, the local variance shows a dip around µ = 0, making a local minimum at µ = 0. While the local variance increases monotonically with µ < 0 for λ = 2, it shows a nonmonotonic dependence on µ < 0 for λ = 3. The difference between λ ≥ 2 and λ < 2 may be attributed to the FIG. 22. (a,b,c) The local variance for λ = 1, 2, and 3, respectively. The red lines denote the values of µ presented in Fig. 2. presence/absence of the jump in the n i histogram.