Exponentially long lifetime of universal quasi-steady states in topological Floquet pumps

Topological phenomena in periodically driven, isolated quantum many body systems are difficult to obtain due to the generic tendency of such systems to heat up towards a maximum entropy state in the long time limit. Here we investigate a mechanism to transiently stabilize topological phenomena over a long-time window, for systems driven at low frequencies. This mechanism provides a route for obtaining long-lived prethermal states with anomalous topological properties, unattainable in equilibrium. We consider the dynamics of interacting particles in a slowly driven lattice with two (or more) bands separated by a bandgap. In the case where the bandgap of the Hamiltonian at any specific time is large we obtain an analytical bound for the rate of change in the number of particles populating these bands. The bound is exponentially small in the ratio between the instantaneous bandgap and the maximum between the driving frequency, interaction strength, and the bandwidth. Within the prethermal time window, this leads to the existence of a quasi-steady state which is characterized by maximum entropy subject to the constraint of fixed number of particles in each band. By initializing the system with a majority of particles in one of the bands, the topological properties of the bands can be measured. For example, this mechanism can be used to obtain quasi-steady states in slowly driven one dimensional systems which carry universal currents, or quasi-steady states of drive-induced Weyl points in 3-dimensional systems.


Introduction
The possibility to dynamically control material properties through a time-periodic driving field has lead to many proposals  and experiments [23][24][25][26] aimed at realizing topological states in periodically driven systems. Topological properties of materials are usually robust and cannot be changed easily using external perturbations. However, Floquet engineering allows for dynamical modification of topological properties using driving fields. These are generally weaker than the bare energy scales of the material system, such as the bandgap in the case of insulators.
A major challenge for Floquet engineering in interacting quantum many body systems is their tendency to absorb energy from the drive and generically heat towards featureless high-entropy states [27][28][29]. In such a state all correlations are trivial, and any topological properties are washed out. Several exceptions to the fate of reaching a featureless state were proposed, such as Floquet integrable or many-body localized systems [30][31][32][33][34][35][36][37][38][39][40].
The featureless high-entropy fate of interacting Floquet systems at long times may not preclude a system from exhibiting topological phenomena at intermediate times. For example, when subjected to high-frequency drives, the heating rates of quantum many body systems are suppressed. In this case, a long prethermal time window can emerge in which the evolution of the system is governed by an effective time-independent Hamiltonian [41][42][43][44][45][46][47][48][49][50]. On the other hand, new types of topological phenomena possible only in periodically driven systems have recently been discovered [3,8,12,17,21,[51][52][53][54][55] which require that the driving frequency is comparable to or smaller than the natural energy scales of the system. We therefore investigate the conditions under which topological phenomena can be exhibited over a long prethermal time window for systems driven at low frequencies [56].
We focus our attention on the dynamics of interacting fermions in slowly driven lattices. We consider a system where the bandstructure of the Hamiltonian at any specific time exhibits two sets of bands separated by a large bandgap. In the limit in which the bandgap is larger than the interaction strength, bandwidths, and the driving frequency, we expect the rates of processes in which particles are scattered across the bandgap to be suppressed relative to intraband scattering processes, due to the large energy mismatch involved in the former. Therefore, the total populations in bands below and above the bandgap are nearly conserved quantities over a long time interval (before interband scattering becomes significant). Due to intraband scattering processes the system quickly attains a high entropy state, subject to the restriction that the initial total band populations are conserved. In Ref. [56], we showed how such a state emerges in a one dimensional system. Moreover, we showed that this state carries a universal quasi-steady state current whose value depends only on the topological character of the Floquet bandstructure and the initial populations of bands.
The goal of this paper is to show that the prethermal regime for low-frequency driving exists across a wide variety of systems in different spatial dimensions, and with different topological properties. To this end we derive a universal and rigorous bound on the interband scattering rates for interacting fermions in slowly driven lattices. The bound demonstrates that the band populations are indeed almost conserved quantities on a timescale which is exponentially long in the ratio between the minimal instantaneous bandgap and the maximum among the bandwidth, driving frequency, and interaction strength.
Our results allow us to go beyond the one dimensional pump analyzed in Ref. [56] and set the stage for finding new types of topological transport phenomena in periodically driven systems with high energy density. An example are driving-induced Weyl points in 3 dimensions [57][58][59]. Furthermore, since the prethermal states obtained in the low driving frequency regime exhibit a homogenous distribution of particles over all momenta within each band, they can also be used to measure the topological character of the Floquet bands (for example, by measuring the average of the Berry curvature to infer the Chern number of the bands [54,55]). This paper is structured as follows: In section 2 we motivate the problem with a generic description of the effects of interactions on topological pumps and describe the setup of the problem. In the main text we restrict ourselves to a 1-dimensional model with 2 bands and on-site interaction. In Appendix G we generalize the calculation to multi-band systems in any number of spatial dimensions, with generic short-range interactions. Here we want to emphasize that our derivation is independent of the topological nature of the system, but applies to all slowly-driven systems with a large bandgap compared compared to driving frequency, bandwidth, and interaction strength. In section 3 we describe a separation of terms of the interaction into nearby and distant terms. Crucial for this are locality properties of the operators which we discuss in section 3.2, together with introducing key quantities and notation for our analytical treatment. By using Lieb-Robinson bounds in section 4.1 we show that the distant terms only yield exponentially small contributions and can be neglected. In section 4.2 we turn to the nearby terms which become the main focus of this paper. We show that the density of excitations grows linearly in time, at a rate which is exponentially suppressed in the ratio of the bandgap over the maximum of driving frequency, bandwidth, and interaction strength. In section 5 we discuss our findings and their implications.

Problem setup and main results
In this section we describe the setup of the problem in this paper and define the basic notations that will be used throughout. Using these definitions, we give a formal statement of the main result.
For simplicity, throughout the main text we will focus on spinless fermions hopping on a one dimensional lattice with two atoms in each unit cell. Our results can be straightforwardly generalized to higher dimensions and larger numbers of bands, as discussed in Appendix G. The hopping of fermions on the lattice is described by the time periodic Bloch Hamiltonian: where A † k and B † k create Bloch states with crystal momentum k on sublattices A and B, respectively, and h AA,k (t), . . . , are periodic functions of time with a period T . A specific lattice Hamiltonian realizing Thouless' one-dimensional pump, which takes the form of Eq. (1), is given in Appendix A. For other models realizing topological pumps, see Figure 1. Plot of the instantaneous bands E i,k (t) of a periodically driven one-dimensional lattice system. Details of the used model are in Appendix A, which is an example of a topological pump. The solid lines are the minimum and maximum over all times t for each momentum k, the shaded areas are the possible instantaneous eigenenergies. We define the energy scale W as the width of the bands, while ∆ denotes the separation of the centers of the bands. For ∆ > W the system is gapped at all times t. In the following we refer to ∆ as the bandgap. [51,[53][54][55][57][58][59][60][61][62][63][64]. The HamiltonianH 0 (t) yields two energy bandsẼ i,k (t), i = 1, 2. We will consider the case where the parameters appearing inH 0 (t) are such that its two instantaneous energy bands are separated by a gap at any time t. We define the maximal bandwidth of the bands as W = max i [max k,t E i,k (t) − min k,t E i,k (t)], and the distance between the bands as ∆ = E 2 − E 1 with the middle of the bands E i = [max k,tẼi,k (t) + min k,tẼi,k (t)]/2. In the following we consider the case ∆ > W for which the system is gapped at all times, and refer to ∆ as the bandgap. These definitions are illustrated in Fig. 1.
Using the time-periodicity of the Hamiltonian in Eq. (1) we find a complete basis of single-particle Floquet states |ψ λ,k (t) = e −i λ,k T |ψ λ,k (t + T ) which solve the Schrödinger equation i∂ t |ψ = H 0 (t)|ψ . Here λ = R, L denotes the two Floquet bands (the choice for the names R and L will be clarified below), and λ,k are the quasi-energies of the states which are visualized in Fig. 2. We use the symbolR † k (t) [L † k (t)] to denote the creation operator for a fermion in the state ψ R,k (t) ψ L,k (t) . We also define the populations . Importantly, we will show that the band populations N R (t) and N L (t) are almost conserved quantities in the prethermal time window that we define below.
We study thermalization of the system in the presence of interactions. For concreteness, Figure 2. Quasienergies for a topological pump with two bands, which are defined modulo ω. For details of the model see Appendix A. The driving frequency ω is smaller than the intrinsic energy scales. This introduces jumps at quasienergy = ± ω 2 . Following the bands throughout the Brillouin zone the periodicity condition gives λ,k+2π = λ,k ± ω, which makes the red band a right-moving R-band and the blue band a left-moving L-band. The red arrows denote a possible intraband scattering process between two right-moving particles, the green arrows an interband process where two right-moving particles are scattered to a right-moving and a left-moving particle.
we consider an on-site two-particle interaction of the form In Appendix G we generalize the interaction to any short-range operator. Our results apply for both repulsive and attractive interactions. In the formulas below we assume U > 0; for U < 0 one should replace U → |U | throughout. Generically such a system will eventually reach a featureless infinite-temperature state where all single-particle Floquet eigenstates are populated equally, implying N R = N L [27]. Here we study the dynamics at intermediate times before reaching the final state, and prove that there exists a long-lived prethermal state in which N R (t) and N L (t) are close to their initial values.
To bound the rate of relaxation of N R (t) and N L (t), we split the interaction term into an intraband partṼ 0 , which preserves the number of particles in each Floquet band, and interband scattering partsṼ + andṼ − which transfer particles from the R-to the L-band and vice versa, respectively: We make the separation in Eq. (3) by expressing the interaction V in Eq. (2) in terms of the L-and R-band creation and annihilation operators,L k (t) andR k (t). The operatorsṼ 0,± are thus time-dependent, due to the time-dependence of the Floquet states encoded inL k (t) andR k (t).
Consider a system in which the interaction term consists only of the band conserving termṼ 0 . In such a system, due to the time dependence of the Hamiltonian, after a short thermalization time τ intra the system reaches a maximal-entropy state subject to the constraint that the band populations N R (t) and N L (t) are conserved. In such a state, all momentum states within each band are equally populated.
We focus on times which are longer than the intraband thermalization time τ intra . For these times the fermions within each band have been fully thermalized. The many-body Floquet eigenstates of the Hamiltonian with the band-conserving part of the interactioñ H 0 (t) +Ṽ 0 (t) are all maximal entropy states indexed by its conserved quantities N R and N L . Therefore the prethermal state established after time τ intra can be represented by any one of these states with the appropriate particle numbers.
When the interband scattering partsṼ ± are included, N R (t) and N L (t) are not strictly conserved. To explore the decay of the quasi-steady prethermal state, we consider an initial state with N R,i and N L,i fermions in the R-and L-band, respectively. We calculate the probability for transferring one particle from the R-to the L-band within Fermi's golden rule, Here T is a long time much larger than the driving period, T T . Our goal in this work is to bound the rate of transfer of particles from the R-to the L-band, τ −1 inter = lim T →∞ P (T )/T . The state |Ψ i (t) in Eq. (4) is a many-body Floquet eigenstate ofH 0 (t) +Ṽ 0 (t), with N R,i (N L,i ) fermions in the R-(L-)band. Likewise, we take the state |Ψ f (t) in Eq. (4) to be a many-body Floquet state ofH 0 (t) +Ṽ 0 (t) with N R − 1 (N L + 1) particles in the R-(L-)band (assuming that N R > N L ).
In Sec. 4.2 we prove that interband transitions are strongly suppressed for ∆ max(ω, W, U ), where ω = 2π/T is the driving angular frequency. Specifically, we show for a general system with short-range interactions in any dimension that the rate of interband transitions is exponentially suppressed, for a positive constant µ. This suppression of interband transitions creates a time window τ intra t τ inter during which the particles within the two bands separately are thermalized, while the densities of particles in the R-and L-bands n R,L = N R,L L remain 6 conserved quantities up to exponentially small corrections. We refer to this state at intermediate times as a prethermal quasi-steady state. The combination of maximal entropy and almost conserved quantities in the quasisteady state leads to the emergence of new universal phenomena. In particular, the homogeneous distribution of particles in the momentum states within each band allows the topological features of the Floquet bandstructure to be manifested in the quasisteady state.
In the main text we consider a generic one-dimensional case with two bands, for which the quadratic part of the Hamiltonian is given by Eq. (1). In Appendix G we generalize our results to systems in any dimension with an arbitrary number of bands. The Floquet bandstructure is obtained by examining the spectrum of the unitary matrix is the matrix appearing in Eq. (1) and P refers to path ordering. In the limit ω/∆ 0 the two Floquet bands may exhibit non-trivial chirality. Each quasienergy band λ,k of U T must be periodic in k. Since the quasienergies are defined in a Floquet-Brillioun zone which is periodic in both k and , the quasienergy bands can wind an integer number of times in the quasienergy direction when k goes from 0 to 2π (we take the lattice constant to be 1). This winding is captured by the integer winding number ν = T 2π 2π 0 dk ∂ ∂k . Bands with non-zero ν are chiral and exhibit a quantized non-zero average group velocity ν/T . Note that the sum of the winding numbers for all the bands must necessarily vanish; for two bands, the winding numbers are equal in magnitude and opposite in sign. The winding number ν is in fact the topological invariant characterizing Thouless' quantized charge pump [3,53,56].
For chiral Floquet bands, the system may carry a nonvanishing current in the quasisteady state. Recall that in the quasi-steady state the distribution of particles is uniform as a function of momentum within each band. Therefore, we expect the contribution to the current carried by each band to be given by the average group velocity of the band times the density of particles occupying this band [56] (setting e = 1 for the charge of the particles). Using the quantization of the averaged group velocity described above, this leads to a simple form for the quasi-steady state current: Remarkably, this result is universal and does not depend on the details of the bandstructure. The mechanism we described here extends beyond the one dimensional example discussed above. In fact, non-trivial topology of the Floquet bands in a variety of slowly driven systems can similarly be reflected in other universal quasi-steady state observables.

Nearby and distant contributions to the excitation rate
In this section we show how to separate the contributions to the excitation probability in Eq. (4) into nearby and distant terms. We define this decomposition in Sec. ??, and in Sec. 3.2 we explicitly show that the intra-and interband operators, V 0 and V ± in Eq. 3, can be written as a sum of local terms. The notations defined in this section will be used throughout the rest of this paper. Our quantity of interest is the rate of change of particle density n R,L (t) in the two bands, Γ = τ −1 inter . Fermi's golden rule, Eq. (4), gives the probability of exciting one particle from the R-to the L-band in time T . Normalizing the probability by 1 LT gives the rate Γ of change of the particle densities n R,L as

Decomposition into nearby and distant terms
The interaction in Eq. (2) is a sum over local terms We rewrite each term in the basis of the Floquet eigenstate operatorsL † k ,R † k . In this basis we can distinguish between terms that conserve the population numbers N R,L and make upṼ 0 , and terms which move particles from the R-to the L-band and formṼ + (or vice versa forṼ − ). For the interband interaction we writẽ The square of the matrix element in Eq. (7) produces a double-sum over x and x , To separate the terms we define a critical distance r * according to which we split the doublesum into L x ,x=1 = L x=1 |x −x|>r * + L x=1 |x −x|≤r * . We refer to the first part as distant and the second part as nearby terms. However, note that due to the projection on the Floquet bands these terms are not strictly local. For this separation to be meaningful we first need to show thatṼ x (t) is quasilocal.

Local representation of the intra-and interband interaction operators
The creation operators for the single-particle Floquet-Bloch states are related to those in the sublattice basis via a unitary transformation, Note that the periodicity condition of the Floquet-Bloch states carries over to the functions that define the change of basis,α k+2π (t) =α k (t) = e i 1,k Tα k (t + T ), and similarly forβ k (t). 8 With this definition we can write the explicit form of the interband scattering interactioñ V + in Eq. (8), wherẽ represents the scattering process of two particles from the R-band to a final state where one particle was moved to the L-band.Ṽ L x (t) describes the interaction of two particles from opposite bands where both end up in the L-band. Here and in the following we only need the terms that scatter exactly one particle from the R-to the L-band. Note that this operatorṼ + is not Hermitian, its adjoint operator isṼ − .
It follows from the k-periodicity of the functionsα,β that at all timesṼ x (t) is localized around x with exponentially small tails. Hence the separation of terms in Eq. (9) yields pairs of operators which are localized nearby each other, or distant from each other.
As mentioned above we want to use Lieb-Robinson bounds to show that the distant terms in Eq. (9) are negligible. The Lieb-Robinson bound can be applied if the Hamiltonian of the system is local at all times t. As mentioned above, time evolution of the states is given byH(t) =H 0 (t) +Ṽ 0 (t). In general the quadratic part in Eq. (1) is a local operator. When rewriting the band population conserving interactionṼ 0 in the band basis (10) we obtain a sum of three types of terms, with the functiong RR ({q i }, t) = −α * q 1 (t)β * q 2 (t)α q 3 (t)β q 4 (t), describes scattering of two particles within the R-band. The band population conserving interaction between one particle in the R-band and one particle in the L-band is described bỹ wherẽ The operatorṼ LL 0 (t, {q i }) contains only creation and annihilation operators in the L-band, . Similarly to our discussion above, periodicity ofα q ,β q ensures that at all timesṼ 0 (t) is local. This property is crucial for the bounds for both the distant and the nearby terms.

Derivation of an upper bound on the excitation rate
With these preliminaries we are prepared to derive a rigorous upper bound on the excitation rate Γ in Eq. (9). To this end we consider the contribution of the distant terms and of the nearby terms separately, as discussed in section 3.1, and derive an upper bound for each.

Lieb-Robinson bound on the contribution of the distant terms
We begin with the distant terms in Eq. (9), where |x − x| > r * . For this we use Lieb-Robinson bounds, which provide a bound on the expectation value of the commutator of two localized operators in the Heisenberg picture, evaluated at different times [65,66]. To obtain a commutator we subtract from Eq. (9) the term This term equals zero because the adjoint operatorṼ x (t) † moves a particle from the Lband to the R-band, however by assumption |Ψ f (t) has one particle more in the Lband than |Ψ i (t) . Next we switch to the Heisenberg picture where the states are timeindependent and the operators evolve with the time-evolution operator of the Hamiltoniañ H(t) =H 0 (t) +Ṽ 0 (t). Then we can sum over all final states, in the subspace of fixed particle numbers this gives f |Ψ f Ψ f | = 1, and rewrite the contribution of the distant terms in (9) as Following the discussion in section 3.2 the HamiltonianH(t) =H 0 (t) +Ṽ 0 (t) is local at all times t, which is a necessary requirement for the application of Lieb-Robinson bounds.
Under this condition, the norm of the commutator of time-dependent local operators in the Heisenberg picture is bounded by where the usual operator norm is defined as Ṽ 2 ≡ sup |ψ|=1 ψ |Ṽ †Ṽ |ψ . For each local term the norm of the interaction is bounded by Ṽ H x ≤ U . C, c, and v are numerical constants, where the latter is referred to as the Lieb-Robinson velocity. After substituting in the bound on the expectation value, we are left with a double time integral and a sum over r = x − x. Converting the sum over r to an integral (setting a = 1 for the lattice spacing), we obtain: For times T < r * v this contribution is exponentially small. If we pick r * large enough the contribution of the distant terms is negligible compared to the nearby terms, which we discuss next.

Bounding the nearby terms
Now we turn to the nearby terms in Eq. (9), As a first step toward bounding the magnitude of this contribution, we use the triangle inequality with show that the contributions of mixed terms with different x and x are bounded by terms with both operators located at the same site x. Hence we can reduce the calculation to same-site terms through the bound where 2r * + 1 is the number of terms with |r| ≤ r * .

Rotating frame transformation space
As discussed in section 2 we know that the instantaneous single-particle eigenstates fall into two bands of width bounded by W , which are separated by a large bandgap ∆ W . For a slowly-driven system with ω ∆ the Floquet eigenstates and the instantaneous eigenstates are almost identical. Most importantly the diagonal elementsẼ 11,k (t),Ẽ 22,k (t) of the Hamiltonian and the instantaneous eigenenergies E 1,k (t), E 2,k (t) agree up to small corrections of order O(ω/∆) 2 [67]. Hence the diagonal elements also fall into a narrow interval, max k,tẼii,k (t) − min k,tẼii,k (t) = W for i = 1, 2. We now perform a rotating frame transformation to shift the diagonal elements to the interval [−W/2, W/2], see illustration in Fig. 3. To this end we define new operators as These operators define Floquet eigenstates which solve the time-dependent Schrödinger equation for the transformed Hamiltonian with The quantitiesẼ ij,k are obtained from the original Hamiltonian in Eq. (1) via the change of basis in Eq. (10). As a consequence, for all k and t the absolute values of these matrix elements are bounded by half the bandwidth, The same applies to the instantaneous eigenenergies in Fig. 3, |E i,k (t)| ≤ W 2 . Furthermore the diagonal elements E 11,k (t), E 22,k (t) of the transformed Hamiltonian have the same periodicity in both k and t as the original matrix elements. When writing the rotated operators in the sublattice basis akin to Eq. (10), the new functions α k (t), β k (t) differ from the ones defined in Eq. (10) only by a fast oscillating phase, It turns out that after this transformation the forms of the intraband scattering operator V 0 in (13) and the interband scattering interaction V + in (8) remain unchanged. For example, the local termṼ x (t) transforms to where in the rotated frame The states, however, obtain fast oscillating phases, because the initial (final) state contains N R,i (N R,f ) creation operators R † k (particles in the R-band) and N L,i (N L,f ) operators L † k (particles in the L-band). Since N L,f = N L,i + 1 and N R,f = N R,i − 1 the matrix element appearing in Fermi's golden rule (23) obtains a phase, For the contribution of the nearby terms in Eq. (23) this means From now on we will only work in the basis of the transformed frame where the diagonal elements of H 0 are of order O(W ), and the system experiences a large effective driving frequency ∆ W . In this way we cast the low-frequency driving problem to a highfrequencey driving one, akin to the situation studied in [46]. However, in the present case an additional driving frequency ω is present in addition to the large frequency ∆.
We will now proceed to derive a bound on the excitation rate from Eq. (31) for ∆ ω, W, U . To this end we use e i∆t = (−i∂ t ) M ∆ −M e i∆t to convert the fast oscillations e i∆t inside the integral in Eq. (31) into a small prefactor controlled by a large denominator 1/∆ M for a large power M (of size to be determined below). Next we perform integration by parts to move the derivative to the matrix element: where and the boundary terms are given by To simplify notation we define the operator O M Note that this condition on a single matrix element does not uniquely define the operator O M x (t). However, in the following it is sufficient to consider one such operator that satisfies Eq. (35). In section 4.2.3 we specify our choice for the operator O M x (t).

Excitation rate space
We apply expression (32) to Eq. (31) and take this as starting point to calculate the excitation probability. This implies taking a square of (32). Using a triangle inequality we bound the cross-terms by the sum of the product of the squares, The boundary terms in Eq. (32), for each x, are oscillatory functions of the total time T . In Appendix B we show that they are bounded by a constant of order O(U/∆) 2 . Due to the 1/T prefactor in the scattering rate the contribution of the boundary terms |G b (T , M, x)| 2 in Eq. (36) decay to zero at long times. The persistent heating comes from the remaining integrals. In the same appendix we show that in leading order they give a T -independent contribution and depend on the norm of the operator O M x defined in Eq. (35). Altogether we obtain as bound on the excitation rate with the usual operator norm, The second term gives a small offset at short times because, by assumption, we have U ∆. At long times T this term decays to zero. The rate of persistent heating is thus controlled by the norm of the operator O M x . In the next section we will show that this norm is exponentially small: where E is the largest of the small energy scales, E = max{K U U, K E W, K ω ω}. The constants K U,E,ω are of order unity and depend on properties of the Floquet single-particle eigenstates.

Bound on the operator norm space
All that remains is to derive a bound on the norm of the operator O M x in Eq. (38). We first consider the effect of taking one derivative of the matrix element, the iterate the procedure. The many-body states |Ψ i,f (t) are defined as solutions of the Schrödinger equation with the Hamiltonian H 0 +V 0 , therefore the Heisenberg equation for time-dependent operators gives with the interband scattering operator Vx in the transformed frame from Eq. (28). There are three different contributions in Eq. (39): a commutator of V x (t) with the quadratic Hamiltonian H 0 (t), a commutator with the band-preserving scattering interaction V 0 (t), and a time-derivative of the interband scattering interaction. We consider these terms individually to show that each yields a small energy scale ω, W, U respectively, together with the prefactor ∆ −1 this gives an operator with a small norm. Then we iterate the process for M derivatives to obtain a factor ((ω, W, U )/∆) M . The calculation for V R x and V L x is very similar in each of these three contributions, therefore we will demonstrate only the calculation for V R x .
Commutator with quadratic Hamiltonian: To calculate the commutators in Eq.
, the remainder are scalar functions. Here we focus on the L † k 1 R † k 2 R k 3 R k 4 term, the other term is computed in the same way.
We first examine the commutator of V x with H 0 , the quadratic part of the Hamiltonian (25). The off-diagonal elements in Eq. (25) correspond to the products R † k L k , L † k R k , which 15 change the occupation numbers of the two bands. However, the initial and final states |Ψ i,f have fixed occupation numbers: |Ψ f has one particle more in the L-band than |Ψ i . The commutator of R † k L k with V x is an operator that preserves the number of particles in each band, while the commutator of L † k R k with V x is an operator that changes the band occupancies by two. In both cases, the resulting matrix element between Ψ f | and |Ψ i vanishes: Therefore we only need to consider commutators involving the diagonal elements of H 0 in (25), corresponding to the operators R † k R k and L † k L k . The commutator gives , each term has a prefactor that is periodic in all k i . In the rotating frame we showed that |E 11,k |, |E 22,k | ≤ W 2 . Hence each term that is generated by the commutator in Eq. (41) has a magnitude that is reduced by a factor of W/2∆ as compared with the original magnitude of V R x (t). The same argument applies to the the commutator with V L x .
Commutator with band-conserving interaction: The commutator of the local operator V R x (t) with the intraband scattering operator from Eq. (13) contains commutations with V RR 0 (t), V LR 0 (t), and V LL 0 (t). As an example, the commutator of V R x (t) with V RR 0 (t) in Eq. (14) is given by After some algebra (see Appendix C for details), the commutator in Eq. (42) is reduced to a sum of 6 terms, each involving a string of 6 creation/annihilation operators (i.e., corresponding to three-particle scattering). One example term is where This resulting term h({k i }, t) is periodic in all k i , as is f R ({k i }, t). Together with the exponential factor e ix(k 1 +k 2 +k 3 −k 4 −k 5 −k 6 ) this ensures localization. The commutators with V LR 0 (t) and V LL 0 (t) yield similar terms as in Eq. (43). Compared with the magnitude of the four-point interaction terms in V x , the magnitudes of the six-point interaction terms generated by the commutator in Eq. (42) are suppressed by a factor U/∆. In Appendix C we show that overall there are 24 such terms. Additionally there are 24 more similar terms coming from the commutator of V 0 with V L x .
Time-derivative of the operator: The third term in Eq. (39) contains the explicit timederivative of the operator V x (t) in Eq. (28). The time-dependent parts therein are the Floquet single-particle creation and annihilation operators R † k (t), L † k (t), R k (t), L k (t) and the functions α k (t), β k (t). Overall we have to calculate the derivative We use the product rule of the derivative to separate this into 4 terms, where in each we take the derivative of terms with one k i only. From the relation between Floquet eigenstates and sublattice basis in Eq. (26) we have, e.g., for k 4 These are products of two functions which by themselves are not periodic in time, but their products are: In the last step we expand this product into its Fourier series. With its help we calculate the derivative in Eq. (46): We rewrite this in terms of the Floquet eigenstates operators in the rotating frame to arrive The L k 4 term in Eq. (49), combined with the operators L † k 1 R † k 2 R k 3 in Eq. (45), yields a final state |Ψ f with the same number of particles in the L-band as |Ψ i . However, as we argued in Eq. (40) the initial and final states |Ψ i,f have fixed occupation numbers and therefore can only be connected by an operator which moves exactly one particle from the R-band to the L-band. Hence the term involving L k 4 in Eq. (49) can be neglected. Thus from the explicit action of the time derivative we obtain a new term which contains the same single-particle operators as V R x (t), with an additional prefactor ω ∆ . We perform analogous calculations for k 1,2,3 in Eq. (45), where each gives a new term involving the same single-particle operators. Thus we overall obtain 4 new terms which are all of the form where for all cases G({k i }, t) is periodic in all k i . In the example from Eq. (49), A similar calculation can be done to obtain the time-derivative of V L x (t). Here we want to note that the functionsα k (t),β k (t) defined in Eq. (10) oscillate rapidly with a frequency ≈ ∆/2. However in the functions α k (t), β k (t) defined after the rotating frame transformation these fast modes are cancelled, see Eq. (26). The strongest modes in the Fourier series in (47) are those with l ≈ 0, and they decay exponentially fast for large l. In [56] it was shown that this decay goes as e −A|l| 3/2 , therefore the sum over all modes l is finite. We will return to this sum later.
Iterated derivatives of the matrix element: So far we discussed the effects of taking one derivative i∂ t of the matrix element in Eq. (39). To calculate the norm of the operator O M x (t) defined in Eq. (35) we need to iterate this step M times. We showed that all three contributions to the derivative yield operators which are similar to V x (t), in the sense that they contain a string of single-particle Floquet mode creation and annihilation operators R † k , L † k , R k , L k which overall move one particle from the R-band to the L-band, and their prefactor is a periodic function in all k i .
When we iterate the derivative, in each step we start with a sum of terms like V x (t) in Eq. (28). From each term we obtain several new terms akin to equations (41), (43) and (51). Starting from one term with n Floquet single-particle operators R † k , L † k , R k , L k the three different contributions from Eq. (39) give • Commutator with H 0 : n terms with n Floquet state operators and factor W/∆, • Commutator with V 0 : 6n terms with n + 2 Floquet state operators and factor U/∆, • Derivative ∆ −1 i∂ t : n terms with n Floquet state operators and factor ω/∆. Therefore the total number of terms which one term with n single-particle operators creates is n + 6n + n = 8n. This process is visualized in Fig. 4. In the m-th step we obtain terms involving different numbers of single-particle operators: the commutator with the intraband scattering interaction V 0 (t) increases that number n by 2, while the commutator with H 0 (t) and the time-derivative leave n unchanged.
As we can see from the discussion above, the longer the operator string (i.e., the larger the number n), the larger the number of additional terms that are generated in each iteration. Hence, to obtain an upper bound on the total number N of terms in O M x (t) we can replace n by the maximal possible length of an operator string generated after m iterations of the derivative, n max (m) = 2m + 2. The iteration begins with two terms involving n = 4 singleparticle operators, V R x and V L x . Hence the total number N of terms after M steps is bounded by In the last step we use the bound (M + 1)! ≤ 5 2 M M , which is valid for all positive M.
Operator norm after iterated derivatives: Now we are in the position to calculate a bound for the norm O M x (t) . From the discussion above we conclude the operator is of the form Here η indexes all terms that were created by the iterated derivative, i.e., it labels the path that was taken on the tree in Fig. 4 but also labels all the different terms on one branch. For each path η we define two integers p, q which count the number of times the commutator is taken with H 0 and V 0 , respectively, out of the M applications of the time derivative in Eq. (39). In Appendix D we show that the norm of this operator can be bounded by The parameter E is the largest of the renormalized energy scales In Appendix D we derive expressions for the dimensionless constants K, K U , K E , K ω . In general they depend on the properties of the Floquet single-particle states and K, K U , K E are of order O(1) for a local Hamiltonian. For ω > W the constant K ω also is of order one, K ω = O(1), while for W > ω it goes as K ω ∼ W ω . In either case the resulting energy scale E in Eq. (55) is proportional to one of the small energy scales U, W, ω.
So far M was a free parameter, Eq. (32) holds for all M . Hence we can freely choose this parameter. The tightest bound is reached for M = µ∆ E with µ = (16eK) −1 , and Eq.
Our main result from this section is the exponential bound on the operator norm in Eq. (56). The most important part is that this norm is exponentially small if ∆ W, U, ω; in the bound, these energy scales are renormalized by constants which depend on properties of the Floquet single-particle eigenstates. To obtain the final result for the excitation rate from the nearby terms we apply the result (56) to the part which is constant in time in Eq. (37): Together with the exponentially suppressed contribution from the distant terms, Eq. (20), this gives an upper bound on the total rate of change of the densities n R,L of particles in the R-and L-band: For times T T r * /v the first term is the dominant contribution. In thermodynamically large systems r * can be chosen large, hence this gives a long time window during which a prethermal quasi-steady state persists.

Summary and discussion
In this paper we showed that in slowly driven systems, assuming the hierarchy of energy scales ∆ U, W, ω, the rate at which particles transition from one band to the other is exponentially small in the large ratio ∆/E. The intraband thermalization rate is not suppressed by ∆/E. Therefore we expect thermalization within each band to occur rapidly compared to the timescale associated with the decay of the imbalance of the band populations. This leads to a prethermal quasi-steady state in which all single-particle eigenstates within one band are populated equally while there is a population imbalance between the two bands. The inverse of the excitation rate in Eq. (57) gives a lower bound on the lifetime of the quasi-steady prethermal state. Although we were motivated by dynamics of topological pumps we stress that our results hold for slowly-driven systems under very general assumptions. We only require that all terms in the Hamiltonian and the interaction are local operators, and the band gap ∆ is the largest energy scale in comparison to the drive frequency ω, the band width W , and the interaction strength U . The dimensionless factor µ which multiplies the ratio of the energy scales in the exponent only depends on the properties of the single-particle eigenstates; in Appendix D we provide an explicit calculation of µ.
The factor µ is defined as µ = (16eK) −1 with the constant K in equations (D. 19) and (D.6). It is calculated by summing over the coefficients of the Floquet-Wannier functions of the states L † k |0 and R † k |0 , see Eq. (D.3) in Appendix D. In an eigenstate with smaller localization length the coefficients of the Floquet-Wannier functions decay faster at large distance, yielding a smaller sum of coefficients. Their localization length depends also depends on the ratios ∆/W and ∆/ω, which implies that the factor µ in the exponent in Eq. (57) also has a dependence on the energy scales. In general the Floquet-Wannier states become more localized if the bandgap increases, i.e. for larger ratio ∆/W or ∆/ω). Hence K decreases as these ratios increase which leads to an increase in µ ∼ K −1 . Therefore there are subleading corrections in the final result in Eq. (57) which enhance the suppression of the transition rate as ∆/E increases.
We obtained an upper bound for the rate Γ of particles being transferred from the Rto the L-band. If the L-band in the initial state |Ψ i is not empty, the opposite process of particles transitioning from the L-to the R-band is also possible. Hence the decay rate of the population imbalance between the R-and the L-bands will be smaller than the rate Γ. If we assume an initial state where the L-band is empty there is no backflow of particles. In this case we can derive a tighter bound on the excitation rate, see Appendix E. This tigher bound is similar to Eq. (57) with the replacement of µ →μ = 2µ. We expect that the rate of exciting particles is largest in the absence of a backflow of particles, i.e. when the L-band is empty. This is in agreement with numerical simulations in [56]. Hence we may expect that the tighter bound also holds when the L-band is partially filled.
In Appendix G we provide generalizations of our bound. First we generalize to any number of spatial dimensions, and second we generalize the interaction. In Eq. (2) we assume a 2-particle interaction which only acts on particles on the same site. In Appendix G we take a general short-range 2-particle interaction and derive a similar exponential bound. In the same spirit it is also possible to assume short-range interactions between multiple particles. Thirdly we generalize to systems with an arbitrary number of bands. In multiband systems we only require that the bands fall into groups which are separated by a large bandgap ∆. Hence our results apply to a wide class of slowly-driven systems.
Throughout the derivation we assumed a system of fermions. A further generalization would be considering a bosonic system. The final form of the operator O M x (t) in Eq. (53) in section 4.2.3 can be derived in the same way, with only minor changes to the final form. When estimating the norm of the operator we used that for fermionic operators This is also true for hardcore bosons. If there is a limit to the number of bosons that can occupy a single site we can derive a similar bound with our method. However in general the norm of a bosonic operator goes as where N is the number of particles in the system. Hence the norm diverges in the thermodynamic limit. However in a steady state we expect the density of particles to be uniform along the system. In Appendix F we show how to use this to derive an upper bound to the excitation rate. This result shows a significant difference to the fermionic case because the exponential suppression is replaced by a stretched exponential, Γ ∼ e −μ(∆/E) 2/3 . While this is the best bound we can obtain using our method, in Appendix F we also argue that in the case of weak interactions, i.e. ∆ w, ω U , the excitation rate is similarly small as in the fermionic case, Γ ∼ e −μ∆/E . It is an interesting question for future research whether there exists an exponential bound also for bosons for moderate interaction strength, or whether in certain systems bosons behave fundamentally different from fermions. Furthermore it is worth investigating under which circumstances a bosonic system quickly approaches this slow rate of heating, and whether an instability can occur before the system reaches a prethermal state.
Following the original argument by Thouless [53], half filling is required in order to obtain a band insulator state in a non-interacting system, and thereby to obtain a quantized pump. This was observed in cold atom experiments with bosons [54] and fermions [55], as well as in quantum dots [61]. Here we show that (see also [56]) in the presence of interparticle interactions, there exist prethermal states with universal values for the currents at arbitrary filling fractions. Thus our results opens new possibilities for observing topological Floquet bandstructures in quantum many body systems.

