Many-body quantum dynamics of initially trapped systems due to a Stark potential --- thermalization vs. Bloch oscillations

We analyze the dynamics of an initially trapped cloud of interacting quantum particles on a lattice under a linear (Stark) potential. We reveal a dichotomy: initially trapped interacting systems possess features typical of both many-body-localized and self-thermalizing systems. We consider both fermions ($t$-$V$ model) and bosons (Bose-Hubbard model). For the zero and infinite interaction limits, both systems are integrable: we provide analytic solutions in terms of the moments of the initial cloud shape, and clarify how the recurrent dynamics (many-body Bloch oscillations) depends on the initial state. Away from the integrable points, we identify and explain the time scale at which Bloch oscillations decohere.


Introduction-
The historical focus of many-body quantum physics has been the low-energy parts of the many-body spectrum. In recent years, the perspective has changed, largely due to experiments with cold atoms [1][2][3][4], which have inspired the study of non-equilibrium situations in isolated quantum systems [5][6][7][8]. In an isolated situation, energy conservation ensures that a system with an initially high energy will never explore the low-energy parts of the spectrum. The quantum dynamics of isolated systems poses new challenging questions, such as whether observables thermalize [5][6][7]9].
In this work, we address the real-time dynamics of an initially trapped interacting lattice system subject to a linear potential. We present a comprehensive study for two representative systems (featuring bosons and fermions), for all interaction regimes. At zero or infinite interaction, each model becomes integrable (can be mapped to free particles). For intermediate interactions, we have an example of many-body localization without disorder [78][79][80][81][82][83][84], where nevertheless a version of thermalization is valid when we focus on the part of the Hilbert space spanned by states in which particles are confined within a connected spatial region, i.e., the subspace explored by initially trapped systems. We show that the dynamics within such a subspace is thermalizing.
At the "free" points, there is perfectly periodic behavior. We provide a series of exact analytical results for the cloud dynamics in these cases. For strongly interacting (hardcore) bosons, we show dynamical generation (and periodic disappearance!) of fragmented condensation of an initial un-condensed cloud. At strong (weak) initial trapping, the dynamics consists primarily of width (position) oscillations. At intermediate trapping, the skewness undergoes unusual dynamics during every period, of which we do not know of an analog in the literature. Near the integrable points, we show and explain beating behavior of the cloud dynamics, with linear dependences on the integrability-breaking parameter. This explains and generalizes the experimental observation of [36].
Models-We consider N p particles on an infinite lattice subjected to a tilt potential. The total Hamiltonian consists of a kinetic term T = −J/2 j a † j a j+1 + h.c. , with a † j the creation operator of a particle in site j and J the hopping amplitude; a potential term E = E j j a † j a j due to a constant tilt strength E; and an interaction term V. We consider two families of models: the Bose-Hubbard model (BHM) We will mostly take the initial state |Φ 0 to be the ground state of the non-tilted system in the presence of a harmonic potential (H 0 = T + V + W j j 2 a † j a j ), parametrized by the dimensionless constantρ = N p W/J [85,86]. The initial condition can be varied form an extended Gaussian-like cloud (small ρ) to a highly packed state at largeρ. We also consider initial states which are product states, e.g., of the form . For bosons at U = ∞ and for fermions at all V = ∞, the ground state has this form at largeρ.
In addition to U, V = 0, in both strong interacting regimes (U, V → ∞), the dynamics is that of a set of non-interacting particles. For U → ∞ (BHM), double occupancy is kinematically forbidden and the finite energy Hilbert space reduces to that of hard core bosons. In this limit the BHM maps to the Ft-VM with V = 0 via a Jordan-Wigner (JW) transformation. The spectrum of the Ft-VM with V → ∞ and L sites can also be shown to map onto that of a Ft-VM with L − N p sites and V = 0 [87]. In all these (effectively) non-interacting cases, the spectrum of the tilted Hamiltonian consists of equally spaced highly degenerate levels, with spacing E. This yields periodic evolution, with period T = 2π/E, for any initial state. In fact, exact analytical solutions can be found for the many-body evolution [87]. Away from these 'free' cases the dynamics is non-integrable: either because the non-tilted model is already so (BHM); or because a finite tilt breaks the integrability present in the E = 0 case (Ft-VM).
Long time behavior and thermalization-Figs. 1(a,a ) show the eigenenergies ε α of H, corresponding to eigenvectors |α , as a function of the interaction strength, color-coded with |c α | 2 , with c α = α |Φ 0 the overlap amplitude with the initial state. Only some eigenstates have a non-negligible overlap with the initial state; the other eigenenergies are not visible. For fixed N p , Increasing the chain length L (with fixed N p ) increases the Hilbert space dimension polynomially, rendering the spectrum dense at L → ∞, but leaves Figs. 1(a,a ) invariant. Density profiles of the many-body eigenstates which have non-negligible |c α | 2 are exponentially localized within a length proportional to 1/E. Therefore the dynamics of an initially confined cloud of atoms is always localized. This can be traced to the fact that a cloud of atoms in an infinite system is always in the dilute density regime; as interactions are short-range, if the cloud (a, a ) Spectra of the Ft-VM and of the BHM as a function of the interaction strength. The color coding corresponds to the overlap-squared, |cα| 2 . of the initial state, which is |.. expands too much the particles cease to interact with each-other. The exponential localization of the manybody eigenstates, and consequently of the dynamics, is thus ensured by the exponential localization of the single particle eigenstates [10,11].
The effective dimensionality of the Hilbert space spanned by the initial state is shown in Figs. 1(b,b ). This quantity is larger for intermediate interactions than near the 'free' points (small or large U , V ). d eff decreases algebraically with E and increases algebraically with the number of particles N p . We now analyze the long time asymptotic behavior of the cloud dynamics in light of these spectral properties. We define the time averaged densityn j = lim T →∞ T −1 T 0 dt n j (t), with n j (t) = a † j (t)a j (t) the site occupancy, andσ 2 j = lim T →∞ T −1 T 0 dt[n j (t)−n j ] 2 , which quantifies the temporal deviations around the average. For a system with a non-degenerate spectrum these quantities are given by their diagonal ensemble [9] valuesn j = α |c α | 2 α| n i |α andσ 2 j = 0.01 1 100 U/J For systems fulfilling the so called eigenstate thermalisation hypothesis (ETH) [9,89,90] the temporal fluctuations of local observables are strongly suppressed, decreasing exponentially with system size. In contrast, for integrable models, ETH does not hold: The decrease is merely polynomial. In the present case the system does not fulfill ETH trivially -there are an infinite number of eigenstates with the same energy but a vanishing overlap with the initial state. Moreover, as all eigenstates are localized throughout the spectrum, the system behaves as a many-body-localized (MBL) one.
Nonetheless, away from the 'free' points, equilibration may still arise for a sufficiently large N p , i.e. large d eff , in the sense that different trapped initial states with roughly the same energy yield the samen i profile and that long time deviations from the average are suppressed σ ∝ 1/d eff . A comparison between Figs. 1(b, b ) and (c, c ) shows that d −1 eff andσ are qualitatively similar and that the values ofσ substantially decrease with the number of particles in the cloud. This supports an equilibration scenario for both fermionic and bosonic systems away from U, V = 0 and U, V = ∞. At these special values the system becomes integrable and the limits U, V → 0, ∞ and t → ∞ do not commute. At these pointsσ is much larger and decreases much slower with particle number.
Dynamics near 'free' points- Fig. 2 shows some BHM time evolutions at finite interaction values near the 'free' points U = 0, ∞.
The center of mass x t = ( j n j (t) j)/N p and the width σ t = of the cloud both generically show a "collapse and revival" or beating behavior. Other cloud characteristics (skewness or kurtosis) show the same effect [87]. To what extent the phenomenon is visible varies with the initial state and the quantity observed, but generically for U/J not too close to 1, a beat is visible. The beat period is seen to have clear linear dependences, ∝ U −1 at small U/J and ∝ U at large U/J, on the interaction. The behavior at small U/J has recently been observed experimentally [36]. We have found the same behavior in the fermionic case as a function of V [87].
This remarkably simple dependence can be explained using the many-body spectrum. At the free points, this spectrum is exactly equally spaced (steps of E) and highly degenerate. As one moves away from these simple points, the degeneracy is lifted, so that the frequencies available for the dynamics are a range of values around E, the range being small compared to E. This explains the beat behavior. A perturbative argument yields an energy level splitting of the order of V ν or U ν with ν = ±1 for weak/strong interactions. The splitting scale provides the beat frequency.
Spectral considerations also explain why there is rapid relaxation behavior without beats in the U, V ∼ J regime. In this regime, the eigenstates mix, destroying the ladder structure, and the chaotic structure of the spectrum leads to relaxation, as we have analyzed above. The present study in terms of the spectrum thus explains the results of the experiments of Ref. [36].
Cloud dynamics at 'free' points-In contrast with the equilibration seen for moderate interactions, at U, V = 0, ∞ there are perfectly periodic oscillations. The long term state is not equilibrated and has strong dependence on the initial condition. Fig. 3 shows time evolution for the JW-related cases V = 0 and U → ∞, respectively labeled by F or B. The cases of an initially spread-out and narrow cloud (small and largeρ) are shown (top and bottom). The density plots show the evolution of the density n j (t) (identical for F and B), and of the momentum occupation numberñ F (k, t) andñ B (k, t). For the bosonic system we also compute the occupation numbers of the natural orbitals λ n (t) [91,92] (with λ 0 ≥ λ 1 ≥ ...), defined as the eigenvalues of the single particle density matrix A macroscopic occupation (i.e., a λ i of order N p ) corresponds to quasi-condensation.
The density profile n(x, t) displays qualitatively different dynamics for small and largeρ: Bloch oscillations consist of mainly position oscillations forρ 1 and mainly width oscillations forρ 1. For largeρ the shape of the initially localized cloud changes considerably within a period, the shape becoming double-peaked when the cloud widens. The oscillation amplitude of the center of mass x t is large forρ 1 and small forρ 1. The cloud width σ t shows the opposite behavior. (Fig. 3 right.) This distinction is analogous to that observed in single-particle Bloch oscillations [12]. Additional shape dynamics appear at intermediateρ -the cloud becomes strongly skewed once every period [87]. The amplitude of skewness oscillations is non-monotonic as a function ofρ, unlike amplitudes of position (width) oscillations which decreases (increases) monotonically withρ [87].
For any Gaussian initial state, the subsequent cloud dynamics (time evolution of moments) can be obtained analytically as a function of the initial moments of correlators [87]. The center of mass has purely si-   nusoidal oscillations, . (The behavior of µ 1 and µ 2 as functions of N p andρ is described in [87].) This allows to compute the amplitudes of oscillation of the moments, e.g, ∆x = max t x t − min t x t , and ∆σ 2 = max t σ 2 t − min t σ 2 t as a function ofρ for different N p . The position oscillation amplitude is ∆x = 2J/E forρ → 0, and at largeρ decreases as ∆x ∝ J/ (Eρ) for N p > 1 [87]. Conversely, ∆σ 2 increases from zero to 2(J/E) 2 asρ is increased (N p > 1) [87].
The momentum distribution of the fermionic system (F) has simple time evolution: n F (k, t) = n F (k − 2π T t, 0), reminiscent of single-particle Bloch oscillations [87]. For B, the momentum distribution n B (k, t) has similar behavior for smallρ, but it is now a sharply peaked distribution that traverses the Brillouin zone periodically, signaling quasi-condensation in the initial state that survives during the oscillations. The natural orbital occupancy accordingly shows a dominant eigenvalue that stays dominant throughout the evolution. The largeρ behavior is more intricate. Although the condensate is initially non-condensed, two well-defined coherence peaks appear. Remarkably, they disappear periodically for a short fraction of the period when returning to the initial state. The set {λ i } now has two dominant occupancies, λ 0 λ n>0 , signaling a fragmented condensate that is dynamically generated [70,93] and persists for almost all times within each period.
Discussion-We have presented a thorough study of many-body Bloch dynamics in two standard lattice models in one dimension, one fermionic and one bosonic. A main result is that generic many-body systems under a tilt potential have a dichotomic nature, possessing both ETH and MBL features. Although their eigenstates are exponentially localized, and an initially trapped cloud has finite overlap only with a zero-measure set of eigenstates within the relevant energy window, the long time dynamics yield a thermalized state within a Hilbert space of effective dimension d eff which increases with the number of particles N p .
The approach to the thermalized state can be seen as the destruction of the many-body Bloch oscillations which are present at the integrable ('free') limits, both for weak and strong coupling. We show that the relevant time-scale grows as U (U −1 ) or V (V −1 ) away from the weak (strong) integrable limit. At the free limits we present several striking features of the cloud dynamics, including a dynamical generation (and periodic disappearance) of fragmented condensation for strong initial trapping.

