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

TL;DR

This paper demonstrates that in driven quantum systems, topological quasi-steady states can have exponentially long lifetimes, maintaining topological features while losing non-universal details, with implications for quantum control.

Contribution

The authors derive an analytical bound showing the lifetime of topological quasi-steady states grows exponentially as driving frequency decreases, revealing a new stabilization mechanism.

Findings

Lifetime of quasi-steady states is exponentially large in inverse frequency.

Topological properties persist while non-universal features are washed out.

States are characterized by maximum entropy with fixed particle number in Floquet bands.

Abstract

We investigate a mechanism to transiently stabilize topological phenomena in long-lived quasi-steady states of isolated quantum many-body systems driven at low frequencies. We obtain an analytical bound for the lifetime of the quasi-steady states which is exponentially large in the inverse driving frequency. Within this lifetime, the quasi-steady state is characterized by maximum entropy subject to the constraint of fixed number of particles in the system's Floquet-Bloch bands. In such a state, all the non-universal properties of these bands are washed out, hence only the topological properties persist.

Equations

$${\tilde{\mathcal{H}}}_0(t)=\sum _k\left(\begin{array}{c}{A}_k^{†}{B}_k^{†}\end{array}\right)\left(\begin{array}{cc}{h}_{\mathrm{AA},k}(t)& h_{\mathrm{AB},k}(t)\\ h_{\mathrm{BA},k}(t)& h_{\mathrm{BB},k}(t)\end{array}\right)\left(\begin{array}{c}{A}_k\\ B_k\end{array}\right),$$

$$\mathcal{V}=U\sum _x{A}_x^{†}{A}_x{B}_x^{†}{B}_x.$$

$$\mathcal{V}={\tilde{\mathcal{V}}}_{intra}(t)+{\tilde{\mathcal{V}}}_{+}(t)+{\tilde{\mathcal{V}}}_{++}(t)+{\tilde{\mathcal{V}}}_{-}(t)+{\tilde{\mathcal{V}}}_{--}(t).$$

$$\Gamma (\mathcal{T})=\frac{1}{L\mathcal{T}}\sum _f{\left|{\int }_0^{\mathcal{T}}dt⟨{\tilde{\Psi }}_f(t)|{\tilde{\mathcal{V}}}_{+}(t)|{\tilde{\Psi }}_i(t)⟩\right|}^2.$$

$$\Gamma \sim \text{exp}\left(-\mu \frac{\Delta }{\text{max}(\omega ,W,U)}\right),$$

$$j=(n_R-n_L)\frac{\nu }{T}.$$

$${\tilde{\mathcal{V}}}_{+}(t)=\sum _x{\tilde{\mathcal{V}}}_x(t),$$

$$\Gamma =\frac{1}{L\mathcal{T}}\sum _f{\iint }_0^{\mathcal{T}}d{t}^{\prime }dt\sum _{x^{\prime },x=1}^L⟨{\tilde{\Psi }}_i(t^{\prime })|{\tilde{\mathcal{V}}}_{x^{\prime }}^{†}(t^{\prime })|{\tilde{\Psi }}_f(t^{\prime })⟩⟨{\tilde{\Psi }}_f(t)|{\tilde{\mathcal{V}}}_x(t)|{\tilde{\Psi }}_i(t)⟩.$$

$$\left(\begin{array}{c}{\tilde{R}}_k^{†}(t)\\ {\tilde{L}}_k^{†}(t)\end{array}\right)=\left(\begin{array}{cc}{\tilde{\alpha }}_k(t)& {\tilde{\beta }}_k(t)\\ -{\tilde{\beta }}_k(t{)}^{\ast }{e}^{-i{\tilde{\theta }}_k(t)}& {\tilde{\alpha }}_k(t{)}^{\ast }{e}^{-i{\tilde{\theta }}_k(t)}\end{array}\right)\left(\begin{array}{c}{A}_k^{†}\\ B_k^{†}\end{array}\right).$$

$${\tilde{\mathcal{V}}}_{+}(t)=\sum _x{\tilde{\mathcal{V}}}_x(t)=\sum _x\left({\tilde{\mathcal{V}}}_x^{\mathrm{RR}\to \mathrm{RL}}(t)+{\tilde{\mathcal{V}}}_x^{\mathrm{RL}\to \mathrm{LL}}(t)\right),$$

$${\tilde{\mathcal{V}}}_x^{\mathrm{RR}\to \mathrm{RL}}(t)$$

$${\tilde{\mathcal{V}}}_x^{\mathrm{RL}\to \mathrm{LL}}(t)$$

$${\tilde{\mathcal{V}}}_{intra}(t)={\tilde{\mathcal{V}}}_{intra}^{\mathrm{RR}}(\left\{q_i\right\},t)+{\tilde{\mathcal{V}}}_{intra}^{\mathrm{LR}}(\left\{q_i\right\},t)+{\tilde{\mathcal{V}}}_{intra}^{\mathrm{LL}}(\left\{q_i\right\},t).$$

$${\tilde{\mathcal{V}}}_{intra}^{\mathrm{RR}}(t,\left\{q_i\right\})=\frac{U}{L^2}\sum _{x=1}^L\sum _{\{q_i\}}{e}^{ix(q_1+q_2-q_3-q_4)}{\tilde{g}}_{\mathrm{RR}}(\{q_i\},t){\tilde{R}}_{q_1}^{†}(t){\tilde{R}}_{q_2}^{†}(t){\tilde{R}}_{q_3}(t){\tilde{R}}_{q_4}(t),$$

$${\tilde{\mathcal{V}}}_{intra}^{\mathrm{LR}}(\left\{q_i\right\},t)=\frac{U}{L^2}\sum _{x=1}^L\sum _{\{q_i\}}{e}^{ix(q_1+q_2-q_3-q_4)}{\tilde{g}}_{\mathrm{LR}}(\{q_i\},t){\tilde{L}}_{q_1}^{†}(t){\tilde{R}}_{q_2}^{†}(t){\tilde{L}}_{q_3}(t){\tilde{R}}_{q_4}(t),$$

$${\tilde{g}}_{\mathrm{LR}}(\{q_i\},t)=-\left({\tilde{\alpha }}_{q_1}(t){\tilde{\alpha }}_{q_2}^{\ast }(t)+{\tilde{\beta }}_{q_1}(t){\tilde{\beta }}_{q_2}^{\ast }(t)\right)\left({\tilde{\alpha }}_{q_3}(t{)}^{\ast }{\tilde{\alpha }}_{q_4}(t)+{\tilde{\beta }}_{q_3}^{\ast }(t){\tilde{\beta }}_{q_4}(t)\right)e^{i({\tilde{\theta }}_{q_1}(t)-{\tilde{\theta }}_{q_3}(t))}.$$

$${\tilde{\mathcal{V}}}_{intra}^{\mathrm{LL}}(t,\left\{q_i\right\})=\frac{U}{L^2}\sum _{x=1}^L\sum _{\{q_i\}}{e}^{ix(q_1+q_2-q_3-q_4)}{\tilde{g}}_{\mathrm{LL}}(\{q_i\},t){\tilde{L}}_{q_1}^{†}(t){\tilde{L}}_{q_2}^{†}(t){\tilde{L}}_{q_3}(t){\tilde{L}}_{q_4}(t),$$

$$⟨{\tilde{\Psi }}_i(t)|{\tilde{\mathcal{V}}}_x(t)|{\tilde{\Psi }}_f(t)⟩⟨{\tilde{\Psi }}_f(t^{\prime })|{\tilde{\mathcal{V}}}_{x^{\prime }}(t^{\prime }{)}^{†}|{\tilde{\Psi }}_i(t^{\prime })⟩=0$$

$${\Gamma }_{dist}=\frac{1}{L\mathcal{T}}{\iint }_0^{\mathcal{T}}d{t}^{\prime }dt\sum _{x=1}^L\sum _{|x^{\prime }-x|>r^{\ast }}⟨{\tilde{\Psi }}_i|[{\tilde{\mathcal{V}}}_{x^{\prime }}^{H†}(t^{\prime }),{\tilde{\mathcal{V}}}_x^H(t)]|{\tilde{\Psi }}_i⟩.$$

$$\Vert [{\tilde{\mathcal{V}}}_{x+r}^{H†}(t^{\prime }),{\tilde{\mathcal{V}}}_x^H(t)]\Vert \le C\Vert {\tilde{\mathcal{V}}}_{x+r}^{H†}\Vert \Vert {\tilde{\mathcal{V}}}_x^H\Vert e^{-c\left(|r|-v|t-t^{\prime }|\right)},$$

$${\Gamma }_{dist}$$

$${\Gamma }_{near}=\frac{1}{L\mathcal{T}}\sum _f\sum _{x=1}^L\sum _{|r|\le r^{\ast }}{\iint }_0^{\mathcal{T}}d{t}^{\prime }dt⟨{\tilde{\Psi }}_i|{\tilde{\mathcal{V}}}_{x+r}^{†}(t^{\prime })|{\tilde{\Psi }}_f⟩⟨{\tilde{\Psi }}_f|{\tilde{\mathcal{V}}}_x(t)|{\tilde{\Psi }}_i⟩.$$

$${\mathcal{I}}_x{\mathcal{I}}_{x^{\prime }}^{\ast }+{\mathcal{I}}_{x^{\prime }}{\mathcal{I}}_x^{\ast }\le |{\mathcal{I}}_x{|}^2+|{\mathcal{I}}_{x^{\prime }}{|}^2,$$

$${\Gamma }_{near}\le \frac{2r^{\ast }+1}{L\mathcal{T}}\sum _f\sum _{x=1}^L{\left|{\int }_0^{\mathcal{T}}dt⟨{\tilde{\Psi }}_f|{\tilde{\mathcal{V}}}_x(t)|{\tilde{\Psi }}_i⟩\right|}^2,$$

$$R_k^{†}(t)={\tilde{R}}_k^{†}(t)e^{-i\Delta t/2};L_k^{†}(t)={\tilde{L}}_k^{†}(t)e^{i\Delta t/2}.$$

$${\mathcal{H}}_0(t)$$

$${\mathcal{H}}_0(t)=W\sum _k\left(\begin{array}{cc}{R}_k^{†}& L_k^{†}\end{array}\right)\left(\begin{array}{cc}{e}_{11,k}& e_{12,k}\\ e_{21,k}& e_{22,k}\end{array}\right)\left(\begin{array}{c}{R}_k\\ L_k\end{array}\right).$$

$$\left(\begin{array}{c}{R}_k^{†}(t)\\ L_k^{†}(t)\end{array}\right)=\left(\begin{array}{cc}{\alpha }_k(t)& {\beta }_k(t)\\ -{\beta }_k(t{)}^{\ast }{e}^{-i{\theta }_k(t)}& {\alpha }_k(t{)}^{\ast }{e}^{-i{\theta }_k(t)}\end{array}\right)\left(\begin{array}{c}{A}_k^{†}\\ B_k^{†}\end{array}\right),$$

$${\alpha }_k(t)={\tilde{\alpha }}_k(t)e^{-i\Delta t/2};{\beta }_k(t)={\tilde{\beta }}_k(t)e^{-i\Delta t/2};{\theta }_k(t)={\tilde{\theta }}_k(t).$$

