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

TL;DR

This paper investigates the complex dynamics of trapped interacting quantum particles under a Stark potential, revealing a mix of localization and thermalization behaviors, with analytic solutions at integrable points and decoherence analysis away from them.

Contribution

It provides new analytic solutions for initial trapped quantum systems at integrable limits and explains the decoherence of Bloch oscillations in non-integrable regimes.

Findings

Analytic solutions for integrable limits of fermions and bosons.

Identification of the decoherence time scale for Bloch oscillations.

Demonstration of mixed localization and thermalization features.

Abstract

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.

Tables

Table 1. Table S.1: Same mapping as in Figure S1 , now for N p = 4 subscript 𝑁 𝑝 4 N_{p}=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.

$N_p=4$, $L=6$, $V=0$	$N_p=4$, $L=9$, $V\to \infty$
0 0 1 1 1 1	0 0 1 0 1 0 1 0 1
0 1 0 1 1 1	0 1 0 0 1 0 1 0 1
0 1 1 0 1 1	0 1 0 1 0 0 1 0 1
0 1 1 1 0 1	0 1 0 1 0 1 0 0 1
0 1 1 1 1 0	0 1 0 1 0 1 0 1 0
1 0 0 1 1 1	1 0 0 0 1 0 1 0 1
1 0 1 0 1 1	1 0 0 1 0 0 1 0 1
1 0 1 1 0 1	1 0 0 1 0 1 0 0 1
1 0 1 1 1 0	1 0 0 1 0 1 0 1 0
1 1 0 0 1 1	1 0 1 0 0 0 1 0 1
1 1 0 1 0 1	1 0 1 0 0 1 0 0 1
1 1 0 1 1 0	1 0 1 0 0 1 0 1 0
1 1 1 0 0 1	1 0 1 0 1 0 0 0 1
1 1 1 0 1 0	1 0 1 0 1 0 0 1 0
1 1 1 1 0 0	1 0 1 0 1 0 1 0 0

Equations

$$H=\mathcal{T}+\mathcal{E}+\mathcal{V},$$

$$\mathcal{T}$$

$$\mathcal{E}$$

$$e^{-iHt}$$

$$\begin{gathered} {⟨c_k^{†}{c}_{k^{\prime }}⟩}_t=e^{-i{\int }_0^{t}d{t}^{\prime }\left[\epsilon \left(k^{\prime }+E{t}^{\prime }\right)-\epsilon \left(k+E{t}^{\prime }\right)\right]} \\ \times {⟨c_{\left(k+Et\right)}^{†}{c}_{\left(k^{\prime }+Et\right)}⟩}_{t=0} \end{gathered}$$

$$\begin{gathered} {⟨c_k^{†}{c}_{k^{\prime }}⟩}_t=e^{2i\frac{J}{E}\sin \left(\frac{Et}{2}\right)\left[\cos \left(k^{\prime }-\frac{tE}{2}\right)-\cos \left(k-\frac{Et}{2}\right)\right]} \\ \times {⟨c_{\left(k+Et\right)}^{†}{c}_{\left(k^{\prime }+Et\right)}⟩}_{t=0} \end{gathered}$$

$$\begin{gathered} {⟨c_x^{†}{c}_{x^{\prime }}⟩}_t=\sum _{y,y^{\prime }}\left\{e^{i\frac{Et}{2}\left(x+y\right)}{I}_{x-y}\left[-2i\frac{J}{E}\sin \left(\frac{Et}{2}\right)\right]\right\}\times \\ {⟨c_y^{†}{c}_{y^{\prime }}⟩}_0\left\{e^{-i\frac{Et}{2}\left(x^{\prime }+y^{\prime }\right)}{I}_{-x^{\prime }+y^{\prime }}\left[2i\frac{J}{E}\sin \left(\frac{Et}{2}\right)\right]\right\} \end{gathered}$$

$$G_a\left(\lambda ,t\right)$$

$${\mu }_{a,m}\left(t\right)$$

$$\begin{gathered} G_a\left(\lambda ,t\right)=e^{iEta}\sum _b{e}^{-i\left(\frac{E}{2}t+\frac{\pi }{2}\right)b} \\ \times I_b\left[4i\frac{J}{E}\sin \left(\frac{Et}{2}\right)\sin \left(\frac{\lambda }{2}\right)\right]G_{a-b}\left(\lambda ,0\right) \end{gathered}$$

