Floquet prethermalization and regimes of heating in a periodically driven, interacting quantum system

We study the regimes of heating in the periodically driven O(N)-model, which is a well established model for interacting quantum many-body systems. By computing the absorbed energy with a non-equilibrium Keldysh Green’s function approach, we establish three dynamical regimes: at short times a single-particle dominated regime, at intermediate times a stable Floquet prethermal regime in which the system ceases to absorb, and at parametrically late times a thermalizing regime. Our simulations suggest that in the thermalizing regime the absorbed energy grows algebraically in time with an exponent that approaches the universal value of 1/2, and is thus significantly slower than linear Joule heating. Our results demonstrate the parametric stability of prethermal states in a many-body system driven at frequencies that are comparable to its microscopic scales. This paves the way for realizing exotic quantum phases, such as time crystals or interacting topological phases, in the prethermal regime of interacting Floquet systems.

The eigenstate thermalization hypothesis (ETH) suggests that generic interacting many-body systems heat up to infinite temperature 26,27 , thus inhibiting the realization of such novel phases. A possible resolution is to stabilize the Floquet states by disorder such that the system becomes many-body localized and ETH does not apply [8][9][10]13 , as recently demonstrated experimentally 14 . However, this restricts the variety of accessible phases. Another route would be to resort to driving frequencies much higher than all other microscopic scales [28][29][30][31] . But in that case Ĥ F becomes quasi-local and cannot possess any exotic phases. A more general approach, is to resort to a transient prethermal regime [32][33][34][35] , which is characterized by the interaction time scale of the Floquet Hamiltonian required to realize exotic phenomena being much shorter than the heating timescale. It is therefore eminent to study the stability of such a Floquet prethermal regime in a general context.
In this work, we investigate the stability of the Floquet prethermal regime and the thermalization time scales in a generic interacting many-body system subject to a periodic drive. In particular, we focus on the quantum O(N)-model with modulated mass. To this end, we employ the 2-particle irreducible (2PI) effective action approach on the closed Keldysh contour including corrections up to next-to-leading order (NLO) in 1/N which allow the system to thermalize. The O(N)-model is a well established model for interacting many-body systems, both in condensed matter and cosmology 34,[36][37][38][39][40][41][42][43][44][45] . In particular, the presence of nontrivial interactions at NLO as well as the bosonic nature of excitations render the O(N)-model useful for studying heating of a driven many-body system to infinite temperature. Based on our numerical simulations, we find a parametrically large regime of Floquet prethermalization, even when the driving frequency is comparable to other microscopic scales of the undriven Hamiltonian so long as the interactions of the system are not too strong; Fig. 1

Model and nonequilibrium Keldysh formalism
We study the quantum O(N)-model of N real scalar fields Φ a , a = 1, … , N with the action 46 We use the abbreviation ∫ ∫ ∫ where the time integration runs over the closed-time Keldysh contour . Furthermore we assume that repeated indices are summed over. In momentum space a finite cut-off Λ is applied to regularize eventual UV divergencies. Consequently, we are effectively discussing a lattice system with a finite quasi-particle bandwidth. The bare mass is driven with amplitude A and frequency Ω, which, in a linear response regime ( Ω  A m , 0 ) creates pairs of excitations. It is convenient to rescale time t → 2t/Ω and the fields Φ a → (2/Ω) 1/2 Φ a , and to introduce the effective coupling constants in the presence of an external drive: The driving amplitude is rescaled by Ω 2 , which is a consequence of the relativistic form of the action. The model (1) displays an equilibrium phase transition to an ordered, symmetry broken phase for < m 0 0 2 and small λ < λ c at low temperatures. The drive destroys the ordered phase already at leading order in 1/N 34 . Hence, as we are interested in the long-time dynamics, we restrict ourselves to initial states in the symmetric phase. Furthermore, in the case of symmetric initial states, we find the same qualitative behavior in all spatial dimensions d = 1, 2, 3, and thus the presented results focus on d = 1. We emphasize that our results represent the thermodynamic limit, and thus should be contrasted to the exact diagonalization of small systems.
In order to simulate the dynamics of the driven system, we use the nonequilibrium Keldysh formalism 47 . The time evolution of the two-point contour ordered Green's function Ĝ is governed by the self-consistent Dyson equation where □ t,x is the d' Alembert operator and the self-energy Σ is given as the functional derivative of the 2PI effective action Γ 2