Supplemental Materials for:
Many-body quantum dynamics of initially trapped systems due to a Stark potentialthermalization vs. Bloch oscillations

S.I. CONTENTS
In these Supplemental Materials, • We describe the four non-interacting ('free') points in the two one-dimensional models we have considered (Section S.II). These are the U = 0, ∞ points for the Bose-Hubbard model (BHM) and the V = 0, ∞ points for the fermionic t-V model (Ft-VM).
• We provide a number of exact results valid at two of the free points, with derivations (Section S.III).
• We discuss the beating behavior and provide further numerical data for the beating dynamics of the interacting systems (Section S.IV).
• We discuss some aspects of the many-body spectrum, including level-spacing statistics (Section S.V). The level-spacing data supports our picture of dual thermalizing and locallizing behavior.

S.II. THE FOUR 'FREE' POINTS
For the bosonic U = 0 and the fermionic V = 0 cases, the Hamiltonian is quadratic, i.e., that of free bosons and free fermions respectively. This implies that the knowledge of all n-point correlators of the initial state allows for an analytical solution of the subsequent evolution of the n-point correlators. If the initial state is Gaussian (so that all n-point correlators are determined by 2-point correlators), then the state remains Gaussian under time evolution. This is particularly straightforward for the V = 0 case where all physically motivated initial states we used are Gaussian, i.e. can be seen as ground states of free particle models.
The non-interacting bosonic model (U = 0) does not admit an initial Gaussian state at zero temperature and finite particle number. Nonetheless, from the knowledge of the correlation matrix (matrix of two-point correlators) in the initial state one can obtain all subsequent two-point correlators. If the two-point correlators of the bosonic U = 0 and of the fermionic V = 0 are the same initially, then in the subsequent dynamics the two-point correlators of the two models continue to be the same.
The dynamics at the bosonic U = ∞ point can be mapped to that of the fermionic V = 0 point by a Jordan-Wigner (JW) transformation that provides a mapping at the operator level. In particular, the density of the fermionic and bosonic systems are the same. (As a side product, this allows us to conclude that the dynamics of the two-point correlators at U = 0 and U → ∞ is also the same if the initial state is the same.) Off-diagonal elements in position basis for the bosonic U = ∞ point and the fermionic V = 0 point are not so simply related, but the hard-core boson correlators can be computed numerically from the V = 0 free fermion dynamics. We have presented such results in the main text.
Finally, although we can prove that the spectral properties of the V = ∞ are that of free particles (next section), no mapping was found at the operator level and thus the dynamics of correlators or densities is not simply related to that of a non-interacting model.
The fermionic t-V model allows for a mapping between the V → ∞ Hamiltonian (in the reduced Hilbert space sector where the energy is finite) and the V = 0 Hamiltonian on a smaller chain.
In the V → ∞ limit, states with particles at nearestneighbour sites are kinematically excluded. For N p fermions in L sites, this means an effective Hilbert space of dimension L − N p + 1 N p , which is the same as the Hilbert space dimension of a system with N p fermions in L − N p + 1 sites without this constraint. We can design a mapping of the Hilbert space for V → ∞ for a chain of size L and N p fermions onto the Hilbert space for V = 0 with N p fermions on L = L−N p +1 sites. The mapping preserves the Hamiltonian. The mapping of the Hilbert spaces is shown in Figure S1 for N p = 2 and 3, by showing the configuration spaces allowed in the two cases. For a 1D system with N p fermionic particles, antisymmetry dictates that the configuration space consists of states with j 1 < j 2 < ... < j Np . Thus, for N p = 2 particles where configurations can be described as (j 1 , j 2 ) pairs, the allowed configurations exclude the diagonal line on the j 1 -j 2 plane. In the V → ∞ case, the states with particles at nearest-neighbour sites j i and j i + 1 are also excluded; hence the next-to-diagonal points on the j 1 -j 2 plane are also excluded. This results in an identical number and topology of allowed configurations in the V = 0 case and V → ∞ cases, when there is one more site in the latter case. This is displayed in the left panels of Figure  S1.
The description is analogous for N p = 3, with the diagonal plane being excluded due to antisymmetry, and one further next-to-diagonal plane being excluded in the case of V → ∞, as shown in the right panels of Figure  S1. The construction is trivially generalized to arbitrary N p < L/2, but difficult to display visually for larger N p . In Table S.1, the mapping is shown to work for a N p = 4 case, by listing all the configurations in the two cases.
The matrix elements of the Hamiltonian between basis states correspond to the hopping of a single particle. Diagonal matrix elements are given by the potential energy of the tilting field. The Hamiltonians are identical in the Mapping between Hilbert spaces of the V = 0 fermionic system and the V → ∞ fermionic system with Np − 1 more sites on the tight-binding chain. Left panels: Np = 2 particles. The allowed configurations of the V = 0 system are shown on the (j1, j2) plane (left top), where j1 and j2 are the positions of the two particles. Lines indicate configurations connected by a single-particle hopping process. Left bottom panel shows that the V → ∞ system with one more lattice site has configuration space of the same size and same hopping topology. Right panels: Np = 3 particles. The corresponding V → ∞ system now has two extra lattice sites compared to the V = 0 system. case of the V → ∞ system and the V = 0 system with the number of sites reduced by N p − 1, except for a possible constant shift on the diagonal terms. (Since the lattice sizes are unequal in the two cases, the definition of zero Stark energy might be chosen independently in the two cases, so this shift is arbitrary, and anyway does not affect the dynamics.) The off-diagonal matrix elements due to single-particle hopping are identical because the mapping preserves the topology of the configuration spaces, i.e., if two configurations of the L-site V → ∞ system are connected by a single-particle hopping process, then the corresponding configurations of the (L − N p + 1)site V = 0 system are also connected by a single-particle hopping process. In Figure S1 for the N p = 2 case, such pairs are joined by lines, and it is visually obvious that the network topology is preserved under the mapping of Hilbert spaces.
The present mapping provides a one-to-one correspondence between basis states and establishes the equality of the Hamiltonians in the two cases. However it does not translate to a simple mapping between creation and annihilation operators, which would have allowed a computation of correlators in the V = ∞ case from a free-   Figure S1, now for Np = 4 particles. The configuration space would be 4-dimensional in the representation used in Figure S1, so we simply list all allowed configurations in the two cases. It is easy to verify that: (1) if a pair of configurations on the left are connected by single-particle hopping, then the corresponding pair on the right is also connected by single-particle hopping; (2) once a zero energy is chosen for the electric field for the two systems, the energy difference between a configuration on the left and the corresponding configuration on the right is the same for all 15 configurations.
particle calculation. We are not aware of a mapping at the operator level that takes us from the V = ∞ model to a non-interacting system.