$${\mu }_{0,1}\left(t\right)=N_p^{-1}⟨x⟩=-\frac{2}{N_p}\frac{J}{E}{\sin }^2\left(\frac{Et}{2}\right){\mu }_{1,0}\left(0\right),$$

$$\begin{gathered} {\mu }_{0,2}\left(t\right)=N_p^{-2}⟨x^2⟩={\mu }_{0,2}\left(0\right) \\ +\frac{2}{N_p^2}{\left(\frac{J}{E}\right)}^2{\sin }^2\left(\frac{Et}{2}\right)\left[1-\cos \left(Et\right){\mu }_{2,0}\left(0\right)\right], \end{gathered}$$

$$\begin{gathered} {\mu }_{0,3}\left(t\right)=N_p^{-3}⟨x^3⟩=-\frac{6}{N_p}\frac{J}{E}{\sin }^2\left(\frac{Et}{2}\right){\mu }_{1,2}\left(0\right) \\ -\frac{1}{N_p^3}{\left(\frac{J}{E}\right)}^3{\sin }^2\left(\frac{Et}{2}\right)\{\left[6{\sin }^2\left(\frac{Et}{2}\right)+\frac{1}{2}{\left(\frac{E}{J}\right)}^2\right] \\ \times {\mu }_{1,0}\left(0\right)+\left[\cos \left(2Et\right)-\cos \left(Et\right)\right]{\mu }_{3,0}\left(0\right)\}. \end{gathered}$$

$$\Delta x=\underset{t}{\max }{⟨x⟩}_t-\underset{t}{\min }{⟨x⟩}_t$$

$$\Delta {\sigma }^2=\underset{t}{\max }{\sigma }_t^2-\underset{t}{\min }{\sigma }_t^2.$$

$$\begin{gathered} \overline{⟨c_j^{†}{c}_{j^{\prime }}⟩}=\underset{t\to \infty }{\lim }\frac{1}{t}{\int }_0^t{⟨c_j^{†}{c}_{j^{\prime }}⟩}_{t^{\prime }}d{t}^{\prime } \\ =\sum _{n,n^{\prime }}{M}_{j{j}^{\prime };n{n}^{\prime }}{⟨c_n^{†}{c}_{n^{\prime }}⟩}_0. \end{gathered}$$

$$\begin{gathered} M_{j{j}^{\prime };n{n}^{\prime }}=\frac{1}{T}{\int }_0^{T}dt\{e^{i\frac{V_{0}t}{2}\left(j+n\right)} \\ \times I_{j-n}\left[-2i\frac{J}{V_0}\sin \left(\frac{t{V}_0}{2}\right)\right] \\ \times e^{-i\frac{V_{0}t}{2}\left(j^{\prime }+n^{\prime }\right)} \\ \times I_{-j^{\prime }+n^{\prime }}\left[2i\frac{J}{V_0}\sin \left(\frac{t{V}_0}{2}\right)\right]\}. \end{gathered}$$

$${\bar{n}}_j=\overline{⟨c_j^{†}{c}_j⟩}=\sum _a{m}_a{⟨c_{j+a}^{†}{c}_{j+a}⟩}_0$$

$$\begin{gathered} m_a=M_{jj;(j+a)(j+a)} \\ =\frac{1}{T}{\int }_0^{T}dt{\left(-1\right)}^a{\left\{I_a\left[-2i\frac{J}{V_0}\sin \left(\frac{t{V}_0}{2}\right)\right]\right\}}^2 \end{gathered}$$

$$\begin{gathered} {\bar{\sigma }}_j^2=\overline{{\left(⟨c_j^{†}{c}_j⟩-\overline{⟨c_j^{†}{c}_{j^{\prime }}⟩}\right)}^2} \\ =\sum _{a,b}{m}_{a,b}{⟨c_{j+a}^{†}{c}_{j+a}⟩}_0{⟨c_{j+b}^{†}{c}_{j+b}⟩}_0-{\bar{n}}_j^2 \end{gathered}$$

$$\begin{gathered} m_{a,b}=\frac{1}{T}{\int }_0^{T}dt{\left(-1\right)}^{a+b}\times \\ {\left\{I_a\left[-2i\frac{J}{V_0}\sin \left(\frac{t{V}_0}{2}\right)\right]I_b\left[-2i\frac{J}{V_0}\sin \left(\frac{t{V}_0}{2}\right)\right]\right\}}^2. \end{gathered}$$

