Probing Transport and Slow Relaxation in the Mass-Imbalanced Fermi-Hubbard Model

Constraints in the dynamics of quantum many-body systems can dramatically alter transport properties and relaxation timescales even in the absence of static disorder. Here, we report on the observation of such constrained dynamics arising from the distinct mobility of two species in the one-dimensional mass-imbalanced Fermi-Hubbard model, realized with ultracold ytterbium atoms in a state-dependent optical lattice. By displacing the trap potential and monitoring the subsequent dynamical response of the system, we identify suppressed transport and slow relaxation with a strong dependence on the mass imbalance and interspecies interaction strength, consistent with eventual thermalization for long times. Our observations demonstrate the potential for quantum simulators to provide insights into unconventional relaxation dynamics arising from constraints.


I. INTRODUCTION
The fundamental understanding of thermalization and its failure in isolated quantum many-body systems has seen remarkable progress in recent years, driven both by novel theoretical concepts [1][2][3][4] and emerging experimental platforms for quantum simulations such as ultracold atoms, trapped ions, and superconducting qubits [5][6][7].The simplest approach to the problem is to posit the two distinct classes of ergodic and nonergodic dynamics.In the former case, fast relaxation to a local equilibrium is followed by slow thermalization of globally conserved quantities according to the laws of hydrodynamics [8][9][10][11].By contrast, nonergodic dynamics can arise in systems with a large number of conserved quantities, as exemplified by many-body localization in strongly disordered quantum systems [12][13][14][15] that retain their memory of an initial state for arbitrarily long times [16][17][18].
Recent experimental and theoretical investigations of nonequilibrium dynamics in isolated quantum manybody systems, however, suggest considerable refinements to this classification.For instance, Rydberg quantum simulators have observed surprising oscillatory dynamics in the blockade regime [19,20], and the slow late-time dynamics of fractonic quantum matter with constrained excitations [21] is described by substantially modified hydrodynamic equations [22][23][24].Furthermore, interacting mixtures of heavy and light particles [see Fig. 1(a)] have even been proposed to feature a dynamical type of manybody localization arising from the heavy particles acting as a form of disorder for the light ones [25][26][27][28].Although subsequent studies have suggested that such heavy-light mixtures are ergodic and do thermalize at late times, the relaxation of initial nonequilibrium states is expected to be unconventionally slow [27,[29][30][31][32].This is a direct consequence of the strongly constrained motion in the presence of mass imbalance and interactions.However, since such heavy-light mixtures are particularly challenging to simulate with classical resources, these theoretical studies do not entirely agree on the exact properties of relaxation [30][31][32], which underlines the importance of experimental evidence.Crucially, this phenomenology connects to the general question of how dynamical constraints can introduce slow equilibration and nonergodicity in quantum many-body systems [33][34][35][36][37][38].
In this work, we experimentally study the nonequilibrium dynamics of a heavy-light mixture with an ultracold quantum gas of 171 Yb atoms in an optical lattice.Our experiment realizes the strongly mass-imbalanced Fermi-Hubbard model in one dimension (1D), as illustrated in Fig. 1(a).In addition to the 1 S 0 ground state (denoted |g ), the alkaline-earth-like atom 171 Yb features the metastable 3 P 0 excited state (denoted |e ), often referred to as the clock state [39].We harness a statedependent optical lattice (SDL) [40] to introduce different timescales for the hopping of atoms in |g and |e taking on the role of light and heavy particles, respectively.Our system is harmonically confined, and we gradually displace the trap minimum with the help of an additional optical potential to probe the transport dynamics of the light species [see Fig. 1(b)].Measuring the dynamics of the atomic cloud, we find relaxation at late times as a signature of ergodicity for all finite interaction parameters.However, this relaxation becomes slow compared to conventional timescales for strong interactions and large mass imbalance.Our observations cannot be explained by a mere reduction of the mobility for the heavy species alone, but rather they can be understood as emergent properties from the constrained many-body dynamics.Numerical simulations based on matrix-product states and exact diagonalization support our findings.

II. EXPERIMENT
Our experiment begins with a Fermi gas of 171 Yb atoms in a balanced mixture of the two nuclear spins |m F = −1/2 ≡ |↓ and |m F = +1/2 ≡ |↑ in the ground state |g .In the initial optical dipole trap, we prepare a total of approximately 10 4 atoms at a temperature of T 0.15T F , with T F the Fermi temperature.The weakly interacting spin mixture is loaded from the optical dipole trap into the ground band of a two-axes, approximately 30E m rec deep, state-independent lattice operated at the magic wavelength (λ m = 759.3nm) [41] and a V L = 6.9(3)E rec deep SDL (λ = 671.5 nm) along the third axis.Here, E m rec = h × 2.0 kHz and E rec = h × 2.6 kHz are the recoil energies of the corresponding lattice photons.The atoms are nonuniformly distributed across the array of decoupled 1D systems (tubes) generated by the perpendicular state-independent lattices.We estimate approximately 300 tubes are considerably filled with a mean atom number of N 18 per spin state and standard deviation ∆N 8.In a typical tube, the atoms are spread over a system size of l ≈ 30 lattice sites (rootmean-square width) [42].
Shortly after loading the lattices, we use a 0.17-ms-long clock laser pulse to selectively drive atoms from |g ↑ to |e ↑ .This suddenly introduces a distinct hopping amplitude t H h × 5 Hz t L for the clock state atoms in |e ↑ ≡ |H (heavy).Here, t L h × 105 Hz is the unaltered hopping amplitude of the remaining ground-state atoms in |g ↓ ≡ |L (light).The ratio of the hopping amplitudes and the associated mass imbalance is determined by the different lattice depths V L and V H = 3.06(4)V L experienced by |L and |H atoms in the SDL [42].We ensure that the mixture initially remains noninteracting by ramping the magnetic field to the zero crossing of an orbital Feshbach resonance at approximately 1533 G [42,43] before the excitation pulse.Finally, we lower the SDL depth to adjust the hopping ratio t H /t L and slowly ramp the magnetic field to 1400-1600 G in accordance with the desired interaction strength.We note that the fraction of doublons as a critical property of the initial state only weakly depends on the chosen interaction parameter [42].After the state preparation procedure, we translate the minimum of the trapping potential along the tubes by slowly turning on a displaced, state-independent dipole trap, which initiates the transport dynamics.Following the subsequent evolution, we rapidly ramp up the SDL to V L ≈ 15E rec , which freezes the motion of light and heavy atoms.We then turn off the magnetic field and record the density of the light atoms with in situ absorption imaging and determine the fraction transported to the new trap minimum [see Fig. 1(b)].Our imaging intrinsically integrates along one axis of the system, which averages the measurement over an ensemble of tubes with different atom numbers [42].
For a completely frozen heavy species (t H = 0), the Hamiltonian Ĥ decouples into single-particle terms for the light species.In this limit, known as the Falicov-Kimball model for κ = 0 [44], the heavy particles only contribute as an effective binary disorder potential giving rise to single-particle localization [45].Previous experiments have employed binary mixtures of bosons to probe such a regime, with disorder formed by a fully localized species [46].In contrast, we study a system of two mobile species with distinct hopping amplitudes t L t H > 0 [see Fig. 1(c)], realizing two strongly differing but finite scales in the Hamiltonian.For finite interaction strength |U | > 0, numerical simulations have suggested that the resulting dynamical constraints can lead to anomalous transport and unconventionally slow thermalization with a nontrivial dependence on the parameters of the model [27,[29][30][31][32].At late times, this behavior is expected to be superseded by diffusive transport, as shown in Fig. 1(c).Interestingly, for large interactions and strongly differing masses, the crossover timescale is expected to become extraordinarily large-a genuine many-body effect [42].In the following measurements, we probe the nature of this transient regime at different times with parameters t H /t L > ∼ 0.1 and |U |/t L < ∼ 10.

IV. LOCALIZED SINGLE-PARTICLE STATES
In a first reference measurement, we characterize the single-particle physics originating from the harmonic confinement κ > 0 in our experiment.It gives rise to singleparticle localized eigenstates at the edge of the trap, a phenomenon studied in Refs.[47][48][49][50][51][52].These Wannier-Stark-localized states occur because of the finite gradient ∂ i Ĥ = κ(i − i 0 )n iα between neighboring lattice sites [53].To experimentally probe this effect, we prepare a clean sample of noninteracting light atoms by employing a short resonant "push" pulse that removes all |g ↑ atoms.Subsequently, we measure the response to a linear translation of the trap minimum as shown in Fig. 1(b).Here, the trap minimum is displaced by ∆x/d = i 1 − i 0 = 47(3) within time 90(5)h/t L and we consider variable κ/t L by adjusting the SDL depth in the range 0.69(3)-11.8(5)Erec .
For all finite lattice depths, we observe a separation of the atoms into two clouds, one close to the initial trap minimum i 0 and one at the final trap minimum i 1 , as visible in the insets of Fig. 2(a).We quantify this by determining the number of atoms in the right half of the system N r , which corresponds to i ≥ i 0 + ∆x/(2d), and compare it to the total number of atoms N .When measuring N r before the displacement of the trap minimum, we find a nearly constant value of approximately 10%, which we attribute to our finite imaging resolution and a small extent of the initial cloud into the counting region.In Fig. 2(a), we plot the fraction of transported atoms, N r /N , which is exponentially suppressed for increasing confinement strength.This significant reduction of mass transport in the system results from the properties of the single-particle eigenstates, which are visual- We add a small linear potential to lift the degeneracy of the localized states.
ized in Fig. 2(b).In general, eigenstates with energies E n ∈ [−2t L , 2t L ] are delocalized across the trap with a nonzero probability density at the center.This behavior changes dramatically for energies E n > 2t L corresponding to localized states, in which atoms cannot efficiently follow the displaced trap minimum.Crucially, the number of states with E n ≤ 2t L is reduced as we increase κ/t L and therefore, a larger number of localized states becomes occupied in each tube.This explains the observed suppression of N r /N , which we also reproduce with a theoretical calculation [42] of the experimental protocol [see Fig. 2(a)].

V. INHIBITED TRANSPORT AT EARLY TIMES
Next, we probe how interactions with the heavy species modify the mobility of the light atoms.For this measurement, we set the interaction strength U/t L ∈ [−20, 5], a range accessible via magnetic-field tuning of the orbital Feshbach resonance, and choose the fixed hopping ratio t H /t L 0.1 and confinement κ/t L 1.7 × 10 −2 by ramping the SDL to 4.7(2)E rec .Here, the confinement κ/t H 10κ/t L = 1.7 × 10 −1 suggests that most heavy atoms occupy single-particle localized states ac- cording to our previous measurement in the noninteracting regime [see Fig. 2(a)].However, in the presence of interactions, several of these states delocalize in analogy to the suggested critical gradient required for localization of the interacting Wannier-Stark ladder [48,[54][55][56][57].
We numerically verify this for our mass-imbalanced case and finite but weak interactions.Here, we find a central trap region of approximately 30 lattice sites, in which the heavy atoms already become mobile during state preparation [42].Thus, a significant fraction of the system will be responsive to the displacement of the trap potential.
In analogy to our first measurement, we again determine the fraction of transported atoms after displacing the trap minimum by 47(3)d within time 92(5)h/t L .As we increase the interaction strength |U |/t L , we find a substantial reduction of the transported fraction N r /N up to a factor of 2 compared to U = 0 (see Fig. 3).The behavior for attractive (U < 0) and repulsive (U > 0) interactions is almost identical, resulting from a similar initial state [42] and the dynamical symmetry of our model in the limit κ = 0 [58].We find that the most significant change in the fraction of transported atoms occurs for interaction energies |U |/t L < 2, whereas the signal saturates and remains nearly constant at a low value for |U |/t L > 4. To support our experimental data, we perform matrix-product-state simulations for the Fermi-Hubbard model described in Eq. ( 1) and a simplified version of the experimental protocol [42].As shown in the inset of Fig. 3, the extracted N r /N curve agrees qualitatively with the experiment.Quantitative disagreement is expected due to the significantly lower atom number in the simulation, N sim = 5 N .In our measurement, the apparent suppression of transport results from the constraints in the dynamics of the light atoms arising from interactions with the heavy species, as explored in Fig. 1(c).We emphasize that this measurement primarily probes the early-time dynamics on the timescale h/t H since the system is effectively frozen and fully separated towards the end of the gradual displacement over the relatively large distance ∆x/d 47, which exceeds the typical system size l ≈ 30 (see bottom panels of Fig. 3).

VI. SLOW RELAXATION AT LATE TIMES
We now turn to the dynamical response of the heavylight mixture and investigate whether light atoms, remaining close to their initial position during the trap translation, relax towards the final trap minimum at much later times.Such dynamics lead to signatures in the in situ density, which we systematically record for a variable hold time.Here, we reduce the displacement to ∆x/d 20 < l while keeping the speed approximately 0.5d(t L /h) unchanged to ensure that the system remains connected and a significant fraction of the heavy atoms are mobile over the traversed distance.The relaxation dynamics can be captured by quantifying the change of the in situ density distribution between the initial state after the translation of the trap minimum at a hold time τ = 0 and later times 0 < τ < ∼ 400h/t L .To this end, we introduce the density deviation observable, which is obtained from the density of the light atoms integrated perpendicularly to the transport direction x and normalized such that dx n(x, τ ) = 1.For a slow relaxation of the density, this observable grows monotonously with increasing τ , whereas a constant δn(τ ) indicates a stationary state of the system.We probe the influence of the interaction strength on the dynamics by measuring the density deviation for t H /t L 0.1 and U/t L 0, −2, −10, as shown in Fig. 4(a).First, we focus on the noninteracting time trace with U/t L 0, which shows transient dynamics of δn(τ ) at early times for both hopping ratios.These large-amplitude transients are quickly damped due to contributions from different eigenstates.For later times τ > ∼ 100h/t L , the density deviation reaches an almost stationary signal, which we also expect from our numerical calculation.
For finite interaction strengths, we find a contrasting behavior with a monotonously increasing density deviation δn(τ ) for τ > 0. In this regime, the integrated atomic densities n(x, τ ) show that the light atoms drift towards the final trap minimum at late times [see panels right of Fig. 4(a)].To extract the timescale of the associated relaxation dynamics, we fit the experimental data to the function δn(τ ) = a(1 − e −τ /τ0 ) + c with the exponential decay constant τ 0 and parameters a, c.For U/t L −2, we find good agreement between the fitted function and the experimental data with τ 0 = 85(17)h/t L .In contrast, for the larger interaction strength U/t L −10, the fit parameters a and τ 0 take unreasonably large values yielding a linear curve.Until τ ≈ 100h/t L , the measured density deviation grows far less compared to the case with U/t L −2.To compare the dynamics for the two interaction strengths, we fit the measured density deviation to a linear curve for early times τ ≤ 101h/t L [42].Here, we find the linear slope 2(2) × 10 −6 t L /h for U/t L −10 compared to 13(2) × 10 −6 t L /h for U/t L −2.The significantly smaller linear slope illustrates the emergence of unconventionally slow dynamics for strong interactions.
To demonstrate that the slow relaxation stems from dynamical constraints induced by mass imbalance, we also measure the density deviation for the significantly larger t H /t L = 0.34 (3).For this value of the hopping ratio, the transient regime of unconventional relaxation is expected to be significantly shorter than for the larger ratio studied previously, as illustrated in Fig. 1(c).By changing the wavelength of the SDL to λ = 690.1 nm such that V H = 1.97(5)VL , we realize parameters closely matching our first measurement apart from the hopping ratio t H /t L [42].Because of this, the confinement strength κ/t L 1.7 × 10 −2 also remains unchanged, ensuring directly comparable systems, in particular, regarding Wannier-Stark-localized states of light atoms.As shown in Fig. 4(b), for U/t L 0 and −2, the time-dependent behavior remains largely similar to the previous measurement.However, the time trace for U/t L −10 now only differs negligibly from the one for weaker interactions.From the numerical fit of δn(τ ), we extract similar decay constants τ 0 = 120(22)h/t L and 138(26)h/t L for the two interaction strengths U/t L = −2 and −10, agreeing with our qualitative observation.The experimental data for t H /t L 0.3 confirm a strong dependence of the relaxation dynamics on the interaction strength and mass imbalance between the heavy and light particles.Hence, the unconventionally slow trans-port emerging in our system can be ascribed to genuine many-body effects.
This observation is further supported by considering the limit of strong interactions and mass imbalances much larger than probed experimentally.In this regime, small-scale numerical calculations predict a pronounced change in the relaxation dynamics of a density modulation, setting in at the timescale τ ∼ hU/t 2 H (valid for U t L ) [31].This timescale arises from the collective motion of a bound state comprised of two heavy particles on neighboring sites and a light one delocalized over these sites [27].However, we find numerical evidence that a noticeable change of dynamics at τ vanishes with increasing system size for moderate hopping ratios t H /t L ≥ 0.01 [42].These findings are consistent with the experimentally observed relaxation, which does not exhibit prominent features at τ .The smooth relaxation dynamics in larger systems could potentially originate from an increased relevance of the many-body medium and a hierarchy of relaxation timescales.Such a hierarchy would emerge from N -body bound states consisting of (N − 1) heavy particles and one light particle, only becoming mobile at the timescale hU t We briefly comment on how experimental imperfections could also drive the relaxation of δn(τ ) in addition to the closed system dynamics of Eq. ( 1).Most notably, we expect the mobility of the light species to increase from the small loss of heavy atoms, which becomes appreciable for τ 100h/t L and reaches up to approximately 20% at τ 400h/t L [42].However, these losses should be mostly independent of U/t L since the dominant dissipation process is the off-resonant scattering of SDL photons.By contrast, our data show a strong interaction dependence.Thus, we conclude that the distinct features of the time traces cannot be explained by dissipation.

VII. DISCUSSION AND OUTLOOK
We have characterized the density dynamics of the mass-imbalanced Fermi-Hubbard model after a gradual change of the external trapping potential, which can be summarized as follows.First, our data show relaxation of the density for all finite interaction parameters and hopping ratios at late times.Second, for the strongest interaction U/t L −10 and smallest hopping ratio t H /t L 0.1 probed in our experiment, the system relaxes much slower as a result of the dynamical constraints induced by the heavy species.Comparing this observation to the much faster relaxation for U/t L −2 shows that the slow emergent timescale must involve many-body effects.Thus, our experimental platform gives access to a particularly interesting class of unconventionally slow dynamics that exist in the absence of static disorder and can be controlled solely by interactions.
In future experiments, variable-wavelength density modulations could provide an ideal setting to directly study anomalous transport in this system.Moreover, it would be particularly interesting to investigate the effect of integrability on transport in the one-dimensional massbalanced Fermi-Hubbard model and explore the weak breaking of integrability by only slightly detuning the hopping amplitudes of the two species.Our work also motivates the experimental study of transport in other strongly constrained states of quantum matter, which can be found in Rydberg-blockaded systems, strongly tilted lattices, and realizations of lattice gauge theories [19,24,37,38,59,60].

Supplemental Material
Probing transport and slow relaxation in the mass-imbalanced Fermi-Hubbard model

S.I. EXPERIMENTAL TECHNIQUES A. State preparation
The details on the preparation of degenerate gases of 171 Yb are described in Ref. [1] and we only briefly outline the procedure here.Our experiment starts by loading approximately 1.4 × 10 6 174 Yb and 1.0 × 10 6 171 Yb atoms from a magneto-optical trap (MOT) into a crossed optical dipole trap.We perform forced evaporation by lowering the optical dipole trap.Towards the end of the cooling sequence, we remove any remaining 174 Yb atoms with a resonant "push" pulse on the broad 1 S 0 → 1 P 1 transition.After evaporative cooling, typically 3-5 × 10 3 171 Yb atoms in each of the two nuclear spin states m F ∈ {−1/2, +1/2} remain in the optical dipole trap at a temperature of T 0.15T F , which is determined by fitting density profiles to in-situ absorption images.Here, T F is the Fermi temperature.We note that the extremely small scattering length All parameters are estimated from the voltages controlling the lattice depth and magnetic field strength in the experiment.Time is given relative to the clock-laser excitation pulse which quenches the hopping ratio from tH /tL = 1 (not visible on the chosen scale) to tH /tL 0.05, indicated with a small circle.
of 171 Yb in the ground state [2] generally leads to a particularly long thermalization time scale.Here, a 0 is the Bohr radius.
The atoms are loaded into the optical lattices sequentially with s-shaped ramps, and we first ramp up the vertical state-independent lattice to ≈ 30 E m rec within 120 ms.Then, the state-independent lattice along the other axis is ramped up to the same depth within 300 ms.After subsequently loading the atoms into the ≈ 7E rec deep state-dependent lattice (SDL) within 300 ms, we proceed with the state preparation (also see main text), which includes changing the SDL depth within 25 ms (linear ramp) and the magnetic field within 70 ms (s-shaped ramp) to set the desired Hubbard parameters.We show the time dependence of these values for an exemplary case in Fig. S1 and also refer to Table S1, which summarizes the experimental parameters for all figures in the main text.
We use a linearly polarized clock laser beam to excite atoms from the 1 S 0 , m F = +1/2 ground state (de-arXiv:2011.12411v2[cond-mat.quant-gas]17 Aug 2022 TABLE S1.Experimental parameters for each figure presented in the main text.The lattice depth is given for the light species [VH = 3.06(4)VL for the heavy species unless noted otherwise, see Section S.I F].We also show next-nearest neighbor and perpendicular hopping amplitudes to illustrate the accuracy of a one-dimensional tight-binding description.19 [8] 19 [9] 17 [7] 18 [8] Lattice a For this measurement, we operate the state-dependent lattice at the wavelength λ 690.1 nm such that V H = 1.97(5)VL .
noted |g ↑ ) to the 3 P 0 , m F = +1/2 clock state (denoted |e ↑ ≡ |H ).This laser beam co-propagates with one of the deep perpendicular lattices to drive the transition in the Lamb-Dicke regime.To overcome the inhomogeneous broadening present in the SDL and to quench the hopping amplitude for the heavy species, we use a high-intensity excitation pulse with Rabi frequency Ω 2π × 3 kHz.We address singly and doubly occupied lattice sites with this single-frequency pulse by applying a large magnetic field 1533 G corresponding to the zero crossing of the orbital Feshbach resonance (see Section S.I E), where the interaction shift of doubly occupied lattice sites vanishes.In general, the fidelity of the excitation pulse is fairly good, with > 95% of the atoms driven to |H .We note that the residual number of remaining |g ↑ atoms only interact with atoms in |H (interaction strength U H,g↑ t L ).However, this is independent of the chosen heavy-light interaction strength and should only have a minor influence on our observations.After quenching the hopping ratio t H /t L of the noninteracting mixture, we slowly ramp the interaction strength towards the desired value by varying the magnetic field accordingly (see Fig. S1).In general, this process can significantly change the initial state, particularly the fraction of doublons as the Hamiltonian is gradually modified during the ramp of the interaction.However, numerical simulations suggest that this fraction only weakly depends on the chosen U parameter [see Section S.II B 1 for details].We also verify this experimentally by determining an estimate of the doublon fraction with clock-line spectroscopy.After the statepreparation protocol (see Fig. S1), we freeze the motion of both species by rapidly ramping up the SDL to V L ≈ 15E rec .Subsequently, we ramp the magnetic field to B 1450 G and employ clock-line spectroscopy to determine the amplitude ratio of the signal originating from singly and doubly occupied lattice sites.To circumvent the inhomogeneous broadening in the SDL and increase the signal, we only consider a central region of the atomic cloud containing ≈ 5% of atoms.In this way, we determine a maximum change of the doublon fraction ≈ 16% between U/t L ≈ −20 and 4, which is qualitatively consistent with the numerical simulations shown in Fig. S8.

B. Lattice loading
For the theoretical description of the experimental procedure (see Section S.II), the distribution of atoms across the different one-dimensional systems (tubes) is of central importance.Since our absorption imaging inherently averages across one axis, we cannot determine this quantity directly.Instead, we estimate it from parameters in the dipole trap, together with the assumption of constant entropy during the process of lattice loading [3].This assumption only approximately holds as we can detect an increase of the temperature (and entropy) by ≈ 30% when loading the atoms from the lattice back into the dipole trap.We also simplify the lattice loading process by assuming the state-independent lattices are ramped up simultaneously (see Section S.I A).
From the total number of light atoms N tot , the trapping frequencies (ω x , ω y , ω z ) = 2π × [38 (1), 35(1), 402(1)] Hz and temperature T 0.15T F in the optical dipole trap, we determine the total entropy S(N tot , T ) N tot × 1.4k B of the initial state.For the final state in the array of tubes formed by the perpendicular lattices, we fit the chemical potential µ 0 and temperature T as parameters of the total entropy, under the constraint of constant entropy, S = S , and conservation of atom number, i,j N ij = N tot .Here, S ij is the entropy and N ij the number of atoms in the (i, j)-th tube, described by a harmonically trapped onedimensional Fermi gas.The chemical potential µ ij accounts for the local offset of the harmonic potential from the optical dipole trap, with m the mass of a 171 Yb atom, d m = λ m /2 the spacing of the state-independent lattice, and ϕ i , ϕ j corresponding to the relative phase of the optical lattice and optical dipole trap.In Fig. S2(a), we plot the resulting distribution of {N ij } for typical parameters in our experiment.From this distribution, we also calculate the effective quantities shown in Table S1 and Fig. S2(b), For our theoretical calculations (see Section S.II A), we distribute {N ij } over six bins such that the effective width of each bin is approximately constant and average over ten relative phases for each ϕ i and ϕ j [see Fig. S2(c)].We only consider the ≈ 300 significantly filled tubes (N ij ≥ 1) in this procedure, which account for more than 98% of the signal.Finally, we calculate all relevant parameters in the state-dependent lattice using the atom numbers {N ij } and entropies {S ij }, which we increase by 30%/2 = 15% to account for half of the temperature increase observed when loading back into the optical dipole trap.From this estimate and the exact diagonalization of the singleparticle Hamiltonian in Eq. (S.16), we find an effective filling n 0.5, fraction of doublons D 0.3, and temperature T 2t L /k B for a typical total number of light atoms N tot = 4.5 × 10 3 and lattice depth V L = 4.7E rec .In addition, we estimate the effective system size l ≈ 30 lattice sites of a typical tube from the root-mean-square width of the density distribution along the one-dimensional system, which is also obtained with exact diagonalization techniques (see Section S.II A).
To verify our modeling of the lattice loading process, we compare the in-situ column density to a theoretical prediction obtained from the atom numbers and entropies in each tube.Here, we find qualitative agreement with deviations most likely originating from the finite imaging resolution and imperfect confinement potentials.

C. Absorption imaging
We perform in-situ absorption imaging of the atomic cloud on the broad 1 S F =1/2 0 transition, which allows us to directly determine the column density of the light atoms.In the experiment, we freeze the atomic density by rapidly increasing the lattice depth to ≈ 15E rec within 2.5 ms.Then, we lower the magnetic field within 30 ms to 2 G and image the light atoms with a 5 µs long high-intensity pulse.To calibrate the atomic crosssection and the absolute atom number, we use a technique analog to the one used for 173 Yb in Ref. [4].
Finally, we briefly comment on the indirect imaging of the heavy atoms, which are detected after the first imaging pulse has removed all light atoms in the 1 S 0 state.We employ a 3 ms long laser pulse resonant with the 3 P transition to repump all clock state (heavy) atoms back to the ground state.Subsequently, these atoms are detected with absorption imaging on the same transition as the light atoms.However, this detection technique does not reliably determine the in-situ density of the heavy atoms for several technical reasons.First, the repumping introduces motional blurring and additional atom loss up to 30% [1], which depends significantly on the local density.Second, on-site heavy-light atom pairs can either be projected onto the repulsive branch or are continuously transferred into the deeply bound molecular state [1] during the relatively quick ramp of the magnetic field before the imaging.The exact process strongly depends on the technical details of the ramp, and our data suggests that heavy atoms transferred to the molecular state are dark for our detection method, thereby introducing additional loss on the order of ∼ 10%.Although the generally ≈ 40% lower number of detected heavy atoms [see Figs.S7(c) and S7(d)] is approximately consistent with our estimates, the measured column densities still contain systematic errors, and we therefore do not use them for any quantitative analysis.

D. Hubbard parameters
We determine the hopping amplitudes t L and t H from two separate tight-binding band structures for light and heavy atoms [5].The on-site interaction U is determined from the overlap of the Wannier functions w α (x) and the magnetic field-dependent s-wave scattering length a (see Section S.I E), Here, r = (x, y, z), m is the mass of a 171 Yb atom, and w ⊥ (y, z) = w y (y)w z (z) is the combined Wannier function along the perpendicular axes y and z.The main experimental input parameters for the band-structure calculations are the lattice wavelength and depth.The former is fixed with high precision by locking the lasers to a wavelength meter using a digital servo.In contrast, the latter is calibrated for each measurement using deep lattices to minimize systematic uncertainties.For this calibration, we employ lattice modulation spectroscopy [6] for one of the state-independent lattice axes and resolved sideband spectroscopy on the clock transition for the remaining two lattice axes (state-independent and state-dependent, also see Section S.I F).From these measurements, we also estimate the uncertainties of the Hubbard parameters.Finally, the value of the confinement parameter κ = mω 2 d 2 is obtained by measuring the trapping frequency ω after suddenly displacing the trap center and fitting the frequency of the resulting center-of-mass oscillations.We perform this measurement with the statedependent lattice turned off as it does not significantly contribute to the confinement.

E. Orbital Feshbach resonance
Details on the orbital Feshbach resonance (OFR) in 171 Yb can be found in Ref. [1].Here, we only briefly discuss how we employ clock-line spectroscopy on the singleparticle and two-particle resonance to measure the magnetic field-dependent s-wave scattering length between |e ↑ ≡ |H and |g ↓ ≡ |L atoms.This measurement is FIG.S3.Scattering length extracted from clock-line spectroscopy of the single-and two-particle resonance in a nearly isotropic three-dimensional state-independent optical lattice with depth V0 = 30.7(8)Em rec .The circles correspond to individual data points and the solid line is a fit to Eq. (S.6).We do not plot data points close to a(B) = 0 since the resonances overlap in this region making the fit of the experimental data prone to systematic errors.In the inset, we show the zero crossing (B0 + ∆) in an isotropic state-independent optical lattice with variable depth.Here, the circle corresponds to the data set in the main panel whereas the squares represent separate measurements.The dashed line shows the expected scaling from Eq. (S.7) relative to the data point V0 = 30.7Em rec .
performed in a nearly isotropic three-dimensional stateindependent optical lattice with depth ≈ 30.7(8)E rec to avoid the inhomogeneous broadening present in the statedependent lattice (SDL).From the transition energies at various magnetic fields B, we extract the corresponding scattering length a(B) with an expression similar to Eq. (S.5).In Fig. S3, we show the experimental data together with a fit to the following general form [7], Here, the fit parameters are a bg = 333(14)a 0 , ∆ = 255(6) G, and B 0 = 1285(5) G with a 0 the Bohr radius.We find a significant dependence of the zero crossing B a=0 = B 0 + ∆ (calculated from the fitted parameters) on the chosen lattice depth (see inset of Fig. S3).The effective shift of B a=0 most likely originates from the energy dependence of the open-channel scattering amplitude [8].For the lattice depth V 0 , the zero-point energy of the ground band is approximately given by and the resulting shift of the zero crossing can be obtained from δB a=0 = 0 /|δu| with δu = −h × 399.0(1)Hz the differential Zeeman shift of the open and closed channel [1,8].In this simple approximation, we neglect any energy dependence of the scattering length itself but still find good agreement between experiment and theory.The situation for the SDL is slightly different, and we have to consider the mean zero-point energy of a light and a heavy atom in combination with different optical lattices.We utilize the symmetry (with respect to U = 0) of the signal in Fig. 3 and benchmark our approach with additional measurements for different depths of the SDL [see Fig. S5(c)].Since we do not know the precise location of the zero crossing without the trap, we extract it from a fit to this data.Finally, we briefly comment on the lifetime of the heavy-light mixture close to the OFR, which we can monitor up to a maximum duration of 400 ms in our experiment.In general, losses over this time scale are relatively small for the regime 1450-1600 G but increase quickly for magnetic fields < ∼ 1450 G. Importantly, these losses should be negligible for most magnetic fields employed for the measurements in the main text [also see Figs.S7(c) and S7(d)].

F. State-dependent lattice
Our state-dependent optical lattice (SDL) is generated by retro-reflecting a monochromatic laser beam with wavelength λ = 671.509(1)nm, almost identical to the realization in Ref. [9].At this particular wavelength, the polarizability α of |g (light) and |e (heavy) atoms differs significantly, giving rise to different lattice depths and hopping amplitudes.We determine the corresponding polarizability ratio p = α e /α g = V H /V L with interleaved sideband spectroscopy on the clock transition.Here, V H and V L are the lattice depths experienced by heavy (|e ) and light (|g ) atoms, respectively.
For this measurement, we use a separate clock laser beam co-aligned with the SDL axis to change the orbital state from |g to |e (or vice versa) and simultaneously the band index from m to m .We perform two sets of measurements, one with the atoms initially prepared in |e, m = 0 and the other with the atoms in |g, m = 0 .By driving the carrier (m = 0 → m = 0) as well as the first sideband transition (m = 0 → m = 1), we find two transition frequencies by fitting the line shapes as shown in Figs.S4(b) and S4(c).Subtracting the two frequencies yields the band gap individually for |g and |e atoms, which allows us to determine the corresponding lattice depths V L and V H from a band-structure calculation.With these results, we find a polarizability ratio p = 3.06(4) from multiple measurements.By employing the same technique, we also determine the polarizability ratio p = 1.97 (5) for the SDL wavelength λ = 690.076(1)nm.
The much larger polarizability α e of the clock state 3 P 0 (|e ) arises since the SDL is operated only ≈ 20 nm detuned from the 3 P 0 → 3 S 1 transition wavelength.Consequently, the scattering of lattice photons becomes significant and heavy atoms off-resonantly excited to the 3 S 1 state quickly decay back to either of the lower-lying 3 P J=0,1,2 states with a ratio of 15:40:45 according to the Clebsch-Gordan coefficients.Thus, the scattering of only two photons is sufficient to almost completely deplete the |e state, and the dominant dissipation in the system is the loss of heavy atoms.More precisely, heavy atoms are converted with almost equal probability to either 3 P 1 atoms with a subsequent quick decay back to |g ↑ or 3 P 2 atoms, which are most likely quickly lost from the trap.We characterize the loss dynamics by monitoring the remaining number of heavy atoms after preparing a pure sample in |e ↑ = |H .For the relevant lattice depths (3-7E rec ), we find 1/e lifetimes ∼ 10 3 h/t L obtained by fitting the experimental data.While these loss rates are negligible for probing the transport directly after the displacement of the trap minimum in Fig. 3, the total loss becomes appreciable with ≈ 20% for the longest observation times 400h/t L of the density dynamics in Fig. 4.

G. Displacement of the trap center
We use an off-centered state-independent dipole trap beam operated at the magic wavelength to displace the minimum of the trap.The spatially-dependent potential resulting from this approximately Gaussian laser beam is given as Here, m is the mass of a 171 Yb atom, ω ∆ ∈ 2π ×[0, 23] Hz is the trapping frequency proportional to the chosen intensity of the laser beam, x ∆ ≈ 50 µm is the displacement of the beam with respect to the initial trap minimum i 0 (see main text), and w 0 ≈ 1.2x ∆ is the estimated waist of the Gaussian laser beam.Neglecting the higher-order corrections to U ∆ (x), the shifted location of the trap center x 0 is given as with ω = 40(1) Hz (see main text).For our measurements discussed in the main text, the intensity of the off-centered dipole trap is linearly increased, such that the trap minimum is displaced accordingly.We note that the terms in Eq. (S.8) also modify the effective trapping frequency ω, which we estimate to be on the order of ∼ 10% and neglect in our theoretical description.We calibrate the final value of ∆x independently by ramping up the dipole trap before loading the statedependent lattice.This allows us to extract the trap center and effective displacement from a Gaussian fit of the atomic column density.In Fig. S5(a), we show the result of such fits for variable intensity of the optical dipole trap.We estimate the statistical uncertainty of ∆x conservatively with ≈ 3d from the observed 2σ-fluctuations of the trap center over typical measurement durations.

H. Transport measurement
In this section, we discuss details and show additional experimental data for the transport measurements.As discussed in the main text, we observe a mostly symmetric signal (see Fig. 3), which only depends on the absolute value of the interaction strength |U |/t L , shown explicitly in Fig. S5(b).Formally, the dynamical symmetry of our model in the limit κ = 0 [10] is broken by the finite confinement.Due to the relatively small energy scale of κ, we nevertheless still expect the dynamics to be comparable for attractive and repulsive interactions, which is confirmed by our numerical simulations shown in Fig. 3 [also see Section S.II B 2].
To determine the optimal speed for the translation of the trap minimum in our measurement, we also characterize the transported fraction N r /N for constant ramp distance ∆x and variable ramp duration τ ramp .For the probed interaction parameters U/t L 0, −2, −10, we find a monotonous increase of the observable, which approximately follows an exponential saturation curve, as shown in Fig. S5(d).For the measurements in the main text, we chose the ramp speed ≈ 0.5dt L /h, which corresponds to approximately twice the 1/e time in the noninteracting case.We note that the transported fraction N r /N remains significantly reduced for all finite interactions and independent of the exact ramp duration.
To complement the data in Fig. 3, we perform an additional measurement for a mass-balanced configuration with t ≡ t L = t H .This is achieved by replacing the SDL with a state-independent lattice operated at the magic wavelength with the same lattice depth for the |L as well as |H state.We emphasize that this measurement only serves as a qualitative comparison since the beam parameters of the two lattices differ strongly, which introduces systematic differences.In general, we observe a much smaller suppression of the transported fraction and a large asymmetry between attractive (U < 0) and repulsive (U > 0) interactions [see Fig. S5(e)].In contrast to the situation in the SDL, the fraction of doublons here strongly depends on the interaction parameter due to the nearly adiabatic state preparation.Thus, the inhibited transport originates most likely from the slow dynamics ∝ t 2 /U of doublons in the system, as similarly observed in Ref. [11].We verify the increased (reduced) doublon fraction for attractive (repulsive) interactions with clockline spectroscopy and find a modulation of up to ≈ 50%in reasonable agreement with the observed asymmetric lineshape.

I. Density dynamics
We characterize the non-equilibrium dynamics after the displacement of the trap minimum for variable hold times by recording the density deviation, , (S.10) as discussed in the main text.The effective weight n(x, τ ) helps to significantly reduce the contribution of noise in the integral.The one-dimensional density n(x, τ ) is calculated from the two-dimensional atomic column density CD(x, y; τ ) by normalization and integration along y, n(x, τ ) = 1 dxdy CD(x, y; τ ) dy CD(x, y; τ ).(S.11) The relaxation of n(x, τ ) shown explicitly in the two right panels of Fig. 4(a) for t H /t L 0.1 and U/t L −10, 0 can also be observed for other interaction parameters and hopping ratios, which are shown in Fig. S6.Compared to smaller interaction strength or larger hopping ratio, the one-dimensional densities for t H /t L 0.1 and U/t L −10 appear to not show significant relaxation for short hold times τ < ∼ 50h/t L .Importantly, the signatures of δn(τ ) discussed in the main text can also be similarly found in the fraction of transported atoms N r /N .However, this observable is particularly sensitive to the exact calibration of the trap minimum i 0 (see Section S.I G) as the atomic cloud does not fully separate for the relatively short transport distance ∆x 20d in this measurement.Additionally, we note that the center of mass determined for each n(x, τ ) shows interaction-dependent dynamics with similar features as δn(x, τ ).
We estimate the uncertainty of δn(τ ) with jackknife resampling.To this end, we exclude one out of M total measurements and determine n(x, τ ) from the average of the remaining M − 1 measurements.This yields M distinct one-dimensional densities n(x, τ ) from which we obtain the density deviation δn(τ ).From the resulting distribution {δn(τ ) (i) } i=1,...,M , we obtain an uncertainty of δn from the jackknife estimate of the standard error [12] where δn(τ ) = M i=1 δn(τ ) (i) /M .Since the observable δn(τ ) rectifies any differences between n(x, 0) and n(x, τ ), it is susceptible to technical noise.To estimate the contribution of noise to our observable, we consider the density variance δn 2 and assume noise n with vanishing mean n = 0, finite variance n 2 > 0, and vanishing third moment n 3 = 0. Here, X denotes the expectation value of the quantity X.From these assumptions and Eq.(S.10), we find that the contribution of n to δn 2 is proportional to n 2 .We therefore expect for uncorrelated n , Here, γ is a constant, M is the number of individual measurements averaged for the calculation of n(x, τ > 0) to the measurement-independent observable δn 2 .The value δn 2 M for M ≤ 5 can be estimated from subsampling since we take between three and five individual measurements for each parameter.This allows us to extract the otherwise unknown value of δn 2 and γ from a fit.To avoid systematic errors, we use non-interacting data points in the stationary regime and find the noise contribution δn M =3...5 (τ ) − δn(τ ) ≈ 0.4 × 10 −3 , which is used to indicate the band of the theory data in Fig. 4. We emphasize that this approach only serves as an estimate since it neglects potential correlations of n .
In Figs.S7(c) and S7(d), we show the number of heavy and light atoms during the measurement of the density dynamics.The initial atom numbers and loss rates of the heavy species are generally quite comparable for the different interaction parameters.
As discussed in the main text, we fit the experimental data shown in Fig. 4 with an exponential saturation curve, given by Here, a, c are dimensionless fit parameters, and τ 0 is the fitted exponential decay constant, which is discussed in the main text.By excluding each point individually in a separate fit routine [jackknife resampling, see Eq. (S.12)], we generate a distribution of τ 0 from which we determine the estimated uncertainty of τ 0 given in the main text.
To compare the two cases for t H /t L 0.1 and finite interaction strength, we fit the measured density deviation to a linear function, Restricting the numerical fit to early times, τ ≤ 101h/t L , we find the fitted parameters m = 13(2) × 10 × 10 −6 t L /h and c = 1.2(2) × 10 −3 t L /h, 1.3(1) × 10 −3 t L /t H for the cases with U/t L −2 and −10, respectively.The same procedure yields m 14 × 10 −6 t L /h and c 0.9 × 10 −3 t L /h for both finite interaction strengths with t H /t L 0.34.

S.II. THEORETICAL DESCRIPTION A. Non-interacting calculation
For the theoretical calculations in the non-interacting limit (see Figs. 2 and 4), we employ the following singleparticle Hamiltonian which we numerically diagonalize for a finite but large system size l = 250 (i ∈ [−l/2, l/2]) to find the eigenenergies k and eigenstates with time evolution |ψ k (τ ) = e −i k τ /h |ψ k .To account for the finite temperature T of the system (see Section S.I B), we calculate the occupation of the different eigenstates from the Fermi-Dirac distribution Here, T is the temperature, and the chemical potential with N L the total number of light atoms in one tube.Finally, we calculate the time-dependent density profile from which we determine the experimentally relevant quantities.
The ramp of the trap minimum in the experiment is simulated by integrating the time-dependent Schrödinger equation corresponding to the Hamiltonian (S. 16) for all eigenstates.For the long evolution times in Fig. 4, we average the calculation over ten equally spaced noninteger offsets ∈ [−0.5, 0.5) of the trap minimum relative to i 0 .This is motivated by the fact that the relative phase of the lattice with respect to the minimum of the harmonic confinement is not fixed between individual shots of the experiment and also varies across the different tubes.To approximately model the contributions of tubes with different atom numbers in the experiment, we use a binned distribution with different weights (see Section S.I B) calculated for the specific total number of light atoms N tot (see Table S1) and estimate our observable from a weighted average of these bins.The theory curve in Fig. 2 systematically overestimates N r /N , which can be corrected with a ≈ 25% increase of the estimated temperature or atom number.We observe a similar deviation for the experimental and theoretical data in Fig. 4. Overall, we attribute the disagreement to the approximative nature of our modeling of the lattice loading process.

B. Many-body calculation
To study the interacting Hamiltonian defined in Eq. ( 1) we employ tensor network methods as we aim to closely predict the dynamics observed in the experiment.To account for the finite temperature of the system, we apply a density matrix-based approach, where the density matrix ρ of the closed system is approximated as a matrix-product operator (MPO).By means of purification, this MPO can be expressed as pure matrix-product state (MPS) |ρ in an enlarged Hilbert space, such that we obtain the density matrix by tracing out the ancilla degrees of freedom ρ = Tr a |ρ ρ|.In this representation, we can make use of established tensor network algorithms developed in the MPS framework [13,14].In particular, we implement our code based on the TenPy package [15].To calculate the thermal equilibrium state ρ ∝ e −β Ĥ we start with the maximally mixed, infinite temperature state ρ0 ∝ 1, which can be represented exactly as a trivial MPO.For our nearest-neighbor Hamiltonian we employ the time evolving block decimation (TEBD) algorithm [16,17]

Initial state
In general, we expect that the details of our initial state preparation in the experiment (see Section S.I A) have an important influence on the dynamics.Therefore, we study the preparation protocol numerically by approximating the time dependence of the Hubbard parameters [see Fig. S1] with a sequence of linear ramps and intermediate hold times.We integrate the von Neumann equation for the time-dependent Hamiltonian by applying the TEBD algorithm with a first-order Trotter decomposition, while updating the time evolution operator within each simulation step.This procedure resembles the simplest possible unitary integration scheme within our MPS approximation.We start with a non-interacting thermal state at T init = 4t L /k B to approximately model the initial state in the ≈ 7E rec deep state-independent lattice, quench the hopping ratio t H /t L = 1 → 0.05, apply the ramps for the confinement strength as well as hopping ratio, and subsequently turn on the interaction by ramping U/t L to the desired value.
In Fig. S8, we show the density profiles niα of the light atoms (α = L) and heavy atoms (α = H) as well as the density of doublons niL niH after calculating the state preparation for an exemplary N L = N H = 10 tube with U/t L = −10 and 10.Here, Ô = Tr(ρ Ô) denotes the expectation value of the corresponding operator Ô.We compare the prepared state to a thermal equilibrium state at fixed temperature T eff for strong interactions and for the non-interacting case.For this purpose, we choose the effective temperature of this state according to T eff = (κ 2 /κ 1 )T init 2.3t L /k B where κ 1 (κ 2 ) is the trapping potential before (after) the state-independent optical lattice is ramped to the final depth (see Fig. S1).As shown in Figs.S8(a) and S8(b), the density profiles of light and heavy atoms for the prepared state are similar to a thermal state at temperature T eff on the repulsive side, while significant differences become visible for the strong attractive case.Furthermore, we put emphasis on the doublon density niL niH and relative fraction Due to the expected exponentially long lifetime [18] and suppressed mobility of doublons [19] (τ D ∝ t 2 /U , for the mass-balanced case t = t L = t H ), we expect the initial fraction of doublons to influence the transport properties in our system critically.Strikingly, the doublon fraction drastically deviates from the thermal state, which has a significantly higher (lower) value on the attractive (repulsive) side [see inset of  and (c),(d) doublons at the end of our state preparation protocol from a numerical simulation for a single tube with NL = NH = 10 atoms.Solid colored lines correspond to the full state preparation protocol for U/tL = −10 (attractive, blue) and 10 (repulsive, red) with typical experimental parameters, while the dashed lines represent the corresponding thermal state at the T eff obtained from the localdensity approximation and at finite interactions U/tL.The black solid lines correspond to a non-interacting thermal state at temperature T eff .In the inset of panel (d), we show the doublon fraction D for variable interaction strength U/tL and similar preparation protocols as in the main panels.Here, the solid black line corresponds to the non-interacting thermal state with an interaction quench 0 → U/tL and a subsequent hold time of 10h/tL.suggest that the system does not thermalize during the state preparation and the resulting state is still close to a non-interacting state at T eff .Moreover, the doublon fraction only weakly depends on the interaction strength, which is consistent with the symmetry between the attractive and repulsive side observed in the experiment [see Fig. 3].

Simulation of the transport measurement
Simulating the full dynamical response of our system even after displacing the trap minimum is generally limited by the rapidly growing operator space entanglement.To numerically access long enough times before truncation errors become significant, we do not perform the full preparation protocol but rather choose the noninteracting thermal state at T eff as an approximate initial state and quench to the desired interaction strength.Comparing the density profiles and the doublon fraction verifies that this approximates the prepared state reasonably well (see black lines in Fig. S8).To further mitigate truncation errors due to growing entanglement, we only consider N L = N H = 5 atoms in this simulation.Fig. S9 shows the transport after displacing the trap minimum over a distance of 45d within 90h/t L for different interaction strength on the attractive side.The suppression of mobility with increasing |U |/t L becomes very clear and can be quantified by measuring the fraction of atoms on the right side of the system with i ≥ 0, in analogy to the measurement of N r /N in the experiment.The results for attractive and repulsive interactions are plotted in the inset of Fig. 3, and qualitatively agree with our experiment.

C. Many-body Wannier-Stark localization
As discussed in the main text, harmonic confinement in the experiment leads to Wannier-Stark localization of single-particle eigenstates with energies above the band edge at 2t L .Hence, transport in the non-interacting case is strongly suppressed, especially for the heavy atoms, as these experience a much stronger effective confinement κ/t H κ/t L compared to the light atoms.For the interacting case, we expect some of these states to be delocalized through interactions [11,20] already in the early time dynamics of our system.Nevertheless, manybody Wannier-Stark localization [21] is still expected for strong enough trapping and has been studied, e.g., in Ref. [22] for the mass-balanced Fermi-Hubbard model (albeit relaxation at extremely late times cannot be fully ruled out).Here, we extend these studies to the mass- imbalanced case.To this end, we study the quench dynamics for an initial staggered product state in a system of 64 sites for a typical hopping ratio of t H /t L = 0.1.We express the initial state as an exact MPS and calculate the time evolution with the TEBD algorithm.Fig. S10 depicts the dynamics of the atomic densities for both species.The light atoms equilibrate on a time scale h/t L , and we observe relaxation for the given system size and confinement κ/t L = 2 × 10 −2 .In contrast, relaxation of heavy atoms starts at a time scale h/t H and occurs in the center of the trap.Close to the edges of the system, the local gradient leads to a confinement of the heavy atoms on the simulated time scales.For the massbalanced Fermi-Hubbard model, the local critical potential has been estimated as ∂ i H/n iα = κ (i crit − i 0 ) 2.8 [22] which agrees well with our observations in the mass-imbalanced case.We therefore estimate that already at short times the heavy atoms are mobile over a central trap region of 2i crit ≈ 30 lattice sites.Our simulations explicitly focus on weak interactions U/t L = −1.0 to ensure that dynamical constraints arising from doublon physics and mass imbalance are weak, and we consequently only probe the influence of harmonic confinement on the mobility of the heavy atoms.

D. Time scales
It has been discussed that strong mass imbalance in small translational-invariant (κ = 0) systems leads to an ergodic regime of extremely slow equilibration [23].For the delocalized atoms in the center of the trap, we expect the dynamics to be close to this κ → 0 limit.This suggests that our finite-size experimental system with large mass imbalance will feature similar physics.To study the involved time scales, we turn to very large mass imbalance and long evolution times.For this purpose, we perform exact diagonalization of the translational-invariant model with the Hamiltonian in Eq. (1) taking the form, subject to periodic boundary conditions, ĉ † lα = ĉ † 0α , for a system of size l = 8 at half filling N L = N H = 4.To estimate the relaxation time of a density modulation, we calculate the dynamical correlator which measures the decay of an inhomogeneous density perturbation nqα = j e −iqj njα for one species (α = L, H).
The decay of long-wavelength modulations is shown in Fig. S11(a) for t H /t L = 10 −3 and indicates the existence of distinct time scales governing the dynamics.These time scales also become apparent in the entanglement dynamics corresponding to plateau-like features of The wavevector q correspond to the longest wavelength λ = 2π/q supported by the respective system size under periodic boundary conditions.Here, the strongly interacting case at U/tL = −10 with the mass ratio tH /tL = 10 −2 is considered.Note that the features of the decay dynamics become less pronounced with increasing system size.
very slow entanglement growth [see Fig. S11(b)].Here, we compute the half-chain entanglement entropy for random initial product states, averaging the results of 100 states.After initial relaxation of the light particles on the time scale ∼ h/t L the dynamics stagnate, and the entanglement entropy stays nearly constant up to time ∼ h/t H .Following further relaxation and an increase of the entanglement entropy, the system enters a new regime at time scale ∼ hU/t 2 H [dashed lines in Figs.S11(a) and S11(b)].This regime features slow decay of non-equilibrium states [23], and we observe an accordingly slow growth of entanglement entropy.Although the phenomenology shows similarities to a dynamical type of many-body localization, the system is still ergodic, and density perturbations eventually decay.From our exact diagonalization data, we extract the time scale τ = C × h|U |/t 2 H for the onset of this regime and estimate C = 0.29(3) from a fit [see Fig. S11(c)].
We note that for increasing t H /t L and increasing system size, the features separating different regimes become less pronounced.In Fig. S12 we illustrate the finite-size scaling for systems of size l = 8, 10, and 12.When the system size is increased new time scales enter the dynamics and additional features in the relaxation dy- namics and entanglement growth become visible.This is exemplified for l = 10 in Fig. S12(b), where the time scale ∼ hU t L /t 3 H is observed.Such a hierarchy of time scales "washes" out distinct plateau-like features, which for large systems leads to featureless decay dynamics.Note that distinct time scales in the relaxation dynamics are only expected for small systems or extremely small mass ratios t H /t L [24].
Few-body bound states constitute a promising explanation for the observed time scales in small systems.A bound state of N H heavy and N L light atoms is stable for t H /t L < 1, and can become mobile at time ∼ h U N L t N H −N L −1 L /t N H H , bottlenecking the overall dynamics.The emergence of such bound states could explain the unconventionally slow relaxation, where a cascade of time scales leads to the observed dynamics.

E. Properties of the transient regime
In our experimentally realized Fermi-Hubbard model, mass imbalance between the two species induces strong dynamical constraints on the light atoms.Contrary to the small size exact diagonalization studies discussed in the previous section, the experiment is conducted at moderate mass-ratios t H /t L ≥ 0.05, and in a significantly larger system than accessible to exact diagonalization.In this case, we expect overall slow relaxation dynamics, but no distinct time scales in form of plateau-like features.
Here, we briefly discuss the transport properties of our system in this transient regime of unconventionally slow relaxation.To this end, we consider the Hamiltonian in Eq. (1) for mass ratios t H /t L ≥ 0.1 with vanishing confinement κ/t L = 0, and moderate attractive interactions U/t L = −4.Within the accessible simulation time, we find the features to be most pronounced for this intermediate interaction strength, while other values of U/t L lead to qualitative similar behavior.To characterize transport properties, we calculate the density-density correlation function C iα (τ ) = niα n0α (τ ) . (S.26) We simulate the Heisenberg time evolution of the operator n0α with the same method as discussed for the density matrix (see Section S.II B).From n0α (τ ) the full spatial correlation profile C iα (τ ) can be obtained due to translational invariance.We note, that for the Heisenberg equation of motion ∂ τ Ô = i[ Ĥ, Ô], the superoperator reads L = − Ĥ ⊗ 1 + 1 ⊗ Ĥ.The simulations are carried out for a system of size l = 201, which has been verified not to exhibit finite-size effects on the accessible simulation times.
Transport properties can be characterized by the scaling of the mean squared displacement which quantifies the spatial width Σ α (τ ) of the correlation profile.For diffusion, we expect the spatial profiles to be well described by Gaussians [25], with the spatial width scaling as Σ α (τ ) ∝ τ 1/2 .Generally, a power law scaling Σ α (τ ) ∝ τ 1/z indicates ballistic and diffusive transport, for z = 1 and z = 2, respectively, superdiffusion for 1 < z < 2, and subdiffusion for z > 2 [25].In Fig. S13(a), spatial profiles of the density-density correlator corresponding to the light species are shown for different mass ratios and at different times.For small mass ratios we find that, within the accessible simulation time, the spatial profiles cannot be described by Gaussians, but rather have a peaked shape in the semilogarithmic plots.This observation indicates anomalous transport in the system, and has previously been found in small systems with exact diagonalization [26].
Moreover, the broadening of the spatial width with time slows down with increasing mass imbalance.To characterize this regime, we extract the dynamic exponent z from the logarithmic derivative of Σ L (τ ) averaged over time windows of fixed size ∆τ = 10.As shown in Fig. S13(b), this reveals subdiffusive behavior z > 2 of our system on the accessible simulation times.However, we observe that the dynamic exponent z is timedependent and decreases with time.The value of z approaches 2 at late times, which indicates a crossover from subdiffusive (z > 2) to diffusive dynamics (z = 2).These findings suggest that our experimental system exhibits a transient regime of slow subdiffusive relaxation dynamics, as indicated in the suggested dynamical phase diagram depicted in Fig. 1(c).

FIG. 1 .
FIG. 1. Transport in the mass-imbalanced Fermi-Hubbard model.(a) Illustration of the mass-imbalanced Fermi-Hubbard model with on-site interaction U and hopping tL tH of the light (blue circles) and heavy (red circles) particles, respectively.We also show the atomic states relevant for the experimental implementation.(b) Schematic representation of the different steps in our transport measurement (top row) and integrated densities of the light atoms (bottom row) from experimental snapshots.(c) Scenario for the transient regime following a quench in the mass-imbalanced Fermi-Hubbard model for variable inverse hold time (vertical axis).Here, the dashed line marks the crossover between unconventional relaxation dynamics at early times (question mark) and diffusion at late times.The orange arrows indicate a trajectory in the regime of our experiment.

FIG. 2 .
FIG. 2. Wannier-Stark localization in the absence of interactions.(a) Fraction of light atoms Nr/N transported to the right half of the system (blue circles) in the absence of heavy atoms and for variable confinement κ/tL determined by the lattice depth VL.Each point is the average of 2-5 measurements, and error bars indicate the uncertainty of κ/tL (partly smaller than the marker size).The solid line corresponds to our numerical calculation [42], and the insets show raw atomic column densities for κ/tL = 1.7(1)×10 −2 and 9(1)×10 −2 with the dashed lines indicating the boundary that determines the atom count Nr in the right half of the system.(b) Probability density |Ψn(x)| 2 and eigenenergy En/tL of the lowest-lying single-particle eigenstates for VL = 9Erec (κ/tL = 4.9×10 −2 ).We add a small linear potential to lift the degeneracy of the localized states.

3 FIG. 3 .
FIG. 3. Constrained early-time dynamics in the interacting heavy-light mixture.Fraction of light atoms Nr/N transported to the right half of the system (blue markers) for variable heavy-light interaction strength U , fixed hopping ratio tH /tL = 0.104(7), and confinement κ/tL = 1.7(1) × 10 −2 .Each point is the average of 3-4 measurements, and error bars (partly smaller than the marker size) denote the standard error of the mean of Nr/N and the uncertainty of U/tL.The bottom panels show raw atomic column densities for U/tL = −4.8(4),0.02(4), and 2.9(2) (from left to right, hexagonal markers) with the dashed lines indicating the boundary that determines the atom count Nr in the right half of the system.The inset of the main panel shows the numerical matrixproduct-state simulations of a single tube with the Hubbard parameters and confinement strength of the experiment and Nsim = 5 atoms of each species.

FIG. 4 .
FIG. 4. Late-time relaxation dynamics.Density deviation δn(τ ) of the light atoms with respect to the initial atom distribution after translating the trap minimum by ∆x [see Eq. (2)], for a variable hold time τ > 0, as illustrated in the bottom-right schematic of panel (b).The time traces correspond to three heavy-light interaction strengths (see legend) with (a) tH /tL = 0.102(6) and (b) 0.34(3).Each data point is calculated from the average of 3-5 atomic densities and an independent measurement of n(x, 0).Error bars (partly smaller than the marker size) denote the uncertainty of τ , determined from the estimated uncertainty of the hopping amplitude tL, and standard error of δn(τ ), determined from jackknife resampling [42].The colored lines show exponential fits, as discussed in the main text, and in panel (a), a three-point moving average for U/tL −10 as a guide to the eye.We show our numerical calculation for U = 0 as a gray band with the height derived from the estimated systematic uncertainty of δn(τ ) [42].The two central panels show integrated atomic densities, which are used to calculate δn(τ = 330h/tL) in (a).

FIG
FIG. S1.Time dependence of the Hubbard parameters tH , U , and κ during the state preparation for the target values tH /tL 0.1 and U/tL −10.All parameters are estimated from the voltages controlling the lattice depth and magnetic field strength in the experiment.Time is given relative to the clock-laser excitation pulse which quenches the hopping ratio from tH /tL = 1 (not visible on the chosen scale) to tH /tL 0.05, indicated with a small circle.

Fig. 2 4
Fig. 2(a) Fig. 3 Fig.4(a) Fig. 4(b) a Total number of light atoms Ntot 10 3 4.2(3) 4.5(3) 3.4(3) 3.8(6) Effective number of light atoms per tube N [∆N ]19[8]   19[9]   17[7]   18[8] FIG. S2.Distribution of atoms across different onedimensional systems (tubes).(a) Typical distribution of Ntot = 4.5 × 10 3 (light) atoms across the tubes generated by the perpendicular state-independent lattices.Each black circle represents a tube with the size indicating the number of atoms (see legend).We only plot a single quadrant due to the symmetry of the distribution.(b) Effective number of light atoms per tube N for variable total light atom numbers Ntot.The blue-shaded area corresponds to the region of N ± ∆N .(c) Approximate distribution of light atoms across different bins for Ntot = 4.5 × 10 3 (yellow star) and their contribution to the measured signal.This distribution is obtained by averaging over multiple relative phases −0.5 ≤ ϕi, ϕj < 0.5 [see Eq. (S.2)].
FIG. S4.Measurement of the polarizability ratio p in the state-dependent lattice (SDL).(a) Band structure of a 22.6(6)Erec deep SDL with the two lowest-lying energy bands of |g atoms (blue) and |e atoms (yellow).The four transitions required to determine the |g and |e band gap are shown as black arrows at exemplary quasi-momenta q.All energies are shown relative to E0, which is the mean energy of the |g ground band.(b),(c) Clock-line spectra of the carrier and sideband transitions starting from (b) |e, 0 and (c) |g, 0 for variable detuning of the clock laser.Here, the circles correspond to data points and the solid lines to individual fits of the line shapes.
FIG. S5.Details on the displacement of the trap center and the transport measurement.(a) Transport distance ∆x extracted from a Gaussian fit of the atomic column density for variable intensity I∆ in the off-centered dipole trap.The dashed line corresponds to a linear fit.(b) Transported fraction Nr/N from Fig. 3 shown for |U |/tL to illustrate the symmetry between attractive (U < 0, red diamonds) and repulsive (U > 0, blue circles) interactions.(c) Transported fraction Nr/N for variable hopping ratio tH /tL and interaction strength U/tL.Note that the tH /tL = 0.104(7) data set (light blue circles) corresponds to the data shown in Fig. 3.The solid lines are fits to A exp[−(U/tL) 2 /(2σ 2 )] + B for extracting the amplitude A and minimum B shown in the inset together with a linear fit (dashed gray line).(d) Transported fraction Nr/N for tH /tL = 0.102(6), ∆x 43(3)d, and variable ramp duration τramp.The solid lines are a guide to the eye.(e) Transported fraction Nr/N for the special case t ≡ tH = tL and variable interaction strength U/t with the transport distance ∆x 44(3)dm and ramp duration 94(3)h/t.
FIG.S6.One-dimensional densities for the hopping ratios (a) tH /tL = 0.102(6)  as well as (b) 0.34(3) and multiple hold times as indicated in the top right legends (units of h/tL).Each row corresponds to either one of the interaction strengths U/tL = −1.9(3)(light blue, top row of each panel) and −10(1) (dark blue, bottom row of each panel).

2 and δn 2 M 2 FIG
FIG. S7.Atom loss during the measurement of the density deviation δn(τ, x).Number of heavy (|H , top panels) and light (|L , bottom panels) atoms during the measurement of the density deviation for (a) tH /tL = 0.102(6) and (b) 0.34(3) presented in Fig. 4. Data points spaced less than 50h/tL are binned to reduce visual clutter, and the interaction parameters are shown in the legend.
to obtain the equilibrium state at T = 1/(k B β), ρ ∝ Tr a e −β/2 Ĥ |ρ 0 ρ 0 | e −β/2, Ĥ , (S.20) by evolving |ρ 0 up to β/2 in imaginary time.The dynamics of our system are governed by the von Neumann equation ∂ τ ρ = −i[ Ĥ, ρ], which we can write in terms of the purified MPS, ∂ τ |ρ = −iL |ρ , (S.21) with the Liouville superoperator L = Ĥ ⊗ 1 − 1 ⊗ Ĥ.The solution |ρ(τ ) = e −iLτ |ρ can be efficiently calculated by propagating the initial state in real time with the TEBD algorithm.The accuracy of our simulations depends critically on the operator-space entanglement, defined as the half-chain von Neumann entanglement entropy corresponding to the normalized MPS |ρ .Throughout our simulations, we dynamically truncate the singular values at each bond, by keeping the discarded weight smaller than a threshold ε trunc , but keeping at most χ max Schmidt values.Note that in order to keep truncation errors negligible at fixed χ max , our simulations are generally restricted to small atom numbers N L = N H ≤ 10.

1 D
FIG. S8.Density profiles of (a) light atoms (α = L), (b) heavy atoms (α = H), and (c),(d) doublons at the end of our state preparation protocol from a numerical simulation for a single tube with NL = NH = 10 atoms.Solid colored lines correspond to the full state preparation protocol for U/tL = −10 (attractive, blue) and 10 (repulsive, red) with typical experimental parameters, while the dashed lines represent the corresponding thermal state at the T eff obtained from the localdensity approximation and at finite interactions U/tL.The black solid lines correspond to a non-interacting thermal state at temperature T eff .In the inset of panel (d), we show the doublon fraction D for variable interaction strength U/tL and similar preparation protocols as in the main panels.Here, the solid black line corresponds to the non-interacting thermal state with an interaction quench 0 → U/tL and a subsequent hold time of 10h/tL.
FIG. S9.Interaction-dependent transport of light atoms from a numerical simulation for a single tube of NL = NH = 5 atoms, transport distance 45d, and transport duration of 90h/tL.The solid gray line corresponds to the initial state located at the initial trap minimum i0 and blue lines represent the evolved density profiles for different interaction strength after translating the trap minimum to i1.We indicate the initial and final harmonic confinement with dotted gray lines.
FIG. S10.Time-dependent quench dynamics of the atomic densities for an initial staggered state in a quarter filled system of size l = 64, shown for (a) light and (b) heavy atoms.The dashed lines in panel (b) indicate the boundary between localized and delocalized states.The Hubbard parameters are close to typical experimental values with κ/tL = 2 × 10 −2 , tH /tL = 0.1, and U/tL = −1.0.We note that the visible asymmetry on the edges of the system originates from the non-inversion-symmetric initial state.

2 FIG
FIG. S11.Extraction of time scales for large mass imbalance.(a) Decay of a long-wavelength excitation and (b) growth of entanglement entropy for initial product states evaluated for a translational-invariant, periodic system of size l = 8 for NL = NH = 4 atoms.We consider the interaction strength U/tL = −10 for the strong mass-imbalanced case of tH /tL = 10 −3 to illustrate the different time scales.The inset shows the time scale τ = C × hU/t 2 H corresponding to the onset of a regime of unconventionally slow relaxation [dashed line in panels (a) and (b)].Here, the constant C is extracted from a fit displayed as solid lines.We indicate the experimentally relevant value for tH /tL = 0.1 with a star-shaped marker.

12 FIG
FIG.S12.System-size dependence of the dynamical correlator Cqα(τ ) for translational-invariant and periodic systems of size (a) l = 8, (b) l = 10, and (c) l = 12.The wavevector q correspond to the longest wavelength λ = 2π/q supported by the respective system size under periodic boundary conditions.Here, the strongly interacting case at U/tL = −10 with the mass ratio tH /tL = 10 −2 is considered.Note that the features of the decay dynamics become less pronounced with increasing system size.

6 FIG
FIG. S13.Transient transport properties of the mass-imbalanced Fermi-Hubbard model.(a) Spatial profiles of the densitydensity correlation function of the light species for different mass-ratios tH /tL = 0.6, 0.3, 0.1 (left to right) at four points in time.(b) Scaling of the spatial width ΣL(τ ) ∝ τ 1/z obtained from the mean squared displacement.The dynamic exponent z is obtained by power law fits to the time dependent spatial width ΣL(τ ) over time windows of size ∆τ = 10/tL.For ordinary diffusion Gaussian, profiles and dynamical exponents z = 2 are expected, while z > 2 indicates subdiffusion.