Disorder-controlled relaxation in a three-dimensional Hubbard model quantum simulator

W. Morong ,1,* S. R. Muleady ,2,3 I. Kimchi,2,3,† W. Xu,1,‡ R. M. Nandkishore,3,4 A. M. Rey ,2,3 and B. DeMarco1,§ 1Department of Physics, IQUIST, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA 2JILA, National Institute of Standards and Technology and Department of Physics, University of Colorado, Boulder, Colorado 80309, USA 3Center for Theory of Quantum Matter, University of Colorado, Boulder, Colorado 80309, USA 4Department of Physics, University of Colorado, Boulder, Colorado 80309, USA


I. INTRODUCTION
Strong disorder and interactions are known to give rise to the celebrated Anderson and Mott metal-insulator transitions. In an Anderson insulator, a random spatial potential localizes noninteracting particles through destructive interference [1,2], while for a unit-filled Mott insulator, strong repulsive interactions create an energy gap that prevents particle motion [3]. The combined presence of disorder and interactions in many physical systems poses the open challenge of understanding the interplay between these vastly different localization mechanisms [4,5]. The development of highly tunable and isolated quantum simulators, such as ultracold atoms trapped in optical lattices, has created new opportunities to experimentally study this long-standing problem using a minimal model combining both elements: the disordered Fermi-Hubbard model (DFHM) [6].
A potential result of combined disorder and interactions in isolated systems is many-body localization (MBL), in which relaxation to thermal equilibrium is prevented by sufficiently strong disorder [7][8][9][10]. Despite a concerted theoretical effort in recent years and several experimental results for one-dimensional chains [11][12][13][14], many questions still remain regarding the nature of this phenomenon. This shortcoming is especially true for systems with more than one spatial dimension, although initial experimental studies with atoms in two-and three-dimensional optical lattices have also observed slow dynamics consistent with MBL [15][16][17].
Another mechanism that can suppress relaxation emerges in the strongly interacting regime of the DFHM: the formation of quasibound doubly occupied sites (i.e., doublons), which slowly decay in a clean lattice through high-order processes that generate many low-energy excitations [18,19]. The effect of disorder on doublon relaxation is largely unexplored. Reconciling the interplay of slow dynamics caused by doublon binding and disorder-induced localization is critical to obtaining a more complete understanding of thermalization in the DFHM, including the possibility of MBL. However, advancing this frontier demands exploring the highly nontrivial regime characterized by comparable disorder and interaction energies.
Here, we investigate how strong interactions and disorder compete and cooperate to affect the far-from-equilibrium dynamics of the three-dimensional (3D) DFHM. Using a quantum simulator of fermionic atoms in an optical lattice, we perform measurements of doublon relaxation following (1)], which involves tunneling with amplitude t, doublon formation with interaction energy U , and a local disordered energy offset with an exponential distribution ρ( ) characterized by disorder strength . (c) While the full dynamics of the system involve the interplay of doublons and single atoms, the case of an isolated doublon is useful for intuition about the dynamical regimes. For weak disorder (compared to U ), the decay of doublons into pairs of single atoms is suppressed by the gap between the doublon and single-particle energies. Moderate disorder can enhance the decay rate by creating resonant pairs of sites with an energy difference comparable to U . However, sufficiently strong disorder can hinder decay by suppressing diffusion to these resonant sites.
an interaction quench and use the resulting decay times to characterize the system behavior. Exploring the parameter regime from interaction-dominated to disorder-dominated behavior, we are able to classify the dynamics in terms of two distinct regimes: disorder-suppressed relaxation at strong disorder, and disorder-enhanced relaxation at weaker disorder. The latter effect has yet to be observed in a quantum simulator platform and may be related to disorder-driven insulator-metal transitions observed in certain correlated materials [20][21][22]. We compare our results to beyond-mean-field numerical simulations of the quantum dynamics, which capture many features of our observed results and suggest the creation of resonances in the lattice by the disorder as the physical origin of this disorder-enhanced regime. We further supplement this picture by developing a simple phenomenological model that incorporates both disorder-enhanced and disorder-suppressed mechanisms for doublon decay.

