A stabilization mechanism for many-body localization in two dimensions

Experiments in cold atom systems see almost identical signatures of many body localization (MBL) in both one-dimensional ($d=1$) and two-dimensional ($d=2$) systems despite the thermal avalanche hypothesis showing that the MBL phase is unstable for $d>1$. Underpinning the thermal avalanche argument is the assumption of exponential localization of local integrals of motion (LIOMs). In this work we demonstrate that addition of a confining potential -- as is typical in experimental setups -- allows a non-interacting disordered system to have super-exponentially (Gaussian) localized wavefunctions, and an interacting disordered system to undergo a localization transition. Moreover, we show that Gaussian localization of MBL LIOMs shifts the quantum avalanche critical dimension from $d=1$ to $d=2$, potentially bridging the divide between the experimental demonstrations of MBL in these systems and existing theoretical arguments that claim that such demonstrations are impossible.


Introduction:
The study of disordered systems has borne rich discussion and novel phenomena ever since Anderson's seminal work [1] and the subsequent theoretical observation that all single particle eigenstates of non-interacting, time-reversal symmetric systems in one and two dimensions are localized in the presence of disorder [2].Of particular interest in recent years is the phenomenon of Many-Body Localization (MBL), wherein strong disorder drives localization of the entire eigenspectrum in the presence of interactions.MBL has since been subject to intense investigation due to both fundamental and practical reasons [3][4][5][6][7].While its existence in 1D is accepted because of good agreement among numerical [8][9][10][11], analytical [12] and experimental [13,14] work, the situation in 2D remains contentious.On one hand, experimental [15][16][17] and numerical [18][19][20][21][22][23] signatures of MBL in 2D are almost identical to those in 1D [13], but on the other the thermal avalanche hypothesis (TAH) [24,25] posits that MBL cannot exist in any system of dimension greater than 1.However, the TAH relies strongly on the exponential localization of local integrals of motion (LIOMs) [26][27][28], an assumption that, as we show below, may be broken on more careful treatment of the disordered potential in these many-body systems.
In this Letter, we show that a confining potential, always present in experiments, affects the MBL transition by stabilizing the localized phase.We argue that this is a consequence of super-exponential localization mediated by the confining potential.* c2ddfcw@nus.edu.sgWe begin with a brief overview of the TAH, noting in particular the main assumption, which if broken then allows MBL to occur in 2D.We then present evidence of super-exponential (Gaussian) localization in the noninteracting picture.Using exact diagonalization (ED) of interacting spinless fermion Hamiltonians in a disordered cosine trap, we show how such a trap promotes localization, in defiance of the TAH.Our work challenges the commonly accepted view that thermal avalanches always destroys MBL in 2D.
Thermal avalanche hypothesis: In a sufficiently large system with uncorrelated disorder, it is inevitable for a region of locally weak disorder to emerge.These rare regions of weak disorder may host "thermal bubbles", regions where the system is well described by the eigenstate thermalization hypothesis (ETH) [29,30].At the interface between localized and thermal regions, interactions between the thermal bubble and individual LIOMs may trigger the thermalization of those LIOMs, thereby incorporating them into the bubble.It has been noted that the situation is highly asymmetric and the reverse, viz., LIOM's localizing the thermal bubble, rarely happens [31].The energy scales governing this thermalization are the bubble-LIOM matrix element Γ = V √ δρ and the bubble level spacing δ, where V is the interaction strength and ρ is the bubble spectral function, with the thermalization of the spin proceeding if Γ δ [24,32].The interaction decays in accordance with the decay (localization) of LIOMs through the relation , where ψ i is the LIOM at site i and ξ is the localization length.
For every additional LIOM incorporated into the thermal bubble, the level spacing roughly halves, and so δ decays exponentially with the number of thermalized spins.The number of thermalized spins itself grows algebraically with the bubble size R, δ , where δ 0 is the bare bubble spacing, A is a positive geometric constant and d is the system dimension.In this way the thermal avalanche is driven by the ever decreasing bubble level spacing while it is limited by the bubble-LIOM interaction strength.Assuming a form V (R) = V 0 exp(−(R/ξ) a ) for the interaction, where ξ is the localization length and a is a shape parameter, and that the bubble spectral function does not change dramatically with R, the general criterion for avalanche propagation at a distance R from a thermal bubble in 2D is exp AR 2 − (R/ξ) a 1 [24,25,33].We omit a dimensionless prefactor involving comparison of energy scales.The original formulation [24] set a = 1 as the authors assumed exponentially localized LIOMs, and therefore an exponentially decaying coupling between the thermal bubble and LIOMs.
The study of the TAH and the MBL-ETH transition has grown beyond the simple argument described above, with a wealth of numerical and analytical studies discussing, for example, Kosterlitz-Thouless scaling near the critical point [32], localization of the critical point itself [34], coexistence of localized and thermal regions [35], and the dynamical and transport properties [36,37].However, these build upon the basic description above with the same assumption of a = 1.
The phase diagram derived from the avalanche condition in 2D is shown in Fig. 1.Due to the R dependence of the avalanche criterion and the implicit assumption of exponential localization, a = 1, it has been argued that thermal avalanches unequivocally destroy MBL in 2D, as the quantity on the left hand side increases without  bound beyond a critical R * for a < 2. However, exponential localization, a = 1, relies on assumptions that may be broken in real experiments.For example, we argue here that the presence of a trap potentials likely alters the localization profiles to be Gaussian, a = 2.Such a change of shape qualitatively changes the behaviour of the criterion in 2D, giving rise to a critical ξ * below which thermal avalanches cannot propagate and therefore cannot destroy MBL.In such a scenario MBL survives in d=2.
Super-exponential localization: The concept of exponentially localized single-particle wavefunctions or LIOMs is rooted in two related arguments; first in the Furstenberg theorem [38], which is most applicable in 1D, and second in the forward scattering approximation to the locator expansion [1,28,[39][40][41] which can be seen as a mean-field approximation that is more accurate in higher dimensions.Both arguments consider the joint distribution of a product of individual, identically distributed (IID) elements, coming to the conclusion that a Lyapunov exponent naturally emerges characterising the exponential localization of a state.
The natural question therefore is what happens when these elements are not IID but are instead correlated, as is the case in experiments where the single-particle on-site energy consists of both disorder and confining potential terms [15][16][17].This observation opens the possibility of new types of localization and thus the question of what effect the confining potential may have on the localization is of vital importance.
To investigate the possibility of super-exponential localization numerically, we solve the 2D Anderson tight 3. Plot of the distribution, P (r), of level spacing ratios, r, for an atomic gas in a cosine trap as trap depth is increased.Without the trap, the system thermalizes despite disorder under the influence of strong, long-range interactions, resulting in a P (r) consistent with the GOE.As trap depth is increased, the system undergoes a localization transition leading to level repulsion and a Poisonnian P (r).
binding model with a confining potential, where c † i (c i ) creates (annihilates) a spinless fermion at site i = (i x , i y ), n i is the number operator, L is the number of sites in the linear dimension, w i is the disordered on-site potential uniformly drawn from [−W, W ] and is the confining potential.Angle brackets denote nearest neighbors (NN), and we impose periodic boundary conditions on both directions.The hopping term sets the energy scale, and we use W = 5.
We set L = 100 giving a Hilbert space of L 2 = 10 4 basis states, and look specifically at the eigenstate closest to the center of the spectrum.For each disorder realization, we fit the wavefunction, ψ, to both exponential and Gaussian forms in the x and y directions.We then calculate the coefficient of determination, R 2 , with values closer to 1 indicating a better fit [42], and average over 400 disorder realizations.Fig. 2 shows that at fixed disorder strength, as the confining potential depth V 0 increases, the wavefunction is better described as having a Gaussian envelope than an exponential one, suggesting that confining potentials may aid in localization by encouraging super-exponential (Gaussian) decay of the wavefunctions.
Trap-mediated MBL: Having demonstrated the change of wavefunction envelope in an Anderson localization context, we now turn on interactions to see how an ordered potential may affect the MBL transition.Owing to the exponential increase in the size of the Hilbert space, we first look at 1D.The Hamiltonian to be discussed is where H 0 is given by the 1D analogue of Eq. ( 1) with w i uniformly drawn from [−W, W ], V i = V0 2 cos 2πi L , and the second term on the right hand side of Eq. (2) describes the interactions, with g the interaction strength.Infinite range interactions are thought to suppress localization [43,44] and so we use these to demonstrate the ability of confining potentials to promote localization.As mentioned, a 1D Hamiltonian is studied to keep the Hilbert space amenable to ED, however, the model is readily generalized to higher dimensions, and we present quantum Monte Carlo (QMC) and ED results for 2D systems in the next section and Supplementary Material respectively.Eq. ( 2) is analogous to the XXZ spin chain with additional longer range diagonal interactions.For the XXZ spin chain with NN interactions, MBL is thought to occur when V 0 = 0, g = 0 at a critical W c ≈ 7 [8,45].We therefore set g = 32, W = 2 to start in the delocalized phase as we investigate the effect of varying V 0 .
We probe the localization transition by using ED to obtain the middle 1% of eigenenergies, repeated over 6400 disorder configurations, and use these to determine the probability distribution P (r) of the ratio of successive energy gaps, r i = min ri , For time-reversal symmetric Hamiltonians such as Eq. ( 2), P (r) is expected to follow the Gaussian orthogonal ensemble (GOE) in the thermal delocalized phase and a Poisson distribution in the MBL phase [46][47][48][49].For V 0 = 0 and V 0 = 4, P (r) closely adheres to that of the GOE before breaking away and transitioning to the Poisson distribution at higher values of V 0 , seen in the evolution of dashed lines in Fig. 3.These results clearly indicate a transition from thermal to localized mediated by the trap depth.
These results bear some resemblance to previous work on Stark-MBL [50][51][52], though we note that we use long range interactions and work with periodic rather than open boundary conditions, and therefore in the W = 0 limit the system does not admit localized solutions.This is in contrast to the Wannier-Stark localization that precedes Stark-MBL [53,54].The localization observed here therefore requires both disorder and appropriate confining potential.
Finite-size effects: An alternate probe to identify the MBL transition is the entanglement entropy, defined as S = −Tr A (ρ A log ρ A ) where ρ A = Tr B |ψ ψ| is the partial density matrix of subsystem A when the entire system is in the state ψ.Subsystem B is the complement of A. The trapping potential, V i , removes the "translationally symmetric" freedom to choose subsystems A and B; here we specifically choose both to span a half period of V i from maxima to minima.Results are shown in Fig. 4, where we average over the middle 1% or 100 eigenstates, whichever is fewer, and 6400 disorder configu- rations.We furthermore normalize by the Page value [55], S P = (L log 2−1)/2 to get a figure of merit between 0 and 1, where 0 indicates MBL and 1 indicates purely thermal behaviour.The normalized entanglement entropy is plotted against the aspect ratio of the trap potential, V 0 /L, for various values of L to get a sense of the transition point and the magnitude of finite size effects.The results are in agreement with those for P (r), indicating a thermal to localized transition on increasing V 0 with a critical aspect ratio of (V 0 /L) c ≈ 0.4 and the transition sharpening as L increases.
Further insight may be obtained by considering the harmonic approximation of the confining potential about its minimum, V i ≈ π 2 V 0 i 2 /L 2 for i near L/2.Identifying this potential with that of the quantum harmonic oscillator gives a trap frequency ω = π √ V 0 /L which scales as 1/ √ L at the critical point, indicating that in the thermodynamic limit, the critical trap frequency tends to 0. These results suggest that the trap-mediated localization persists in the thermodynamic limit.
2D numerics: We use the stochastic series expansion (SSE) QMC technique to investigate the 2D analog of Eq. ( 2) with NN interactions, with w i uniformly distributed in [−W, W ] and . QMC algorithms are restricted to probing only ground state properties of a Hamiltonian, H.This hurdle may be mitigated by using the Eigenstate to Hamiltonian construction [56] to find a mapped Hamiltonian, H, with the same interaction strength and trap depth but different disorder configuration, hosting the ground state of the parent Hamiltonian as an excited state of H.The ground state properties of H, estimated with SSE QMC,therefore describe the excited state properties of H of the same form [23,57,58]. Simulating a 2D system of 10 × 10 lattice sites, we calculate the participation entropy, where ψ is the many-body ground state and φ i is the i-th Fock basis state.The value S ∞ / ln N , with N the Hilbert space dimension, denotes the multifractal dimension of the state in the Fock space, whose finiteness is a characteristic signature of the MBL phase, while it is equal to 1 in the ergodic phase.We see in Fig. 5 that the multifractal dimension decreases as V 0 increases, indicating increasingly nonergodic behavior induced by the confining potential, thus confirming the 1D results.Discussion and conclusions: We have revisited the TAH and noted in particular that the underlying assumption of exponentially localized LIOMs may be broken under application of trap potentials, relevant for real experiments [15][16][17], and that the localization decay can be made Gaussian.This observation directly challenges the assertion that MBL is generically unstable to thermal avalanches in dimension greater than 1 and opens up the possibility for stable MBL in 2D, as has previously been reported in experimental [15][16][17] and numerical [18][19][20][21][22][23] studies.
We have further demonstrated that the addition of such a potential term to an MBL Hamiltonian may trigger an MBL transition in a parameter regime that would otherwise host a thermal phase, and that such a transition persists in the thermodynamic limit.Nevertheless, the reader may be concerned about whether the observed stabilization of MBL in 2D are truly the result of Gaussian localization or might have some other physical ori-gin brought on by the addition of the confining potential.Two alternative possibilities are fragmenting of the Hilbert space through emergence of new conserved quantities, and the rescaling of effective hopping amplitudes as the confining potential impedes particle motion.We show in the Supplementary Material that neither of these explanations account for the observed localization.Our results point to external trap potentials and the concomitant Gaussian localization being the most likely stabilization mechanism for MBL in 2D and in principle able to overcome the TAH.
Potential avenues of further work include the study of the effect of trap potentials on the dynamics of disordered systems, and experiments to verify stability of the MBL phase against thermal avalanches by changing the trap shape or depth.For example, a 2D system confined with a cosine, or locally quadratic, potential in one direction and a uniform square well trap [59] in the orthogonal direction could give an anisotropic avalanche propagation.
FIG. 6. Plot of the distribution, P (r), of level spacing ratios, r, for a 2D atomic gas described by Eq. ( 5) as trap depth is increased.Data gathered using ED.Without the trap, the system thermalizes despite disorder under the influence of strong, long-range interactions, resulting in a P (r) consistent with the GOE.As trap depth is increased, the system undergoes a localization transition leading to absence of level repulsion and a Poisonnian P (r).
Exact diagonalization results in 2D: Given the main argument that Gaussian localization may overcome thermalization in 2D, it is prudent to present evidence of MBL in a 2D system, even if the system size accessible to us is small.A straightforward extension of Eq. ( 2) gives with w i uniformly drawn from [−W, W ] and , which we exactly diagonalize in analogy to Fig. 3 for a 4 × 4 system, all other physical and numerical parameters being equal.Results are shown in Fig. 6, where the same GOE to Poissonian evolution of P (r) on increasing V 0 is seen as in the 1D study of the main text, though the trap depth required here is much higher.This is due to the increased coordination of the hopping and interaction terms promoting thermalization.However, due to the small system size, the lattice model cannot resolve the cosine confining potential and at L = 4 it looks like a sawtooth wave.There is therefore a need for alternate techniques such as the QMC -EHC approach that we use in our letter to properly investigate 2D systems.
Multiple wells: There is no need a priori to identify the confining potential wavelength with the system size, and we could for example make the replacement L → kL for the system size with k ∈ Z + in Eq. ( 1) and maintain the periodic boundary conditions.The quantity k would then count the number of wavelengths spanned by the system.The results obtained for k > 1 and finite disorder had the wavefunction localized about a single point in space with no obvious periodicity, and were indistinguishable from k = 1, so we choose k = 1 for computational expedience.On the contrary, at zero disorder, the eigenstates were periodic in L as expected from Bloch's theorem.We therefore conclude that for periodic systems, the breaking of discrete translational symmetry by inclusion of disorder is necessary for localization about a single point.Such concerns have thus far not been addressed in studies of disorder-free Stark-MBL [50][51][52]54] where open boundary conditions are used.
Other potentials: Further investigation into the importance of Gaussian localization envelopes may be conducted by varying the form of the confining potential present in Eqns.( 1) and (2).We choose two simple periodic forms, where Θ(x) is the Heaviside step function and the 2D versions are readily generalized as V i = 1 2 (V ix + V iy ).The sawtooth wave in particular is reminiscent of the  6)) and square (lower, Eq. ( 7)) potential respectively.2 sets of data are shown in each case to show the fit done in 2 dimensions.Wavefunction fitting done on a noninteracting system.
constant gradient Stark field seen in Stark-MBL studies [50][51][52] but we repeat that the previous studies were performed with open boundary conditions and that in the non-interacting disorder free limit, while the open system admits exponentially localized Wannier-Stark orbitals, the periodic system in contrast does not admit localized solutions as a consequence of Bloch's Theorem.
In the non-interacting case, we solve for the lowest energy free state as in Fig. 2 to determine the shape of the localized wavefunctions.The results in Fig. 7 show that the sawtooth potential admits wavefunctions of indeterminate character while the square potential hosts wavefunctions that are definitively better described as exponential rather than Gaussian.
Turning on interactions, we see from Fig. 8 that the system under the influence of the sawtooth potential is localized to a similar degree as that of the cosine potential, while the square potential does not appear to promote localization to the same degree, in agreement with our assertion on the importance of Gaussian localization.
Hilbert space fragmentation: Concomitant with application of a square potential in particular is the emergence Plot of the distribution of level spacing ratios r for an atomic gas under a cosine (red, Eq. ( 2)), sawtooth (yellow, Eq. ( 6)) and square (green, Eq. ( 7)) at fixed trap depth, and under no potential but with effective (reduced) hopping amplitude (blue, Eq. ( 10)).Only the cosine and sawtooth potentials show appreciable agreement with the localized Poissonian prediction.
of a new conserved quantity, that is the proportion of particles in the low (or high) potential regions, fragmenting the Hilbert space into weakly connected regions that become truly separate in the infinite V 0 limit.One must take care therefore when gathering level statistics as aggregating samples from disjoint sectors, that individually adhere to the GOE prediction, may result in a total sample more reminiscent of the Poisson distribution and thus a spurious indication of localization [60].We see in Fig. 8 that this has not happened, with the level statistics of a system under a square potential not being well described by the Poissonian prediction, in particular having no evidence of strong level repulsion.This is due to the new conserved quantity having a strong energy dependence and so by sampling in the middle of the spectrum, we exclusively sample the sector of Hilbert space that has half of the particles in the low potential region and half in the high potential region.Such emergent conserved quantities are not present in continuously varying potentials and so we conclude that the level statistics seen there are not due to spurious sampling of a fragmented Hilbert space.
Effective hopping: Another natural question to ask is whether the observed localization is a consequence of the potential energy term rescaling the kinetic hopping term, which was set to unity in Eq. ( 2).This idea is pertinent to the highly excited states we consider, where the particle(s) may be classically free and the (positive) difference between total and potential energy may be interpreted as kinetic energy.We start by considering the non-interacting Schrödinger equation and attempt to obtain a space-varying effective hopping t(x) that encodes the effect of the applied potential V (x) while maintaining the eigenstate ψ n (x) where n is just some (set of) good quantum number(s) to label the state.The result is and for V (x) everywhere negative and focusing on classically free states, E n > 0, the effective hopping element is bounded as 0 < t(x) < 1, indicating a reduced kinetic energy and thus, naively, an increased tendency for localization.We furthermore see that since t(x) depends on E n and thus on the state labels n, there is in general no effective hopping model that preserves both the spectrum and eigenstates of the tight-binding model with confining potential.We can however select a particular E n to exactly recover one eigenstate and approximately recover others close in energy.To investigate this, we perform the same ED study on the effective-hopping Hamiltonian where we now add interactions, E mid is the energy of the exact eigenstate at the middle of the spectrum and the form of t i comes from Eq. ( 9).For the system parameters given, E mid < 0 generically for any disorder configuration, and so with V i ≥ 0, we have 0 ≤ t i ≤ 1.This would be expected to push the system towards localization, however as seen in the blue line of Fig. 8, P (r) remains consistent with the GOE prediction even at high effective V 0 and so a simple rescaling of the hopping terms consistent with the applied potential V i is insufficient to trigger a transition to the MBL phase.

