Universal transport in periodically driven systems without long-lived quasiparticles

An intriguing regime of universal charge transport at high entropy density has been proposed for periodically driven interacting one-dimensional systems with Bloch bands separated by a large single-particle band gap. For weak interactions, a simple picture based on well-defined Floquet quasiparticles suggests that the system should host a quasisteady state current that depends only on the populations of the system's Floquet-Bloch bands and their associated quasienergy winding numbers. Here we show that such topological transport persists into the strongly interacting regime where the single-particle lifetime becomes shorter than the drive period. Analytically, we show that the value of the current is insensitive to interaction-induced band renormalizations and lifetime broadening when certain conditions are met by the system's non-equilibrium distribution function. We show that these conditions correspond to a quasisteady state. We support these predictions through numerical simulation of a system of strongly interacting fermions in a periodically-modulated chain of Sachdev-Ye-Kitaev dots. Our work establishes universal transport at high entropy density as a robust far from equilibrium topological phenomenon, which can be readily realized with cold atoms in optical lattices.

The realization of topological pumps in metallic interacting systems is challenging due to an interplay of inter-particle interactions and the non-adiabatic evolu-tion stimulated by the periodic drive.Such an interplay often results in an extensive generation of entropy and incessant heating up of the system to a featureless, highentropy state [96][97][98].The heating can be significantly slowed down in the high driving frequency regime or under special conditions, giving rise to a long-lived prethermal state [99][100][101][102][103][104][105][106][107][108][109][110][111].Recently it was shown that a slowly driven topological pump in the weakly interacting limit can form a quasi steady state [112][113][114].In this limit, the quasi-steady state can be understood heuristically on the level of free dynamics and weak scattering of particles in well-defined Floquet-Bloch bands.Notably, the quasi steady state hosts a universal current that depends only on the populations of the system's Floquet-Bloch bands and their associated quasienergy winding numbers.
Here we show that the quasi-steady state persists into the strongly interacting regime where the single-particle scattering lifetime is shorter than the drive period.Furthermore, in this regime, the current exhibits a similar universal value as in the weakly interacting case, despite the absence of long-lived single-particle Floquet states.We support our analytical predictions through numerical simulations of a system of strongly interacting fermions in a periodically-modulated chain of Sachdev-Ye-Kitaev (SYK) [115][116][117] dots.Our work demonstrates a new approach for numerically exact simulations of driven, strongly interacting chains of many sites, by specializing to SYK-type interactions.This method outperforms the conventional exact diagonalization methods that can be applied to significantly smaller systems.In turn, approximate methods such as Hilbert space decimation (i.e., the time-dependent density matrix renormalization group [118,119]), can not be applied here, because thermalization dynamics generates long-ranged correlations.

II. DEFINITION OF THE PROBLEM
In this work, we consider a one-dimensional bipartite chain of L unit cells with periodic boundary conditions, hosting N flavors of otherwise spinless fermions (see Fig. 1a).We label the two sublattices by A and B, and denote the lattice constant by a.For simplicity of notation, throughout we set = k B = e = 1.
The chemical potential µ 0 sets the average density of the fermions in the chain (see Fig. 1b).For t < 0, the system is assumed to be in an equilibrium state with respect to the Hamiltonian Ĥ(t = 0), at inverse temperature β 0 .The interparticle interactions are described by the Hamiltonian In our analytical study, we assume a single flavor, α = 1, and generic short-ranged interactions of characteristic strength U αβγδ (j, j ) = U χ(|j − j |), where χ(|j − j |) is a rapidly decaying dimensionless function of its argument, with χ(1) = 1.(Note that for a single species of fermions, the on-site interaction terms, j = j , do not contribute.) In the numerical study, we consider the limit of a large number of flavors, N , with N → ∞, see Fig. 1a for an illustration.Particles of different flavors can locally interact through an SYK-type on-site interaction term, where we consider random and constant in space interactions with U αβγδ (j, j ) = 0 and U 2 αβγδ (j, j ) = δ jj U 2 /N 3 , such that the system preserves invariance to translations for every realization of disordered couplings [121,122].

