We experimentally observe many-body localization of interacting fermions in a one-dimensional quasi-random optical lattice. We identify the many-body localization transition through the relaxation dynamics of an initially-prepared charge density wave. For sufficiently weak disorder the time evolution appears ergodic and thermalizing, erasing all remnants of the initial order. In contrast, above a critical disorder strength a significant portion of the initial ordering persists, thereby serving as an effective order parameter for localization. The stationary density wave order and the critical disorder value show a distinctive dependence on the interaction strength, in agreement with numerical simulations. We connect this dependence to the ubiquitous logarithmic growth of entanglement entropy characterizing the generic many-body localized phase.

Observation of many-body localization of interacting fermions in a quasi-random optical lattice

Michael Schreiber, Sean S. Hodgman, Pranjal Bordia, Henrik P. Lüschen, Mark H. Fischer, Ronen Vosk, Ehud Altman, Ulrich Schneider and Immanuel Bloch

Fakultät für Physik, Ludwig-Maximilians-Universität München, Schellingstr. 4, 80799 Munich, Germany

Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Str. 1, 85748 Garching, Germany

Department of Condensed Matter Physics, Weizmann Institute of Science, Rehovot 76100, Israel

## Introduction

The ergodic hypothesis is one of the central principles of statistical physics. In ergodic time evolution of a quantum many-body system, local degrees of freedom become fully entangled with the rest of the system, leading to an effectively classical hydrodynamic evolution of the remaining slow observables [1]. Hence, ergodicity is responsible for the demise of observable quantum correlations in the dynamics of large many-body systems and forms the basis for the emergence of local thermodynamic equilibrium in isolated quantum systems [2, 3, 4]. It is therefore of fundamental interest to investigate how ergodicity breaks down and search for alternative, genuinely quantum paradigms in the dynamics, and to understand the long-time stationary states that ensue in the absence of ergodicity.

One path to breaking ergodicity is provided by the study of integrable models, where thermalization is prevented due to the constraints imposed on the dynamics by an infinite set of conservation rules. Such models have been realized and studied in a number of experiments with ultracold atomic gases [5, 6, 7]. However, integrable models represent very special and fine-tuned situations, making it difficult to extract general underlying principles.