48
, see the methods section for details. The advantages of this approach are that it operates in the thermodynamic limit and respects the conservation laws associated with the global symmetries of the microscopic action, such as energy or momentum conservation. We decompose the contour ordered Green's function as , with the Keldysh or statistical correlation function F , that is symmetric under a permutation of arguments, and the spectral function ρ, that is antisymmetric when permuting the arguments.
We employ a 1/N fluctuation expansion to the real-time effective action Γ 2 to next-to-leading order (NLO) 49,50 . While in the symmetric phase only a single diagram contributes at leading-order (LO), at NLO an infinite series of diagrams has to be summed. The self-energy up to NLO can be schematically represented by the following diagrammatic series The energy density ε(t) exhibits three distinct dynamical regimes: (I) At short times t, single-particle rearrangements lead to a fast increase of the energy density, up to times t int at which interactions become relevant. (II) At intermediate times, t int < t < t th , a stable Floquet prethermal regime occurs in which the interacting system ceases to absorb. (III) At late times, beyond the thermalization time scale t > t th , interactions between a large number of generated quasiparticle excitations cause strong heating. In that regime, the energy density displays an algebraic growth, ε α t t ( ) , with an exponent that approaches α~1/2 for strong interactions. The data is shown for drive frequency Ω = 2.3 and for three different values of the interaction strength u, see legend.
where lines represent full Green's functions G and dots vertices, each of which comes with a factor ~λ/N. In this scheme, the LO (first diagram) is equivalent to a self-consistent Hartree-Fock approximation and thus results in a time-local self-energy that solely renormalizes the bare mass; see methods section. A LO analysis is thus not sufficient to answer the question of whether a prethermal state can be stabilized, as it eliminates the possibility of infinite heating from the beginning. Only at NLO [all other diagrams in Eq. (4)] the self-energy contains parts which are non-local in time and lead to scattering and memory effects that ultimately enable thermalization.
The NLO evolution equations are integrated numerically for times up to 3.18 ⋅ 10 4 driving cycles. The momentum cutoff is set to Λ = π. As initial condition we use the LO groundstate of the O(N)-model for given interaction u and fixed renormalized mass = m 1 eff 2 , i.e. the bare mass m 0 2 gets adjusted accordingly. We have chosen this convention, since the physically relevant observable quantity is the renormalized mass m eff 2 , which has to be fixed to get comparable results. Furthermore, we set the drive amplitude to g = 1/4 and scan the interaction u and drive frequency Ω.