Appendix A. Lattice model for Thouless' pump
To obtain the graphs of the band structure we use the same model for a 1-dimensional topological charge pump as in reference [56]. The tight-binding Hamiltonian is This is a tight-binding Hamiltonian with alternating nearest-neighbor hopping J(t) and J (t), and alternating on-site potential ±µ(t). For the figures we set J(t) = J 0 + δJ(t) and J (t) = J 0 − δJ(t) with J 0 = 1 and δJ(t) = 1.2 sin(ωt). The on-site potential is chosen as V (t) = 2.5 cos(ωt). In Fig. 2 the frequency is ω = 0.15.

Appendix B. Constant excitation rate and offset
Here we show that from Eq. (36) we obtain an excitation rate which is constant in time plus a decreasing offset. Eq. (36) can also be written as The boundary terms in the first line give a function which is oscillating in T but not increasing. Every derivative i∂ t of the matrix element corresponds to an energy scale ω, U, W .
Together with the small factor 1/∆ we see that the main contribution in the offset is the m = 0 term which has the lowest power of (ω, U, W )/∆. Justification of this qualitative argument follows from the derivation in section 4.2.3. Hence for the boundary terms we calculate In the first step we apply another triangle inequality, so that we only obtain terms with the same time. In the second step we use f |Ψ f (t) Ψ f (t)| = 1 and for the norm V x † (t)V x (t) = U 2 , as in section 4.1. Therefore the offset to the excitation probability goes, in leading order, as (U/∆) 2 which is small for a large bandgap and weak interactions. Hence, the dominant contribution to heating and eventual thermalization is the part which comes from the second line in Eq. (B.1). Here we want to show that this part is independent of time T , i.e. the system has a constant heating rate. To simplify notation we define the operator O M x (t) implicitely by This definition is not unique, however in the following we only need to consider one such In the last line we set E l,m,n = f − i − ∆ + (l − m − n)ω, and sinc(x) = sin(x) x is the cardinal sine function. It is sharply peaked at x = 0 and has a decaying envelope x −1 . Following the standard derivation of Fermi's golden rule for a continuum of final states we argue that the cardinal sine is negligible outside the interval x ∈ [−π, π], which are the two central zeroes. Hence we only need to consider E l,m,n ∈ [− 2π T , 2π T ]. For T > 2T not for all quasi-energies f ∈ [0, ω] exist l, m, n so that E l,m,n falls into this interval. In the thermodynamic limit there is a continuum of final states. By assumption they are thermalized, which means they are indistinguishable and obtain a constant density of states, ρ( f ) = ρ for f ∈ [0, ω]. [27] Hence we can convert the sum over final states into an integral, f → d f ρ( f ), with a constant density of states ρ( f ) = ρ for f ∈ [0, ω]. Thus the fraction of final states that we need to consider is 2T T of all states, and we continue: Here we used translational invariance, i.e. the sum over all sites becomes L x=1 → L and we consider the operator at one example point x 0 . Because the states are indistinguishable the local operator O M x (t) acts on all of them in a similar way, and the matrix element does not depend on the particular choice of Ψ f . To obtain the matrix element and the density of states we calculate in two ways: Combined these show that the matrix element times the density of states is bounded by the norm of the operator. Applying this to Eq. (B.6) we arrive at

