Interaction-enhanced many body localization in a generalized Aubry-Andre model

We study the many-body localization (MBL) transition in a generalized Aubry-Andre model (also known as the GPD model) introduced in Phys. Rev. Lett. 114, 146601 (2015). In contrast to MBL in other disordered or quasiperiodic models, the interaction seems to unexpectedly enhance MBL in the GPD model in some parameter ranges. To understand this counter-intuitive result, we demonstrate that the highest-energy single-particle band in the GPD model is unstable against even infinitesimal disorder, which leads to this surprising MBL phenomenon in the interacting model. We develop a mean-field theory description to understand the coupling between extended and localized states, which we validate using extensive exact diagonalization and DMRG-X numerical results.

Introduction.-Interactingmany-body systems should manifest eigenstate thermalization hypothesis (ETH) with long-time thermalization as the system becomes ergodic [2][3][4].ETH, however, seems to fail generically in disordered or quasiperiodic 1D systems, at least for finite-size systems, although the situation is unclear in the thermodynamic limit.Such a phenomenon is known as many-body localization (MBL) [5][6][7].MBL has been numerically and experimentally verified in numerous 1D disordered and quasiperiodic interacting systems with the generic finding that the interacting system undergoes an MBL transition at large disorder for a fixed interaction strength.In general, the disorder strength necessary for inducing MBL is much larger in the interacting system (and increases with increasing interaction) than in the corresponding noninteracting system.For example, the noninteracting 1D Anderson model is localized for any finite disorder whereas the corresponding interacting Anderson model undergoes MBL for large disorder [8].This is expected as interaction tends to thermalize the system by sharing energy and information among the constituents, leading to ergodicity as postulated in ETH (all noninteracting systems are trivially nonergodic whereas the nonergodicity in MBL is nontrivial violating ETH).Another example is the well-known Aubry-Andre (AA) quasiperiodic model which has all states localized or extended for a critical "disorder" (i.e., quasiperiodic potential strength) larger or smaller than 1 for the noninteracting model [using convention in Eq. ( 2) below], whereas the corresponding interacting model exhibits MBL for disorder larger than 1.7 [9][10][11][12].
However, an important question is whether exceptions to the above scenario can be constructed.To address this question, here we describe and analyze a surprising counter-intuitive situation arising in the GPD model [1], where finite interactions may lead to enhanced MBL in the sense that MBL occurs "early" with the critical disorder strength for the MBL in the interacting GPD model being lower than that in the noninteracting model.This is unexpected as interaction is thought to oppose MBL, and not induce it.Since the GPD model has already been studied experimentally [13], our predictions are directly verifiable in the laboratory.Our work provides important new insights into the competition between localized and extended degrees of freedom in interacting many-body dynamics.
The GPD model.-Westart with the GPD model, introduced in Ref. [1]: where t is the hopping strength, which serves as the energy unit in this work.The on-site potential is where q = (3 − √ 5)/2 (equivalent to the golden ratio), V is the potential strength, α ∈ (−1, 1), and ϕ ∈ (0, 2π) is a random phase.
For α = 0, the GPD model reduces to the AA model, which carries no single-particle mobility edge (SPME).However, for nonzero α, the GPD model is known to carry an exact SPME given by αE = 2(1 − V ) [1].To show this, we calculate the fractal dimension of the singleparticle eigenstates in Fig. 1 for α = −0.8.The fractal dimension is defined by Γ = − ln j |ψ j | 4 / ln L, where ψ j is the wave function.Consequently, we have Γ → 1 for extended states and Γ → 0 for localized states.As shown in Fig. 1, the extended and localized states are separated by the exact SPME, and the system is fully localized at V c = 2.
An "early" MBL.-We now consider the spinless interacting GPD model at half-filling [14][15][16][17].Specifically, the interacting GPD Hamiltonian is given by H = H GPD + U j n j n j+1 , where U is the interaction strength.Throughout this work we will take ϕ = 0, α = −0.8, and adopt the open boundary condition, unless specified otherwise.
We use two standard diagnostics to obtain the MBL phase diagram in this model: the entanglement entropy (EE) and the mean gap ratio.Here, we use the half-chain second Renyi entropy given by S = − ln tr L ρ 2 L and ρ L = tr R |ψ⟩⟨ψ|, where tr L , tr R denote respectively the partial trace over the left and the right half of the system.In addition, the mean gap ratio ⟨r⟩ is defined as the average value of r i = min{δE i , δE i+1 }/ max{δE i , δE i+1 }, where δE i = E i+1 − E i is the energy gap between two adjacent energies.In the thermal phase, the EE of the eigenstate approaches the Page value S T = (L ln 2 − 1)/2 [18] and ⟨r⟩ = 0.53 for the Gaussian orthogonal ensemble.In the MBL phase we have S/S T → 0 and ⟨r⟩ = 0.38 for the Poisson distribution.
In Fig. 2(a) and 2(b), we calculate the EE and the mean gap ratio averaged over the whole spectrum in an L = 16 system.In Fig. 2(c) and 2(d), we obtain the same quantities for U = 1 in systems of varying sizes.One prominent feature in these results is that the system seems to have an "early" MBL transition.To see this, we focus on the 0.85 < V < 2 region in Fig. 2(a) and 2(b).For very weak interactions (U ≪ 1), the mean gap ratio is always 0.38, while the EE scales as S ∝ ln L in the presence of extended states.Such features are consistent with the single-particle properties of the GPD model.However, as the interaction U increases, the system becomes localized rather than thermalized, as shown by the drastic decease of the EE.This interaction-induced localization is also confirmed by the mean gap ratio.To further verify this "early" MBL transition, we take U = 1 and analyze the finite-size effect in Fig. 2(c) and 2(d).Both of them indicate an MBL transition around V = 0.75, which is much smaller than the full single-particle localization transition at V c = 2.We thus see that the MBL transition in this model happens consistently at a critical disorder strength well below the single-particle localization point.Such a feature is in stark contrast with the MBL transition found in all other models, where MBL happens at a critical disorder substantially larger than that for the noninteracting case [5-8, 11, 12].
Fragile flat bands.-Toexplain this surprising "early" MBL transition, we take a closer look at the singleparticle properties of the GPD model.From Fig. 1 we see that for V ∈ (0.85, 2) most eigenstates in the model are already localized, and only the highest-energy band remains extended.Moreover, we find numerically that this flatband contains approximately qL states, so nearly 40% of the states are extended for V ∈ (0.85, 2).Moreover, compared to other bands, the width of this extended band is very small.It turns out that this extended band is very fragile against external perturbations, which distinguishes the GPD model from some other extensively studied models, such as the AA model [19].
We demonstrate this feature by adding to the singleparticle GPD model an additional weak random disorder, where We find that, as a result of the above perturbation, the highest-energy band in the GPD model becomes completely localized, causing the full localization transition point to decrease from V c = 2 to below V = 1, as shown in Fig. 3(a).Intuitively, this observation can be understood by noticing that the width of the flatband w is much smaller than the band gap ∆ between this flatband and the other states, making this flatband susceptible to perturbations.In particular, as long as the disorder satisfies w ≪ δV ≪ ∆, first-order perturbation theory leads to the following effective Hamiltonian, where P is the projection operator of the flatband.Therefore, the system tends to localize the flatband states in the presence of any disorder, as shown in Fig. 3(b).These localized eigenstates constitute a localized basis of the flatband, so they can be regarded as the Wannier functions for this flatband.Keep in mind, however, that there is no standard definition of Wannier function in a quasiperiodic system.Consequently, the shape of the Wannier functions is highly dependent on the perturbation we select.Nonetheless, the localization centers of any two Wannier functions are separated by at least two lattice sites, regardless of the choice of perturbation.To quantify the extent of localization, we extract the localization length ξ of the Wannier functions by fitting the wave function to ψ j ∝ exp(−|j − j 0 |/ξ), as shown in Fig. 3(b).Our results in Fig. 3(c) show that all the Wannier functions are deeply localized states with ξ ∼ 1.Hence, the single-particle GPD model with α = −0.8resembles a deeply localized system, disguised by weak tunnelling that can be destroyed by disorder.
In fact, the fragile flatband is only prominent for large α (i.e., |α| close to 1) and small q, and therefore our discussion below is strictly carried out in this limit.
A mean-field description for the transition.-The"early" MBL transition suggests that the extended orbitals in the flatband, amounting to 40% of all the orbitals, are localized by the other localized singleparticle states.Based on our discovery of the fragile flatband in the GPD model, an intuitive argument for the "early" MBL transition is that quantum fluctuations coming from the localized orbitals serve as the additional disorder and localize the extended flatband as the interaction is turned on.This heuristic viewpoint can be theoretically formulated by a mean-field (MF) theory.The MF theory is essentially a set of nonlinear selfconsistent equations, and the physics can be explained by the MF theory if the solution of the self-consistent equations agrees with the numerically calculated manybody eigenstates.
Let us construct a complete and deep-localized basis utilizing the Wannier functions together with the other localized orbitals.Each lattice site can be associated with a unique basis orbital localized on this site.Under this basis, the model of Eq. ( 1) can be rewritten as where f k annihilates the basis orbital localized on site j, and ñj = f † j f j is the particle number of the orbital j.In addition, j ∈ flat denotes the Wannier functions of the flatband, and ′ ijkl denotes that three of the four indices ijkl are different.The first term in Eq. ( 5) is the energy of the orbital, and the second is the diagonal part of the interaction in this basis, serving as the additional disorder.These two terms contribute the dominant part of the Hamiltonian.The third term comes from the weak tunnelling between the Wannier functions, and we estimate that t ij ∼ w/2 ∼ O(10 −1 ).Finally, the last term is the off-diagonal part of the interaction, and thus Hence, the last two off-diagonal terms are perturbative, suggesting that the MF theory should work well.
Although we choose a specific basis in this argument, the self-consistent equations are basis-independent, and so are the MF solutions.Our numerical results below show that the precision of the MF solution is much better than O 10 −1 .The above argument also applies to MBL in other models, and we specifically show for the AA model [19].
Accuracy of the MF solution.
-We now demonstrate numerically the accuracy of the MF theory.We fix V = 1.5, which is smaller than the single-particle localization transition point.The standard MF theory is targeting the low-temperature physics, and thus we need to calculate the ground state of the self-consistent equations.
However, to analyze the "early" MBL transition, we must obtain the highly excited solution of the self-consistent equations, which is notoriously difficult in generic models.In the present problem, when the system is in the MBL phase, we are able to construct an efficient algorithm to solve the MF equations for the highly excited states, as we now explain.To obtain the ground state, one generally starts from an initial product state, replaces the nonlinear terms with the solution in the previous round of iteration, and minimizes the energy in each round of iteration.To derive the highly excited states, in contrast, we start from the product state of the localized orbitals, replace the nonlinear terms, and maximize the overlap between the next state and the current state, a procedure resembling the DMRG-X algorithm [20,21].If the iteration converges, we then generate an MF state from the initial Fock state.We emphasize that despite the similarity between DMRG-X and our modified MF theory, these two algorithms focus on very different aspects.First, although DMRG-X is expected to provide a more quantitatively accurate approximation, the MF theory provides a more intuitive understanding of the physics behind the numerical results.Second, the MF theory is much faster numerically than the DMRG-X algorithm.Finally, the MF theory can actually provide more accurate initial states for the DMRG-X algorithm.
To illustrate the power our theory, we first study a particular MF solution, the MF Néel state |Z 2 ⟩ generated by the Néel-like Fock state In Fig. 4, we compute the fidelity of the MF Néel state F (t) = |⟨Z 2 |Z 2 (t)⟩| in an L = 18 system using exact diagonalization (ED) and L = 26 system using the kernel polynomial method (KPM) [22].The fidelity in both systems remains surprisingly high (> 0.997) and does not decay after 1000 tunnelling times.This implies that the MF Néel state is fairly close to one of the manybody eigenstates.To compare the MF Néel state with the eigenstates, we use ED to obtain all eigenstates in the L = 18 system and choose the eigenstate with the highest overlap (= 0.9994) with the MF solution.As shown in Fig. 4, the density profile of the exact eigenstate and that of the MF Néel state are almost indistinguishable.For the L = 26 system, we utilize the DMRG-X algorithm with the MF Néel state being the initial state to generate the eigenstate.The DMRG-X algorithm gives us a rather accurate result with the energy uncertainty of H 2 − ⟨H⟩ 2 = 2.9 × 10 −9 .Similar to the result in the smaller system, one can barely differentiate the MF and the DMRG-X results.We also The probability density f (δn) of δn in the GPD model, and the inset zooms into δn ∈ [0, 0.05].For L = 14, 16 and 18, we consider all the sites and all the eigenstates obtained by ED.For L = 26, we consider all the sites and use DMRG-X to generate 2000 random eigenstates.Here, we take U = 1.
check that the overlap between the DMRG-X and the MF results is |⟨ψ MF |ψ DMRG-X ⟩| = 0.9992.Note that in general the density profile in the L = 18 chain is quite similar to in the L = 26 chain, which shows that system is deep in the MBL phase.
MF theory for generic To further validate the MF we now investigate the agreement between generic eigenstates and their corresponding MF Given we are focusing on the MBL regime, it is more appropriate to characterize the quality of the approximation using local quantities, such as particle density, rather than global quantities, such as the overlap.Particularly, we plot f (δn), the probability density of δn in Fig. 5 with normalization ´f (δ) dδn = 1.This quantity is defined as the absolute value of the particle number difference on a random site between a random eigenstate and its corresponding MF state.To efficiently obtain the corresponding MF solutions of a given eigenstate, we use the one-body reduced density matrix of the eigenstate ⟨E|c i c † j |E⟩ as the initial state for the iteration of the self-consistent equations.In practice, we find that about 2% of the states do not converge at all, and about 10% of them converge to the MF solution corresponding to the other eigenstates.Even including these 12% of solutions, we find in Fig. 5 that δn is highly concentrated around zero.What is more, the statistics of the quantity manifests almost no finite-size effects, which can be seen by comparison with the ED results in the smaller systems and the DMRG-X results in the larger system.Hence, the numerical results verify that the MF theory is a general explanation for the "early" MBL transition in the GPD model.
Conclusion.-In this work we studied a surprising interaction enhanced MBL, which is a special feature of the GPD model.This phenomenon is explained by a MF theory construction, which is validated by comparison with numerical results obtained by ED and DMRG-X algorithms.Given the recent experimental implementation of the GPD model [13], our predictions can be experimentally verified.
We now show that the existence of an extended and narrow band of eigenstates in the GPD model is an intrinsic and unusual property of this model.We illustrate this point by comparing the GPD model with the AA model and the Anderson model.
The pristine model considered here is where V j is given in the main text.We still take α = −0.8 for the clean GPD model, while the clean AA model is the α = 0 case.Their corresponding disordered models are obtained by adding δV introduced in in the main text.Further, we also study the weak disordered Anderson model as a benchmark, Although infinitesimal disorder can localized 1D systems, weak disorder will result in a very long localization length, and the system still seems extended in a finite small system as shown in Fig. S1.This is in clear contrast to the disorder GPD model where the system is obviously localized even in a small system.Hence, the instability should be an intrinsic property of quasiperiodic systems.