$$\mathcal{V}_{x}(t)=\frac{U}{L^{2}}\sum_{\{k_{i}\}}e^{ix(k_{1}+k_{2}-k_{3}-k_{4})}\Big{(}f_{\rm R}(\{k_{i}\},t)L_{k_{1}}^{\dagger}R_{k_{2}}^{\dagger}R_{k_{3}}R_{k_{4}}+f_{\rm L}(\{k_{i}\},t)L_{k_{1}}^{\dagger}L_{k_{2}}^{\dagger}L_{k_{3}}R_{k_{4}}\Big{)},$$

$$|{\Psi }_i(t)⟩=e^{i(N_{L,i}-N_{R,i})\Delta t/2}|{\tilde{\Psi }}_i(t)⟩;|{\Psi }_f(t)⟩=e^{i(N_{L,f}-N_{R,f})\Delta t/2}|{\tilde{\Psi }}_f(t)⟩,$$

Full text

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

Tobias Gulden1,2∗, Mark S. Rudner3, Erez Berg4,5, and Netanel H. Lindner1

1 Physics Department, Technion, 3200003 Haifa, Israel

2 IST Austria, 3400 Klosterneuburg, Austria

3 Niels Bohr International Academy and Center for Quantum Devices, University of Copenhagen, 2100 Copenhagen, Denmark

4 Physics Department, University of Chicago, IL 60637, USA

5 Faculty of Physics, Weizman Institute of Science, 7610001 Rehovot, Israel

• tgulden@ist.ac.at

Abstract

**We investigate a mechanism to transiently stabilize topological phenomena in long-lived quasi-steady states of isolated quantum many-body systems driven at low frequencies. We obtain an analytical bound for the lifetime of the quasi-steady states which is exponentially large in the inverse driving frequency. Within this lifetime, the quasi-steady state is characterized by maximum entropy subject to the constraint of fixed number of particles in the system’s Floquet-Bloch bands. In such a state, all the non-universal properties of these bands are washed out, hence only the topological properties persist.

**

Contents

1. 1 Introduction

2. 2 Problem setup and main results

3. 2.1 Interaction induced quasi-steady states

4. 2.2 Universality in the quasi-steady state

5. 3 Locality structure of the contributions to the excitation rate

6. 3.1 Decomposition into nearby and distant terms

7. 3.2 Local representation of the intra- and interband interaction operators

8. 4 Derivation of an upper bound on the excitation rate

9. 4.1 Lieb-Robinson bound on the contribution of the distant terms

10. 4.2 Bounding the nearby terms

11. 4.2.1 Rotating frame transformation

12. 4.2.2 Excitation rate

13. 4.2.3 Bound on the operator norm

14. 5 Summary and discussion

15. A Lattice model for Thouless’ pump

16. B Constant excitation rate and offset

17. C Details for the calculations in section 4.2.3

18. C.1 Commutator of the interaction with the intraband operator

19. C.2 Explicit time-derivative

20. C.3 Bound on the operator norm

21. D Initial state with empty L-band

22. E Generalization to arbitrary dimensions and multi-band systems

1 Introduction

The possibility to dynamically control material properties through a time-periodic driving field has lead to many proposals [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22] 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 fields 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, 56] 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 prethermal states can be stabilized over a long time window for systems driven at low frequencies [57]. This includes, but is not limited, to states with topological properties.

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 of the individual bands, 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. [57], we showed how such a state emerges in a particular 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 derive a universal and rigorous bound on the interband scattering rates for interacting fermions in slowly driven lattices. The bound holds for a large class of systems in different spatial dimensions, and with different topological properties. The bound implies that the prethermal regime for low-frequency driving is indeed a general phenomenom which can be found across this wide variety of systems. The bound demonstrates that the band populations are indeed almost conserved quantities on a timescale that 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. [57] and set the stage for finding new types of topological transport phenomena in periodically driven systems with high energy density. Such extensions include, for example, driving-induced Weyl points in 3 dimensions [58, 59, 60]. 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 in slowly-driven systems, in particular in 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 E 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 to driving frequency, instantaneous bandwidth, and interaction strength. In Section 3 we differentiate between two qualitatively different classes of contributions to the excitation rate which we treat separate thereafter. Crucial for this distinction 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 one of the classes above only yields exponentially small contributions and can be neglected. In section 4.2 we bound the remaining contributions. We show that the density of excitations grows at a rate that is exponentially small 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.

2 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 which then will be derived in the subsequent sections.

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 or degrees of freedom, as discussed in appendix E. The dynamics of fermions on the lattice is described by a time-periodic Bloch Hamiltonian

[TABLE]

where $A_{k}^{\dagger}$ and $B_{k}^{\dagger}$ create Bloch states with crystal momentum $k$ on sublattices $A$ and $B$, respectively, and $h_{{\rm AA},k}(t),\ldots,$ 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 [53, 61, 62, 54, 55, 51, 63, 64, 65, 58, 59, 60]. The Hamiltonian $\tilde{\mathcal{H}}_{0}(t)$ yields two instantaneous energy bands $\tilde{E}_{i,k}(t),i=1,2$. We will consider the case where the parameters appearing in $\tilde{\mathcal{H}}_{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}\left[\max_{k,t}\tilde{E}_{i,k}(t)-\min_{k,t}\tilde{E}_{i,k}(t)\right]$, and the distance between the bands as $\Delta=\overline{E}_{2}-\overline{E}_{1}$ with the middle of the bands $\overline{E}_{i}=[\max_{k,t}\tilde{E}_{i,k}(t)+\min_{k,t}\tilde{E}_{i,k}(t)]/2$. In the following we consider the case $\Delta>W$ for which the system is gapped at all times, and refer to $\Delta$ 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 $|\tilde{\psi}_{\lambda,k}(t)\rangle=e^{-i\epsilon_{\lambda,k}T}|\tilde{\psi}_{\lambda,k}(t+T)\rangle$ which solve the Schrödinger equation $i\partial_{t}|\tilde{\psi}\rangle=\tilde{\mathcal{H}}_{0}(t)|\tilde{\psi}\rangle$. Here $\lambda=R,L$ denotes the two Floquet bands (the choice for the names $R$ and $L$ will be clarified below), and $\epsilon_{\lambda,k}$ are the quasi-energies of the Floquet states. When $\tilde{\mathcal{H}}_{0}(t)$ realizes a topological pump, the quasi-energies wind around the Floquet-Brillouin zone in the adiabatic limit (see section 2.2 below for more details). Such a spectrum is visualized in Fig. 2. We therefore use the symbol $\tilde{R}_{k}^{\dagger}(t)$ and $\tilde{L}_{k}^{\dagger}(t)$ to denote creation operators for fermions in the states $|\tilde{\psi}_{R,k}(t)\rangle$ which corresponds to a “right mover” with positive winding, and $|\tilde{\psi}_{L,k}(t)\rangle$ which corresponds to a “left mover” with negative winding. There are actually small gaps from avoided crossings between the right- and left-moving bands in Fig. 2. However these gaps are exponentially small in $\Delta/\omega$ (cf. [57]), hence they can be ignored over exponentially long time scales. We also define the populations $N_{R}(t)=\sum_{k}\tilde{R}_{k}^{\dagger}(t)\tilde{R}_{k}(t)$ and $N_{L}(t)=\sum_{k}\tilde{L}_{k}^{\dagger}(t)\tilde{L}_{k}(t)$. 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.

2.1 Interaction induced quasi-steady states

We study thermalization of the system in the presence of interactions. For concreteness, we consider a two-particle interaction which acts within one unit cell,

[TABLE]

In appendix E 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\to|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)$ remain close to their initial values.

To bound the rate of relaxation of $N_{R}(t)$ and $N_{L}(t)$ we split the interaction term (2) into five different parts, dependent on the change of the number of particles in the two Floquet bands:

[TABLE]

The intraband part $\tilde{\mathcal{V}}_{intra}$ preserves the number of particles in each Floquet band, while the interband scattering parts $\tilde{\mathcal{V}}_{+}$ ($\tilde{\mathcal{V}}_{++}$) and $\tilde{\mathcal{V}}_{-}$ ($\tilde{\mathcal{V}}_{--}$) transfer one (two) particles from the $R$- to the $L$-band and vice versa, respectively. We make this separation by expressing the interaction $\mathcal{V}$ in Eq. (2) in terms of the $L$- and $R$-band creation and annihilation operators, $\{\tilde{L}_{k}(t)^{\dagger},\tilde{L}_{k}(t),\tilde{R}_{k}(t)^{\dagger},\tilde{R}_{k}(t)\}$. The operators $\tilde{\mathcal{V}}_{intra}(t)$, $\tilde{\mathcal{V}}_{\pm}(t)$, $\tilde{\mathcal{V}}_{++}(t)$ and $\tilde{\mathcal{V}}_{--}(t)$ are thus time-dependent, due to the time-dependence of the Floquet states encoded in $\{\tilde{L}_{k}(t)^{\dagger},\tilde{L}_{k}(t),\tilde{R}_{k}(t)^{\dagger},\tilde{R}_{k}(t)\}$.