S.III. EXACT SOLUTION AT NON-INTERACTING POINTS
In the main text, we highlighted some results for the fermionic V = 0 and bosonic U → ∞ systems. We now provide some more details and explicit expressions for time evolution and asymptotic behaviors.

S.III.A. Analytic expressions for moments of the cloud
We present the derivations for free fermions. Since this concerns occupancies, the U → ∞ bosonic system is described by the same equations.
The Hamiltonian can be written as H = T + E with where the operators in real and momentum space obey the usual relations c n = π −π dk 2π e ikn c k , c k = n e −ikn c n . For nearest-neighbor hoppings, ε (k) = −J cos (k), however the following argument holds for a generic dispersion relation. Using the fermionic commutation relations, the evolution operator can be written as Applying the evolution operator in this form to the single-particle density matrix in momentum space yields where ... t denotes the mean value taken at time t. For ε (k) = −J cos (k), Eq.(S.4) simplifies to so that the real-space correlators are found to be using the identity e z cos(θ) = ∞ n=−∞ I n (z) e inθ where I n (z) is the modified Bessel function.
In order to compute the moments of the cloud we define the generalized characteristic function For a = 0, the µ a,m are simply the moments of the cloud shape: x m t = N m p µ 0,m (t). For a = 0, they may be regarded as moments of two-point correlators. In particular, the quantities µ 1 and µ 2 defined in the main text equal to µ 1,0 (t = 0) and µ 2,0 (t = 0) respectively. In case the initial state is the ground-state of an harmonic trap one has µ −a,m (0) = µ a,m (0) ∈ R for m even and µ a,m (0) = 0 for m odd. The first non-trivial generalized moments µ a,m (t) of the trap ground state are shown in Fig.S2 as function ofρ, for different values of N p . Using Eq.(S.6), the characteristic function can be written as Taylor expanding the previous expression and identifying the powers of λ in both sides of Eqs. (S.7) and (S.9) for the case a = 0, one obtains explicit expressions for the first 3 moments of the cloud shape: as a function of the initial values µ a,m (0) of the generalized moments. In the cases of interest here (starting with trap ground states), these initial values are displayed in Fig.S2.