III. NON-INTERACTING DYNAMICS
Before studying the interacting model, we briefly summarize the dynamics of the non-interacting topological pump [2,68,75,112,123,124].We initialize the pump in an equilibrium state of H(0), at inverse temperature β 0 and a chemical potential µ 0 that fixes the average density of particles at n 0 , see Fig. 1b.The spectral function is initially periodic in k, following the spectrum of the static Hamiltonian, H(0), as is demonstrated in Fig. 1b.
After switching on the drive (i.e., for t > 0), the dynamics of the system follows the time-dependent Hamiltonian H(t).Shortly after the quench, the bands of high intensity in the spectral function develop a pronounced structure of sidebands, spaced by the drive frequency Ω, and furthermore obtain nonvanishing net slopes [2] as a function of k (see Fig. 1c and attached video).The peaks of the spectral function at each k correspond to quasienergies associated with the single-particle Floquet state solutions [31,125] of the time-periodic Rice-Mele problem, [i∂ t − H 0 (k, t)]|Ψ ν (k, t) = 0, with the multiple values across the different sideband peaks capturing the indeterminacy of quasienergy up to integer multiples of the drive frequency, Ω.The single particle Floquet states [126,127] are given by |Ψ ν (k, t) = e −iεν (k)t m e −iΩmt |φ m ν (k) , where {|φ m ν (k) } are timeindependent states.Here, ε ν (k) is the corresponding quasienergy of the single particle state with crystal momentum k, and chirality (or Floquet band) index ν = {L, R} for the net left-and right-moving bands, respectively.We denote the bandwidth of the Floquet bands by W F and the gap between them by ∆.The net chiralities of the bands are determined by the topological index [2,31] , where the + (−) sign corresponds to the R (L) band.The chiralities of the bands are exhibited by the spectral function shown in Fig. 1c.As the momentum changes from k = −π/a to k = π/a, the peaks of the spectral function of the R (L) band shift in frequency from ω to ω ± Ω.

IV. TIME-EVOLUTION TOWARDS THE QUASISTEADY STATE
We now study the dynamics of the system in the interacting case.In particular, we focus on the evolution of the two-point Green's functions, providing information about expectation values of one-body operators, such as particle densities and the current (see below).The time-evolution of the interacting system's two-point Green's functions is described by the Kadanoff-Baym equations [128] where for brevity we have suppressed the crystal momentum and time indices on the right hand sides of these equations, and the • symbol indicates a convolution over time and matrix product in the sublattice indices.In Eq. ( 3), the (flavor-averaged) retarded and lesser Green's functions are defined as In these expressions, the expectation values are calculated with respect to the initial state described by the density matrix ρ 0 , describing an equilibrium state with respect to H(0) with temperature β 0 and average density of particles n 0 .The bar denotes averaging over the random interaction strength (in the case of the SYK interactions), "{, }" denotes an anticommutator, and θ(t) is the Heaviside step function.Throughout, we omit the sublattice indices s, s , leaving the 2 × 2 matrix structure of the Green's functions implicit.The retarded and lesser components of self-energy are denoted by Σ R (k; t, t ) and Σ < (k; t, t ), respectively (see Fig. 2a and Appendix C for technical details).

A. Single-particle spectral function
The renormalized single-particle spectrum of the nonequilibrium system is encoded in the retarded Green's function, whose time-evolution is given by Eq. (3a).In order to facilitate the separation of intraband and interband scattering processes below, we write the full retarded Green's function as a sum of R-and L-band projected Green's functions: The bandresolved Green's functions are defined through the Dyson series shown in Fig. 2a, corresponding to Dyson's equation, where, as in Eqs.(3a) and (3b), we suppress the crystal momentum and time indices on the last term for brevity.
Here, g R ν (k; t, t ) is the non-interacting retarded Green's function projected to band ν, see Appendix A for more details.
We define the renormalized band-resolved spectral functions as The renormalized bandresolved spectral function broadens due to interactions, with tails extending into the gap that decay approximately as e −ω/ξ .In what follows, we focus on the limit ξ ∆, where the band-resolved spectral functions are well separated in frequency, see Figs. 2b,d.This separation of the renormalized Floquet bands in the frequency domain is crucial for obtaining a long-lived quasi steady state in the system (see below).
The renormalization of the single-particle spectral function is caused by the dressing of the non-interacting Floquet bands by virtual electron-hole pair creation and annihilation processes.Our analytical estimate near the quasi steady state (see Appendix C) suggests that the broadening of G where f 0 L and f 0 R are the occupation probabilities of the Floquet bands (see below) and f 0 ν = 1 − f 0 ν .

