Universal chiral quasi-steady states in periodically driven many-body systems

We investigate many-body dynamics in a one-dimensional interacting periodically driven system, based on a partially-filled version of Thouless's topologically quantized adiabatic pump. The corresponding single particle Floquet bands are chiral, with the Floquet spectrum realizing nontrivial cycles around the quasienergy Brillouin zone. For generic filling, with either bosons or fermions, the system is gapless; here the driving cannot be adiabatic and the system is expected to rapidly absorb energy from the driving field. We identify parameter regimes where scattering between Floquet bands of opposite chirality is exponentially suppressed, opening a long time window where the many-body dynamics separately conserves the occupations of the two chiral bands. Within this intermediate time regime we predict that the system reaches a chiral quasi-steady state. This state is universal in the sense that the current it carries is determined solely by the density of particles in each band and the topological winding numbers of the Floquet bands. This remarkable behavior, which holds for both bosons and fermions, may be readily studied experimentally in recently developed cold atom systems.


I. INTRODUCTION
Topological transport has garnered great attention in condensed matter physics, ever since the discovery of the quantized Hall effects [1,2]. Traditionally, robust features of transport have been linked to topological properties of the ground states of many-body systems. Recently, a new paradigm has emerged for altering the topological properties of band structures using external driving . These ideas have sparked a variety of proposals and initial experiments aimed at realizing so-called Floquet topological insulators in a range of solid state, atomic, and optical systems [24][25][26][27][28]. While the prospect of dynamically controlling band topology is exciting, establishing how and to what extent topological phenomena can be observed in such systems raises many crucial and fundamental questions about many-body dynamics in periodically driven systems.
The long-time behavior of periodically-driven manybody systems is a fascinating question of current investigation. Several recent theoretical works suggest that closed, interacting, driven many-body systems generically absorb energy indefinitely from the driving field, tending to infinite-temperature-like states in the long time limit [29][30][31]. In such a state, all correlations are trivial, indicating in particular that any topological features of the underlying Floquet spectrum are expected to be washed out. On the other hand, several interesting exceptions to the infinite temperature fate have been proposed [32][33][34]. Particularly, heating may be circumvented by local conservation laws, e.g., in integrable [35] or many-body localized systems [36][37][38][39][40][41][42][43].
A system's featureless fate at long times may not preclude it from exhibiting topological phenomena transiently, on timescales shorter than that of its unbounded heating. Indeed, under some conditions such as high fre- The dimerized hopping amplitudes and on-site potentials are changed adiabatically as depicted by the magenta curve. The origin δJ = δV = 0 is a degeneracy point. b) The Floquet spectrum for ω = 0.2J0 and λ = 1, see text below Eq. (1). The right (left) moving Floquet band is colored green (yellow). c) The period-averaged current vs. time obtained from a numerical simulation of a system of length L = 16 unit cells with N = 8 particles. The initial value of the current depends on the initial state, but after a few driving periods the system reaches a quasi-steady state featuring a current J ≈ ρ/T , where ρ = N/L is the density of particles. d) At long times, particle scattering from the right to the left moving band leads to a decay of the current. The parameters used for panels (c) and (d) are ω = 0.33J0, λ = 0.66 and U = 2J0. quency driving [44][45][46][47], exponentially-slow heating rates provide long time windows in which interesting "prethermal" behavior may be observed [47][48][49][50][51][52][53][54][55].
In this paper we show that a non-trivial long lived quasi-steady state can be stabilized in an interacting periodically driven system, in the low frequency limit. We focus our attention to one dimension (1d), where Thouless showed that cyclic adiabatic driving in filledband/gapped many-body systems leads to topologicallyquantized charge pumping [56]; this phenomenon was recently observed in cold atoms experiments [57,58]. An example of a time-dependent tight-binding model with this property is illustrated in Fig. 1a [see Eq. (1) below]. In the adiabatic limit, the system's single-particle Floquet spectrum exhibits a non-trivial winding of each band [5]: the quasi energy changes by ω ≡ 2π/T (where T is the driving period) as the crystal momentum changes from 0 to 2π/a, where a is the lattice constant [5,59]. A characteristic Floquet spectrum with such winding (or "chiral") bands is shown in Fig. 1b.
The average current carried by the system over one driving period is quantized when one of the Floquet bands is completely filled with fermions, and the other is completely empty. What happens when the bands are initially only partially filled, or if the system is comprised of bosons? Here the system lacks a many-body gap, and the evolution cannot be considered to be adiabatic. In the absence of interactions any value of the pumped charge per cycle is possible, depending on the details of the Hamiltonian and the initial state. However, in the interacting case we identify a parameter regime where heating naturally drives the system to a quasi-steady state with universal properties, reflecting the topological nature of the Floquet band structure.
We study dynamics in the situation where one of the chiral Floquet bands is initially partially filled, while the other is empty. We assume that the minimum gap between the two bands of the single particle Hamiltonian (minimized over all times and quasi-momenta), ∆, is much larger than max[ ω, U, W ], where U and W are the interaction strength and the maximum bandwidth, respectively (see below for a precise definition of W ). Based on an analysis of low and high order scattering rates, we argue that that (i) regardless of details, the system equilibrates on a short time scale τ intra to a chiral infinite-temperature-like state in which all momentum states in one of the chiral bands are equally populated, while the other band remains nearly empty; (ii) in this state, the average pumped charge per period is approximately wρ, where w is the winding number of the partially filled band (see below) and ρ is the density of particles [60]; and (iii) the quasi-steady state persists for a long time scale τ inter , that can be larger than τ intra by many orders of magnitude. At times t τ intra , the system relaxes to a state in which the values of all local observables are as in an infinite temperature state. In particular, the current tends to zero. This physical picture is supported by exact numerical simulations of finite systems, see Figs. 1c,d.
Model. -For concreteness, we consider a 1d lattice with two sites per unit cell, labeled A and B. The hopping matrix elements and on-site potentials are time dependent (see Fig. 1a). The lattice is populated by a finite density ρ of identical fermions or bosons per unit cell.
We write the Hamiltonian as H(t) = H 0 (t) + H int , where the single particle part H 0 (t) is given by Here,ĉ † j,A (ĉ † j,B ) creates a particle of type A (B) in unit cell j. To avoid later ambiguity, we use hats to denote creation and annihilation operators throughout this work. Below we take J(t) = J 0 +δJ(t), J (t) = J 0 −δJ(t), and V (t) = V 0 + δV (t), with δJ(t) = λJ 0 cos ωt and δV (t) = V 1 sin ωt. For simplicity we fix V 1 = 3λJ 0 in all simulations presented below.
We consider a local interaction The intra-unit-cell form of the interaction in Eq. (2) is convenient for the analysis below. However, we do not expect our conclusions to depend on this choice. The single particle Floquet spectrum corresponding to H 0 (t) is found by seeking solutions to the Schrödinger equation which satisfy [61]: |Ψ 1P (t) = e −iε1Pt |Φ 1P (t) , with |Φ 1P (t + T ) = |Φ 1P (t) . We decompose the periodic function |Φ 1P (t) in terms of an infinite set of (nonnormalized) discrete Fourier modes {|ϕ The full time-dependent evolution of |Ψ 1P (t) is specified by the quasienergy ε 1P and a vector of Fourier coefficients Throughout this work we choose the quasi-energies of all single particle states to lie within a fundamental Floquet-Brillouin zone, 0 ≤ ε 1P < 2π/T . The Fourier coefficients comprising |Ψ(t) are determined by an eigenvalue equation εϕ = H 0 ϕ, where the "extended Hamiltonian" H 0 is constructed from the Fourier decomposition of H 0 (t) (see Appendix A) and ϕ without a ket symbol stands for a column vector of Fourier modes as in Eq. (4). We use calligraphic symbols for matrices in the space of Fourier coefficients. For harmonic driving, H 0 (t) = H dc +Λe iωt +Λ † e −iωt , the matrix H takes a simple block tri-diagonal form in harmonic (m) space: (H 0 ) mm = (H dc + mω)δ mm + (Λδ m,m −1 + Λ † δ m,m +1 ). The single particle Floquet spectrum is shown in Fig. 1b for parameters specified in the caption. For each of the two bands, we assign a winding number w, such that the quasi-energy band winds by wω when the quasi-momentum k changes from 0 to 2π/a. In our case, the two bands have winding numbers w = +1 and w = −1, and we refer to them as the right-moving (R) and left-moving (L) bands, respectively. More generally, a non-trivial winding is achieved in the adiabatic limit when the curve (δJ(t), δV (t)) encircles the origin (as in Fig. 1a).

II. MANY-BODY DYNAMICS
We now turn to the many-body dynamics of this system. We consider the situation where the system is initialized with a finite density of particles in the net rightmoving (R) Floquet band, shown in green in Fig. 1b. The initial momenta of the particles are arbitrary.
To investigate the timescales for intra-band and interband scattering, we develop a perturbative analysis of the many-body dynamics of the system. The perturbation series is organized in terms of powers of the interaction strength U . Crucially, we work in a Floquet picture where the time-dependent driving is first taken into account exactly, to all orders in the driving. As above, we work in the extended space of Fourier coefficients, where the many-body Floquet eigenstates are described by the eigenvectors of the extended Hamiltonian, The extended Hamiltonian defines a (static) eigenvalue problem that yields the Fourier coefficients describing Floquet eigenstates. One may also use the extended Hamiltonian to generate an effective evolution in the extended space, in an auxiliary time variable τ , via i∂ τ ϕ(τ ) = Hϕ(τ ). For the stroboscopic times τ = nT , where n is an integer and T is the driving period of the original problem, the "evolved" vector of Fourier coefficients ϕ(nT ) precisely captures the state of the system in the physical Hilbert space at the corresponding time t = nT (see Appendix B). Using this mapping we obtain transition rates for the stroboscopic evolution by employing standard Green's function techniques to the auxiliary evolution problem in the extended space.
Our aim is to calculate the rate at which particles are scattered into the left-moving (L) band. For weak interactions and short times, it is natural to view this process in terms of a perturbation series in the interaction U. We express the auxiliary-time evolution of the Fourier vector ϕ(τ ) in terms of its Fourier transform, ϕ(τ ) = ∞ −∞ dτ e iΩτφ (Ω). In terms of the extended Green's function G 0 (Ω) = (Ω − H 0 + iδ) −1 and T-matrix where ϕ 0 is the Fourier vector corresponding to the "free" initial state in which all particles are initialized in singleparticle Floquet eigenstates in the right-moving (R) band of the non-interacting system (see Appendix C for details of the construction of ϕ 0 ).

A. Born approximation
As a first step, we investigate the scattering rates to leading order in U, i.e., in the Born approximation T (Ω) ≈ U. This approximation captures the leadingorder behavior of two-particle scattering, which one may expect to be relevant for weak interactions and low densities (see discussion below and Refs. [62][63][64]).
Within the Born approximation, the transition rate is given by Γ This is Fermi's golden rule, adapted for a Floquet system [65]. Here f labels all final "free" Floquet eigenstates, and {E f } are their eigenenergies (with respect to the extended Hamiltonian H 0 ). We break Γ into two parts, Γ = Γ intra +Γ inter , corresponding to intra-band and interband scattering, respectively. The latter scattering processes transfer one or more particles from the R to the L band.
Figures 2a,b show the two-particle scattering rates for a pair of bosons in the R band, with momenta {k 1 , k 2 }, to scatter either to states within the R band (processes we denote as RR→RR, Fig. 2a), or to states with one particle in the R band and one in the L band (RR→RL, Fig. 2b). The two-particle scattering rates are given by where ϕ 2P,0 and ϕ 2P,f are the wavefunctions of the initial and final two-particle Floquet states (see Appendix C), ∆v f is the difference of group velocities of the two outgoing particles, and the summation is over all final states that satisfy quasi-energy and quasi-momentum conservation. The factor 1/|∆v f | comes from the density of outgoing states.
Our key observation is that the inter-band scattering rates are strongly suppressed compared with the intraband ones. For the parameters chosen, the mean RR → RL scattering rate is down by a factor ∼ 10 −12 compared to the average RR → RR rate. The average RR → LL rate (not shown) is suppressed by a factor of ∼ 10 −10 relative to RR → RL. The bright lines visible in the interband scattering rate, Fig. 2b, are associated with single particle resonances; such resonances yield significant interband scattering rates only in exponentially small regions of phase space (see Appendix A).
The suppression of inter-band relative to intra-band scattering originates from the matrix element in Eq. (6). Since U is diagonal in Fourier harmonics, the matrix element for inter-band scattering is suppressed if the initial and final states have support in different regions in harmonic space. In the adiabatic limit ω ∆, this is indeed the case. Figure 2c shows two representative singleparticle Floquet wavefunctions with quasi-momenta k 1 , Intensity,  6), with forward scattering contribution removed. Here k1 and k2 are the momenta of two incident particles, and we take λ = 0.56, ω = 0.2J0, and U = 0.67J0. b) Same as above, for interband scattering in which one particle is scattered from the right to the left moving band. c) Single particle Floquet states. On the left, the shaded region spans the energies of the instantaneous bands as a function of time. The significant Fourier components |ϕ 1P | 2 of a Floquet state α with momentum k fall between the minimum and maximum values of the instantaneous energies E α,k (t). d) Fourier components |ϕ 2P | 2 of incoming and outgoing two particle Floquet states. For intraband scattering (bottom) the Fourier components overlap. For interband scattering (top), the overlap is strongly suppressed leading to a suppression of the interband scattering rate. k 2 , one from each band. The support of each state in harmonic space corresponds to the energy window spanned by the instantaneous energy, E α,k (t) (the eigenvalue of H 0 (t) for band α = R,L), with 0 ≤ t < T [66]. This energy window is bounded by W , which we define as Fig. 2c. Outside of this window, the Floquet wavefunctions decay rapidly. The separation of the Floquet states of the two bands in harmonic space can be derived by mapping the Floquet problem to a Zener tunneling problem in a weak electric field (see Appendix A). Figure 2d shows representative two-particle states that participate in either inter-or intra-band scattering. The two-particle states are constructed as convolutions of two single-particle states. For intra-band scattering (RR → RR), the initial and final states occupy the same region in harmonic space. In contrast, for inter-band scattering (RR → RL) the initial and final states are separated in harmonic space by a gap of order ∆/ω. (Note that this requires the energy spread of the single particle states, of order W , to be smaller than ∆.) Hence, at least within the Born approximation, inter-band scattering is strongly suppressed with respect to intra-band scattering.