Appendix C. Commutator of the interaction with the intraband operator
In this section we show the explicit calculation of the commutator of V x (t) with the intraband operator V 0 (t) in section 4.2.3. We start with Eq. (42), the commutator of V R x with V RR 0 . It contains 2 terms of eight Floquet single-particle operators each, in different order. By using the canonical fermionic anticommutation relations of the Floquet single-particle operators, we bring the operators in the same order to cancel the terms with eight operators. This requires 6 non-trivial commutations according to (C.1), in the end we obtain The result are 6 terms with six operators each. This we apply to Eq. (42) and calculate for one example term: In the first step we use the definitions of the functions f R ({k i }, t) and g RR ({q i }, t) from equations (28) and (13). In the second step we solve the sum over x which yields a δfunction over the momenta q i . In the third step we solve the sums over q 1 and k 2 by using the δ-functions. All 6 terms from Eq. (C.2) have a similar shape. Next we turn to the commutator with V LR 0 . Commuting the single-particle Floquet operators gives These are 4 additional terms, however the prefactor g LR in Eq. (15) is a sum of 4 terms. Therefore this gives a total of 16 new terms which are all of similar structure as Eq. (C.3).
The last contribution is the commutator with V LL 0 , for which we calculate These are 2 more terms with the same structure as Eq. (C.3). Hence in total the commutator of V R x with the intraband interaction V 0 generates 24 new terms. The same calculation with V L x yields the same result, the commutator generates 24 more terms.