B. The kinetic equation and population dynamics
To study the formation and properties of the quasisteady state, we define occupation probabilities by parametrizing the lesser Green's function as where crystal momentum and time indices are suppressed on the right hand side.By analogy to the case of thermodynamic equilibrium, where the (Fourier-transformed) lesser Green's function and the spectral function are related via the Fermi occupation function, in Eq. ( 5) the Hermitian matrices f R (k; t, t ) and f L (k; t, t ) play the roles of distribution functions for the two bands.As we will discuss further below, this interpretation is most meaningful when f R (k; t, t ) and f L (k; t, t ) take simple forms in terms of their matrix and time (or frequency) structures.We will see that such a simple form naturally emerges in the quasisteady state of the system.
To assess the dynamics of the distribution functions we derive a kinetic equation for ∂tf L (k; ω, t), where f L (k; ω, t) is obtained via the Wigner transform of the two-time function, f L (k; t, t ).A similar approach is often employed in studies of weakly interacting dynamical systems [129][130][131].Here, we generalize this approach to the strongly interacting case, without imposing the "onshell" approximation [132].Crucially, such an approximation would not correctly capture interband scattering, which occurs through high-order processes in U .
The kinetic equation is obtained by combining the Wigner transforms of Eqs. ( 5) and (3b) and subtracting terms proportional to ∂tG R (k; ω, t), which itself is described by Eq. (3a) (see Appendix D for the full derivation).Generically, the kinetic equation for f L can be written as [133] where the collision integral I L is a functional of f R and f L with a matrix structure in the sublattice indices.Notably, for values of ω where G ∆ L (k; ω, t) has significant weight, i.e., near the lower band (see Fig. 2d), the net collision integral is exponentially small if its arguments are independent of k and ω [see appendix D]: , where f 0 L and f 0 R are constants and 1 denotes the identity matrix in sublattice space.In particular, we estimate , characterized by uniform occupation within each band, exhibits exponentially slow population dynamics and thus describes a long-lived quasisteady state of the system.

C. Universal value of the current
Using the form of the quasisteady state found above in terms of two-point Green's functions, we can now characterize observables in the quasisteady state.In particular we focus on the value of the time-averaged current, which was previously conjectured to take a universal value based on a weak-coupling picture and evidence from numerical simulations on modestly-sized systems.
The instantaneous current averaged along the chain in a generic translation-invariant state described by G < reads where the momentum integral is performed over the first Brillouin zone.Next, we evaluate Eq. ( 6) in the quasisteady state given by G 5) and the following discussion].At equal times, G < (qs) can be further simplified using G ∆ ν (k; t, t + 0 + ) = g ∆ ν (k; t, t + 0 + ), due to the fact that the time-integrals in Eq. ( 4) vanish for t = t .Substituting the resulting quasisteady state form of G < (qs) into Eq.( 6), we obtain is the current carried by Floquet band ν of the system in the absence of interactions, when fully filled (see Appendix A).As defined in Sec.III above, |Ψ ν (k, t) denotes the single-particle Floquet state with crystal momentum k in band ν of system in the absence of interactions.In the adiabatic limit, the period averaged current J 0 T .In a system where the upper (R) band is initially empty and the lower (L) band has fractional filling ν 0 = an 0 , we therefore expect f 0 L = ν 0 and f 0 R = 0, such that the current in the quasisteady state is equal to J (qs) = ν0 T + O(e −∆/ξ ), where the correction captures the deviation of the quasi steady state from the maximal entropy state in which the upper band is empty.FIG. 3. Relaxation to the quasi steady state.a "Distance" to equilibrium, Deq(t), defined in Eq. ( 7), as a function of time.Lines of different colors in this and following panels correspond to the different interaction strengths indicated on the top of the figure, in units of J0 = 1.Shaded areas indicate error bars due to the extrapolation procedure, see text.At t = 0, the extrapolation yields non-physical, negative values of Deq(t) due to the quench, we thus cut-off these values in the plot.b Effective temperature of the lower band minimizing Deq(t) in units of J0 = 1, as a function of time.Blue dashed line indicates the point where T eff = WF.c Entropy density of the system's one-body reduced density matrix as a function of time, defined in Eq. (8).A dashed line indicates maximal entropy for a quarter-filled system, when all the particles occupy the lower band.d Period averaged current normalized by the total density of particles and driving frequency, J (t)(T /ν0) as a function of time.
This is a remarkable result: even in a limit where the single particle scattering lifetime may be short compared with the driving period, where the single particle Floquet states and associated spectrum are not well-resolved or defined, the current still attains a universal value associated with the nontrivial topology of the system's singleparticle Floquet spectrum in the absence of interactions.The universal value of current holds, up to an exponentially small correction, provided that the scattering rate (and associated level broadening, captured by ξ) remains small compared with the single particle band gap, ∆.

