## Bosons refusing to thermalize in 2D

Messy, interacting quantum-mechanical systems are difficult to analyze theoretically. In a single spatial dimension, the calculations are still tractable, and experiments have recently confirmed the prediction that sufficiently strong disorder can disrupt the transport of interacting particles. In two dimensions, however, the theoretical blueprint is missing. Choi *et al.* used single-site imaging of cold ^{87}Rb atoms in an optical lattice to show that similar localization occurs in two-dimensional (2D) systems. The study highlights the power of quantum simulation to solve problems that are currently inaccessible to classical computing techniques.

*Science*, this issue p. 1547

## Abstract

A fundamental assumption in statistical physics is that generic closed quantum many-body systems thermalize under their own dynamics. Recently, the emergence of many-body localized systems has questioned this concept and challenged our understanding of the connection between statistical physics and quantum mechanics. Here we report on the observation of a many-body localization transition between thermal and localized phases for bosons in a two-dimensional disordered optical lattice. With our single-site–resolved measurements, we track the relaxation dynamics of an initially prepared out-of-equilibrium density pattern and find strong evidence for a diverging length scale when approaching the localization transition. Our experiments represent a demonstration and in-depth characterization of many-body localization in a regime not accessible with state-of-the-art simulations on classical computers.

In his seminal work on localization in quantum mechanical systems, Philip Anderson emphasized the implications of localization on the thermodynamics of closed quantum systems (*1*). Recently, perturbative arguments suggested the existence of nonthermalizing, many-body localized systems at low energy (*2*, *3*). Soon thereafter, these arguments were extended to all interaction strengths and energy densities for systems with a bounded spectrum (*4*, *5*). The implication—nothing less than a breakdown of equilibrium statistical mechanics for certain generic macroscopic systems—triggered tremendous theoretical efforts (*6*–*8*). Furthermore, the breakdown of the eigenstate thermalization hypothesis (*9*–*12*) caused by the failure of these systems to act as their own heat bath implies the persistence of initial state information, which might serve as a useful resource for quantum information technologies (*13*). Several other notable features of many-body localization (MBL) have been uncovered, such as the description of fully localized systems by coupled localized integrals of motion (*14*, *15*). This underlies the absence of particle transport but allows the transport of phase correlations, leading to a characteristic logarithmic growth of the entanglement entropy in the case of short-range interactions (*16*–*20*). Another distinctive feature of many-body localized systems, as compared with noninteracting low-dimensional systems, is the requirement of a nonzero disorder strength for the localized phase to appear (*21*, *22*).

Recently, the absence of thermalization due to MBL in a quasi-disordered one-dimensional (1D) Fermi lattice has been reported (*23*, *24*). These studies explored the system’s behavior at long times and high energy density, as opposed to earlier experiments with noninteracting systems (*25*–*30*) or interacting ultracold atoms in lower-energy states (*31*–*36*). A recent experiment with 3D disordered lattice fermions provided evidence for the absence of particle transport, even at elevated temperatures (*37*). Indications for localization in Fock space, one characteristic property of MBL (*2*), have been reported in short ion chains (*38*), and MBL has been suggested as one possible explanation for the recently observed vanishing conductance in disordered superconductors at nonzero temperature (*39*). However, despite intensive theoretical and experimental efforts, some aspects of MBL, such as the details of the localization transition, including the identification of diverging length scales, are still not fully understood. Whereas in one dimension the localization transition is rather well studied (*21*, *22*, *40*, *41*), the nature of MBL in higher dimensions is an open question.

Here we address the open question of the nature of a MBL transition in two dimensions, which we observe experimentally and characterize. We report on the single-site–resolved study of thermalization and transport in a disordered 2D bosonic optical lattice, starting from a high–energy density initial state far from equilibrium. By tracking the time evolution of an initially prepared density domain wall for variable disorder strengths, we reveal the fairly sharp onset of nonthermalizing behavior above a critical value of disorder strength. The observed localization transition is found when the disorder, single-particle bandwidth, onsite interaction, and energy density are all of comparable strength. This regime of parameters, being far from the disorder-dominated (classical) and the interaction-dominated limits, is highly nontrivial from a theoretical perspective. Precise and direct characterization of the projected disorder potential by site-resolved spectroscopy allows for a direct comparison of our results with a numerical prediction for the noninteracting model, highlighting the dramatic differences and pointing out the key role of interactions. Furthermore, by locally comparing the observed density profiles to a thermalizing reference measured without disorder, we obtain strong evidence for a diverging length scale when approaching the transition from the localized side.