II. MF THEORY IN THE AA MODEL
The MF argument in the main text is based on the existence of a set of deep-localized orbitals and perturbative off-arXiv:2305.20090v1[cond-mat.dis-nn]31 May 2023 diagonal terms, which is generally satisfied in the deep MBL phase of various models.To test the universality of the MF theory, we study the MBL phase in the AA model.The MF theory in the AA model is much simpler than that in the GPD model.First, as the AA model has no mobility edge, it is impossible that the localized states serve as disorder.Second, even if there are external random disorders, the AA model is stable against weak random disorders as shown in Fig. S1.Hence, the MBL transition in the AA model is much later than the single particle transition, meaning that all single particle orbitals are deep localized states and that we do not need to construct Wannier functions.
Here, we take V = 2.5 because the AA model with U = 1 is known to be localized around V ≈ 1.7.In Fig. S2(a) and S2(b), we calculate the MF Néel state in the AA model, and compare them with the ED and the DMRG-X.The high overlap (0.9260 for L = 18, and 0.9644 for L = 26) and the agreement of the density profile validate the MF argument in the AA model.Note that similar to the GPD model in the main text, the particle density on the first 16 sites in the L = 26 system resembles that in the L = 18 system.This is a direct consequence of the localization because the particles on the left can hardly feel the effect of the additional patch from site 17 to site 26 on the right.Nonetheless, the region close to site 17 can still feel the additional patch, so the density close to site 17 is general different as shown by Fig. S2(a) and S2(b).In fact, it is a coincidence that the MF Néel states in the GPD model even agrees well around site 17.Furthermore, we also calculate the δn in Fig. S2(c).The result is quantitatively similar to Fig. ??, manifesting an excellent precision of the MF solutions.