S.III.B. Cloud shape dynamics
Using the exact solution above, one can describe the position and shape oscillations of the cloud during the periodic evolution.
We consider the amplitude of variation of the center of mass and the amplitude of cloud width oscillations Theρ-dependence (for both ∆x(E/J) and ∆σ 2 (E/J) 2 ) are very similar for all N p > 1, converging rapidly to the large-N p limit. The single-particle (N p = 1) behavior differs significantly.
For N p > 2, ∆x(E/J) decreases from 2 to 0, while ∆σ 2 (E/J) 2 increases from 0 to 2. This reflects the physics that the Bloch oscillations are primarily position oscillations for smallρ and primarily width oscillations for largeρ. (For the single-particle case, the large-ρ limit is different.) From the exact solutions, ∆x(E/J) and ∆σ 2 can be expressed in terms of µ 1 = µ 1,0 (0) and µ 2 = µ 2,0 (0). It is easy to see that ∆x(E/J) = 2µ 1 . The expression for ∆σ 2 is more complicated: < 1 and ∆σ 2 E 2 /J 2 = ξ otherwise. Here ξ = 2 µ 2 − 2µ 2 1 + 1 . We now consider the third moment, which gives the skewness of the cloud. In Fig.(S4) we show the skewness computed at t = T /2, where the cloud typically shows a larger deformation with respect to its initial shape. The skewness it seen to have a non-monotonic dependence onρ. s (T /2) vanishes for bothρ = 0 andρ → ∞. Forρ = 0 this is due to an almost undeformed cloud evolution. In theρ → ∞ case, while there is significant deformation, the cloud remains symmetric throughout the whole oscillation period. For a fixed tilt strength E and a large number of particles (see Fig.(S4) lower panel for E = 0.1J and N p = 50) this quantity passes by an N p -dependent minimum. Fig. S4 upper panel shows the the density profile of the atomic cloud for t = 0 and t = T /2 for the points marked (with arrows) in the lower panel. The minimum skewness point corresponds to a highly asymmetric cloud shape having a shock-wave-like form (e.g., theρ = 0.53 panel in Fig. S4). For smallish particle numbers (N p 10), thẽ ρ-dependence is more intricate -there is both a positive maximum and a negative minimum of s(T /2). The shape of the distorted cloud when having a positive maximum is exemplified in theρ = 1.72 panel. For N p 20 there is a unique (negative) minimum that shifts to largerρ with increasing N p .