Appendix D. Bound on the operator norm
In this appendix we derive the final bound on the operator O M x (t) in section 4.2.3. We start with the expression in Eq. (53), is a periodic function for all k i because it is a product of periodic functions α k , β k , E 11,k , and E 22,k . That implies the Fourier coefficients in the expansion the coefficients go to zero at large |y|. We apply these expansion to Eq. (D.1) and rewrite the Floquet state operators in the local basis {A x i , B x i } according to definition (10): The momenta k i only appear in the exponent thus in the second step we can sum over all of them and obtain Kronecker-δs. We use these to solve the sums over all sites x i . Note that in this expression any operator far away from x has an exponentially small prefactor because either F η or the α, β are small in that case. This is a manifestation of locality of the operator O M x (t). The norm for each on-site operator is A x i = B x i = 1, so we can estimate the norm of O M x (t) by taking absolute values of each individual factor, and we see that the sums over {y i , z i } decouple: As we argued above, the absolute value coefficients α y , β y , E y in the Fourier series (D.3) go to zero exponentially fast for large |y|. Hence the sums over their absolute value are finite and we define the constants Similarly we know that the coefficients F η ({z i }, t) decay exponentially at large |z i | and therefore the sum {z i } |F η ({z i }, t)| converges. We calculate it iteratively. The main idea hereby is that each step in the iterated derivative causes the same change to this sum.
In the initial step M = 0 we have N = 2 operators with . For now we focus on the terms from f R ({k i }, t). The summation over {z i } is bounded by In the main text in section 4.2.3 we discussed the three different kinds of terms that appear from taking the derivative of the matrix element, and there are multiple terms of each kind. The index η labels all these terms including which branch we follow during the iterated derivatives. Next we calculate the change to the sum over the coefficients F η ({z i }) in Eq. (D.5) after taking one derivative.
Commutator with the quadratic Hamiltonian: After one commutation with the quadratic Hamiltonian we obtain four terms where in each the function F η ({k i }, t) is of the form To express the Fourier coefficients F η ({z i }) we use the Fourier series in (D.3). From this we derive The momentum k 4 appears twice in the expression for F η ({k i }, t) which manifests itself in a convolution of the indices z 4 and z. However this has no effect on the summation over the 28 absolute values: . Compared with the bound in Eq. (D.7) this only differs by an additional factor K E = z |E z |. So in every iteration of taking the commutator with the quadratic Hamiltonian the bound on the sum over coefficients |F η ({z i }, t)| is multiplied by a factor K E .
Commutator with the intraband operator: In a very similar fashion the commutator with the intraband term creates a convolution in the Fourier series but with more terms involved. The calculation follows the same lines. Taking the example term in Eq. (43) we have In the Fourier series this gives a convolution of three momenta k 2 , k 3 , k 4 and the five terms which contain them. The coefficients are and summation over all terms gives Compared with the bound in Eq. (D.7) there is an additional factor K 2 α K 2 β . However, commutation with the intraband operator also increases the number of Floquet single-particle operators by 2. As we have seen in Eq. (D.5) these sums decouple and the contribution to the norm from each pair of Floquet operators R † k R k is (K α + K β ) 2 . Hence in every iteration of commuting with the intraband operator the bound on the sum over coefficients is multiplied by Time derivative of the operator: The function F η for one example term when taking the time derivative is calculated in Eq. (49): 14) The Fourier series in k 1 , k 2 , k 3 is simple, in k 4 this gives a convolution of three terms. We obtain the Fourier coefficients as The summation over all {z i } gives As discussed in the main text, the product α k (t)β * k (t) is periodic in t (Eq. (47)). Therefore the absolute value of the coefficients in the Fourier series over time t, (α k β * k ) (l) , goes to zero exponentially fast at large l. Starting from a Zener tunneling problem [56] shows that these coefficients decay as e −A|l| 3/2 . Therefore the sum in Eq. (D.16) is finite and with the definition of a constant we get This differs from the result in Eq. (D.7) by a factor κ(K 2 α + K 2 β ). Every iteration of the application of the time derivative produces this additional factor.
To summarize this: after M iterations of taking the derivative of the matrix element there are 4 from Eq. (D.7) and the initial 4 Floquet operators as an overall prefactor.
• a factor K p E from p commutators with the quadratic Hamiltonian.
In the last step we defined the energy scale E as the maximum of all the involved renormalized energy scales, The expression (D. 19) holds for all M and we can apply it to Eq. (B.8) in the main text.
Appendix E. Initial state with empty L-band In section 5 we mention that we obtain a tighter bound if in the initial state |Ψ i the Lband is empty. This is the case at short times T if the system is initialized in the ground state of the static Hamiltonian before switching on the drive. There are two changes to the derivation, which mostly affect the counting of terms in Eq. (52). First we can neglect the intraband interaction that requires an L-particle in the initial state, V L x . This operator can't act on the state |Ψ i . Secondly, when calculating the commutator of V x with the intraband term V 0 in section 4.2.3 and Appendix C the counting of terms changes. The commutator [V R x , V RR 0 ] is unchanged, from there we obtain 6 new terms. Next we commute with V LR 0 defined in Eq. (15). The matrix element, reduced to the operators, is The part of the commutator where the term to the right first acts on the initial state is zero because by assumption the L-band is empty and L q 3 removes a particle from that band. In the second part L † k 1 creates the only particle in the L-band, hence L q 3 has to act on this particle which gives δ k 1 ,q 3 . With this we arrive at the same shape as in Eq. (C.3), however the prefactor function g LR ({q i }, t) contains 4 pieces, hence this gives 4 additional terms similar to (C.3). The commutator with V LL 0 does not yield any contribution because V x |Ψ i has only one particle in the L-band and V LL 0 can only act on states with at least 2 L-particles. Thus we overall obtain 10 new terms. When iterating these steps the commutator with V LR 0 always only gives 4 terms, independent of the number n of single-particle Floquet operators. The commutator with V RR 0 gives 2n − 2 terms, which in total makes 2n + 2 terms from the commutator.
This changes the counting of operators in section 4.2.3 and the bound on the total number in Eq. (52). We now obtain at most n + (2n + 2) + n = 4n + 2 terms at each step of taking the iterated derivatives. The total number of terms after M steps therefore is bounded bỹ This is in comparison to the bound N ≤ 2(16M ) M in the main text. This changes the final result from Eq. (57) to Hence the overall prefactor is divided by 4, and the factorμ in the exponent is doubled which leads to stronger suppression.