Theoretical studies over the last decade point to many-body localization (MBL) in a disordered isolated quantum system as a more generic alternative to thermalization dynamics. In his original paper on single-particle localization, Anderson already speculated that interacting many-body systems subject to sufficiently strong disorder would also fail to thermalize [8]. Only recently, however, have convincing theoretical arguments been put forward that Anderson localization remains stable under the addition of moderate interactions, even in highly excited many-body states [9, 10, 11]. Further theoretical studies have established the many-body localized state as a distinct dynamical phase of matter that exhibits novel universal behavior [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. In particular, the relaxation of local observables does not follow the conventional paradigm of thermalization and is expected to show explicit breaking of ergodicity. In many ways, the MBL transition is fundamentally different from all other known transitions [23, 24]. On one side of the transition ergodicity prevails and quantum effects decay at long times, whereas on the other side quantum correlations persist indefinitely. Hence the MBL transition sets a sharp boundary between a macroscopic world showing quantum phenomena and one governed by classical physics.

While Anderson localization of non-interacting particles has been experimentally observed in a variety of systems, including light scattering from semiconductor powders in 3D [25], photonic lattices in 1D [26] and 2D [27] and cold atoms in 1D and 3D random [28, 29, 30] and quasi-random [31] disorder, the interacting case has proven more elusive. Initial experiments with interacting systems have focused on the superfluid [32, 33] or metal [34] to insulator transition in the ground state. Evidence for inhibited macroscopic mass transport was reported even at elevated temperatures [34], but is hard to distinguish from exponentially slow motion expected from conventional activated transport or effects stemming from the inhomogeneity of the cloud. Until now conclusive experimental evidence for many-body localization at finite energy density has thus been lacking.

In this paper we report the first experimental observation of ergodicity breaking due to many-body localization. Our experiments are performed in a one-dimensional system of ultracold fermions in a bi-chromatic, quasi-randomly disordered lattice potential. We identify the many-body localized phase by monitoring the time evolution of local observables following a quench of system parameters. Specifically, we prepare a high-energy initial state with strong charge density wave (CDW) order (as shown in Fig. 1A) and measure the relaxation of this charge density wave in the ensuing unitary evolution. Our main observable is the imbalance between the respective atom numbers on even () and odd () sites

(0) |

which directly measures the CDW order. While the initial CDW () will quickly relax to zero in the thermalizing case, this is not true in a localized system, where ergodicity is broken and the system cannot act as its own heat bath (Fig. 1B) [35]. Intuitively, if the system is strongly localized, all particles will stay close to their original positions during time evolution, thus only smearing out the CDW a little. A longer localization length corresponds to more extended states and will lead to a lower steady state value of the CDW. The long-time stationary value thus effectively serves as an order parameter of the MBL phase and allows us to map the phase boundary between the ergodic and non-ergodic phases in the parameter space of interaction versus disorder strength. In particular, in the non-interacting system the CDW vanishes asymptotically as [36]. In contrast to previous experiments, which studied the effect of disorder on the global expansion dynamics [28, 31, 32, 34, 33], the CDW order parameter acts as a purely local probe, directly capturing the ergodicity breaking.

Our system can be described by the one-dimensional fermionic Aubry-André model [37] with interactions [35], given by the Hamiltonian

(0) |

Here, is the tunneling matrix element between neighboring lattice sites and () denotes the creation (annihilation) operator for a fermion in spin state on site . The second term describes the quasi-random disorder, i.e. the shift of the on-site energy due to an additional incommensurate lattice, characterized by the ratio of lattice periodicities , disorder strength and phase offset . Lastly, represents the on-site interaction energy and is the local number operator (see Fig. 1C).

This quasi-random model is special in that, for almost all irrational [36], all single particle states become localized at the same critical disorder strength [37]. For larger disorder strengths the localization length decreases monotonically. Such a transition was indeed observed experimentally in a non-interacting bosonic gas [31]. In contrast, truly random disorder will lead to single-particle localization in one dimension already for arbitrarily small disorder strengths. Previous numerical work indicates many-body localization in quasi-random systems to be similar to that obtained for a truly random potential [35].

## Experiment

We experimentally realize the Aubry-André model by superimposing on the primary, short lattice () a second, incommensurate disorder lattice with (thus ) and control , and via lattice depths and relative phase between the two lattices [36]. The interactions () between atoms in the two different spin states and are tuned via a magnetic Feshbach resonance [36]. In total, this provides independent control of , and and enables us to continuously tune the system from an Anderson insulator in the non-interacting case to the MBL regime for interacting particles.

An additional long lattice () forms a period-two superlattice [38, 39] together with the short lattice and is employed during the preparation of the initial CDW state, and during detection [36]. Deep lattices along the orthogonal directions ( and ), create an array of decoupled 1D tubes. Here, denotes the recoil energy, with being Planck’s constant, the mass of the atoms and the respective wavelength of the lattice lasers.

We employ a two component degenerate Fermi gas of K atoms, consisting of an equal mixture of 25-30 atoms in each of the two lowest hyperfine states and , at an initial temperature of 0.24(2) , where is the Fermi temperature. The atoms are initially prepared in a finite temperature band insulating state [40] in the long and orthogonal lattices. We then split each lattice site by ramping up the short lattice in a tilted configuration [36] and subsequently ramp down the long lattice. This creates a charge density wave, where there are no atoms on odd lattice sites but zero, one or two atoms on each even site [39, 41]. This initial CDW is then allowed to evolve for a given time in the 8.0(2) deep short lattice at a specific interaction strength in the presence of disorder . In a final step, we detect the number of atoms on even and odd lattice sites by employing a band-mapping technique which maps them to different bands of the superlattice [41, 36]. This allows us to directly measure the imbalance , as defined in Eq. (Introduction).

## Results

We track the time evolution of the imbalance for various interactions and disorder strengths (see Fig. 2). At short times the imbalance exhibits some dynamics consisting of a fast decay followed by a few damped oscillations. After a few tunneling times the imbalance approaches a stationary value. In a clean system () and for weak disorder, the stationary value of the imbalance approaches zero. For stronger disorder, however, this behaviour changes dramatically and the imbalance attains a non-vanishing stationary value that persists for all observation times. Since the imbalance must decay to zero on approaching thermal equilibrium at these high energies, the non-vanishing stationary value of directly indicates non-ergodic dynamics. Deep in the localized phase, where unbiased numerical Density-Matrix Renormalisation Group (DMRG) calculations are feasible due to the slow entanglement growth, we find the stationary value obtained in the simulations to be in very good agreement with the experimental result. These simulations were performed for a single homogeneous tube without any trapping potentials [36]. The stronger damping of oscillations observed in the experiment can be attributed to a dephasing caused by variations in between different 1D tubes [41, 36].

We experimentally observe an additional very slow decay of on a timescale of several hundred tunneling times for all interaction strengths, which we attribute to the fact that our system is not perfectly closed due to small losses, technical heating and photon scattering [42, 36]. Another potential mechanism for delocalization at long times is related to the intrinsic SU(2) spin symmetry in our system [43]. However, for the relevant observation times our numerical simulations do not indicate the presence of such a thermalization process.

To characterize the dependence of the localization transition on and , we focus on the stationary value of , plotted in Fig. 3 for non-interacting atoms and in Fig. 4 for interacting atoms. For non-interacting atoms (Fig. 3), the measured imbalance is compatible with extended states within the finite, trapped system for . Above the critical point of the homogeneous Aubry-André model at [37], however, the measured imbalance strongly increases as the single-particle eigenstates become more and more localized. The observed transition agrees well with our theoretical modeling including the harmonic trap [36].

The addition of moderate interactions slightly reduces the degree of localization compared to the non-interacting case, i.e. they decrease the imbalance and hence increase the critical value of necessary to cross the delocalization-localization transition (Fig. 4A and B). Importantly, we find that localization persists for all interaction strengths. For a given disorder, the imbalance decreases up to a value of before increasing again. For large , the system even becomes more localized than in the non-interacting case. This can be understood qualitatively by considering an initial state consisting purely of empty sites and sites with two atoms (doublons): for sufficiently strong interactions, isolated doublons represent stable quasiparticles as the two atoms cannot separate and hence only tunnel with an effective second-order tunneling rate of [44, 45]. This strongly increases the effective disorder and promotes localization. In the experiment, the initial doublon fraction is well below one [36] and the density is finite, such that we observe a weaker effect. We find the localization dynamics and the resulting stationary values to be symmetric around , highlighting the dynamical symmetry of the Hubbard Hamiltonian for initially localized atoms [46]. The effect of interactions can be seen in the contour lines (Fig. 4A, dotted white lines) as well as directly in the characteristic ‘W’ shape of the imbalance at constant disorder (Fig. 4B), demonstrating the re-entrant behaviour of the system [22]. This behaviour extends to our best estimate of the localization transition, which is shown in Fig. 4A as the solid white line.

We can gain additional insight into how localization changes with interaction strength by computing the growth of the entanglement entropy between the two halves of the system during the dynamics, as shown in Fig. 5A. For long times, we observe a logarithmic growth of the entanglement entropy with time as , which is characteristic of the MBL phase [12, 13]. The slope is proportional to the bare localization length , which in a weakly interacting system in the localized phase corresponds to the single particle localization length. In general, is the characteristic length over which the effective interactions between the conserved local densities decay [17, 18] and connects to the many-body localization length deep in the localized phase. In contrast to , however, is expected to remain finite at the transition [23]. We find to exhibit a broad maximum for intermediate interaction strengths (Fig. 5B), corresponding to a maximum in the thus inferred localization length. This maximum in turn leads to a minimum in the CDW value. The characteristic ‘W’ shape in the imbalance is thus directly connected to the maximum in the entanglement entropy slope, as both are consequences of the maximum in localization length. Equivalent information on the localization properties as obtained from the entanglement entropy can be gained in experiments by monitoring the temporal decay of fluctuations around the stationary value of the CDW [36]. While we do not have sufficient sensitivity to measure these fluctuations in the current experiment, we expect them to be accessible to experiments with single site resolution [47, 48].

To systematically study the effect of the initial energy density on the MBL phase, we load the lattice using either attractive, vanishing or repulsive interactions (Fig. 6), thereby predominantly changing the number of doublons in the initial state [36]. Since the initial state consists of fully localized particles only, the local energy density is directly given by the product of interaction strength and doublon density. We find that for an interaction strength during the evolution of the energy density does not significantly affect the localization properties, proving that MBL persists over a wide energy range. For , localization properties depend significantly on the doublon fraction because of the second emerging energy scale , as discussed above. For the case of repulsive loading, i.e. a low fraction of doubly occupied sites, the imbalance for and strong interactions match within error. Indeed, a rigorous mapping can be made between the non-interacting system and the dynamics in the doublon free subspace at strong interactions [36]. At very large interactions and high doublon fractions, the additional long timescales start to also compete with heating and loss processes, rendering the definition of stationary states challenging.

## Conclusion

We have created and characterized the many-body localized phase of a system of interacting ultracold fermions in a quasi-random optical lattice. The dependence of the MBL phase-transition point on interactions was measured by studying the evolution of an artificially created charge density wave. Moderate interactions have a delocalizing effect and increase the critical disorder strength. We have also shown the MBL phase to be stable over a wide range of energy densities.

Our experimental demonstration of ergodicity breaking due to many-body localization paves the way for many further investigations. An interesting extension would be to use ‘true’ random disorder created by e.g. an optical speckle pattern, as has been used to study Anderson localization [28]. Another important next step is extending the present study to higher dimensions. Additional insight can also be gained by analyzing the full relaxation dynamics of local observables [19, 20, 21], in an experimental setup featuring single site resolution [47, 48]. For instance, the decay of fluctuations of with time could be directly measured, providing an even more direct connection to the entanglement entropy. Another important direction for future investigation is the effect of opening the system in a controlled way. This could be done, e.g. by adding a near-resonant laser to deliberately enhance photon scattering or by employing a Bose-Fermi mixture, in which excitations of the BEC form a well controlled bath for the fermions. This will allow a systematic study of the critical dynamics associated with the MBL phase transition, where the bath relaxation time now provides the only scale. Such a study would also allow the MBL phase to be clearly distinguished from classical glassy dynamics. The latter, unlike MBL, is insensitive to coupling of the system to an external bath.

## References:

- [1] J. Lux, J. Müller, A. Mitra, A. Rosch, Hydrodynamic long-time tails after a quantum quench. Phys. Rev. A 89, 053608 (2014).
- [2] J. M. Deutsch, Quantum statistical mechanics in a closed system. Phys. Rev. A 43, 2046 (1991).
- [3] M. Srednicki, Chaos and quantum thermalization. Phys. Rev. E 50, 888 (1994).
- [4] M. Rigol, V. Dunjko, M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems. Nature 452, 854 (2008).
- [5] B. Paredes, A. Widera, V. Murg, O. Mandel, S. Fölling, I. Cirac, G. V. Shlyapnikov, T. W. Hänsch, I. Bloch, Tonks-Girardeau gas of ultracold atoms in an optical lattice. Nature 429, 277 (2004).
- [6] T. Kinoshita, T. Wenger, D. S. Weiss, A quantum Newton’s cradle. Nature 440, 900 (2006).
- [7] M. Gring, M. Kuhnert, T. Langen, T. Kitagawa, B. Rauer, M. Schreitl, I. Mazets, D. Adu Smith, E. Demler, J. Schmiedmayer, Relaxation and prethermalization in an isolated quantum system. Science 337, 1318 (2012).
- [8] P. W. Anderson, Absence of Diffusion in Certain Random Lattices. Phys. Rev. 109, 1492 (1958).
- [9] D. M. Basko, I. L. Aleiner, B. L. Altschuler, Metal-insulator transition in a weakly interacting many-electron system with localized single-particle states. Ann. Phys. 321, 1126 (2006).
- [10] I. V. Gornyi, A. D. Mirlin, D. G. Polyakov, Interacting electrons in disordered wires: Anderson localization and low-T transport. Phys. Rev. Lett. 95, 206603 (2005).
- [11] J. Z. Imbrie, On many-body localization for quantum spin chains. arXiv:1403.7837 (2014).
- [12] M. Žnidarič, T. Prosen, P. Prelovšek, Many-body localization in the Heisenberg XXZ magnet in a random field. Phys. Rev. B 77, 064426 (2008).
- [13] J. H. Bardarson, F. Pollmann, J. E. Moore, Unbounded growth of entanglement in models of many-body localization. Phys. Rev. Lett. 109, 017202 (2012).
- [14] B. Bauer, C. Nayak, Area laws in a many-body localized state and its implications for topological order. J. Stat. Mech. Theor. Exp. 09, P09005 (2013).
- [15] R. Vosk, E. Altman, Many-body localization in one dimension as a dynamical renormalization group fixed point. Phys. Rev. Lett. 110, 067204 (2013).
- [16] M. Serbyn, Z. Papić, D. A. Abanin, Universal slow growth of entanglement in interacting strongly disordered systems. Phys. Rev. Lett. 110, 260601 (2013).
- [17] M. Serbyn, Z. Papić, D. A. Abanin, Local conservation laws and the structure of the many-body localized states. Phys. Rev. Lett. 111, 127201 (2013).
- [18] D. A. Huse, R. Nandkishore, V. Oganesyan, Phenomenology of fully many-body-localized systems. Phys. Rev. B. 90, 174202 (2014).
- [19] F. Andraschko, T. Enss, J. Sirker, Purification and many-body localization in cold atomic gases. Phys. Rev. Lett. 113, 217201 (2014).
- [20] R. Vasseur, S. A. Parameswaran, J. E. Moore, Quantum revivals and many-body localization. arXiv:1407.4476 (2014).
- [21] M. Serbyn, Z. Papić, D. A. Abanin, Quantum quenches in the many-body localized phase. Phys. Rev. B. 90, 174302 (2014).
- [22] Y. Bar Lev, G. Cohen, D. R. Reichman, Absence of diffusion in an interacting system of spinless fermions on a one-dimensional disordered lattice. arXiv:1407.7535 (2014).
- [23] R. Vosk, D. A. Huse and E. Altman, Theory of the many-body localization transition in one dimensional systems. arXiv:1412.3117 (2014).
- [24] A. C. Potter, R. Vasseur, S. A. Parameswaran, Theory of phase transitions from quantum glasses to thermal fluids. arXiv:1501.03501 (2015).
- [25] D. S. Wiersma, P. Bartolini, A. Lagendijk, R. Righini, Localization of light in a disordered medium. Nature 390, 671 (1997).
- [26] Y. Lahini, A. Avidan, F. Pozzi, M. Sorel, R. Morandotti, D. N. Christodoulides, Y. Silberberg, Anderson localization and nonlinearity in one-dimensional disordered photonic lattices. Phys. Rev. Lett. 100, 013906 (2008).
- [27] T. Schwartz, G. Bartal, S. Fishman, M. Segev, Transport and Anderson localization in disordered two-dimensional photonic lattices. Nature 446, 52 (2007).
- [28] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, D. Clément, L. Sanchez-Palencia, P. Bouyer, A. Aspect, Direct observation of Anderson localization of matter waves in a controlled disorder. Nature 453, 891 (2008).
- [29] S. S. Kondov, W. R. McGehee, J. J. Zirbel, B. DeMarco, Three-Dimensional Anderson Localization of Ultracold Matter. Science 334, 66 (2011).
- [30] F. Jendrzejewski, A. Bernard, K. Müller, P. Cheinet, V. Josse, M. Piraud, L. Pezzé, L. Sanchez-Palencia, A. Aspect, P. Bouyer, Three-dimensional localization of ultracold atoms in an optical disordered potential. Nat. Phys. 8, 398 (2012).
- [31] G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort, M. Zaccanti, G. Modugno, M. Modugno, M. Inguscio, Anderson localization of a non-interacting Bose-Einstein condensate. Nature 453, 895 (2008).
- [32] B. Deissler, M. Zaccanti, G. Roati, C. D’Errico, M. Fattori, M. Modugno, G. Modugno, M. Inguscio, Delocalization of a disordered bosonic system by repulsive interactions. Nat. Phys. 6, 354 (2010).
- [33] C. D’Errico, E. Lucioni, L. Tanzi, L. Gori, G. Roux, I. P. McCulloch, T. Giamarchi, M. Inguscio, G. Modugno, Observation of a Disordered Bosonic Insulator from Weak to Strong Interactions, Phys. Rev. Lett. 113, 095301 (2014).
- [34] S. S. Kondov, W. R. McGehee, W. Xu, B. DeMarco, Disorder-induced localization in a strongly correlated atomic Hubbard gas. arXiv:1305.6072 (2013).
- [35] S. Iyer, V. Oganesyan, G. Refael, D. A. Huse, Many-body localization in a quasiperiodic system. Phys. Rev. B 87, 134202 (2013).
- [36] See Supplementary Material for details.
- [37] S. Aubry, G. André, Analyticity breaking and Anderson localization in incommensurate lattices. Ann. Israel Phys. Soc. 3, 133 (1980).
- [38] J. Sebby-Strabley, M. Anderlini, P. S. Jessen, J. V. Porto, Lattice of double wells for manipulating pairs of cold atoms. Phys. Rev. A 73, 033605 (2006).
- [39] S. Fölling, S. Trotzky, P. Cheinet, M. Feld, R. Saers, A. Widera, T. Müller, I. Bloch, Direct observation of second-order atom tunnelling. Nature 448, 1029 (2007).
- [40] U. Schneider, L. Hackermüller, S. Will, Th. Best, I. Bloch, T. A. Costi, R. W. Helmes, D. Rasch, A. Rosch, Metallic and insulating phases of repulsively interacting fermions in a 3D optical lattice. Science 322, 1520 (2008).
- [41] S. Trotzky, Y-A. Chen, A. Flesch, I. P. McCulloch, U. Schollwöck, J. Eisert, I. Bloch, Probing the relaxation towards equilibrium in an isolated strongly correlated one-dimensional Bose gas. Nat. Phys. 8, 325 (2012).
- [42] H. Pichler, A. J. Daley, P. Zoller, Nonequilibrium dynamics of bosonic atoms in optical lattices: decoherence of many-body states due to spontaneous emission. Phys. Rev. A 82, 063605 (2010).
- [43] R. Vasseur, A. C. Potter, S.A. Parameswaran, Dynamics of hot random quantum spin chains: from anyons to Heisenberg spins. arXiv:1410.6165 (2014).
- [44] K. Winkler, G. Thalhammer, F. Lang, R. Grimm, J. Hecker Denschlag, A. J. Daley, A. Kantian, H. P. Büchler, P. Zoller, Repulsively bound atom pairs in an optical lattice. Nature 441, 853 (2006).
- [45] S. Trotzky, P. Cheinet, S. Fölling, M. Feld, U. Schnorrberger, A. M. Rey, A. Polkovnikov, E. A. Demler, M. D. Lukin, I. Bloch, Time-resolved observation and control of superexchange interactions with ultracold atoms in optical lattices. Science 319, 295 (2008).
- [46] U. Schneider, L. Hackermüller, J. P. Ronzheimer, S. Will, S. Braun, T. Best, I. Bloch, E. Demler, S. Mandt, D. Rasch, A. Rosch, Fermionic transport and out-of-equilibrium dynamics in a homogeneous Hubbard model with ultracold atoms. Nat. Phys. 8, 213 (2012).
- [47] W. S. Bakr, J. I. Gillen, A. Peng, S. Fölling, M. Greiner, A quantum gas microscope for detecting single atoms in a Hubbard-regime optical lattice. Nature 462, 74 (2009).
- [48] J. F. Sherson, C. Weitenberg, M. Endres, M. Cheneau, I. Bloch, S. Kuhr, Single-atom-resolved fluorescence imaging of an atomic Mott insulator. Nature 467, 68 (2010).

Acknowledgments: We acknowledge technical assistance by D. Garbe and F. Görg during the setup of the experiment. We acknowledge financial support by the Deutsche Forschungsgemeinschaft (FOR801, Deutsch-Israelisches Kooperationsprojekt Quantum phases of ultracold atoms in optical lattices), the European Commision (UQUAM, AQuS), the US Defense Advanced Research Projects Agency (Quantum Emulations of New Materials Using Ultracold Atoms), the Minerva Foundation, ISF grant # 1594/11 and Nanosystems Initiative Munich (NIM).

Supplementary Material:

Supplementary Text

## General sequence

The experimental sequence for producing degenerate Fermi gases of K via sympathetic cooling with Rb in a magnetic quadrupole and optical dipole trap followed by evaporative cooling has been described elsewhere [40]. Final clouds contain typically 50-60k K atoms at a temperature of 0.24(2) times the Fermi temperature , with 50(3)% of atoms in each of the and spin states. We tune the interactions between the two spin states via a Feshbach resonance at 202.1G [50]. All experiments described in this work employ absorption imaging after time of flight (TOF) along the axis ( axis is oriented along the vertical direction). This is orthogonal to the axis along which the superlattice is aligned and all experiments take place. For the experiments shown in Fig. 3, where we employ a single spin state, we set the magnetic field early in the evaporation sequence to , which is the centre of a Feshbach resonance between K atoms and Rb atoms in the state [51]. This results in all atoms being lost (along with some Rb atoms). By adjusting the cooling sequence slightly we are able to prepare a sample with 98(5)% of atoms in the state with the same temperature and atom number as in the usual sequence.

## Lattice Details:

### Lattice setup

All lattice potentials result from standing waves created by retro-reflected single-frequency laser beams along the , and axes. All beams have a Gaussian profile and are focussed down to waists of at the positions of the atoms. As our cloud has a FWHM of in the horizontal plane and in the vertical direction, this results in slight variations of the lattice parameters , and across the cloud (see below).

### Superlattice lock

In our superlattice setup the laser (long lattice) serves as a master oscillator and its frequency is locked to a thermally stabilized and acoustically isolated Fabry-Pérot cavity. The absolute short term stability of the master oscillator is over and is characterized through the residual locking error. The short lattice laser at is offset locked relative to the second harmonic of the long lattice laser. The offset lock frequency sets the relative phase of the superlattice potential [39]. We obtain a relative line width of the offset lock of , which corresponds to a phase fluctuation of . This also dominates the absolute short term stability of the short lattice laser.

### Phase adjustment of the superlattice

In order to calibrate the phase of the superlattice, we prepare a spinpolarized band insulating cloud in a combination of long and perpendicular lattices and subsequently split each well non-adiabatically by adding the short lattice. Time of flight pictures of such arrays of double wells show the typical double-slit interference pattern only for phase settings of the superlattice that result in a symmetric splitting. By determining the difference in offset frequencies necessary to reach neighboring symmetric configurations, we can precisely calibrate the phase of the superlattice.

### Disorder lattice

We use a Coherent MBR-110 Ti:Sa laser operated at to create the disorder potential and the perpendicular lattices. This laser is stabilized via its internal reference cavity, which can be externally tuned to change the laser frequency. An additional external temperature stabilized and vibration isolated Fabry-Pérot cavity is used to measure the linewidth and to calibrate the frequency changes necessary to set the phase of the disorder lattice. The phase is always set before the lattice sequence starts and kept constant throughout the experiment.

Note that, strictly speaking, the incommensurate ratio should be an irrational Diophantine number [52] for there to be a sharp localization transition in the Aubry-André model. However, for a realistic finite system, it is sufficient for the actual period of the combined disorder and regular lattices to exceed the system size [53].

### Disorder strength

The disorder strength in the tight-binding limit depends on the lattice depth of the primary lattice , the disorder lattice depth and the wavelength ratio . It is given by [53]

(0) |

where denotes the Wannier function of the unperturbed short lattice, , and and are the depths in terms of of the short and disorder lattices, respectively. Here, is given in units of . We use a primary lattice depth of yielding an undisturbed tunnel coupling of . Combining these numbers, we obtain a conversion factor for

(0) |

## Long-Term evolution

For extremely long hold times (ms or ) we observe an additional decay of the imbalance, along with the total particle number. This decay can be attributed to unwanted additional effects, such as off-resonant light scattering and collisions with background gas atoms, as well as technical imperfections including laser noise, frequency instabilities or vibrations and other pointing instabilities. All of these effects violate the condition of a perfectly closed system and will therefore lead to a vanishing imbalance in the long-time limit. Since the relevant timescale is, however, more than an order of magnitude larger than the longest times used in the experiment, these effects have a negligible impact on the reported observations.

## Trap effects on non-interacting dynamics

Compared to the homogeneous Aubry-André model, the experimental results are modified due to the presence of the 3D harmonic trap. All experimental results are automatically averaged over many inequivalent 1D tubes and, furthermore, the atoms in each tube experience a confining harmonic potential.

In the presence of a longitudinal harmonic trapping potential (trapping frequency Hz), the sharp transition at expected in the homogeneous case [37] becomes smeared out (see Fig. 3), as even in a weak trap the localization length is bounded from above by the finite cloud size. In a sufficiently strong trap, particles situated far away from the trap center will experience offset energies that are bigger than and thereby become localized to their respective side of the trap, even in the absence of any disorder. Nonetheless, although the trap causes some degree of localization even in the absence of disorder, we still observe a clear crossover from trap dominated behavior to disorder-induced localization.

The finite-size of the used lattice beams (waists ) results in an additional small but finite variation of between the central and off-center tubes which increases with the distance from the beam center. More than of the atoms see a less than change in the hopping term. The Rayleigh length of the lattice beams is long enough that there is no appreciable change in tunneling along the tubes. The above variation in mainly gives rise to a dephasing between the dynamics in individual tubes, and is responsible for the more pronounced damping of imbalance oscillations compared to the homogeneous case (see Fig. 12). In the ideal homogeneous system with zero disorder and no variation in , the imbalance would decay according to a Bessel function [54].

Lastly, since the quasi-random potential is produced by superimposing a disorder lattice on top of the primary lattice, this additional potential gives rise to a longitudinal modulation of the hopping rate. This can be modeled as a sinusoidal variation of the same periodicity as the quasi-random potential (i.e. as , where is the phase difference between the modulation of on-site energy and tunneling. For , the amplitude of this modulation is about of the hopping term. This modulation does not change the average hopping in a tube and hence does not change the average imbalance. However, it further increases the damping of the imbalance oscillations.

## CDW preparation sequence

The preparation sequence is illustrated in Fig. 7. To prepare the initial CDW state we start with a 3D band insulating state in a deep long lattice and deep orthogonal lattices. The short lattice is then ramped up in to with the superlattice phase set to . This creates an array of tilted double wells, where the tilt offsets are large enough for all atoms to be loaded into the even wells. This tilt additionally suppresses all tunnelling while we set the interaction strength by ramping the Feshbach field to the desired value in . The long lattice is then switched off in while the disorder lattice is simultaneously ramped up. Finally, the short lattice is ramped down to in , thereby enabling tunneling and initiating the dynamics. The transverse lattices remain deep to decouple individual tubes.

This preparation sequence yields 96(2)% of atoms on even sites in the short lattice. We attribute the small residual number of atoms sitting on odd sites to non-adiabaticities in the splitting process, tunneling events during the preparation sequence, higher band occupancies and detection imperfections.

## Side-resolved detection sequence

The detection sequence to determine the number of atoms on even and odd lattice sites (and hence ) is shown in Fig. 7. Once the evolution in the short lattice with disorder is complete, we increase the short lattice in to to minimize the tunneling rate. The long lattice is then ramped up to in , while the disorder lattice is simultaneously ramped down. With the spatial distribution of atoms now frozen in the tilted superlattice, we ramp the Feshbach field to the non-interacting () point in , to ensure that there are no unwanted interaction effects during the subsequent band-mapping and imaging process. Next, the long lattice depth is further increased to in , transferring atoms from the second band into the third band of the superlattice via an avoided crossing [41]. Any atoms on odd sites of the short lattice (the higher side of the tilted superlattice double wells) are now in the 3rd band, while atoms from even sites of the short lattice remain in the lowest band of the superlattice. After of hold time, the short lattice is ramped to in , which is slow enough to be adiabatic with respect to inter-band transitions. Since the ordering of the bands remains unchanged during this ramp, atoms originally on even (odd) sites in the short lattice end up in the 1st (3rd) band of the long lattice. Finally, we implement a bandmapping sequence followed by TOF and absorption imaging, from which we can count the number of atoms in each band directly, providing us with and . A few atoms remaining in the second band due to an imperfect transfer are counted together with the third band.

## Doublon fractions

In order to characterize our initial state, we measure the fraction of atoms forming doubly occupied sites, termed the doublon fraction, by converting doublons into molecules by crossing the Feshbach resonance [40]. Experimentally this involves loading into a 20 deep long lattice and deep orthogonal lattices as per the standard sequence with loading interactions given by , then ramping the long lattice to 40(1) in s to supress hopping in all three dimensions, which preserves the atomic distribution for the duration of the doublon detection sequence. The Feshbach field is then ramped down to either 206.3G (= -117(2), still above the resonance) or 198.3G on the molecular side of the resonance. When the Feshbach field ramp crosses the resonance, any pairs of atoms on the same lattice site will be converted into Feshbach molecules [55]. Since the molecular transition is out of resonance of the imaging light, we can determine the doublon fraction directly by comparing the number of atoms detected at the two different Feshbach ramp endpoints. The number of atoms is measured using standard absorbtion imaging after 6ms TOF. Fig. 8 shows the doublon fraction plotted against the loading interaction during the lattice ramp up. Error bars show the standard deviation. We quote the loading interactions in units of rather than , as the lattice configuration at the end of the loading differs from the one during the evolution. In the lattice configuration used for time evolution of the CDW, a scattering length of corresponds to .

We also investigate the effect of the initial temperature on the doublon fraction. The initial temperature is changed by reducing the duration of the final sympathetic cooling stage and increasing the final dipole evaporation depth. This has the effect of making the cooling less efficient while keeping the final atom number relatively constant. Fig. 9 shows the doublon fraction plotted as a function of the temperature of the initial state, expressed in units of the Fermi temperature , for three different loading interactions. In the case of loading with attractive interactions or in the non-interacting case, we observe higher temperatures to reduce the resulting doublon fraction, since the increased entropy and average kinetic energy per particle reduce the average density and render the formation of doublons less favourable. For atoms which are loaded repulsively, there are very few doublons for any temperature and thus no real dependence is observable within our error bars.

## Imbalance vs energy density

We investigate the effects of both loading interactions and initial temperature on the imbalance in the many-body localized regime at , close to the localization transition. While changing the interactions during loading predominantly affects the number of doublons for a given temperature [56], changing the temperature of the initial state affects both doublon fraction and average density. Fig. 10 shows the effect of loading interactions on the imbalance for three different interactions during evolution.

For vanishing or weak interactions during evolution, we observe no influence of loading interactions, highlighting the robustness of many-body localization in a lattice with respect to energy density. However, for strong repulsive interactions during evolution, we see a significant effect of the loading interactions on the imbalance . Increasing the number of doublons by loading attractively increases the imbalance, since the additional energy scale of the doublons brings the system much deeper into the localized phase. In contrast, for repulsive loading, doublons are strongly suppressed. In the limit of strong interactions (during evolution) and the absence of doublons, interacting fermions can be exactly mapped to spinless (non-interacting) fermions (see below). This is precisely what we observe in our experiment for strong repulsive loading interactions and thus vanishing doublon fraction: the imbalance observed for strong interactions during evolution converges towards the case of a non-interacting evolution.

For intermediate interactions , we cannot detect any notable variation of the imbalance within errors for all loading interactions. Additionally, in Fig. 11 we investigate the effect of loading temperatures on . Fig. 11A shows the case of attractive loading at , while Fig. 11B shows the case of repulsive loading interactions . In all cases, increased temperatures reduce the effect of interactions due to the decreasing average density, connecting the interacting cases asymptotically to the behaviour of non-interacting fermions. In the case of attractive loading, higher temperatures in addition reduce the number of doublons, which for strong repulsive interactions during the evolution reduces the localization. For intermediate interactions, the effect of initial temperature is less pronounced, but for both loading interactions increasing initial temperatures connect the minimum of localization (Fig. 6) to the case of no interactions , see Fig. 11.

## Fitting of the imbalance vs interaction curves

We extract the minima of the imbalance as a function of interactions by employing the following phenomenological fitting function:

(0) |

Here , , , and offset are free fitting parameters. We have implemented the symmetry to the function by fixing the center of both Gaussians to zero. The positive amplitude Gaussian with a smaller width captures the reduction of the imbalance up to the minimum, while the decreasing contribution of the negative amplitude Gaussian resembles the increase for strong interactions. The minima of this function are used to generate the green dot-dashed lines in Fig. 4A.

## Time evolution of the CDW in a non-interacting model

We consider an initial state with density wave order, i.e.

(0) |

where is the total number of sites, () is the number of fermions on even (odd) sites, and represents a trace over the initial density matrix. This state is evolved with a (spinless) non-interacting fermionic Hamiltonian

(0) |

The disorder could represent true randomness or a quasi-periodic potential as in the experiment. To obtain the density wave order at later times we need to compute the expectation values of the local densities at later times, i.e. . We can write the density operators explicitly as

(0) |

where is the non-interacting propagator from point to point and is the first-quantized Hamiltonian. Taking the trace over the density matrix leaves only the terms with yielding

(0) |

where we have assumed an initial state with particles only on even sites, i.e. .

Therefore the CDW or (imbalance) at time is

(0) |

which can be calculated numerically

## Details of system modelling via exact diagonalisation

In order to calculate the red line in Fig. 3, we employ Eq. (Time evolution of the CDW in a non-interacting model) to calculate the time evolution of the initial CDW in a single tube in the presence of the harmonic trap. The initial CDW is modeled as a collection of randomly distributed localized particles, whose distribution follows a Gaussian with and matches the initial imbalance of observed in the experiment. We extract the stationary imbalance by averaging the imbalance in the evolution time interval from to . In addition, the calculations are averaged over twelve disorder realizations corresponding to twelve equally spaced phases in Eq. (Introduction).

In order to correctly model the dynamic evolution shown in Fig. 12, we additionally integrate the results over many individual 1D tubes with tube-dependent hoppings in order to fully model the experimental situation. The widths and amplitudes of the corresponding Gaussian distributions were chosen in accordance with the experimental cloud sizes of and that were obtained via in-situ imaging. We included both the hopping perturbation resulting from the disorder potential as well as the effect of the finite waist of the lattice beams. The strength and phase difference of the disorder induced hopping modulation are extracted from fits to the exact position dependent hopping rates.

## Stationary imbalance value in the localized phase

We can give an analytic estimate for the steady state value at long times in the localized regime. Time averaging the wave function of a particle originating from a single site should lead to a form which is exponentially decaying over a localization length . Hence, we can approximate , where the prefactor ensures that the initial site of the wave function is an even site. Note that this form holds for , i.e. a localization length much larger than a lattice constant . Inserting this expression into Eq. (Time evolution of the CDW in a non-interacting model) for the density wave order yields

(0) |

where in the last equality we changed to a sum over relative and center of mass coordinates . The second term vanishes and we find (for large )

(0) |

This formula establishes the connection between the CDW stationary value and the localization length close to the transition in a non-interacting system.

## Mapping from hard-core to non-interacting fermions

We consider the dynamics of the (spin-) Hubbard model in the limit of infinite (on-site) when the initial state has no doubly-occupied sites. In this case, the time evolution can be mapped to that of non-interacting spinless fermions.

In the limit the Hilbert sub-space with no doubly-occupied sites is closed under unitary time evolution. Double occupancies cannot be generated in the dynamics while conserving energy because they incur an infinite energy cost. Hence, we can replace the interaction with a hard-core constraint implemented by the operator , which projects out the doubly-occupied subspace,

(0) |

Let us denote the ordering of the fermion spins from left to right on the chain with , where is the total fermion number. is an invariant of the dynamics because in order to exchange the spins of two particles with opposite spin orientation, the particles need to meet on the same site, which is prohibited by the hard-core constraint. Now, suppose we start the time evolution with an initial state having a definite spin ordering . Since this is an invariant of the dynamics, we can describe the time evolution within this sector with a Hamiltonian in which the spin ordering enters as a parameter. The Hilbert space in which this Hamiltonian acts can be faithfully represented by the Fock space of spinless fermions , because given the fermion positions their spins are uniquely determined by the configuration which parameterizes the space. Hence we can write the Hamiltonian which acts within the subspace as a Hamiltonian of spinless fermions

(0) |

Note that the projection of doubly-occupied sites is automatically taken into account using the representation with spinless fermions.

The mapping of the time evolution within a given spin-ordering sector can also be understood as a result of a non-local unitary change of basis

(0) |

where denotes the particle number to the left of point , and

(0) |

In other words, the transformation counts the number of particles to the left of point to find its order in the chain . It then rotates the spin by around the axis if in the configuration the -th spin is . The transformed fermions are then all and we can drop the spin index.

The dynamics is thus described by the same Hamiltonian of non-interacting spinless fermions regardless of the initial ordering of the spins . In particular, the time dependence of the CDW order can be generated by unitary evolution of a spinless non-interacting fermion model and is independent of the initial spin ordering.

## Time evolution of the CDW in an interacting model

We consider the following ensemble of many-body states as initial states for the numerical calculation:

(0) |

where describes the vacuum and . In order to generate the ensemble of initial states, the occupations are chosen from an appropriate joint probability distribution such that the average site occupation and the average doublon occupation are fixed. This state is time evolved with the interacting Hamiltonian

(0) |

using the time-evolving block decimation (TEBD) [57, 58] on matrix product states implemented with the ITensor Library [59].

For our numerical simulations we use a quasi-periodic potential with as in the experiment. However, we do not take into account the periodic trap, but rather consider a system with sites and open boundary conditions. When simulating the time evolution for a given set of system parameters, i.e. , , site occupation, and doublon density, we average over both the initial state and the disorder realization (through sampling of the random phase ). The total number of runs for a given set of parameters is , i.e. about 10 initial states for every disorder realization. For the singular value decomposition, we use a cutoff of and allow for a maximum bond dimension of . Using a time step of , this leads to a maximum truncation error .

## Fluctuations of the imbalance

In the main text we presented TEBD (i.e. DMRG) calculations showing logarithmic growth of the entanglement entropy with time, which is characteristic of the many-body localized state. The entanglement entropy, however, is a non-local quantity and therefore not directly measurable. In this section we show that the same information about the many-body localized state can be inferred from the temporal fluctuations of the measured imbalance (CDW order parameter) around its stationary value.

We focus on the system deep in the localized phase, where it can be described using an effective model of local integrals of motion. For simplicity we initialize the system in a CDW state with exactly one particle on every even site. We now partition the lattice into two-site dimers and, as a starting point, take the dimers to be disconnected. We define and , which within the accessible Hilbert space act like Pauli spin operators. In terms of these operators the imbalance (or CDW order) reads

(0) |

For a system with strongly localized particles (i.e. ) we will continue to use an effective model of decoupled dimers, where from now on denotes the number of dimers. We take into account the tail of the single particle wave-function, which extends outside of the dimers, only through the interaction it induces between the local spin operators. Hence an effective model is

(0) |

where , are the integrals of motion and the interaction is given by .

For the non-interacting case, only the first term of the effective model is non-zero, and we directly find the expectation value of the spin operator on site

(0) |

where and . Thus, we find the time-averaged stationary value of the imbalance is . The temporal fluctuations around this average are

(0) |

where represents averaging over a long time window . We conclude that without interactions the fluctuations of the imbalance remain on average constant in time with .

When interactions are present the local contribution to the density wave is affected by the time evolution of the other spins. The pseudo spin gradually becomes entangled with a growing number of other spins. Specifically, because the interactions decay exponentially with distance, the number of spins that become significantly entangled with a given spin at time is . We will see below how this affects the fluctuations of and ultimately contributes to the temporal fluctuations of the order parameter .

Let us calculate the reduced density matrix of the pseudo spin at site . The system starts the time evolution in a product state

(0) |

with and . We can formally write the time-dependent density matrix as

(0) |

Here, . We can now trace out all but one site to obtain the reduced density matrix of a single pseudo-spin at site in the basis of the eigenvalues of :

(0) | |||||

(0) | |||||

(0) |

where . and represent states of all other sites except site . The number of frequencies involved is roughly , i.e. all possible interactions between the spins that are significantly entangled with the observed spin at time . More precisely, the number of frequencies measures the size of the entangled Hilbert space at the observation time.

To find the time dependence of the CDW we have to transform back from the basis of eigenvalues to the basis of . This is achieved with a rotation around the axis by

(0) |

Hence, the fluctuations of the local imbalance between an even site and neighboring odd site behaves as

(0) |

while the fluctuation of the global imbalance (CDW order) is further suppressed by a factor of

(0) |

We see that the decay of the fluctuations is intimately connected to the growth of the entanglement entropy, as it also ‘measures’ the number of spins coupled through the interactions after time .

The above analysis is of a simplified effective model. However, the main conclusions are supported by direct simulation of the Hubbard model on the quasi-periodic lattice using time dependent DMRG. Fig. 13A shows the temporal noise of the imbalance as a function of the time for different values of the interaction strength . The fluctuations are measured by averaging them over a time window of around the time . The results fit well to a power law decay. Fig. 13B compares the fitted exponent to the slope of the logarithmic growth of the entanglement entropy showing the direct correspondence between the two.

Thus, we conclude that measurements of the temporal fluctuations of the CDW order provide a viable experimental route to determine the bare localization length and distinguish the many-body localized state from an Anderson localized state of non-interacting particles. In the present experiment, this is, however, not possible, as the unavoidable averaging over many tubes suppresses the fluctuations below the detection limit after only few oscillations. For future experiments, a single tube, or even single-site resolution would be desirable to overcome this limitation. We also note that the temporal fluctuations of the expectation value are different from shot-to-shot fluctuations at a given time which reflect the quantum uncertainty of the observable and would be finite even in the infinite-time limit.

## References:

- [49]
- [50] C. A. Regal, C. Ticknor, J. L. Bohn, D. S. Jin, Creation of ultracold molecules from a Fermi gas of atoms. Nature 424, 47 (2003).
- [51] C. Klempt, T. Henninger, O. Topic, J. Will, W. Ertmer, E. Tiemann, J. Arlt, K-Rb Feshbach resonances: modeling the interatomic potential. Phys. Rev. A 76, 020701(R) (2007).
- [52] S. Y. Jitomirskaya, Metal-insulator transition for the almost Mathieu operator. Ann. Math. 150, 1159 (1999).
- [53] M. Modugno, Exponential localization in one-dimensional quasi-periodic optical lattices. New J. Phys. 11, 033023 (2009).
- [54] D. Pertot, A. Sheikhan, E. Cocchi, L. A. Miller, J. E. Bohn, M. Koschorreck, M. Köhl, C. Kollath, Relaxation dynamics of a Fermi gas in an optical superlattice. Phys. Rev. Lett. 113, 170403 (2014).
- [55] T. Köhler, K. Góral, P. S. Julienne, Production of cold molecules via magnetically tunable Feshbach resonances. Rev. Mod. Phys. 78, 1311 (2006).
- [56] L. Hackermüller et. al., Anomalous Expansion of Attractively Interacting Fermionic Atoms in an Optical Lattice, Science 327, 1621 (2010).
- [57] G. Vidal, Efficient Simulation of One-Dimensional Quantum Many-Body Systems, Phys. Rev. Lett. 93, 040502 (2004).
- [58] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96–192 (2011).
- [59] http://itensor.org