Consider a system in which the interaction consists only of the band conserving term $\tilde{\mathcal{V}}_{intra}(t)$. Generically, due to the time dependence of the Hamiltonian, after a short thermalization time $\tau_{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 $\tau_{intra}$, when the fermions within each band have been fully thermalized. The many-body Floquet eigenstates $|\tilde{\Psi}(t)\rangle$ of the Hamiltonian with the band-conserving part of the interaction $\tilde{\mathcal{H}}_{0}(t)+\tilde{\mathcal{V}}_{intra}(t)$ are all maximal entropy states indexed by its conserved quantities $N_{R}$ and $N_{L}$. We will represent prethermal states established after time $\tau_{intra}$ by any one of these states with the appropriate particle numbers.

When the interband scattering parts $\tilde{\mathcal{V}}_{\pm}$ 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 rate of transferring particles from the $R$- to the $L$-band within Fermi’s golden rule,

[TABLE]

Here $\Gamma({\mathcal{T}})$ captures the rate of change of the densities of particles in the $R$- and $L$-band, given by $N_{R,i}/L$ and $N_{L,i}/L$ where $L$ is the length of the system, and $\mathcal{T}$ is much larger than the driving period, $\mathcal{T}\gg T$. The state $|\tilde{\Psi}_{i}(t)\rangle$ in Eq. (4) is a many-body Floquet eigenstate of $\tilde{\mathcal{H}}_{0}(t)+\tilde{\mathcal{V}}_{intra}(t)$, with $N_{R,i}$ ($N_{L,i}$) fermions in the $R$- ($L$-)band. Likewise, we take the state $|\tilde{\Psi}_{f}(t)\rangle$ in Eq. (4) to be a many-body Floquet eigenstate of $\tilde{\mathcal{H}}_{0}(t)+\tilde{\mathcal{V}}_{intra}(t)$ with $N_{R,i}-1$ ($N_{L,i}+1$) particles in the $R$- ($L$-)band (assuming that $N_{R}>N_{L}$). The term with the operator operator $\mathcal{V}_{++}(t)$ in Fermi’s golden rule also transfers particles from the $R$- to the $L$-band, however as we argue below the rates of such processes are exponentially suppressed compared to the rate in Eq. (4).

Our goal in this work is to bound the rate of transfer of particles from the $R$- to the $L$-band, $\tau_{inter}^{-1}=\lim_{\mathcal{T}\to\infty}\Gamma(\mathcal{T})$. In section 4.2 we prove that interband transitions are strongly suppressed for $\Delta\gg\text{max}(\omega,W,U)$, where $\omega=2\pi/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,

[TABLE]

for a positive constant $\mu$. We note that the process corresponding to the operator $\tilde{\mathcal{V}}_{++}(t)$, in which two particles must traverse the gap $\Delta$ has a contribution to the rate which scales as $\text{exp}\left(-2\mu\frac{\Delta}{\text{max}(\omega,W,U)}\right)$ and is therefore negligible in comparison with $\Gamma(\mathcal{T})$ in Eq. (5). This suppression of interband transitions creates a time window $\tau_{intra}\ll t\ll\tau_{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}=\frac{N_{R,L}}{L}$ remain 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 quasi-steady 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.

2.2 Universality in the quasi-steady 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). The Floquet bandstructure is obtained by examining the spectrum of the unitary matrix $U_{T}=\mathcal{P}e^{-i\int_{0}^{T}dt\,\hat{h}_{0,k}(t)}$, where $\hat{h}_{0,k}(t)$ is the matrix appearing in Eq. (1) and $\mathcal{P}$ refers to path ordering. In the limit $\omega/\Delta\ll 1$ the two Floquet bands may exhibit non-trivial chirality. Each quasienergy band $\epsilon_{\lambda,k}$ of $U_{T}$ must be periodic in $k$. Since the quasienergies are defined in a Floquet-Brillouin zone which is periodic in both $k$ and $\epsilon$, the quasienergy bands can wind an integer number of times in the quasienergy direction when $k$ goes from [math] to $2\pi$ (we take the lattice constant to be $1$). This winding is captured by the integer winding number $\nu=\frac{T}{2\pi}\int_{0}^{2\pi}dk\frac{\partial\epsilon}{\partial k}$. Bands with non-zero $\nu$ are chiral and exhibit a quantized non-zero average group velocity $\nu/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. For any non-zero value of $\omega$, avoided crossings in the Floquet spectrum open at the intersections between counter-propagating Floquet bands. As a result, when these avoided crossings are resolved, all Floquet bands have zero winding number [57, 66]. However, in the limit $\omega/\Delta\ll 1$, the corresponding gaps at these avoided crossings are exponentially suppressed in $\Delta/\omega$. Therefore, chiral Floquet bands with non-zero values of $\nu$ can be defined in the limit $\omega/\Delta\ll 1$, in which the avoided crossings can be ignored, see Fig. 2. The winding number $\nu$ is in fact the topological invariant characterizing Thouless’ quantized charge pump [53, 3, 57].

