Research ArticlesPhysics

Exploring the many-body localization transition in two dimensions

See allHide authors and affiliations

Science  24 Jun 2016:
Vol. 352, Issue 6293, pp. 1547-1552
DOI: 10.1126/science.aaf8834

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 87Rb 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


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 (68). Furthermore, the breakdown of the eigenstate thermalization hypothesis (912) 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 (1620). 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 (2530) or interacting ultracold atoms in lower-energy states (3136). 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 87Rb atoms in a single plane of a cubic optical lattice with lattice spacing alat = 532 nm. To ensure that the initial state was separable, we froze out any motion at a lattice depth of 40Er in the x and y directions, where Embedded Image 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.6alat 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 12Er 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).

Fig. 1 Schematics of the experiment and raw images.

(A) A 2D random disorder potential is imaged onto a single atomic plane in an optical lattice. The disorder is controlled by a digital mirror device (DMD), which converts a Gaussian laser intensity profile into a 2D random intensity distribution with spatially uniform mean light intensity (bottom image). The limited numerical aperture (NA = 0.68) of the microscope objective introduces a finite correlation length and leads to a smoothing of the disorder distribution. The histogram at bottom right (red bars) is the measured disorder distribution and its asymmetric Gaussian fit curve (red solid line), where Δ is the full width at half maximum of the disorder distribution. Distinct to the other two images showing the original (bottom) and smoothed (middle) light intensity distributions, the top image displays the local disorder potential determined by in situ spectroscopy (42). The yellow circles on the lower images indicate the spectroscopically calibrated region. (B) Raw fluorescence images (the red-to-yellow color scale corresponds to increasing detected light level) showing the evolution of the initial density step without disorder. The left column shows single images (isolated red dots are individual atoms) of the parity-projected atomic distribution for the indicated evolution times. The right column displays the mean density distribution averaged over Embedded Image different disorder potentials. The top left image depicts the initial state for which the analysis region (dx × dy = 5 × 31) is indicated by the white box. For the high-disorder case shown in (C) the detected initial-state filling is slightly lower, which is an artifact of the parity projection (42). In contrast to (B), traces of the initial state remain at all times in the disordered case. The white circles in the averaged density profiles after t = 249τ highlight the difference. a.u., arbitrary units.

