Disorder-controlled relaxation in a 3D Hubbard model quantum simulator
W. Morong, S.R. Muleady, I. Kimchi, W. Xu, R.M. Nandkishore, A.M. Rey, B. DeMarco
DDisorder-controlled relaxation in a 3D Hubbard model quantum simulator
W. Morong,
1, 2, ∗ S. R. Muleady,
3, 4, ∗ I. Kimchi,
3, 4
W. Xu,
1, 5
R. M. Nandkishore,
4, 6
A. M. Rey,
3, 4 and B. DeMarco Department of Physics, IQUIST, University of Illinois, IL 61801, USA. Present address: Joint Quantum Institute and Joint Center for Quantum Information andComputer Science, University of Maryland and NIST, College Park, Maryland 20742, USA JILA, National Institute of Standards and Technology and Departmentof Physics, University of Colorado, Boulder, CO, 80309, USA Center for Theory of Quantum Matter, University of Colorado, Boulder, CO, 80309, USA Present address: Massachusetts Institute of Technology, Department of Physics, Cambridge, MA 02139 Department of Physics, University of Colorado, Boulder, CO 80309, USA
Understanding the collective behavior of strongly correlated electrons in materials remains a cen-tral problem in many-particle quantum physics. A minimal description of these systems is providedby the disordered Fermi-Hubbard model (DFHM), which incorporates the interplay of motion in adisordered lattice with local inter-particle interactions. Despite its minimal elements, many dynam-ical properties of the DFHM are not well understood, owing to the complexity of systems combiningout-of-equilibrium behavior, interactions, and disorder in higher spatial dimensions. Here, we studythe relaxation dynamics of doubly occupied lattice sites in the three-dimensional (3D) DFHM usinginteraction-quench measurements on a quantum simulator composed of fermionic atoms confined inan optical lattice. In addition to observing the widely studied effect of disorder inhibiting relaxation,we find that the cooperation between strong interactions and disorder also leads to the emergenceof a dynamical regime characterized by disorder-enhanced relaxation. To support these results, wedevelop an approximate numerical method and a phenomenological model that each capture theessential physics of the decay dynamics. Our results provide a theoretical framework for a previ-ously inaccessible regime of the DFHM and demonstrate the ability of quantum simulators to enableunderstanding of complex many-body systems through minimal models.
I. INTRODUCTION
Strong disorder and interactions are known to give riseto the celebrated Anderson and Mott metal–insulatortransitions. In an Anderson insulator, a random spa-tial potential localizes non-interacting particles throughdestructive interference [1, 2], while for a unit-filled Mottinsulator, strong repulsive interactions create an energygap that prevents particle motion [3]. The combinedpresence of disorder and interactions in many physicalsystems poses the open challenge of understanding the in-terplay between these vastly different localization mech-anisms [4, 5]. The development of highly tunable andisolated quantum simulators, such as ultracold atomstrapped in optical lattices, has created new opportuni-ties to experimentally study this long-standing problemusing a minimal model combining both elements: thedisordered Fermi-Hubbard model (DFHM) [6].A potential result of combined disorder and inter-actions in isolated systems is many-body localization(MBL), in which relaxation to thermal equilibrium is pre-vented by sufficiently strong disorder [7–10]. Despite aconcerted theoretical effort in recent years and severalexperimental results for one-dimensional chains [11–14],many questions still remain regarding the nature of thisphenomenon. This shortcoming is especially true for sys-tems with more than one spatial dimension, althoughinitial experimental studies with atoms in two and three- ∗ Authors W.M. and S.R.M. contributed equally to this work. dimensional optical lattices have also observed slow dy-namics consistent with MBL [15–17].Another mechanism that can suppress relaxationemerges in the strongly interacting regime of the DFHM:the formation of quasi-bound doubly occupied sites (i.e.,doublons), which slowly decay in a clean lattice throughhigh-order processes that generate many low-energy exci-tations [18, 19]. The effect of disorder on doublon relax-ation is largely unexplored. Reconciling the interplay ofslow dynamics caused by doublon binding and disorder-induced localization is critical to obtaining a more com-plete understanding of thermalization in the DFHM, in-cluding the possibility of MBL. However, advancing thisfrontier demands exploring the highly non-trivial regimecharacterized by comparable disorder and interaction en-ergies.Here, we investigate how strong interactions and dis-order compete and cooperate to affect the far-from-equilibrium dynamics of the three-dimensional DFHM.Using a quantum simulator of fermionic atoms in an op-tical lattice, we perform measurements of doublon relax-ation following an interaction quench and use the result-ing decay times to characterize the system behavior. Ex-ploring the parameter regime from interaction-dominatedto disorder-dominated behavior, we are able to classifythe dynamics in terms of two distinct regimes: disorder-suppressed relaxation at strong disorder, and disorder-enhanced relaxation at weaker disorder. The latter effecthas not previously been observed in a quantum simulatorplatform and may be related to disorder-driven insulator-metal transitions observed in certain correlated materials[20, 21]. We compare our results to beyond-mean-field a r X i v : . [ c ond - m a t . qu a n t - g a s ] J un FIG. 1. Experimental setup. i: A nonequilibrium doublon population for atoms in two spin states is prepared in a 3D cubiclattice using an interaction quench. ii: The atoms are described by the disordered Fermi-Hubbard model (Eq. 1), whichinvolves tunneling with amplitude t , doublon formation with interaction energy U , and a local disordered energy offset (cid:15) withan exponential distribution ρ ( (cid:15) ) characterized by disorder strength ∆. iii: While the full dynamics of the system involve theinterplay of doublons and single atoms, the case of an isolated doublon is useful for intuition about the dynamical regimes. Forweak disorder (compared to U ), the decay of doublons into pairs of single atoms is suppressed by the gap between the doublonand single-particle energies. Moderate disorder can enhance the decay rate by creating resonant pairs of sites with an energydifference comparable to U . However, sufficiently strong disorder can hinder decay by suppressing diffusion to these resonantsites. numerical simulations of the quantum dynamics, whichcapture many features of our observed results and suggestthe creation of resonances in the lattice by the disorder asthe physical origin of this disorder-enhanced regime. Wefurther supplement this picture by developing a simplephenomenological model that incorporates both disorder-enhanced and disorder-suppressed mechanisms for dou-blon decay. II. MEASURING AND SIMULATINGRELAXATION
We employ an optical lattice experimental platform(Fig. 1) described in previous work [15]. Two spin states(denoted |↑(cid:105) and |↓(cid:105) ) of fermionic K atoms are trappedin a cubic lattice superimposed with 532 nm opticalspeckle disorder. This system realizes the DFHM withconfinement (Fig. 1): H = (cid:88) (cid:104) i,j (cid:105) ,σ (cid:16) − t ij ˆ c † jσ c iσ + h.c. (cid:17) + (cid:88) i U i n i ↓ n i ↑ + (cid:88) i,σ (cid:18) (cid:15) i + 12 mω r i (cid:19) n iσ . (1)Here, σ indexes the two spin states, t ij is the tunnel-ing energy between sites i and j (restricted to nearest- neighbors, indicated by (cid:104) i, j (cid:105) ), U i is the on-site interac-tion energy, (cid:15) i is the local disordered energy offset, ω isthe harmonic confinement, m is the atomic mass, and r i is the distance of site i from the trap center. Thesingle-particle bandwidth is 12 t in the absence of disor-der. The applied speckle potential creates disorder in the t , U , and (cid:15) terms, leading to distributions of these Hub-bard parameters with widths that depend on the opticalpower [22, 23]. We characterize the disorder strength by∆, which is approximately equal to the standard devi-ation of the (cid:15) distribution (Fig. 1(ii)). The influenceof spatial correlations in the disorder potential is weak,since the correlation length of the speckle field is smallerthan two lattice spacings along every lattice direction.Because of the harmonic confinement, all measurementsare averaged over a density profile that varies from anestimated occupancy (cid:104) n (cid:105) = 0 . |↓(cid:105) atom in each molecule to an ancillary spinstate using an rf sweep (see Appendix A 1). Alterna-tively, we can selectively transfer and image only atomsfrom singles, which allows us to separate doublon decayfrom overall number loss.In all regimes, the doublon population decreases fol-lowing the quench with a rate sensitive to the disorderstrength. Typical results at different disorder strengthsare shown in Fig. 2. To quantify the decay, we fit thedata to a model that describes exponential doublon decaywith a time constant τ and that includes overall parti-cle number loss (see Appendix A 1). While we find thatthis fit provides a reasonable characterization of the de-cay timescale, the functional form of the relaxation is un-known beyond the clean limit. We therefore turn towardsnumerics to provide an interpretation of the timescalewith disorder present.The large scale and dimensionality of our system, aswell as the far-from-equilibrium nature of the dynam-ics, preclude exact numerical studies and commonly em-ployed approximate techniques, such as DMRG or diag-onalization methods. Furthermore, the fermionic signproblem forbids use of quantum Monte Carlo meth-ods. We therefore develop a numerical method—ageneralized discrete truncated Wigner approximation(GDTWA) [24–26]—to simulate the relaxation process.The GDTWA approach invokes a factorization of thedensity matrix ˆ ρ of the many-body system over individ-ual lattice sites i , ˆ ρ = (cid:78) i ˆ ρ i . Furthermore, random sam-pling of initial conditions from a discrete semiclassicalphase space accounts for quantum noise.More specifically, we can fully describe the reducedstate ˆ ρ i at lattice site i (with a four-dimensional localHilbert space spanned by basis states {|↑↓(cid:105) , |↑(cid:105) , |↓(cid:105) , | (cid:105)} )by a 16-dimensional vector (cid:126)λ i , so that ˆ ρ i = (cid:80) α λ αi ˆΛ αi fora complete basis of local observables { ˆΛ αi } . Insertingthe product state ansatz into the von Neumann equa-tion ∂ t ˆ ρ = − i (cid:126) [ ˆ H, ˆ ρ ] results in a set nonlinear differentialequations describing the evolution of each (cid:126)λ i , and whichmay be numerically integrated to obtain the mean fielddynamics. While this ansatz retains full information re-garding the strong on-site Hubbard correlations, such asolution describes a product state at all instances in timeand therefore neglects the important cross-site quantumcorrelations responsible for doublon decay. In fact, foran initial product state of definite particle number andspin, the mean field solution results in no dynamics forthis system.To capture the buildup of quantum correlations be-tween sites, we instead describe the initial value of each (cid:126)λ i as a probability distribution over a discrete phasespace that samples over all allowed values for each ob-servable in our local basis (see Appendix B). Averagingindependently evolved sets of randomly sampled initial conditions from this distribution yields a highly non-trivial solution for the dynamics (owing to the nonlin-ear nature of the dynamical equations) and describes theevolution to a correlated state exhibiting entanglement[24, 26]. This approximation, while only rigorously validat short times, has demonstrated the ability to provideaccurate results in generic spin models at longer timesand properly capture quantum thermalization of local ob-servables [24, 27, 28]. Here, for the first time, we adaptthe GDTWA to model DFHM dynamics and we providebenchmarks for its applicability in Appendix B.The GDTWA results, shown as traces in Fig. 2, ap-proximately capture the observed changes in doublonpopulation decay with applied disorder. The GDTWAdynamics also allow us to determine, within the assump-tions of the technique, the extent to which τ accuratelycharacterizes the relaxation timescale. As illustrated inFig. 2, we observe qualitative differences between themeasurements and simulations, which tend to exhibit afaster initial decay (over timescales (cid:126) /t ) before transi-tioning to a slower decay on longer timescales that is wellfit by an exponential decay. Nonetheless, by constrain-ing the initial doublon fraction in the fitted exponentialdecay to the measured value, we are able to obtain a rea-sonable exponential fit of the entire dynamical curve thatis used to extract a consistent relaxation timescale. Theresulting effective timescale incorporates both the initialnon-exponential features and the subsequent exponentialdecay (see Appendix B, Fig. 9).To account for the spatially varying trap density andvariations in the initial preparation, which affects theequilibrium doublon density, we choose initial singlesdensities for the simulations that yield best fits to theexperimental data. These best-fit singles densities in thesimulations are generally consistent with the experimen-tally measured single density at the trap center, exceptat small disorders where GDTWA does not fully capturethe decay processes associated with the clean system andthus results in a much larger equilibrium doublon densityfor all initial singles densities, as evidenced in Fig. 2(i)(see also Appendix B). However, we have also confirmedthat the relaxation timescale is not strongly influenced bythe average density nor the spatial profile of the gas in theregime of interest (see Appendix B, Fig. 8). These ob-servations support characterizing the doublon relaxationdynamics by a single parameter τ : τ primarily dependson the interplay of disorder with interactions and tunnel-ing and is insensitive to parameters that may fluctuateor 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 andnumerics. Strikingly, we find that applying disorder firstcauses the relaxation time τ to rapidly decrease. For theclean system, τ is substantially longer than the single- FIG. 2. Examples of experimental and numerical doublon relaxation data. Sample traces of the doublon population vs. timefor different disorder strengths are shown for experimental measurements (points); the error bars give the statistical standarderror of the mean for averaging over multiple measurements. The fits used to extract τ are shown using solid lines. Thefinite offset value at long times reflects the equilibrium doublon population; see Appendix A 1 and Fig. 5. The correspondingGDTWA simulation (see Appendix B) is displayed as a shaded region that shows the standard error of the mean from disorderand trajectory averaging. The horizontal axis for panel iv (the strongest disorder) is compressed to bring the slow decay intoview, and all data are taken at U/ t = 1 . E R lattice depth in the experiment). particle tunneling time (cid:126) /t , which is consistent withprevious studies that identified doublons as repulsivelybound pairs [18, 19, 29–31]. However, disorder causes τ to decrease to a minimum value comparable to thetunneling time (cid:126) /t at a disorder value near ∆ ∼ U/ U , τ eventually increases,growing by over two orders of magnitude at the strongestdisorder we can apply in the experiment. The separationinto two dynamical regimes (distinguished by the slope of τ with ∆) combined with the crossover at ∆ ∼ U suggest that the dynamics are controlled by competingmechanisms arising from interactions and disorder.We can understand these dynamical regimes usinga minimal model of diffusing doublons in a disorderedenvironment. In this model, the interchange betweendoublon–hole pairs and pairs of opposite-spin singles iscontrolled by a set of reaction–diffusion (R-D) equations.The R-D model (see Appendix C) augments the classi-cal continuum diffusion equation for each particle specieswith a source term that converts a doublon–hole to asingle–single combination, only when the local parame-ters allow this process to be resonant. In particular, ourmodel requires the local energy difference arising fromthe speckle disorder to lie within a window of width FIG. 3. Dynamical regimes of doublon relaxation. The doublon lifetime τ shows a strong non-monotonic dependence ondisorder, first decreasing (green region, labeled 1) and then increasing (blue region, labeled 2) with larger disorder strength.The measurements (points) are compared to numerical GDTWA simulations (crosses). The error bars for the experimentaldata provide the standard error of the fit used to determine τ , while the error bars for the GDTWA numerical data are smallerthan the symbols. The solid line represents a fit of the measurements to the analytical reaction–diffusion model with a diffusionconstant that decreases monotonically with disorder. The lifetime for atom number loss in the experiment τ loss (dash–dottedline) is independent of disorder and well separated from the doublon dynamics timescale. The labels i.–iv. indicate the valuesof ∆ /U corresponding to panels i.–iv. of Fig. 2. ∼ t around U (see Fig. 1). The doublon diffusion co-efficients are taken to decrease monotonically with disor-der strength, as increasingly large local energy differenceswill inhibit doublon transport throughout the lattice.The R-D model gives the decay rate 1 /τ as a productof the effective diffusion rate D eff and the probability P reaction for a conversion between a doublon–hole pairand two singles ( ↑ and ↓ ):1 /τ = P reaction × D eff . (2)The probability of a reaction per site encountered (de-rived in Appendix C) is modeled as P reaction = p + exp ( − U/ ∆) sinh (cid:16) √ zt/ ∆ (cid:17) , (3)where p is a small (i.e., of order 0.01) parameter cor-responding to reactions in the clean limit, and z is thelattice coordination number. The effective diffusion co-efficient D eff gives the rate (linear in time) at which newsites 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 asymptot-ically localized phase predicted at large disorder (wheredoublon diffusion vanishes), D eff should be supplementedby an additional small velocity that allows sampling ofsites within the finite localization length.The exact dependence of diffusion D eff on disorder∆ /U is not qualitatively important as long as D eff de-creases monotonically to a very small or vanishing value.In Fig. 3, we choose a diffusion coefficient that decreasesexponentially with disorder. For this functional form,the best-fit value for D eff is of order t/ (cid:126) at low disor-der (4 ± t/ (cid:126) ) and is reduced by a factor of over 100(to 0.02 ± t/ (cid:126) ) at the highest ∆. We find that analgebraically decaying diffusion coefficient yields quan-titatively similar decay rates at all measured disorderstrengths. The combined expression for the decay rate1 /τ within this reaction–diffusion model explains the tworelaxation regimes as follows.In regime 1 (∆ (cid:46) U ), increasing disorder leads to de-creasing decay times. This regime is controlled by thereaction probability P reaction . In the clean limit, the re-action rate is greatly suppressed by the strong interactionenergy U , and doublon decay occurs slowly. However, asthe disorder strength increases, the tail of the local energydistribution allows for a finite probability of an adjacentsite energy difference of order U . By compensating theinteraction energy U through this mechanism, the disor-der produces a quantum resonance that allows the decayreaction to proceed with high probability, thereby result-ing in fast doublon relaxation.In regime 2 (∆ (cid:38) U ), increasing disorder leads to in-creasing decay times. This regime is controlled by the dif-fusion constant D eff , which becomes heavily suppressedby localization effects produced by the large disorder.To a lesser extent, the gradual suppression of the reac-tion probability due to fewer resonances also contributesto the behavior in this regime. Recent theoretical worksuggests that, aside from rare region fluctuations, diffu-sion should vanish at sufficiently large disorder and signalan asymptotic many-body localized phase [10]. The diffu-sion in the 3D interacting system is an unknown functionof the disorder ratio to the bandwidth ∆ / t and inter-actions ∆ /U , but dimensional considerations suggest areduction when the ratios become of order unity, consis-tent with the observations Fig. 3. While the experimen-tal data cannot distinguish between D eff = 0 and D eff saturating to a finite small value, we can robustly con-clude that diffusion must be highly inhibited in regime2. IV. CONCLUSION
Quantum many-body systems involving interactionsand disorder can exhibit emergent behavior that, espe-cially in two and three dimensions, is challenging to pre-dict from first principles. By experimentally construct-ing a DFHM quantum simulator in a 3D optical lattice,we are able to study the dynamical out-of-equilibriumbehavior of doublons across vastly different scales of dis-order and interactions. Remarkably, we find that con-structing a simple model for this extraordinarily compli-cated system is possible using reaction–diffusion equa-tions that are analytically solvable. We devise a numeri-cal approach that captures the same physical effects andshows quantitative agreement with the experiment.Intriguing open questions remain to be explored. Inthe intermediate disorder regime (∆ ∼ U ), the observedfast doublon relaxation may be related to the behaviorof “bad metals” [33, 34], which can be characterized bya lack of conserved excitations [35]. Because the rapiddoublon decay can be understood as a consequence ofdisorder disrupting the gap in the single-particle spec-trum, spectroscopic measurements in this regime couldidentify a disorder-created pseudogap in the density ofstates [20, 21]. As a practical tool, the fast disorder-mediated thermalization that we observe also suggests anapproach to avoid challenges in the adiabatic preparationof strongly correlated atomic gases [36]. In the strong dis- order regime, our measurements imply a suppression ofparticle transport as disorder is increased. However, weare unable to distinguish whether this behavior is a man-ifestation of slow diffusion or a signature of asymptoticmany-body localization. Further experiments could clar-ify this difference and extend the results of Ref. [15] tothe U/ t > ACKNOWLEDGEMENTS
We acknowledge helpful discussions with GretchenCampbell and Jason Iaconis. This work is supportedby the AFOSR grant FA9550-18-1-0319; the ARO singleinvestigator award W911NF-19-1-0210; the NSF PHY-1820885 and NSF JILA-PFC PHY-1734006 grants; andNIST. It is also supported in part by the AFOSR undergrant number FA9550-17-1-0183 (R.M.N.).
Correspondence and requests for materialsshould be addressed to B. DeMarco (email: [email protected]).
Appendix A: Experimental details1. Parameters t loss (s-1) s ( E R ) D ( E R ) FIG. 4. Experimental atom number loss rate. The measuredloss rate 1 /τ loss is shown for various lattice depths s and dis-order strengths. Both the lattice depth and disorder are mea-sured in units of the recoil energy of the lattice ( E R ). Theloss rate scales linearly with lattice depth and is independentof disorder strength. This behavior is consistent with off-resonant scattering of lattice light as the dominant loss mech-anism. The slope of a linear fit (solid line, constrained to haveno intercept) is 0 . ± . · E R ) − . We use this fit todetermine the value of τ loss used in the fitting model (Eqs.A1-A4) for data sets in which only the doublon population ismeasured. Nd,eq/Ntot D / U FIG. 5. Disorder dependence of steady-state doublon frac-tion ( N d,eq /N tot ), corresponding to the data in Fig. 3. Thesteady-state doublon population N d,eq is extracted from thefits to the doublon data, and corrected for finite imaging effi-ciency, while the total number in one spin N tot is measured inseparate experimental runs. Error bars are the standard errorof the fit. The dashed line is an atomic limit calculation forequilibrium at the initial temperature and number of atoms. We use a balanced mixture of the | F = 9 / , m F = − / (cid:105) and | / , − / (cid:105) K spinstates, which are denoted as |↑(cid:105) and |↓(cid:105) . Typical startingconditions are 100,000 atoms at
T /T F = 0 . T = 100nK). The atoms are initially confined in a 1064 nmoptical trap with the interactions tuned to be attractive(scattering length a s = − . a , where a is the Bohrradius) using the Feshbach resonance near 202.1 G [37].We load the atoms into a 3D cubic optical lattice withlattice depth s = 12 E R , where E R = h / (2 mλ ) is therecoil energy ( m is the atomic mass and λ = 782 . µ s, after which a s = 61 . a and U/ t = 1 . . E R . Thedisorder is generated by a 532 nm speckle laser beam cre-ated using a holographic diffuser [15, 38]. The strengthof the disorder ∆ is tuned by controlling the 532 nmoptical power and is equal to the standard deviation ofthe disorder potential. The disorder ramp rate is cho-sen to be adiabatic relative to the lattice band gap, butfast compared to the doublon decay rate in the absenceof disorder. The rapid turn-on of the disorder may in-duce disorder-dependent heating of the atoms. Using anindirect thermometry method, in which we measure thesize of the gas after adiabatic turnoff of the lattice anddisorder potential, we infer that this effect changes thetemperature by less than 20%.After a variable hold time t hold , the lattice potential depth is quickly ramped to 30 E R to freeze the densitydistribution, and the disorder is removed by turning offthe 532 nm light. The doublon population is measuredby ramping the magnetic field across the Feshbach res-onance to associate doublons into molecules [39], selec-tively transferring the |↓(cid:105) atom in each molecule to theancillary | / , − / (cid:105) state with an rf pulse, and imagingonly the atoms in the ancillary state. Before imaging,we remove all atoms in the |↓(cid:105) and |↑(cid:105) spin states us-ing a combination of rf sweeps and resonant light pulses.Then, the atoms in the ancillary state are transferred tothe | ↑(cid:105) state to be imaged on a closed transition. Wecan also image only the |↓(cid:105) atom singles instead of dou-blons by changing the rf pulse to be resonant with thespin transition for free atoms instead of molecules.
2. Fitting model
We extract τ using a fit to a rate model:˙ N d = ˙ N d,neq + ˙ N d,eq , (A1)˙ N d,neq = − τ N d,neq − τ loss N d,neq , (A2)˙ N d,eq = − τ loss N d,eq , (A3)˙ N s ↓ = − τ loss N s ↓ + 1 τ N d,neq . (A4)This model describes a non-equilibrium population ofdoublons N d,neq that dissociate to create N s ↓ singles |↓(cid:105) (along with undetected |↑(cid:105) atoms) at a rate 1 /τ , a steady-state population of doublons N d,eq , and overall numberloss at rate 1 /τ loss . For the longest doublon lifetimessampled in this work, we perform a simultaneous fit tothe changes in the doublon and singles populations toextract τ , while for the shorter lifetimes we find that thesingles provide no additional constraint. Measurementsof τ loss , which primarily depends on the lattice depth,are shown in Fig. 4.The fit values for N d,eq are shown in Fig. 5. Unlike themeasurements of τ presented in the main text, these areonly physically meaningful given knowledge of the imag-ing efficiency. Experimentally, we find that this efficiencyis lower (by approximately a factor of 2) for measurementof doublons than for singles, which we attribute primar-ily to inelastic decay of Feshbach molecules during theimaging procedure. Therefore, to arrive at an equilib-rium doublon fraction, we scale the measured number bythe initial total atom number (in one spin state) and cor-rect for the unequal imaging efficiencies of doublons andsingles by comparing their respective population changesduring decay. The resulting steady-state doublon frac-tion is comparable to the result from an atomic limitcalculation at equilibrium (dashed line) [40]. This cal-culation neglects the entropy generated by the quench,which is expected to be substantially smaller than theentropy present initially from finite temperature.
3. Dependence of τ on U/ t By changing the lattice depth, we measure the depen-dence of the relaxation on interaction strength, allowingus to map out an experimental relaxation phase diagram.We vary s from 12–15 E R , deep into the doped Mott insu-lator regime at equilibrium [41]. For lower lattice depths,the effect of the quench becomes too small to accuratelyextract τ . The resulting data, along with the correspond-ing numerical GDTWA calculations, are shown in Fig. 6(we omit the reaction-diffusion model, as it provides amodel and fit, rather than an independent prediction).For all values of U we observe the initial disorder-driven decrease in relaxation time, both in experimentand numerics. In both experiment and numerics, we alsosee a strong dependence on disorder at low U , while athigh U the differences become washed out. However, dif-ferences between experiment and numerics are also ev-ident. For high U the regime of increasing relaxationtime appears to recede beyond accessible values of ∆ inexperiment.For the strongest interactions ( U/ t = 2 . , . U/ t ) are one and nine seconds, respectively. However,the numerics also do not capture this limit, in which thedecay mechanism is via higher-order processes [19]. Moregenerally, the experimental data may exhibit many-bodyeffects due to the finite filling and also finite temperatureeffects, which are interesting subjects for future work. Appendix B: Generalized Discrete TruncatedWigner Approximation
To numerically simulate doublon dynamics, we adaptthe generalized discrete truncated Wigner approximation(GDTWA) to our system [24]. We begin by choosing atensor product structure for the Hilbert space; for theDFHM, we let each subspace correspond to the four-dimensional Hilbert space for each lattice site with basis {|↑↓(cid:105) , |↑(cid:105) , |↓(cid:105) , | (cid:105)} . We then invoke a factorization of thedensity matrix as a function of time, ˆ ρ ( t ), along this ten-sor product structure, resulting inˆ ρ ( t ) = (cid:79) j ˆ ρ j ( t ) . (B1)This ansatz allows us to retain information regardingthe strong on-site correlations between spin species re-sponsible for doublon formation, while treating tunnel-ing induced correlations between sites in a semiclassicalmanner. However, this method also implicitly invokes areplacement of the anti-commuting fermionic operatorswith commuting hardcore boson operators. Nonetheless,we expect this to be a good approximation, given that the doublon population relaxation dynamics are expected tobe highly insensitive to quantum statistics in the diluteregime explored in the experiment.In principle, it is possible to increase the size of eachHilbert space in our product state factorization ansatz(e.g., by including multiple lattice sites within each den-sity matrix factor) and thus presumably increase the ac-curacy of the simulated dynamics. However, the com-putational cost associated with randomly sampling overdisorder distributions precludes the use of this strategyin this work. Already, a cluster of only two sites requiresa description via 256 phase space variables per cluster.We have performed limited simulations for select param-eters using this clustering method, and generally foundthat the resulting decay times are robust, though theequilibrium doublon values are altered slightly.The factorization approximation directly results in aset of mean-field-like equations for variables λ αj ( t ) ≡(cid:104) ˆΛ αj ( t ) (cid:105) , where the Hermitian operators ˆΛ αj correspondto the set of generalized Gell-Mann matrices (GGM) forSU(4) (along with the identity matrix) with index α foreach lattice site j , which form a complete orthonormalbasis under a Hilbert-Schmidt norm for the observableson j . The resulting mean-field equations are given by: dλ αj d t = i (cid:88) β M αβj λ βj ( t ) + (cid:88) β,β (cid:48) ,j (cid:48) C αββ (cid:48) jj (cid:48) λ βj ( t ) λ β (cid:48) j (cid:48) ( t ) (B2)where, for a generic Hamiltonian written as ˆ H = (cid:80) j ˆ H j + (cid:80) jj (cid:48) ˆ H jj (cid:48) , we have M αβj = Tr (cid:104) ˆΛ βj [ ˆ H j , ˆΛ αj ] (cid:105) and C αββ (cid:48) jj (cid:48) =Tr (cid:104) ˆΛ βj ˆΛ β (cid:48) j (cid:48) [ ˆ H jj (cid:48) , ˆΛ αj ] (cid:105) . We use the equations arising fromthis approximation (which lead to trivial mean-field dy-namics) to dynamically evolve multiple trajectories ran-domly sampled from a phase space [24], with a samplingdistribution determined by the initial product state of thesystem. Specifically, if each observable basis element haseigendecomposition ˆΛ αj = (cid:80) a λ α, [ a ] j ˆ P α, [ a ] j for eigenvalue λ α, [ a ] and corresponding eigenspace projector ˆ P α, [ a ] j , thenfor each trajectory we set λ αj (0) = λ α, [ a ] j with probability p αj ( a ) = Tr (cid:104) ˆ ρ j (0) ˆ P α, [ a ] j (cid:105) .In our simulations, we neglect the relatively small dis-order in the tunnelling and interaction parameters. Wealso neglect the presence of the harmonic confining trap:assuming diffusive motion, the doublons in the trap cen-ter will not travel an appreciable distance towards reso-nant regions in the trap edges before decaying. We ex-pect both of these features to result in small shifts to theeffective disorder strength and diffusion coefficients, butnot to affect the main features of our results.As we are interested in only the total doublon popula-tion, we let each sampled trajectory correspond to a dif-ferent quenched disorder realization and a different ran-dom initial product state configuration respecting the lat-tice filling (see comparisons with exact results in Fig. 7). U / 1 2 t
G D T W Al o g ( t t / (cid:1) ) D /12t U / 1 2 t
FIG. 6. Dependence of measured τ on U/ t . Left: experimental data. The lattice depth is tuned from 12–15 E R to vary U/ t by a factor of approximately 2.5. The crosses represent experimental points, which are interpolated to create the colormap.Right: GDTWA numerical results. (a) (b)(d) (e) (f)(c) FIG. 7. GDTWA benchmarking in 1D. A comparison of the total doublon number dynamics obtained from exact diagonalizationof the DFHM(a,d) and GDTWA (b,e) is shown for a periodic 1D chain of eight lattice sites. These simulations were carried outwith U/ t = 2 .
5, and the results were averaged over 100 disorder realizations for each disorder strength. The top row (a-c)corresponds to an initial loadout of two doublons and no singles, while the bottom row (d-f) corresponds to an initial loadoutof two doublons and one single of each flavor. Panels (c) and (f) give the doublon population time-averaged over the last halfof the time interval shown. For system sizes accessible by exact simulation, the initial doublon population quickly jumps to thesteady state, without any further relaxation features observable on longer timescales.
Each dynamical curve in Fig. 2 and the dynamical curves corresponding to the decay times in Fig. 3 results from0 n s = 0 . 1 0 n s = 0 . 1 9 n s = 0 . 2 8 n s = 0 . 3 8 t t/ (cid:1) D / U FIG. 8. GDTWA variation of decay times with singles den-sity. A comparison of GDTWA decay times at selected valuesof ∆ /U for a 6 × × U/ t = 1 .
8, andthe initial doublon density remains fixed at 11%. While thesteady-state population depends strongly on the initial sin-gles density, the robustness of the decay time to this parame-ter allows comparison with experiment, despite the spatiallyvarying density in the trap. averaging 500-1,000 such trajectories in a 4 × × .
11, and thetotal particle density to (cid:104) n (cid:105) = 0 . t t/ (cid:1) D / U FIG. 9. GDTWA relaxation times from biased fits. A com-parison between experimental (black circles) and GDTWAdecay times (blue squares) is shown. The GDTWA decaytimes are extracted utilizing a (disorder-dependent) samplingfunction consisting of delta peaks at the experimental holdtimes, weighted by the appropriate experimental uncertainty.The error bars show the fit standard error, obtained from alinearized curvature analysis. while the second is optimized for comparison to the ex-periment. In the first method (displayed in Fig. 3), wefit the dynamical curves to an exponential decay with anonzero steady-state, and we constrain this function tostart at the exact initial doublon number. This techniqueleads to decay times that incorporate the effect of theserapid early time decay dynamics (see Fig. 2). An uncon-strained fit virtually ignores the early time regime and ismore heavily influenced by the long period of slow, expo-nential relaxation at late times, leading to an underesti-mate of the decay time and a significant underestimate ofthe initial doublon number. In the second method, whichis displayed in Fig. 9, we only use simulation data attimes corresponding to the experimental measurements.This latter method takes into account experimental biasfrom the selected hold times, thus providing a consis-tency check between simulation and experiment. Thistechnique produces generally better quantitative agree-ment in the crossover regime.In contrast to the decay time, the steady-state doublonfraction is strongly dependent on the singles density andgenerally increases as more initial singles are added tothe system. To compare with the dynamical data fromthe experiment in Fig. 2, for each disorder strength weperform simulations for various singles densities (each dif-fering by ≈ .
03) and select the simulation correspondingto the best weighted least-squares fit to the experiment.Keeping the initial doublon density fixed at 0 .
09, we findbest-fit total particle densities (cid:104) n (cid:105) = 0 .
22, 0 .
53, 0 .
53, and0 .
41 for Fig. 2(i), (ii), (iii), and (iv), respectively. Whilethis best-fit singles density is generally representative ofthe initial densities measured in the experiment over awide range of disorders, for small disorders the best-fit1singles density are well below reasonable values for thesingles density in the experiment, to the point where avanishingly small initial singles density cannot accommo-date the measured change in doublon density. This is theresult of various decay processes in the clean system notbeing captured by GDTWA, resulting in a higher equi-librium doublon fraction. However, we still find that theextracted decay timescales are generally consistent withexperiment in this regime.
Appendix C: Reaction-Diffusion Model
In this model, we let each particle species obey a R-D equation, e.g. for the doublon number n d ( r , t ), as afunction of space r and time t , ∂ t n d = D d ∇ n d + S ( r )( n s ↑ n s ↓ − n d n h ) . (C1)The source term S ( r ) has contributions from quantumfluctuations and disorder. Here we focus on the lattercontribution, relevant when disorder is not too small. S ( r ) is modeled by a small constant clean-system reac-tion probability p , together with a Boolean step functionsignifying resonance with local energy difference µ ( r ), S ( r ) = p + Θ ( γ − | µ ( r ) − U | ) , (C2)where γ is the width of the resonance and µ ( r ) is thedifference of two independent random local energies fromthe speckle distribution, ℘ ( µ ) = 12∆ exp ( −| µ | / ∆) . (C3) The resulting S ( r ) distribution yields a reaction ( S ≈
1) with probability P ( { S = 1 } ) = p + (cid:90) U + γU − γ ℘ ( µ ) dµ = p + exp( − U/ ∆) sinh( γ/ ∆) , (C4)which is also the average of S ( r ). We expect γ ≈ √ zt forlattice coordination number z ( z = 6 for a cubic lattice in3D), though the exact form is unimportant for the generalconclusions we wish to draw. We also find that the formof S ( r ) is qualitatively similar to the form obtained byconvolving the speckle distribution with a heavy-tailedLorentzian resonance peak.To proceed, let us work to linear order in n d (con-trolled in the dilute limit) and take t (cid:28) ∆, arriving (with D = D d ) at a linear differential equation with an inho-mogenous random source term, ∂ t n d ( r , t ) = D ∇ n d ( r , t ) − S ( r ) n d ( r , t ) . (C5)The homogeneous solution to the heat equation of apoint particle in d = z/ n ( r , t ) =(4 πD t ) − d/ exp (cid:0) − r / D t (cid:1) . The linear diffusion rangegrows as √ dD t , which is also proportional to the numberof sites sampled when d = 1. In 3D, however, the numberof sites sampled grows linearly in time as N S = αD t with α ≈ .
32. This yields a decay constantΓ = αD (∆) (cid:104) p + exp( − U/ ∆) sinh (cid:16) √ zt/ ∆ (cid:17)(cid:105) . (C6)For a diffusion constant of the form D (∆) = D exp ( − ∆ / ∆ ), we fit the model to the experimen-tal data with unknowns p , D , and ∆ , resulting inthe line shown in Fig. 3. The best fit parameters are p = 0 . ± . (cid:126) D /t = 4 ±
1, and ∆ /U = 0 . ± . [1] P. W. Anderson, Phys. Rev. , 1492 (1958).[2] P. A. Lee, Rev. Mod. Phys. , 287 (1985).[3] N. F. Mott, Proc. Phys. Soc. Sect. A , 416 (1949).[4] D. Belitz and T. R. Kirkpatrick, Rev. Mod. Phys. ,261 (1994).[5] K. Byczuk, W. Hofstetter, and D. Vollhardt, Int. J. Mod.Phys. B , 1727 (2010).[6] L. Sanchez-Palencia and M. Lewenstein, Nature Physics , 87 (2010).[7] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Phys.Rev. Lett. , 206603 (2005).[8] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Ann.Phys. , 1126 (2006).[9] R. Nandkishore and D. A. Huse, Annu. Rev. Condens.Matter Phys. , 15 (2015).[10] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev.Mod. Phys. , 021001 (2019).[11] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Luschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, andI. Bloch, Science , 842 (2015).[12] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess,P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Nat.Phys. , 907 (2016).[13] P. Roushan, C. Neill, J. Tangpanitanon, V. M. Bastidas,A. Megrant, R. Barends, Y. Chen, Z. Chen, B. Chiaro,A. Dunsworth, A. Fowler, B. Foxen, M. Giustina, E. Jef-frey, J. Kelly, E. Lucero, J. Mutus, M. Neeley, C. Quin-tana, D. Sank, A. Vainsencher, J. Wenner, T. White,H. Neven, D. G. Angelakis, and J. Martinis, Science , 1175 (2017).[14] H. P. L¨uschen, P. Bordia, S. Scherg, F. Alet, E. Altman,U. Schneider, and I. Bloch, Phys. Rev. Lett. , 260401(2017).[15] S. S. Kondov, W. R. McGehee, W. Xu, and B. DeMarco,Phys. Rev. Lett. , 083002 (2015).[16] J.-y. Choi, S. Hild, J. Zeiher, P. Schauss, A. Rubio- Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch,and C. Gross, Science , 1547 (2016).[17] P. Bordia, H. L¨uschen, S. Scherg, S. Gopalakrishnan,M. Knap, U. Schneider, and I. Bloch, Phys. Rev. X ,041047 (2017).[18] N. Strohmaier, D. Greif, R. J¨ordens, L. Tarruell,H. Moritz, T. Esslinger, R. Sensarma, D. Pekker, E. Alt-man, and E. Demler, Phys. Rev. Lett. , 080401(2010).[19] R. Sensarma, D. Pekker, E. Altman, E. Demler,N. Strohmaier, D. Greif, R. J¨ordens, L. Tarruell,H. Moritz, and T. Esslinger, Phys. Rev. B , 224302(2010).[20] E. Lahoud, O. N. Meetei, K. B. Chaska, A. Kanigel, andN. Trivedi, Phys. Rev. Lett. , 206402 (2014).[21] Z. Wang, Y. Okada, J. O’Neal, W. Zhou, D. Walkup,C. Dhital, T. Hogan, P. Clancy, Y.-J. Kim, Y. F. Hu,L. H. Santos, S. D. Wilson, N. Trivedi, and V. Madha-van, Proc. Natl. Acad. Sci. , 11198 (2018).[22] S. Q. Zhou and D. M. Ceperley, Phys. Rev. A , 013402(2010).[23] M. White, M. Pasienski, D. McKay, S. Q. Zhou,D. Ceperley, and B. DeMarco, Phys. Rev. Lett. ,055301 (2009).[24] B. Zhu, A. M. Rey, and J. Schachenmayer, New. J. Phys. , 082001 (2019).[25] J. Schachenmayer, A. Pikovski, and A. M. Rey, Phys.Rev. X , 011022 (2015).[26] O. L. Acevedo, A. Safavi-Naini, J. Schachenmayer, M. L.Wall, R. Nandkishore, and A. M. Rey, Phys. Rev. A ,033604 (2017).[27] S. Lepoutre, J. Schachenmayer, L. Gabardos, B. H.Zhu, B. Naylor, E. Mar´echal, O. Gorceix, A. M. Rey,L. Vernac, and B. Laburthe-Tolra, Nature Communica- tions , 1714 (2019).[28] A. Patscheider, B. Zhu, L. Chomaz, D. Petter, S. Baier,A.-M. Rey, F. Ferlaino, and M. J. Mark, Phys. Rev.Research , 023050 (2020).[29] K. Winkler, G. Thalhammer, F. Lang, R. Grimm,J. Hecker Denschlag, A. J. Daley, A. Kantian, H. P.B¨uchler, and P. Zoller, Nature , 853 (2006).[30] D. Hansen, E. Perepelitsky, and B. S. Shastry, Phys.Rev. B , 205134 (2011).[31] J. P. Covey, S. A. Moses, M. Grttner, A. Safavi-Naini,M. T. Miecnikowski, Z. Fu, J. Schachenmayer, P. S. Juli-enne, A. M. Rey, D. S. Jin, and J. Ye, Nat. Comm. ,11279 (2016).[32] B. Deissler, M. Zaccanti, G. Roati, C. D’Errico, M. Fat-tori, M. Modugno, G. Modugno, and M. Inguscio, Na-ture Physics , 354 (2010).[33] W. Xu, W. R. McGehee, W. N. Morong, and B. De-Marco, Nat. Commun. , 1588 (2019).[34] P. T. Brown, D. Mitra, E. Guardado-Sanchez,R. Nourafkan, A. Reymbaut, C.-D. H´ebert, S. Bergeron,A.-M. S. Tremblay, J. Kokalj, D. A. Huse, P. Schauß,and W. S. Bakr, Science , 379 (2019).[35] S. A. Hartnoll, Nature Physics , 54 (2014).[36] C. L. Hung, X. Zhang, N. Gemelke, and C. Chin, Phys-ical Review Letters , 1 (2010).[37] T. Loftus, C. Regal, C. Ticknor, J. L. Bohn, and D. S.Jin, Phys. Rev. Lett. , 173201 (2002).[38] S. S. Kondov, W. R. McGehee, J. J. Zirbel, and B. De-Marco, Science , 66 (2011).[39] R. J¨ordens, N. Strohmaier, K. G¨unter, H. Moritz, andT. Esslinger, Nature , 204 (2008).[40] L. De Leo, J.-S. Bernier, C. Kollath, A. Georges, andV. W. Scarola, Physical Review A , 023606 (2011).[41] D. Semmler, J. Wernsdorfer, U. Bissbort, K. Byczuk,and W. Hofstetter, Phys. Rev. B82