## Experimental setup

We began by preparing a 2D approximately unity filled Mott insulator of bosonic ^{87}Rb atoms in a single plane of a cubic optical lattice with lattice spacing *a*_{lat} = 532 nm. To ensure that the initial state was separable, we froze out any motion at a lattice depth of 40*E*_{r} in the *x* and *y* directions, where is the recoil energy, *h* the Planck constant, and *m* is the atomic mass. Next, we used a digital mirror device–based spatial light modulator (*42*) to optically remove the right half of the atomic population such that about *N* = 125(11) atoms remained in the lattice sites located at *x* < 0, thus preparing a sharp density domain wall. For each experimental realization, a new computer-generated random disorder pattern drawn from a uniform distribution was displayed on the light modulator and subsequently projected onto the atoms (Fig. 1A). The optical projection results in a slightly asymmetric disorder distribution and a convolution of the disorder, with the point spread function of the imaging system leading to a finite disorder correlation length 0.6*a*_{lat} and a narrowing of the disorder distribution to a width Δ (Fig. 1A, inset) (*43*). We used single-site–resolved spectroscopy of the atomic sample to directly characterize the disorder. The topmost disorder picture in Fig. 1A displays the result of one such spectroscopic measurement for the configuration shown in the two images below. We then initiated the dynamics of the initial domain wall by lowering the *x* and *y* lattice depths to 12*E*_{r} in 5 ms, which is less than one tunneling time. Next, we allowed the system to evolve for a variable time *t* in the disorder potential, after which we used fluorescence imaging to measure the local parity-projected density (*44*).

During the dynamics, the system is described by a 2D Bose-Hubbard Hamiltonian with onsite disorderHere, is the bosonic creation (annihilation) operator, is the local density operator on site , and the first sum includes all neighboring sites. A harmonic trapping potential with frequencies Hz in the *x* and *y* direction confines the atoms around the trap minimum. The nearest-neighbor hopping strength at a lattice depth of 12*E*_{r} is *J*/*h* = 24.8 Hz, corresponding to a tunneling time of τ = *h*/2π*J* = 6.4 ms. Longer-range hopping terms are neglected, as they are exponentially suppressed. The onsite interaction strength is *U* = 24.4*J*, and δ** _{i}** denotes the onsite disorder potential. For these parameters, in the absence of disorder, the system’s ground state is in the Mott insulating phase, albeit with strong particle-hole fluctuations (

*45*).

## MBL in two dimensions

