Observing non-ergodicity due to kinetic constraints in tilted Fermi-Hubbard chains

The thermalization of isolated quantum many-body systems is deeply related to fundamental questions of quantum information theory. While integrable or many-body localized systems display non-ergodic behavior due to extensively many conserved quantities, recent theoretical studies have identified a rich variety of more exotic phenomena in between these two extreme limits. The tilted one-dimensional Fermi-Hubbard model, which is readily accessible in experiments with ultracold atoms, emerged as an intriguing playground to study non-ergodic behavior in a clean disorder-free system. While non-ergodic behavior was established theoretically in certain limiting cases, there is no complete understanding of the complex thermalization properties of this model. In this work, we experimentally study the relaxation of an initial charge-density wave and find a remarkably long-lived initial-state memory over a wide range of parameters. Our observations are well reproduced by numerical simulations of a clean system. Using analytical calculations we further provide a detailed microscopic understanding of this behavior, which can be attributed to emergent kinetic constraints.

U nderstanding the complex out-of-equilibrium dynamics of quantum many-body systems is central to a number of research areas ranging from statistical physics to quantum information theory [1][2][3] . State-of-the-art experimental platforms are now able to test novel theoretical concepts and approximate descriptions based on experimental observations. Important experimental results were obtained in particular with integrable 4 or many-body localized (MBL) [5][6][7] systems. Both phenomena emerge due to the existence of extensively many conserved quantities and have been of considerable interest, because they break the eigenstate thermalization hypothesis, which assumes that each individual eigenstate behaves locally like a thermal ensemble and is believed to hold for generic ergodic systems [8][9][10] .
In between the two extreme limits of ergodic and localizing dynamics there exists a rich variety of more complex thermalizing behavior. Models with many-body scar states, e.g., host a vanishing fraction of non-thermal eigenstates embedded within an otherwise thermal spectrum [11][12][13][14][15][16] . They exhibit a weak form of ergodicity-breaking, that strongly depends on the initial state, as has been observed with Rydberg atoms 14,17,18 . More recently, a whole new class of models has been suggested, where the presence of only few conserved quantities, in particular dipole conservation, results in non-ergodic dynamics due to an emergent fragmentation of the Hilbert space into exponentially many disconnected subspaces [19][20][21][22] . Fragmented models offer an alternative view on a central open question, namely if many-body localization can occur in translationally-invariant models without disorder [23][24][25][26][27][28][29] .
In this work, we study non-ergodic behavior in the disorderfree tilted one-dimensional (1D) Fermi-Hubbard model (Fig. 1a), which lies at the interface of MBL and Hilbert-space fragmentation. In the presence of additional weak disorder or harmonic confinement, theoretical studies have found characteristic MBL phenomenology, known as Stark MBL [30][31][32][33][34] . This, however, does not hold for a clean system with pure linear potential 30,33 . While conventional MBL predicts localization for any typical initial state, we do not expect this to hold for our system, where resonances can occur between interaction and tilt energies (regime ① in Fig. 1b). Intriguingly, it has been predicted, that in the limit of large tilts, Δ ≫ J, |U|, non-ergodicity may still occur despite the absence of disorder. In this regime, the large tilt energy imposes kinetic constraints, which result in an emergent dipole conservation 19,20,22,31,33 . This emergent behavior is in fact governed by a fragmented Hamiltonian resulting in non-ergodic dynamics. Starting from an initial charge-density wave (CDW) of singlons (singly-occupied site), we study relaxation dynamics in the tilted 1D Fermi-Hubbard model for a large range of interaction strengths and moderate values of the tilt (Δ < 4J), where none of the two mechanisms described above should apply and where naively one may expect the system to thermalize 35,36 . At short times we observe coherent dynamics due to Bloch oscillations, whose amplitude strongly depends on the Hubbard interactions. Surprisingly we find that after intermediate times and even close to resonance (regime ①), the evolution converges to a steady-state, that persists for long evolution times up to 700 tunneling times, signaling a robust memory of the initial CDW throughout.
Using numerical calculations we show that the observed nonergodicity cannot be explained by the phenomenon of Stark-MBL, i.e., the robust memory is not due to experimental imperfections, such as residual harmonic confinement or disorder, and the bipartite entanglement entropy does not exhibit the characteristic behavior of MBL systems 30,37 (Supplementary Fig. 5). Hence, non-ergodicity appears to have a different origin, despite similar experimental signatures. This raises the question about the origin of the observed non-ergodicity. We construct effective Hamiltonians in two distinct regimes (① and ②, Fig. 1b) by taking the large tilt limit and find strongly-fragmented Hamiltonians in both cases (Supplementary Note 3). While these models are only expected to describe the dynamics at large tilt values and for intermediate times (on the order of a few tens of tunneling times), they allow us to identify the microscopic processes that initiate dynamics at short times (Fig. 1b). In both regimes these are correlated tunneling processes, which result in the formation of doublons (doubly-occupied sites), either resonantly (regime ①) or detuned by the Hubbard interaction energy U (regime ②). Higher-order terms are expected to eventually drive the system towards thermalization 19 . However, we are able to show that energy penalties for the second-or higher-order tunneling processes, which occur naturally in the model, render these dynamics inefficient. This results in extremely slow relaxation (Supplementary Note 3), which appears stable for > 10 4 tunneling times in our exact diagonalization studies of small systems, in agreement with our experimental observations ( Supplementary Fig. 4).
In order to characterize the dynamics across the whole parameter regime studied experimentally, we compute the finite-time connectivity of our initial CDW state C ϵ ¼ dimðN ϵ Þ= dimðHÞ, which is defined by the fraction of states that participate in the time evolution up to a finite time T N ; here N ϵ denotes the subspace in the complete Hilbert-space H, which is defined, such that the residual overlap of the time-evolved state ψðtÞ outside of N ϵ is at most ϵ at any time t ≤ T N (Methods). The value of ϵ is typically chosen between 1 and 10%. The finite-time connectivity can be understood as a measure of non-ergodicity, similar to the more conventional return probability or other multifractality measures 38 . While effective Hamiltonians can only be derived explicitly in certain limits, the numerical construction is applicable in the whole parameter regime probed in this work (Fig. 1c). We find that the finite-time connectivity vanishes in the thermodynamic limit for all parameters, suggesting that only a small fraction of the states participates in the dynamics, signaling nonergodic behavior. Our results suggest that the emergent kinetic constraints result in transient non-ergodic behavior across the whole parameter range studied in this work. We further show analytically that the relevant microscopic constraints in the resonant ① regime give rise to Hilbert-space fragmentation in the large tilt limit (Supplementary Note 4). b Dominant resonant tunneling processes for different regimes. c Finite-time connectivity C ϵ (for a cut-off ϵ = 10%) defined as the fraction of states that participate in the dynamics up to an evolution time T N ¼ 1000τ (main text, "Methods"). The calculation was performed for a Néel-ordered singlon CDW initial state, using exact diagonalization (ED) with system size L = 13 and Δ ↑ = Δ ↓ ≡ Δ. In the large-tilt limit, Δ/J → ∞, we find emergent strongly-fragmented effective Hamiltonians for regime ① and ② (see Supplementary Note 3 and 4).

