Signatures of Many-Body Localization in a Controlled Open Quantum System

In the presence of disorder, an interacting closed quantum system can undergo many-body localization (MBL) and fail to thermalize. However, over long times even weak couplings to any thermal environment will necessarily thermalize the system and erase all signatures of MBL. This presents a challenge for experimental investigations of MBL, since no realistic system can ever be fully closed. In this work, we experimentally explore the thermalization dynamics of a localized system in the presence of controlled dissipation. Specifically, we find that photon scattering results in a stretched exponential decay of an initial density pattern with a rate that depends linearly on the scattering rate. We find that the resulting susceptibility increases significantly close to the phase transition point. In this regime, which is inaccessible to current numerical studies, we also find a strong dependence on interactions. Our work provides a basis for systematic studies of MBL in open systems and opens a route towards extrapolation of closed system properties from experiments.


I. INTRODUCTION
In a perfectly closed system, many-body localization (MBL) presents a novel paradigm of time evolution, in which quantum correlations can persist locally to arbitrarily long times [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. It thereby provides a robust alternative to conventional thermalizing dynamics. However, experiments are invariably coupled, at least weakly, to their environment, and thus cannot realize a strictly closed system. It is therefore crucial to understand how such dissipative couplings affect the many-body dynamics. Recently, experiments have observed MBL through the persistence of initially prepared density or spin patterns on intermediate to long timescales [13][14][15][16][17]. At even longer times, however, the patterns vanish and the systems become thermal due to residual couplings to the environment. These couplings not only set a timescale for thermalization, but are also expected to broaden the localization transition into a crossover ( Fig. 1(a)), in which the dynamics smoothly interpolate between those of an ergodic and an MBL system. Features of the critical point could then be encoded into a universal dependence of the relaxation curves on the dissipation rate. This is similar to the role of temperature in ground state quantum phase transitions, where measuring universal temperature dependencies allows for a characterization of the critical point [18]. Systematically varying the strength of the dissipative couplings promises the analogous possibility of extrapolating to closed systems.
Several theoretical works have recently addressed different aspects of MBL systems coupled to an external bath [8,[22][23][24][25][26][27][28][29][30][31][32][33][34]. In particular, recent analytical and numerical studies considered the relaxation of almost local integrals of motion associated with MBL under the influ-ence of a weak coupling to a photon bath [29][30][31][32][33]. This constitutes a Markovian heat bath at infinite temperature operating mainly through two dissipation channels ( Fig. 1(b)): (i) effectively measuring the position of the particles and (ii) particle loss. Position measurements affect the system by dephasing coherent superpositions of Wannier states, and therefore show a much stronger effect on systems with finite or longer range quantum coherences as compared to e.g. certain glasses that may behave classically already at the scale of the inter-particle distance.
Here, we make use of the exceptional microscopic understanding and control over dissipative processes that is possible with cold atoms in optical lattices [35,36] to explore MBL in an open quantum system. Specifically, we investigate how photon scattering affects the dynamics of a density pattern in the many-body localized regime. Deep in the localized phase, we find a stretched exponential relaxation of the density pattern in good agreement with theoretical studies [29][30][31][32][33]. In the experiment, we are furthermore able to study the intriguing regime close to the MBL transition that is not accessible to current numerical studies. There we find that dissipation has an increasingly strong effect that is significantly enhanced by inter-particle interactions.  (a) Coupling a disordered system to a thermal environment destroys MBL on a timescale inversely proportional to the coupling strength γ. Nonetheless, for sufficiently weak couplings the characteristics of both the MBL (i.e. persisting density pattern) and the ergodic Griffiths phase (i.e. power law decay of density pattern [19][20][21]) survive sufficiently long to enable their experimental characterization. The shaded area represents the regime of weak dissipation in which the intrinsic dynamics of the phases can be discerned. Due to diverging timescales at the critical point, the respective signatures become increasingly difficult to observe in its proximity. A sharp transition is only expected in the closed system limit (γ = 0), while finite dissipative couplings are expected to produce a crossover regime between the two characteristic dynamics, similar to the effect of temperature in ground state quantum phase transitions. Black arrows indicate the regime that is considered in this work. (b) Schematic of a scattering event. An atom in an initial localized superposition of Wannier states (blue) becomes localized on the length scale of the scattered photon's wavelength λ/2π, which, in our system, is less than the lattice constant d. This dephases superpositions to incoherent mixtures of single Wannier states in the ground band (yellow), as well as producing a small population in higher bands (faint yellow), which can be seen as the result of position measurements with sub lattice site resolution. Band excitations can then lead to atom loss, since in most excited bands atoms are not trapped.

II. EXPERIMENT
We start the experiment by cooling a gas of 130 × 10 3 40 K atoms in a dipole trap to a temperature of 0.15 T F , where T F denotes the Fermi temperature. The gas is then loaded into a three dimensional optical lattice consisting of two deep (λ o = 738.2 nm) orthogonal lattices, which create an array of one dimensional tubes, and a primary (λ p ≈ 532.2 nm) lattice with lattice constant d = λ p /2 along the tubes. We superimpose the primary lattice with an incommensurate (λ d ≈ 738.2 nm) disorder lattice to implement the interacting Aubry-André model [13,37] in the individual tubes. This model describes spinful fermions on a tight-binding lattice with on-site interaction U and nearest-neighbor tunneling amplitude J ≈ h · 500 Hz, subject to a quasi-periodic potential ∆ cos(2παi + φ) with amplitude ∆. Here, i ∈ Z numbers the lattice sites, α = λ p /λ d is the disorder periodicity and φ is the relative phase between the primary and the incommensurable disorder lattice. In the absence of interactions this model exhibits a localization transition at ∆ = 2 J [37], and it has been shown to be many-body localized at U = 0 above a parameter dependent critical disorder strength [13].
Using a period-two superlattice, we artificially create a charge-density wave state with an initial imbalance denotes the number of atoms on even (odd) sites. After the desired evolution time, we extract the remaining imbalance using a superlattice band-mapping technique [38]. For ergodic systems, the imbalance decays to zero during the evolution, while a persisting imbalance signals localization. Further details of the system, the preparation and the readout sequence can be found in references [13,14,38].
The atoms are prepared in an equal mixture of the two lowest spin states in the lower ground state manifold with hyperfine quantum number F=9/2. Photon scattering is introduced via a dedicated π-polarized plane wave laser beam at a detuning of 1.3 GHz below the D 2 line (see Appendix A). Starting from the F=9/2 manifold, the absorption and reemission of a photon can leave an atom back in its original state, but may also, with a probability of ≈ 33%, excite it to the upper F=7/2 ground state manifold. The detuning of the scattering beam is chosen such that the light is essentially resonant for atoms in the F=7/2 manifold, resulting in a quick transfer back to the lower manifold by resonant optical pumping, scattering typically one to five additional photons. Such 'scattering bursts', which start and end in the lower hyperfine manifold, happen on a timescale that is much shorter than the tunneling time τ = /J. Consequently, we can consider these bursts as single effective scattering events happening at a total rate of N · γ, where N denotes the atom number and γ the resulting single-particle scattering rate, which sets the effective coupling strength to the bath. We focus on the weak scattering regime ( γ J), where atoms can freely time evolve under the closed system Hamiltonian between successive scattering bursts. This is in stark contrast to the strong scattering limit, where a quantum Zeno effect would result in the localization of atoms [35,36,39]. Since we are close to the Paschen-Back regime, optical processes couple only weakly to the magnetic quantum number, such that a scattering burst will leave the spin state of the atom mostly unchanged. Details of the scattering bursts as well as the scattering beam are discussed in the Appendices A,B and C.
III. RESULTS Fig. 2 shows sample time traces of imbalance I and atom number N for various scattering rates γ at U = 2 J and moderate disorder ∆ = 4 J. As observed previously [13], the imbalance settles to a plateau at finite imbalance within a few tunneling times. In a perfectly isolated system, this finite imbalance would persist for all times, as is indicated by the numerical simulations. However, residual couplings between neighboring tubes [14], as well as off-resonant scattering of lattice photons limit the imbalance and atom number lifetimes to O(10 3 τ ) at the chosen parameters. Note that deeper in the localized phase we have observed significantly longer lifetimes [14].
For finite values of γ we observe a faster relaxation of the imbalance and an increased atom loss. This can be understood from a microscopic picture, in which scattering a photon results in the measurement of an atom's position on the length scale of the photon's wavelength λ ( Fig. 1(b)) [40]. Because of the relative size of the wavelength and the lattice constant d, in our experiment each scattering event can be interpreted as projecting the affected atom onto a single lattice site (λ/(2πd) ≈ 0.46) [41]. In this process, the probability for finding the atom on a specific final lattice site after the scattering event is given by the squared wavefunction overlap of the corresponding Wannier state with the atom's original state. This measurement effectively turns any coherent superposition of Wannier states into an incoherent mixture and can be described as dephasing the coherence terms in the initial single-particle density matrix at a rate γ dp = p dp ·γ, without altering the occupations [29,30,32]. Here p dp gives the probability of a scattering burst resulting in a dephasing event, where an atom remains in the lowest band (see Appendices B,C). Crucially, in the weak scattering limit considered in this work, time evolution under the closed system's Hamiltonian allows atoms to evolve into new coherent superpositions between successive scattering events. Since the new superpositions can be centered around a different lattice site than the original superposition, this effectively re-introduces hopping processes.
In addition, the induced measurement of the atom's position on a length scale λ/2 π implies a position measurement also within the lattice site, which can excite population to higher Bloch bands at a rate of γ ex = (1 − p dp ) · γ. These excitations ultimately result in atom loss, since weak trapping and strong tunnel couplings in higher excited bands allow the atoms to quickly tunnel out of the system. Note, however, that in our system atoms in the lowest longitudinally excited band remain trapped but are delocalized due to the higher tunneling rate (see Appendix D). Hence, in the presence of interactions, band excitations can contribute to the imbalance decay through both a complex rearrangement of the ground band wavefunction when an atom is excited, as well as through interactions of ground band atoms with delocalized atoms in higher bands.
The rates of the dephasing and band excitation processes sum up to the total scattering rate γ, which is controlled by the intensity of the scattering beam. In our system, the ratio γ dp /γ ex ≈ 2.3 is set by the lattice parameters as well as the wavelength of the scattering beam and is fixed throughout this work. We obtain its value from an ab-initio calculation of a scattering burst, which is discussed in detail in the Appendices B and C.
We quantify the imbalance relaxation via fits to a heuristic fit function, which for long times t decays by a stretched exponential of the form e −(Γ I t) β [29,30]. The inset shows measured imbalance decay rates ΓI as a function of the scattering rate γ at ∆ = 3 J. We observe a linear behavior, the slope of which is directly given by χ. We compare the data to the predictions of a rate model [30], indicated by the brown line, which is parallel to the fit through the experimental data. Hence experiment and theory give the same susceptibility. The offset is caused by the constant background decay Γ bg in the experiment.
The stretched exponential form arises naturally within a model with a spatial distribution of local relaxation rates [30]. We find that the global imbalance relaxation rate Γ I increases linearly with γ, i.e. (Γ I − Γ bg ) ∝ γ, in all parameter regimes (inset of Fig. 3 and Appendix H). This is consistent with an incoherent sum of two independent decay processes, namely the previously studied constant background decay Γ bg [14] and the effects of photon scattering (∝ γ). Motivated by the above observation we parameterize the imbalance relaxation rate as Γ I = χ · γ + Γ bg . For our system the susceptibility to photon scattering χ depends on the intrinsic system parameters (U, ∆), as well as the ratio of dephasing to excitation processes. Intuitively, 1/χ is a measure of the stability of a localized system to the effects of photon scattering and directly relates to the color gradients in Fig. 1a, with higher values of χ corresponding to a steeper gradient along γ. We first analyze the non-interacting case as it is exactly solvable and can be used to calibrate the experiment, before proceeding to the interacting case.

A. Non-interacting case
In the absence of interactions, we expect excitations of atoms to higher bands to have no effect on the imbal-ance, as they occur on even and odd sites with identical rates, and cannot affect the remaining atoms in the non-interacting case. Fig. 3 shows the non-interacting susceptibility χ as a function of disorder strength in the single-particle localized regime (∆ > 2 J). The susceptibility strongly decreases for increasing disorder strength, which can be understood by considering a single particle localized around a site i: Deep in the localized phase, its time averaged density distribution will be almost identical to that of the Wannier state on site i with almost no weight on neighboring sites. In this limit, a photon scattering event has negligible probability of moving the particle away from site i, resulting in a vanishing susceptibility. At weaker disorder strength, single-particle eigenstates are less localized and have finite overlap with the Wannier states of the neighboring sites. Hence, there is now a finite probability of scattering induced hopping transferring the particle to a neighboring site and thereby relaxing the imbalance, giving rise to a finite susceptibility. This intuitive idea is also at the heart of a recently proposed rate model [30], which we compare to our data (Fig. 3). Since the rate model describes only dephasing events, its scattering rate has been rescaled by p dp to take their finite probability into account. We find very good agreement between experiment and theory. This demonstrates that atom losses and excited band populations cannot affect the imbalance in the absence of interactions.
Our observable does not allow us to characterize the susceptibility at disorder strengths below ∆ 3 J, since close to the phase transition point the localization length becomes too large and the stationary imbalance of the closed system is already close to zero. However, we can derive a simple upper bound for the susceptibility based on the rate equation model: When the localization length diverges, each dephasing event has equal probability to project the atom onto an even or odd site, thereby canceling its contribution to the imbalance. In this limit, the imbalance thus decays with the rate γ dp , giving an upper bound to the susceptibility of χ ≤ p dp .

B. Interacting case
In the interacting case we expect a higher susceptibility, since now any dephasing event can also affect particles close by. Additionally, atom losses will now perturb the surrounding atoms and thereby further increase the susceptibility [30]. On top of these purely dissipative effects, also any delocalized atoms in excited bands can interact with the ground band. Hence, the non-interacting limit χ ≤ p dp no longer applies. In analogy to divergent susceptibilities at other phase transitions, one might in fact expect a divergent behavior at the MBL transition. As a consequence, even infinitesimally small couplings would dominate the dynamics close to the critical point, as is indicated in Fig. 1(a).
Measured susceptibilities at different disorder strengths versus interactions U . At finite interaction strengths, we compare our results to numerical TEBD simulations that do not include particle loss (triangles). The theoretical values at U = 0 are calculated from the rate model discussed earlier (squares). We observe a strong interaction dependence of the experimental susceptibilities close to the phase transition (∆ = 4 J), but only a weak effect deep in the localized phase (∆ = 6 J). Errorbars indicate the fit uncertainty. The solid lines are guides to the eye and the gray shaded region indicates the statistical uncertainty of the TEBD simulations. Fig. 4 shows the measured susceptibilities, as well as the results of TEBD simulations as a function of interaction strength. As in the non-interacting case, the numerical simulation does not implement atom loss, and hence its scattering rate has been rescaled accordingly. In the presence of interactions we do expect this to result in deviations between the experimental and numerical susceptibility, since excitations to higher bands will affect the imbalance. Deep in the localized phase, at ∆ = 6 J, we observe only a weak effect of interactions consistent with earlier works suggesting that interactions become less important at very strong disorder strengths [13]. In this regime we also find good agreement between theory and experiment, suggesting that atom losses only marginally affect the imbalance. However, at ∆ = 4 J we experimentally observe a strongly increasing susceptibility for growing interaction strengths, a trend that we expect to saturate at even larger U . This is in stark contrast to the TEBD result, which again approaches its non-interacting value at large U/J. While the TEBD simulations at U = 4 J suffer from large truncation errors (see Appendix J) and hence need to be considered with care, returning to the non-interacting value is the behavior expected for hardcore fermions, due to an exact mapping between the respective Hamiltonians [13]. However, this mapping breaks down when particle numbers are not conserved, suggesting that the difference between experiment and the TEBD simulations are most likely due to the effects of particles being excited to higher bands. Ex-perimentally disentangling the respective contributions of dephasing and particle excitations is not possible in our setup due to the fixed ratio of γ dp /γ ex .
An additional challenge is to unravel the effects of pure particle loss from the effects of the trapped but delocalized atoms accumulating in the first excited band. These atoms present an interesting field for future work, since they implement a 'small' bath, the properties of which might be strongly influenced by the back-action from the MBL system in the ground band [42][43][44].

IV. CONCLUSION
We have realized a controlled open MBL system by introducing dissipation via photon scattering and have found a stretched exponential decay of an initially imprinted charge-density wave in qualitative agreement with recent numerical studies [29][30][31]. Systematically varying the scattering rate γ enabled us to characterize the robustness of the MBL system via the definition of a susceptibility χ, which we found to be essentially independent of the interaction strength deep in the localized phase. Furthermore, we were able to experimentally study the interesting regime close to the MBL transition that is not accessible to current numerical studies, and have found an increasing susceptibility upon approaching the critical point. For the non-interacting system we derive an upper bound of χ ≤ p dp . However, we have found that interactions dramatically increase the system's susceptibility and speculate that they might even cause it to diverge at the MBL transition point, such that even infinitesimally small couplings would dominate the dynamics.
Our study paves the way towards a systematic characterization of the critical point by extrapolating the dynamics at finite coupling to the closed system limit. A complementary study in the ergodic Griffiths regime, where power-law decays of the imbalance are expected [19][20][21], would give insight into the delocalized side of the MBL transition. Furthermore, the applied scheme of implementing open quantum systems via controlled photon scattering is rather general and can straightforwardly be generalized to interesting delocalized states such as e.g. superfluids or topological insulators, where controlled dissipation appears to be essential to change the effective Chern number of a state [45,46]. ported in part by the National Science Foundation under Grant No. NSF PHY11-25915. MHF acknowledges additional support from the Swiss Society of Friends of the Weizmann Institute of Science.

Appendix A: Level scheme and scattering bursts
We prepare the 40 K atoms in our system in the lower lying F = 9/2 hyperfine manifold of their electronic ground state 4 2 S 1/2 , in an equal mixture of m F = −9/2 (|↓ ) and m F = −7/2 (|↑ ). We control the interaction strength U between our spins |↓ and |↑ using a Feshbach resonance centered around 202.1 G [47]. At these magnetic fields the level structure is close to the Paschen-Back regime, where m j and m i become good quantum numbers. This suppresses transitions between different m i states due to optical transitions to below 10% (See Appendix G). Hence, we can restrict the discussion of optical transitions to the quantum number m j .
A level scheme illustrating all levels and transitions important to the scattering of photons from our dedicated scattering beam is illustrated in Fig. 5. As is indicated by the red arrows, the scattering beam is π-polarized and its frequency is chosen such that atoms in the F = 9/2 ground state manifold see a detuning of roughly 1.3 GHz to the D 2 -line, while the upper ground state manifold (F = 7/2) is coupled resonantly to the excited states.
We define a scattering burst as a series of absorption and reemission processes of photons, where an atom both starts and ends in the lower lying F = 9/2 manifold. After absorption of an initial photon from the scattering beam, an atom can, via reemission, directly decay back into the lower lying F = 9/2 manifold, thereby ending the scattering burst, or, with a 33% probability, decay into the upper F = 7/2 manifold. Atoms that decayed to the F = 7/2 manifold experience resonant light and will be excited again. They will therefore quickly scatter multiple photons until, with 33% probability per scattering, they decay back into the F = 9/2 manifold, which ends the scattering burst. Hence, a scattering burst typically involves around 1-5 scattered photons.
Since the excitation from the upper F = 7/2 manifold is resonant, scattering rates from this state are much higher than tunneling rates in the lattice, and hence the total duration of a scattering burst is much shorter than a tunneling time. This means that atoms effectively remain frozen during a scattering burst. Therefore, considering the measurement, or dephasing, effect of a photon scattering, we can treat a full scattering burst as a single effective scattering event. Since the probabilities of exciting an atom into higher bands of the lattice increases with the number of scattered photons in a burst we use an average band excitation probability for the scattering bursts (see Appendix C).
The total rate of scattering bursts is controlled by the rate of absorbing the first photon from the lower lying F = 9/2 manifold. On this transition, the detuning and Here F labels the manifold that the state is adiabatically connected to at low magnetic fields. The quantum number mi is, to good approximation, not coupled to optical transitions as the system is close to the Paschen-Back regime. The two spin states |mI = −3 and |−4 adiabatically connect to the |mF = −7/2 and |−9/2 states at low magnetic fields. The dedicated scattering light is shown as red arrows, spontaneous emission processes are indicated as wavy lines, along with their branching ratios.
the low light intensities of below 3.6 µW/cm 2 result in scattering rates of only a few scattering bursts per atom per 100 tunneling times.

Appendix B: Single photon band excitation probabilities
Scattering photons gives rise to two processes: dephasing of atoms in the ground band and excitation of atoms to higher bands. Understanding the relative rates of ground band dephasing and excitations associated with the scattering bursts discussed in this work requires a detailed understanding of the processes associated with a single photon. We can calculate the single photon rates by calculating the band excitation probabilities of a stimulated absorption followed by a spontaneous emission. In this picture, dephasing will be associated with atoms remaining in the ground band. While in highly excited bands atoms will quickly be lost from the trap, the first excited band is an exception, as atoms in this band are trapped (see Appendix D).
The calculations are analogous to previous work on heating of atoms in dipole traps [40,41,[48][49][50], and are performed on the combination of three lattices along the three spatial directions for our system parameters. The primary lattice along the x-axis with λ p = 532. along the y-and z-axis have a depth of 36 E o R . Here E i R = h 2 /2mλ 2 i is the recoil energy corresponding to the wavelength of the lattice laser λ i and atomic mass m. For this calculation, we will neglect the weak (< 1E d R ) disorder lattice, assuming that it only marginally influences the bandstructure.
We calculate the excitation probabilities for atoms starting in a Wannier state of band (i x , i y , i z ). Stimulated absorption provides a momentum kick of k along the longitudinal direction, which is the direction of our scattering beam. Here k is the momentum of a photon from the dedicated scattering beam. Afterwards we act with another momentum kick of k, along an arbitrary direction to model spontaneous emission. The results of extracting the final excitation probabilities from band (i x , i y , i z ) into band (j x , j y , j z ), averaged over all emission directions, are shown in Table I.
Note that the excitation probabilities into the excited bands of the orthogonal lattices are equal due to symmetry and much lower than the excitation probabilities in the x-direction. This is due to the orthogonal lattices being deeper, and the momentum kicks from absorption of photons being along the x-direction due to the direction of travel of our dedicated scattering beam.
All estimations for the scattering bursts are based on the results of this calculation. Note that we are making small simplifications by using a Wannier state as the starting state, instead of the actual localized wavefunctions, which are a superposition of few Wannier states. Thereby we neglect the potential small effects from coherences. Also, we neglect return processes from 'higher' bands that are not explicitly considered, since atoms in these bands are not trapped. (see Appendix A) consists of a first photon absorbed from the F = 9/2 manifold, followed by a small number n ∈ [0, 1, 2, 3, ...] of photons scattered from the F = 7/2 transition, before returning to the F = 9/2 manifold. The band excitation probabilities of such a burst will depend heavily on the number of individual photons involved. We calculate the average band excitation probability of a scattering burst by averaging over all possible realizations of a burst, characterized by the number of cycled photons n. Specifically, we sum over the band excitation probabilities after scattering n photons (which steadily increases), weighted by the probability of scattering n photons in a burst, which is given by the geometric series P (n) = (1/3) · (2/3) n−1 (which quickly converges to zero). This sum is plotted in Fig. 6(a) as a function of the maximum number of photons considered in a burst.
We observe a quickly converging behavior of all populations after only a few photons. The limiting (n → ∞) values give the average excitation probabilities of a scattering burst.
Appendix D: Atom loss mechanism and finite excited band population In this work, we distinguish between two effects of a photon scattering burst. While dephasing is associated with events where atoms stay in the ground band, atom loss occurs due to particles being excited to higher bands and tunneling out of the system. Since atoms are mainly excited into higher bands of the longitudinal (x) lattice, we expect tunneling along this direction to constitute the main loss mechanism. Fig. 7 illustrates the energies of the ground, 1st and 2nd excited band of the longitudinal lattice as a function of real space position. The Gaussian shape trapping potential stems from the Gaussian beam shapes of the dipole trap. The bandwidth of the bands increases away from the trap center, because the orthogonal lattice beams also have a Gaussian shape. At a distance of 200 µm from the trap center, these beams have zero intensity and hence the atoms only experience the x-lattice. Vertical dashed lines mark the width of the atom cloud in the ground band, based on an in-situ measurement of the cloud. Indicated is the full 1/e 2 width of a Gaussian fit, which corresponds to approximately 200 lattice sites. Since the ground band is localized via disorder and photon assisted hopping gives only slow, diffusive spreading, we expect the cloud size to remain essentially constant during the dynamics.
In order to enable tunneling out of the system, a band needs to be i) delocalized (2 J band ≤ ∆ band ), where ∆ band is the disorder strength felt by atoms in the respective band, and ii) untrapped (4 J band > V trap , where V trap is the trap depth). Due to higher tunnel couplings, the first criterion is true for all longitudinally excited bands. A graphical visualization of the second condition is illustrated for the 1st and 2nd longitudinally excited band: A horizontal line from the upper band edge must not cross the lower band edge. This criterion is fulfilled for the 2nd longitudinally and higher excited bands and hence atoms in these bands can be lost from the system.
In the case of the 1st excited band, however, the line crosses the lower band edge, marking a finite size that atoms in the second band will expand to. We have checked these predictions by i) measuring the size that the 1st excited band expands to by deliberately loading atoms into the 1st excited band using the superlattice, letting them time evolve and imaging the cloud in-situ, as well as ii) directly measuring the lifetime of the 1st excited band. We obtained good agreement with the predicted size and found a lifetime similar to the lifetime of the ground band in the absence of photon scattering.
Furthermore, calculations of the bandstructure along the orthogonal (y, z) direction shows that, due to the deeper lattices, both the 1st and 2nd excited band are trapped. While photon scattering only excites a few atoms into these bands, they might nonetheless become relevant, since at the spatial edges of the system the first excited bands along y and z are resonant with the first excited band along x, enabling transfer of atoms between the bands. We have experimentally checked this by preparing atoms in the first excited band and found that atoms indeed distribute between the 1st excited bands in all three directions. The 1st excited bands being trapped will result in a finite population in these bands building up. Fig. 6(b) shows the band populations versus the number of scattering bursts. The values are calculated using rate equations based on the average excitation probabilities of a scattering burst. The total population in the 1st excited band quickly builds up to about 15% of the initial atoms before slowly decaying. Due to the decay of the ground band population, the 1st excited band population quickly reaches a significant portion of the ground band population, namely approximately 30% after 4 scattering bursts.