II. MEASURING AND SIMULATING RELAXATION
We employ an optical lattice experimental platform ( Fig. 1) described in previous work [15]. Two spin states (denoted | ↑ and | ↓ ) of fermionic 40 K atoms are trapped in a cubic lattice superimposed with 532 nm optical speckle disorder. This system realizes the DFHM with confinement ( Fig. 1): Here, σ indexes the two spin states, t i j is the tunneling energy between sites i and j (restricted to nearest-neighbors, indicated by i, j ), U i is the on-site interaction energy, i is the local disordered energy offset, ω is the harmonic confinement, m is the atomic mass, and r i is the distance of site i from the trap center. The single-particle bandwidth is 12t in the absence of disorder. The applied speckle potential creates disorder in the t, U , and terms, leading to distributions of these Hubbard parameters with widths that depend on the optical power [23,24]. We characterize the disorder strength by , which is approximately equal to the standard deviation of the distribution [ Fig. 1(b)]. The influence of spatial correlations in the disorder potential is weak, since the correlation length of the speckle field is smaller than two lattice spacings along every lattice direction. Because of the harmonic confinement, all measurements are averaged over a density profile that varies from an estimated occupancy n = 0.5 at the center of the trap to zero at the edges of the system.
To probe far-from-equilibrium doublon dynamics, we measure the population of atoms in doubly occupied sites following a quench in which the interactions are reversed from attractive to repulsive using a Feshbach resonance. Before the quench, the gas is in equilibrium with an energetically favorable doublon population. The equilibrium temperature is estimated (from measurements of the entropy before the lattice and disorder are turned on) as T /12t = 5.3, resulting in an initial state that is at high temperature with respect to the ground band. After the quench, the doublons become excitations that can decay by breaking apart into a pair of single atoms ("singles"). The atomic doublon population is allowed to evolve in a disordered lattice by turning on the optical speckle field following the interaction quench. After a variable time, the doublon population is measured by mapping each doublon to a tightly bound Feshbach molecule and selectively transferring the | ↓ atom in each molecule to an ancillary spin state using an rf sweep [25]. Alternatively, we can selectively transfer and image only atoms from singles, which allows us to separate doublon decay from overall number loss.
In all regimes, the doublon population decreases following the quench with a rate sensitive to the disorder strength. Typical results at different disorder strengths are shown in Fig. 2. To quantify the decay, we fit the data to a model that describes exponential doublon decay with a time constant τ towards an equilibrium value, which itself slowly decays because of overall particle loss with τ loss = 2100h/t (see the Supplemental Material [25] for model details). While we find that this fit provides a reasonable characterization of the decay timescale, the functional form of the relaxation is unknown beyond the clean limit. We therefore turn towards numerics to provide an interpretation of the timescale with disorder present.
The large scale and dimensionality of our system, as well as the far-from-equilibrium nature of the dynamics, preclude exact numerical studies and commonly employed approximate techniques, such as density matrix renormalization group (DMRG) or diagonalization methods. Furthermore, the fermionic sign problem forbids use of quantum Monte Carlo methods. We therefore develop a numerical method-a generalized discrete truncated Wigner approximation (GDTWA) [26][27][28][29]-to simulate the relaxation process. The GDTWA All data are taken at U/12t = 1.8 (corresponding to a 12E R lattice depth in the experiment). Sample traces of the doublon population vs time for different disorder strengths are shown for experimental measurements (points); the error bars give the statistical standard error of the mean for averaging over multiple measurements. The fits used to extract τ are shown using solid lines. The dashed lines denote the fit values for the equilibrium doublon population, which is typically 10% of the total atom number [25] and which decays over time due to overall number loss. The corresponding GDTWA simulations [25] are displayed as shaded regions that show the standard error of the mean from disorder and trajectory averaging. The horizontal axis for (d) (the strongest disorder case) is compressed to bring the slow decay into view. At the highest disorder the doublon lifetime τ becomes closest to the overall number loss τ loss , potentially affecting the fitted value for the equilibrium doublon fraction [25].
approach invokes a factorization of the density matrixρ of the many-body system over individual lattice sites i,ρ = iρ i . Furthermore, random sampling of initial conditions from a discrete semiclassical phase space accounts for quantum noise. GDTWA thus provides a numerically tractable method for examining the role of quantum correlations in the system dynamics.
More specifically, we can fully describe the reduced statê ρ i at lattice site i (with a four-dimensional local Hilbert space spanned by basis states {|↑↓ , |↑ , |↓ , |0 }) by a 16dimensional vector λ i , so thatρ i = α λ α iˆ α i for a complete basis of local observables {ˆ α i }. Inserting the product state ansatz into the von Neumann equation ∂ tρ = −ih[Ĥ ,ρ] results in a set of nonlinear differential equations describing the evolution of each λ i , and which may be numerically integrated to obtain the mean-field dynamics. While this ansatz retains full information regarding the strong on-site Hubbard correlations, such a solution describes a product state at all instances in time and therefore neglects the important crosssite quantum correlations responsible for doublon decay. In fact, for an initial product state of definite particle number and spin, the mean-field solution results in no dynamics for this system, highlighting the importance of quantum correlations in enabling doublon decay.
To capture the buildup of quantum correlations between sites, we instead describe the initial value of each λ i as a probability distribution over a discrete phase space that samples over all allowed values for each observable in our local basis [25]. Averaging independently evolved sets of randomly sampled initial conditions from this distribution yields a highly nontrivial solution for the dynamics (owing to the nonlinear nature of the dynamical equations) and describes the evolution to a correlated state exhibiting entanglement [27,29]. This approximation, while only rigorously valid at short times, has demonstrated the ability to provide accurate results in generic spin models at longer times and properly capture quantum thermalization of local observables [27,30,31]. Here, we adapt the GDTWA to model DFHM dynamics and we provide benchmarks for its applicability in the Supplemental Material [25].
The GDTWA results, shown as traces in Fig. 2, approximately capture the observed changes in doublon population decay with applied disorder. With the notable exception of /U = 0.067, the GDTWA real-time dynamics generally reside within experimental error bars. For small disorder where doublon decay is expected to proceed via a large number of high-order processes [18,19], we expect that the nature of our approximation will be unable to fully account for all involved decay channels. The GDTWA dynamics also allow us to determine, within the assumptions of the technique, the extent to which τ accurately characterizes the relaxation timescale. As illustrated in Fig. 2, we observe qualitative differences between the measurements and simulations, which tend to exhibit a faster initial decay (over timescalesh/t) before transitioning to a slower decay on longer timescales that is well fit by an exponential decay. Nonetheless, by constraining the initial doublon fraction in the fitted exponential decay to the measured value, we are able to obtain a reasonable exponential fit of the entire dynamical curve that is used to extract a consistent relaxation timescale. The resulting effective timescale incorporates both the initial nonexponential features and the subsequent exponential decay [25].
To account for the spatially varying trap density and variations in the initial preparation, which affects the equilibrium doublon density, we choose initial single densities for the simulations that yield best fits to the experimental data. These best-fit single densities in the simulations are generally consistent with the experimentally measured single density at the trap center, except at small disorders where GDTWA does not fully capture the decay processes associated with the clean system and thus results in a much larger equilibrium doublon density for all initial single densities, as evidenced in Fig. 2(a) (see also the Supplemental Material [25]). However, we have also confirmed that the relaxation timescale is not strongly influenced by the average density nor the spatial profile of the gas in the regime of interest [25]. These observations support characterizing the doublon relaxation dynamics by a single parameter τ : τ primarily depends on the interplay of disorder with interactions and tunneling and is insensitive to parameters that may fluctuate or have some uncertainty in the experiment.