For chiral Floquet bands, the system may carry a nonvanishing current in the quasi-steady 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 $n_{R,L}$ occupying this band [57] (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:

[TABLE]

Remarkably, this result is universal and does not depend on the details of the bandstructure. The mechanism we describe 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. In appendix E we generalize our results to systems in any dimension with an arbitrary number of bands to prove the general existence of quasi-steady states.

3 Locality structure of the contributions to the excitation rate

In this section we show how to separate the contributions to the excitation rate in Eq. (4) into two types of terms which we call nearby and distant terms. We define this decomposition in section 3.1. In section 3.2 we explicitly show that the intra- and interband interaction operators, $\tilde{\mathcal{V}}_{intra}$, $\tilde{\mathcal{V}}_{\pm}$, and $\tilde{\mathcal{V}}_{\pm\pm}$ in Eq. (3), can each be written as a sum of quasi-local terms (see definition below), which is crucial for this separation. The notations defined in this section will be used throughout the rest of this paper.

3.1 Decomposition into nearby and distant terms

The interaction in Eq. (2) is a sum over local terms $A_{x}^{\dagger}A_{x}B_{x}^{\dagger}B_{x}$. We rewrite each term in the basis of the Floquet eigenstate operators $\tilde{L}_{k}^{\dagger},\tilde{R}_{k}^{\dagger}$. In this basis we can distinguish between terms that conserve the population numbers $N_{R,L}$, which make up $\tilde{\mathcal{V}}_{intra}$, and terms which move particles from the $R$- to the $L$- band forming $\tilde{\mathcal{V}}_{+}$ and $\tilde{\mathcal{V}}_{++}$ (or vice versa for $\tilde{\mathcal{V}}_{-}$ and $\tilde{\mathcal{V}}_{--}$). For the interband interaction operator which moves one particle we write

[TABLE]

where each $\tilde{\mathcal{V}}_{x}(t)$ is the corresponding component of $A_{x}^{\dagger}A_{x}B_{x}^{\dagger}B_{x}$. The square of the matrix element in Eq. (4) produces a double-sum over $x$ and $x^{\prime}$ and double-integral over $t$ and $t^{\prime}$,

[TABLE]

To separate the terms by distance we define a fixed radius $r^{*}$ according to which we split thedouble-sum into $\sum_{x^{\prime},x=1}^{L}=\sum_{x=1}^{L}\sum_{|x^{\prime}-x|>r^{*}}+\sum_{x=1}^{L}\sum_{|x^{\prime}-x|\leq 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 the operators $\tilde{\mathcal{V}}_{x}(t)$ are not strictly local. For this separation to be meaningful we first need to show that $\tilde{\mathcal{V}}_{x}(t)$ is quasi-local. To define quasi-locality around the site $x$ we express an operator $\tilde{\mathcal{O}}_{x}$ in terms of annihilation and creation operators in the site basis, i.e., the operators $A^{\dagger}_{x^{\prime}},B^{\dagger}_{x^{\prime}},A_{x^{\prime}},B_{x^{\prime}}$. We call the operator $\tilde{\mathcal{O}}_{x}$ quasi-local around $x$ if the coefficients in this expansion decay at least exponentially with the distance $r=|x-x^{\prime}|$ of any contributing creation/annihilation operator from the site $x$.

3.2 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, which in general can be written as

[TABLE]

Note that the periodicity condition of the Floquet-Bloch states carries over to the basis transformation functions, $\tilde{\alpha}_{k+2\pi}(t)=\tilde{\alpha}_{k}(t)=e^{i\epsilon_{R,k}T}\tilde{\alpha}_{k}(t+T)$, and similarly for $\tilde{\beta}_{k}(t)$. With this definition we can write the explicit form of the interband interaction operator $\tilde{\mathcal{V}}_{+}$ in Eq. (7),

[TABLE]

where

[TABLE]

with the two functions $\tilde{f}_{\rm R}(\{k_{i}\},t)=\left(\tilde{\alpha}_{k_{1}}(t)\tilde{\alpha}_{k_{2}}^{*}(t)+\tilde{\beta}_{k_{1}}(t)\tilde{\beta}_{k_{2}}^{*}(t)\right)\tilde{\alpha}_{k_{3}}(t)\tilde{\beta}_{k_{4}}(t)e^{i\tilde{\theta}_{k_{1}}(t)}$ and $\tilde{f}_{\rm L}(\{k_{i}\},t)=\tilde{\alpha}_{k_{1}}(t)\tilde{\beta}_{k_{2}}(t)e^{i(\tilde{\theta}_{k_{1}}(t)+\tilde{\theta}_{k_{2}}(t)-\tilde{\theta}_{k_{3}}(t))}\left(\tilde{\alpha}_{k_{3}}(t)^{*}\tilde{\alpha}_{k_{4}}(t)+\tilde{\beta}_{k_{3}}^{*}(t)\tilde{\beta}_{k_{4}}(t)\right)$. $\tilde{\mathcal{V}}_{x}^{\rm{RR}\to\rm{RL}}(t)$ represents the scattering process of two particles from the $R$-band to one $R$- and one $L$-particle. $\tilde{\mathcal{V}}_{x}^{\rm{RL}\to\rm{LL}}(t)$ describes a scattering process of two particles from opposite bands when both end up in the $L$-band. Note that the operator $\tilde{\mathcal{V}}_{+}$ is not Hermitian, its adjoint operator is $\tilde{\mathcal{V}}_{-}$.

The functions $\tilde{\alpha}_{k},\tilde{\beta}_{k}$ are periodic in $k$, and in general their Fourier series only contains low modes in $k$. From this it follows that at all times $\tilde{\mathcal{V}}_{x}(t)$ is quasi-local around $x$, as defined above. Hence we classify each term in Eq. (8), which involves operators $\tilde{\mathcal{V}}_{x}(t)$ and $\tilde{\mathcal{V}}_{x^{\prime}}(t^{\prime})$, as “nearby” if $|x-x^{\prime}|\leq r^{*}$ or “distant” if $|x-x^{\prime}|>r^{*}$.

In section 4.1 we will use Lieb-Robinson bounds to show that the distant terms in Eq. (8) are negligible. Lieb-Robinson bounds can be applied if the Hamiltonian of the system is a sum of quasi-local terms at all times $t$. As mentioned in section 2, time evolution of the states $|\tilde{\Psi}_{i/f}(t)\rangle$ is given by the Hamiltonian $\tilde{\mathcal{H}}(t)=\tilde{\mathcal{H}}_{0}(t)+\tilde{\mathcal{V}}_{intra}(t)$. In general the quadratic part, Eq. (1), is a sum of local operators. To analyze locality properties of the band population conserving interaction $\tilde{\mathcal{V}}_{intra}$ we rewrite it in the Floquet band basis defined in Eq. (9). This gives a sum of three types of terms,

[TABLE]

Here

[TABLE]

with the function $\tilde{g}_{\rm RR}(\{q_{i}\},t)=-\tilde{\alpha}_{q_{1}}^{*}(t)\tilde{\beta}_{q_{2}}^{*}(t)\tilde{\alpha}_{q_{3}}(t)\tilde{\beta}_{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 by

[TABLE]

where

[TABLE]

The operator $\tilde{\mathcal{V}}^{LL}_{0}(t,\left\{q_{i}\right\})$ contains only creation and annihilation operators in the $L$-band,

[TABLE]

with $\tilde{g}_{\rm LL}(\{q_{i}\},t)=-\tilde{\beta}_{q_{1}}(t)\tilde{\alpha}_{q_{2}}(t)\tilde{\beta}_{q_{3}}^{*}(t)\tilde{\alpha}_{q_{4}}^{*}(t)$. Using similar reasoning to the argument demonstrating quasi-locality of $\tilde{\mathcal{V}}_{x}$ above, the periodicity of $\tilde{\alpha}_{q},\tilde{\beta}_{q}$ can be used to show that at all times $\tilde{\mathcal{V}}_{intra}(t)$ is a sum of quasi-local terms. This property is crucial for the bounds for both the distant and the nearby terms which we will derive in the next section.

4 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 $\Gamma$ in Eq. (8). To this end we consider separately the contribution of the distant terms and of the nearby terms, as discussed in section 3.1, and derive an upper bound for each.

4.1 Lieb-Robinson bound on the contribution of the distant terms

We begin with the distant terms in Eq. (8), where $|x^{\prime}-x|>r^{*}$. Lieb-Robinson bounds [67, 68] provide a bound on the expectation value of the commutator of two quasi-local operators, evaluated at different times in the Heisenberg picture. We aim to rewrite Eq. (8) with a commutator. Therefore we subtract the term

[TABLE]

from the argument of the sum over $|x-x^{\prime}|>r^{*}$. This term equals zero because the adjoint operator $\tilde{\mathcal{V}}_{x^{\prime}}(t^{\prime})^{\dagger}$ moves a particle from the $L$-band to the $R$-band, however by assumption $|\tilde{\Psi}_{f}(t^{\prime})\rangle$ has one particle more in the $L$-band than $|\tilde{\Psi}_{i}(t^{\prime})\rangle$. Next we switch to the Heisenberg picture where the states are time-independent and the operators evolve with the time-evolution operator corresponding to the Hamiltonian $\tilde{\mathcal{H}}(t)=\tilde{\mathcal{H}}_{0}(t)+\tilde{\mathcal{V}}_{intra}(t)$. Then we can sum over all final states, which in the subspace of fixed particle numbers in each band gives $\sum_{f}|\tilde{\Psi}_{f}\rangle\langle\tilde{\Psi}_{f}|=\mathbbm{1}$, and rewrite the contribution of the distant terms in Eq. (8) as

[TABLE]

Here $\tilde{\mathcal{V}}^{H}_{x}(t)$ is the equivalent operator in the Heisenberg picture.

Following the discussion in section 3.2, the Hamiltonian $\tilde{\mathcal{H}}(t)=\tilde{\mathcal{H}}_{0}(t)+\tilde{\mathcal{V}}_{intra}(t)$ is a sum of quasi-local terms at all times $t$, which is a necessary requirement for the application of the Lieb-Robinson bound. Under this condition, the norm of the commutator of time-dependent quasi-local operators in the Heisenberg picture is bounded by [67, 68]

[TABLE]

where the usual operator norm is defined as $\|\tilde{\mathcal{V}}\|^{2}\equiv sup_{|\tilde{\psi}|=1}\langle\tilde{\psi}|\tilde{\mathcal{V}}^{\dagger}\tilde{\mathcal{V}}|\tilde{\psi}\rangle$. For each quasi-local term the norm of the interaction is bounded by $\|\tilde{\mathcal{V}}_{x}^{H}\|\leq U$. In Eq. (18) $C$, $c$, and $v$ are numerical constants, where $v$ is referred to as the Lieb-Robinson velocity. After substituting the bound (18) into Eq. (17) and using translational invariance we are left with a double time integral and a sum over $r=|x^{\prime}-x|$. Converting the sum over $r$ into an integral (with $a=1$ for the lattice spacing), we obtain:

[TABLE]

For times $\mathcal{T}<\frac{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.

4.2 Bounding the nearby terms

We now turn to the nearby terms in Eq. (8),

[TABLE]

We proceed in two steps. First we move to a rotating frame in which all instantaneous eigenstates of $\tilde{\mathcal{H}}_{0}$ lie within a narrow energy window of width $W$. Through this transformation the gap $\Delta$ turns into a high-frequency drive, in addition to the low-frequency drive $\omega$. In a system driven by a high frequency, the excitation rate is strongly suppressed, as shown in Ref. [46]. We use this intuition to derive a bound in the presence of both a high- and a low-frequency drive.

As a first step toward bounding the magnitude of the contribution in Eq. (20), we use the triangle inequality

[TABLE]

with $\mathcal{I}_{x}=\int_{0}^{\mathcal{T}}dt\langle\tilde{\Psi}_{f}|\tilde{\mathcal{V}}_{x}(t)|\tilde{\Psi}_{i}\rangle$, to show that the contributions of mixed terms with different $x$ and $x^{\prime}$ are bounded by terms with both operators centered at the same site $x$. Hence we can reduce the calculation to same-site terms through the bound

[TABLE]

where $2r^{*}+1$ is the number of terms with $|r|\leq r^{*}$.

4.2.1 Rotating frame transformation

space

As discussed in section 2 we assume that at all times $t$ the instantaneous single-particle eigenstates fall into two bands of width bounded by $W$, which are separated by a large bandgap $\Delta\gg W$. For a slowly-driven system with $\omega\ll\Delta$ the Floquet eigenstates and the instantaneous eigenstates are almost identical. Most importantly, the diagonal elements $\tilde{E}_{11,k}(t),\tilde{E}_{22,k}(t)$ of the Hamiltonian in the basis $\{\tilde{L}_{k}^{\dagger}(T),\tilde{R}_{k}^{\dagger}(t)\}$ and the instantaneous eigenenergies $E_{1,k}(t),E_{2,k}(t)$ agree up to small corrections of order $\mathcal{O}(\omega/\Delta)^{2}$ [69]. Hence the diagonal elements also fall into a narrow interval, $\max_{k,t}\tilde{E}_{ii,k}(t)-\min_{k,t}\tilde{E}_{ii,k}(t)\leq 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

[TABLE]

These operators define Floquet eigenstates which solve the time-dependent Schrödinger equation for the transformed Hamiltonian

[TABLE]

The quantities $\tilde{E}_{ij,k}(t)$ are obtained from the original Hamiltonian in Eq. (1) via the change of basis in Eq. (9). As a consequence, for all $k$ and $t$ the absolute values of the new diagonal elements are bounded by half the bandwidth, $|E_{11,k}(t)|,|E_{22,k}(t)|\leq\frac{W}{2}$. The same applies to the instantaneous eigenenergies in Fig. 3, $|E_{i,k}(t)|\leq\frac{W}{2}$. We define dimensionless quantities $e_{ij,k}(t)\equiv E_{ij,k}(t)/W$ and rewrite the Hamiltonian in Eq. (4.2.1) as

[TABLE]

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. (9),

[TABLE]

the new functions $\alpha_{k}(t),\beta_{k}(t)$ differ from the ones defined in Eq. (9) only by a fast oscillating phase,

[TABLE]

It turns out that after this transformation the forms of the intraband interaction operator $\mathcal{V}_{intra}$ in Eq. (12) and the interband interaction $\mathcal{V}_{+}$ in Eq. (7) remain unchanged. For example, the quasi-local term $\tilde{\mathcal{V}}_{x}(t)$ transforms to

[TABLE]

where in the rotated frame $f_{\rm R}(\{k_{i}\},t)=\left(\alpha_{k_{1}}(t)\alpha_{k_{2}}^{*}(t)+\beta_{k_{1}}(t)\beta_{k_{2}}^{*}(t)\right)\alpha_{k_{3}}(t)\beta_{k_{4}}(t)e^{i\theta_{k_{1}}(t)}$ and $f_{\rm L}(\{k_{i}\},t)=\alpha_{k_{1}}(t)\beta_{k_{2}}(t)e^{i(\theta_{k_{1}}(t)+\theta_{k_{2}}(t)-\theta_{k_{3}}(t))}\left(\alpha_{k_{3}}^{*}(t)\alpha_{k_{4}}(t)+\beta_{k_{3}}^{*}(t)\beta_{k_{4}}(t)\right)$. The states, however, obtain fast oscillating phases,

[TABLE]

because the initial (final) state contains $N_{R,i}$ ($N_{R,f}$) creation operators $R_{k}^{\dagger}$ (particles in the $R$-band) and $N_{L,i}$ ($N_{L,f}$) operators $L_{k}^{\dagger}$ (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 (22) obtains a phase,

[TABLE]

For the contribution of the nearby terms in Eq. (22) this means

[TABLE]

From now on we will only work in the basis of the transformed frame where the diagonal elements of $\mathcal{H}_{0}$ are of order $\mathcal{O}(W)$, and the system is time-dependent with a large effective effective frequency $\Delta\gg W$. In this way we cast the low-frequency driving problem to a high-frequencey driving one, akin to the situation studied in [46]. However, in the present case an additional driving frequency $\omega$ is present in addition to the large frequency $\Delta$.

We will now proceed to derive a bound on the excitation rate from Eq. (31) for $\Delta\gg\omega,W,U$. To this end we use the identity $e^{i\Delta t}=(-i\partial_{t})^{M}\Delta^{-M}e^{i\Delta t}$ to convert the fast oscillations $e^{i\Delta t}$ inside the integral in Eq. (31) into a small prefactor controlled by a large denominator $1/\Delta^{M}$ for a large power M (the value of $M$ will be determined in Eq. (50)). Using this identity and then performing integration by parts to move the derivative to the matrix element, we obtain

[TABLE]

where

[TABLE]

and the boundary terms are given by

[TABLE]

To simplify notation we define the operator $\mathcal{O}_{x}^{M}(t)$ via

[TABLE]

Note that this condition on a single matrix element does not uniquely define the operator $\mathcal{O}_{x}^{M}(t)$. However, in the following it is sufficient to consider one such operator that satisfies Eq. (35) for Floquet eigenstates $|\Psi_{i/f}(t)\rangle$. In section 4.2.3 we specify our choice for the operator $\mathcal{O}_{x}^{M}(t)$.

4.2.2 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 squares,

[TABLE]

The boundary terms in Eq. (34), for each $x$, are oscillatory functions of the total time $\mathcal{T}$. In appendix B we show that they are bounded by a constant of order $\mathcal{O}(U/\Delta)^{2}$. Due to the $1/\mathcal{T}$ prefactor in the scattering rate the contribution of the boundary terms $|\mathcal{G}_{b}(\mathcal{T},M,x)|^{2}$ in Eq. (36) decay to zero at long times. The persistent heating comes from the remaining integral in Eq. (33). In the same appendix we show that in leading order these integrals give a $\mathcal{T}$-independent contribution in Eq. (36) and depend on the norm of the operator $\mathcal{O}_{x}^{M}$ defined in Eq. (35). Altogether we obtain as bound on the excitation rate

[TABLE]

with the usual operator norm, $\|\mathcal{O}_{x}^{M}\|^{2}\equiv sup_{|\Psi|=1}\left\langle\Psi\right|{\mathcal{O}_{x}^{M}}^{\dagger}\mathcal{O}_{x}^{M}\left|\Psi\right\rangle$. The second term gives a small offset at short times because, by assumption, we have $U\ll\Delta$. At long times $\mathcal{T}$ this term decays to zero. The rate of persistent heating is thus controlled by the norm of the operator $\mathcal{O}_{x}^{M}$. In the next section we will show that this norm is exponentially small:

[TABLE]

where $K$ is a dimensionless constant and $\mathcal{E}$ is the largest of the small energy scales, $\mathcal{E}=\text{max}\{K_{U}U,K_{E}W,K_{\omega}\omega\}$. The constants $K_{U,E,\omega}$ are of order unity and depend on properties of the Floquet single-particle eigenstates.

4.2.3 Bound on the operator norm

space

All that remains is to derive a bound on the norm of the operator $\mathcal{O}_{x}^{M}$ in Eq. (35). We first consider the effect of taking one derivative of the matrix element, then we iterate the procedure. The many-body states $\left|\Psi_{i/f}(t)\right\rangle$ are defined as solutions of the Schrödinger equation with the Hamiltonian $\mathcal{H}_{0}(t)+\mathcal{V}_{intra}(t)$, therefore the matrix elements $\left\langle\Psi_{f}(t)\right|\mathcal{V}_{x}(t)\left|\Psi_{i}(t)\right\rangle$ obey

[TABLE]

with the interband scattering operator $\mathcal{V}_{x}(t)$ in the transformed frame from Eq. (28). There are three different contributions in Eq. (39): a commutator of $\mathcal{V}_{x}(t)$ with the quadratic Hamiltonian $\mathcal{H}_{0}(t)$, a commutator with the band-preserving interaction $\mathcal{V}_{intra}(t)$, and a time-derivative of the interband interaction operator $\mathcal{V}_{x}(t)$. We consider these terms individually to show that each yields one of the energy scales $\omega,W,U$, together with the prefactor $\Delta^{-1}$ this gives an operator with a small norm. Then we iterate the process for $M$ derivatives to obtain a factor $(U^{q}W^{p}\omega^{M-p-q})/\Delta)^{M}$, where $p$ and $q$ are positive integers. The calculation for $\mathcal{V}_{x}^{\rm RR\to RL}$ and $\mathcal{V}_{x}^{\rm RL\to LL}$ is very similar in each of these three contributions, therefore we will demonstrate only the calculation for $\mathcal{V}_{x}^{\rm RR\to RL}$.

Commutator with quadratic Hamiltonian:

To calculate the commutators in Eq. (39) the relevant parts of $\mathcal{V}_{x}$ in Eq. (28) are the operators $L_{k_{1}}^{\dagger}R_{k_{2}}^{\dagger}R_{k_{3}}R_{k_{4}}$ and $L_{k_{1}}^{\dagger}L_{k_{2}}^{\dagger}L_{k_{3}}R_{k_{4}}$; the remainder are scalar functions. Here we focus on the $L_{k_{1}}^{\dagger}R_{k_{2}}^{\dagger}R_{k_{3}}R_{k_{4}}$ term. The second term is computed in the same way.

We first examine the commutator of $\mathcal{V}_{x}$ with $\mathcal{H}_{0}$, the quadratic single-particle Hamiltonian (4.2.1). The off-diagonal elements in Eq. (4.2.1), which contain $R_{k}^{\dagger}L_{k}$ or $L_{k}^{\dagger}R_{k}$, change the occupation numbers of the two bands. When commuted with $\mathcal{V}_{x}$, the $R_{k}^{\dagger}L_{k}$-term gives an operator that preserves the occupation numbers, and the $L_{k}^{\dagger}R_{k}$-term gives an operator that changes the occupation numbers by two. However, since the occupation numbers of $\left|\Psi_{i}\right\rangle$ and $\left|\Psi_{f}\right\rangle$ differ by exactly one, $N_{L,f}-N_{L,i}=1$, the resulting matrix element vanishes:

[TABLE]

Therefore we only need to consider commutators involving the diagonal elements of $\mathcal{H}_{0}$ in (4.2.1). The commutator gives

[TABLE]

This is a sum of four terms which all, following from the periodicity of $f_{\rm R}(\{k_{i}\},t)$ and $e_{11,k}(t)$ [$e_{22,k}(t)$], have a prefactor that is periodic in all $k_{i}$ and oscillates with frequency $\omega$. In the rotating frame we showed that $|E_{11,k}|,|E_{22,k}|\leq\frac{W}{2}$, i.e. $|e_{11,k}|,|e_{22,k}|\leq\frac{1}{2}$. Hence the magnitude of each term that is generated by the commutator in Eq. (41) is bounded by $UW/2\Delta$, a reduction by a factor of $W/2\Delta$ as compared to the original magnitude of $\mathcal{V}_{x}^{\rm RR\to RL}$. The same argument applies to the commutator of $\mathcal{V}_{x}^{\rm RL\to LL}$ with $\mathcal{H}_{0}$.

Commutator with band-conserving interaction:

The commutator of the quasi-local operator $\mathcal{V}_{x}^{\rm RR\to RL}(t)$ with the intraband interaction operator from Eq. (12) contains commutations with $\mathcal{V}_{intra}^{\rm RR}(t)$, $\mathcal{V}_{intra}^{\rm LR}(t)$, and $\mathcal{V}_{intra}^{\rm LL}(t)$. As an example, the commutator of $V^{\rm RR\to RL}_{x}(t)$ with $\mathcal{V}_{intra}^{\rm RR}(t)$ in Eq. (13) is given by

[TABLE]

After some algebra (see appendix C.1 for details), the commutator in Eq. (4.2.3) is reduced to a sum of six terms, each involving a string of six creation/annihilation operators (i.e., corresponding to three-particle interaction). One example term is

[TABLE]

where

[TABLE]

This resulting term $h(\{k_{i}\},t)$ is periodic in all $k_{i}$, as is $f_{\rm 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 the resulting operator is a sum of quasi-local operators. The commutators with $\mathcal{V}_{intra}^{\rm LR}(t)$ and $\mathcal{V}_{intra}^{\rm LL}(t)$ yield similar terms as in Eq. (43). Compared with the magnitude of the four-point interaction terms in $\mathcal{V}_{x}$, the magnitudes of the six-point interaction terms generated by the commutator in Eq. (4.2.3) are suppressed by a factor $U/\Delta$. In appendix C.1 we show that overall there are 24 such terms. Additionally there are 24 more similar terms coming from the commutator of $\mathcal{V}_{x}^{\rm RL\to LL}$ with $\mathcal{V}_{intra}$.

Time-derivative of the operator:

The third term in Eq. (39) contains the explicit time-derivative of the operator $\mathcal{V}_{x}(t)$ in Eq. (28). The time-dependent parts therein are the Floquet single-particle creation and annihilation operators in the rotated frame, $R_{k}^{\dagger}(t),L_{k}^{\dagger}(t),R_{k}(t),L_{k}(t)$, and the functions $\alpha_{k}(t),\beta_{k}(t)$. In appendix C.2 we show that the derivative of $\mathcal{V}_{x}^{\rm RR\to RL}(t)$ is a sum of four terms of the form

[TABLE]

where in all terms $G(\{k_{i}\},t)$ is periodic in all $k_{i}$. The function $G(\{k_{i}\},t)$ is a product of functions $\alpha_{k}(t),\beta_{k}(t)$ and their derivatives, hence it is periodic in $t$ with frequency $\omega$. We further show in appendix C.2 that its Fourier transformation with respect to time has the largest coefficients at small frequencies much less than $\Delta$, and the Fourier coefficients decay exponentially beyond that regime. We will return to this fact later. Note that Eq. (45) contains the same single-particle Floquet operators as $\mathcal{V}_{x}^{\rm RR\to RL}(t)$. A similar calculation can be done to obtain the time-derivative of $\mathcal{V}_{x}^{\rm RL\to LL}(t)$.

Iterated derivatives of the matrix element:

So far we discussed the effects of taking one derivative $i\partial_{t}$ of the matrix element in Eq. (39). To calculate the norm of the operator $\mathcal{O}_{x}^{M}(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 $\mathcal{V}_{x}(t)$ in Eq. (28), in the sense that they contain a string of creation and annihilation operators $R_{k}^{\dagger},L_{k}^{\dagger},R_{k},L_{k}$, they overall move one particle from the $R$-band to the $L$-band, and their prefactors are periodic functions in all $k_{i}$.

When we iterate the derivative, in each step we start with a sum of terms like $\mathcal{V}_{x}(t)$ in Eq. (28). From each term we obtain several new terms akin to equations (41), (43) and (45). Starting from one term with $n$ single-particle Floquet state creation/annihilation operators $R_{k}^{\dagger},L_{k}^{\dagger},R_{k},L_{k}$, the three different contributions from Eq. (39) give:

• •

Commutator with $\mathcal{H}_{0}$: $n$ terms with $n$ creation/annihilation operators, each multiplied by an additional factor $W/\Delta$,

• •

Commutator with $\mathcal{V}_{intra}$: $6n$ terms with $n+2$ creation/annihilation operators, each multiplied by an additional factor $U/\Delta$,

• •

Derivative $\Delta^{-1}i\partial_{t}$: $n$ terms with $n$ creation/annihilation operators, and each multiplied by an additional factor $\omega/\Delta$

Therefore the total number of terms created by one term with $n$ single-particle operators is $n+6n+n=8n$. This process is visualized in Fig. 4. In the $m$-th step we obtain terms involving strings with different numbers of single-particle operators: the commutator with the intraband interaction $\mathcal{V}_{intra}(t)$ increases that number $n$ by 2, while the commutator with $\mathcal{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 $\mathcal{N}$ of terms in $\mathcal{O}_{x}^{M}(t)$ we can replace $n$ in the $m$-th iteration by the maximal possible length of an operator string generated after $m-1$ iterations of the derivative, $n_{max}(m)=2m+2$. The iteration begins with two quartic terms $\mathcal{V}_{x}^{\rm RR\to RL}$ and $\mathcal{V}_{x}^{\rm RL\to LL}$, each involving $n=4$ creation/annihilation operators. Hence the total number $\mathcal{N}$ of terms after $M$ steps is bounded by

[TABLE]

In the last step we use the bound $(M+1)!\leq\frac{5}{2}M^{M}$, which is valid for all $M>0$.

Operator norm after iterated derivatives:

Now we are in the position to write down the explicit form of the operator $\mathcal{O}_{x}^{M}(t)$ in Eq. (35) and calculate a bound for the norm $\|\mathcal{O}_{x}^{M}(t)\|$ in Eq. (38). After $M$ iterations of taking derivatives $\mathcal{O}_{x}^{M}(t)$ is a sum of $\mathcal{N}$ terms. Each contains a prefactor $U^{q}W^{p}\omega^{M-p-q}/\Delta^{M}$. The integers $p,q$ depend on the path that was taken in Fig. 4, i.e. $p$, $q$ count the number of times the commutator is taken with $\mathcal{H}_{0}$ and $\mathcal{V}_{intra}$, respectively. Note that each path leads to multiple terms. We use the label $\eta$ to index all of the terms arising from all of the different paths. Each term contains an additional dimensionless prefactor $\mathcal{F}_{\eta}(\{k_{i}\},t)$ which is periodic in all $k_{i}$. Hence we write the operator $\|\mathcal{O}_{x}^{M}(t)\|$ in the form

[TABLE]

where in the above sum, $p$ and $q$ are determined by the path corresponding to each term $\eta$. We discuss the properties of $\mathcal{F}_{\eta}(\{k_{i}\},t)$ in detail in appendix C.3, and show that the norm of this operator is bounded by

[TABLE]

The parameter $\mathcal{E}$ is the largest of the renormalized energy scales

[TABLE]

Here $K,K_{U},K_{E},K_{\omega}$ are dimensionless constants. They depend on the properties of the Floquet single-particle states. In general $K,K_{U},K_{E}$ are of order $\mathcal{O}(1)$ for a Hamiltonian which is a sum of local terms. For $\omega>W$ the constant $K_{\omega}$ also is of order 1, $K_{\omega}=\mathcal{O}(1)$, while for $W>\omega$ it goes as $K_{\omega}=\mathcal{O}(\frac{W}{\omega})$. In either case the resulting energy scale $\mathcal{E}$ in Eq. (49) is proportional to one of the small energy scales $U,W,\omega$. In appendix C.3 we derive explicit expressions for these constants.

So far $M$ was a free parameter; Eq. (32) holds for all $M$. Hence we can freely choose this parameter. By differentiating we see that the tightest bound is reached for $M=\frac{\mu\Delta}{\mathcal{E}}$ with $\mu=(16e)^{-1}$. With this optimal choice of $M$, Eq. (48) becomes

[TABLE]

Our main result from this section is the exponential bound on the operator norm in Eq. (50). The main message here is that this norm is exponentially small if $\Delta\gg W,U,\omega$. To obtain the final result for the excitation rate from the nearby terms we apply the result (50) to Eq. (37):

[TABLE]

Together with the exponentially suppressed contribution from the distant terms, Eq. (4.1), 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 $\Gamma=\Gamma_{near}+\Gamma_{dist}$. For times $T\ll\mathcal{T}\ll r^{*}/v$ the first term in Eq. (51) is the dominant contribution, and the effect from the distant terms in Eq. (4.1) is negligible. In thermodynamically large systems $r^{*}$ can be chosen sufficiently large to have a long window during which $T\ll\mathcal{T}\ll r^{*}/v$. In this time interval the rate of creating excitations is exponentially small. Hence, during this long time window a prethermal quasi-steady state persists.

5 Summary and discussion

In this paper we showed that in slowly driven systems, assuming the hierarchy of energy scales $\Delta\gg U,W,\omega$, the rate at which particles transition from one band to the other is exponentially small in the large ratio $\Delta/\mathcal{E}$ where $\mathcal{E}=\text{max}(U,W,\omega)$. The intraband thermalization rate is not suppressed by $\Delta/\mathcal{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 Floquet 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. (51) 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 $\Delta$ is the largest energy scale in comparison to the drive frequency $\omega$, 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 Floquet eigenstates; in appendix C.3 we provide an explicit calculation of this factor. Generally, a small localization length of the Floquet-Wannier functions leads to a large factor in the exponent and strong suppression of the excitation rate. This localization length also depends on the ratios $\Delta/W$ and $\Delta/\omega$. Therefore there are subleading corrections in the final result in Eq. (51) which enhance the suppression of the transition rate as $\Delta/\mathcal{E}$ increases.

We obtained an upper bound for the rate $\Gamma$ of particles being transferred from the $R$- to the $L$-band. If the $L$-band in the initial state $\left|\Psi_{i}\right\rangle$ 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 bound we calculated. 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 D. This tighter bound is similar to Eq. (51) with the replacement of $\mu\to\tilde{\mu}=2\mu$. 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 [57]. Hence we may expect that the tighter bound also holds when the $L$-band is partially filled.

In appendix E 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 E 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. Third we generalize to systems with an arbitrary number of bands. In multi-band systems we only require that the bands fall into groups which are separated by a large bandgap $\Delta$. All in all we can apply our results to a wide class of slowly-driven systems. Another interesting question is whether a similar bound holds for a system of bosons, or whether bunching effects lead to qualitatively different behavior.

Following the original argument by Thouless [53], half filling is required in order to obtain a gapped state 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 [62]. Here we show that in the presence of inter-particle interactions, there exist exponentially long-lived chiral prethermal states at arbitrary filling fractions (see also [57]). Thus our results opens new possibilities for observing topological Floquet bandstructures in quantum many body systems.

Acknowledgements

N.L., T.G. and E.B. acknowledge support from the European Research Council (ERC) under the European Union Horizon 2020 Research and Innovation Programme (Grant Agreement No. 639172). T.G. was in part supported by an Aly Kaufman Fellowship at the Technion. T.G. acknowledges funding from the Institute of Science and Technology (IST) Austria, and from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 754411. N.L. acknowledges support from the People Programme (Marie Curie Actions) of the European Unions Seventh Framework 546 Programme (FP7/20072013), under REA Grant Agreement No. 631696, and by the Israeli Center of Research Excellence (I-CORE) Circle of Light funded by the Israel Science Foundation (Grant No. 1802/12). M.R. gratefully acknowledges the support of the European Research Council (ERC) under the European Union Horizon 2020 Research and Innovation Programme (Grant Agreement No. 678862). M.R. acknowledges the support of the Villum Foundation. M.R. and E.B. acknowledge support from CRC 183 of the Deutsche Forschungsgemeinschaft.

Appendix A Lattice model for Thouless’ pump

To obtain the plots of the band structure we use the same model for a 1-dimensional topological charge pump as in Ref. [57]. The tight-binding Hamiltonian is

[TABLE]

This is a tight-binding Hamiltonian with alternating nearest-neighbor hopping $J(t)$ and $J^{\prime}(t)$, and alternating on-site potential $\pm\mu(t)$. For the figures we set $J(t)=J_{0}+\delta J(t)$ and $J^{\prime}(t)=J_{0}-\delta J(t)$ with $J_{0}=1$ and $\delta J(t)=1.2\sin(\omega t)$. The on-site potential is chosen as $\mu(t)=2.5\cos(\omega t)$. In Fig. 2 the frequency is $\omega=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.

The full form of Eq. (36) is

[TABLE]

The boundary terms in the first line give a function which is oscillating in $\mathcal{T}$ but not increasing. The prefactor $1/\mathcal{T}$ ensures that this term vanishes at long times. Every derivative $i\partial_{t}$ of the matrix element yields an energy $\omega,U,W$. Together with the small factor $1/\Delta$ we see that these terms give a small offset in which the main contribution is the $m=0$ term. Full justification of this qualitative argument follows from the derivation in section 4.2.3. Hence for the boundary terms we estimate

[TABLE]

In the first step we apply a triangle inequality, so that we only obtain terms with the same time. In the second step we use $\sum_{f}\left|\Psi_{f}(t)\right\rangle\left\langle\Psi_{f}(t)\right|=\mathds{1}$ within the subspace of fixed particle numbers $N_{R,f}$ and $N_{L,f}$, and the norm is $\|{\mathcal{V}_{x}}^{\dagger}(t)\mathcal{V}_{x}(t)\|=U^{2}$ as in section 4.1. Therefore the offset to the excitation probability goes, in leading order, as $(U/\Delta)^{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. (53). Here we want to show that this part is independent of time $\mathcal{T}$, i.e. the system has a constant heating rate. To simplify notation we define the operator $\mathcal{O}_{x}^{M}(t)$ implicitely by (cf. Eq. (35) in the main text)

[TABLE]

This definition is not unique, however in the following we only need to consider one operator $\mathcal{O}_{x}^{M}(t)$ which fulfills Eq. (55). From the periodicity of $\mathcal{V}_{x}(t)$ follows that the operator $\mathcal{O}_{x}^{M}(t)$ is periodic in time, and we decompose it into its Fourier modes ${\mathcal{O}_{x}^{M}}^{m}$. Similarly $\left|\Psi_{i/f}(t)\right\rangle$ are Floquet eigenstates of the Hamiltonian $\mathcal{H}_{0}+\mathcal{V}_{intra}$ with quasienergies $\epsilon_{i/f}$, so we can write them in their Fourier modes,

[TABLE]

Rewriting the second line in Eq. (53) in terms of these Fourier modes gives

[TABLE]

In the last line we set $E_{l,m,n}=\epsilon_{f}-\epsilon_{i}-\Delta+(l-m-n)\omega$, and $\operatorname{sinc}(x)=\frac{\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\in[-\pi,\pi]$, which are the two central zeroes. Hence we only need to consider $E_{l,m,n}\in[-\frac{2\pi}{\mathcal{T}},\frac{2\pi}{\mathcal{T}}]$. For $\mathcal{T}>2T$ not for all quasi-energies $\epsilon_{f}\in[0,\omega]$ 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, $\rho(\epsilon_{f})=\rho$ for $\epsilon_{f}\in[0,\omega]$.[27] Hence we can convert the sum over final states into an integral, $\sum_{f}\to\int d\epsilon_{f}\rho(\epsilon_{f})$, with a constant density of states $\rho(\epsilon_{f})=\rho$ for $\epsilon_{f}\in[0,\omega]$. Thus the fraction of final states that we need to consider is $\frac{2T}{\mathcal{T}}$ of all states, and we continue:

[TABLE]

Here we used translational invariance, i.e. the sum over all sites becomes $\sum_{x=1}^{L}\to L$ and we consider the operator at one example point $x_{0}$. Because the states are indistinguishable the quasi-local operator $\mathcal{O}_{x}^{M}(t)$ acts on all of them in a similar way, and the matrix element does not depend on the particular choice of $\Psi_{f}$. To obtain the matrix element and the density of states we calculate in two ways:

[TABLE]

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) we arrive at

[TABLE]

Appendix C Details for the calculations in section 4.2.3

In this appendix we summarize the technical details for the calculations in section 4.2.3. These are required to derive the operator $\mathcal{O}_{x}^{M}(t)$ in Eq. (4.2.3), and subsequently derive a bound on its norm.

C.1 Commutator of the interaction with the intraband operator

In this section we show the explicit calculation of the commutator of $\mathcal{V}_{x}(t)$ with the intraband operator $\mathcal{V}_{intra}(t)$ in section 4.2.3. We start with Eq. (4.2.3), the commutator of $\mathcal{V}_{x}^{\rm RR\to RL}$ with $\mathcal{V}_{intra}^{\rm RR}$. It contains two terms of eight Floquet single-particle operators each, in different order. By using the canonical fermionic anticommutation relations of the Floquet single-particle operators,

[TABLE]

we bring the operators in the same order to cancel the terms with eight operators. This requires six non-trivial commutations according to (61). In the end we arrive at

[TABLE]

The result are six terms with six operators each. This we apply to Eq. (4.2.3) and calculate for one example term:

[TABLE]

In the first step we use the definitions of the functions $f_{\rm R}(\{k_{i}\},t)$ and $g_{\rm RR}(\{q_{i}\},t)$ from equations (28) and (12). In the second step we solve the sum over $x^{\prime}$ which yields a $\delta$-function over the momenta $q_{i}$. In the third step we solve the sums over $q_{1}$ and $k_{2}$ by using the $\delta$-functions. All six terms from Eq. (62) have a similar shape.

Next we turn to the commutator with $\mathcal{V}_{intra}^{\rm LR}$. Commuting the single-particle Floquet operators gives

[TABLE]

These are four additional terms, however the prefactor $g_{\rm LR}$ in Eq. (14) is a sum of 4 terms. Therefore this gives a total of 16 new terms which are all of similar structure as Eq. (63).

The last contribution is the commutator with $\mathcal{V}_{intra}^{\rm LL}$, for which we calculate

[TABLE]

These are 2 more terms with the same structure as Eq. (63). Hence in total the commutator of $\mathcal{V}_{x}^{\rm RR\to RL}$ with the intraband interaction $\mathcal{V}_{intra}$ generates 24 new terms. The same calculation with $\mathcal{V}_{x}^{\rm RL\to LL}$ yields the same result, the commutator generates 24 more terms.

C.2 Explicit time-derivative

This subsection deals with the explicit time-derivative of the operator $\mathcal{V}_{x}(t)$. To take the derivative of $\mathcal{V}_{x}^{\rm RR\to RL}(t)$ we need to calculate

[TABLE]

We use the product rule of the derivative to separate this into four terms, where in each we take the derivative of terms with one $k_{i}$ only. We have e.g. for $k_{4}$

[TABLE]

where we used the relation between Floquet eigenstates and sublattice basis, Eq. (26). These are products of two functions which by themselves are not periodic in time, but their products are:

[TABLE]

In the last step we expand this product into its temporal Fourier series. With its help we calculate the derivative in Eq. (67):

[TABLE]

We rewrite this in terms of the Floquet eigenstate operators in the rotating frame to arrive at

[TABLE]

The $L_{k_{4}}$ term in Eq. (C.2), combined with the operators $L_{k_{1}}^{\dagger}R_{k_{2}}^{\dagger}R_{k_{3}}$ in Eq. (66), yields an operator which preserves the band occupation numbers. Since the states $\left|\Psi_{i}\right\rangle$ and $\left|\Psi_{f}\right\rangle$ have different occupation numbers this matrix element is zero, similar to the argument in Eq. (40). Hence the term involving $L_{k_{4}}$ in Eq. (C.2) 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 $\mathcal{V}_{x}^{\rm RR\to RL}(t)$, with an additional prefactor $\frac{\omega}{\Delta}$. We perform analogous calculations for the terms with $k_{1},k_{2},k_{3}$ in Eq. (66), where each gives a new term involving the same single-particle operators. Overall we obtain four new terms which are all of the form

[TABLE]

where for all cases $G(\{k_{i}\},t)$ is periodic in all $k_{i}$. In the example from Eq. (C.2),

[TABLE]

The functions $\tilde{\alpha}_{k}(t),\tilde{\beta}_{k}(t)$ describing the single particle Floquet states in the original (non-rotating) frame, defined in Eq. (9), have most of their spectral weight in the Fourier modes with frequency $l\omega\approx\Delta/2$. In the rotating frame, the spectral weights of $\alpha_{k}(t),\beta_{k}(t)$ are shifted to the slow modes around $l\approx 0$ with decay for large $l$ as $e^{-A|l|^{3/2}}$ [57]. Therefore, the sum over the absolute value of all modes $l$ is finite. We will return to this sum in section C.3.

C.3 Bound on the operator norm

In this appendix we derive the final bound on the operator $\mathcal{O}_{x}^{M}(t)$ in section 4.2.3. We start with the expression in Eq. (4.2.3),

[TABLE]

$\mathcal{F}_{\eta}(\{k_{i}\},t)$ is a periodic function for all $k_{i}$ because it is a product of periodic functions $\alpha_{k}$, $\beta_{k}$, $e_{11,k}$, and $e_{22,k}$. That implies the Fourier coefficients in the expansion

[TABLE]

go to zero exponentially fast for large $|z_{i}|$. The same applies to the Fourier series of the other periodic functions,

[TABLE]

the coefficients go to zero at large $|y|$. We apply these expansion to Eq. (C.3) and rewrite the Floquet state operators in the site basis $\{A_{x_{i}},B_{x_{i}}\}$ according to definition (9):

[TABLE]

The momenta $k_{i}$ only appear in the exponent thus in the second step we can sum over all of them and obtain Kronecker-$\delta$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 $\overline{\mathcal{F}}_{\eta}$ or the $\overline{\alpha},\overline{\beta}$ are small in that case. This is a manifestation of locality of the operator $\mathcal{O}_{x}^{M}(t)$. The norm for each on-site operator is $\|A_{x_{i}}\|=\|B_{x_{i}}\|=1$, so we can estimate the norm of $\mathcal{O}_{x}^{M}(t)$ by taking absolute values of each individual factor. As a result the sums over $\{y_{i},z_{i}\}$ decouple:

[TABLE]

As we argued above, the absolute values of the coefficients $\overline{\alpha}_{y},\overline{\beta}_{y},\overline{e}_{y}$ in the Fourier series (75) go to zero exponentially fast for large $|y|$. Hence the sums over their absolute value are finite and we define the dimensionless constants

[TABLE]

Similarly we know that the coefficients $\overline{\mathcal{F}}_{\eta}(\{z_{i}\},t)$ decay exponentially at large $|z_{i}|$ and therefore the sum $\sum_{\{z_{i}\}}|\overline{\mathcal{F}}_{\eta}(\{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 $\mathcal{N}=2$ operators with $\mathcal{F}_{1}(\{k_{i}\},t)=f_{\rm R}(\{k_{i}\},t)=(\alpha_{k_{1}}\alpha_{k_{2}}^{*}+\beta_{k_{1}}\beta_{k_{2}}^{*})\alpha_{k_{3}}\beta_{k_{4}}e^{i\theta_{k_{1}}}$ and $\mathcal{F}_{2}(\{k_{i}\},t)=f_{\rm L}(\{k_{i}\},t)$. For now we focus on the terms from $f_{\rm R}(\{k_{i}\},t)$. The summation over $\{z_{i}\}$ is bounded by

[TABLE]

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 $\eta$ labels all these terms including which branch we follow during the iterated derivatives. In the following we calculate the change to the sum over the coefficients $\overline{\mathcal{F}}_{\eta}(\{z_{i}\})$ in Eq. (77) 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 $\mathcal{F}_{\eta}(\{k_{i}\},t)$ is of the form

[TABLE]

To express the Fourier coefficients $\overline{\mathcal{F}}_{\eta}(\{z_{i}\})$ we use the Fourier series in (75). From this we derive

[TABLE]

The momentum $k_{4}$ appears twice in the expression for $\mathcal{F}_{\eta}(\{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 absolute values:

[TABLE]

Compared with the bound in Eq. (79) this only differs by an additional factor $K_{E}=\sum_{z}|\overline{e}_{z}|$. So in every iteration of taking the commutator with the quadratic Hamiltonian the bound on the sum over coefficients $|\overline{\mathcal{F}}_{\eta}(\{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

[TABLE]

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

[TABLE]

and summation over all terms gives

[TABLE]

Compared with the bound in Eq. (79) there is an additional factor $K_{\alpha}^{2}K_{\beta}^{2}$. However, commutation with the intraband operator also increases the number of Floquet single-particle operators by 2. As we have seen in Eq. (77) these sums decouple and the contribution to the norm from each pair of Floquet operators $R_{k}^{\dagger}R_{k^{\prime}}$ is $(K_{\alpha}+K_{\beta})^{2}$. Hence in every iteration of commuting with the intraband operator the bound on the sum over coefficients is multiplied by $K_{\alpha}^{2}K_{\beta}^{2}(K_{\alpha}+K_{\beta})^{2}$.

Time derivative of the operator:

The function $\mathcal{F}_{\eta}$ for one example term when taking the time derivative is calculated in Eq. (C.2):

[TABLE]

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

[TABLE]

The summation over all $\{z_{i}\}$ gives

[TABLE]

As discussed in the main text, the product $\alpha_{k}(t)\beta_{k}^{*}(t)$ is periodic in $t$ (Eq. (C.2)). Therefore the absolute value of the coefficients in the Fourier series over time $t$, $(\alpha_{k}\beta_{k}^{*})^{(l)}$, goes to zero exponentially fast at large $l$. Starting from a Zener tunneling problem, [57] shows that these coefficients decay as $e^{-A|l|^{3/2}}$. Therefore the sum in Eq. (C.3) is finite and with the definition of a constant

[TABLE]

we get

[TABLE]

This differs from the result in Eq. (79) by a factor $\kappa(K_{\alpha}^{2}+K_{\beta}^{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

• •

the initial factor $K\equiv K_{\alpha}K_{\beta}(K_{\alpha}^{2}+K_{\beta}^{2})(K_{\alpha}+K_{\beta})^{4}$ from Eq. (79) and the initial 4 Floquet operators as an overall prefactor.

• •

a factor $K_{E}^{p}$ from $p$ commutators with the quadratic Hamiltonian.

• •

a factor $K_{U}^{q}\equiv K_{\alpha}^{2q}K_{\beta}^{2q}(K_{\alpha}+K_{\beta})^{2q}$ from $q$ commutators with the intraband operator.

• •

a factor $K_{\omega}^{M-p-q}\equiv\kappa^{M-p-q}(K_{\alpha}^{2}+K_{\beta}^{2})^{M-p-q}$ from $M-p-q$ applications of the time-derivative to the operator.

Now we are ready to apply this to the norm in Eq. (77):

[TABLE]

In the last step we defined the energy scale $\mathcal{E}$ as the maximum of all the involved renormalized energy scales,

[TABLE]

The expression (91) holds for all $M$ and we can apply it to Eq. (4.2.3) in the main text.

Appendix D Initial state with empty L-band

In section 5 we mention that we obtain a tighter bound if in the initial state $\left|\Psi_{i}\right\rangle$ the $L$-band is empty. This is the case, up to corrections of order $\mathcal{O}(\omega/\Delta)^{2}$ [69], at short times $\mathcal{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. (46). First we can neglect the intraband interaction that requires an $L$-particle in the initial state, $\mathcal{V}_{x}^{\rm RL\to LL}$. This operator can’t act on the state $\left|\Psi_{i}\right\rangle$. Secondly, when calculating the commutator of $\mathcal{V}_{x}$ with the intraband term $\mathcal{V}_{intra}$ in section 4.2.3 and appendix C.1 the counting of terms changes. The commutator $[\mathcal{V}_{x}^{\rm RR\to RL},\mathcal{V}_{intra}^{\rm RR}]$ is unchanged, from there we obtain six new terms. Next we commute with $\mathcal{V}_{intra}^{\rm LR}$ defined in Eq. (14). The matrix element, reduced to the operators, is

[TABLE]

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}}^{\dagger}$ creates the only particle in the $L$-band, hence $L_{q_{3}}$ has to act on this particle which gives $\delta_{k_{1},q_{3}}$. With this we arrive at the same shape as in Eq. (63), however the prefactor function $g_{\rm LR}(\{q_{i}\},t)$ only contains four pieces, hence this gives four additional terms similar to (63). The commutator with $\mathcal{V}_{intra}^{\rm LL}$ does not yield any contribution because $\mathcal{V}_{x}\left|\Psi_{i}\right\rangle$ has only one particle in the $L$-band and $\mathcal{V}_{intra}^{\rm LL}$ can only act on states with at least 2 $L$-particles. Thus we overall obtain ten new terms. When iterating these steps the commutator with $\mathcal{V}_{intra}^{\rm LR}$ always only gives four terms, independent of the number $n$ of single-particle Floquet operators. The commutator with $\mathcal{V}_{intra}^{\rm RR}$ 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. (46). 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 by

[TABLE]

This is in comparison to the bound $\mathcal{N}\leq 2(16M)^{M}$ in the main text. This changes the final result from Eq. (51) to

[TABLE]

with $\tilde{\mu}=(8e)^{-1}=2\mu$. Hence the overall prefactor is divided by 4, and the factor $\tilde{\mu}$ in the exponent is doubled which leads to stronger suppression.

Appendix E Generalization to arbitrary dimensions and multi-band systems

In this appendix we generalize the calculation of the excitation rate in section 4.2 to multi-band systems in any number of dimensions, and to interactions that do not only act within a single unit cell. We consider a non-interacting Hamiltonian with $N$ bands in the basis of Floquet single-particle eigenstates $\tilde{\psi}_{\nu}(\mathbb{k})$,

[TABLE]

Here and in the following we use Einstein summation for greek indices with $\nu,\tau=1,...,N$ running over all bands, and we use bold font $\mathbb{k}$ for vectors in $d$ dimensions. The Floquet states are related to the sublattice basis by a unitary transformation akin to Eq. (9),

[TABLE]

where $\tilde{\alpha}(\mathbb{k},t)$ is a unitary $N\times N$ matrix. The eigenstates and the matrix elements in the Hamiltonian follow the same periodicity properties as in Eq. (9) in the main text,

[TABLE]

for any inverse lattice vector $\mathbb{G}$. We assume the instantaneous eigenstates split into two groups of bands which are separated by a large bandgap $\Delta$. In analogy with the main text we define the two groups as $L$-and $R$-bands, respectively, i.e. $\tilde{\psi}=(\tilde{R},\tilde{L})$. The unitary transformation for the single-particle operators in this notation is

[TABLE]

In the following the indices $\nu$ and $\tau$ always run over all bands, while $\rho$ indexes the $R$- and $\lambda$ the $L$-bands only, respectively. A generic two-particle interaction between all bands can be written as

[TABLE]

We assume a quasi-local interaction which means that $|\chi_{\nu_{1},\nu_{2}}^{\nu_{3},\nu_{4}}(\{\mathbb{d}_{i}\})|$ decays exponentially at large $\|\mathbb{d}_{i}\|$, and we choose $\chi$ and $U$ so that $\sum_{\{\mathbb{d}_{i}\},\{\nu_{i}\}}|\chi_{\nu_{1},\nu_{2}}^{\nu_{3},\nu_{4}}(\{\mathbb{d}_{i}\})|=1$. The intraband interaction does not change the population of the $R$-bands, but as in the main text it splits into parts which act on two particles on the $R$-bands, two particles in the $L$-bands, or one particle in each group of bands. In momentum space the $RR$-term is

[TABLE]

with the implicit definition of $(\tilde{g}_{\rm RR})_{\rho_{1},\rho_{2}}^{\rho_{3},\rho_{4}}(\{\mathbb{k}_{i}\},t)$. Similarly we obtain the $LL$- and the $LR$-terms where the latter has a prefactor $\tilde{g}_{\rm LR}$ which consists of 4 terms, as in Eq. (14) in the main text. $V$ is the ($d$-dimensional) volume of the system.

The interband part of the interaction which excites one particle from the $R$- to the $L$-bands splits into two pieces, one part $\tilde{\mathcal{V}}_{+}^{\rm RR\to RL}(t)$ that acts on two particles in the $R$-bands and one part $\tilde{\mathcal{V}}_{+}^{\rm RL\to LL}(t)$ that acts on one particle in the $R$- and $L$-bands each. The first part becomes

[TABLE]

and similarly for $\tilde{\mathcal{V}}_{+}^{\rm RL\to LL}(t)$. We want to calculate the rate at which this interband interaction changes the particle density in the $L$-bands $n_{L}=N_{L}/V$:

[TABLE]

We start by separating nearby and distant terms as in section 4.1, determined by a critical distance $r^{*}$:

[TABLE]

After changing to the Heisenberg picture and completing the commutator as in Eq. (17) we apply the Lieb-Robinson bounds and obtain

[TABLE]

Here $S_{d}$ is the surface of the $d$-dimensional unit sphere, i.e. $S_{1}=2,S_{2}=2\pi,S_{3}=4\pi$. It is worth noting that the Lieb-Robinson velocity $v$ also depends on the localization length of the intraband interaction $\tilde{\mathcal{V}}_{\mathbb{x}}(t)$ which includes the localization length of $\chi_{\nu_{1},\nu_{2}}^{\nu_{3},\nu_{4}}(\{\mathbb{d}_{i}\})$ in (100).

Turning to the nearby terms we apply a rotating frame transformation akin to section 4.2.1:

[TABLE]

This implies for the entries in the unitary matrix $\alpha$ that

[TABLE]

In the rotated Hamiltonian

[TABLE]

all matrix elements in the diagonal blocks shift to a narrow interval of width $W$ (cf. Eq. (4.2.1)),

[TABLE]

In analogy with Eq. (25) we define the dimensionless quantities $e_{\nu,\tau}(\mathbb{k},t)\equiv E_{\nu,\tau}(\mathbb{k},t)/W$ to write the Hamiltonian as

[TABLE]

After applying a triangle inequality similar to Eq. (22) the excitation probability in the rotated frame becomes

[TABLE]

Here $V^{*}$ is the volume of a ball of radius $r^{*}$, and in the second step we introduced derivatives $e^{i\Delta t}=(-i\partial_{t})^{M}\Delta^{-M}e^{i\Delta t}$ and performed integration by parts. From there we derive a similar result for the excitation probability as in Eq. (37),

[TABLE]

in which the rate of continued excitations only depends on the norm of the operator $\|\mathcal{O}_{\mathbb{x}}^{M}(t)\|$ which is implicitely defined by

[TABLE]

Note that, as Eq. (35) in the main text, Eq. (113) does not uniquely define the operator $\|\mathcal{O}_{\mathbb{x}}^{M}(t)\|$, but we only need to consider one operator which fulfills the condition (113).

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 equation (39),

[TABLE]

From now on we only consider the operator $\mathcal{V}_{\mathbb{x}}^{\rm RR\to RL}(t)$, the contribution from $\mathcal{V}_{\mathbb{x}}^{\rm RL\to LL}(t)$ is calculated analogously.

First we consider the commutator of $\mathcal{V}_{\mathbb{x}}^{\rm RR\to RL}(t)$ with the quadratic Hamiltonian. The elements in the off-diagonal blocks change the occupation of the $L$-and $R$-band, therefore

[TABLE]

Hence we only need to account for the block on the diagonal of $\mathcal{H}_{0}$ in (4.2.1) with $R^{\dagger}R$ and $L^{\dagger}L$. The commutator gives

[TABLE]

As in the case with two bands this gives 4 new terms, however each of them contains a sum over all $L$- or all $R$-bands.

Next we calculate the commutator of the interaction with the intraband operator $\mathcal{V}_{intra}(t)$. For example, in the multi-band system the commutator with the $RR$-part $\mathcal{V}_{intra}^{\rm RR}(t)$ takes the form

[TABLE]

where the $\gamma_{i}$ index the $R$-bands only. As in the main text, the commutator with the full intraband operator gives 24 terms of the form

[TABLE]

where the prefactor function for one example term is

[TABLE]

To calculate the explicit time derivative of $\mathcal{V}_{\mathbb{x}}^{\rm RR\to RL}(t)$ we separate the time-dependent terms by their momenta $\mathbb{k}_{i}$, as in Eq. (67). From Eq. (E) we see that the functions $\alpha_{\nu,\tau}(\mathbb{k},t)$ are not periodic in time, however as in the case of two bands in Eq. (C.2) the relevant products are. We calculate in a similar manner:

[TABLE]

The part with the operator $L_{\lambda}(\mathbb{k}_{4},t)$ reduces the number of particles in the $L$-band, therefore it does not connect to the final state in the matrix. We can drop this part and only keep the operator $R_{\rho}(\mathbb{k}_{4},t)$. Performing this calculation for all four momenta $\mathbb{k}_{i}$ gives four new terms of the form

[TABLE]

where in the example from Eq. (E) the resulting prefactor function is

[TABLE]

From this discussion we see that each individual step gives several terms which have the same shape as the initial term $\mathcal{V}_{\mathbb{x}}(t)$. Therefore we can easily iterate taking derivatives. Furthermore the number of terms that are generated does not change compared to the special case in the main text. Thus the total number $\mathcal{N}$ of final terms after taking $M$ derivatives is identical to the result in Eq. (46) in the main text, including the terms from $\mathcal{V}_{\mathbb{x}}^{\rm RR\to RL}(t)$ and $\mathcal{V}_{\mathbb{x}}^{\rm RL\to LL}(t)$ gives

[TABLE]

note that here each term contains sums over all bands, however this does not affect the rest of the calculation.

Similarly to Eq. (4.2.3) in the main text we can write the operator $\mathcal{O}_{\mathbb{x}}^{M}(t)$ in Eq. (113) as

[TABLE]

We rewrite the operators $R_{\rho}(\mathbb{k},t),\,L_{\lambda}(\mathbb{k},t)$ in the sublattice basis according to Eq. (99). The functions $\alpha_{\nu,\nu^{\prime}}(\mathbb{k})$ and the functions $\mathcal{F}_{\eta,\lambda_{1},\{\rho_{i}\}}(\{\mathbb{k}_{i}\},t)$ are periodic in all momenta $\mathbb{k}_{i}$. We write them in their Fourier series akin to equations (75) and (74). Summing over all momenta and using $\|A_{\nu}(\mathbb{x})\|=1$ for the norm of the on-site operators we find

[TABLE]

The absolute values of $\overline{\mathcal{F}}$ and $\overline{\alpha}$ decay exponentially at large $\|\mathbb{y}_{i},\mathbb{z}_{i}\|$ because of the $\mathbb{k}$-periodicity of $\mathcal{F}$ and $\alpha$. Thus the sums over $\mathbb{y}_{i},\mathbb{z}_{i}$ in Eq. (125) converge to a finite value. In the initial term with $M=0$ derivatives we have $\mathcal{F}_{1}=f_{\rm R}$ and $\mathcal{F}_{2}=f_{\rm L}$. We bound the constant by

[TABLE]

As in appendix C.3 we now consider the changes to this constant $K$ that come from taking one commutator or one derivative. In all cases the norm of $\mathcal{F}$ only obtains an additional factor. The commutator with the quadratic Hamiltonian $\mathcal{H}_{0}$ changes

[TABLE]

in the expression of $|\overline{\mathcal{F}}_{\eta}|$. We define the constant $K_{E}$ analog to the definition in Eq. (C.3) so that

[TABLE]

The change to Eq. (126) due to the commutator with the intraband interaction $\mathcal{V}_{intra}$ is

[TABLE]

After the smoke clears this merely multiplies the constant $K$ by a numerical factor. We define this factor as $K_{U}$, in analogy with $K_{E}$ in Eq. (128). The explicit time-derivative of the operator yields the change

[TABLE]

We define the additional contribution to the constant in Eq. (126) as $K_{\omega}$. We want to stress again that even though these constants all contain sums over all coefficients in the Fourier series and over all bands, due to the periodicity of the functions the Fourier coefficients go to zero exponentially fast and the constants are all finite. We obtain the same result for the operator norm as in Eq. (50) in the main text,

[TABLE]

where the constant $K$ is defined in Eq. (126). In the last step we again choose $M=\frac{\mu\Delta}{\mathcal{E}}$ with $\mu=(16e)^{-1}$, and we use the renormalized energy scales $\mathcal{E}=\text{max}\{K_{U}U,K_{E}W,K_{\omega}\omega\}$. This is identical to the previous result, only the constants $K_{U,E,\omega}$ that renormalize the energy scales change. Applying this to the excitation probability in Eq. (E) and adding the contribution from the distant terms in Eq. (E) gives

[TABLE]

This final result is the direct generalization for the special case of a two-band system in one dimension, Eq. (51). In any dimension the excitation rate grows linearly in time with a heating rate that is exponentially small in $\frac{\Delta}{\mathcal{E}}$. The rate of change of the particle densities in the two bands is independent of the system volume which is expected from the thermodynamic limit. The tighter bound for a system with initially empty $L$-bands in appendix D also generalizes straightforwardly to the multiband system.