Full text

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 Bloch et al. (2008); Morsch and Oberthaler (2006); Kinoshita et al. (2006); Greiner et al. (2002), which have inspired the study of non-equilibrium situations in isolated quantum systems D’Alessio et al. (2016); Eisert et al. (2015); Polkovnikov et al. (2011); Dziarmaga (2010). 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 D’Alessio et al. (2016); Eisert et al. (2015); Polkovnikov et al. (2011); Rigol et al. (2008).

A well-known example for which isolation leads to drastically different dynamics is the phenomenon of Bloch oscillations Bloch (1929); Zener (1934). Particles in a tight-binding lattice subject to a linear potential, e.g., due to gravity or an electric field, do not accelerate toward lower potentials, but rather undergo local oscillations. For a single particle, the shape and/or position of the particle wavefunction oscillates, perfectly periodically Bloch (1929); Zener (1934); Hartmann et al. (2004); Thommen et al. (2004a). Long after its prediction, Bloch oscillations were observed in semiconductor super-lattices Waschke et al. (1993); Lyssenko et al. (1997), in cold atoms Ben Dahan et al. (1996); Niu et al. (1996); Raizen et al. (1997); Anderson and Kasevich (1998), and in periodic photonic structures Morandotti et al. (1999); Pertsch et al. (1999); Sapienza et al. (2003). In cold atom experiments, Bloch oscillations have by now been observed many times Ben Dahan et al. (1996); Niu et al. (1996); Raizen et al. (1997); Anderson and Kasevich (1998); Morsch et al. (2001); Cristiani et al. (2002); Roati et al. (2004); Battesti et al. (2004); Ferrari et al. (2006); Morsch and Oberthaler (2006); Cladé et al. (2006); Drenkelforth et al. (2008); Gustavsson et al. (2008); Haller et al. (2010); Poli et al. (2011); Tarruell et al. (2012); Tarallo et al. (2012); Genske et al. (2013); Meinert et al. (2014); Preiss et al. (2015); Keßler et al. (2016); Georges et al. (2017); Geiger et al. (2018), and are used widely as a measurement tool, e.g., for metrological applications Battesti et al. (2004); Cladé et al. (2006); Poli et al. (2011); Tarallo et al. (2012) to detect Dirac points in optical lattices Tarruell et al. (2012), etc. Some experiments have also explored the effect of inter-particle interactions on Bloch oscillations Morsch et al. (2001); Roati et al. (2004); Gustavsson et al. (2008); Drenkelforth et al. (2008); Haller et al. (2010); Meinert et al. (2014). Theoretical treatments of Bloch oscillations have addressed a variety of single-particle situations Kolovsky et al. (2002); Hartmann et al. (2004); Thommen et al. (2004a); Breid et al. (2006, 2007); Thommen et al. (2004b); Walter et al. (2010); Longhi (2010); Collura et al. (2012); Kolovsky and Bulgakov (2013); Khomeriki and Flach (2016); Kartashov et al. (2016); Longhi (2016); Zheng et al. (2017), interacting few-particle systems Dias et al. (2007); Krimer et al. (2009); Dias et al. (2010); Longhi and Della Valle (2012); Wiater et al. (2017), and interacting many-body systems Buchleitner and Kolovsky (2003); Kolovsky (2003); Kolovsky and Korsch (2004); Breid et al. (2007); Schulte et al. (2008); Salerno et al. (2008); Kolovsky et al. (2009, 2010); Mierzejewski and Prelovšek (2010); Walter et al. (2010); Eckstein et al. (2010); Eckstein and Werner (2011); Witthaut et al. (2011); Rubbo et al. (2011); Cai et al. (2011); Mandt et al. (2011); Gaul et al. (2011); Longhi (2012); Carrasquilla et al. (2013); Mandt (2014); Mahmud et al. (2014); Driben et al. (2017). Interactions have been treated both in mean-field (e.g., Gross-Pitaevskii) regimes Kolovsky and Korsch (2004); Breid et al. (2007); Salerno et al. (2008); Schulte et al. (2008); Kolovsky et al. (2009); Walter et al. (2010); Kolovsky et al. (2010); Witthaut et al. (2011); Gaul et al. (2011); Driben et al. (2017) and beyond the mean-field regime Buchleitner and Kolovsky (2003); Mierzejewski and Prelovšek (2010); Eckstein et al. (2010); Eckstein and Werner (2011); Rubbo et al. (2011); Carrasquilla et al. (2013); Mahmud et al. (2014).