S.III.C. Recurrent occupancies of the natural orbitals
Here we present some more details of the time evolution of the natural orbital occupancies, λ n . Figure S5 shows the λ n as a function of n, corresponding to the two cases presented in Fig.3 in the main text, through density plots in addition ot snapshots.
For theρ = 0.1 case the lowest natural orbital has a substantial occupation already in the initial state -this is a single-mode quasi-condensate. During the course of the evolution the distribution of the λ n 's does not get substantially modified and the initially quasi-condensed state is observed to remain stable throughout the time evolution.
On the contrary, largeρ induces a Mott insulator state as the initial condition for which the occupation of the natural orbits is given by λ n = 1 for n < N p and λ n = 0 otherwise. In Figure S5 (right panels), we emphasize that this state is realized periodically but only for times in a small vicinity of the multiples of the period t = mT with m ∈ Z. During most of the evolution the two lowest modes get substantially occupied giving rise to a periodically regenerated bimodal quasi-condensate, as described in the main text. Since the evolution is periodic, in the quantities M jj ;nn the integral can be taken over a period of the evolution This derivation is completely general for an initial state correlation matrix c † n c n 0 . For simplicity, let us concentrate on the asymptotic form of the density, for which j = j , for the special case of initial condition with c † n c n 0 = δ nn c † n c n 0 , i.e., product states as initial states. For this particular case the long time averaged density yields In the same way, the fluctuation around the average are given bȳ Note that heren j andσ j were computed performing the integral over time explicitly. Away from the integrable points, the same quantities given in the main text were obtained using the diagonal ensamble expression:n j = α |c α | 2 α| n j |α andσ 2 j = α =α |c α | 2 |c α | 2 | α| n j |α | 2 . Diagonal ensemble values are only valid for a system with non-degenerate energy levels which is the case for finite interactions. Therefore the limits U, V → 0, ∞ of the quantities obtained in the main text do not coincide with those computed here. One way of understanding this is that, in the computation of the asymptotic long time averages, the limits t → ∞ and U, V → 0, ∞ do not commute.