V. NUMERICAL ANALYSIS
To benchmark our analytical results, we numerically simulated the model given in Eqs. ( 1) and ( 2), with SYK interactions.In addition to the terms described above, we also included a weak random quadratic term ĤSYK−2 = j,αβ K αβ ĉ † jα ĉjβ + h.c., where K αβ = 0 and K 2 αβ = K 2 /N , with K = 0.05.This additional term is essential for stabilizing the numerics in the weakly interacting regime.
The unique structure of the SYK interactions allows us to simulate considerably larger systems compared to exact diagonalization methods applied to systems with conventional interactions.Here, we simulated the time evolution of a chain of 100 SYK dots arranged into L = 50 unit cells.We used the Kadanoff-Baym equations [given in Eq. ( 3)] to evolve the Keldysh-ordered Green's functions in time [110,[134][135][136][137]; for further details see Appendix E. The system is initialized in an equilibrium state of Ĥ(0) with temperature β −1 0 = 0.1 and µ 0 = −2.93,which is set to fix the density of electrons approximately at quarter filling, see Fig. 1b.For the model itself, we select the parameter values: J 0 = 1, J 1 = 0.85, v 0 = 2.55, Ω = 0.5, see Eq. ( 1) and surrounding text for the definitions of the Hamiltonian and its parameters.We note that all the energies and frequencies are given in units in which J 0 = 1.
The time-evolution algorithm is based on the discretization of the time and frequency domains, with small steps δt and δω, respectively.We performed the evolution for several values of steps in the range 0.08 ≤ δt ≤ 0.16 and 0.015 ≤ δω ≤ 0.04, and performed a two-dimensional linear extrapolation to δt = δω = 0.In the numerical results we present the extrapolated values with error bars indicating the uncertainty of the extrapolation procedure (defined as the difference between the extrapolated value and the closest numerically-determined point).

A. Formation of the quasisteady state
We first analyze the formation of the quasisteady state, wherein the distribution functions f L and f R for the two bands become independent of crystal momentum and frequency, while the total populations of the two bands remain approximately constant.As a means of characterizing the nonequilibrium state of the system, and to enable the extraction of an effective temperature (when it is appropriate to do so), we define the "distance" of the distribution from a thermal equilibrium state as where ω c is the center of the spectrum and the Green's functions are evaluated in the Wigner transformed representation and averaged over one period: −T /2 ds G < (k; ω, t + s) and similarly for G ∆ (k; ω, t).
In Figs.3a,b we plot the evolution of the distance to equilibrium D eq ( t) and the extracted effective temperature T eff = β −1 min , corresponding to the minimization in Eq. (7).When the drive is switched on at t = 0, D eq rapidly grows, indicating evolution into a far from equilibrium state.Following the rapid rise, the system relaxes to the quasisteady state, which is manifested by the decay of D eq .In parallel, the effective temperature grows approximately exponentially with a rate Γ intra : T eff ∼ β −1 0 e Γintrat (Fig. 3b).The quasisteady state observed in the simulation approximately realizes the conditions discussed in Sec.IV B, once the effective temperature exceeds the width of the single-particle Floquet bands, T eff W F , indicated by blue dashed line (and for D eq 1).The curves of different colors in Figs.3a,b correspond to different interaction strengths, U ; the time to reach the quasisteady state rapidly decreases with interaction strength, U .
To track the system's evolution towards a high entropy density state, we calculated the average von-Neumann entropy density of the system's one-body reduced density matrix where Ḡ> = −i Ḡ∆ + Ḡ< and the Green's functions are evaluated at equal times t = t = t, see Fig. 3c.The value of S( t) for a maximal entropy density state in a quarterfilled system, subject to the constraint that all the particles occupy the lower band, is given by S max L = log(2)/a.As can be seen in Fig. 3c, the entropy density stabilizes slightly above this value due to a small population excited to the upper band at t = 0.After stabilizing near S(t) ≈ S max L , the entropy slowly grows further due to interband transitions.In the infinite time limit, we expect S → S max ≈ 1.12/a corresponding to one quarter filling of the entire system.
In Fig. 3d, we extracted the period-averaged current normalized by the filling, J /ν 0 [see Eq. ( 6)].As follows from the discussion below Eq. ( 6), we expect an approximately quantized value in the units of T −1 for the normalized current in the quasisteady state.Fig. 4b shows the period-averaged current normalized by ν 0 as a function of the stroboscopic time.The gray strips indicate the uncertainty intervals associated with the extrapolation to infinitesimal grid spacing in the simulations, as described above.In the regime of strong interactions, the average current rapidly increases on a timescale set by ∼ Γ intra .When the quasisteady state is reached, the current obtains the expected universal value to within the uncertainties of the numerical simulation, as expected.For later times, the current slowly decays with the rate Γ inter due to interband heating.In the weakly interacting case, the normalized current remains non-universal for much longer times; for these cases, the quasisteady state was not reached within the time window that we were able to simulate.Interestingly, slowly-driven Fermi liquids have been shown to persist in non-thermal states for parametrically long timescales [110].The connection between the slow intraband heating observed here and the mechanism in Ref. a Density of electrons excited to the upper band as a function of time, in the interacting and non-interacting cases.Circles represent the period averaged density.In the interacting case, the density of excitations increases with an approximately constant rate Γinter, while in the non-interacting case the average charge is constant, following an initial jump at t = 0. b Intraband (Γintra) and interband (Γinter) equilibration rates as a function of the interaction strength.The intraband equilibration rate is extracted from the temperature growth [see Fig. 3b].The interband equilibration rate is extracted from the slope of the period-averaged excitation density, in a.
gate in future work.