Appendix E: Effects of the finite excited band populations
Trapped atoms in the excited bands can affect the imbalance in the ground band in multiple ways. The most direct way of influence is due to the imaging procedure not being able to distinguish fully between higher and ground band atoms, which directly affects the measured imbalance. However, we believe this effect to be small, since the atoms distribute among bands in all directions such that the individual populations are too small and vanish in the noise. Furthermore, we find good agreement with theory in the non-interacting case.
In the presence of interactions one can envisage another possible channel of influence, as atoms in the higher bands, which are delocalized, can act as a bath for atoms in the ground band. While the coupling to this bath should be rather weak due to the bigger spatial size of the 1st longitudinally excited band, it might be significant in certain regimes. While models describing such a two band behavior have been studied theoretically [42][43][44], those studies have been limited to very small system sizes. A detailed study of the effects of higher band population on the ground band would constitute a particularly interesting future direction for this work.

Appendix F: Calibration of the scattering rates
In the experiment, we vary the amount of scattering light by controlling the intensity of the scattering beam via an acousto-optic modulator and stabilize the total power using a calibrated photodiode. We calibrate the photodiode via the intensity profile of the scattering beam by imaging it at the position of the atoms and comparing it to an in-situ image of the atomic cloud. From these images we can obtain the average intensity I at the position of the atoms. Finally, the scattering rate can be calculated as Here ω D2 and Γ D2 denote the transition frequency and the decay rate of the D 2 line, respectively. The detuning δ sc refers to the detuning seen by atoms in the lower F = 9/2 hyperfine manifold of the ground state, since the absorption from this state controls the rate of scattering bursts. Due to the detuning being δ sc ≈ 1.3 GHz, we can neglect the effects of the D 1 line, which is much further away, and assume that we do not resolve the hyperfine levels of the excited state, allowing us to use this simple formula.