Results
Dynamics of the energy density. The central observable to study heating in any driven system is the where V is the system volume. In our scheme the expectation value of the Hamiltonian is directly available from the Keldysh Green's function F. Calculating the expectation of the quadratic part of the Hamiltonian is straightforward, whereas for the quartic term, we use Heisenberg's equations of motion to express it in terms of higher order time derivatives of the Keldysh Green's function. We obtain Typical plots of ε(t) are shown in Fig. 1. We can divide the heating of the system into three regimes: (I) At short times, up to the interaction timescale t int , the dynamics is dominated by single-particle rearrangements, leading to exponentially fast heating. In that regime, a LO approximation is sufficient to describe the dynamics and scattering of quasi-particles is essentially irrelevant. We define, the interaction timescale t int as the time at which the LO and NLO results starts to deviate, which characterizes the time at which non-local contributions to the self-energy become important. (II) After this initial stage of heating, the system quickly enters a prethermal plateau with low absorption. This Floquet prethermal state persists up to the thermalization time t th and can span several decades in time, thus providing a solid regime for Floquet engineering. (III) At late times  t t th heating becomes significant and we expect the system to approach the infinite temperature state. In that regime our data suggests a power-law growth of the energy density ε α t t ( ) . In the following, we discuss these regimes and where possible provide analytical arguments for the observed behavior. Short time dynamics. At short times, NLO corrections are essentially irrelevant for the dynamics, as confirmed explicitly by comparing LO and NLO results; inset in Fig. 2(a). At LO, the system is equivalent to a multi-dimensional anharmonic oscillator with periodically modulated frequencies (see details in the methods section). We can understand the dynamics in terms of a parametric resonance with the resonance condition set by 2 1/2 is the initial dispersion relation of excitations. The momentum-mode ⁎ p grows exponentially and the fastest growing observable will be γ F t t p e ( , , * ) * t 2 p . Consequently, using (5), the energy density will also grow exponentially in time ε γ t e ( ) * t 2 p . In the Gaussian limit, u = 0, this exponential growth would last indefinitely, but for finite u the self-consistently determined effective mass grows simultaneously with ⁎ F t t p ( , , ), breaking the resonance condition at a certain time, and preventing any further energy-absorption in a LO approximation 34 . Note that since we fix the initial effective mass, , the growth rate γ p is independent of u. This is due to the fact, that the parametric resonance only depends on the frequency and amplitude of the drive as well as the initial quasi-particle spectrum of the system.
Taking into account NLO corrections quasi-particle excitations interact with each other which will eventually lead to heating. We estimate the validity of the LO calculation by  Fig. 2(a). Deviations from the logarithmic scaling exist for  u 1, as in the strong interaction regime NLO processes are important already at initial times, which renders the interaction time scale ill-defined.
In order to validate that the scattering of quasi-particle excitations is the reason for the deviation of the LO and NLO results, we derive a Floquet Fermi's Golden rule (FFGR) 51 , which formally considers NLO diagrams with the lowest number of interaction vertices (sunset diagram); see methods section. We find perfect agreement between the interaction time t int evaluated with the full NLO calculation and the FFGR, respectively, which demonstrates that scattering of created excitations is responsible for the deviations between the leading and next-to-leading order time evolution. This explains why the system can heat up further: Once scattering is taken into account, not only pairs of quasi-particles can be created but the energy can also be distributed over many excitations.
Floquet prethermalization. Once the parametric resonance regime is left, heating becomes extremely slow and the prethermal plateau is entered. In that regime the number of quasi-particles is small and hence the multi-particle scattering, which is enabling further energy absorption, is much slower than pair creation. The number of quasi-particle excitations is directly related to the equal time Keldysh Green's function F, which due to the self-consistent feedback continues to grow. As the thermalization timescale t th is reached, the higher order loop diagrams [ Eq. (4) ] that allow for multi-particle scattering start to dominate. Thus, heating becomes significant and the Floquet prethermal state breaks down.
To quantitatively understand the thermalization time scale t th , we study it as a function of the interaction strength u and driving frequency Ω; Fig. 2(b). The thermalization timescale and thus the lifetime of the prethermal plateau decreases with increasing u, however the functional dependence cannot be described by an exponential or power-law. The dependence is quite strong, with t th changing over one order of magnitude as u varies in the interval [0. 5,15] and Ω = 2.3. Fixing the interaction u, we find that t th decreases with increasing Ω. This is a consequence of all chosen frequencies lying within the initial bandwidth of quasi-particle pairs, <Ω< Λ + 2 2 1 2 , as illustrated in the inset of Fig. 2(b). With increasing Ω, more momentum modes participate in the parametric resonance and consequently the Keldysh component F already ends up being larger as t int is reached; see methods section. Based on our previous arguments on the quasi-particle density, the system thus will be earlier driven out of the prethermal plateau.
Our results do not contradict refs 28-30, which predict that heating is exponentially suppressed at large drive frequencies, as our results are all for small drive frequencies within the two-particle bandwidth. When increasing the drive frequency Ω in our model beyond the two particle bandwidth, the energy absorption becomes very slow. In that case, the system is far away from a parametric resonance and hardly responds to the drive at all.
Even though heating is slow within the prethermal regime t int < t < t th it remains finite and the system does not become fully stationary. Nevertheless, in this regime the Green's function only depends extremely weakly on the stroboscopic center-of-mass time T n = (t + t′ )/2 = 2πn/Ω, where n is an integer. Thus, this extremely slow center-of-mass time dependence should not affect the much faster microscopic processes, that are required to realize novel prethermal states.
Thermalization. At times  t t th , the system is driven out of the prethermal regime and the absorption increases. Our numerical simulations suggest that the energy density grows as a powerlaw ε α t t ( ) (Fig. 1), which can persist for several decades. We show the exponent α as a function of the interaction strength u for different driving frequencies Ω in Fig. 3. With increasing interaction u and drive frequency Ω the exponent approaches 1/2, which appears as a lower bound. In the limit of large u and Ω, the thermalization time scale is smallest and hence, given the fixed maximum time that we can reach in our simulations, the accessible thermalization regime is largest for these parameters. This suggests that the powerlaw exponent might slowly creep to the universal value 1/2 for any interaction u and drive frequency Ω in the asymptotic limit, t → ∞ . In contrast, we found linear heating at late times in the O(N)-model subject to colored noise; see methods section. Moreover, our results suggest that the driven O(N)-model heats to infinite temperature following the well defined prethermal plateau.