S.IV. BEATING VS EQUILIBRATION BEHAVIORS
In this section we expand on the results reported in the main text concerning many-body Bloch dynamics slightly away from the non-interacting points, i.e., the Ft-VM at V J and V J, and the BHM at U J and U J.

S.IV.A. Spectral explanation of beating frequency
For the BHM, the beating frequency is ∝ U −1 for U J and ∝ U for U J, as shown in Figure 2 of the main text by plotting the beat period against U . For the Ft-VM, the beating frequency will similarly be ∝ V −1 for V J and ∝ V for V J, near the two non-interacting points.
As announced in the main text, a perturbative argument starting from the corresponding 'free' point (U, V = 0, ∞) explains this behavior. This is illustrated in Figure  S7 (and expanded below) for the case of large U . Exactly the same argument holds for small U and for the Ft-VM. One simply has to replace 1/U by the corresponding perturbative parameter (U , V or 1/V ), as appropriate.
The spectrum of a single particle in a Stark ladder has equally spaced non-degenerate levels with spacing E. Note that, each eigenstate corresponds to localization around a particular site. For a non-interacting ) Left: The single-particle spectrum is non-degenerate and equally spaced. Center: The many-body spectrum for the non-interacting case is obtained by filling up the singleparticle levels in all possible ways, hence it is also equally spaced, but each many-body eigenenergy is massively degenerate. Right: moving away from the free point, the degeneracies get lifted, so that each energy level is broadened. In the perturbative regime, the broadening is O(U −1 ). Dynamics of three bosons starting at an initial product-state configuration ...0011100..., with E = 0.35 and U = 20. The dynamics of the cloud is shown through the time evolution of first four moments of the occupancy distribution: the center of mass < x >, the r.m.s. width σ, the skewness s and the kurtosis. Beating behavior with the same beating period is seen in all these quantities. many-body system, the many-body spectrum can be constructed out of the single-particle spectrum by filling the single-particle levels with various numbers of particles. In this case, the many-body eigenenergies are sums of singleparticle eigenenergies. Hence the possible values of the many-body eigen-energies are also equally spaced with spacing E. However, these levels are now highly degenerate, as many different combinations of single-particle eigenstates can lead to the same total eigenenergy. For a fixed number N of particles, the degeneracy of a level far from the spectral edges would scale with the number of sites as O(L N −1 ) (ignoring the fact that, for any finite L, the levels are not exactly equally spaced and so the degeneracies are not exact). In the limit L → ∞ that we are interested in, each level is infinitely degenerate for any N > 1, and the degeneracies are exact.
Since the levels are all equally spaced with spacing E, the time evolution of any observable quantity will be perfectly periodic with period 2π/E. Now we consider perturbing the system by moving away slightly from the 'free' points. The perturbation parameter is 1/U for the BHM system near the hardcore limit. The perturbation will lift the degeneracies, and the splitting is proportional to the 1/U . The time evolution of any observable quantity under this Hamiltonian can be written (when the initial state is expanded in the energy eigenstate basis) as a sum of oscillating terms, with the oscillation frequencies being the energy differences between eigenstates. Now, because of the splitting as shown in Figure S7 (right panel), the frequencies are not all equal to E, rather they are clustered around the value E with frequency difference of order 1/U . This explains why the beat periods are inversely proportional to the perturbation parameter. In the main text, we presented time evolution data displaying beating behavior in the center of mass and width of the cloud, for the BHM, and reported that the same behavior can be seen in other observables, such as the skewness and kurtosis of the cloud. This is shown in Fig. S7. The beating behavior is also visible in the time evolution of the site occupancy or double occupancy (not shown).  FIG. S10. Fermionic t-V model, dynamics of 3 fermions, E = 0.3, V = 0.1. The initial state is changed from left to right. The product initial state (...001010100...) used for the leftmost panels is a proxy for box-like (tightly-trapped) initial state while the rightmost column corresponds to weakest trapping.
In the main text, we also reported that the same beating behavior appeared in the fermionic system. This is shown in Figure S9.