During the dynamics, the system is described by a 2D Bose-Hubbard Hamiltonian with onsite disorderEmbedded ImageHere, Embedded Image Embedded Image is the bosonic creation (annihilation) operator, Embedded Image is the local density operator on site Embedded Image, and the first sum includes all neighboring sites. A harmonic trapping potential Embedded Image with frequencies Embedded ImageHz in the x and y direction confines the atoms around the trap minimum. The nearest-neighbor hopping strength at a lattice depth of 12Er 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.4J, 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 12Er (that is, assuming decoupled sites), we extracted an upper limit of the temperature of T < 0.54(1)U/kB, where kB 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 Er/N = 0.58(1)U might be overestimated by up to 0.15U, due to the finite band width, and agrees reasonably well with the expectation for a thermalized state, assuming that the lattice ramp to 12Er was adiabatic with respect to the band gap. Here, the energy density of the initial out-of-equilibrium state contributes with E0/N = 0.28(3)U (where E0 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 [EH/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 (NL) and right (NR) atom number imbalance Embedded Image, 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 Embedded Image. 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 ts, 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 ~Embedded Image within one relaxation time ts. Small couplings with the environment might lead to relaxation of the quasi–steady state on time scales much longer than our experimental time scale (4650).

Fig. 2 Relaxation dynamics of a density domain wall.

The evolution of the imbalance Embedded Image shown for five different disorder strengths [Δ/J = 0 (dark green), 3 (medium green), 4 (light green), 8 (light blue), and 13 (dark blue)] displays a saturation behavior toward a quasi–steady state for all disorder strengths. For low disorders (green curves), the asymptotic value of the imbalance is vanishing, whereas a finite imbalance remains for higher disorder (blue). Solid curves are fits to the data with Embedded Image, of which the decay time ts is plotted versus disorder strength Δ in the inset. Error bars represent 1 SD of the mean in the main figure and Embedded Image confidence bounds of the fit parameters in the inset. ħ, Planck’s constant h divided by 2π.

Identifying the MBL transition

The transition from zero to nonzero imbalance Embedded Image 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 Embedded Image, which indicates that the transition is taking place in the system, where we extracted the critical disorder from a simple double linear fit Embedded Image with Embedded Image. In the vicinity of the transition, slow subdiffusive transport has been predicted (21, 22, 5153), 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 deviationEmbedded Imageof the vertically (y direction) averaged reference profile Embedded Image and the finite disorder profiles ni(Δ). We observe that the profiles start to deviate at Δcn = 5.3(2)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 8J are comparable.

Fig. 3 Identifying the MBL transition.

(A) Disorder dependence of the imbalance after equilibration of the dynamics [constant evolution time of 187τ, also in (B)]. The data show a sharp onset of nonzero quasi–steady-state imbalance Embedded Image at a disorder strength of Embedded Image, and the solid line is a double linear fit (described in the text) to extract the critical disorder. The inset illustrates the sharpening of the transition versus time and demonstrates the saturation of its shape by showing the imbalance extracted from the fits in Fig. 2 at different evolution times. (B) Deviation from the zero-disorder thermal profile, as measured by the root mean square density difference δn at various disorder strengths. The relaxed density profiles differ by more than the random position-induced measurement noise from the thermal profile abruptly above the critical disorder strength Δcn = 5.3(2)J. The inset shows the averaged reference density profile for zero disorder (green) and the averaged profile at a high disorder of Δ = 13J (blue). Error bars represent 1 SD of the mean.

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 ni(Δ) to the thermal profile ni(0) by calculatingEmbedded Image (Fig. 4A). The observed decay is found to be well captured by an exponential fitEmbedded Image 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 Δ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 Embedded Image. We emphasize, that all three methods of extracting the critical disorder strength of the MBL transition in the experiments agree within the fit errors.

Fig. 4 Diverging density decay length at the localization transition.

(A) The spatial dependence of the normalized average density n(Δ)/n(0) in the initially empty region is fitted by an exponentially decaying model (solid lines). The level of blue brightness encodes the disorder strength increasing from light to dark: Δ/J = 4, 7, 13, and 17. (B) Fitted decay constant λ as a function of disorder strength Δ. The solid light blue line is a double linear fit (described in the text) to locate the transition point Δcn = 5.3(2)J. The inset shows the diverging decay length ξ = 1/λ near the critical disorder strength. Evolution time 187τ is for all data shown here. Error bars represent 1 SD of the mean in (A) and Embedded Image confidence bounds of the fit parameters in (B).

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 (E0/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 Embedded Image displays a sharp transition behavior at a smaller critical disorder of Embedded Image. 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.

Fig. 5 Interaction dependence of the localization transition.

Quasi–steady-state imbalance Embedded Image versus disorder strength Δ for different initial densities (evolution time 187τ). The interaction effects are reduced by lowering the initial filling to 0.23, which is 25% of the value previously discussed in Fig. 3A (light blue). The clear difference in critical disorder strengths highlights the strong influence of interactions on the localization. The right two insets show representative fluorescence images of the initial density distribution for each case. The left inset is a zoomed-out view of the main figure where we added the results of exact diagonalization numerics for the noninteracting system with the same experimental conditions (black open squares). Here, the horizontal error bars denote the systematic uncertainty in the disorder strength. Vertical error bars represent 1 SD of the mean.


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, 5153). 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

Supplementary Text

Figs. S1 to S4

Reference (58)

References and Notes

  1. See supplementary materials on Science Online.
  2. 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.).
View Abstract

Stay Connected to Science

Navigate This Article