Probing many-body localization in the presence of a quantum bath

Closed generic quantum many-body systems may fail to thermalize under certain conditions even after long times, a phenomenon called many-body localization (MBL). Numerous studies have left no doubt about the stability of the MBL phase in strongly disordered one-dimensional systems. However, the situation is much less clear when a small part of the system is ergodic, a scenario which also has important implications for the existence of many-body localization in higher dimensions. Here we address this question experimentally using a large-scale quantum simulator of ultracold bosons in a two-dimensional optical lattice. We prepare two-component mixtures of varying relative population and implement a disorder potential which is only experienced by one of the components. The second non-disordered"clean"component plays the role of a small quantum bath of adjustable size that is collisionally coupled to the"dirty"component. Our experiments provide evidence that MBL is remarkably stable in 2D and can even survive in the presence of a small bath for our observed long timescales. A too large bath, however, ultimately leads to the destruction of MBL.

Typical quantum many-body systems evolve into a locally thermal state after driven out of equilibrium by a global quench [1]. This quantum version of thermalization is explained by the eigenstate thermalization hypothesis, which postulates that small subsystems are described by a thermal density matrix even for individual many-body eigenstates of the global system [2][3][4][5]. Quantum thermalization can, however, fail generically in systems exhibiting quenched disorder [6][7][8] when the strength of the disorder is large enough to prevent efficient spreading of entanglement [9]. Non-thermalizing behaviour and strong indication for the existence of a many-body-localized phase have been observed experimentally in several systems, in one [10][11][12][13][14] as well as in two dimensions [15,16]. The stability of MBL is, however, challenged in certain models and parameter regimes by the emergence of partially thermalizing features. In Fock space, eigenstates in parts of the energy spectrum may obey the eigenstate thermalization hypothesis while other parts of the spectrum remain localized [17][18][19]. In real space, rare regions of low disorder may form locally ergodic inclusions [20,21] that destabilize MBL. A system which shares many similarities with these two situations and that can provide insight into such MBLdelocalizing phenomena is an interacting two-component mixture composed of a "dirty" component in a random potential, and a "clean" component insensitive to the disorder [22,23]. In such a hybrid system, the clean component, which on its own would quantum thermalize, can be viewed as a bath with few degrees of freedom [24][25][26], a very different scenario compared to the theoretically [27][28][29][30] and experimentally studied [31] coupling of an MBL system to a classical bath at infinite coupling bandwidth. From the coupling of the two components, three possi-ble outcomes might emerge: 1) the full thermalization of the system, when the clean component acts as a "good" quantum bath which delocalizes the dirty component. 2) localization survives in the dirty component, while the clean one is localized as well through intercomponent interactions (an MBL proximity effect [22,23]) or 3) coexistence of a localized phase in the disordered component and an ergodic phase in the clean one. Here we report on the experimental study of such a two-component system. By tuning the relative population between the clean and dirty component we observe a qualitative change in behaviour, from localized for few clean atoms, to apparently ergodic for many. These results are remarkable in both regimes, since on the one hand they reveal a surprising robustness of MBL in two dimensions, while on the other hand the delocalization of the full system due to mere collisions with the clean component and no additional energy scales demonstrates that localization is present in a highly non-trivial regime.
Our experimental system consists of a two-component bosonic mixture of rubidium-87 atoms in a twodimensional optical lattice with lattice spacing a lat = 532 nm. The mixture is formed by populating the |d = |F = 2, m F = −2 and |c = |F = 1, m F = −1 hyperfine states, which feature almost equal inter-and intraspecies interactions U dc = U cc = U dd ≡ U [32]. We optically induce a state-dependent on-site disorder potential δ i , which only affects the dirty |d component, thereby breaking SU(2) symmetry (see Supplementary Information [33] whereâ i,σ ,â † i,σ andn i,σ denote the annihilation, creation and number operators for a particle in state σ ∈ {c, d} at site i of our 2D system [i = (i x , i y )]. The first sum indicates the hopping between nearest-neighbour sites i, j with a state-independent tunneling amplitude J, and V i characterizes the harmonic trapping potential.
To prepare an out-of-equilibrium initial state, we start with a unity-filling Mott insulator in the atomic limit and remove the atoms on every second column such that N = 124(12) atoms remain (see Fig. 1b). We then prepare a certain fraction η of the atoms in state |c via a microwave pulse before imposing the disorder potential and quenching the optical-lattice tunneling down to J/ = 2π × 24.8 Hz (see Fig. 1 and Supplementary Information for details [33]). The random disorder potential is chosen to be different for each experimental realization and its distribution is approximately Gaussian with a full width at half maximum of ∆ = 28 J [15]. The interactions amount to U = 24.4 J and after the lattice quench we allow the system to evolve for up to t 1000 τ , where τ = /J is the characteristic tunneling time. After the evolution, the lattices are increased to their maximum depth, in order to freeze the spatial distribution and image the atomic occupations on each individual lattice site [34].
We first focus on the dynamics of a system composed only of the dirty component for which we tracked the persistence of the initially prepared density distribution. This is measured by the imbalance I = Ne−No Ne+No , where N e,o are the occupation numbers on sites at even and odd columns, which is a probe for localization at short length scales. When only the dirty component is present, we observe the initial imbalance decaying over a few hundreds of tunneling times until a quasi-steady state is reached with a finite imbalance I d ≈ 0.13 -a signature of MBL. The decay of the imbalance towards this steady state can be well described by two exponentials with vastly different time constants. We identify the initial time scale, with a characteristic decay time of 0.6(2) τ , in which interactions are negligible since the atoms expand into empty sites. This initial dynamics is followed by a much slower relaxation with a decay time of 105(5) τ . Additional measurements of the number of doubly-occupied This picture is also supported by small-system numerical simulations (see Supplementary Information [33]). In addition to the two described time scales, the quasi-steadystate observed for long times also exhibits a very slow decay with a characteristic decay time of at least 2500 τ . We attribute this residual decay to the finite coupling of the system to the environment, as we observe a decay of the total number of particles in our system on a similar time scale (see Supplementary Information [33]).
We now turn to the main objective of this work, the study of a two-component system. In a first experiment we measured the dynamics of the dirty-component imbalance I d in the presence of a variable-size clean component. To ensure that the detection is only sensitive to atoms in state |d , we removed all |c atoms by a resonant light pulse prior to detection (see Supplementary Information [33]). The results are shown in Fig. 3, where we In the case of a small quantum bath (η = 0.15), a finite imbalance remains in a quasi-steady state, which is -up to an offset -indistinguishable from the η = 0 case. The inset shows the steady-state imbalance (measured after t = 625 τ ) as a function of the clean fraction η, which indicates a qualitative change from a localized to a delocalized system with increasing η. The horizontal dashed gray line indicates the typical statistical threshold at which the imbalance is compatible with zero. The error bars indicate one standard deviation of the mean.
included the η = 0 case as reference. Surprisingly, up to a global reduction of the imbalance, the long-time dynamics is similar for the η = 0 and the η = 0.15 case, where 15% of the atoms are in the clean state |c . We thus conclude that MBL persists in the presence of this small quantum bath on our observed timescales. The stability of MBL in contact with a weak bath has been studied recently theoretically [22,23]. In the simulation of two coupled clean and dirty 1D systems of equal population but different mobilities, localization was found to persist when the hopping strength of the clean component is reduced below a threshold value. Such a reduction on the hopping effectively decreases the rate at which it can couple spatially distant points of the dirty component. A decrease in population of the clean-component population η, as performed in our experiment, should also lower the effective strength of such a mediated coupling, leading to a similar outcome. When increasing η in the experiment, the imbalance I d keeps decreasing, until it vanishes at around η ≈ 0.4. In this situation, the localization length has grown well beyond a few sites, such that the initial short-scale density pattern becomes invisible. Notably, the disappearance of the imbalance is caused only by collisional interactions with the bath, whose en-ergy scale is comparable to the one set by the intracomponent interactions in the system. This underlines that the localization observed for lower η is not trivial at all.
The observed delocalization of the dirty component due to the bath coupling reflects that interactions generally tend to drive a system towards thermalization. Nonetheless, the back action exerted by the localized component may not be neglected in our experiment, and could instead lead to a localization of the bath. To probe the dynamics of the quantum bath in our system, we performed experiments in which we removed the dirty component |d prior to detection, in order to detect the clean component |c only (Fig. 4). The imbalance of these atoms I c relaxes very fast, on the time scale of a few tunneling times τ , and is seemingly independent of the clean-component fraction η. The vanishing I c after few tens of tunneling times τ does not, however, imply the thermalization of the clean component. The localization length of the bath might be too long to be detectable by our short-distance probing of the density wave. This argument is supported by the small remaining imbalance of the dirty component I d = 0.07(2) in the relevant η = 0.15 case, indicating localization on long length scales. Additionally, the occurrence of an MBL proximity effect requires an interaction-induced effective disorder, which, even for strong localization of the dirty component, would be inhibited due to the lack of randomness in the initially prepared state. We can therefore not rule out a localization of the bath by this measurement.
The experiments reported in this work shed first light on MBL systems in contact with a quantum bath of variable size. The results indicate a surprising robustness of localization for small baths, even in two dimensions. Follow-up experiments can explore different initial states and disorder regimes to settle the question of proximityinduced localization [22,23], demonstrating the localization of a system due to and not despite of interactions. Furthermore, the debated question on the stability of MBL in the presence of thermal inclusions can be directly addressed in a similar experiment with engineered low-disorder regions [20,21].

Initial state preparation
The experiment began with the preparation of a unityfilling rubidium-87 Mott insulator in a two-dimensional square optical lattice with a lattice spacing of a lat = 532 nm [34]. The typically prepared system consists of 250 atoms, and we keep the depth of the two lattices in the atomic plane at 40 E r , where E r = h 2 /8ma 2 lat is the recoil energy, h the Planck constant and m the atomic mass. At such a lattice depth, virtually no tunneling takes place, and thus the prepared state is separable.
We proceeded by removing all the atoms on the odd lattice sites along the x axis, thereby preparing a chargedensity wave (CDW) along one direction (see schematic in Fig. 1 of the main text). To do so, we use a spatially modulated laser beam with σ − polarization and a wavelength of 787.55 nm, which induces a differential lightshift of h × 10 kHz between the two hyperfine states |c = |F = 1, m F = −1 and |d = |F = 2, m F = −2 [35]. We then apply a microwave sweep to transfer the illuminated atoms to the hyperfine state |d and remove them by shining resonant light on the cycling transition of the D 2 line. The initial state fidelity is characterized by the imbalance between the even N e and odd N o lattice sites I = (N e − N o )/(N e + N o ) = 0.91 (1), and the remaining atom number is N = 124 (12).
To prepare an admixture of the two hyperfine states |d and |c , we afterwards apply a resonant microwave pulse of a certain length, which generates a state |Ψ i = √ 1 − η |d i + √ η |c i in each lattice site. After less than one tunneling time, inhomogeneities due to the disorder potential will lead to dephasing between the two spin states, and hence the whole system can be treated as a statistical spin mixture with a fraction η in the |c state.

Disorder potential
After the preparation of a CDW, we quenched the system by ramping up a projected disorder potential and lowering the depth of the in-plane lattices from 40 E r to 12 E r in less than 5 ms. The disorder potential is generated by the spatial modulation of a laser beam using a digital micromirror device (DMD), such that each lattice site of our 2D system features an individually programmable light shift [15]. The DMD consists of a 1024 × 768 micromirror array with a 13.7 µm micromirror pitch, and approximately 7 × 7 of these mirrors oversample the point spread function (Gaussian with σ = 0.48(1) a lat ) of our system. To image the disorder onto the atoms we use a high-resolution objective with a numerical aperture of NA = 0.68 [34]. A pseudorandom number generator is used to produce a 2D random pattern, which is taken different for every experimental realization. We resolved the microwave resonances spatially, which is equivalent to extracting the local light shift [15]. The histogram of all the local shifts displays a Gaussian distribution with full width at half maximum ∆, which characterizes the strength of the disorder. Our imaging system introduces a finite correlation to the disordered potential due to its finite resolution, for which we measure a correlation length of 0.63(1) a lat [15].
The disorder beam is tuned to the so-called 'tuneout' wavelength of the |F = 1, m F = −1 state, such that species-dependent potentials can be tailored [36]. This allows for tuning to a configuration where the |c component experiences a vanishing light shift, while the light shift for the |d state leads to an attractive potential. Aside from the programmed on-site disorder potential δ i , all atoms are equally sensitive to an overall harmonic trapping potential V i = ma 2 lat (ω 2 x x 2 i + ω 2 y y 2 i )/2 with frequencies (ω x , ω y ) = 2π × (51(2), 55(2))Hz.

Measuring the occupation number
At the end of each run we measure the atomic occupation in each lattice site. We first freeze the tunneling motion by increasing the lattice depth to 60 E r in less than half a millisecond. To selectively measure only one of the two hyperfine components, we then remove all the atoms in the state that we do not wish to measure. To image the |c state, we push the |d -state atoms using a resonant D 2 light pulse, while for detecting the |d state, we apply a microwave sweep to swap the two hyperfine spin states before the optical push-out pulse, thus removing the atoms which were originally in the |c state.
After the state selection, all lattices are ramped to their maximum depth and an optical molasses is used to scatter fluorescence photons and simultaneously cool the atoms. We expose an EMCCD camera for 1 s, thereby obtaining single-site-resolved images of the parity-projected density distribution [34].

Atom loss and doublon generation
An estimate of the total atom number in the system cannot be directly obtained from the in-situ fluorescence imaging, since parity projection prevents us from distinguishing empty lattice sites from doubly-occupied ones. When discarding local information we can circumvent this issue. Performing a short time of flight in the 2D plane right before imaging dilutes the atom density, such that essentially no doubly-occupied sites remain in the final atomic distribution. We used this technique to measure the atom number of an initially prepared Mott insulator by turning both in-plane lattices off for 1 ms. By measuring the total atom number N (t) for different holding times, we were able to estimate the atom losses. This is shown in fig. S1, in which the purple points represent the normalized atom number, n(t) = N (t)/N (0), and the purple line is a linear fit with an atom number decay of 3300(300) τ . This means that after 600 τ , approximately 15 % of the initially prepared atoms have been lost, for which the main loss mechanism is induced by technical fluctuations of the lattice-beam intensities [15]. We estimate the photon scattering rate from the disordered potential to be below γ = 7 · 10 −5 τ −1 . Though our atom loss is comparable to previous work in other quantum-gas experiments [31], the measured imbalance decay seems little affected by it.
By comparing the parity-projection-free atom number with the one obtained from direct in-situ measure- ments of the CDW relaxation ( fig. 2 in the main text), we can additionally estimate the fraction of dynamically generated doublons even for long times. We define the doublon fraction as p d (t) = 2·N dou (t) , where N dou (t) is the number of doublons, and N CDW (t) and n CDW (t) are respectively the absolute and normalized parity-projected atom number. For the data measured at ∆ = 28 J, we find that the doublon fraction remains approximately constant after the fast initial formation (see Fig. S1). We did not take into account the effect of triply-occupied sites, which at the current experimental settings is estimated to be negligible.

Numerical simulations
To get more insight on the slow relaxation observed for the dirty-component imbalance dynamics (Fig. 2 in the main text) we have performed numerical simulations based on exact diagonalization for a small disordered-Bose-Hubbard system, for which we used the Quspin package [37]. The considered system is a ladder of 2 × 6 lattice sites with periodic boundary conditions, populated with 5 particles in a CDW-like pattern. We choose the system parameters close to the experimental ones, with a disorder distribution given by a Gaussian with full width at half maximum of ∆ = 25 J and we have considered both a non-interacting (U = 0) and an interacting case (U = 25 J).
Without interactions, the imbalance dynamics shows a very fast relaxation to I = 0.7, which takes less than 10 τ . In contrast, for the interacting system a steady state is only reached after almost 100 τ , and one can distinguish a very fast initial decay to I = 0.6 in few tunneling times from a much slower relaxation afterwards, quite similar to our experimental findings. Interestingly, such a slow relaxation is not observed in one-dimensional numerical simulations with the same interaction strength, not even for an intermediate disorder strength ∆. The slow decay stops at a steady-state imbalance of I ≈ 0.5, which is notably higher than the one measured in our experiment (I ≈ 0.15). Such a discrepancy can be explained by the small number of atoms in the simulation, which will diminish the effect of interactions in comparison to the real system. We have also obtained the doublon-fraction dynamics, which clearly indicate a fast doublon formation during the initial time scale, followed by a much slower increase, also in a similar fashion to the experimental results.