Recent experiments Meinert et al. (2014) have found, by tuning bosonic on-site repulsion, the collapse and revival of the oscillation of the cloud position, with the revival period proportional to interaction strength. In addition, sufficiently far from the non-interacting point, the atom cloud was found to have ‘chaotic’ behavior leading to rapid relaxation.

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 Gavish and Castin (2005); Yao et al. (2016); Schiulaz et al. (2015); Antipov et al. (2016); Smith et al. (2017); Mondaini and Cai (2017); Schulz et al. (2018), 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 Meinert et al. (2014).

Models— We consider $N_{p}$ particles on an infinite lattice subjected to a tilt potential. The total Hamiltonian

[TABLE]

consists of a kinetic term $\mathcal{T}=-J/2\sum_{j}\left(a_{j}^{\dagger}a_{j+1}+\text{h.c.}\right)$, with $a_{j}^{\dagger}$ the creation operator of a particle in site $j$ and $J$ the hopping amplitude; a potential term $\mathcal{E}=E\sum_{j}j\,a_{j}^{\dagger}a_{j}$ due to a constant tilt strength $E$; and an interaction term $\mathcal{V}$. We consider two families of models: the Bose-Hubbard model (BHM) for which $\mathcal{V}=U/2\sum_{j}b_{j}^{\dagger}b_{j}^{\dagger}b_{j}b_{j}$, in which case the particles are bosons $a_{j}=b_{j}$; and the $t$-$V$ model (Ft-VM) with $\mathcal{V}=V\sum_{j}c_{j}^{\dagger}c_{j+1}^{\dagger}c_{j+1}c_{j}$, featuring interacting spinless fermions $a_{j}=c_{j}$.

We will mostly take the initial state $\left|\Phi_{0}\right\rangle$ to be the ground state of the non-tilted system in the presence of a harmonic potential ($H_{0}=\mathcal{T}+\mathcal{V}+W\sum_{j}j^{2}\,a_{j}^{\dagger}a_{j}$), parametrized by the dimensionless constant $\tilde{\rho}=N_{p}\sqrt{W/J}$ Rigol and Muramatsu (2004a, b). The initial condition can be varied form an extended Gaussian-like cloud (small $\tilde{\rho}$) to a highly packed state at large $\tilde{\rho}$. We also consider initial states which are product states, e.g., of the form $\left|\Phi_{0}\right\rangle=a_{i+1}^{\dagger}a_{i+2}^{\dagger}\ldots a_{i+{N_{p}}}^{\dagger}\left|0\right\rangle$. For bosons at $U=\infty$ and for fermions at all $V\neq\infty$, the ground state has this form at large $\tilde{\rho}$.

