Moire Localization in Two Dimensional Quasi-Periodic Systems
MMoire Localization in Two Dimensional Quasi-Periodic Systems
Biao Huang ∗ and W. Vincent Liu
1, 2, † Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh PA 15260, USA Wilczek Quantum Center, School of Physics and Astronomy and T. D. Lee Institute,Shanghai Jiao Tong University, Shanghai 200240, China (Dated: October 17, 2019)We discuss a two-dimensional system under the perturbation of a Moire potential, which takes the samegeometry and lattice constant as the underlying lattices but mismatches up to relative rotation. Such a self-dualmodel belongs to the orthogonal class of a quasi-periodic system whose features have been evasive in previousstudies. We find that such systems enjoy the same scaling exponent as the one-dimensional Aubry-Andre model ν ≈
1, which saturates the Harris bound ν > / d = ff erent from typical one-dimensional situation where only a few or no step-like contours show up. An experimental scheme based onoptical lattices is discussed. It allows for using lasers of arbitrary wavelengths and therefore is more applicablethan the one-dimensional situations requiring laser wavelengths close to certain incommensurate ratios. I. INTRODUCTION
Overlaying lattices with mismatches leads to the so-calledMoire pattern, which has long been an intriguing subject intwo-dimensional heterostructures [1–3]. Usually, a commen-surate twist angle between two layers are considered, such thatthe bilayer system still possesses an enlarged periodicity. Theinterlayer hopping and interactions would then induce a re-construction of the original Bloch waves in each layer intomultiple bands. That provides a new knob to tune the prop-erty of systems with various twisting angles, leading to novelexamples including the flat-band induced superconductivityfor bilayer graphene twisted at magic angle [4, 5]. More re-cently, the experimental advancement has made it possible tostabilize the bilayer graphene system at incommensurate largetwist angle [6]. In such situations, the combined bilayer sys-tem breaks any translation symmetry and the Bloch waves ineach layer are destroyed completely. Then, it is of interestto ask what may be the general phenomena expected in thesecircumstances.The incommensurately twisted bilayer system resemblesquasi-crystals, in the sense that only rotation but not trans-lation symmetry is preserved therein. In the case of quasi-crystals, localization of particles has been a focus of studyfor decades [7–10]. One could relate the two systems usingthe following schematic reasoning. Consider a bilayer, two-dimensional system with interlayer interactions, H bi = (cid:88) (cid:104) i , j (cid:105) ( a † i a j + b † i b j ) + (cid:88) i , j V i j n Ai n Bj , (1)where a † i , b † i create particles in two layers respectively, and n Ai = a † i a i , n Bi = b † i b i . Now, suppose one could take asa starting point the decomposition µ Ai = (cid:80) j V i j (cid:104) n Bj (cid:105) , and ∗ [email protected] † [email protected] µ Bj = (cid:80) i V i j (cid:104) n Ai (cid:105) , then the above model reduces to H e ff = (cid:88) (cid:104) i , j (cid:105) a † i a j + (cid:88) i µ Ai n Ai + ( A ↔ B ) , (2)which describes a decoupled bilayer under onsite chemicalpotentials µ A , Bi respectively. With reasonable interactions,i.e. not infinite-range interaction with identical strength, onewould expect an incommensurate twist angle to result in anincommensurate, quasi-periodic pattern for µ A , Bi . That givesan analogous scenario as in quasi-crystals, with the di ff erencethat here, it is the chemical potentials rather than lattice sitelocations that break translation symmetry.The inter-sample interaction induced slow-relaxation(“quasi many-body localization”) has previously been demon-strated for a one-dimensional translation invariant ladder byYao et. al. [11]. Specifically, given a highly non-uniform ini-tial density distribution and strong interactions, the particlesin one chain would serve as random chemical potential forthe other, similar to the case mentioned above. But in con-trast, due to the incommensurate twist angles breaking trans-lation symmetry here, it is expected a genuine localizationwithout initial state dependence would occur for incommen-surate Moire bilayers. In particular, even if the particles areuniformly distributed in one layer, they could still function asdisorder potentials for the other through quasi-periodic V i j .Thus, we are motivated to investigate the “Moire localiza-tion” possibly held in a two-dimensional quasi-periodic sys-tem. As a first step, we focus on an e ff ective scenario inEq. (2), namely, a regular single-layer lattice under the per-turbation of quasi-periodic potentials. This would serve as agood starting point for further taking into account fluctuationsin Eq. (1). Also, although the localization of quasi-periodicsystems in one-dimension have been thoroughly studied fordecades [12–17], since the renowned work by Aubry and An-dre [18] in 1980, its generalization to higher dimensions hasjust started quite recently [19–21]. The theoretical studies onthis regard have remained chiefly in the single-particle level,with many crucial aspects still open to discussions. In partic-ular, the scaling exponents for orthogonal class of models in two-dimensions have been left out in Ref. [20, 21] due to var- a r X i v : . [ c ond - m a t . qu a n t - g a s ] O c t ious di ffi culties. But this class of models are most relevant forexperiments on genuine two-dimensional quasi-periodic sys-tem with real hopping, and also for the theoretical analysisbased on Eqs. (1) and (2). Therefore, it is useful to clarifythe relevant single-particle physics before considering the fullmany-body interactions.Apart from preparing for analysis of interaction-induced lo-calization, generalizing the paradigmatic Aubry-Andre modelto higher dimensions itself may yield interesting results. Asshown by Devakul and Huse in Ref. [20], for orthogonalclass of models in three dimensions, the scaling exponentsfor (single-particle) localization transitions are the same forboth quasi-periodic and purely random potentials. This is incontrast to the situation in one-dimension, where the quasi-periodic Aubry-Andre model gives rise to the scaling expo-nent strongly violating the Harris bound [22] formulated forpurely random models. Such a violation in one-dimensionis indicated in Ref. [14] to persist into the many-body local-ization scenario. Experimentally, a higher-dimensional quasi-periodic potential is actually more suitable for cold atom ex-periment. This is because the property of the system is mostsensitive to the relative rotation angle between the Moire per-turbing potential and the main lattice potential, rather than therelative lattice constants. That means one could use lasers ofany frequency available in the laboratory to engineer the de-sirable system. Finally, having a good understanding of thesingle particle physics for quasi-periodic systems would pavethe way for possible perturbative treatment [23] of many-bodylocalization therein in the future.In this work, we study a two-dimensional square lattice un-der the perturbation of Moire-type of potential, which takesexactly the same lattice geometry and constants as the under-lying lattice. A commensuration condition for the relative ro-tation angles between the two lattice potentials is provided.Such a condition also helps identifying “better” rotation an-gles that are further from a commensurate one, resembling the“better” ratio of lattice constants in a one-dimensional quasi-periodic system, that is usually taken to be the golden ratio.Two important features of such a system is revealed. First,for the inverse-participation-ratios in the eigenstate-disorderplane, in contrast to the one-dimensional situation where onlya few step-like jumps exist, our two-dimensional model ex-hibits a rapid and continuous change. Such a character indi-cates the unusual mobility edges in two dimensions. Second,the critical exponent for localization length is extracted fromthe multifractal analysis, with the value ν ≈ ν > / d = II. MODELS AND THE COMMENSURATIONCONDITIONS
Consider a square lattice under perturbations of Moire-typeof potentials with the same geometry as the underlying ma-jor one. There are three types of controlling parameters for such a perturbing potential. (a) Stretching: it could possessdi ff erent lattice constants. (b) Rotation: it may be rotated byan angle θ with respect to the underlying lattice. (c) Trans-lation: there may be a global relative translation betweenthese two lattices. Both (a) and (b) could generate an in-commensurate potential, but as in the experiment Ref. [19],without interactions, (a) alone would result in a separable V ( x , y ) = V ( x ) + V ( y ) such that the model is reduced to twoorthogonal lower-dimensional ones. Meanwhile, (c) does notchange the universal properties of the system in thermody-namic limits, but only a ff ects microscopic details for a finitesize system, i.e. the wave functions in the boundary. Thus, theglobal rotation angle of perturbing potential is the most impor-tant factor characterizing a genuine two-dimensional quasi-periodic system. It was also found in previous works that thesystem’s properties are most sensitive to rotation angles of theperturbing potentials [20, 21].Thus, we focus on the following model with two overlap-ping square lattices, where the stronger one gives rise to thetight-binding approximation, and the weaker one generatesthe Moire type of onsite chemical potentials, H = − (cid:88) (cid:104) i , j (cid:105) ( c † i c j + c † j c i ) + V d (cid:88) i µ i c † i c i ,µ i = sin (cid:2) π u i − ϕ (cid:3) + sin (cid:2) π v i − ϕ (cid:3) , u i = x i cos θ − y i sin θ, v i = x i sin θ + y i cos θ,ϕ = π ( a cos θ − b sin θ ) , ϕ = π ( a sin θ + b cos θ ) . (3)Here i = ( x i , y i ) denotes the square lattice sites, with x i , y i ∈ Z .The Moire superlattice potential µ i , compared with the mainlattice, is rotated by an angle θ at the origin ( x , y ) = (0 ,
0) andthen translated by ( x , y ) = ( a , b ). The stretching (a) wouldbe neglected in most parts below, because it does not leadto extra features and only complicates the analysis. For in-stance, either stretching the lattice constant by a factor √
2, ora rotation of π/ ff erent rotation angles θ , and average over various transla-tions ( a , b ) (or equivalently, ϕ , ϕ ) when extracting universalproperties from finite size numerical results. T T (a) Commensurate (b) Incommensurate FIG. 1. The Moire potential pattern µ i for (a) the commensurateangle θ = arccos(60 /
61) ( q = , p = θ = π/
5. Each pixel denotes one lattice site, and the colorindicates the strength of the disorder potential therein.
The first question to ask is under what condition would therotation angle θ result in a commensurate perturbing potential µ i . Similar to the analysis in honeycomb lattices [6, 24, 25],the commensuration condition is prescribed by solving a Dio-phantine equation specified by lattice geometries. For squarelattices, the solution is readily provided by the Pontryagin’striples [26] (see Appendix A for derivations) θ c = m π ± arccos q − p q + p , m ∈ Z (4)where q , p are positive, coprime integers. This is the com-mensurate condition for Moire rotation angles of rectangularlattices with rational aspect ratios γ . Specific to the squarelattice γ =
1, due to four-fold rotation symmetry and mirrorsymmetry along ˆ x ± ˆ y directions, one only needs to consider θ c = arccos q − p q + p ∈ (0 , π/
4) and θ (cid:48) c = π/ − θ c = arccos qpq + p ∈ (0 , π/ q > p . The resulting Moire pattern hasperiodic structures, with the new Moire-Bravais lattice vec-tors T , T given below. When p or q is an even number (theycannot both be even due to the coprime requirement), T = q ˆ x − p ˆ y , T = p ˆ x + q ˆ y , (5)and when ( q , p ) are both odd numbers, T = q + p x + q − p y , T = q − p x − q + p y . (6)In both cases, the two vectors have identical length | T | = | T | = (cid:112) q + p (when q or p is even) or (cid:112) ( q + p ) / q , p are both odd), and are orthogonal to each other T · T =
0. An example of commensurate Moire pattern isdrawn in Fig. 1(a), where q = , p =
1, and the Bravaisvectors are T = x + y , T = x − y . In contrast, an in-commensurate example is shown in Fig. 1(b) where only thefour-fold rotation, but not any translation symmetry, exists.Here, we are interested in localization induced by incom-mensurate potentials, so it is worthwhile to choose among dif-ferent incommensurate angles for further numerical analysis.Strictly speaking, any angle other than those given by Eq. (4),i.e. any θ = arccos( A ) with A (or A module π/
2) being ir-rational numbers, would render incommensurate systems be-longing to the same universal class. But when one is trying toextract information from numerical results in a finite-size sys-tem, there is a subtle question of “how incommensurate” thepotentials are. Mathematically, it means if one takes a seriesof commensurate angles { θ n , n = , , . . . } , whose n → ∞ gives an incommensurate angle θ ∞ , then a “more incommen-surate” angle would require each cos θ n being given by larger( q , p ) in Eq. (4). Equivalently put, since the length of Moire-Bravais vectors ∼ (cid:112) q + p , that means given a system size(and therefore an upper limit of ( q , p ), since for larger ( q , p ),the system cannot cover even one Moire period), “more in-commensurate” θ ∞ ’s are farther away from commensurate an-gles. This is in the same spirit as in one-dimension, where the“golden-ratio” ( √ − / π / π / π / π / π / FIG. 2. The density of commensurate angles (blue polar lines), for0 ≤ p ≤ q ≤
20. Due to mirror symmetry along x = y , the angles aresymmetric along θ = π/ To illustrate this point more clearly, we plot the commensu-rate angles for maximal ( q , p ) ≤
20 in Fig. 2. We readily notethat for all small angles θ → θ = π/ θ = π/ ff erence be-tween bilayer heterostructure and optical lattice systems. Forbilayer ones, there is a competition between interlayer attrac-tive van der Waals force and elastic force. The former onetends to align the two layers by distorting lattices so as tominimize interlayer distance, while the latter one would favorkeeping lattice shapes and therefore the interlayer twisting.For small twist angles θ →
0, the van der Waals force over-whelms and the system could always be taken as a commen-surate one [27], with the enlarged periodicity approximatelygiven by 1 /θ . Although such an issue does not exist in opti-cal lattices where lattice shapes are fixed, for generality of theresults in practice, we consider the large twist angle π/ β , there still exists the need to pick outtwist angles further away from “almost-commensurate” oneso as to minimize finite size e ff ects. The process there wouldconsists of first finding a sequence of rational lattice constantsclose to β , and then for each member in this sequence, findthe density of commensurate angles. We would not digress tosuch a situation in this work. III. MOIRE LOCALIZATION
For the one-dimensional Aubry-Andre model [18], its self-duality gives rise to the unique localization transition (lackof mobility edge) and the uniform localization length. Butif the incommensurate potential is given by relative twistsrather than stretch, the two-dimensional eigenstates cannotbe decomposed into orthogonal one-dimensional components.Then we would see that such a higher-dimensional general-ization breaks the unique localization length and transition, inaccordance with the three-dimensional results in Ref. [20]. Itcan be traced back to that the Thouless formula adopted toprove these two points [18], which we review in Appendix D,no longer holds for a higher-dimensional system. We wouldalso illustrate the mobility edge clearly in numerics later on.
A. What does self duality imply in the Moire model?
For later convenience, we rewrite Eq. (3) as H = − (cid:88) x , y ( c † x + , y + c † x − , y + c † x , y + + c † x , y − ) c x , y − V d (cid:88) x , y µ x , y c † x , y c x , y ,µ x , y = cos (cid:104) π u x , y + ϕ (cid:105) + cos (cid:104) π v x , y + ϕ (cid:105) u x , y = x cos θ − y sin θ, v x , y = x sin θ + y cos θϕ = π ( a cos θ − b sin θ ) , ϕ = π ( a sin θ + b cos θ )(7)which is equivalent to Eq. (3) up to a constant. That alsomakes the model particle-hole symmetric after averaging overthe phases ϕ , ϕ . The single-particle eigenfunction φ x , y = (cid:104) x , y | φ (cid:105) = (cid:104) | c x , y | φ (cid:105) satisfies − ( φ x + , y + φ x − , y + φ x , y + + φ x , y − ) − V d µ x , y φ x , y = E φ x , y , (8)where H | φ (cid:105) = E | φ (cid:105) and | (cid:105) is vacuum. Consider the Fouriertransformation φ x , y = e π i ( x ϕ + y ϕ ) ∞ (cid:88) m , n = −∞ ψ m , n e π i [ n ( u x , y + ϕ ) + m ( v x , y + ϕ ) ] , (9)where m , n ∈ Z . Then, for incommensurate rotations where( m , n ) sin θ or ( m , n ) cos θ are never integers and therefore dif-ferent Fourier modes do not couple, the wave functions ψ m , n satisfy the dual equation − µ m , n ψ m , n − V d ψ m + , n + ψ m − , n + ψ m , n + + ψ m , n − ) = E ψ m , n . (10)Compare Eqs. (8) and (10), we see that the self-duality maps V d ↔ V d , and E ↔ EV d . Then the invariant disorder strengthunder duality mapping can be extracted as V ( D ) d = . (11)The dual Fourier transformation Eq. (9) has the propertythat if | ψ m , n | is a localized wave function, in the sense that (cid:80) m , n | ψ m , n | and | ψ m , n || m or n →∞ is finite ( | ψ m , n | is exponentiallylocalized around certain ( m , n )), we would have a delocal-ized φ x , y as (cid:80) x , y | φ x , y | → ∞ [28]. The boundary condition (orthe normalization condition) for the wave function is that themaximal value of | ψ m , n | at the localized site ( m , n ) is takento be a constant. As such, the duality mapping implies a one-to-one correspondence between a localized eigenfunction atdisorder V d and energy E , and a delocalized one at 16 / V d and energy 4 E / V d . It does not specify which one is the localizedsolution.So far, the analysis is completely in parallel with that for theone-dimensional Aubry-Andre model. However, to prove alleigenstates are localized for V d > V ( D ) d = V d < V ( D ) d ), that is, a unique localization transition forall eigenstates, Ref. [18] invoked the Thouless formula [29].Such a formula relates the localization length ξ E at energy E with the density of states D ( ε ) as 1 /ξ E = (cid:82) ∞−∞ (ln | E − ε | ) D ( ε ) d ε .Together with the duality mapping for density of states, onecould show that all eigenstates are localized at the V d > V ( D ) d side with uniform localization length 1 /ξ = ln( V d /
4) (see Ap-pendix D). However, the Thouless formula, which requires thenearest-neighbor hopping in one-dimension (i.e. Eq. (D4) and(D5)), no longer holds in higher dimensions. Then it is ex-pected the mobility edge would generally exist for a higherdimensional model.In the following subsections, we perform exact diagonaliza-tion to further reveal the nature of the localization transitionsin our Moire model.
B. Level statistics
We first compute the energy level statistics for the ratio ofneighboring gaps [30, 31]. Arranging the energy levels { E m } in the order E m < E m + , and defining the gap between neigh-boring levels in a finite-size system as δ m = E m + − E m , theratio reads r m = min( δ m , δ m + ) / max( δ m , δ m + ). The distribu-tion function of r m ’s, P ( r ), would approach a Poisson limit P ( r ) = / (1 + r ) with (cid:104) r (cid:105) = (cid:82) dr ( rP ( r )) ≈ .
386 for a fullylocalized system with complete sets of integrals of motion. Incontrast, for delocalized systems, the neighboring energy levelrepulsions would lead to the Gaussian orthogonal ensemble(GOE) with (cid:104) r (cid:105) ≈ . V ( D ) d with (cid:104) r (cid:105) approaching the Gaussian ensembles. It is of in-terest to see whether such an intermediate phase, which wouldindicate the existence of mobility edge, persist to our two-dimensional situation. Also, we check the change of (cid:104) r (cid:105) asthe system size increases so as to indicate the behavior in thethermodynamic limit.The qualitative behavior of (cid:104) r (cid:105) as a function of V d can beexpected as follows. At V d =
0, the translationally-invariantlimit, momentums are conserved and E k = − k x + cos k y ).Thus, it can be regarded as a localized system in momentumspace. Deviating away from V d =
0, the momentum conserva-tion is immediately broken by the Moire perturbing potential,while real-space localization has not been established. Then,due to the lack of complete local integrals of motion, the levelrepulsion would lead to the GOE distribution for (cid:104) r (cid:105) . As V d increases, more and more eigenstates are localized, and whenthe full localization for all eigenstates occur, (cid:104) r (cid:105) would onceagain approach the Poisson limit with real-space positions asgood quantum numbers to denote localized states.The results for θ = π/ ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■■■■■■■■■■■■■■■■■■■■■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ GOEPoisson ● L = ■ L = ◆ L = ▲ L = ▼ L = V d 〈 r 〉 FIG. 3. The average value for neighboring level gap ratios. Here θ = π/
5. The numbers of phases ϕ , averaged over are 4000 / L = /
500 ( L = /
100 ( L = /
50 ( L = /
20 ( L = V d within / out-of the range [3 , reasoning given above. Here, we avoid the subtle regime withsmall V d (cid:28)
1. This is because at V d =
0, the system possessesmany symmetries, such as translation, reflection, and four-fold rotation symmetries. That renders the Hamiltonian block-diagonalized into di ff erent symmetry sectors and a meaning-ful level statistics requires the decomposition of Hilbert spaceinto these sectors [32–35]. Close to V d =
0, those symme-tries are only weakly broken and a severe finite-size e ff ect isexpected. We readily see from Fig. 3 that at small V d , (cid:104) r (cid:105) de-creases with irregular behaviors. On the other hand, the tran-sition near V ( D ) d = V d = V ( D ) d = (cid:104) r (cid:105) keeps approaching the GOE limit as system sizes L increase.Passing V d = (cid:104) r (cid:105) starts to drop. We note that for relativelysmaller system size (i.e. L = , , V d ≈ . L . However, for larger L ’s such a “crossing” appears to be replaced by a convergenceof curves for di ff erent L ’s. We have also verified several othertwisting angles and similar results emerge. Thus, the levelstatistics indicates the existence of mobility edge; only when V d (cid:39) . (cid:104) r (cid:105) ap-proaches the Poisson limit.The average value (cid:104) r (cid:105) gives an averaged behavior for alleigenstates. We next look at the eigenstate-resolved inverseparticipation ratio for more details regarding the mobilityedge. C. Inverse Participation Ratio and Mobility Edge
The inverse participation ratio (IPR) is defined asIPR m = L (cid:88) x , y = | φ ( m ) x , y | , (12)for a normalized eigenstate m : (cid:80) x , y | φ ( m ) x , y | =
1. The IPRcrosses from 0 in the strongly delocalized (i.e. take | φ r | = / L )regime to 1 in the localized situation (i.e. take | φ r | = δ r , r (cid:48) ). Seeresults in Fig. 4. As mentioned previously, we average overvarious ϕ , ϕ ’s in Eq. (7) so as to reduce finite size e ff ectsconcerning microscopic details. (a) L =
40 (1000 ϕ , ’s) (b) L =
80 (100 ϕ , ’s) FIG. 4. IPR for di ff erent eigenstates (arranged by E m < E m + ), aver-aged over phases ϕ , for each eigenstate. To check possible finite-size e ff ects, we compute two system sizes L =
40 in (a) and L = θ = π/
5. The green dashed line is thedual point V ( D ) d = From the results, there does seem to be mobility edges aswe see that the transition of IPR for di ff erent eigenstates ap-pears at di ff erent critical V ( c ) d ’s. Since the localization lengthsdo not need to be uniform, it is worthwhile to check whethersome eigenstates have much longer localization lengths thanothers such that they appear like an extended one for a finite-size sample. For this purpose, we compute systems of dif-ferent sizes for the same parameter in Fig. 4. Should sucha scenario occur, those localized states with larger localiza-tion legnths would take larger IPR values as one increases thesystem size (brighter colors should occupy larger areas). Butthere appears no notable finite-size e ff ects as Fig. 4 (a) and(b) for sizes L = ,
80 look almost identical. Thus, the IPRresults indeed suggest the existence of mobility edges for thetwo-dimensional generalization of Aubry-Andre model.Further, we note that the IPR feature here is quite di ff er-ent from what is typically seen in its one-dimensional coun-terparts. Here, the transition of IPR for di ff erent eigenstateschanges continuously, constrained only by the particle-holesymmetry (after ϕ , averages, or in L → ∞ limits) around E m = m = L /
2. That is, there ap-pear to be no “plateaus” of constant V ( c ) d for certain energywindows. This is to be compared with the one-dimensionalscenario where only a few V ( c ) d exits, related by a “step”-liketransition at certain E m ’s (see, e.g. Ref [36]).We trace such a di ff erence back to the lack of gaps in thedensity of states for our two-dimensional models. In the one-dimensional system with mobility edge, large gaps (compara-ble with band width) typically exist in the density of states.Usually, an abrupt change of V ( c ) d occurs when the eigenstatesgo from one gapped band to another, while V ( c ) d remains al-most constant within one band [36]. In our two-dimensionalcase, we note that for small size systems, such as the θ = π/ L ≤
20, gaps in the density of states do ap-pear, together with “plateaus” of V ( c ) d for E m within one band.However, as the system size L increases, those gaps in the den-sity of states shrink, and eventually vanish together with the“plateaus” of V ( c ) d . Then, a smooth change of IPR a ff ected bydensity of states at E m takes over, with typical features shownin Fig. 4. After this point, the IPR configurations no longerchange notably with the increase of L . Such a scenario showsup for all cases we have tried, including various rotation an-gles θ and stretched lattice constants for perturbing potentials.Thus, we expect that for a generic two-dimensional, quasi-periodic system in the thermodynamic limit, no gapped struc-ture with “plateaus” of V ( c ) d should appear. The lack of gaps inthe density of states is also shown for three-dimensional quasi-periodic models in Ref. [20]. It is interesting to note that sucha “mini-band” type of finite-size e ff ects also exists in certainstrongly interacting systems [37]. As such, for all results dis-cussed in this work, we have made sure that the system size iswell above the limit for band gaps to appear. D. Critical Exponent
Near the critical point V ( c ) d , the localization length is ex-pected to diverge ξ ∼ | V d − V ( c ) d | − ν (13)with the scaling exponent ν . For the two-dimensional quasi-periodic model in the orthogonal class, such an exponent hasbeen evasive in the previous studies [20, 21]. Here, by sam-pling a large number of systems with di ff erent ϕ , ’s basedon the “multifractal” analysis, we obtain such an exponent asshown in Fig. 5. One advancement of our work is that the sizes L considered are systematically larger than those in previousworks. For instance, the system sizes in Ref. [20] ranges over L = ∼ L = ∼
90 in the ex-act diagonalizations for two-dimensional systems. We noticedthat in our example of θ = π/
5, for size L ∼ , ,
60, thereare severe finite-size e ff ects as ˜ α defined below in Eq. (14)shows irregular crossovers for di ff erent system sizes. As wetested, such finite size e ff ects are much more severe than theone-dimensional Aubry-Andre model or a factorizable two-dimensional one ( θ = L already provide good scal-ing behaviors. In the following, let us describe the procedureof extracting ν in more detail.The method of multifractal finite size scaling (MFSS) wasintroduced in Refs. [38–40]. The steps can be summarizedas follows. (1) Partition the system of size L × L into boxesof size l × l , each containing ( L / l ) lattice sites. In our case, L / l =
10, so there are totally 100 boxes. (2) Obtain the real-space eigenstate φ x , y with energy close to certain value, i.e. E = α ,where we have avoided boxes on the boundary so as to reducefinite-size e ff ects,˜ α = L / l ) − (cid:80) ( L / l ) − a , b = ln A a , b ln( l / L ) , A a , b = al (cid:88) x = ( a − l + bl (cid:88) y = ( b − l + | φ x , y | . (14)Here, A a , b is the wave function amplitudes within the box in-dexed by ( a , b ). For each system size L , we need to diagonal-ize a large number of Hamiltonians with di ff erent ϕ , ϕ . Each V d ( c ) ≈ ± ν≈ ± V d - V d ( c ) ν L L = = = = = = V d α ˜ FIG. 5. The scaling of ˜ α at energy E =
0. The inset is the datacollapse of ˜ α as a function of L /ξ ∝ | V d − V ( c ) d | ν L for various systemsizes L . The number of ϕ , ϕ averaged over are 4 × for L = , , ,
180 and 2 × for L = ,
220 respectively. Thestable fit is given by expansion of Eq. (15) to the 5-th order. Thereare totally N D =
144 data points, N p = χ = p -value 0 .
52. (SeeAppendix E 1 for more fitting details.) of such a diagonalization (using Lanczos method) would ren-der one eigenstate φ x , y closest to E = α . Those ˜ α ’s for the same system sizes and V d ’s are to beaveraged over to provide one data point in Fig. 5. As empha-sized in Ref. [40], only one sample state φ x , y should be usedin each diagonalization.The qualitative behavior of ˜ α can be expected as fol-lows. When the system is completely delocalized, i.e. wavefunction amplitudes are homogeneously distributed among allsites / boxes, | φ x , y | ∼ / L d , A a , b ∼ ( l / L ) d , so ˜ α → d = d is the spatial dimension. In contrast, when approach-ing strongly localized regime with most of the boxes empty A a , b →
0, ˜ α → ∞ . Thus, near the critical regime when thesystem crosses from metallic to insulating behavior at certainenergy E , ˜ α is a monotonically increasing function of V d .The single-parameter scaling form can be written as˜ α = f ( L /ξ ) = g (cid:16) ( V d − V ( c ) d ) L /ν (cid:17) , (15)which suggests the data collapse of ˜ α in terms of L /ξ ∝ | V d − V ( c ) d | ν L . Near the critical point, we expand ˜ α ≈ (cid:80) Nn = a n w n ,where w = ( V d − V ( c ) d ) L /ν , and fit the data for di ff erent V d and L with the polynomial. The quality of fit is evaluated bythe χ statistics and the p -value, and the uncertainty range isgenerated by a Monte-Carlo type of procedure using syntheticdata sets having the same mean and standard deviation as theoriginal data. Detailed fitting procedures are described in Ap-pendix E 1.The results are shown in Fig. 5, where we see that thescaling exponent ν ≈ ν > / d = θ = π/ ν ≈ ν is simi-lar to the one-dimensional Aubry-Andre one with analyticallyobtainable ν =
1, where mobility edge is absent. As such, wecan summarize the critical behavior for the self-dual modelsin various dimensions in Table I.
TABLE I. The critical behavior for self-dual, orthogonal class ofmodels in a quasi-periodic potential.Dimension Mobility edge? Exponent ν Harris bound satisfied?1 [18] No 1 No2 Yes ≈ ≈ . IV. EXPERIMENTAL SIGNATURES
In principle, one could surely produce exactly the modelin Eq. (3) by superposing a standard square optical lattice to-gether with the Moire perturbing potential given by a digitalmirror device, like in the experiment Ref. [17]. Even with-out using the digital mirror device, a standard superlatticetechnique could realize a generalized version of the previousmodel. For the generality and broadness of our results, wewould focus on such a generalized model in the following. λ = λ = FIG. 6. The experimental scheme for the Moire lattice. Black arrowsand crosses denote the polarization directions such that the four pairsof laser beams do not interfere with each other. Here θ = π/
3. Therandom phases are ϕ = , ϕ = In the superlattice scheme, we need two identical and mutu-ally independent lattice potentials (without interference) over-lapping with each other. Note that two square lattices involve4 pairs of laser beams of the same wavelengths. Since thereare only 3 perpendicular directions for polarizations, at leasttwo pairs of the beams would end up interfering with eachother and therefore the two sets of square lattice will not bemutually independent. To circumvent this di ffi culty, we con-sider a generalization to two identical rectangular lattices,where laser beams of very di ff erent frequencies will not in-terfere with each other when their e ff ects are averaged over t / J N odd / N V d = V d = V d = V d = V d = V d = FIG. 7. The portion of particles remaining in odd x = , , . . . .Here θ = π/
3, and L = time. That means the Moire potential becomes˜ µ x , y = sin [ π ˜ u x , y − ϕ ] + sin [( π/γ )˜ v x , y − ϕ ] , ˜ u x , y = x cos θ − y γ sin θ, ˜ v x , y = x sin θ + y γ cos θ, (16)with x , y ∈ Z still denoting site indices. The self-dualitymapping can be similarly performed by replacing u x , y , v x , y → ˜ u x , y , ˜ v x , y /γ in Eq. (9), and therefore the self-dual point is still V ( D ) d =
4. In the following, we choose the aspect ratio of γ = λ = nm and λ = nm as in many experiments. Aschematic plot for the lasers and their polarization directionsare shown in Fig. 6. Note that the commensurate conditionEq. (4) still holds for rectangular lattices with rational aspectratios, as discussed in Appendix A.Consider a spin-polarized fermion system. The initial statecan be chosen similar to the previous experiments as a charge-density-wave type [16], where even sites x = , , , . . . alonga direction are empty, while odd sites x = , , , . . . are oc-cupied. Fig. 7 shows the numerical simulation for the timeevolution, where N odd are particles at sites x = , , . . . forarbitrary y , and N is the total particle number. We readily seethat close to self-duality point V ( D ) d =
4, i.e. V d =
3, there al-ready appear a significant portion of localized states such that N odd / N remain slightly above 0 .
5, which signals the mobilityedge in such a system.
V. CONCLUSION
We have provided, as the first step, a study on e ff ec-tive models of Moire localization in two-dimensional quasi-periodic systems. The commensuration conditions discussedin Sec. II can be used as a guide to engineer bilayer ormulti-layer systems made of simple rectangular / square lat-tices, which is of primary interest to cold atom experiments.Such a model also fills up the gap for the understandingof two-dimensional quasi-periodic system in the orthogonalclass (time-reversal invariant). We find that the mobility edgehere, shown by the IPR, appears quite di ff erently from that ina one-dimensional quasi-periodic system studied before. Inour case, there exists a rapid and continuous change with-out “plateaus” for IPR in the eigenstate-disorder plane. Wetrace such a di ff erence to the lack of band gaps for the two-dimensional spectrum. Note that such a feature is most no-tably revealed in the IPR calculation rather than the levelstatistics, as the latter one would inevitably involve an averageover many nearby eigenstates. On the other hand, level statis-tics reveals a cross-point for relatively small system sizes butappear to converge for larger L ’s, with the region well aboveduality point V d = ffi cient for this class of model, as we found thatchoosing su ffi ciently large system sizes would yield satisfac-tory scaling behaviors. The exponent for our two-dimensionalmodel, ν ≈
1, turns out to be the same as its one dimensionalcousin, while the di ff erence is that in two-dimensions such anexponent saturates the Harris bound ν > / d = ν ≈ .
38 [20].Finally, we have provided a scheme for cold atom experi-ments relying solely on standard superlattice techniques, with-out using the digital mirror device. Since the system is onlysensitive to the relative rotation angles rather than lattice con-stants between main and perturbing potentials, our results in-dicate that the two dimensional system is actually more viablefor cold atom experiments as there are no requirements forparticular laser wavelengths.
VI. ACKNOWLEDGMENT
This work is supported by AFOSR Grant No. FA9550-16-1-0006, MURI-ARO Grant No. W911NF-17-1-0323, andNSF of China Overseas Scholar Collaborative Program GrantNo. 11429402 sponsored by Peking University. This workused the Bridges system, which is supported by NSF awardnumber PHY190004P at the Pittsburgh Supercomputing Cen-ter (PSC).
Appendix A: Commensurate rotations for rectangular Lattices
Take an arbitrary lattice site as the center of rotation, thenthe commensurate rotations should satisfy m ˆ x + γ m ˆ y = n ˆ x (cid:48) + γ n ˆ y (cid:48) , (A1)where ˆ x (cid:48) = ˆ x cos θ + ˆ y sin θ , ˆ y (cid:48) = − ˆ x sin θ + ˆ y cos θ , and γ is theaspect ratio for the rectangular lattices. The following deriva-tions are valid for a rational γ , i.e. γ = m γ m = cos θ − sin θ sin θ cos θ n γ n . (A2) To have integer solutions ( m , m ) , ( n , n ), it is necessary tohave cos θ = k / k , sin θ = k / k , k , k , k ∈ Z . (A3)The integers k , k , k satisfy k + k = k . (A4)The above diophantine equation corresponds to finding Pon-tryagin triples whose solutions are well-known [24, 26, 41].One can enumerate them by the following procedure. Notethat Eqs. (A3) and (A4) mean finding the rational points onthe unit circle cos θ + sin θ =
1. One such point can bepinned down as (0 , ,
1) and ( q / p , q , p ∈ Z . Combining (A4)and k / k = − ( q / p ) k / k , we can set k = q − p , k = pq , k = q + p , p , q ∈ Z . (A5)Then we have Eq. (4) in the main text. Note here k , k areinterchangeable, which correspond to angles θ and π/ − θ .Note that the commensuration condition is the same for ar-bitrary aspect ratios γ . But the Bravais lattice vectors are dif-ferent, as we discuss below. Appendix B: Moire Bravais lattice vectors at commensuraterotation angles
Define the matrices S = (cid:112) q + p q p − p q = q (cid:112) q + p I + p (cid:112) q + p i σ y , S − = S T , (B1) A = q + p q − p − qp qp q − p = q − p q + p I − qpq + p i σ y (B2)where T is transpose, and A is the rotation matrix in Eq. (A2)parametrized by Eq. (A3) and (A5). They satisfy the relation S AS = I . Then the vectors m , n can be written as n γ n = α q − p + β pq , (B3) m γ m = α qp + β − pq (B4)where α, β are some rational numbers. Thus, the Moire Bra-vais vectors for a square lattice γ = T = q − p , T = pq (B5)We aim to find the shortest vectors T , as it is possible thatsome linear combinations could yield α i T + β i T = N i T (cid:48) i (B6)where N i ∈ Z and T (cid:48) i is also a Moire Bravais vector. Note here( α i , β i ) are coprime numbers, for otherwise the common fac-tors can be canceled by N i . It turns out that the only exceptionis when ( q , p ) are both odd numbers, in which case T (cid:48) = q + pq − p , T (cid:48) = q − p − q − p (B7)For a rectangular lattice with rational aspect ratio γ , oneshould apply the general conditions (B3) to find the linearindependent vectors T = ( n , n ) such that α q + β p and( − α p + β q ) /γ are integers, and q , p are coprime. For our pur-pose, we give the lattice vectors for γ =
2: when q , p are bothodd, T = ( q + p , ( q − p ) / T , T = ( − ( q − p ) , ( q + p ) / T ; when q is odd and p is even, T = ( q , − p / T , T = (2 p , q ) T ; finally,when q is even and p is odd, T = (2 q , − p ) T , T = ( p , q / Appendix C: Numerical algorithm for free fermion dynamics
Consider H = Ψ † H Ψ , Ψ = ( c , . . . , c N ) T , (C1)with c i the fermion operator at site i and H an N × N matrix.For an observable that can also be expressed in the bilinearform (such as the density) A = Ψ † A Ψ , A ( t ) = e iHt Ae − iHt , (C2)note [ AB , C ] = A { B , C } − { A , C } B , e iHt Ψ µ e − iHt = ∞ (cid:88) n = ( it ) n n ! [( Ψ † α H αβ Ψ β ) ( n ) , Ψ µ ] = ( e − i H t ) µβ Ψ β , (C3)and therefore A ( t ) = Ψ † (cid:16) e i H t A e − i H t (cid:17) Ψ . (C4)Now, if the Hamiltonian matrix can be diagonalized by thematrix T , T † HT = diag { ε , . . . , ε N } , (C5)we have A ( t ) = N (cid:88) α,β = Γ † α T † α A αβ T β Γ β e i ( ε α − ε β ) t , Γ = T † Ψ . (C6)For our purposes, consider the initial state where half of thesystem is filled, | ψ ini (cid:105) = c † . . . c † N / | (cid:105) . Then (cid:104) ψ ini | A ( t ) | ψ ini (cid:105) = N (cid:88) i = n (0) i T i α e i ε α t (cid:16) T † AT (cid:17) αβ e − i ε β T † β i , (C7)where α, β are summed over all eigenstates, and n (0) i is theinitial particle number at site i . The above formula applies toboth free bosons and fermions when the initial state is a Fockone, with the restrictions that for fermions, n (0) i ≤ Appendix D: Review of Aubry-Andr´e model in one-dimension
To be self-contained and to help trace the origin for mo-bility edge in two-dimensional models, we review some fea-tures of the one dimensional Aubry-Andre model. Thereare two necessary conditions for having uniform localizationlengths: (1) A one-dimensional model with nearest-neighborhopping; (2) The duality between localization / delocalizationin real / momentum spaces. Generalizing to higher-dimensionsclearly violate condition (1) and therefore do not exhibitunique localization length. The following derivation is basedon the discussions in Refs. [18, 29, 42, 43].The Aubry-Andre model is H = − t N (cid:88) m = ( a † m + a m + a † m a m + ) + N (cid:88) m = µ m a † m a m , (D1) µ m = V d π qm + ϕ ) , V ( D ) d = . (D2)Performing a similar analysis as in Sec. III A, we could simi-larly obtain the duality mapping with critical disorder strength V ( D ) d . However, due to the one-dimensional nearest neighborfeature, one could related the localization length with densityof states through the Thouless formula as shown below.Consider the real space matrix for Green’s function G mn ( E ) = (cid:32) E − H (cid:33) mn ≡ (cid:104) | a m ( E − H ) − a † n | (cid:105) , (D3)where m , n = , . . . , N denote lattice sites, | (cid:105) is the vacuum.For an open boundary system, the matrix E − H = E − µ − t . . . − t E − µ − t . . . − t E − µ − t . . . . . . . . . . . . . . . . . . . (D4) Due to the one-dimensional nearest-neighbor feature , the cor-ner element for the inverse matrix can be easily acquired asthe cofactor G N = E − H ) ( − t ) N − ( − N − = t N − (cid:81) N α = ( E − E α ) , (D5)where E α are eigenvalues of H . On the other hand, accord-ing to the definition, one can first expand the Green’s func-tion in the energy eigenbasis G ( E ) = ( E − H ) − = (cid:80) N α = ( E − E α ) − | α (cid:105)(cid:104) α | , and then the real space matrix element is G N ( E ) = N (cid:88) α = ψ α (1) ψ ∗ α ( N )( E − E α ) − , (D6)where ψ α ( m ) = (cid:104) | a m | E α (cid:105) is the real-space eigenfunctions. Si-multaneously using the above two equations to compute theresidual of G N ( E ) at a certain E β ,lim E → E β ( E − E β ) G N ( E ) = ψ β (1) ψ ∗ β ( N ) = t N − (cid:81) α,α (cid:44) β ( E β − E α )(D7)0Suppose the eigenstate | E β (cid:105) is localized centering at some site m , with the ansats ψ β (1) = A e − λ ( E β )( m − , ψ β ( N ) = A N e − λ ( E β )( N − m ) (D8)where λ is the exponent of interest, and A m is the multipli-cation of some constant and non-decaying oscillating factorsat site m . Combing the above two equations and taking loga-rithm, we have λ ( E β ) = N − ln( A A N ) + (cid:88) α,α (cid:44) β ln | E β − E α | − ln t . (D9)Taking the N → ∞ limit, the ln( A A N ) term vanishes. Thesecond term can be converted into an integration. Then λ ( E β ) = − ln( t ) + (cid:90) ∞−∞ (ln | E β − ε | ) dN ( ε ) , (D10)where dN ( ε ) = D ( ε ) d ε and D ( ε ) is the density of state at ε . This is the Thouless formula for the localization length ξ = /λ . Usually, the hopping constant t =
1, so the first termdrops out. This formula can be easily generalized to inhomo-geneous nearest neighbor hopping situations [29]. But notethat once the hopping goes beyond nearest neighbor, i.e. forsecond-neighbor hopping or in higher dimensions, Eq. (D4)no longer takes the same form, and Eq. (D5) cannot be ob-tained. Then the Thouless formula would not hold generally.Due to the duality mapping, one can obtain the relation N t , V d ( ε ) = N t , t Vd ( 4 t ε V d ) (D11)for an incommensurate system. Then, applying Thouless for-mula, λ t , V d ( E β ) = λ t , t Vd ( 4 tE β V d ) + ln V d t . (D12)Since λ ≥ V d > t , λ t , V d ( E β ) ≥ ln( V d / t ) >
0. That means all eigenstates | E β (cid:105) are exponentially local-ized. From the transformation between eigenfunctions in twodual representations, all eigenfunctions in the dual space arethen delocalized (localization / delocalization are opposite inthe two spaces related by Fourier transformation), meaning λ t , t / V d (4 tE β / V d ) =
0. Then, we have the same critical expo-nent (inverse of localization length ξ ) λ t , V d = ξ = ln V d t (D13)for all eigenstates when V d / t > Appendix E: Additional details for scaling exponents1. Fitting procedure
The procedure for fitting the scaling of ˜ α ( V , L ) consists ofthree steps. First, the average (cid:104) ˜ α (cid:105) ( V , L ) and standard deviation σ ( V , L )of ˜ α at each point ( V , L ) is determined. There are data of theorder 10 at each ( V , L ) calculated by using di ff erent ϕ , ϕ .Note that although the average value is straightforward to un-derstand, the “standard deviation” must be taken with care. Inall regimes, the ˜ α takes values according to a broad distribu-tion function (see, i.e. Ref. [40] Fig. 3). We reproduced sucha character for our model in Fig. 8. It is in fact the maximalpoint of the distribution function that gives the scaling behav-ior. Only when averaged over large numbers of ϕ , , the meanvalue of ˜ α would approach the maximal probability point.Thus, at each ( V , L ), we divide all data into 10 bins, obtain therespective averaged ˜ α , and compute the standard deviationusing the 10 bin-averaged ˜ α ’s. L = = L = = L = = FIG. 8. The (normalized) histogram for ˜ α at V d = . V d = . V d = . ff erent ϕ , . Forlarger systems L , the distribution inclines towards smaller ( V d in themetal limit) or larger ( V d in the insulator limit) values, while in thecritical V d they overlap. Second, we fit (cid:104) ˜ α (cid:105) ( V , L ) into the scaling function g ( L /ξ )mentioned in the main text by minimizing χ , χ = (cid:80) V , L ( (cid:104) ˜ α (cid:105) ( V , L ) − g ( V , L )) /σ ( V , L ). For a reasonable fit, the χ value should be similar or less than the number of degreesof freedom, k = N D − N p , where N D is the number of (cid:104) ˜ α (cid:105) ( V , L )to be fit for di ff erent ( V , L ), and N p is the number of fitting pa-rameters. Whether or not a fit is acceptable is decided by the p -value, p = Γ ( k / , χ / Γ ( k / , (E1)where Γ ( a , z ) = (cid:82) ∞ a t z − e − t dt is the generalized Euler-Gammafunction, with Γ ( z ) = Γ (0 , z ). A larger p -value signals a betterfit, and we take p ≥ . p -value, we also make sure the fitting stability bychecking that increasing the expansion order would result in V ( c ) d , ν within the uncertainty range (to be discussed below).Under such criterions, the expansion order is kept as low aspossible.Finally, once a best fit g ( V , L ) is obtained for (cid:104) ˜ α (cid:105) ( V , L )and σ ( V , L ), we generate 1000 synthetic data α ( syn ) i ( V , L ) , i = , . . . , V , L ). They distribute according tothe same mean and standard deviations discussed previously.Then, we fit each set of { α ( syn ) i ( V , L ) |∀ V , L } with the same ex-pansion as in the best fit g ( V , L ), and obtain a histogram of V ( c ) d , ν corresponding to i = , . . . , .
5% of data, the uncertaintyrange is obtained. The histogram of V ( c ) d , ν for θ = π/ V d ( c ) = % range = [ ] ν = % range = [ ] FIG. 9. The histogram for fitted V ( c ) d and ν at the angle θ = π/ .
5% of maximal / minimal data for V ( c ) d , ν on the twosides.
2. Scaling at another angle
We perform the same analysis for the angles θ = π/ E =
0. But here we used a wider range of system lengths L ,and it is clear from Fig. 10 that there is a shift in the “crossing”point with the increase of L . Such a “shift” diminishes as L becomes larger, which means there is an additional scaling-irrelevant term entering ˜ α [40]. Therefore, we include such aterm in the scaling function,˜ α = g (cid:16) ( V d − V ( c ) d ) L /ν (cid:17) + L −| y | g (cid:16) ( V d − V ( c ) d ) L /ν (cid:17) . (E2) Compared with Eq. (15), the L −| y | term describes the first orderexpansion of scaling-irrelevant term. With the fitting param-eter −| y | <
0, this term vanishes in the thermodynamic limit.Then we expand g , similarly as before, g = (cid:80) N m = a m w m , g = (cid:80) N m = b m w m , where w = ( V d − V ( c ) d ) L /ν . The fit-ting parameters are a m , b m , y , ν, V ( c ) d , with the total number N p = N + N + V d ( c ) ≈ + (- )(+ ) ν≈ + (- )(+ ) V d - V d ( c ) ν L L = = = = = = = =
240 L = = = V d α ˜ FIG. 10. The scaling of ˜ α for θ = π/ E =
0. The inset showsthe scaling of ˜ α in terms of L /ξ ∝ | V d − V ( c ) d | ν L after subtracting thescaling-irrelevant term (the second term in Eq. (E2)). The numbersof ϕ , being averaged over are 10 for L ≤
200 and 4 × for L ≥ / parameters has N D = , N p =
16 (with N = N = χ = , p = .
98. The uncertainty range is determined by synthetic data set asdescribed previously.[1] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov,and A. K. Geim, Rev. Mod. Phys. , 109 (2009).[2] C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao,J. Katoch, M. Ishigami, P. Moon, M. Koshino, et al. , Nature , 598 (2013).[3] B. Hunt, J. Sanchez-Yamagishi, A. Young, M. Yankowitz, B. J.LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino,P. Jarillo-Herrero, et al. , Science , 1427 (2013).[4] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxi-ras, and P. Jarillo-Herrero, Nature , 43 (2018).[5] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe,T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science , 1059 (2019).[6] W. Yao, E. Wang, C. Bao, Y. Zhang, K. Zhang, K. Bao, C. K.Chan, C. Chen, J. Avila, M. C. Asensio, J. Zhu, and S. Zhou,Proc. Natl. Acad. Sci. U.S.A. , 6928 (2018).[7] J. B. Sokolo ff , Phys. Rev. Lett. , 2223 (1986).[8] M. Krajci and J. Hafner, J. Non-Cryst. Solids , 338 (1995).[9] A. D. Villa, S. Enoch, G. Tayeb, F. Capolino, V. Pierro, andV. Galdi, Opt. Express , 10021 (2006).[10] M. J. Park, H. S. Kim, and S. Lee, 1812.09170v1.[11] N. Y. Yao, C. R. Laumann, J. I. Cirac, M. D. Lukin, and J. E.Moore, Phys. Rev. Lett. , 240601 (2016).[12] R. Nandkishore and D. A. Huse, Annu. Rev. Condens. MatterPhys. , 15 (2015). [13] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, ArXiv:1804.11065 .[14] V. Khemani, D. N. Sheng, and D. A. Huse, Phys. Rev. Lett. , 075702 (2017).[15] G. Roati, C. DErrico, L. Fallani, M. Fattori, C. Fort, M. Za-ccanti, G. Modugno, M. Modugno, and M. Inguscio, Nature , 895 (2008).[16] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. L¨uschen, M. H.Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Sci-ence , 842 (2015).[17] J.-Y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal,T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross,Science , 1547 (2016).[18] S. Aubry and G. Andr´e, Ann. Israel Phys. Soc. , 133 (1980).[19] P. Bordia, H. L¨uschen, S. Scherg, S. Gopalakrishnan, M. Knap,U. Schneider, and I. Bloch, Phys. Rev. X , 041047 (2017).[20] T. Devakul and D. A. Huse, Phys. Rev. B , 214201 (2017).[21] J. Sutradhar, S. Mukerjee, R. Pandit, and S. Banerjee, arXiv:1810.12931 .[22] A. B. Harris, J. Phys. C: Solid State Phys. , 1671 (1974).[23] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Ann. Phys. , 1126 (2006).[24] S. Shallcross, S. Sharma, E. Kandelaki, and O. A. Pankratov,Phys. Rev. B , 165105 (2010).[25] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Cas-tro Neto, Phys. Rev. B , 155449 (2012). [26] J. Steuding, Diophantine Analysis (Chapman and Hall / CRC,London, 2005).[27] C. R. Woods, L. Britnell, A. Eckmann, R. S. Ma, J. C. Lu, H. M.Guo, X. Lin, G. L. Yu, Y. Cao, R. V. Gorbachev, A. V. Kretinin,J. Park, L. A. Ponomarenko, M. I. Katsnelson, Y. N. Gornos-tyrev, K. Watanabe, T. Taniguchi, C. Casiraghi, H.-J. Gao, A. K.Geim, and K. S. Novoselov, Nat. Phys. , 451 (2014).[28] Specifically, (cid:80) x , y | φ x , y | = (cid:80) x , y (cid:80) mnkl ψ ∗ m , n ψ k , l A xmnkl B ymnkl C mnkl ,where A xmnkl = e π ix [( l − n )cos θ + ( k − m )sin θ ] , B ymnkl = e π iy [ − ( l − n )sin θ + ( k − m )cos θ ] , and C mnkl = e π i [( l − n ) ϕ + ( k − m ) ϕ ] . For m (cid:44) k or n (cid:44) l , one could use (cid:80) Lx = − L e π i β x = sin[ πβ (2 L + / sin[ πβ ]and show that D mnkl ≡ C mnkl (cid:80) x , y A xmnkl B ymnkl is finite as L → ∞ . Then, (cid:80) x , y | φ x , y | = (cid:80) m (cid:44) k or n (cid:44) l ψ ∗ m , n ψ k , l D mnkl + (cid:80) mn | ψ m , n | (cid:80) x , y → ∞ , where the first term is finite due to thelocalization of ψ m , n .[29] D. Thouless, J. Phys. C: Solid State Phys. , 77 (1972).[30] V. Oganesyan and D. A. Huse, Phys. Rev. B , 155111 (2007).[31] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Phys. Rev.Lett. , 084101 (2013).[32] L. F. Santos and M. Rigol, Phys. Rev. E , 036206 (2010). [33] L. F. Santos and M. Rigol, Phys. Rev. E , 031130 (2010).[34] A. W. Sandvik, AIP Conference Proceedings , 135 (2010).[35] L. D’Alessio and M. Rigol, Phys. Rev. X , 041048 (2014).[36] X. Li, X. Li, and S. Das Sarma, Phys. Rev. B , 085119(2017).[37] Z. Papi´c, E. M. Stoudenmire, and D. A. Abanin, Ann. Phys. , 714 (2015).[38] A. Rodriguez, L. J. Vasquez, and R. A. R¨omer, Phys. Rev. Lett. , 106406 (2009).[39] A. Rodriguez, L. J. Vasquez, K. Slevin, and R. A. R¨omer, Phys.Rev. Lett. , 046403 (2010).[40] A. Rodriguez, L. J. Vasquez, K. Slevin, and R. A. R¨omer, Phys.Rev. B , 134209 (2011).[41] L. E. Dickson, History of the Theory of Numbers, Volume II:Diophantine Analysis, pp.165-170 (Dover Books on Mathemat-ics, 2005).[42] D. Herbert and R. Jones, J. Phys. C: Solid State Phys. , 1145(1971).[43] I. Y. Goldsheid and B. A. Khoruzhenko, Israel Journal of Math-ematics148