Conclusion and Outlook
Our results demonstrate, that a prethermal Floquet state can be stabilized in a periodically-driven quantum many-body system, despite strong interactions and despite the driving frequency being comparable to microscopic energy scales of the system. This opens the possibility of realizing exotic states in the Floquet prethermal regime, such as time crystals or other novel symmetry protected topologically phases. Furthermore, our study suggests a algebraic heating at late times of the form εt t ( ) , which is significantly slower than the linear Joule heating. We attribute this peculiar form of heating to the strong interactions between the dynamically generated quasi particles. How such a sublinear growth can be reconciled with the eigenstate thermalization hypothesis is an important open question. A future study based on a Floquet Boltzmann type approach might provide further insights into this behavior. . The functional Γ 2 [φ, G] is the two-particle irreducible (2PI) effective action, which is given by the sum of all 2PI vacuum diagrams of  . With increasing interaction strength u and driving frequency Ω (but still within the single-particle band), the exponent α quickly approaches 1/2, suggesting that in the asymptotic long-time limit, t → ∞ , the heating rate universally scales as εt t ( ) .

Two-Particle Irreducible Effective Action
G t x t y F t x t y i t t t x t y Using this parametrization, the causal Kadanoff-Baym equations follow directly from Eq. (3). We further make use of the fact, that the driving is uniform in space and that it conserves the O (N)-symmetry, i.e., G ab = Gδ ab to obtain The self-energy is decomposed as . Here, Σ (0) (t), is the local contribution to the self-energy, which leads to a mass renormalization = + Σ m m . By contrast, Σ ′ t t p ( , , ), is nonlocal in time and splits into spectral and statistical components, analogously to the Green's function.symmetric phase, the VEV To make use of Eq. (10), we need the 2PI effective action Γ 2 [G]. As the exact form of Γ 2 [G] is unknown for our interacting model, we employ a large-N approximation scheme 49,50 . At next-to-leading order (NLO), the self-energy becomes

t y I t x t y t x t y I t x t y t x t y N F t x t y I t x t y t x t y I t x t
The functions I F , I ρ are often referred to as summation functions and obey the integral equations . The set of Equations (10)(11)(12) have to be solved simultaneously, starting from t = t′ = 0. To this end, we discretize the system in momentum space and sample 46 points, which we have checked to be large enough to describe the thermodynamic limit. The equations of motion are then integrated numerically using a leap-frog method.