FIG. 1 .
FIG. 1.The thermal avalanche hypothesis phase diagram in two dimensions for the many body localization to eigenstate thermalization hypothesis transition as a function of localization length, ξ and the local integral of motion (LIOM) shape parameter a, where a = 2 is critical.LIOM shapes (shaded blue) at both a = 1 and a = 2 are shown on the left, with sample localization potentials (black lines) producing those LIOMs shown below each.The confining potential is shown along with the total potential as a guide to the eye.The overlap of adjacent states (shaded red) highlight the qualitative difference between Gaussian and exponential localization of the LIOMs.

18 W = 2 , g = 32 FIG. 4 .
FIG.4.Entanglement entropy of half the system against trap aspect ratio with increasing system size.Shallow traps permit a thermal phase while deep traps promote localization, with a crossover aspect ratio of (V0/L)c ≈ 0.4, corresponding to a vanishing trap frequency ω = π √ V0/L in the thermodynamic limit.The errors of the curves are within 1% and are omitted for clarity.

1 FIG. 5 .
FIG. 5. Participation entropy of the ground state of 2D NNinteracting spinless fermions with on-site disorder in a cosine confining potential, calculated by SSE QMC.The participation entropy decreases with increasing trap depth, indicating increased nonergodicity.The black dashed line at S∞/ ln N = 1 indicates ergodic behaviour.

FIG. 7 .
FIG.7.Plot of the coefficient of determination R 2 when fitting the absolute value of the lowest energy free state |ψ| to a Gaussian (blue) or an exponential (red) as the depth of a confining potential V0 is increased for a sawtooth (upper, Eq. (6)) and square (lower, Eq. (7)) potential respectively.2 sets of data are shown in each case to show the fit done in 2 dimensions.Wavefunction fitting done on a noninteracting system.