B. Higher order contributions
Next, we consider higher order contributions to the inter-band scattering rate [67]. Importantly, the strong suppression of the Born-level interband scattering rate arises from the exponentially small overlap of harmonicspace wave functions of the initial and final states, along with the fact that the time-independent interaction conserves the harmonic index. As we now show, this exponential factor may be avoided at higher orders in perturbation theory, trading the small matrix element for higher powers of the interaction, U . Optimizing over the competition between small U and small overlaps, we find an optimal order N min that dominates the scattering amplitude. In this way we argue that the scattering rate should be a power law in U , with an exponential suppression in the inverse of the driving frequency, 1/ω.
To analyze the scattering amplitudes at higher orders in the interaction, we return to the expansion of the Tmatrix, Recall that the operators T (Ω), G 0 (Ω) and U are defined in the "extended space" of Fourier harmonics. As we will see, the additional state-space dimensions of the extended space play an important role in describing energy (photon) absorption from the drive.
To facilitate the T-matrix analysis, we define a basis of non-interacting eigenstates of the extended space Hamiltonian H 0 . For each non-interacting N particle Floquet state we associate a label ξ = {k n , α n }, n = 1 . . . N , denoting the quasimomenta and the band indices α n = {R, L} of all particles. As in the single particle case, cf. Eq. (4), the N particle non-interacting Floquet state |Ψ ξ (t) is represented by a vector of Fourier coefficients which we denote by ϕ ξ . We use the convention that the many-body quasi-energy is given by the sum of single particle quasi-energies, ε ξ = n ε 1P,knαn , with 0 ≤ ε 1P,knαn < 2π/T as above. Under this convention, the quasi-energies ε ξ generically fall outside of the fundamental Floquet zone. The detailed construction of the many-body Fourier vectors ϕ ξ is given in Appendix C.
In the extended space representation this spectrum is repeated over and over again, shifted by integer multiples of ω, with the corresponding eigenvectors likewise shifted in harmonic space. Therefore, with each Floquet state ϕ ξ defined above, we can associate a whole family of eigenstates {ϕ ξ,ν } with components ϕ states" we write the extended Green's function G 0 (Ω) as Crucially, the "copy states" defined above play an important role as off-shell virtual intermediate states in the perturbation theory. At each order, the scattering amplitude involves a product of many-body matrix elements of the form ϕ † ξ,ν Uϕ ξ ,ν . Such a matrix element can be expressed in the time-domain as is the periodic part of the (non-interacting) many-body Floquet wave function.
Importantly, H int is a two-body operator. Therefore, computation of the matrix element in Eq. (9) can be reduced to the problem of evaluating matrix elements between two-particle Floquet states, M . By further understanding the properties of these two-particle matrix elements, we may thus characterize the various terms in the high-order perturbation theory.
For simplicity we focus on the case W ω ∆, where the spread δm of the two particle Floquet wave-functions in Fourier space is of order δm = O(1), see Fig. 2c,d, and Fig. 3. The dependence of the two-particle matrix elements on ν − ν is distinguished by the Floquet band indices of the incoming and outgoing states. In particular, as seen for the special case of the on-shell process depicted in Fig. 2d, for an RR→ RR (intraband) transition the harmonic-space wave functions give an order-1 contribution (no suppression) to Eq. (9) for |ν − ν | δm. For an RR→ RL matrix element, order-1 overlap of the harmonic-space wave functions is achieved for ∆ ω − δm |ν − ν | ∆ ω + δm. Note that this nonsuppressed interband scattering matrix element necessarily describes an off-shell process. We neglect "double interband" processes RR → LL, the rates of which are heavily suppressed compared with those of the primary RR → RL decay process.
Using the rules above, we see that at high order in perturbation theory it is possible to construct an interband scattering amplitude that avoids the (super) exponential suppression arising from the tiny overlap of harmonic space wave functions that appears in the Born approximation. At each step in the perturbation theory, U may connect virtual states with |ν − ν | ∼ δm. Therefore, as indicated in Fig. 3b, there is an order N min ∼ ∆ δmω at which the initial (RR) and final (RL) two-particle states can be connected with no suppression due the harmonic space wave functions.
Viewing the construction of the N min -order scattering amplitude as a sequence of steps, each application of U reduces the harmonic-space separation between the virtual intermediate and the final (on-shell) state. These successive off-shell intermediates bring energy denominators of larger and larger magnitudes. After n steps, one particle is transferred from the R to the L band, giving a jump in the sequence of energy denominators. The terms in this sequence, in units of δmω, are of order −1, −2, . . . , −n, N min − n, N min − n − 1, . . . , 1, see Fig. 3c. We write the corresponding transition amplitude as A (n) inter ≈ (−1) n a n n!(N min − n)!(δmω) Nmin .
Here, a n ∝ U Nmin contains a sum of matrix elements of the interaction between the intermediate states. This gives the following rough estimate for the inter-band scattering rate (see Appendix D for details): Here, α is an O(1) numerical constant. Physically, in the regime W ω ∆ that we have considered, the dominant inter-band scattering processes are associated with N min scattering events, each one involving the absorption of δm energy quanta of ω from the driving field, such that the total absorbed energy is ∆. We expect such processes to dominate as long as ω is not too small compared to W . If W ω, the typical change in m in every scattering event scales with the width of the single-particle Floquet wavefunctions in harmonic space, δm ∼ W/ω. Hence we would have N min ∼ ∆ W in Eq. (11), and for low enough frequencies we expect Γ inter to saturate.