Estimating the relative dephasing rate
Comparing the experimental data to theory, which only includes the effects of dephasing, requires an estimation of the fraction of scattering bursts resulting only in dephasing γ dp /γ = p dp . Ignoring back-transfer processes from the 1st excited to the ground band (which would change the ground band population by ∼ 1%), this is equal to the probability of staying in the ground band during an average scattering burst, which was calculated earlier. This gives a relative dephasing rate of γ dp /γ ≈ 70%.

Estimating the relative loss rate
In order to check our calibration of atom loss and the calculations on band excitations we estimate the expected loss rate and compare it to the experimentally measured atom number decay. Since the 1st excited band along x and the higher bands along the orthogonal directions are trapped, the loss rate should be equal to the rate at which atoms are excited to the 2nd excited band along x.
By summing the exact probabilities of an atom being excited to the 2nd excited band after n scattering bursts ( Fig. 6(b)), weighted by the probability of scattering n photons in time t, which is given by a Poisson distribution we can calculate the probability of staying in the system (the probability of staying in the ground or 1st excited bands) until time t. While this is not strictly an exponential decay, it can be approximated as such, allowing us to extract an effective loss rate γ el ≈ 0.175γ ( Fig. 8(a)). Using this rate, the linear relationship between atom loss rate and scattering rate plotted in Fig. 9(b), allows us to define an atom number susceptibility (F3) Fig. 8(b) shows the atom number susceptibility for the non-interacting case for various ∆. We observe a noisy behavior consistent with no trend along ∆. Perfect agreement with our model would be indicated by χ N = 1.
We observe values slightly below one, indicating that our model describes excitation processes reasonably well.
One possible explanation for the observed minor differences is the experimental extraction of the atom number lifetime. Since the time traces were only taken up to times where the imbalance reaches zero, the atom numbers have often not fully decayed yet, rendering the exponential fit unreliable.