S.IV.C. ∆x and ∆σ as function of trap strength
For the non-interacting system (V = 0 for Ft-VM and U = ∞ for BHM), we have shown that Bloch oscillations are primarily position oscillations (large ∆x, small ∆σ) if the trapping is weak and the initial cloud shape is gaussian-like, while the oscillations are predominantly width oscillations (small ∆x, large ∆σ) if the trapping is strong and the initial cloud is 'box'-shaped.
In Figure S10, we demonstrate that some traces of this phenomenon survive, at least initially, when interactions are added. The example shown is for the fermionic model (Ft-VM), for a small interaction. The dynamics of x and σ are shown for three different initial states. On the left column, the initial state is ...001010100..., which may be considered as the analog of a box-like initial state (ρ = ∞) for fermions. The center and right columns correspond to finite trap ground states, with the rightmost column corresponding to weaker traps.
Going from the leftmost to rightmost columns, the amplitude of center-of-mass oscillations gets larger, while the amplitued of width oscillations becomes smaller. Thus the intuition of which type of oscillation (position vs width) dominates, which we have gained from the noninteracting systems, continues to be valid for interacting systems.

S.V. VIEWING THE MANY-BODY SPECTRUM
We have made several arguments about the many-body dynamics based on the many-body spectrum for fixed N p and infinite L. Unfortunately, it is not possible to numerically calculate the many-body spectrum explicitly S8 for L → ∞. Nevertheless, we can still use the results of finite-size numerical diagonalizations to infer relevant features of the many-body spectrum.
One issue arising with finite-L data is that, even for a single particle, the spectrum is not exactly equally spaced. The deviation is more severe near the edges of the spectrum, i.e., for eigenstates with a significant occupancy near the edges of the finite lattice. For eigenstates localized far from the edges, the corresponding singleparticle eigenvalues are nearly equally spaced.
In the main text, we employed the trick of intensitycoding (color-coding) the many-body eigenstates by the overlap with a many-particle state trapped near the center of the lattice. This ensures that the edges of the lattice plays no role, so that we obtain an L-independent picture.
In Figure S11, we use a complementary procedure to provide another view of the many-body spectrum. The numerically calculated eigenstates are filtered so that any eigenstate with a significant occupancy at one of the edge sites is not shown. A view of the central part of the spectrum is shown for various values of the cutoff. It seems reasonable to presume that, for the infinite chain, the many-body spectrum in any slice of energy is similar to the picture obtained with cutoff 0.1, with a more dense spectrum at intermediate U .
In the bottom panel, we have analyzed the level statistics of the spectrum obtained for L = 21 with cut-off 0.1. For intermediate U , the value approaches that expected for a chaotic system (GOE value). In the infinite-size limit, it is expected that the spectrum is still not chaotic and will be Poissonian at all U due to many-body localization. (Eigenstates localized in different regions can be expected not to interact with each other; hence have no level repulsion.) However, our cutoff procedure is biased toward eigenstates trapped near the center of the lattice, and we have shown that if one restricts to the part of the Hilbert space localized or trapped in a particular region then these systems behave like thermalizing (ETH-obeying or chaotic) systems. This is visible in our analysis of level statistics in the bottom panel of Figure  S11.