III. NUMERICAL SIMULATIONS
To further study the many-body dynamics we have performed exact numerical simulations of the dynamics of finite-size systems. To minimize the Hilbert space dimension we consider fermions in these simulations. We believe that the qualitative behavior does not depend on the quantum statistics, but leave the detailed study of the bosonic case to future investigations.
In each simulation, we initialize the system in a Slater determinant state where N momentum states in the right-moving (R) Floquet band are occupied, and the left-moving (L) band is empty. The results do not depend sensitively on which states in the R band are initially occupied, except at very short times. The particles move on a lattice of L unit cells (with two sites each), with periodic boundary conditions. The largest system we studied contained 8 particles with L = 16 unit cells. Figures 1c,d show the period-averaged current, as a function of n T , the number of periods elapsed (see figure caption for model parameters). At time t = 0, the system carries a current that depends on the initial state. Over a timescale τ intra ∼ 10 T the current relaxes to a value J ≈ ρ/T , see Fig. 1c (in this simulation, ρ = 0.5). We ascribe this time scale to relaxation within the R band, while the L band remains almost unpopulated. Examining the occupation numbers of momentum states confirms this interpretation (see Appendix E); for t τ intra , the occupation numbers in the R band are approximately uniform and all close to ρ. The average group velocity of states in the R band is v g = a/T . Therefore, the average current J = ρ a v g observed in Fig. 1c is as expected for a uniform particle distribution in the right moving band.
At longer times the current undergoes a much slower decay process toward zero, with a time scale τ inter of several hundred periods (Fig. 1d). During this process the population of the L band gradually increases. For t → ∞ the occupation numbers of all the momentum states in both bands tend toward ρ/2, corresponding to a maximally randomized infinite-temperature-like state.
For times up to several times τ intra , we find that the average (total) occupation number in the L band increases approximately linearly with time. The slope of the linear growth, which we define as Γ inter ≡ 1/τ inter , is found to be only weakly dependent on system size for L between 10 and 16 (see Appendix E).   Fig. 4 in a log-log plot, showing a power-law dependence on U . The power depends linearly on 1/ω (left panel, inset). The right panel shows the dependence on ω at fixed U . Clearly, Γ inter scales exponentially with 1/ω. For the lowest frequencies we studied, corresponding to ∆/ω ≈ 16 (where ∆ is the minimum instantaneous gap), the rate reaches ∼ 10 −7 in units of J 0 , indicating a very long-lived quasi-steady state where only the R band is populated. We also find that the behaviour of Γ inter shown in Fig 4 is not sensitive to the details of the band structure, and persists throughout a wide range of values of the band structure parameter λ.
Interestingly, the numerically obtained Γ inter (U, ω) is consistent with the form of Eq. (11). This suggests that, within our model, the inter-band relaxation process is dominated by "multi-photon assisted" scattering events, where many energy quanta are absorbed from the driving field.