Appendix F. Bosonic pump
In this appendix we consider a system of bosons in contrast to a fermionic system in the main text. We can perform a similar derivation of the operator O M x and arrive at a result akin Eq. (53). The main difference appears when calculating the norm of the operator. In Eq. (D.5) we use for the local norm of a fermionic operator A x = B x = 1. However there can be arbitrarily many bosons on a single site, hence the norm of a bosonic operator In other words, we do not consider the effect of the operator O M x 0 on an arbitrary state, but on the initial state |Ψ i . By assumption this state is prethermal and all A, B sites are equally populated. The occupation number of each site is bounded by the particle density per unit cell,n = N L particles, where N is the total number of particles in the system. The operator O M x in Eq. (53) is a sum of terms which all contain a product of n/2 creation operators and n/2 annihilation operators. n depends on the index η and is bounded by n max = 2M + 4. So there are at most M + 2 annihilation and creation operators. Application of an annihilation operator gives a numerical factor √ m with m the number of particles on that particular site and reduces the number of particles to m − 1. Hence iterated application of the same annihilation operator (i.e. in the sum over {y i , z i } terms with the same location, x + y i + z i = x + y j + z j ) gives a reducing numerical factor, so the effect of each annihilation operator acting on |Ψ i can be bounded by √n . The creation operators are different, they give a numerical factor m + 1 and increase the number of particles m by 1. Terms in Eq. (53) contain up to M + 2 times the same creation operator, yielding an increasing numerical factor. The total numerical factor hence is bounded by M +2 ϕ=1 √n √n + ϕ = √n M +2 Γ(n + M + 3) where Γ is the factorial gamma function, Γ(x + 1) = x! for x ∈ N. Here we can bound where γ is a constant which only depends on the density of particlesn. To calculate the final bound for the bosonic operator this factor multiplies the norm of the fermionic operator in Eq. (D. 19): In the last step we chose M =μ∆/E withμ = (16e √ γK √n ) 2/3 , similar to the last step in Eq. (56) in the main text. The excitation rate therefore is bounded by The main difference to the fermionic case is the appearance of a stretched exponential with (∆/E) 2/3 . The reason for this lies in the additional numerical factor (F.3) that comes from the properties of the bosonic operators, which grows with the order M in the bosonic case.
However, more can be said. In Eq. (39) in section 4.2.3 there are three different contributions to the operator norm which we discussed separately in the following subsections. There we showed that only the commutator with the interband interaction in Eq. (43) leads to an increase in the number of single-particle operators. On the other hand the commutator with the quadratic Hamiltonian and the time-derivative of the operator, equations (41) and (51), leave the number of single-particle operators unchanged (see also the discussion in section 4.2.3). Therefore we argue that if the interaction strength is the weakest of the energy scales, i.e. in the regime ∆ ω, W U , the contribution of the terms with large number of single-particle operators is negligible. Hence in this regime we conjecture the excitation rate is as small as for fermions, whereĒ = max{K E W, K ω ω}. However we can not prove this statement rigorously. In physical terms this means that only the boson-boson interaction can lead to an excitation rate which is not exponentially small. We suppose this is due to possible bunching effects of many bosons on a single site. The low-frequency drive and two-particle scattering alone always lead to an exponentially small excitation rate. It is an interesting question to show this rigorously, and to analyze under which conditions a bunching effect of bosons may or may not happen and cause a larger excitation rate. particle in each group of bands. In momentum space the RR-term is V RR 0 (t) = U x,{d i },{k i } χ ν 3 ,ν 4 ν 1 ,ν 2 ({d i })e i(k 1 (x+d 1 )+k 2 (x+d 2 )−k 3 (x+d 3 )−k 4 x) (G.6) α * ρ 1 ,ν 1 (k 1 , t)α * ρ 2 ,ν 2 (k 2 , t)α ρ 3 ,ν 3 (k 3 , t)α ρ 4 ,ν 4 (k 4 , t) and similarly we obtain the LL-and the LR-terms where the latter has a prefactorg LR which consists of 4 terms, as in Eq. (15) in the main text. The interband part of the interaction which excites one particle from the R-to the L-bands splits into two pieces, one partṼ R + (t) that acts on two particles in the R-band and one part V L + (t) that acts on one particle in the R-and L-band each. The first part becomes V R and similarly forṼ L + (t). We want to calculate the rate at which this interband interaction changes the particle density in the upper band, (G.8) Here the particle density in the L-bands is n L = N L /V . We start by separating nearby and distant terms as in section 4.1, determined by a critical distance r * : x,x To find a bound on the norm of this operator we follow the same steps as in section 4.2.3. We consider iterative application of the M -fold derivative, where in each step the effect is given by the Heisenberg Eq. (39), From now on we only consider the operator V R x (t), the contribution from V L x (t) is calculated analogously.