Probing Slow Relaxation and Many-Body Localization in Two-Dimensional Quasi-Periodic Systems

In a many-body localized (MBL) quantum system, the ergodic hypothesis breaks down completely, giving rise to a fundamentally new many-body phase. Whether and under which conditions MBL can occur in higher dimensions remains an outstanding challenge both for experiments and theory. Here, we experimentally explore the relaxation dynamics of an interacting gas of fermionic potassium atoms loaded in a two-dimensional optical lattice with different quasi-periodic potentials along the two directions. We observe a dramatic slowing down of the relaxation for intermediate disorder strengths and attribute this partially to configurational rare-region effects. Beyond a critical disorder strength, we see negligible relaxation on experimentally accessible timescales, indicating a possible transition into a two-dimensional MBL phase. Our experiments reveal a distinct interplay of interactions, disorder, and dimensionality and provide insights into regimes where controlled theoretical approaches are scarce.


I. INTRODUCTION
The ergodic hypothesis underlies quantum statistical mechanics, linking reversible microscopic dynamics to irreversible macroscopic behavior. In an ergodic system, local degrees of freedom get rapidly entangled with one another, and local quantum correlations are rapidly erased [1][2][3][4]. Nonergodic many-body localized (MBL) [5][6][7][8][9][10] systems, however, defy this ubiquitous behavior and show persistent local quantum correlations [11][12][13][14]. Furthermore, MBL systems are believed to be robust to small, local perturbations and form a distinct, nonergodic phase of matter. The phase transition from the ergodic phase to the MBL phase appears to be a highly unusual critical phenomenon; as ergodicity breaks down in the MBL phase, its description lies beyond the scope of thermodynamics and traditional statistical physics [8,9].
Because of limitations of the available numerical methods, most theoretical explorations of MBL concentrate on one dimension. Whether, and under which conditions, MBL can occur in higher-dimensional systems remains a challenging question for both theory and experiment. While the initial theoretical work in Ref. [5] on MBL does not depend strongly on dimensionality, it was recently argued that rare, locally thermal regions [15] in systems with true random disorder can destabilize the MBL phase in two dimensions. It is presently unclear if such arguments also hold for systems with deterministic disorder such as quasiperiodic potentials. At the same time, initial experiments provided evidence for a MBL phase in higher dimensions by measuring global transport [16,17]. Moreover, the nature of a possible MBL transition in higher dimensions might itself be very different compared to the one-dimensional transition [18][19][20][21][22][23][24]; for example, a subdiffusive phase as a precursor to localization in one dimension [18,[24][25][26][27][28] might not exist in higher dimensions [29] (but see Ref. [30]). Given the apparent conflict of available theoretical results [15,29,31,32] and infeasibility of reliable numerical simulations, experiments stand to play an important role in elucidating these regimes [10,16,17,33,34].
Ultracold atoms in optical lattices provide a particularly well-suited platform to explore these phenomena, as they combine almost ideal isolation from the environment with individual experimental control of all microscopic parameters. In this work, we employ ultracold fermions in a quasiperiodic optical lattice to experimentally investigate the appearance of a nonergodic many-body phase in two dimensions by directly tuning the strength of a quasiperiodic potential. By quantifying the dynamical relaxation of an imprinted striped density pattern, we find evidence for three dynamical regimes: a regime of fast relaxation at weak disorders consistent with thermalization; a regime of slow relaxation at intermediate disorders, resembling the relaxation expected in a Griffiths regime [23]; and, finally, a strong-disorder regime with negligible relaxation, consistent with the appearance of a MBL phase. The slow relaxation regime begins only once the single-particle states are already strongly localized, highlighting that the slow dynamics is an inherent interaction effect. Compared to one dimension [27], the slow relaxation regime is observed to much stronger disorders in two dimensions, revealing an important role of dimensionality. Furthermore, tracking the relaxation dynamics appears to be useful in locating the many-body localization transition, even in the presence of weak couplings to the environment [27,35].

A. Experiment and model
Our system is composed of a degenerate 40 K Fermi gas prepared in an equal two-component spin mixture of its two lowest hyperfine states. The spinful fermions hop on a square lattice, and the two species interact via on-site interactions that are tunable by a Feshbach resonance. Two quasiperiodic potentials with different incommensurabilities are created along the x and y directions of the lattice and form a quasiperiodic two-dimensional disorder potential; see Fig. 1. Our system is described by the following Hamiltonian: Here,ĉ † i;σ ðĉ i;σ Þ is the creation (annihilation) operator of a fermion with spin σ ∈ fj↑i; j↓ig on a lattice site i ¼ ðm; nÞ, characterized by the Cartesian coordinates ðm; nÞ, andn i;σ ¼ĉ † i;σĉ i;σ is the particle number operator. In the first term, the angle brackets h; i restrict the sum over nearest-neighbor sites. The tunneling matrix element is set to J ≈ h × 300 Hz (h is Planck's constant), and U denotes the on-site interspecies interaction strength. The disorder potential is characterized by the strength Δ and the incommensurable wavelength ratios β x ≈ 0.721 and β y ≈ 0.693 [36]. In the absence of interactions, this system is separable along the two directions and admits an Aubry-André-type metal-insulator transition at a critical disorder strength of Δ U¼0 c ¼ 2J [37]. To probe the many-body dynamics of this system, we prepare a far-from-equilibrium initial state where atoms are selectively loaded only on the even stripes; see Fig. 1(a).
In an ergodic time evolution, this density-wave pattern will quickly vanish as the dynamics erase the microscopic details of the initial conditions. In contrast, a persistent pattern indicates a memory of the initial state and hence nonergodic behavior. This can be captured by the normalized atom number difference between the even N e and odd N o stripes, defined as the imbalance I ¼ ðN e − N o Þ=ðN e þ N o Þ, which serves as our dynamical order parameter. Such an observable has several key advantages. Whereas mass transport is a slow process even in clean ergodic systems [38], the imbalance relaxes within a few hopping times [10,39]. Since all experiments are limited to finite observation times, such a local measurement allows us to clearly identify any longer relaxation time scales emerging because of the disorder. Furthermore, the dynamical time evolution of the imbalance could capture eventual microscopic Griffiths-type effects, even in higher dimensions, where mass transport might not be sensitive to them [23]. The system is initialized in a striped density-wave pattern of fermionic 40 K atoms in a random mixture of two spin states (red and blue) in a square lattice with tunneling matrix element J, quasiperiodic potential of strength Δ, and tunable on-site interactions U between the different spins. The largest realized 2D system is composed of approximately 200 × 100 sites with several thousand atoms. (b) For weak disorder strength, the system thermalizes quickly (green area), whereas at strong disorder it is likely to exhibit a many-body localized regime (blue). Close to the transition (red dot, Δ c ), a regime of slow relaxation is observed (red area), potentially caused by locally insulating regions.

B. Identifying slow relaxation
for varying disorder strengths Δ; see Fig. 2. In the initial state, almost all the atoms occupy even stripes, such that the imbalance at zero evolution time is close to unity [see Fig. 2(a)]. For low disorder strength (Δ ¼ 1J), we observe a quick relaxation, and the imbalance vanishes within a few tunneling times. However, upon increasing the disorder, relaxation slows down dramatically (Δ ¼ 4J) and essentially comes to a full stop for strong disorder To quantitatively analyze this slow relaxation, the time dependence of imbalance is modeled as IðtÞ 1 IðtÞ × fðtÞ. Here,ĨðtÞ is the closed-system imbalance describing the dynamics of a perfectly isolated system, and fðtÞ represents a weak coupling to the environment. Such couplings are present in all real systems and will always thermalize any system at long enough times [40,41]. In our experiment, this weak coupling is dominated by a small but nonzero hopping rate between multiple two-dimensional planes along the z direction, with a rate J z ≈ J=10 3 [40]. We model the resulting imbalance relaxation due to this weak coupling with a stretched exponential fðtÞ ¼ e −ðΓtÞ β , with the decay rate Γ ¼ Γ exp ¼ 10 −3 τ −1 and the stretching exponent β ¼ 0.6 measured independently in a previous experiment [40].
The resultingĨðtÞ is shown in Fig. 3 (a) for short (10τ) and long (100τ) evolution times and fixed interaction strength U ¼ 5J as a function of the disorder strength. We can identify three dynamical regimes. For weak disorders Δ ≲ 2J (1), we observe vanishing values of short and long time imbalances, signaling the presence of a rapidly thermalizing system. Upon increasing the disorder strength, for 2J ≲ Δ ≲ 9J, we find a regime (2) characterized by a nonvanishing imbalance and clearly visible differences between the short and long time closedsystem imbalances. This indicates that, in this regime, the system relaxes much slower than the microscopic time scales. For Δ > 2J, all single-particle states are localized, but in many regions of the system, interactions with nearby atoms can still result in local thermal equilibrium. However, in some rare regions with anomalously low density or large spin imbalances (see below), this thermalization mechanism could be largely ineffective. Such regions can be thermalized only by their greater surroundings, which are thermal, but to which they couple in a significantly weaker fashion. Thus, the overall thermalization time scale grows. As the disorder is increased, these surroundings themselves gradually become more localized and less effective thermal baths. For strong disorder Δ ≳ 9J, we identify regime (3), where the values ofĨðtÞ at short and long times are both large and, within the experimental uncertainty, identical. This is consistent with the system being many-body localized.

C. Relaxation exponents and noninteracting inclusions
Identifying a suitable model to analyze the slow relaxation in regime (2) is challenging, as the underlying dynamics in two dimensions at intermediate disorder is theoretically unknown. In one-dimensional models with random potentials, anomalously strongly disordered regions have been argued to give rise to a subdiffusive phase via Griffiths effects [18,20,21,23]. Because our system contains quasiperiodic rather than random potentials, it should not contain such anomalously disordered regions. One possible mechanism explaining the slow relaxation at the observed time scales could be the randomness in the initial state, which can contain rare configurations, where some regions or inclusions are effectively noninteracting and therefore might appear localized at intermediate times, provided that the single-particle states are localized (Δ > 2J). In particular, thermalizing collisions cannot occur in a region where all atoms have the same spin. These regions can relax on longer time scales by coupling to the rest of the system, provided the disorder is not too high, by two possible mechanisms: First, particles inside the inclusion can tunnel out of the inclusion, at a rate that is exponentially small in their distance from the inclusion's edge, and then become thermalized by the thermal surroundings. Second, the spin and density imbalances can relax through slow nonlinear diffusion processes so that the inclusion evaporates from outside. While this could dominate at late times, we do not expect these diffusion processes to becomes relevant on the probed time scales [36].
We define inclusions of size L or larger as those containing a site from which any path of L steps encounters only atoms of one spin species or vacancies. The time to thermalize such a region through tunneling increases exponentially with L, as tðLÞ ∼ τ Ã e 2L=ξ , where ξ is the single-particle localization length of the noninteracting Aubry-André model. The factor of 2 in the exponent results from squaring the matrix element in Fermi's golden rule. A simple perturbative estimate of the prefactor τ Ã suggests that τ Ã ∼ max½U=J 2 ; ðΔ=JÞ 2 =U. The density of such inclusions is exponentially small in their volume, i.e., nðLÞ ∼ p L d in d dimensions, where p is the probability of a given site belonging to such an inclusion. Combining these expressions for nðLÞ and tðLÞ, we find that the density of inclusions that have not relaxed at time t goes as e −ζ log d ðtÞ [23], where ζ is the relaxation exponent quantifying the relaxation. While this ansatz results in a power-law relaxation in one dimension [24], in two dimensions, d ¼ 2, it results in a slightly faster than power-law relaxation of lognormal form:Ĩ ðtÞ ∼ e −ζ log 2 ðtÞ ∼ t −ζ logðtÞ : We find that this model is consistent with the experimental time traces in Fig. 2(b) (solid lines) in the entire intermediate regime (2). Because of the limited dynamical range of the data, we cannot unambiguously rule out other functional forms of the decay.
In particular, we cannot unambiguously distinguish between a power-law and a log-normal relaxation, the two main possibilities: A scenario giving rise to power laws is outlined and analyzed in Ref. [36]. A statistical χ 2 analysis slightly favors the log-normal form of Eq. (2) [36]. Therefore, we have employed it here to further analyze the data.
We find the extracted relaxation exponents ζ to continuously decrease with increasing disorder strength, as shown in Fig. 3(b), demonstrating that the system takes increasingly longer to relax. Furthermore, the fitted exponent ζ appears to vanish completely beyond a critical disorder Δ c . A simple power-law fit near the critical region ζ ∝ jΔ − Δ c j ν for Δ < Δ c and 0 otherwise yields a critical disorder strength Δ c ¼ ð9 AE 0.5ÞJ and a critical exponent ν ≈ 0.9. As the size of the critical region where scaling is expected to hold is unknown, the uncertainty in extracting such a critical exponent is also unknown. (2) decreases continuously with increasing disorder until a threshold value Δ c , at which it vanishes (green points). A piecewise fit ζ ∝ fjΔ − Δ c j ν ; 0g for ζ in the regime Δ ∈ f5J; 12Jg yields Δ c ¼ 9J AE 0.5J as the critical disorder strength of the possible MBL transition. The yellow line is an approximate upper bound for the relaxation exponent. The inset shows the same analysis for IðtÞ, i.e., neglecting the weak coupling to the environment. Apart from a small offset (dashed line) consistent with the known background coupling [40], we still find a sharp change of the extracted exponents at Δ c . Error bars denote the fit error and are typically smaller than the symbol size.
If we interpret a vanishing relaxation exponent as a signal for the MBL transition (as in one dimension [27]), our observations would support a transition into a MBL phase in two dimensions in the isolated limit. The extracted critical disorder strength of Δ 2D c ≈ 9J is significantly larger than the corresponding critical disorder strength in one dimension Δ 1D c ≈ 4J [27], even though the noninteracting critical disorder Δ U¼0 c ¼ 2J is identical in both cases. To demonstrate that these main results are essentially independent of our description for the environmental coupling, a corresponding analysis of IðtÞ, i.e., ignoring the weak coupling to the environment [27,36], is shown in the inset of Fig. 3(b). The resulting relaxation exponents also continuously decrease and exhibit a sharp kink before becoming constant at the same critical disorder strength Δ c , thereby highlighting the robustness of our result. The finite offset is consistent with the known coupling to the environment (dashed line) [27,40]. We emphasize that we find the same behavior near the extracted critical disorder strength Δ c also for power-law fits [36]. Therefore, the obtained results do not crucially depend on the choice of the fit function but rather illustrate that the dramatic slowing-down of relaxation dynamics on approaching the MBL transition from below is a generic feature of the two-dimensional quasiperiodic models.

D. Estimating the relaxation exponents
We can employ the above-mentioned model of counting the expected noninteracting inclusions to obtain a simple theoretical estimate of the exponents ζ. To this end, we assume that the surroundings of the inclusions act as a good thermal bath. To be specific, we focus on the central part of the system, in which even rows essentially reach unit filling. Here, a site in an occupied row has probability 1=2 of hosting, say, a spin-up atom; meanwhile, sites in empty rows are automatically part of the "inclusion." Accounting for the geometry of the inclusion, we estimate the exponent ζ th ¼ ðξ=2Þ 2 logð2Þ [36]. Note that this estimate does not account for imperfections of the initial state preparation [36], the inhomogeneity of the experimental system, or for the possibility that an inclusion can contain some small density of particles in the other spin state without causing it to thermalize. This theoretical estimate of ζ [as indicated by the yellow line in Fig. 3(b)] could thus be regarded as an upper bound on the relaxation exponent. It carries qualitatively similar features as the experimental data and, except close to Δ c , is also of the same order of magnitude as the experimentally extracted exponent. Near the transition, the experimentally measured relaxation exponent becomes significantly smaller than our theoretical estimate. In this regime, the typical surroundings of an inclusion themselves become increasingly localized and thus cease to act as a good bath, which invalidates the assumptions of our theoretical model. Furthermore, we estimate this simple model to become invalid for Δ ≲ 3J, where other relaxation mechanisms can become dominant already at the experimental time scales [36].

E. Interaction effects and energy density dependence
To highlight that the observed slow relaxation at intermediate disorders is driven by interactions and is completely absent in the noninteracting system, we compare it to the case of vanishing interactions U ≈ 0 in Fig. 4(a). We find indeed that the noninteracting system is strongly localized, showing a saturation of the imbalance at a finite value. In contrast, the interacting system relaxes slowly. This supports the fact that the slow relaxation is completely interaction driven.
To probe the impact of changing the energy density on the relaxation, we additionally prepare initial states with FIG. 4. Interaction and energy density dependence. (a) Time evolution of the imbalance for fixed disorder strength Δ ¼ 5J and vanishing U ≈ 0 and finite interactions U ¼ 5J. While there is negligible relaxation in the noninteracting case, the interacting time trace shows slow relaxation. This indicates that the relaxation is inherently due to interactions. We also show the normalized total atom number N to show that atom loss is minimal on the time scales of the experiment [40]. (b) We measure the imbalance for fixed disorder strength Δ ¼ 6J but starting from initial states with two different doublon fractions D ≈ 5%, 30%. In the large interaction limit, we observe a significantly higher imbalance for the larger doublon fraction. This can be qualitatively understood by the reduced mobility of a bound doublon. Each point is averaged over six individual experimental realizations, and error bars denote the error of the mean. Solid lines are guides to the eye. finite fraction D of particles on doubly occupied lattice sites (doublons), which changes the energy density of the initial state by approximately D × U=2. Figure 4(b) shows the resulting imbalance after an evolution time of 50τ. Because of the dynamical symmetry of the Hamiltonian, we expect the evolution to be the same for both repulsive and attractive interactions [10,38]. Hence, we concentrate on the attractive side of the Feshbach resonance [42] to be able to access very strong interactions. For small and moderate interactions (jUj ≲ 10J), we measure a strong decrease in the imbalance, highlighting again the delocalizing effect of interactions. Additionally, the imbalance is identical for both energy densities. For stronger interactions, however, we find a higher imbalance for the state with more doublons, indicating a stronger localization. Qualitatively, this effect can be understood by considering the reduced hopping rate of doublons J D ∼ J 2 =U at large interactions, which should result in stronger localization [10]. In one dimension, the hard-core limit of infinitely strong interactions in the absence of doublons can be mapped back onto a noninteracting model [10,43]. Consequently, the imbalance for a low doublon fraction is identical for the two extreme cases of vanishing and hard-core interactions [10]. In contrast, such a mapping does not exist in two dimensions and, accordingly, the imbalance is significantly different for the two extremes, again a striking qualitative difference in comparison to one dimension.

F. Conclusion
We observe an extended regime of exceedingly slow relaxation of an imprinted density pattern, in which the relaxation becomes progressively slower for increasing disorder strength. After accounting for a known weak coupling to the environment, the relaxation vanishes above a critical disorder strength, thereby indicating the existence of a MBL phase in two-dimensional quasiperiodic systems. We describe a simple model based on configurational inhomogeneities in high-temperature states that captures the qualitative trend of the experiment. However, a full quantitative description of the regime of critically slow relaxation and the apparent MBL transition goes substantially beyond the currently known theory, thereby underlining the importance of experimental results in resolving such regimes. Already in one dimension, recent numerical calculations have supported the fact that subdiffusion can arise, at least to intermediate times, even in systems with no rare regions in the underlying potential, such as quasiperiodic systems [44,45], and have shown that the emerging relaxation can appear very similar to the case of systems with a truly disordered potential [24,35]. Another possibility in two dimensions would be an intermediate critical phase between the fully ergodic and the fully MBL phase, which could correspond to this extremely broad regime with a markedly slow relaxation.
While our results already provide important insights into MBL in higher dimensions, revealing the universal critical properties of the transition remains a challenging task for future experiments. This would require identifying the regime where critical scaling holds and understanding the role of rare regions in quasiperiodic vs truly random systems in determining the critical properties [35,45,46]. Further isolating the system would allow us to experimentally access even longer time scales and understand the coupling to an external bath [27,40,41], as well as to probe the interplay between the spin and charge sectors arising from SU(2) symmetry [47][48][49], along with plaquette resonances in two dimensions [50]. Supplementing such results with frequency-resolved spectroscopy should further provide fundamental insights into the MBL transition [51].

ACKNOWLEDGMENTS
We are grateful to David Huse and Soonwon Choi for multiple stimulating discussions. We also thank Jae-yoon Choi