Appendix G: Spin flip probabilities
At the magnetic fields of around 200 G used in our experiment, the 2 P 3/2 excited state manifold is deep in the Paschen-Back regime. However, the 2 S 1/2 ground state manifold still has a weak coupling between the nuclear and electronic spins, causing a finite probability of changing m i by scattering a photon.
We can calculate the probability of a spin flip by including nuclear spin in our calculation of scattering rates and branching ratios. We find probabilities of 4% for |m F = −9/2 and 10% for |m F = −7/2 . Note that most of the spin flips will simply convert atoms between |m F = −9/2 and |m F = −7/2 . However also |m F = −5/2 states can be created, which would have a different interaction strength. But given the minimal excitation probabilities, we expect any effects due to these additional spin states to be negligible.
Appendix H: Decay rate scaling with scattering rate As discussed in the main paper, the imbalance decay rate shows a linear behavior with the scattering rate. Fig. 9(a) shows further exemplary data for various interaction strengths at ∆ = 4 J. The susceptibilities are extracted via a linear fit to this data. The errorbars plotted in Figs. 3,4 are calculated as the maximum of i) the square root of the covariance error of the fit and ii) the results of linear fits through the imbalance decay rates plus or minus their respective errorbars. Fig. 9(b) shows experimental data for the atom number decay rate as a function of the calculated band excitation rate γ ex = (1 − p dp ) · γ. We observe the expected linear trend, but with a slope of ≈ 0.5, smaller than one. This indicates that not all atoms excited to higher bands are lost. Indeed we find that atoms in the first excited band of the longitudinal lattice remain trapped (see Appendix D).