In addition to $U,V=0$, in both strong interacting regimes ($U,V\to\infty$), the dynamics is that of a set of non-interacting particles. For $U\to\infty$ (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\to\infty$ and $L$ sites can also be shown to map onto that of a Ft-VM with $L-N_{p}$ sites and $V=0$ sup . 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\pi/E$, for any initial state. In fact, exact analytical solutions can be found for the many-body evolution sup . 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^{\prime}$) show the eigenenergies $\varepsilon_{\alpha}$ of $H$, corresponding to eigenvectors $\left|\alpha\right\rangle$, as a function of the interaction strength, color-coded with $\left|c_{\alpha}\right|^{2}$, with $c_{\alpha}=\left\langle\alpha\right.\left|\Phi_{0}\right\rangle$ 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\to\infty$, but leaves Figs. 1($a$,$a^{\prime}$) invariant. Density profiles of the many-body eigenstates which have non-negligible $\left|c_{\alpha}\right|^{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 expands too much the particles cease to interact with each-other. The exponential localization of the many-body eigenstates, and consequently of the dynamics, is thus ensured by the exponential localization of the single particle eigenstates Bloch (1929); Zener (1934).

The effective dimensionality of the Hilbert space spanned by the initial state is $d_{\text{eff}}=\left(\sum_{\alpha}\left|c_{\alpha}\right|^{4}\right)^{-1}$ Popescu et al. (2006), shown in Figs. 1($b$,$b^{\prime}$). This quantity is larger for intermediate interactions than near the ‘free’ points (small or large $U$, $V$). $d_{\text{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 density $\bar{n}_{j}=\lim_{T\to\infty}T^{-1}\int_{0}^{T}dt\,n_{j}(t)$, with $n_{j}(t)=\left\langle a_{j}^{\dagger}(t)a_{j}(t)\right\rangle$ the site occupancy, and $\bar{\sigma}_{j}^{2}=\lim_{T\to\infty}T^{-1}\int_{0}^{T}dt[n_{j}(t)-\bar{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 Rigol et al. (2008) values $\bar{n}_{j}=\sum_{\alpha}\left|c_{\alpha}\right|^{2}\left\langle\alpha\right|n_{i}\left|\alpha\right\rangle$ and $\bar{\sigma}_{j}^{2}=\sum_{\alpha\neq\alpha^{\prime}}\left|c_{\alpha}\right|^{2}\left|c_{\alpha^{\prime}}\right|^{2}\left|\left\langle\alpha\right|n_{i}\left|\alpha^{\prime}\right\rangle\right|^{2}$. Some representative density profiles and $\bar{\sigma}=(\sum_{j}\sigma_{j}^{2})^{1/2}$ for different values of $U$ and $V$ are depicted in Figs. 1($c$,$c^{\prime}$).

For systems fulfilling the so called eigenstate thermalisation hypothesis (ETH) Deutsch (1991); Srednicki (1994); Rigol et al. (2008) 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_{\text{eff}}$, in the sense that different trapped initial states with roughly the same energy yield the same $\bar{n}_{i}$ profile and that long time deviations from the average are suppressed $\bar{\sigma}\propto 1/d_{\text{eff}}$. A comparison between Figs. 1$(b,b^{\prime})$ and $(c,c^{\prime})$ shows that $d_{\text{eff}}^{-1}$ and $\bar{\sigma}$ are qualitatively similar and that the values of $\bar{\sigma}$ 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=\infty$. At these special values the system becomes integrable and the limits $U,V\to 0,\infty$ and $t\to\infty$ do not commute. At these points $\bar{\sigma}$ 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,\infty$. The center of mass $\left\langle x\right\rangle_{t}=(\sum_{j}n_{j}\left(t\right)j)/N_{p}$ and the width $\sigma_{t}=\left[\sum_{j}\left(j-\left\langle x\right\rangle_{t}\right)^{2}n_{j}\left(t\right)/N_{p}\right]^{1/2}$ of the cloud both generically show a “collapse and revival” or beating behavior. Other cloud characteristics (skewness or kurtosis) show the same effect sup . 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, $\propto{U}^{-1}$ at small $U/J$ and $\propto{U}$ at large $U/J$, on the interaction. The behavior at small $U/J$ has recently been observed experimentally Meinert et al. (2014). We have found the same behavior in the fermionic case as a function of $V$ sup .

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^{\nu}$ or $U^{\nu}$ with $\nu=\pm 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\sim 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. Meinert et al. (2014).

Cloud dynamics at ‘free’ points— In contrast with the equilibration seen for moderate interactions, at $U,V=0,\infty$ 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\to\infty$, respectively labeled by F or B. The cases of an initially spread-out and narrow cloud (small and large $\tilde{\rho}$) 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 $\tilde{n}_{F}\left(k,t\right)$ and $\tilde{n}_{B}\left(k,t\right)$. For the bosonic system we also compute the occupation numbers of the natural orbitals $\lambda_{n}(t)$ Penrose and Onsager (1956); Leggett (2001) (with $\lambda_{0}\geq\lambda_{1}\geq...$), defined as the eigenvalues of the single particle density matrix ${\rho_{\text{B}}}_{i,j}(t)=\left\langle b^{\dagger}_{i}(t)b_{j}(t)\right\rangle$. A macroscopic occupation (i.e., a $\lambda_{i}$ of order $\sqrt{N_{p}}$) corresponds to quasi-condensation.

The density profile $n(x,t)$ displays qualitatively different dynamics for small and large $\tilde{\rho}$: Bloch oscillations consist of mainly position oscillations for $\tilde{\rho}\ll 1$ and mainly width oscillations for $\tilde{\rho}\gg 1$. For large $\tilde{\rho}$ 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 $\left\langle x\right\rangle_{t}$ is large for $\tilde{\rho}\ll 1$ and small for $\tilde{\rho}\gg 1$. The cloud width $\sigma_{t}$ shows the opposite behavior. (Fig. 3 right.) This distinction is analogous to that observed in single-particle Bloch oscillations Hartmann et al. (2004). Additional shape dynamics appear at intermediate $\tilde{\rho}$ — the cloud becomes strongly skewed once every period sup . The amplitude of skewness oscillations is non-monotonic as a function of $\tilde{\rho}$, unlike amplitudes of position (width) oscillations which decreases (increases) monotonically with $\tilde{\rho}$ sup .

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 sup . The center of mass has purely sinusoidal oscillations, $\left\langle x\right\rangle_{t}=-\frac{2J}{E}\sin^{2}\left(\frac{tE}{2}\right)\mu_{1}$. The width dynamics is more complicated: $\sigma_{t}^{2}-\sigma_{t=0}^{2}=-\frac{4J^{2}}{E^{2}}\sin^{4}\left(\frac{tE}{2}\right)\mu_{1}^{2}+2\left(\frac{J}{E}\right)^{2}\sin^{2}\left(\frac{E}{2}t\right)\left[1-\cos\left(Et\right)\mu_{2}\right]$. Here $\mu_{a}=\frac{1}{2N_{p}}\sum_{y}\left\langle c_{y}^{\dagger}c_{y-a}+c_{y-a}^{\dagger}c_{y}\right\rangle_{t=0}$. (The behavior of $\mu_{1}$ and $\mu_{2}$ as functions of $N_{p}$ and $\tilde{\rho}$ is described in sup .) This allows to compute the amplitudes of oscillation of the moments, e.g, $\Delta x=\max_{t}\left\langle x\right\rangle_{t}-\min_{t}\left\langle x\right\rangle_{t}$, and $\Delta\sigma^{2}=\max_{t}\sigma_{t}^{2}-\min_{t}\sigma_{t}^{2}$ as a function of $\tilde{\rho}$ for different $N_{p}$. The position oscillation amplitude is $\Delta x=2J/E$ for $\tilde{\rho}\to 0$, and at large $\tilde{\rho}$ decreases as $\Delta x\propto J/\left(E\tilde{\rho}\right)$ for $N_{p}>1$ sup . Conversely, $\Delta\sigma^{2}$ increases from zero to $2(J/E)^{2}$ as $\tilde{\rho}$ is increased ($N_{p}>1$) sup .

The momentum distribution of the fermionic system (F) has simple time evolution: $n_{F}(k,t)=n_{F}(k-\frac{2\pi}{T}t,0)$, reminiscent of single-particle Bloch oscillations sup . For B, the momentum distribution $n_{B}(k,t)$ has similar behavior for small $\tilde{\rho}$, 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 $\tilde{\rho}$ 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 $\{\lambda_{i}\}$ now has two dominant occupancies, $\lambda_{0}\gg\lambda_{n>0}$, signaling a fragmented condensate that is dynamically generated Rigol and Muramatsu (2004c); Cai et al. (2011) 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_{\text{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.

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,\infty$ points for the Bose-Hubbard model (BHM) and the $V=0,\infty$ 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=\infty$ 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\to\infty$ is also the same if the initial state is the same.) Off-diagonal elements in position basis for the bosonic $U=\infty$ 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=\infty$ 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.

S.II.A $V\to\infty$ to $V=0$ mapping for the $t$-$V$ model

The fermionic $t$-$V$ model allows for a mapping between the $V\to\infty$ Hamiltonian (in the reduced Hilbert space sector where the energy is finite) and the $V=0$ Hamiltonian on a smaller chain.

In the $V\to\infty$ limit, states with particles at nearest-neighbour sites are kinematically excluded. For $N_{p}$ fermions in $L$ sites, this means an effective Hilbert space of dimension $\begin{pmatrix}L-N_{p}+1\\ N_{p}\end{pmatrix}$, 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\to\infty$ for a chain of size $L$ and $N_{p}$ fermions onto the Hilbert space for $V=0$ with $N_{p}$ fermions on $L^{\prime}=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, anti-symmetry dictates that the configuration space consists of states with $j_{1}<j_{2}<...<j_{N_{p}}$. 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\to\infty$ 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\to\infty$ 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\to\infty$, 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 case of the $V\to\infty$ 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\to\infty$ 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=\infty$ case from a free-particle calculation. We are not aware of a mapping at the operator level that takes us from the $V=\infty$ 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\to\infty$ 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\to\infty$ bosonic system is described by the same equations.

The Hamiltonian can be written as $H=\mathcal{T}+\mathcal{E}$ with

[TABLE]

where the operators in real and momentum space obey the usual relations $c_{n}=\int_{-\pi}^{\pi}\frac{dk}{2\pi}e^{ikn}c_{k}$, $c_{k}=\sum_{n}e^{-ikn}c_{n}$. For nearest-neighbor hoppings, $\varepsilon\left(k\right)=-J\cos\left(k\right)$, however the following argument holds for a generic dispersion relation. Using the fermionic commutation relations, the evolution operator can be written as

[TABLE]

Applying the evolution operator in this form to the single-particle density matrix in momentum space yields

[TABLE]

where $\left\langle...\right\rangle_{t}$ denotes the mean value taken at time $t$. For $\varepsilon\left(k\right)=-J\cos\left(k\right)$, Eq.(S.4) simplifies to

[TABLE]

so that the real-space correlators are found to be

[TABLE]

using the identity $e^{z\cos\left(\theta\right)}=\sum_{n=-\infty}^{\infty}I_{n}\left(z\right)e^{in\theta}$ where $I_{n}\left(z\right)$ is the modified Bessel function.

In order to compute the moments of the cloud we define the generalized characteristic function

[TABLE]

with

[TABLE]

the generalized moments.

For $a=0$, the $\mu_{a,m}$ are simply the moments of the cloud shape: $\left\langle x^{m}\right\rangle_{t}=N_{p}^{m}\mu_{0,m}\left(t\right)$. For $a\neq 0$, they may be regarded as moments of two-point correlators. In particular, the quantities $\mu_{1}$ and $\mu_{2}$ defined in the main text equal to $\mu_{1,0}\left(t=0\right)$ and $\mu_{2,0}\left(t=0\right)$ respectively. In case the initial state is the ground-state of an harmonic trap one has $\mu_{-a,m}\left(0\right)=\mu_{a,m}\left(0\right)\in\mathbb{R}$ for $m$ even and $\mu_{a,m}\left(0\right)=0$ for $m$ odd. The first non-trivial generalized moments $\mu_{a,m}\left(t\right)$ of the trap ground state are shown in Fig.S2 as function of $\tilde{\rho}$, for different values of $N_{p}$.

Using Eq.(S.6), the characteristic function can be written as

[TABLE]

Taylor expanding the previous expression and identifying the powers of $\lambda$ 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:

[TABLE]

[TABLE]

and

[TABLE]

as a function of the initial values $\mu_{a,m}\left(0\right)$ 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

[TABLE]

and the amplitude of cloud width oscillations

[TABLE]

Fig. S3 shows these two quantities as a function of $\tilde{\rho}$ for different numbers of particles. The quantities plotted, $\Delta{x}(E/J)$ and $\Delta\sigma^{2}(E/J)^{2}$, are scaled to be unit-less and independent of the tilt $E$.

The $\tilde{\rho}$-dependence (for both $\Delta{x}(E/J)$ and $\Delta\sigma^{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$, $\Delta{x}(E/J)$ decreases from 2 to 0, while $\Delta\sigma^{2}(E/J)^{2}$ increases from 0 to 2. This reflects the physics that the Bloch oscillations are primarily position oscillations for small $\tilde{\rho}$ and primarily width oscillations for large $\tilde{\rho}$. (For the single-particle case, the large-$\tilde{\rho}$ limit is different.)

From the exact solutions, $\Delta{x}(E/J)$ and $\Delta\sigma^{2}$ can be expressed in terms of $\mu_{1}=\mu_{1,0}(0)$ and $\mu_{2}=\mu_{2,0}(0)$. It is easy to see that $\Delta{x}(E/J)=2\mu_{1}$. The expression for $\Delta\sigma^{2}$ is more complicated: $\Delta\sigma^{2}E^{2}/J^{2}=\max\left\{\frac{1}{4}\frac{\left(1-\mu_{2}\right){}^{2}}{\left|\mu_{1}^{2}-\mu_{2}\right|},\xi,\frac{\left(3\mu_{2}-4\mu_{1}^{2}+1\right){}^{2}}{4\left|\mu_{1}^{2}-\mu_{2}\right|}\right\}$ if $\left|\frac{3\mu_{2}-4\mu_{1}^{2}+1}{2\left(\mu_{1}^{2}-\mu_{2}\right)}\right|<1$ and $\Delta\sigma^{2}E^{2}/J^{2}=\xi$ otherwise. Here $\xi=2\left|\mu_{2}-2\mu_{1}^{2}+1\right|$.

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 $\tilde{\rho}$. $s\left(T/2\right)$ vanishes for both $\tilde{\rho}=0$ and $\tilde{\rho}\to\infty$. For $\tilde{\rho}=0$ this is due to an almost undeformed cloud evolution. In the $\tilde{\rho}\to\infty$ 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 $\tilde{\rho}=0.53$ panel in Fig. S4). For smallish particle numbers ($N_{p}\lesssim 10$), the $\tilde{\rho}$-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 $\tilde{\rho}=1.72$ panel. For $N_{p}\gtrsim 20$ there is a unique (negative) minimum that shifts to larger $\tilde{\rho}$ 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, $\lambda_{n}$. Figure S5 shows the $\lambda_{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 $\tilde{\rho}=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 $\lambda_{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 $\tilde{\rho}$ induces a Mott insulator state as the initial condition for which the occupation of the natural orbits is given by $\lambda_{n}=1$ for $n<N_{p}$ and $\lambda_{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\in\mathbb{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.

S.III.D Asymptotic behaviors — density profile and fluctuations

The exact solution also allows to compute the long time asymptotic dynamics. From Eq.(S.6), the infinite time average value of the two-point functions is given by

[TABLE]

Since the evolution is periodic, in the quantities $M_{jj^{\prime};nn^{\prime}}$ the integral can be taken over a period of the evolution

[TABLE]

This derivation is completely general for an initial state correlation matrix $\left\langle c_{n}^{\dagger}c_{n^{\prime}}\right\rangle_{0}$.

For simplicity, let us concentrate on the asymptotic form of the density, for which $j=j^{\prime}$, for the special case of initial condition with $\left\langle c_{n}^{\dagger}c_{n^{\prime}}\right\rangle_{0}=\delta_{nn^{\prime}}\left\langle c_{n}^{\dagger}c_{n}\right\rangle_{0}$, i.e., product states as initial states. For this particular case the long time averaged density yields

[TABLE]

with

[TABLE]

In the same way, the fluctuation around the average are given by

[TABLE]

with

[TABLE]

As an example, in Fig.(S6) we present the asymptotic time average starting from the states $\left|\Phi_{0}\right\rangle=\left|...01110...\right\rangle,\left|...0111110...\right\rangle$ and $\left|...0101010...\right\rangle$ for different values of the applied tilt $E$. Note that these quantities are the same for the cases $U=0$, $U\to\infty$ and $V=0$. (This solution does not apply to $V=\infty$, as explained in Section S.II.)

Note that here $\bar{n}_{j}$ and $\bar{\sigma}_{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: $\bar{n}_{j}=\sum_{\alpha}\left|c_{\alpha}\right|^{2}\left\langle\alpha\right|n_{j}\left|\alpha\right\rangle$ and $\bar{\sigma}_{j}^{2}=\sum_{\alpha\neq\alpha^{\prime}}\left|c_{\alpha}\right|^{2}\left|c_{\alpha^{\prime}}\right|^{2}\left|\left\langle\alpha\right|n_{j}\left|\alpha^{\prime}\right\rangle\right|^{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\to 0,\infty$ 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\to\infty$ and $U,V\to 0,\infty$ 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\ll J$ and $V\gg J$, and the BHM at $U\ll J$ and $U\gg J$.

S.IV.A Spectral explanation of beating frequency

For the BHM, the beating frequency is $\propto U^{-1}$ for $U\gg J$ and $\propto U$ for $U\ll 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 $\propto V^{-1}$ for $V\gg J$ and $\propto V$ for $V\ll 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,\infty$) 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 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 single-particle 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\to\infty$ 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\pi/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 hard-core 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.

S.IV.B Beating in different quantities, in both BHM and Ft-VM

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).

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 $\Delta{x}$ and $\Delta{\sigma}$ as function of trap strength

For the non-interacting system ($V=0$ for Ft-VM and $U=\infty$ for BHM), we have shown that Bloch oscillations are primarily position oscillations (large $\Delta{x}$, small $\Delta{\sigma}$) if the trapping is weak and the initial cloud shape is gaussian-like, while the oscillations are predominantly width oscillations (small $\Delta{x}$, large $\Delta{\sigma}$) 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 $\langle{x}\rangle$ and $\sigma$ 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 ($\tilde{\rho}=\infty$) 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 non-interacting 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 for $L\to\infty$. 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 single-particle eigenvalues are nearly equally spaced.

In the main text, we employed the trick of intensity-coding (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.