IV. DISCUSSION
In this work we studied the dynamics of a periodicallydriven many-body system where the driving frequency is much smaller than the instantaneous inter-band gap throughout the driving period. When the system is prepared such that one of the bands is initially empty, while the other is partially occupied, a very long intermediate time window emerges, in which a universal quasi-steady state is realized. The quasi-steady state carries a robust current whose value depends solely on the density of par-ticles and the topological winding number of the underlying single-particle Floquet spectrum. In this sense, the combination of periodic driving and interactions leads to a non-equilibrium "prethermalized" state with topological properties.
Beyond the demonstration that such non trivial longlived quasi-steady states are possible in slowly driven many-body systems, our results may be directly applicable to recent experiments on cold atomic fermions [58] and bosons [57], where quantized pumping has recently been demonstrated. For the parameter regime stated above, initializing the system with a fractional filling of one of the bands should result in a current carrying quasisteady state.
The mechanism presented here is expected to apply to any many-body periodically driven system where the driving frequency is much smaller than the instantaneous gap to some subset of excitations. Under such conditions, the high-energy excitations can remain "frozen" for a long intermediate time window, leading to a quasisteady state with non-trivial properties. This work opens many intriguing theoretical and experimental questions about many-body dynamics and topology in driven systems, which will be interesting to explore in future work. e −imωt into the time-dependent Schrödinger equation, we get that the quasi-energy ε k and the Fourier components |ϕ (m) k satisfy: The operator on the left hand side of this equation is the Floquet "extended zone" operator, or "extended Hamiltonian," H 0,k . Note that H 0.k has a block-tridiagonal form, analogous to that of a nearest-neighbor tight binding Hamiltonian. In this picture, the terms proportional to ω on the diagonal correspond to a linear potential on the m-lattice.
In the adiabatic limit, ω ∆, we first solve the problem with ω = 0. The diagonal "linear potential" is later introduced as a perturbation. For ω = 0, the problem is translationally invariant in harmonic (m) space; we can solve it by Fourier transforming along the m-direction. This amounts to transforming from frequency back to time domain (where t plays the role of "momentum" for the m-space tight-binding problem). The secular equation then becomes which is nothing but the Schrödinger equation for the instantaneous eigenstates and eigenenergies. Here, α = R, L is the band index. The two bands are separated by a gap, |E L,k (t) − E R,k (t)| ≥ ∆. Next, we consider the effect of the linear potential term. The situation is illustrated in Fig. 5a. The problem is equivalent to the well-known Zener tunneling problem of a two-band semiconductor in an electric field. The states in both bands become localized in the "classically allowed region" in harmonic space, whose characteristic size is of order Within this picture, states in the R and L bands that are close in quasienergy (within ∼ ω of each other) are separated in harmonic space by a spacing of the order of ∆/ω. The tails of the wavefunctions decay rapidly into the classically forbidden region; in the limit ω W α , there is a broad region where the wavefunctions decay as e −A|m−m0| 3/2 (see Fig. 5b), where A is a band structure dependent dimensionless constant and m 0 is the boundary of the classically forbidden region. (This form is expected from a Wentzel-Kramers-Brillouin treatment of the problem, where m is treated as a continuous variable). At asymptotically long distances, larger than several times W α /ω, the form of the wavefunction crosses over to e −B|m−m0| log |m−m0| , where B is another dimensionless constant.
The tunneling matrix element between the two bands is hence strongly suppressed in the adiabatic limit. The hybridization between the two bands is expected to be very small, unless their energies are very close to each other. As k varies, the bands nearly cross at a set of k points (see Fig. 1b in the main text); at these crossing points the two bands hybridize, and there is an avoided crossing gap whose magnitude is exponentially small in the adiabatic limit. For a generic band structure, the hybridization between the bands is significant only within exponentially small regions in k space around the nearcrossing points. These hybridizations are responsible for the bright lines appearing in Fig. 2b in the main text. The extended zone Hamiltonian H in harmonic space (shown in Eq. (A1) for the single-particle case) is designed such that its spectrum and eigenstates correspond to the quasi-energies and Floquet states, respectively. It can also be used to generate the stroboscopic dynamics at times t = nT , where n is an integer, starting from an arbitrary initial state.
In (We use Dirac bra-ket notation for states in the physical Hilbert space, whereas vectors in the extended space are written without Dirac notation.) We relate the "evolved" state in the extended space to a state in the physical Hilbert space at stroboscopic times t = τ = N T by summing over the harmonic components of ϕ ext (τ ): |ϕ ext (t = N T ) = e −iεN T m |ϕ (m) . At these stroboscopic times, we find that |ψ(t = nT ) and |ϕ ext (t = nT ) coincide.
The reasoning above can be extended to the time evolution of an arbitrary initial state, |ψ 0 . This is done by expanding |ψ 0 in terms of Floquet eigenstates, considering the time evolution with respect to either H(t) or H, and comparing the time-evolved wavefunctions at times t = nT .
Appendix C: Construction of non-interacting ("free") many-body Floquet states We build up the non-interacting N -particle state ϕ ξ one harmonic at a time. First, we define a single particle Floquet state ϕ 1P,ξi , as in Eq. (4), for each ξ i = {k i , α i }, with i = 1 . . . N . As used throughout the main text, the single particle quasienergies ε 1P,ξi are taken in the interval 0 ≤ ε 1P,ξi < 2π/T . Within this convention we define a set of creation operatorsφ The m-th Fourier component of the many-body state ϕ ξ is determined by a convolution over the single particle harmonics: In the above, the sum extends over all sets of N integers m i , for i = 1, .., N . Note that the quasi-energy of the many-body state ϕ ξ is ε ξ = N i=1 ε 1P,ξi , with the convention 0 ≤ ε 1P,ξi < 2π/T specified above. Thus ε ξ is uniquely specified, and generically falls outside the fundamental Floquet-Brillouin zone.
Appendix D: Estimate of the interband scattering rate Here we describe the steps leading to the estimate of the interband scattering rate [Eq. (6) of the main text]. We begin with the term in the expression for the amplitude of order N min , in which a particle is scattered from the L to the R band after n steps, with 0 ≤ n ≤ N min −1. This amplitude (in the limit W ω) is given by Using Stirling's formula for n, N min 1, we may replace . The function f (x) is bounded and f (x) = O(1) for 0 ≤ x ≤ 1. Therefore we approximate the factor e −Nminf (n/Nmin) by β −Nmin 1 to leading exponential accuracy in the limit of large N min , where β 1 = O(1) is some constant.
In order to get a rough estimate of the interband scattering rate Γ inter , we must examine the factor a n in Eq. (6), which contains a sequence of matrix elements of the interaction term U between the initial, intermediate, and final states. The dominant dependence of a n on N min is expected to be of the form a n ∼ (β 2 U ) Nmin , where β 2 = O(1). Hence, summing the amplitudes A (n) inter over n, taking the modulus of the square, and using N min ∼ ∆ δmω , we get to the estimate for Γ inter quoted in Eq. (6) of the main text. While this estimate is rather rough, our numerical results for Γ inter show excellent agreement with the form predicted in Eq. (6).