III. ANALYSIS OF DYNAMICAL REGIMES
The dependence of τ on disorder is shown in Fig. 3; the qualitative features are shared by experiment and numerics. Strikingly, we find that applying disorder first causes the relaxation time τ to rapidly decrease. For the clean system, τ is substantially longer than the single-particle tunneling timeh/t, which is consistent with previous studies that identified doublons as repulsively bound pairs [18,19,[32][33][34]. However, disorder causes τ to decrease to a minimum value comparable to the tunneling timeh/t at a disorder value near ∼ U/2. While the complementary phenomenon of interaction-driven delocalization has been observed in similar systems [11,15,35], here we observe a disorder-driven increase in relaxation in a quantum simulator with a high degree of isolation and tunability.
As is increased beyond U , τ eventually increases, growing by over two orders of magnitude at the strongest disorder we can apply in the experiment. The separation into two dynamical regimes (distinguished by the slope of τ with ) combined with the crossover at ∼ U suggest that the dynamics are controlled by competing mechanisms arising from interactions and disorder.
We can understand these dynamical regimes using a minimal model of diffusing doublons in a disordered environment. In this model, the interchange between doublon-hole pairs and pairs of opposite-spin singles is controlled by a set of reaction-diffusion (RD) equations. The RD model [25] augments the classical continuum diffusion equation for each particle species with a source term that converts a doublonhole to a single-single combination, only when the local parameters allow this process to be resonant. In particular, our model requires the local energy difference arising from the speckle disorder to lie within a window of width ∼t around U (see Fig. 1). The doublon diffusion coefficients are taken to decrease monotonically with disorder strength, as increasingly large local energy differences will inhibit doublon transport throughout the lattice.
The RD model gives the decay rate 1/τ as a product of the effective diffusion rate D eff and the probability P reaction for a conversion between a doublon-hole pair and two singles (↑ and ↓): The probability of a reaction per site encountered (derived in the Supplemental Material [25]) is then expressed via a product of the exponential speckle disorder amplitude at the conversion energy U , and a term corresponding to the width of a conversion resonance, as follows, where the additional constant p is a small (i.e., of order 0.01) parameter denoting the probability of reactions in the clean limit, and z is the lattice coordination number. The effective diffusion coefficient D eff gives the rate (linear in time) at which new sites are sampled by any particular doublon. We expect D eff to decrease rapidly with disorder since it is controlled by the doublon diffusion constant. Even in the asymptotically localized phase predicted at large disorder (where doublon diffusion vanishes), D eff should be supplemented by an additional small velocity that allows sampling of sites within the finite localization length. The exact dependence of diffusion D eff on disorder /U is not qualitatively important as long as D eff decreases monotonically to a very small or vanishing value. In Fig. 3, we choose a diffusion coefficient that decreases exponentially with disorder. For this functional form, the best-fit value for D eff is of order t/h at low disorder (4 ± 1 t/h) and is reduced by a factor of over 100 (to 0.02 ± 0.008 t/h) at the highest . We find that an algebraically decaying diffusion coefficient yields quantitatively similar decay rates at all measured disorder strengths. The combined expression for the decay rate 1/τ within this reaction-diffusion model explains the two relaxation regimes as follows.
In regime 1 ( U ), increasing disorder leads to decreasing decay times. This regime is controlled by the reaction probability P reaction . In the clean limit, the reaction rate is greatly suppressed by the strong interaction energy U , and doublon decay occurs slowly. However, as the disorder strength increases, the tail of the local energy distribution allows for a finite probability of an adjacent site energy difference of order U . By compensating the interaction energy U through this mechanism, the disorder produces a quantum resonance that allows the decay reaction to proceed with high probability, thereby resulting in fast doublon relaxation.
In regime 2 ( U ), increasing disorder leads to increasing decay times. This regime is controlled by the diffusion constant D eff , which becomes heavily suppressed by localization effects produced by the large disorder. To a lesser extent, the gradual suppression of the reaction probability due to fewer resonances also contributes to the behavior in this regime. Recent theoretical work suggests that, aside from rare region fluctuations, diffusion should vanish at sufficiently large disorder and signal an asymptotic many-body localized phase [10]. The diffusion in the 3D interacting system is an unknown function of the disorder ratio to the bandwidth /12t and interactions /U , but dimensional considerations suggest a reduction when the ratios become of order unity, consistent with the observations Fig. 3. While the experimental data cannot distinguish between D eff = 0 and D eff saturating to a finite small value, we can robustly conclude that diffusion must be highly inhibited in regime 2.

IV. CONCLUSION
Quantum many-body systems involving interactions and disorder can exhibit emergent behavior that, especially in two and three dimensions, is challenging to predict from first principles. By experimentally constructing a DFHM quantum simulator in a 3D optical lattice, we are able to study the dynamical out-of-equilibrium behavior of doublons across vastly different scales of disorder and interactions. Remarkably, we find that constructing a simple model for this extraordinarily complicated system is possible using reaction-diffusion equations that are analytically solvable. We devise a numerical approach that captures the same physical effects and shows quantitative agreement with the experiment.
Intriguing open questions remain to be explored. In the intermediate disorder regime ( ∼ U ), the observed fast doublon relaxation may be related to the behavior of "bad metals" [36,37], which can be characterized by a lack of conserved excitations [38]. Because the rapid doublon decay can be understood as a consequence of disorder disrupting the gap in the single-particle spectrum, spectroscopic measurements in this regime could identify a disorder-created pseudogap in the density of states [20,21]. As a practical tool, the fast disordermediated thermalization that we observe also suggests an approach to avoid challenges in the adiabatic preparation of strongly correlated atomic gases [39]. In the strong disorder regime, our measurements imply a suppression of particle transport as disorder is increased. However, we are unable to distinguish whether this behavior is a manifestation of slow diffusion or a signature of asymptotic many-body localization. Further experiments could clarify this difference and extend the results of Ref. [15] to the U/12t > 1 regime.