Appendix I: Fitfunction for the Imbalance
The decay of the imbalance is fitted by an initial oscillation that decays exponentially (with amplitude A, frequency ω, lifetime τ ) to an offset value o, which corresponds to the stationary imbalance of the closed system. This function is multiplied by a stretched exponential decay to zero [14], that models the long term decay studied here. . Imbalance and atom number decay rate as a function of the scattering rate: (a) Imbalance decay rate for various interaction strengths at ∆ = 4 J. Our data is consistent with a linear behavior for all parameter values. The imbalance decay rates additionally show a background decay rate at γ dp = 0, which heavily depends on both the interaction and the disorder strength. (b) Atom number decay rates as a function of the band excitation rate. Again we find a linear scaling.
Here Γ I gives the imbalance decay rate and β the stretching exponent. We find that a stretched exponential fit describes our data much better than a simple exponential decay. We do not perform a systematical analysis of the stretching exponents β, since their values depend heavily on the behavior after long times, where the systems response is heavily affected by the effects of trapped atom's in higher bands, a strongly reduced density and the creation of additional spin states. Fit values for β scatter between typical values of 0.5 ≤ β ≤ 1 and show large errorbars.
whereĉ † i,s (ĉ i,s ) creates (annihilates) a particle at site i with spin s andn i,s =ĉ † i,sĉi,s . The disorder potential for the Aubry-André model is given by V i = ∆ cos(2παi + φ) with α the ratio of the lattice periodicities and φ is a random phase.
To simulate the time evolution of the open system, we introduce the density matrix, for which the time evolution is given by the Lindblad equatioṅ Here, the first term describes the unitary time evolution and the second term the coupling of the system to the environment. The jump operators L µ denote the system operators directly coupled to the bath, which in our case are given by local density measurements, i.e., L µ =n i,s . We then simulate the time evolution of the quantum Lindblad equation (J2) for system size S = 20 and 30 disorder realizations using the time-evolving block decimation (TEBD) scheme for matrix product operators [51].
We note that the results of the numerical calculations, in particular close to the transition and for intermediate interactions, should be treated with some care. Since the local Hilbert-space dimension of the density matrix is d = 16, we could not increase the matrix product state bond dimension to more than 100. For such bond dimension, the truncation error grows rapidly to ≈ 0.1 per bond before decreasing again. The actual error, however, can not be deduced from this truncation error. The error bars in Fig. 4 are the statistical errors from averaging over the random phases φ.