Appendix E: Further details of the numerical simulations
In our numerical simulations, the many-body wavefunction |Ψ(t) is time-evolved numerically, using the Hamiltonian H(t) [Eqs. (1,2)], for up to 1000 driving periods. The simulations are performed using a finite time step, ∆t, and using a Trotter-Suzuki decomposition of the evolution operator within each time step. Most of the simulations were done with ∆t = T /500. For ω < 0.2J 0 , we used ∆t = T /850. We have verified the results do not change upon decreasing ∆t further, even for the longest times simulated.
To extract the inter-band rates plotted in Fig. 3 of the main text, we computed the distribution of particles in the different single-particle Floquet states at times which correspond to integer multiples of the driving period. We thus calculate N (α) k (t = mT ) = Ψ(mT )|ψ † k,αψ k,α |Ψ(mT ) whereψ † k,α is the creation operator for the Floquet state |ψ(0) k,α with momentum k, and the index α =R, L indicates the right or left moving Floquet band.
Typical particle distributions in the Floquet bands are plotted in Fig. 6. States within the right moving band quickly become nearly equally populated, with probabilities close to 0.5 (corresponding to the density ρ = N/L = 0.5, taken in the simulation). This indicates the establishment of a quasi-infinite-temperature state, restricted to the right moving (R) band. After an initial transient of a few driving periods, the population in the left moving (L) band increases linearly with time, with a small rate, while the population in the right moving band decreases with the same rate. The rates shown in Fig. 3 of the main text were obtained by considering the rate of increase of the average population in the left moving band (slope of the yellow line in the right panel of Fig. 6), Γ inter = The largest system that we could reach with our numerical simulations included 8 particles on 32 sites (i.e., L = 16 unit cells and density ρ = 0.5). To verify that the rates reported in Fig. 3 of the main text do not suffer from substantial finite size effects, we studied the size dependence of the interband scattering rate, normalized to the length of the system, Γ inter /(J 0 L), while keeping the density fixed. For the parameter regimes plotted in Fig. 3 of the main text, we found that the size-normalized rate is only weakly size dependent between L = 10 and L = 16, and appears to be saturating by L = 16. This indicates that at L = 16 finite size effects are small and do not change the results qualitatively. A representative plot of the finite size dependence of the size-normalized rate can be found in Fig. 7.