For reference, we first tracked the time evolution without any disorder potential applied. Even from the bare images (Fig. 1B), it is clear that the initially prepared density step is smeared out after a few tens of tunneling times τ, and after a longer time, no information about the initial density step remains. The observed density distribution appears thermal, and neglecting quantum fluctuations at 12*E*_{r} (that is, assuming decoupled sites), we extracted an upper limit of the temperature of *T* < 0.54(1)*U*/*k*_{B}, where *k*_{B} is the Boltzmann constant. This temperature estimate was obtained by a global fit to the radial density profile, assuming a grand canonical ensemble for each site with the chemical potential given by the local density approximation (*44*). The corresponding energy per particle of *E*_{r}/*N* = 0.58(1)*U* might be overestimated by up to 0.15*U*, due to the finite band width, and agrees reasonably well with the expectation for a thermalized state, assuming that the lattice ramp to 12*E*_{r} was adiabatic with respect to the band gap. Here, the energy density of the initial out-of-equilibrium state contributes with *E*_{0}/*N* = 0.28(3)*U* (where *E*_{0} is the sum of the initial thermal and the potential energy due to the trap), determined by the initial thermal energy and the harmonic trap with frequency ω* _{x}*. An additional energy increase is present due to small heating during the 2-s evolution time [

*E*

_{H}/

*N*= 0.18(6)

*U*] used in this measurement (

*43*). In contrast, upon repeating the measurement with strong disorder, traces of the initial state remain, and the system does not relax to a spatially symmetric density distribution expected for thermal state (Fig. 1C).

A direct and model-free quantity that can be used to identify a nonthermalized state is the density asymmetry quantified by a nonzero left (*N*_{L}) and right (*N*_{R}) atom number imbalance , which we analyze, as with all other extracted quantities, in a central region of interest extending over five lattice sites in the *y* direction. The zero line (*x* = 0) separating left and right is defined by the position of the initial density step and was precisely aligned to the closest lattice site to the trap center, resulting in an offset of up to ±1 lattice sites, corresponding to . The evolution of the imbalance (Fig. 2) confirms that, for all disorder strengths, the system reaches a quasi–steady state within ~150τ. For small disorder strengths, we find a vanishing imbalance, whereas for large disorder, a nonzero imbalance remains, even for long evolution times. We interpret this latter regime as the many-body localized phase, where the observed quasi–steady state is clearly nonthermal and transport through the system is blocked. The relaxation time *t _{s}*, extracted from an exponential fit to the data, increases with disorder strength and hints at a saturation in the nonthermal regime (Fig. 2, inset). We now turn to a series of measurements in which we fix the evolution time to ~190τ, which is well in the quasi–steady-state regime but short enough to keep the effects of atom number loss and noise-induced coupling to higher-energy bands negligible. On this time scale, we also expect the effects of low-frequency noise on the disordered system to be small. Considering the measured heating rate in the nondisordered system as an upper bound for the energy increase, the energy per particle would change by only ~ within one relaxation time

*t*. Small couplings with the environment might lead to relaxation of the quasi–steady state on time scales much longer than our experimental time scale (

_{s}*46*–

*50*).

## Identifying the MBL transition

The transition from zero to nonzero imbalance for large disorder indicates the presence of a thermalizing phase for low disorder strengths and an apparent transition to a localized phase at higher disorder. To locate an MBL transition, we recorded a series of measurements with fixed evolution time in the quasi–steady-state regime and scanned the disorder around the critical value (Fig. 3). We find a fairly sharp onset of nonzero–steady-state imbalance at a disorder strength , which indicates that the transition is taking place in the system, where we extracted the critical disorder from a simple double linear fit with . In the vicinity of the transition, slow subdiffusive transport has been predicted (*21*, *22*, *51*–*53*), which suggests that our measurements might not have reached a true steady state in this regime. Because of this, it is possible that the resulting transition might move toward stronger disorder if we were able to study much longer times. Our high-resolution detection allows for a local comparison of the measured density profiles with the equilibrated thermal profile observed at vanishing disorder. As shown in Fig. 3B, we use this method as a more sensitive probe to detect deviations from the thermal profile by calculating the root mean square deviationof the vertically (*y* direction) averaged reference profile and the finite disorder profiles *n _{i}*(Δ). We observe that the profiles start to deviate at Δ

_{c}_{,δ}

*= 5.3(2)*

_{n}*J*, with a fairly sharp kink signaling the transition, which is quantitatively consistent with the imbalance measurements. Whereas localization is predicted for vanishingly small disorder strength in the noninteracting case, the finite interaction

*U*in our system promotes thermalization. This is consistent with our observation of a thermal behavior at nonzero Δ. However, as the disorder strength is increased above a critical value, the localization is restored, which is particularly notable because this transition takes place in the regime where the disorder Δ, the interaction

*U*, and the single-particle bandwidth 8

*J*are comparable.

The MBL phase transition is expected to be a continuous transition (*21*, *22*, *40*, *41*, *54*) for which one expects a characteristic diverging length scale when approaching the transition. With our experiments, we have direct access to the steady-state decay length ξ(Δ) of the initially prepared density step, which is directly related to a density-density correlation length. To minimize the influence of the external trap, we reference the density profile *n _{i}*(Δ) to the thermal profile

*n*(0) by calculating (Fig. 4A). The observed decay is found to be well captured by an exponential fit with a disorder-dependent decay constant λ(Δ) = 1/ξ(Δ). A priori, we would have expected a crossover from an interaction-dominated decay to a single-particle–dominated decay at low densities present at the outer edge of our sample. However, within the experimental uncertainty we cannot identify such an effect. We directly observe the diverging behavior of the decay length ξ(Δ) when approaching the transition from the localized side (Fig. 4B, inset). The related decay constant λ(Δ) also shows a very sharp kink at Δ

_{i}

_{c}_{,λ}= 5.3(3)

*J*, marking the onset of the localized region. We empirically obtain the critical disorder Δ

_{c}_{,λ}using the bilinear fit function, λ(Δ) = λ

_{0}+ λ

_{1}∙ max[(Δ – Δ

_{c}_{,λ}),0] with . We emphasize, that all three methods of extracting the critical disorder strength of the MBL transition in the experiments agree within the fit errors.

## Interaction dependence

To experimentally verify that the observed behavior of the system is induced by interactions, we reduced the initial density and, therefore, the effects of the interaction on the dynamics while keeping the initial energy per particle (*E*_{0}/*N*) constant. Lower atomic densities are obtained by uniformly transferring a given fraction of the atoms by applying a microwave pulse to another hyperfine state and then optically removing the transferred population. When we reduced the density by factor of 4, we observed a clear shift of the localization transition toward lower disorder strength, as was expected (Fig. 5). Here, the left-right imbalance displays a sharp transition behavior at a smaller critical disorder of . Furthermore, the numerical prediction for a noninteracting system obtained by exact diagonalization (*43*) is fully incompatible and strongly differs from both measurements, showing that the interactions facilitate thermalization for low disorder strengths. Reducing the system size by the preparation of a smaller initial Mott insulator did not affect the observed critical disorder (*43*). From these observations, we conclude that the nonzero critical disorder strength is due to the interactions and that the measured critical disorder value is not strongly influenced by the small external driving caused by laser fluctuations discussed above.

## Conclusions

Our experiments provide evidence for MBL in two dimensions by the observation of a transition from a thermalizing phase to a localized phase of interacting bosons in a disordered optical lattice. The system size analyzed in our experiment is far beyond numerically accessible scales, demonstrating a nontrivial quantum realization of the MBL system that challenges both analytical and numerical theory. Furthermore, we supplemented our observation of an MBL transition with the demonstration of a clear shift of the transition point for effectively smaller interaction energy. Even though it is difficult to distinguish a true phase transition from a sharp crossover within experiments, our results mark a first step in understanding MBL in more than one dimension and can be extended to obtain detailed information about the nature of the MBL transition, such as its dynamical critical exponent (*21*, *22*). Furthermore, supplementing transport experiments with density-density correlation measurements, or even measurements of the entanglement entropy (*55*, *56*), offers a promising possibility to demonstrate ongoing phase dynamics in the localized phase while the density or charge transport is frozen (*38*, *57*). Via detailed studies of the dynamics of transport, it should also be possible to study Griffiths effects and subdiffusive transport in the vicinity of the transition (*21*, *22*, *51*–*53*). In addition to dynamical properties, our technique might allow the investigation of many-body eigenstate–related properties, such as the local integrals of motion (*14*, *15*) defined via local operators (*54*).

## Supplementary Materials

## References and Notes

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵See supplementary materials on
*Science*Online. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
**Acknowledgments:**We thank M. Fischer, E. Altman, M. Knap, P. Bordia, H. Lüschen, A. Rosch, E. Demler, and S. Sondhi for discussions. D.A.H. is the A. and H. Broitman Member at the Institute for Advanced Study. We acknowledge funding by Max-Planck-Gesellschaft, Deutsche Forschungsgemeinschaft, the European Union (UQUAM, Marie Curie Fellowship to J.C.).