III. THE INEFFECTIVENESS OF THE NN INTERACTION
In the main text, we demonstrate the localization effect caused by the localized single-particle orbitals.However, if the state only occupies the extended orbitals, then there is no effective disorder to localize the system.Note that this phenomenon never appears at half filling, because the flatband only contains 38% of the orbitals and some localized orbitals must be occupied.Moreover, as the localization centers of the Wannier functions are at least separated from each other by two sites, the NN interaction has very weak effect within the flatband.Hence, we anticipate that the state remains extend in the absence of the localized orbitals.
To numerically confirm it, we study the system of length L = 21, which possesses 8 extended orbitals in the flatband, and calculate the following imbalance dynamics, where ñi (t) is the particle number operator ñi in the Heisenberg picture.In Fig. S3, we calculate the time-averaged imbalance I(t) = t −1 ´t 0 I(s)ds.Both the noninteracting case and the NN interaction case relax from I = 1 to I = 0.2 at late time.Moreover, they manifest similar dynamics, validating the weak effect of the NN interaction within the flatband.

FIG. 1 .
FIG.1.The fractal dimension of the single-particle GPD model (U = 0).Here, we take L = 610 and use the periodic boundary condition to avoid localized edge states.Here (and in all the figures) the α of Eq. (2) is −0.8, and the line defining SPME is given by E = 5(V − 1)/2.
FIG. 2. (a) and (b) show the entanglement entropy and the mean gap ratio averaged over the whole spectrum in an L = 16 system as a function of U and V .(c) and (d) show the same two quantities averaged over the middle quarter of the spectrum for U = 1 in systems of various sizes.The results in (a) and (b) are averaged over 6 random phase realizations, and the results in (c) and (d) are averaged over 1000, 200 and 10 random phase realizations for L = 14, 16, and 18, respectively.
FIG. 3.(a) The fractal dimension of all eigenstates in a single-particle GPD model under the perturbation of Eq. (3).(b) The wave function of one Wannier function and its fitting.(c) The localization length of the localized eigenstates and the Wannier functions.Here, the system size is L = 610.

FIG. 4 .
FIG. 4. (a) and (b) show the fidelity F (t) = |⟨ψ|ψ(t)⟩| of the MF Néel state in an L = 18 (using ED) and an L = 26 system (using KPM), respectively.(c) and (d) show the expectation of the particle number in the MF Néel states and their corresponding eigenstates derived by ED or DMRG-X.The system size is L = 18 in (c) and L = 26 in (d).For the DMRG-X algorithm, the energy uncertainty is H 2 −⟨H⟩ 2 = 2.9 × 10 −9 .Here, we take U = 1.

) 7 FIG
FIG. S1.NPR of single-particle eigenstates in the AA model, the GPD model, and the Anderson model.Here, the system size is L = 2584.
FIG.S3.Time-averaged imbalance for the noninteracting and the NN interaction cases.Here, we take L = 21 and N = 4 for halffilled flatband, and U = 1 for the NN interaction.