Results
The experimental setup consists of a degenerate Fermi gas of 50 (5) × 10 3 40 K atoms that is prepared in an equal mixture of two spin components " in the F = 9/2 ground-state hyperfine manifold. The atoms are loaded into a 3D optical lattice with lattice constant d s = 266 nm along the x direction and deep transverse lattices, with constant d ⊥ = 369 nm, to isolate the 1D chains along x ("Methods"). The central 1D chains have a length of about 290 lattice sites. The residual coupling along the transverse directions is less than 3 × 10 −4 J. The dynamics along x is described by the tilted 1D Fermi-Hubbard model whereĉ y iσ (ĉ iσ ) is the fermionic creation (annihilation) operator andn i;σ ¼ĉ y iσĉiσ . The on-site interaction strength U is controlled by a Feshbach resonance centered at 202.1 G and a magnetic field gradient is used to create the tilt Δ σ , with Δ ↑ ≃ 0.9Δ ↓ . The weak spin-dependence arises due to the different m F quantum numbers (Supplementary Note 9 and 11). The initial state for all subsequent measurements is a CDW of singlons on even sites, which is prepared using a bichromatic optical superlattice (Supplementary Note 7). The initial state can be described as an incoherent mixture of site-localized particles with random spin configuration ("Methods"). The subsequent evolution is monitored by extracting the spin-resolved imbalance here N σ eðoÞ denotes the total number of spinσ atoms on even (odd) sites and N σ ¼ N σ e þ N σ o . A non-zero steady-state imbalance signals a memory of the initial state, where I σ ðt ¼ 0Þ ¼ 1.
In a first set of measurements we study the effect of interactions on the coherent short-time dynamics. In a tilted lattice an initially localized particle exhibits Bloch oscillations 39 , with a characteristic period T σ = h/Δ σ , set by the spin-dependent tilt. In the presence of interactions, Bloch oscillations persist, showing a rich variety of dynamics, such as interaction-induced dephasing and amplitude modulation [40][41][42][43][44][45] . Here, we use the spin-resolved imbalance to probe real-space Bloch oscillations in a parityprojected manner. In the non-interacting limit the timedependence can be computed analytically: which enables a precise calibration of the model parameters Δ σ and J (Fig. 2a) at short times. Here, J 0 denotes the 0th-order Bessel function of the first kind. The dephasing of the oscillations is caused by a residual harmonic confinement that results in a weak local variation δT σ of the Bloch oscillation period T σ between adjacent sites. An upper bound for the trap frequency ω h /(2π) = 39 Hz was extracted from independent measurements (Supplementary Note 11) and corresponds to δT σ /T σ ≪ 10 −3 .
Since the imbalance dynamics for both spin components is very similar (see Supplementary Fig. 10), we focus on one component I # .
For weak tilt values, Δ ↓ = 1.2J, we find that the dynamics of the interacting spin-mixture (U = 3J) exhibits the same dominant frequency components as the non-interacting Bloch oscillations, while the dephasing is strongly enhanced. This can be seen more directly by calculating the power spectral density (PSD) of the imbalance jĨ σ ðνÞj 2 (inset of Fig. 2a). We find three distinct peaks in the spectrum, the Bloch frequency Δ ↓ and an admixture of two higher harmonics with the largest spectral weight in the second harmonic at ν 1 = 2Δ ↓ /h. For U = 3J its weight is decreased by 70% compared to the non-interacting case. The higher-order harmonics originate from the real-space evolution within one Bloch cycle and are determined by the Bloch oscillation amplitude A σ /d s = 4J/Δ σ . We anticipate frequency components at integer multiples of Δ σ , with an upper bound determined by A σ /d s , in agreement with our data. Interaction effects are expected to be less relevant once the Bloch oscillation amplitude is smaller than one site, resulting in negligible overlap between neighboring particles for our CDW initial state. In Fig. 2b we show the PSD of the coherent shorttime dynamics for Δ ↓ = 3.0J. While the largest spectral weight of the PSD is now contained in the Bloch frequency ν 2 = Δ ↓ /h, the reduction is still about 50% compared to the non-interacting case. Indeed, the spectral weight is a sensitive measure of the interaction-induced dephasing. Moreover, the on-site interactions lift the degeneracy of the energy levels in the Wannier-Stark spectrum, which results in additional frequency components in the PSD. For our parameters (Fig. 2b) they occur at ≈ ν 2 ± 0.5Δ ↓ /h in the time-evolving block decimation (TEBD) simulations [46][47][48] , which is consistent with our data.
The sensitivity of the coherent short-time dynamics on the interaction strength is further highlighted by the strong interaction-dependence of the peak power spectral density (PPSD) jĨ ðν j Þj 2 of the respective dominant frequency components ν j , j = {1, 2} (Fig. 2c, d). We find a sharp decrease of the PPSD by about 40% already for small interaction strength U =± 0.5J for Δ σ = 1.2J. After reaching a global minimum at intermediate interaction strength, it slowly recovers to the noninteracting value in the limit of large interactions.
For long enough evolution times, the coherent Bloch oscillations are dephased and a finite steady-state imbalance develops in the non-interacting limit (Fig. 3a). Note that, if the dephasing was solely due to residual harmonic confinement, we would expect a coherent revival of the oscillations, which is suppressed in our experiment by additional dephasing mechanisms and ensemble averaging. The observed finite steady-state imbalance is caused by Wannier-Stark localization and can be computed analytically by time averaging the short-time dynamics: Excellent agreement between our data and the analytical result provides strong evidence that the effect of the harmonic confinement is negligible for the late-time steady-state imbalance, in contrast to previous fermionic transport experiments 35,36 . This is further supported by the data in Fig. 3b, where the steady-state value is probed for a larger range of tilt values, even reproducing the non-monotonous behavior that is found for small values of the tilt. Note, that the vanishing imbalance, as observed for Δ ↓ ≈ 1.5J (dashed line in Fig. 3b), does not indicate delocalization. It results from localized Wannier-Stark orbitals with equal weight on even and odd sites.
In the presence of weak interactions localization was predicted to survive in the limit of small additional disorder or harmonic confinement, signaled by a finite steady-state imbalance 30,31 .
Here, we find that after a small decay at intermediate times a plateau of the imbalance develops, which persists for long evolution times up to 700τ (Fig. 3a) in the strongly-interacting regime. A comparison with ED simulations (inset Fig. 3a) in a clean system without spin-dependent tilt and without harmonic confinement for a Néel-ordered initial CDW (as opposed to the random-spin initial state realized in the experiment) further highlights that this non-ergodic behavior is not due to experimental imperfections at least for the experimentally relevant observation times (see Supplementary Fig. 5 for a systematic finite-size scaling analysis). Moreover, this robust steady-state value survives over a wide range of parameters (Fig. 3b). As a function of the tilt it qualitatively follows the behavior of the noninteracting system, but shows consistently lower steady-state values.
The persistence of non-ergodicity down to very small values of the tilt is surprising at first sight. One may expect that for large Bloch-oscillation amplitudes the interactions between particles result in a dephasing of the coherent dynamics that give rise to Wannier-Stark localization in the non-interacting limit and hence cause ergodic behavior 35,36,[40][41][42][43] . We study the plateau value for Δ ↓ = 1.1J and find that it is largely independent of interactions (Fig. 3c). In a numerical analysis of this regime for a Néel-ordered singlon CDW we indeed find that the imbalance decays to zero for evolution times on the order of 10 4 τ ( Supplementary Fig. 5), which further agrees with the finite imbalance measured at~200 τ. The observed inversion of the spin-resolved imbalance I # < I " after long evolution times (although Δ ↓ > Δ ↑ ) is explained by the non-monotonic dependence of the stationary imbalance on the tilt for Δ σ < 2J as shown in Fig. 3b.
For intermediate values of the tilt Δ/J ≃ 3 on the other hand, we find a surprisingly robust steady-state imbalance, in agreement with numerical calculations, with a clear interaction dependence (Fig. 3d). The behavior is similar for both spin components and well reproduced by numerical simulations. The deviation between experiment and numerical simulations at larger interaction strengths is most likely due to the finite coupling between 1D chains, which plays a larger role for increased interactions 49 . The steady-state imbalance is symmetric around U = 0 due to a dynamical symmetry [for (Δ ↓ − Δ ↑ ) ≪ J] between attractive and repulsive interactions (Supplementary Note 2), similar to the  homogeneous Fermi-Hubbard model 50,51 . The curve displays a global minimum for intermediate interactions, which we identify with resonant processes at |U| ≃ 2Δ, where two singlons separated by two lattice sites form a doublon. This coincides with regime ① in Fig. 1c, where the largest connectivities were found. The precise value of the resonance is slightly shifted, U res ' 2Δ À 8J 2 =ð3ΔÞ, due to perturbative corrections for finite J/Δ, in agreement with our data (dashed line in Fig. 3e). For large interactions and weak spin-dependence (Δ ↓ − Δ ↑ ) ≪ J, we expect the system to recover the non-interacting regime (Supplementary Note 3).
In order to gain additional insights into the observed non-ergodic behavior, we study the properties of our model perturbatively in the large tilt limit for the two distinct regimes ① and ② (Fig. 1c). In regime ②, Δ ≫ J, |U|, an effective Hamiltonian can be derived in powers of λ = J/Δ. As predicted 19,20,22,31,33 , we find an emergent dipole-conserving HamiltonianĤ  19,20 , seemingly consistent with the observed nonergodic behavior. Yet, higher-order processes Oðλ 4 Þ, relevant for Δ ≃ 3J, are expected to melt the CDW within the experimentally studied timescales 19 . These higher-order processes as well as the dominant off-diagonal contribution, however, require the production of doublons, which is penalized by the on-site interaction U. We numerically show that this leads to a significant slowdown of the dynamics (Supplementary Note 5), which explains the robustness of the steady-state value observed in the experiment. Thus, for large values of the tilt, the doublon number is effectively conserved as well, as suggested in Ref. 31 .
On resonance, |U| ≃ 2Δ (regime ① in Fig. 1c), doublons can be formed without energy penalties, possibly leading to faster dynamics. Indeed, after an initial faster dynamics, we find a lower steady-state imbalance, which cannot be solely explained by the second-order resonant tunneling process shown in Fig. 1c, because it leaves the imbalance invariant. In this regime, we derive an effective HamiltonianĤ res eff [Supplementary Eq. (19)] up to second order in λ (the third order vanishes), conserving the dipole moment, the doublon number or the sum of the two (∑ i;σ in i;σ þ 2∑ ini;"ni;# ). The corresponding symmetry sector exhibits strong fragmentation and results in a finite steady-state imbalance (Supplementary Note 4). In Fig. 4c we show the dominant second-order tunneling terms for our initial state, illustrating the importance of doublon-assisted tunneling processes for the reduction of the steady-state imbalance. For finite λ or longer evolution times, higher-order hopping processes Oðλ 4 Þ enable additional dynamics. These processes are expected to eventually melt the CDW completely, although the required timescales may be very large. In the experiment, we find robust steady-state values even for rather low values of the tilt (Δ ≃ 3J) up to evolution times of about 700τ (Fig. 3a).
In order to connect the large-tilt limit described byĤ res eff to the experimental parameter regime, we investigate the states within the explored subspace N ϵ , which we denote numerical fragment in analogy to the phenomenon of Hilbert-space fragmentation. For simplicity, we study a clean system (ω h = 0, Δ ↑ = Δ ↓ ≡ Δ) and a Néel-ordered CDW initial state. In Fig. 4a, we show the density of states in the Hilbert space H and compare it to the density of states in the numerical fragment N ϵ for different values of the cut-off ϵ. Centered around the energy of the initial state, the density of states acquires a finite width within the numerical fragments, that is approximately set by the many-body bandwidth ± 2JN (dashed line in Fig. 4a), where N = N ↑ + N ↓ denotes the total number of atoms. In stark contrast to thermal systems, the low finite-time connectivity indicates that only a small number of states is relevant for the dynamics. Moreover, it vanishes exponentially in the thermodynamic limit for finite evolution times up to 1000 τ (Fig. 4b). Since the perturbative HamiltonianĤ res eff is only valid in the limit of large tilts, the intersection between the numerically constructed fragment and the analytical one K res ("Methods"), which was derived using the perturbative HamiltonianĤ res eff up to third order in λ = J/Δ, is small for our experimental parameters Δ = 3J and U = 5J (Fig. 4c). We expect, however, that the two subsectors coincide for λ → 0. Indeed the normalized intersection saturates to one, although only for Δ/J ≫ 20. For this comparison the cut-off value ϵðK res Þ is chosen such that dim ðN ϵðK res Þ Þ ¼ dim ðK res Þ, since generally, N ϵ contains a much larger number of states. Despite the large value of λ realized in the experiment, we find strong evidence that the slow dynamics is due to kinetic constraints and that the energetically allowed microscopic processes give rise to the phenomenon of Hilbert-space fragmentation in the large tilt limit, as demonstrated for the two regimes (① and ②). This is further supported by the resonance feature that is shown in the inset of Fig. 4c for the resonant regime ①.

Discussion
In conclusion, we have demonstrated both experimentally and numerically non-ergodic behavior in the tilted 1D Fermi-Hubbard model over a wide range of parameters and have provided a microscopic understanding based on perturbative analytical calculations. For future studies it would be interesting to study the limit of large tilts, where strongly-fragmented effective Hamiltonians were identified and to investigate the initial-state depen-  19,20,22 . Although experimentally challenging due to finite evolution times, it would be interesting to reconcile the phenomenon of Stark MBL and Hilbert-space fragmentation, by studying the impact of weak disorder or residual harmonic confinement on the long-time dynamics 33 . Adding periodic modulation as an additional ingredient, other strongly-fragmented models, scarred models and time crystals could be engineered [52][53][54] or drive-induced localization could be investigated 55,56 . By tuning the direction of the tilt in a 2D lattice, dipole-and higher-moment conserving models could be realized 20,57 enabling studies beyond the hydrodynamic regime 58 . Moreover, it will be interesting to explore the connection between lattice gauge theories and the phenomenon of Hilbert-space fragmentation 21,27,29,[59][60][61] , which could be addressed experimentally in a similar model 62 .

Methods
Experimental sequence. Our sequence begins with loading a degenerate Fermi gas with temperature T/T F = 0.15 (1), where T F is the Fermi temperature, into a threedimensional (3D) optical lattice. The wavelength is λ l = 1064 nm along the x direction and λ ⊥ = 738 nm in the transverse directions. Repulsive interactions during loading in combination with a short, off-resonant light pulse after loading ensure an initial state free of double occupancies (Supplementary Note 7). By adding a short lattice with wavelength λ s = λ l /2 along the x direction, we generate a CDW initial state consisting of singlons (Supplementary Note 7). Holding the gas in this deep 3D lattice with a tilted, bichromatic superlattice along the x direction, dephases remaining correlations between neighboring sites and suppresses any residual dynamics, while ramping up a magnetic field gradient and adjusting the interaction strength. The lattice depths are 18 E rs for the short lattice, 20 E rl for the long lattice and 55 E r⊥ for the transverse lattices. The depths are given in the respective recoil energies, E rj ¼ _ 2 k 2 j =ð2mÞ, with j ∈ {l, s, ⊥}, k j = 2π/λ j the corresponding wave vector, m the mass of 40 K and and ℏ = h/(2π) the reduced Planck constant. The deep transverse lattices decouple the 1D chains aligned along x and generate a 2D array of nearly independent 1D systems. The residual coupling along the transverse directions is typically less then 0.03 % of the coupling J along x. The dynamics to probe the tilted 1D Fermi-Hubbard model described by the Hamiltonian in Eq. (1) is initiated by suddenly switching off the long lattice and quenching the short lattice to depths between 6 E rs and 8 E rs . Simultaneously, the strength of the dipole trap is adjusted in order to compensate the anti-confining harmonic potential introduced by the lattice (Supplementary Note 8). After a variable evolution time t the on-site population is frozen by suddenly ramping up the longitudinal lattices to 18 E rs and 20 E rl respectively. Subsequently, we extract the spin-resolved imbalance I σ , by using a bandmapping technique 63,64 in conjunction with Stern-Gerlach resolved absorption imaging. Note, that the imbalance is defined as a charge imbalance between even and odd lattice sites in our system and does not probe spin imbalances. The magnetization of the systems is conserved during the evolution and it is equal to zero at all times.
Initial state. The initial state in all experiments consists of a CDW of singlons, where " and # states are randomly distributed on even lattice sites and odd lattice sites are empty. We work with an equal mixture of both states (N ↑ = N ↓ ) such that the total magnetization is zero. The fraction of residual holes on even lattice sites, due to imperfections in the loading sequence and due to removed doublons is expected to be about 10% 65 . Excellent agreement between the data and numerical simulations, which do not consider residual holes on even sites, indicates, that the hole fraction has a negligible effect on our dynamics. The initial state can be modelled as incoherent mixture within the zero magnetization sector with density matrixρ ¼ 1 N ∑ fσgj∑ i σ i ¼0 ψ 0 ðfσgÞ ψ 0 ðfσgÞ , where each product state ψ 0 ðfσgÞ , is given by a CDW of singlons and where the sum runs over all N possible permutations of spin configurations {σ}. The product state ψ 0 ðfσgÞ is defined as ψ 0 ðfσgÞ ¼ Q i¼even2trapĉ y i" n i"ĉ y i# n i# 0 j i, whereĉ y iσ is the fermionic creation operator, n iσ ∈ {0, 1}, σ ∈ {↑, ↓}, n i = n i↑ + n i↓ ≤ 1 and i is the lattice-site index along x.
Note, that in the clean translationally-invariant Fermi-Hubbard model the initial CDW corresponds to an infinite temperature state and a finite imbalance value is a hallmark signature of localization. In the tilted model, the spectrum is superextensive complicating a meaningful definition of temperature. This is overcome by transforming into the interaction picture with respect to the tilt potential, which leaves all density observables invariant and allows us to establish the imbalance as a good probe for ergodicity breaking (Supplementary Note 1) Details of numerical calculations. The numerical computations that are compared with the experiment in Figs. 2 and 3 of the main text were performed using ED or TEBD. The parameters J, Δ ↑ and Δ ↓ used in the computations were obtained as fit parameters from the corresponding non-interacting data. Additionally, the effect of harmonic confinement present in the experiment was simulated by scaling the trap frequency by a factor ffiffiffiffiffiffi L exp L q where L exp ¼ 290 is the system size in the experiment and L is the system size used in the numerical calculation. This is done to appropriately simulate the collapse and revival dynamics in the Bloch oscillations induced by the harmonic confinement (Supplementary Note 11).
We use TEBD for short-time dynamics (Fig. 2 of the main text) and ED for long-time dynamics (Fig. 3 of the main text). In ED, we consider the Hilbert space as a tensor product H " H # where H σ is the Hilbert space of spin-σ atoms. In order to efficiently compute the time dynamics, we decompose each time step in the dynamics into three unitary propagators. One each corresponding to the hopping of the two spin components and the third one corresponding to the onsite potential and interactions. We use a Trotter-Suzuki approximation in this decomposition (see Supplementary Note 13 for details and error analysis). In Fig. 3a, d, we use L = 16, N ↑ = N ↓ = 4. In order to effectively model a mixed CDW initial state, in Fig. 3a, this computation is averaged over 20 randomly chosen pure CDW states. In Fig. 3d we use a superposition of pure CDW product states as we are concerned only with time-averaged steady-state value. The parameters J, Δ σ and the harmonic confinement are fixed by fitting to the corresponding noninteracting data.
In Fig. 2, we use TEBD calculations with L = 100 and bond-dimension χ = 120. The truncation error was less than 10 −2 . In Fig. 2b, c, we compare the experimental and numerical data in Fourier space. If the two data sets have different number of samplings in the time domain, we scale the numerical data appropriately after the fast Fourier transform.
Construction of the Krylov subspace. The Krylov subspace (corresponding to the fragment K res ) is constructed by using the effective Hamiltonian on resonancê H res eff in Supplementary Eq. (19). This Hamiltonian is then interpreted as an adjacency matrix in the Wannier basis and the Krylov subspace consists of all states, which are connected to the Néel-ordered CDW initial state. The Krylov subspace K res is closed under time-evolution generated by the effective Hamil-tonianĤ res eff . Starting from initial states within the Krylov subspace K res and including higher-order terms Oðλ 4 Þ, the dynamics is captured only approximately (Supplementary Note 4). An improvement is obtained by further rotating the diagonal basis in which the effective Hamiltonian becomes fragmented with the unitary transformation obtained in powers of λ (as given by the Schrieffer-Wolff perturbative expansion (Supplementary Note 3)). This results in a rotated Krylov subspace.
Construction of the numerical fragment. We define the numerical fragment N ϵ as the span of a subset B ϵ of the number basis B of H, where H is restricted to quarter filling and zero magnetization. We define the set B ϵ via its complement, B ϵ ¼ BnB c ϵ , where B c ϵ would be ideally defined as the largest subset of B satisfying max t < T N ∑ n c 2B c ϵ n c jψðtÞ 2 < ϵ. Here T N defines a time window for the evolution of the initial state ψðt ¼ 0Þ . Equivalently, one could define the subset B ϵ as the smallest one, satisfying min t<T N ∑ n2B ϵ njψðtÞ 2 ≥ 1 À ϵ. We work with the complement, because it is easier to implement numerically. This inequality condition for the complement would ensure that the residual overlap of ψðtÞ outside of N ϵ at any time t ≤ T N is bounded by ϵ. Constructing this B c ϵ , however, involves a search in the powerset of B, which is exponential in the dimension of H. This is intractable even for relatively small system sizes such as L = 7. It follows from the inequality max t < T N ∑ n c n c jψðtÞ 2 ≤ ∑ n c max t < T N n c jψðtÞ 2 that keeping the latter sum smaller than ϵ will ensure that the former sum is also bounded by ϵ. Moreover, the latter sum is computationally easier to handle and therefore, we use it to define the fragment. We construct the numerical fragment N ϵ using a B c ϵ , defined such that ∑ n c max t < T N n c jψðtÞ 2 < ϵ. The gap in the inequality max t < T N ∑ n c n c jψðtÞ 2 ≤ ∑ n c max t < T N n c jψðtÞ 2 loosely depends on the sum ∑ n2B max t < T N njψðtÞ 2 , which is in general, not normalized. Although this sum can be as large as the dimension of H, in the examples that we study, it remains small, i.e., < 10 for L < 20, and grows logarithmically in the dimension of H.

Data availability
All data files are available from the corresponding author upon reasonable request. Source data are provided with this paper.

Code availability
The code that supports the plots within this paper are available from the corresponding author upon reasonable request.