B. Interband heating and universal current
To investigate the interband scattering processes and measure their rates, in Fig. 4a we extracted the density of excitations in the (renormalized) R band from our simulations.We define the excitation density as n ex ( t) = −i ∞ ωc dω 2π dk 2π Tr {G < (k; ω, t)}.The periodaveraged density of excitations jumps at t = 0, when the drive is switched on [124], and then gradually increases with an approximately constant rate Γ inter .The interband equilibration rate, Γ inter , compared to the intraband equilibration rate, Γ intra , is shown as a function of the interaction strength in Fig. 4b.Fig. 2c shows the period averaged lesser Green's function as a function of momentum and frequency, after 30 periods of the drive.The excited population can be seen as a pale strip at the location of the upper band.Note that the features are heavily broadened due to fast intraband scattering, such that for the selected parameters the Floquet harmonic side bands are nearly completely washed out.

VI. DISCUSSION AND OUTLOOK
Periodically driven systems can host Floquet-Bloch bands with unique topological properties that cannot be obtained in equilibrium systems.In the presence of interactions, it is natural to wonder if rapid scattering on timescales comparable to the driving period might mask any dynamical features expected to arise from the singleparticle Floquet states.In this work we showed that this need not be the case: even in a state with high entropy density and rapid scattering, universal transport associated with the topological properties of the system's Floquet Bloch bands persists.We demonstrated this phenomenon in the context of a topological pump with non-integer filling, which exhibits a long-lived quasisteady state with maximal entropy density (subject to the constraint of fixed particle number in each band).We derived conditions under which the quasisteady state hosts quantized transport (in units of the particle density), up to an exponentially small correction in the ratio of the system's band gap to its renormalized band width.
To support these arguments, we studied this phenomenon numerically in an SYK-type chain.This setup enabled us to examine the dynamics in a regime of strong scattering, in system sizes much larger than could be accessed by exact evolution.This advantage is gained from the fact that the SYK system can be solved by time evolution of the Kadanoff-Baym equations.Our numerical results allowed us to study the dynamics leading to the formation of the quasisteady state.Importantly, we showed that quantized transport persists even when quasiparticles are short-lived due to fast intraband scattering and the Floquet sidebands are hence not well resolved.
Appendix B: Definition of the Keldysh-Floquet Green's functions in the gauge-invariant form Due to the non-equilibrium nature of the Floquet-Keldysh Green's functions, the energy in the collision processes is only conserved modulo Ω.This property complicates the calculations using the Keldysh formalism, as multiple photon absorption/emission processes have to be take in account in each collision.Here, we present a gauge invariant definition of the Green's functions in which the ω index is conserved in the collision processes, and show how convolution and product of the the two-time Green's functions are defined with this gauge choice.Such a definition, allows us to operate the Keldysh-Floquet Green's functions as equilibrium Keldysh Green's functions with additional matrix structure in the Floquet harmonics.
Given the Wigner-transformed Green's function G(ω, t) [see Eq. (A5) for definition], we define the Green's function in the harmonic basis (with indices m, n ∈ Z) as This definition is invariant under the transformation G m,n (ω + lΩ) = G m+l,n+l (ω) for any integer l.A convolution C(t, t ) = dsA(t, s)B(s, t), reads Similarly, a product of two same-time functions, P (t, t ) = A(t, t )B(t, t ) is given by Appendix C: Evaluation of the self-energy and renormalization of the spectral function Here, we estimate the broadening ξ of the renormalized bandwidth of the Floquet bands, discussed in Sec.IV A. In all the expressions in this section, we assume the Green's functions in the frequency domain are defined as matrices in the harmonic basis [as defined in Eq. (B1)], and implicitly contract the harmonic indices following the rules given in Eqs.(B2) and (B3).The renormalization of the spectral width can be understood from the definition of the renormalized Green's function in terms of the bare one, following from Eq. (3a).Focusing on the values of ω at the gap of the bare function, i.e., where g ∆ (ω) = 0, the renormalized spectral function reads Therefore, to estimate the broadening of the spectral function, we need to estimate Σ ∆ (ω).In what follows, we estimate the lesser and greater components of the self-energy Σ ≶ (ω), constituting the spectral component, Σ To estimate the self-energy, we need to sum over all irreducible diagrams allowed by the interaction term (first four terms in the expansion are demonstrated Fig. 2a).We begin by estimating the greater component of the selfenergy, Σ > (ω).Consider a generic irreducible diagram in this sum, corresponding to the order U p in the interaction strength, see Fig. 5.Such a diagram contains p interaction vertices evaluated at the times t 1 , t 2 ,.., t p and positions on the Keldysh contour i 1 , .., i p , where i p = ± corresponding to positive (+) and negative (−) branches on the Keldysh contour.The vertices are connected by the non-interacting propagators g i,i (t, t ), with the convention g +− = g < , g −+ = g > , g ++ = g T and g −− = g T .For a specific combination of i p , it is useful to arrange the diagram such that all the vertices on the positive branch are on the left side and all the vertices on the negative branch are on the right side, see Fig. 5. Notice that by the definition of Σ > the incoming vertex belongs to the negative Keldysh branch and the outgoing vertex belongs to the positive branch.
Next, we transform the expression for the self energy following the transformation given in Eq. (B1).For convenience, we enumerate the frequencies of the l propagators going from the left to the right side of the diagram by ω 1 , ω 2 ,..., ω l , and l − 1 propagators going from the right side to the left side by ω l+1 , ω l+2 ,..., ω 2l−1 .The maximal value of l is limited by p.We define the sum of the frequencies of the propagators crossing the center of the diagram with opposite signs by ∆ω = ω 1 + ... + ω l − ω l+1 − ... − ω 2l−1 .From the conservation of frequency, ∆ω = ω, where ω is the frequency associated with self-energy.
By construction, the propagators directed from the left to right correspond to the greater component of the Green's function, g > (ω 1 ),..., g > (ω l ) while the ones directed oppositely correspond to the lesser components of the Green's function, g < (ω l+1 ),..., g < (ω 2l−1 ).For simplicity, we suppressed the momentum dependence and the Floquet harmonic indices, as they are not important for the qualitative result.Therefore, the greater component of the self-energy of a single diagram to the order U p is given by The function F p includes the contribution from all the propagators that do not cross the center of the diagram, which include the time-ordered and anti time-ordered components.It is analytical and non vanishing, therefore can not change the support of the convolution in the ω domain, yet it may change the weight of the function.
Our goal is to estimate Σ > (ω) for ω away from the support of the non-interacting density of states.The dominant contribution to the self-energy would arise from ladder-shaped diagram, as is shown in Fig. 5c.Such a diagram has maximal number of lines crossing the center of the diagram, with minimal constraints on the intermediate frequencies.
The unique topology of this diagram allows to write it as a recursive relation with the initial condition For simplicity, we assume a nearly quasi steady state in which the lesser and greater functions can be approximated by g > (ω The dominant contribution arises from diagrams where each pair of propagators in Eqs.(C4) and (C5) corresponds to the same band, i.e., g > (ω 1 )g < (ω 2 ) → (f L fL + f R fR )g ∆ 0 (ω 1 )g ∆ 0 (ω 2 ), where g ∆ 0 (ω), is the density of states of one of the bands shifted to the center of the energy.We also approximate g T (ω)g T (ω) ≈ 1/ω 2 .
The solution to Eq. (C5), can be written as where g∆ ν (d, ω) is a function of width d centered around the ν-th band.Applying Eq. (C4), we obtain Similarly, after p iterations, we arrive at Σ > 2p (ω) We separate the self energy to Σ For a distance |ω| from the left band, the leading order in p is proportional to p ∼ |ω/W F |. Therefore, the self energy to this order reads Σ > L (ω) . For ω of the order of ∆, we arrive at Σ > L (ω) ∝ iU fL e −|ω|/ξ , where Similarly, the lesser component of the self energy reads Σ < L (ω) ∝ −iU f L e −|ω|/ξ .Therefore, Σ ∆ = iΣ > − iΣ < ∝ U e −|ω|/ξ .Using Eq. (C2), we obtain A similar calculation near the upper band leads to the same energy scale for the broadening.

Appendix D: Kinetic equation
In this section, we derive the kinetic equation for the occupation probabilities, defined in Eq. ( 5), and demonstrate that this equation is solved by constant occupations of the bands f L = f 0 L and f R = f 0 R , up to terms proportional to O(δf e −∆/ξ ), where δf = f 0 L − f 0 R .This means that the fixed point of the kinetic equation to the order e −∆/ξ , corresponds to an infinite temperature distribution in each of the bands.If these terms are included, the fixed point is a global infinite temperature state, in which δf = 0.
To derive the kinetic equation , we substitute Eq. ( 5), in Eq. (3b).Before, performing the substitution, we rewrite Eq. (3b) in the frequency-momentum domain by performing the Wigner transformation of the timefrequency domain, yielding Here, "•" denotes the Moyal and matrix product and To the first order in the derivatives and commutators the Moyal commutator reads, Our goal is to derive an equation for ḟL and ḟR without imposing the "onshell" approximation, which otherwise would not include the interband transitions occurring off-shell.Eq. (D1) describes the time-evolution of G < , which includes timeevolutions of both the spectral function and the occupations [see Eq. ( 5)].On the other hand, the evolution of G ∆ alone can be derived from Eq. (3a), and reads To separate the kinetic equation for f L from the kinetic equation for G ∆ , we define Evaluating the l.h.s. of Eqs.(D1) and (D2), we obtain To simplify, we rewrite Eq. ( 5) as To simplify even further, we assume a close to steady state, such that we can keep only the first order terms in the derivatives and commutators of f ν .With this assumption δG < is already given to the leading order and its Moyal commutator will be of higher order and thus can be neglected.In addition, focusing on ω ≈ −∆/2, the spectral function scales as G ∆ R ≈ e −∆/ξ [see Eq. (C7)].Therefore, this term can be neglected compared to f L • G ∆ L .With these approximations and explicitly computing the Moyal commutator of f L in Eq. (D3) to the leading order, we obtain, Eq. (D4) is essentially the l.h.s. of the Boltzmannlike equation (up to the band-renormalization terms discussed below).Similarly, we evaluate I L using the r.h.s. of Eqs.(D1) and (D2), yielding the collision term (and bandrenormalization terms).An explicit calculation yields, Evaluation of Eq. (D5) in a generic state is complex and is performed in the numerical part.As a first order check, we will verify that an infinite temperature state for each of the bands f ν = f 0 ν , with f 0 L = f 0 R almost nullifies the collision integral and estimate the timescale for the full thermalization of the system to an infinite temperature state (in which f 0 R = f 0 L ).Under the assumption of constant occupations, Eq. (D5) simplifies to Here, we denote by Σ L (ω) the self energy near ω = −∆/2, and similarly Σ R (ω), denotes the self-energy near ω = ∆/2.We also used Σ Importantly, all the terms in I L are proportional to ∆f and to either G ∆ R or Σ ∆ R , that are exponentially small near ω = −∆/2.This exponentially small value of I L is proportional to the rate of thermalization to the infinite temperature state, in which δf = 0.In this state Eq.(D6) becomes identically zero, I L = 0.

Appendix E: Details of the numerical simulation
Here, we present the details of the numerical simulation of the time-evolution of the Kadanoff-Baym equations [see Eq. ( 3)].The time-evolution is performed with respect to the SYK Hamiltonian, given in Eq. ( 2), with N → ∞.To stabilize the numerics in the weaklyinteracting limit, we added a weak random quadratic term where K αβ = 0 and K 2 αβ = K 2 /N .Instead of evolving in time the G R and G < functions, as appears in Eq. ( 3), we found it more convenient to evolve the retarded G R , and Keldysh components G K .The latter is defined as The Kadanoff-Baym equations for these two components read The disorder-averaged self-energy for the chain of SYK dots with SYK-4 [Eq.( 2)] and SYK-2 [Eq.(E1)] interactions can be written in the self-consistent form.The diagrammatic structure of the self-energy to the leading order in N is shown in Fig. 6

(E4)
Here, x is obtained from the Fourier transform of the momentum, k, index.The retarded and Keldysh Green's functions G R and G K [appearing in Eq. (E3)], are related to G ≷ via G R ss (k; t, t ) = θ(t − t )[G > ss (k; t, t ) − G < ss (k; t, t )] and Eq.(E2).In turn, the inverse relatons read G ≷ = 1 2 G K ± (G R − G A ) , where G A ss (k; t, t ) = G R s s (k; t , t) † .The relations for the self energy are similar, with G replaced by Σ.

Equilibrium solution
Equations (E3) and (E4) constitute a set of integrodifferential equations which determine the time evolution of the Green's functions.Initial conditions for this timeevolution are set by the state ρ 0 , corresponding to equi-librium with an inverse temperature β 0 and Hamiltonian Ĥ(0).Due to invariance to time-translations in equilibrium, the equilibrium Green's functions depend only on the time-difference ∆t = t − t .To find the equilibrium solution, we we evaluate Eq. (E3a) for H 0 (t) at t = 0 and transform ∆t to the frequency space ω, giving rise to [ω − H 0 (0)]G R (ω) = 1 + Σ R (ω)G R (ω).

Time evolution
Having found an equilibrium solution, G R eq (k; t − t ) and G K eq (k; t − t ), we rearrange the vectors into matrices [G R eq ] ss (k; t, t ) and [G K eq ] ss (k; t, t ) of size N t × N t in the time domain and 2×2 in the sublattice space, for a vector of crystal momenta of size L. In our simulations, we used N t = 3000 (smaller values of N t are used to vary δω and δt), and L = 50.We used the equilibrium solution as the starting point of the simulation to propagate the Green's functions by one time step δt in each iteration according to Eq. (E3) [134], see Fig. 7.In particular, given G R (t 0 , t ), we evolve G R according to: for t ≤ t 0 ; G R (t, t 0 + δt) = 0 for t ≤ t 0 ; and G R (t 0 + δt, t 0 + δt) = G R (t 0 , t 0 ).Here we defined U 0 (t) = e −iδtH0(t) and I R (t ) = Σ R (t 0 , s)G R (s, t )ds.
To optimize the efficiency of the simulation, we keep the overall size of the matrices constants.Therefore, for each new element of the Green's function calculated in the future, we erase one element in the past.Such a truncation of the Green's function fixes the required memory of the simulation and significantly reduces the computational resources used.

FIG. 1 .
FIG.1.The model and the single-particle spectrum.a Illustration of the system.A chain of atoms with timemodulated hopping and staggered potential, following the driving protocol realizing an adiabatic pump, as shown in the right bottom corner of panel a (we denote ∆d(t) = d+(0, t) − d−(0, t)).In the numerical simulation, each site is an SYK dot with N → ∞ orbitals that interact through random interactions [see Eq. (2)].b The single-particle spectrum of the system at t < 0, before the drive is switched on.The dashed line represents the Fermi level and the colored section represents initially occupied states, constituting a quarter of all the states in the system.The shaded gray areas show the instantaneous levels of the driven system in one period.The band structure is calculated for J0 = 1, J1 = 0.85, v0 = 2.55, and Ω = 0.5 [see Eq. (1) for the definition of the single-particle Hamiltonian].c The period averaged spectral function of the non-interacting system, ḡ∆ (k; ω).Energy and frequency in panels b and c are plotted in units of J0 = 1.d The period-averaged non-interacting density of states (DOS).

FIG. 2 .
FIG. 2. Renormalization of the spectrum by interactions.a First row: Dyson's expansion for the band resolved Green's function.The double and single blue lines indicate the band-resolved renormalized and bare functions respectively.Black lines denote the full bare Green's function.Second row: Diagrammatic expansion of the self-energy, used in the analytical model (shown up to order U 2 ).In the numerics, we considered SYK interactions in the infinite N limit, in which all the diagrams with odd number of interaction vertices or with crossed lines are averaged to zero (i.e., only the last diagram in a contributes).b Period-averaged spectral function, Ḡ∆ (k; ω), renormalized by the SYK interactions with U = 1.c Period-averaged lesser function, i Ḡ< (k; ω), indicating occupation after 30 periods of the drive, for the same parameters as in b.Note the weak, yet finite occupation of the upper band produced by the interband processes.d The renormalized period-averaged single-particle DOS.In the presence of interaction, the non-interacting DOS appearing in Fig. 1d, broadens obtaining exponential tails ξ.The broadening creates an overlap between the upper and lower bands, forming an interband heating channel (indicated by the wiggly arrow).Frequencies are plotted in units of J0 = 1.

5 FIG. 4 .
FIG.4.Heating and current in the quasi steady state.a Density of electrons excited to the upper band as a function of time, in the interacting and non-interacting cases.Circles represent the period averaged density.In the interacting case, the density of excitations increases with an approximately constant rate Γinter, while in the non-interacting case the average charge is constant, following an initial jump at t = 0. b Intraband (Γintra) and interband (Γinter) equilibration rates as a function of the interaction strength.The intraband equilibration rate is extracted from the temperature growth [see Fig.3b].The interband equilibration rate is extracted from the slope of the period-averaged excitation density, in a.

FIG. 5 .
FIG.5.Examples of diagrams contributing to the self energy, Σ > (ω), arranged such that the vertices on the positive Keldysh branch are at the left side and the vertices on the negative Kelysh branch are at the right side.Notice that the left diagram would not contribute in the case of the SYK interactions, since it is subleading in the number of the SYK flavours.
FIG.6.Self consistent relations for a Thouless pump with SYK interactions.Single lines represent the non-interacting Green's function g ss (k; t, t ).Double lines represent the renormalized Green's functions G ss (k; t, t ).
[117].As follows from ′ G R (t, t ′ ), G K (t, t ′ ) Time evolution of the retarded and Keldysh Green's functions according to the Kadanoff-Baym equations [Eqs.(E3)].The square in the bottom left corner represents the initial conditions for the Green's function, Geq(t, t ).At each step of the evolution a new row, column and a cell on the diagonal are added to the matrix, corresponding to G(t0 + δt, t), G(t, t0 + δt) and G(t0 + δt, t0 + δt).this diagram, the greater and lesser components of the self-energy are given by Σ ≷ ss (x; t, t ) =K 2 G ≶ s s (−x; t , t).