I t x t y t x t y d t I t x t y t x t y dt I t x t y t x t y I t x t y t x t y d t I t x t y t x t y
Dynamics at Leading Order. To leading order the self-energy is time-local and the evolution equations The equations (13) describe coupled anharmonic parametric oscillators (one oscillator for each t′ and p) with initial eigenfrequencies ω = + p p m ( ) (0) 0 2 eff 2 . Let us first discuss the entirely noninteracting case λ = 0, in which Eq. (13) are independent Mathieu equations. It is known from classical mechanics, that the modes satisfying the resonance condition 2ω 0 (p) = nΩ with n = 1, 2, … experience a parametric resonance and will grow exponentially in time. As there is no feedback on the spectrum of the system for λ = 0 this exponential growth in the resonant modes continues forever.
For finite λ, the exponential growth of the statistical correlation function F(t, t′ , p) for momenta p satisfying the resonance condition leads to an exponential growth of the effective mass  34 . Therefore, the quasi-particle bandwidth will at a certain time lie entirely in between the parametric resonances and the system cannot absorb energy anymore, Fig. 4.
Scientific RepoRts | 7:45382 | DOI: 10.1038/srep45382 The failure of the system to absorb further energy can be traced back to the fact, that the LO self-energy is local in time and only leads to a renormalization of the quasi-particle dispersion. Except for this renormalization the quasi-particles remain sharp excitations and there is no mechanism present, that allows energy-transfer between them. Consequently, there is only energy absorption from the drive when the driving frequency hits the sharp resonance for the creation of quasi-particle pairs and the heating stops as soon as the resonance condition cannot be fulfilled anymore.
Floquet Fermi's Golden Rule. The simplest diagram leading to scattering between quasi-particles is the "sunset" diagram, see Fig. 5. This diagrams includes interactions of only two quasi-particles. By contrast, higher loop diagrams would include scattering events of more than two particles. As we discuss in the main text, these higher-order events become relevant only at later times.
We obtain the following expression for the FFGR self-energy  Expressing Σ FFGR in this way, we see that the Floquet Fermi's golden rule analysis amounts to replacing the summation function I in the expression for the NLO self-energy, Eq. (11), with the polarization bubble Π .
We calculate the energy-density of the system resulting from the FFGR and find very good agreement with the NLO results for t int , see Fig. 2(a). However, we emphasize that the FFGR self-energy, Eq. (15), does not correspond to a conserving expansion of the 2PI effective action Γ 2 and hence is bound to fail for long times, as it eventually becomes divergent. Figure 4. Illustration of the parametric resonance. When a part of the initial bandwidth (red shaded region) for quasi-particle pairs lies in a parametric resonance (blue shaded regions), the corresponding modes grow exponentially, which leads to an exponentially growth of the effective mass m eff (t). Accordingly the dispersion relation ω t (p) is shifted toward higher energies (light red shaded region). Once the quasi-particle bandwidth lies completely in between two regions of parametric resonance, the system stops absorbing energy in the leading order approximation. The parametric resonance at nΩ, n = 1, 2, … is only sharp in the limit of vanishing drive amplitude A and smears out with increasing A. For small A, the width of the resonance grows linearly with A. Figure 5. Floquet Fermi's golden rule. The "sunset" diagram takes into account scattering of two quasiparticles beyond a simple renormalization of the quasi-particle mass. It is non-local in time and leads to memory effects.
Introducing noise ξ(t) is expected to mimic, at least very crudely, the effect of scattering. Therefore, the system is expected to heat to infinite temperature even with the leading order self-energy.
We explore two cases for the random process, which are white noise and correlated noise, respectively. In the case of white noise, ξ w (t) reduce to Gaussian random variables with vanishing mean, 〈 ξ w (t)〉 = 0 and auto-correlation 〈 ξ w (t)ξ w (t′ )〉 = γ 2 δ(t − t′ ). By contrast, the correlated noise ξ c (t) obeys the stochastic differential equation of the Ornstein-Uhlenbeck process and 〈 ξ c (t)〉 = 0. Note that white noise is recovered in the limit τ → 0, σ → ∞ , keeping στ = γ fixed. White noise is completely uncorrelated, while colored noise has exponentially decaying correlations in time. Hence, one expects that the system thermalizes faster when it is subject to white noise. This is what we find by numerically solving Eq. (13). Moreover, we find that the energy-density grows according to a power-law ε α t t ( ) ; Fig. 6. We exploit the similarity of Eq. (13) and an anharmonic oscillator, for which it has been shown that the energy grows quadratically in time for white noise (α = 2), whereas colored noise leads to a linear growth ε~t (α = 1) 52,53 . The dynamical evolution in our system, Eq. (16), confirms these expectations. Therefore, the heating due to either white (ε~t 2 ) or colored (ε~t) noise is substantially faster than the asymptotic heating we observe when solving the equations of motion self-consistently up to NLO (ε~t ). We attribute the slow heating obtained from the full solution up to NLO to the strong interactions between quasi-particles which cannot be simply mimicked by multiplicative noise. With multiplicative noise the system absorbs energy indefinitely even at leading order. The energy density grows with a powerlaw ε α t t ( ) for late times for arbitrary system parameters. Depending on whether the noise is white, i.e., completely uncorrelated in time or colored, i.e., correlated in time, the growth is quadratic or linear, respectively. From that we deduce that correlations in time slow down heating in the system. The data is shown for driving frequency Ω = 2.3 and interaction strength u = 1.0. The strength of the white noise is γ = 2.0/Ω 2 whereas for colored noise we have chosen σ = 2/Ω 2 and τ = 20/Ω.