Thermalization and Prethermalization in Periodically Kicked Quantum Spin Chains

We study the dynamics of periodically-kicked many-body systems away from the high-frequency regime, and discuss a family of Floquet systems where the notion of prethermalization can be naturally extended to intermediate and low driving frequencies. We investigate numerically the dynamics of both integrable and nonintegrable systems, and report on the formation of a long-lived prethermal plateau, akin to the high-frequency limit, where the system thermalizes with respect to an effective Hamiltonian captured by the inverse-frequency expansion (IFE). Unlike the high-frequency regime, we find that the relevant heating times are model dependent: we analyze the stability of the prethermal plateau to small perturbations in the drive period, and show that, in a spin chain whose IFE is intractable, the plateau duration is insensitive to the perturbation strength, in contrast to a chain where the IFE admits the resummation of an entire subseries. Infinitesimal perturbations are enough to restore the ergodic properties of the system, and decrease residual finite-size effects. Although the regime where the Floquet system leaves the prethermal plateau and starts heating up to infinite temperature is not captured by the IFE, we provide evidence that the evolved subsystem is described well by a thermal state w.r.t.~the IFE Hamiltonian, with a gradually changing temperature, in accord with the Eigenstate Thermalization Hypothesis.

The cornerstone of the theory of periodically-driven systems, modeled by a Hamiltonian H(t) = H(t + T ), is Floquet's theorem. It states that, at times integermultiple of the drive period (i.e., stroboscopically), the evolution operator U ( T, 0) is generated by the time- * christoph.fleckenstein@physik.uni-wuerzburg.de † mgbukov@phys.uni-sofia.bg independent Floquet Hamiltonian H F : This comes in stark contrast to generic time-dependent Hamiltonians, where the evolution operator is defined through a complicated time-ordered exponential. Yet, periodically driven systems do not conserve energy. Theoretically, Floquet systems are of conceptual importance, since they feature a nontrivial controllable equilibrium limit: at infinite drive frequencies, energy conservation is restored, and the Floquet Hamiltonian is a static local operator whose properties are indistinguishable from those of static many-body systems. Hence, Floquet systems provide a systematic approach to understand and analyze nonequilibrium behavior. Away from the infinite-frequency limit, energy absorption may occur, and generic local many-body Floquet systems are currently believed to heat up to infinite temperature at infinite times [35][36][37][38] [but see also Refs. [39][40][41][42][43]].
The existence of the infinite-frequency limit affects significantly the dynamics of Floquet systems: when the drive frequency is much larger than the typical single-particle energy scales of the non-driven Hamiltonian, following a quick constrained thermalization stage, fast-driven systems enter an exponentially long-lived prethermal plateau, before unconstrained thermalization brings the system to a featureless infinite-temperature state [44,45].
The physics in the prethermal plateau is well captured by the inverse-frequency expansion (IFE) [4] for the effective approximate local Hamiltonian H eff ≈ H F [46]. This equilibrium-like regime facilitates significantly the analysis of Floquet systems. Moreover, it provides a play-arXiv:2101.04372v1 [cond-mat.stat-mech] 12 Jan 2021 ground for the ideas of Floquet engineering. It opens up a long prethermal time window which, under suitable conditions, supports phases of matter inaccessible in static systems [15,42,[47][48][49][50][51]. Curiously, the physics of the prethermal plateau has been found to exist in isolated (semi-)classical Floquet systems, which suggests that it is not caused by quantum mechanical processes [52][53][54][55][56][57][58].
Recently, it was shown that a similar prethermal plateau exists for random dipolar driving [59], the periodicallydriven SYK model [60], and for quasi-periodically driven systems where its duration is controlled by a stretched exponential [50,59,61]; the latter have also been shown to exhibit topological phenomena [62,63].
At lower drive frequencies, the system starts absorbing increased amounts of energy via a proliferation of Floquet many-body resonances [64]. In order for a Floquet system to absorb energy from the periodic drive, two conditions must be met: (i) the existence of many-body eigenstates in the non-driven system whose energies differ by an integer multiple of the drive frequency (the socalled spectrum folding criterion), and (ii) a finite transition matrix element between these states when exposed to the periodic drive. The prethermal plateau shrinks gradually with decreasing the drive frequency until it disappears completely when the drive frequency becomes of the order of the single-particle energy scales in the non-driven Hamiltonian. This is correlated with a progressively more nonlocal operator structure of the exact Floquet 'Hamiltonian', whose inverse-frequency approximation breaks down as an asymptotic series with the onset of infinite-temperature heating [45,64]. Recently, techniques have been developed to find approximations to the Floquet Hamiltonian in the intermediate and lowfrequency regimes, based on the empty-lattice-type approximation [65], the Flow equation approach [66,67], Floquet perturbation theory [68], and the Replica expansion [69].
In this paper, we discuss in detail an extension of the notion of prethermalization to the intermediate and lowfrequency regime, introduced in Ref. [70]. We investigate various step-driven integrable and nonintegrable Hamiltonians in the vicinity of commensurate driving periods T * k , for which energy conservation is restored exactly. We also provide detailed numerical evidence for the prethermal behaviour of observable quantities.
We find a rich thermalization behavior: unconstrained thermalization to infinite temperature is suppressed with drive-dependent heating rates following both powerlaw and non-powerlaw behavior as a function of the distance ε to the commensurate point T * k . We believe that this suppressed thermalization originates from the suppressed magnitude of matrix elements between resonant manybody states (the folding criterion being readily satisfied close to T * k ). The intermediate-to-low frequency regime enhances the ergodic properties of the dynamics, and allows us to obtain clean data already at moderate system sizes. Moreover, the heating behaviour caused by the Floquet drive is resilient against small perturbations in the driving period. We further consider the evolution of both pure states and thermal ensembles, and demonstrate that prethermalized Floquet systems evolve into a featureless infinite-temperature state by gradually changing their temperature with respect to a local effective Floquet Hamiltonian, although the latter is computed with the help of an asymptotic IFE. Thus, the IFE can provide a useful static description even outside the prethermal plateau.
The paper is structured as follows. In Sec. II, we introduce a class of Floquet systems which later on allows us to extend the notion of Floquet prethermalization to intermediate and low drive frequencies. In Sec. III, we analyze the thermalization dynamics generated by a generic nonintegrable Hamiltonian -the driven mixed-field Ising model -starting from a pure initial state (Secs. III A and III B); we define and discuss the qualitative behavior of heating rates (Sec. III C), as well as its robustness to perturbations in the driving protocol (Sec. III D). Next, in Sec. III E, we investigate the prethermal properties, starting from a thermal initial ensemble. In the second part of the study (Sec. IV), we investigate Floquet dynamics generated by two integrable Hamiltonians: the transverse-field Ising model , and the Ising model without quantum fluctuations. Finally, in Sec. V we conclude and summarize our results. Additional data, including finite-size scaling, and the complete replica derivation of the effective Hamiltonian, are shown in the Appendix.

II. REALIZING PRETHERMAL BEHAVIOR AWAY FROM THE HIGH-FREQUENCY LIMIT
Following Ref. [70], we consider the family of Floquet unitaries: where T = 2π/Ω is the drive period with the associated frequency of switching (henceforth called the drive frequency). The operator V is required to have a commensurate spectrum; H is an arbitrary local many-body Hamiltonian, such that the average Hamiltonian H ave = H + V is nonintegrable (i.e. does not possess an extensive number of local conserved integrals of motion). Since H ave is the leading-order term in the IFE, we assume that the nonintegrability of H ave implies the nonintegrability of H F [71]. Note that the commensurability condition on V is not excessively restrictive, since merely all shortrange interaction terms of density-density type readily satisfy it. For simplicity, in this paper we consider spin-1/2 systems, where is a global magnetic field of strength γ along the x-direction; the Pauli matrices obey [σ α i , σ β j ] = 2iδ ij αβγ σ γ j . We shall discuss both integrable and nonintegrable drive Hamiltonians H.
The choice of a symmetric drive in Eq. (2): {T /4, T /2, T /4}, as compared to {T /2, T /2}, results in a real-valued generator H F of stroboscopic dynamics. We verified that it does not affect our results and conclusions. We mention in passing that, although Eq. (2) bears a formal resemblance with Floquet time crystals [47][48][49]72], investigating time-crystalline behavior is beyond the scope of the present study.
Due to the commensurate structure in the spectrum of V , there exists a sequence of drive periods T = 2T * k = 2πk/γ with k ∈ N (for L even and V from Eq. (3)), and associated frequencies Thus, for the class of Floquet unitaries under consideration, we have Hence, at the special points T * k , the dynamics of the kicked system in Eq. (2), reduces to a quench problem to the static local Hamiltonian H. Therefore, by construction, at T = 2T * k , energy is conserved and the system is prevented from heating up. In this work, we discuss the behavior of this class of systems in the vicinity of T * k . Note that, for k = 0, we recover the familiar infinitefrequency point, surrounded by an interval of large but finite frequencies for which the Floquet system exhibits a prethermal plateau [44,45]. Thus, the setup in Eq. (2) provides a natural candidate to extend the Floquet prethermal physics to finite frequencies. Observe that, for k > 0, the drive frequency Ω can also happen to be in the intermediate-to-low frequency regime [compared to a typical single-particle energy scale in H]. Therefore, the present construction allows to induce stable isolated points Ω * k on the frequency axis, by means of inhibiting Floquet resonances.
In this study, we focus on the vicinity of the stable points T = 2(T * k + ε) (ε T * k ) where resonances are expected to be suppressed, and The problem reduces to that of the system H subject to small periodic kicks V of strength ε [39,40].
To get an intuition for why prethermal behavior can be expected even at low frequencies in this class of systems, consider the (oversimplied) model H 0 = 2γ j σ z j+1 σ z j + σ z . Notice that, in this case, the spectra of both V and H 0 are commensurate with the same T * k ; yet, H ave = H 0 +V is the mixed-field Ising model which is a nonintegrable Hamiltonian; hence, the Floquet system is expected to display thermalizing dynamics and heat up at intermediate to low frequencies. However, for this choice of H 0 and V , it is easy to see that U F (2(T * k +ε)) = U F (ε) for all k, and the period axis compactifies to a circle. Therefore, despite T * k corresponding to a low drive frequency Ω at large k, the behavior of the system around higher-order commensurate points T * k is exactly the same  [73]. Both panels show the formation of a prethermal plateau over a few decades of driving cycles, whose duration increases parametrically out to infinity as ε → 0 at the commensurate point T * k . The two dashed horizontal black lines correspond to the right-hand-side in Eq. (12) for O = Q and O = Sent, respectively. The purple dashed curve in (a) highlights one curve to better compare it to its counterpart shown in Fig. 6. We choose 15 logarithmically spaced ε values, ε ∈ [3 × 10 −4 , 3 × 10 −1 ] (the interval limits including). The parameters are hz/J = 0.809, hx/J = 0.9045, γ/J = 1, and k = 2; the frequency of switching is Ω * 2 /J = 1/2. The system and subsystem sizes are L = 20 and LA = 10, respectively. We show both quantities against + 1 to make the initial value = 0 visible on the log scale. We display a logarithmically decreasing number of data points at large . as around infinite frequency (i.e. k = 0). In particular, it follows that the dynamics features an exponentially long prethermal plateau for T ≈ 2T * k , while energy conservation is restored exactly (and thus the plateau lifetime becomes infinite) for T = 2T * k . This toy model showcases that, in order for a Floquet system to heat up to infinite temperature, it must have a finite matrix element between the states of the non-driven Hamiltonian whose energies differ by integer multiples of Ω; in other words, the criterion for folding the spectrum is a necessary but not a sufficient condition.
Throughout this paper, we consider integrable and nonintegrable Hamiltonians H j with non-commensurate spectra, where the description of the behavior around T * k is not immediately obvious. Specifically, we attempt to answer the following questions: (i) Under what conditions can there exist a prethermal plateau in the vicinity of T * k ? Notice that for classical systems, Nekhoroshev's estimate w.r.t. breaking energy conservation in the vicinity of T * k postulates that integrals of motion are conserved up to exponentially long times in ε −1 [74]. However, these estimates carry a system-size dependence and, to the best of our knowledge, there is no formal proof which holds in the thermodynamic limit. Quantum mechanically, Fermi's Golden Rule (FGR) postulates that the system should start absorbing energy for infinitesimally small ε, but we do not have expressions for how the magnitude of the transition matrix elements depends on ε [for some models, recent results indicate that the latter is captured by ETH [75]]. (ii) What are qualitative differences between the prethermal plateaus at k = 0 (infinite frequency) and k > 0 (moderate to low frequencies)? (iii) Is there an effective approximate analytical description for the dynamics of the system in the vicinity of T * k , similar to the IFE? (iv) Does the thermalization dynamics depend on whether H is integrable or nonintegrable? (v) Is the state of the Floquet system, after the initial transient is over, fully thermal, or are there any drive-induced synchronization effects [54]?
Along the way, we also investigate the following hypothesis: if the thermalization dynamics of a pure state subject to a Floquet drive exhibits a prethermal plateau, then the subsequent approach to the infinite temperature state, caused by unconstrained thermalization, is a quasistatic process; in particular, a subsystem remains in a thermal state of gradually changing temperature. However, when the system heats up to infinite temperature without going through a prethermal plateau, equilibration is first reached at energy densities corresponding to infinite temperature.

III. NONINTEGRABLE DRIVES
Consider first the kicked system (2), with V is given by Eq. (3), for the nonintegrable spin-1/2 mixed-field Ising model The high-frequency (k = 0) behavior in Floquet systems is distinguished by a long-lived prethermal plateau, and our first goal is to investigate the behavior of the kicked Floquet system close to the commensurate point T * k for k > 0. To this end, we prepare the system in the domain wall state P| ↑ . . . ↑↓ . . . ↓ , where P is the projector onto the zero-momentum sector of positive parity. Ordered pure states are of particular importance in view of recent progress in Floquet engineering [1][2][3][4]. That said, we verified that the conclusions laid out below, do not depend on the choice of the initial pure state [although thermalization and equilibration timescales typically do].
We compute the exact evolution of the system numerically up to 5 × 10 4 driving cycles, and do measurements of the energy density E( ) = ψ( )|H ave |ψ( ) /L in the time-evolved state |ψ( ) = U F |ψ i at stroboscopic times T . Let us define the rescaled energy where H ave β=0 ≈ 0 is the infinite-temperature expectation value of the average Hamiltonian. The quantity Q( ) measures energy absorption relative to the energy of the initial state.  Figure 1a shows that a qualitatively very similar behavior to the familiar infinite-frequency point (k = 0), occurs in the vicinity of the commensurate points T * k with k > 0 [cf. also App. A 2]. Because, k > 0 falls in the low-frequency driving regime, one can potentially make use of this parametrically long-lived stable regime to extend ideas from Floquet engineering to the low-frequency regime.
In the limit ε → 0, H ave is close to the exact Floquet Hamiltonian H F which is conserved. Therefore, to guarantee that the observed prethermal behavior is not a property of the energy observable H ave , we also show the time evolution of the entanglement entropy density. Denoting a chain subsystem by A, and the corresponding reduced density matrix by ρ A = trĀ|ψ( ) ψ( )|, we define As anticipated, the prethermal plateau is also clearly visible in Fig. 1b. In particular, we observe the same four stages of thermalization in the vicinity of the commensurate points, familiar from the high-frequency regime: (I) a transient of constrained thermalization precedes (II) a prethermal plateau, followed by (III) a second transient of unconstrained thermalization leading eventually to (IV) a featureless infinite-temperature state. Observing the prethermal physics in S ent , we anticipate that this behavior is generic, i.e. it applies to all local observables; we confirm this numerically in Fig. 2.

B. Local Equilibration and Subsystem Thermalization
Consider the quench problem of preparing a system in some initial state, and then evolving it under a generic Hamiltonian. A defining prediction of ETH is that a subsystem, evolving under a nonintegrable Hamiltonian, thermalizes at a temperature, corresponding to the energy density of the initial state [76,77]. In short, the reduced density matrix ρ A is expected to evolve into the thermal state ρ A th [78,79]. In Floquet systems, it has been established that the prethermal plateau around k = 0 is well described by an effective Hamiltonian H eff obtained using the IFE [41,44,45,54,64]. We now study numerically the applicability of ETH in the vicinity of commensurate points T * k for k > 0. Our objective is to investigate whether, under unitary evolution of the full system, subsystem A evolves into the mixed Gibbs state Here H A eff is the effective Hamiltonian restricted to subsystem A, and β(E i ) is the temperature, corresponding to the energy density of the initial state.
Because we do not have a handy analytical expression for H eff in the nonintegrable mixed-field Ising model, Eq. (6), [cf. Sec. IV B for a model amenable to the IFE], we work to leading order in ε; coincidentally, this provides a sufficient description for the range of ε values that exhibit a prethermal plateau. Thus, close to T * k , we have H eff ≈ H ave + O(ε) which is nonintegrable by construction, and hence we expect ETH to apply w.r.t. H eff . One can, of course, add higher-order corrections whenever they can be computed [cf. Sec. IV B]. We demonstrate the applicability of ETH in two steps. First, we compute the energy density of the initial pure state E i = ψ i |H eff |ψ i /L, defined on the full system of L sites. We can associate an inverse temperature β(E i ) to the initial energy density, by solving the implicit equation for β. This provides us with a theoretically predicted reference value for the inverse temperature of the prethermal plateau. Independently, as a second step, we also extract a value for β from our exact numerical simulations. To do this, we first construct an approximation to the density matrix of the diagonal ensemble ρ d [80] empirically from the time series of the evolved state |ψ( ) : where m = 1, . . . , M are consecutive stroboscopic times. Note that, unlike the exact definition ρ d = n | ψ i |n F | 2 |n F n F |, (i) the empirical definition in Eq. (10) gives the diagonal density matrix in the computational basis and hence it does not involve/require knowledge of the exact Floquet eigenstates |n F . Moreover, (ii) ensemble averages using ρ d ( ) correspond in a natural way to experimental measurements in the system [81]. (iii) the time-or -dependence of ρ d ( ) allows us to monitor the time evolution of the diagonal ensemble. In practice, we use M = 20 consecutive stroboscopic states to do the time-average, but one should be careful that the system does not deviate from its steady state physically during this window, e.g., by monitoring the values of local observables.
To extract a numerical value for β, we first compute the reduced diagonal ensemble density matrix After that, we plot the spectrum of ρ A d against the spec- Using a least-squares fit to extract the slope of this straight line gives the temperature of the system in the thermal state. Figure 3a indicates that, starting from a pure state on the full system, the subsystem evolves into a thermal state, whose temperature matches well the value predicted by ETH w.r.t. H eff . In particular, for ε 10 −3 , the long-lived prethermal plateau appears to be well described by a thermal density matrix with inverse temperature β(E i ). Likewise, Fig. 3b shows the dynamically extracted values for the inverse temperature β as a function of the energy conservation breaking parameter ε; the error bars show the least square fit uncertainty. This curve depends on the time at which the diagonal ensemble is constructed, since all states are expected to reach infinite temperature at sufficiently long times in the thermodynamics limit.
Our data is fully consistent with ETH predictions for the prethermal plateau [ Fig. 3b, solid orange line]; however, it contains more information. Fig. 3c shows the dynamically extracted values of β at different times during the evolution. The dashed line marks the solution to Eq. (9), where we replaced E i by its value at a later time E( ). Although heating processes cause the system to leave the prethermal plateau, the state of the system at subsequent times is still well-described (to a good approximation) by a thermal state w.r.t. the approximate H eff [dashed lines in Fig. 3(b-c)]. Notice that, although thermal states are universal, in the sense that they maximize the thermodynamic entropy, at a finite temperature they are only well-defined if the Hamiltonian is known, w.r.t. which the state is thermal; this is highly non-trivial in time-dependent systems.
In Fig. 4 we demonstrate that the observed behavior is generic: we present the same comparison between ETHpredicted and fitted inverse temperature but for a few different observables; this is performed using Eq. (9) and replacing H A eff and E( ) with the corresponding local observable and its expectation value, respectively. Initially, as the system is not thermal, large deviations appear between ETH predicted and fitted value for non-energy related quantities. Yet, as soon as the system evolves into a thermal state, ETH predicts the expectation value of local observables, given their instantaneous expectation value.
The above finding may come as a surprise, since the heating processes that drive the system out of the prethermal state are the same which cause the failure of the IFE to converge, and which have been shown to arise from non-analytic terms (in Ω −1 ) present in H F but not in H eff to any order [64]. We find that, although the IFE fails to predict the exact value of the energy density E( ) in 'unconstrained thermalization' stage (III) of the dynamics, given E( ) and H eff one can reconstruct the thermal state that characterizes the system at that point of time. This is reminiscent of the observation that the IFE describes well ensemble expectation values in classical many-body Floquet systems, but not the precise dynamics of observables in isolated evolved configurations (due to classical chaos) [55]. This result is remarkable, because it hints at the existence of a simple hydrodynamic effective description for the dynamics of closed many-body Floquet systems all the way up to the infinite-temperature state at sufficiently long times [82]; attempts have already been made in static open systems [83,84]; recently, experimental protocols to measure temperature in systems undergoing a slowlychanging equilibrium were also proposed [85]. To the best of our knowledge, the law that governs the timedependence of β( ) is currently unknown.
Finally, note that the data in Fig. 3c provides numerical evidence in favor of the Hypothesis we stated at the end of Sec. II.

C. Qualitative Heating Rates
We now turn our attention to the heating rates in the vicinity of the commensurate points T * k . In generic Floquet systems, heating in the vicinity of the infinitefrequency point (k = 0) is exponentially suppressed in the drive frequency for both classical and quantum systems [45,54,55,86]. In contrast, here we show that for We define the heating rate Γ(ε) empirically, as the inverse time at which the value of an observable (or the entanglement entropy) drops to half of its prethermal plateau value. Conversely, we call Γ −1 that heating 'time', which solves the equation This definition allows us to extract the ε-dependence of Γ −1 from the numerical data, up to a pre-factor which depends on the model parameters. Figure 5 demonstrates a power-law scaling Γ −1 ∝ ε −α of the heating rates for k = 2. The data is fully consistent with applying Fermi's Golden Rule (FGR) to the periodically-kicked problem (5), which predicts α = 2 [75]. This represents a major difference compared to the k = 0 point, where Γ −1 ∝ exp(ε/ξ) [54]. Thus, heating close to T * k is only power-law suppressed for H 1 , which explains the relatively small values of ε required for a prethermal plateau to form. Note that the power-law scaling of the heating rate is universally seen in the dynamics of all observables and the entanglement entropy. We also checked that this behavior appears for all k > 0, not just k = 2 (see App. A 2).
We mention in passing that observing an exponentially-suppressed heating for k > 0 is likely possible for H 1 if γ/J 1 is large enough. This becomes plausible when the first commensurate point T * 1 at k = 1 corresponds to a drive frequency larger than the single-particle energy scale of the problem. In this case, however, the system falls outside the low-frequency driving regime for k = 1. Moreover, even in such a case, for a large enough k, we expect a power-law scaling of the heating rates.

D. Robustness to Drive Noise
While the observed quadratic scaling provides evidence that FGR underlies the heating behavior for k > 0, our simulations show that the dynamics of the periodically kicked system may not be fully ergodic out to very long times. This can be seen by noticing that for some values of ε (e.g., ε = 0.016, 0.025) the curves showing the time-dependence of observables get stuck before reaching their infinite temperature values [ Fig. (1)]. Moreover, the presence of a similar feature in the entanglement entropy curves [ Fig. (1)b], which by ETH is related to the thermal entropy, suggests that the system does not explore ergodically the entire available Hilbert space. This secondary plateau occurs at high energy densities long after the system has left the prethermal plateau. The phenomenon appears in the behavior of merely all quantities of interest, and is puzzling because H eff is a completely ergodic, nonintegrable Hamiltonian. In time-independent systems, lack of ergodicity typically suggests the existence of hidden (left-over) conservation laws; these are, however, ruled out both for H eff and for the exact Floquet Hamiltonian generated by Eq. (6). Therefore, we look for an explanation related to the nonequilibrium dynamics of the system.
To investigate this non-ergodic feature in detail, we perturb the periodicity of the drive: we keep the strength ε of the small kick fixed, while adding a small random number δ ∈ [0, 0.05T ] to the duration T /2 the Hamiltonian H is applied for: We consider the regime of small perturbations δ T /2, irrespective of the value of ε [which itself is a perturbation around T * k that controls breaking of energy conservation]. Since this procedure destroys the perfect periodicity of the Floquet drive, any drive-induced synchronization effects [54] would be destroyed as well. Therefore, by comparing the perturbed and perturbation-free Floquet dynamics, we can infer whether synchronization effects occur in our system. This is intimately related to the Markovian properties of the Floquet dynamics which tells if the latter retains memory of its evolution.
In Fig. 6a we show the time evolution of the energy in the kicked system subject to small perturbation strength δ/T = 0.005. Comparing the curves to the perturbationfree case [compare especially the colored dashed lines in Fig. 6a and Fig. 1a], we clearly see that the random perturbation in the drive period helps restore ergodicity: all curves in the noise-perturbed dynamics approach the infinite-temperature value at sufficiently long times. Moreover, we also find that adding the noiseperturbation does not change the time it takes for the system to leave the prethermal plateau: in Fig. 6b, we show the heating time curves Γ −1 (ε).
Finite δ breaks the periodicity of the drive, and hence H eff changes from period to period. Naively, one may render Floquet theory inapplicable. However, for kicked systems we can still apply the more general Baker-Campbell-Hausdorff formula. Hence, in the present case, where the thermalization dynamics is mainly driven by the leading order H eff = H 1 /2 + O(ε), the finite perturbation strength results in an additive correction of the order δ/(2T ): As a result, small value of δ/T have negligible effects on the observed prethermalization. The more prominent effect of adding the perturbation is that all crossover values of ε shift to align perfectly on the straight line with increasing δ/T [ Fig. 6b].
The simulation data demonstrates that the periodically-driven system (δ = 0) is not fully ergodic at the finite system sizes L within the reach of reliable simulations. Finite-size scaling indicates that ergodicity is restored as the system size approaches the thermodynamic limit even in the periodically-driven system [App. A 1]. Yet, at finite system size, adding the perturbation provides a useful technique to simulate ergodic behavior. It is currently an open question what mechanism causes this drive-induced synchronization at long times for finite system sizes, and whether this can be interpreted as a collective phenomenon induced by the Floquet drive in a finite-size system.

E. Thermal Initial Ensemble and Dependence on the Energy Density of the Initial State
The discussion on ergodicity in Sec. III D raises the question whether the state of the periodically-driven system in the prethermal plateau is fully thermal w.r.t. H eff . While the results in Sec. III B already provide a strong indication for this claim, they do show small deviations from the expected thermal behavior. To settle this question, and to show that the results from the previous sections are not sensitive to the energy density of the initial state, we simulate the dynamics of a thermal ensemble, and compare the behavior of the time-evolved thermal state to that of the evolved pure state [Sec. III B].
Simulating exactly the dynamics of a thermal ensemble amounts to solving the von Neumann equation for the density matrix of the full system, starting from a thermal initial state. Unfortunately, with the computational power at our disposal, this proves to be infeasible for spin chains of size L = 20. The reasons for this are the exponentially large (in L) Hilbert space size, and the long evolution times required in our study. Therefore, we resort to an approximate approach, based on quantum typicality [87][88][89][90]. Typicality, which is unrelated to integrability, states that the trace of an operator O defined on a Hilbert space H can be approximated as where |r n are Haar-random states. The approximation becomes exact in the limit N → ∞. Hence, thermal expectation values w.r.t. H eff in the full system, at temperature β −1 , can be approximated as [91] O Interpreting the expression on the right-hand side as an ensemble average, the thermal density matrix is approximately equal to where the subscript β denotes the inverse temperature of the thermal state. Notice that this definition requires N pure states |ψ n = e − β 2 H eff |r n . Therefore, to compute the time evolution of the thermal ensemble, by linearity of the ensemble average, it suffices to evolve each state |ψ n separately and then build: Note the difference of this approximate thermal ensemble ρ β ( ) to the empirical diagonal ensemble ρ d ( ) we introduced in Eq. (10): we construct the diagonal ensemble out of a time series of quantum states [at sufficiently long times when the initial transients have died out], starting from a single initial state. In contrast, in the approximate thermal ensemble, we have a set of initial states which we evolve up to some time before taking a measurement. While the two ensembles may seem different, for dynamics governed by nonintegrable Hamiltonians, they become equivalent in the thermodynamic limit: in fact, it is within the sense of the diagonal ensemble, that thermal expectation values, as defined in statistical mechanics, are to be carried out in practice [since experimentalists typically do not have many copies of the many-body system to build a proper statistical ensemble].
The random states |r n are defined in the full Hilbert space; in practice, we decompose the simulation over various symmetry sectors of H eff for efficiency. In order to avoid building and diagonalizing the 2 L × 2 L matrix ρ β , we first reduce the evolved pure states to subsystem A. Noting that the partial trace and the ensemble average are mutually commuting linear operations, we obtain This computation of ρ A β ( ) can be trivially parallelized over the ensemble to gain speed. Figure 7a shows the time evolution of the energy density of the time-evolved approximate thermal ensemble. Because the system starts already in a thermal state w.r.t. H eff , its dynamics does not feature the initial constrained thermalization stage (I), in contrast to starting from a pure state. Repeating the simulation for a few different initial inverse temperatures β, we see that the Floquet dynamics is insensitive to the energy density of the initial state (provided the assumptions of ETH for H eff are satisfied, see inset in Fig. 7a). Figure 7b shows snapshots of the spectrum of ρ A β ( ) at a few different times . Once initialized in a thermal state, the system remains thermal throughout the time evolution. This is a trivial observation but it rules out the possibility for the system to enter a nonequilibrium state during the evolution before reaching infinitetemperature at long times; instead, we see that the state can be described by a thermal ensemble with a slowly varying temperature. This provides additional evidence for the Hypothesis laid out in Sec. II: indeed, irrespective of the energy density of the initial state, once a nonintegrable Floquet system enters a thermal state, it will remain thermal under continued exposure to the periodic drive; its temperature increases slowly as the system heats up to infinite temperature. Similar to Sec. III B, we find once again that the state of the system during the evolution is (approximately) thermal w.r.t. H eff = H F even at times past the prethermal plateau; this confirms that, although H eff is insufficient to capture the heating dynamics of the system, given the energy density of the evolved state at some time , H eff = H F contains the necessary information to describe the thermal state the Floquet system is in, at any point during the evolution. However, what H eff misses, are the very processes that cause the system to heat up in the first place.
Finally, in Figure 7c we show the evolution of the inverse temperature β, and compare it to the theoretical prediction according to Eq. (9). We attribute the mismatch at short times to finite-size effects, finite-ε corrections to the effective Hamiltonian H eff , and to the relatively small number of states N = 20 used for the ensemble average. Moreover, note also that the boundary conditions that we select for the effective Hamiltonian of the subsystem have an influence on the ETH-predicted value for β. In fact, despite being physically unrealistic, periodic boundary conditions (PBCs) seem to describe better the fitted values of β as compared to (the more natural) open boundary conditions (OBCs). However, such differences are not expected to persist in the thermodynamic limit. This can be already observed on the scale of systemsizes we investigate (see App. A, Fig. 13).
Based on the data in Fig 7, we conclude that, besides the initial constrained thermalization transient, there is no difference (within the limits of finite-size simulations) in the later stages of the Floquet evolution of a pure, as compared to a thermal initial state. In both cases, we observe a qualitatively and quantitatively similar behavior.

A. Transverse-field Ising Model
Let us now turn to integrable models. An integrable limit of the Hamiltonian H 1 is given by the transversefield Ising model where we set h z /J = 0.9045. We drive the system according to the protocol of Eqs. (2) and (3). The initial state is the domain wall pure state in the z-basis, projected to the zero momentum sector of positive parity. The corresponding effective Hamiltonian to leading order in ε is . Despite the integrability of H 2 , and similar to H 1 , it is infeasible to obtain a closed-form analytical expression for the higher-order correction terms to the Floquet Hamiltonian.
Investigating the heating behaviour of Eq. (19) around T * k is particularly interesting from two perspectives: (i) unlike the nonintegrable Ising model, where ε breaks only the remaining energy conservation law, here the same ε also breaks integrability. Recently, it was proposed that quantum chaotic behavior, set out by infinitesimal integrability breaking, can be sensitively detected using adiabatic gauge potentials [92]. Exactly at T * k , integrability is restored and, in the vicinity of these points, we can study how the integrability breaking parameter ε influences the thermalizing dynamics. (ii) despite being an integrable model, H 2 possesses a non-commensurate spectrum, which can be found by virtue of the Jordan-Wigner mapping to free fermions. This is a prerequisite for the proliferation of resonances, once the drive is turned on; away from T * k , resonances are expected to facilitate thermalization.
In general, integrable models do not obey ETH. Instead, they are believed to thermalize to a generalized Gibbs ensemble with a Lagrange multiplier associated to each conserved quantity of the system [76]. Nonetheless, quenches from specific initial states may occasionally lead to thermalization in integrable models [77]. Figure 8a shows that, in the regime of small ε, the dynamics of the periodically kicked system H 2 forms a prethermal plateau; however the expectation values of observables in the plateau are marked by large fluctuations, whose origin can be traced back to the integrable character of H 2 . Figure 8b shows that the system prethermalizes approximately for small values of ε, as becomes evident from the eigenvalues of the reduced density matrix. Moreover, for the given initial energy density, we find that ETH is satisfied, provided PBCs are applied to the subsystem effective Hamiltonian H A eff [ Fig. 8c] (note that significant deviations appear when OBCs are applied, yet they cannot survive in the thermodynamic limit). Similar as compared to the nonintegrable drive generated by H 1 , also here our findings reach beyond those of ETH as the inverse temperature follows the theoretical prediction obtained from Eq. (9). However, in contrast to the nonintegrable drive H 1 , in the transition regime, ε ∼ 10 −3 ÷ 10 −2 , before reaching infinite temperature, the eigenvalues of the reduced density matrix show significant deviations from the expected exponential dependence, which survive with increasing the subsystem size [cf. App. B and Fig. 17 upper row]. A plausible explanation for this behavior is that the state of the system is not fully thermal in this ε-regime. This implies that, for the Floquet dynamics generated by H 2 , a substantial number of states in the Hilbert space are restrained from participating in the thermalization process even for L = 20 spins [cf. App. B, Fig. 15].
Interestingly, the thermal character of the state in the transition regime can be restored by adding a small noise δ to the driving period which breaks periodicity [cf. Sec. III] [cf. App. B, Fig. 17 lower panel]. Importantly, δ > 0 results in a thermal state well before the system reaches infinite temperature: the fit values for β shift systematically towards the ones obtained from Eq. (9) using the instantaneous energy densities [cf. App. B, Fig. 18(b-d)]. This corroborates our Hypothesis also for integrable Hamiltonians [cf. Fig. 18(a-d)]. Remarkably, finite noise restores ergodicity only in the unconstrained thermalization stage (III) of the dynamics between the prethermal regime and the featureless infinitetemperature state at long times; it hardly affects the prethermal properties of the dynamics, e.g. the expectation values of observables, and the time required to leave the prethermal plateau [ Fig. 8d]. In turn, this implies that the relevant effective Hamiltonian is not drastically affected by the addition of small noise.

B. Ising Model without Quantum Fluctuations
Last, let us discuss an Ising model without quantum fluctuations, described by the Hamiltonian where J = 1.0 and h z = 0.809. Quantum fluctuations in the driven system are introduced by the kicks V , cf. Eq. (3), so that the leading-order approximation to H F (ε) is non-integrable [93]. In this section, we consider a thermal initial state at β( = 0) = 0.8 [cf. Sec. III E].
The results below can be summarized in the following two points: (i) we provide numerical evidence that the drive generated by the Hamiltonian of Eq. (20) does not obey a Fermi Golden Rule scaling for the heating times. Instead the heating times cross over from a power-law scaling with an anomalous exponent α at large perturbation strength [cf. Sec. III D] to an exponential scaling for infinitesimal noise strengths. (ii) an IFE based on the Replica trick allows us to compute higher-order corrections to the effective Hamiltonian. We demonstrate that these higher order corrections are important to capture the physics in the vicinity of T * k . Fig. 9a displays the time evolution of the entanglement entropy for the dynamics generated by H 3 using the kicks from Eq. (5). Already from this figure it becomes evident that H 3 behaves quite different as compared to H 1 and H 2 : Instead of showing one stable prethermal level for different ε [as is the case for H 1 , cf. Fig. 1], here, different ε values result in different saturation levels at prethermal times (we checked that these plateaus are not a finite-size effects, see App. D, Fig. 19). Moreover, the required ε values to observe prethermal dynamics for H 3 , ε ∼ 10 −1 , are about two orders of magnitude larger compared to the previous two drives H 1 , H 2 , which correlates with the lack of quantum fluctuations in H 3 .
The varying saturation levels of the prethermal plateau complicate extracting the heating times. Yet, it is easy to recognize that the dynamics features two times scales: the first one captures the time needed to reach the prethermal plateau; it carries a clear dependence on ε as evident from Fig. 9a (black dots connected by solid  line). The second describes the time required to heat up to infinite temperature [ Fig. 9b]. For the periodic perturbation-free dynamics (δ = 0), the second time scales is intractable within the evolution time range of our simulations, which points to a much longer heating time. Interestingly, as opposed to former cases, finite periodicity-breaking noise leads to a significant reduction of the heating times so that a clear pattern becomes tractable.
To separate well the two timescales from each other, we apply Eq. (12) iteratively: first, we replace O prethermal by the initial expectation value and O β=0 by the prethermal value. The solution to the corresponding equation provides an estimate p for the time required to reach the prethermal plateau, i.e. the timescale of constrained thermalization, cf. the time required to cross the dasheddotted horizontal line in Fig. 9b. Independently, we attempt to solve Eq. (12) once again, yet this time we replace O prethermal with O( p ), which yields the time scale for unconstrained thermalization, cf. the time required to cross the dashed horizontal line in Fig. 9b. The results of this analysis are depicted in Fig. 9(c-d). The constrained thermalization time is clearly described by a power-law which survives a finite weak noise strength δ > 0, see Fig. 9c, inset. In contrast, the unconstrained heating time follows an exponential law over two decades for small δ/T , i.e., Γ −1 ∼ exp(−ξ ). We note that the exponential scaling clearly cannot persist as ε → 0, since this would imply finite heating times at ε = 0, where heating is inhibited by restored energy conservation. Increasing the noise strength δ leads to increasingly shorter heating times until eventually the prethermal plateau disappears and, therefore, the scaling crosses over to the power-law scaling of the unconstrained thermalization timescale.
These findings are intriguing, because they imply the existence of refined estimates for the scaling of the heat- ing rates with ε [so far, Floquet systems have mostly been treated on equal footing to derive a generic upper bound on the hating rate [44,45]]. One can even speculate about the existence of a wider class of Floquet models with suppressed heating behavior.

Replica Resummation and Thermalization
The significant change of the prethermal plateau level within the range of ε values we investigate [ Fig. 9a] implies that higher order corrections (in ε) to the effective Hamiltonian are of increased importance for understanding the dynamics of the system (as opposed to the models discussed in Secs. III and IV A). The drive H 3 was chosen to allow the analytical treatment of the leading-order correction to the average Hamiltonian using the replica expansion [69], cf. App. C. Note that the replica trick is needed at the commensurate points k > 0 for the kicked system, where higher-order nested commutator terms in the Baker-Campbell-Hausdorff series also contain terms to first order in ε, and hence one is required to re-sum an infinite subseries to correctly identify the first-order correction.
To facilitate the analytical computation, we switch back to a two-step protocol: Then, for k > 0, the replica expansion gives [cf. App. C]: where f χ (T ) = χ cot(χT /2). Note that f χ carries the only k-dependence up to first order via its argument T = 2(T * k +ε). Since it constitutes a periodic function, f χ (T + 2nπ/χ) = f χ (T ), in the regime π/γ ≈ nπ/χ, this implies a very weak dependence of the effective Hamiltonian on the value of k to first order in ε, consistent with our numerical observations.
Note that the cotangent function present in the firstorder terms, can lead to divergences, which likely persist also in higher-order terms [69]. Indeed, for the range of ε values we discuss, it turns out that J = 1.0 might occasionally lead to a divergence as we vary ε. Thus, subsequently we set J = 0.6.
More importantly, notice that H is nonintegrable, and hence the exact Floquet Hamiltonian is (most likely) also nonintegrable, which explains the origin of slow thermalization. Indeed, we observed that a number of pure initial states do not thermalize well within the accessible system sizes and evolution cycle, which is why we choose to prepare the system in a thermal initial state. Figure 10 shows the importance of the first-order corrections to properly capture the dynamics of the prethermal plateau for the Hamiltonian H 3 . The system is prepared in a thermal initial state w.r.t. H (0) eff , and two quantum quenches are performed: a quench to H (0+1) eff results in dynamics which saturates closer to the prethermal plateau obtained from exact time evolution, as compared to a quench to H (0) eff (which obeys no dynamics by construction).
The above check exhibits only the most obvious implication of the first-order correction. More subtly, higherorder corrections to the effective Hamiltonian also improve the description of the thermal state at all stages of the time evolution. To demonstrate this, we investigate the dynamics of two different thermal initial states: (i) an initial state, thermal w.r.t. H [cf. Fig. 11a]. In the vicinity of T * k , we analyze the thermal properties of the associated reduced density matrix along the way up to infinite temperature. Because we want to compare the two leading-order corrections, we now hold two ways of computing an associated inverse temperature: we can either use H , we consistently obtain a very good agreement of the fitted and the ETHpredicted inverse temperature for various values of ε and different times [ Fig. 11(b-c) lower panel]. By contrast, using H (0) eff produces sizable deviations in the fitted values vs. the ETH predictions, until the system is close enough to the infinite-temperature state [ Fig. 11(b-c) upper panel]. We note in passing that the deviations we observe between the fitted and the ETH-predicted inverse temperature can be observable-dependent [cf. App. D].
Considered altogether, the above analysis implies that H (0) eff is not the correct effective Hamiltonian the system thermalizes to, in both the prethermal plateau and the subsequent unconstrained thermalization regime. Instead, higher-order corrections are indeed required to properly capture the thermalizing Floquet dynamics generated by H 3 , and the Replica trick provides a useful expansion to compute them. Moreover, despite the different scaling of the heating times, even for nonintegrable effective Hamiltonians which obey ETH, we see that the IFE, supplemented with the instantaneous energy density, can provide a good description of the thermalizing dynamics of the evolved state throughout the entire evolution cycle all the way up to infinite temperature.
Interestingly, the Replica expansion also proves useful to better understand the sensitivity of the heating times to the perturbation strength in the case of noisy drives. As for H 1 and H 2 , the leading-order effective Hamiltonian H (0) eff only acquires an overall multiplicative factor Given the applicability of ETH, this translates into small energy density, and with it to temperature fluctuations of the order of δ/(2T ). However, the first order correction H eff is subject to more drastic changes. To first order in δ, we find Besides the fact that χδ/sin(χT ) can become large, finite δ/T also changes the relative weights of the different terms appearing in H eff , cf. Eq. (22). This results in applying a significantly different effective Hamiltonian for each period of the drive; thus, the resulting effect on the dynamics cannot be interpreted as small temperature fluctuations since also the eigenspectrum of H eff is subject to nonnegligible changes from one period to the next. Hence, due to lack of periodicity in the drive, we can no longer define prethermalization w.r.t. H

(0+1) eff
; as a result, we observe an increased sensitivity to the noise strength δ.

V. CONCLUSION
In summary, we investigated a class of step-driven Floquet systems, with the help of which it is possible to extend the high-frequency prethermal physics to low drive frequencies on the order of the Hamiltonian couplings. The systems in this class are distinguished by a frequency axis which, by construction, contains isolated stable points, where energy is conserved exactly. We demonstrated that these points come with prethermal windows, whose width as a function of the deviation ε from the commensurate point, varies between a power-law and an exponential, depending on the drive. Surprisingly, Fermi's Golden Rule is not universally applicable. Intriguing open questions are whether one can bridge these windows to enhance stability, and whether there exist models with a single continuous stable window all the way down to the low frequency regime.
Throughout the paper, we studied the thermalization behavior of three integrable and nonintegrable drives H j : the mixed-field Ising model and the transverse-field Ising model generate Floquet dynamics with quadratic in ε heating rates that follow Fermi's Golden Rule; in contrast, the Ising model without quantum fluctuations, exhibits a more refined heating behaviour that violates the Fermi Golden Rule scaling. Thus, intuition based on the equilibrium integrable-nonintegrable classification does not carry over to non-equilibrium systems in a straightforward manner; instead, we find an interesting correlation between whether Floquet heating is nonquadratically or quadratically suppressed in ε on one hand, and whether it is feasible to re-sum a subseries of the inverse-frequency expansion, on the other.
During the study, we introduced a new technique to facilitate Floquet systems to explore the entire underlying Hilbert space: we apply small random perturbation/noise in the duration of the Hamiltonian H j . We showed that this procedure minimizes finite-size effects and allows to look for a proper parametric dependence of the heating rates in the curves for the time evolution of physical quantities. Additionally, we observed that the k > 0 commensurate points also enhance the ergodic properties of Floquet systems, as compared to the infinite-frequency point k = 0, since the smaller drive frequencies at k > 0 provide the required spectrum folding for Floquet many-body resonances to occur. These technical advances allowed us to obtain clean scaling of the numerical data required for a proper study of thermalization in finite-size many-body systems, cf. App. A 1. As long as the prethermal physics is described by the leading-order in ε effective Hamiltonian, finite noise strengths do not reduce the duration of the associated prethermal plateau, as they effectively translate (in accord with ETH) into small temperature fluctuations. In contrast, when higher order corrections in ε become important for the thermalization dynamics, the prethermal physics becomes sensitive to the noise strengths.
We also provided numerical evidence in favor of the following conjecture: Consider a pure state subject to a periodic drive generated by a nonintegrable Floquet Hamiltonian, whose dynamics features a prethermal plateau. We observe that, upon leaving this plateau, a subsystem of the original system remains thermal w.r.t an effective Hamiltonian as defined by the IFE [to the order it can be computed/defined]. Although the expansion represents a divergent asymptotic series and does not capture the heating process itself, given the momentary value of the energy density, H eff is sufficient to construct a thermal ensemble which captures the dynamics of the system as it continues to heat up. The temperature of the subsystem increases gradually with time, and can be obtained from the energy density w.r.t. H eff . In contrast to this behavior, whenever the system heats up straight to infinite temperature (i.e. no prethermal plateau can form), we distinguish two scenarios: (i) in the high driving frequency limit (k = 0) the system is not in a thermal state until it reaches infinite temperature. (ii) for finite intermediate frequencies (k > 0), the system can still thermalize w.r.t. the ergodic H j but only if the corresponding intrinsic thermalization timescale for H j is smaller than T * k . It remains open whether the crossover as a function of ε between a (pre-)thermal state w.r.t. H eff , and a nonthermal state can become a sharp transition, and what conditions would be required for this to happen. Such a behavior could be detectable in the behavior of the thermodynamic entropy of the system which is maximal for a thermal state and smaller for any other state.
Our study also bears relevance to experiments. Recently, Floquet prethermal physics has been observed in a driven cold atomic system [86,94]. A straightforward application of our analysis is to use periodic drives to continuously fine-tune the temperature of cold-atomic systems in time. This could be used, e.g. to trigger temperature-driven phase transitions, such as the Mott insulator transition or the Kosterlitz-Thouelss transition. Related ideas can potentially prove useful to design a new temperature knob in quantum simulators which are well isolated from external reservoirs. Finally, we mention that it may soon be within the scope of present day quantum gas microscopes to reconstruct the density matrix of a subsystem via density matrix tomography, and verify or negate the conclusions and predictions of our work. The main quantity -the reduced diagonal density matrix -defines a natural measurement ensemble in experiments [81].
Emergent Phenomena in Quantum Systems initiative of the Gordon and Betty Moore Foundation, the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Quantum Algorithm Teams Program, and the Bulgarian National Science Fund within National Science Program VIHREN, contract number KP-06-DV-5. This research was supported in part by the International Centre for Theoretical Sciences (ICTS) during a visit for participating in the program -Thermalization, Many body localization and Hydrodynamics (Code: ICTS/hydrodynamics2019/11). We used Quspin for simulating the dynamics of the quantum systems [95,96]. The authors are pleased to acknowledge that the computational work reported on in this paper was performed on the Shared Computing Cluster which is administered by Boston University's Research Computing Services and on the Wuerzburg HPC cluster.  state. Here, finite size scaling effects are nearly completely absent. Typicality calculations are thus expected to be rather independent of the system size.

Frequency or (k-) Dependence
The analysis so far provides a strong indication for the qualitatively similar behavior of the k = 0 and k > 0 points. Here, we show a direct comparison between different commensurate points T * k . Unexpectedly, we find that there is no scaling of thermalization behaviour with increasing k (Fig. 14). This is surprising as the frequency drops with increasing k (keeping the duration of the kick, i. e. ε, constant) and the folding window of the Floquet spectrum becomes smaller so more interaction is expected. The lack of a k-dependence once more manifests that folding is only a necessary yet insufficient criterion and that the strength of the matrix elements is crucial to determine heating rates. In this Appendix we show the finite size scaling for the transverse-field Ising drive H 2 given in Eq. (20). In Fig. 15, we depict the dependence of the pure state dynamics as a function of the system size L for a noise-free and a noisy drive. In Fig. 16, we show the k-dependence of the dynamics. Qualitatively, we find an overall similar scaling behaviour as found in the nonintegrable case in App. A. Finally, in Figs. 17 and 18 we display the thermalization behaviour by investigating the spectrum of the reduced density matrix as a function of the subsystem size as well as the associated inverse temperature. We find that the thermal character of the state is lost in the transition regime (i.e. for ε ∼ 10 −2 ) [ Fig. 17 upper row]. Adding small noise δ to the driving protocol activates the inactive parts of the Hilbert space and in turn reintroduces thermalization even in the transition regime [ Fig. 17 lower panel, Fig. 18].   eventually obtain a closed form expression makes use of the replica expansion [69].
It is easy to verify that the following identity holds We aim to find an expression resummed in orders of the small expansion parameter ε. Using a Taylor expansion of U F in ε yields the generic expression (d) β( ) for the noise-free (lightblue) and noisy (darkblue) cases. The parameters are the same as in Fig. 1 of the main text.
Interchanging sum and limit yields an algebraic expression for H eff as a series in ε where we defined Next, we can insert the piece-wise constant step drive U F = exp(−iT /2H 3 ) exp(−iεV ). For the case of r = 1, simple algebraic manipulations lead to To obtain an expansion with a resummed subseries in ε, we need to evaluate the object U m 0 consists of two terms: (i) a single particle term, which essentially defines a single particle rotation around the z-axis and (ii) a many-body term also along the z-axis: Thus, we are allowed to apply each of the two terms separately, where the order does not matter. Let us start with the many-body rotation. Straightforward manipulations yield Next, we apply the single-particle z rotation and get To evaluate the sum over m in the replica resummation, a mode expansion of Eq. (C9) is required. Collecting terms with different exponents we find j (σ z j σ y j+1 + σ y j σ z j+1 ) It is possible, though tedious, to evaluate the leading higher-order terms numerically [69]. In this Appendix we provide more detailed data regarding the thermalization behaviour of analytically tractable drives, using dynamics generated by H 3 .
We begin by investigating the finite size scaling of the dynamics generated by H 3 . In the main text, we observed that intermediate values of ε can lead to a late-time saturation in the time-evolution curve of given observables away from their infinite-temperature value (see Sec. III D). Yet, a careful investigation of the dynamics driven by H 1 leads to the conclusion that this constitutes a finite size effect, which might be removed by the addition of noise in the drive protocol at finite system sizes [cf. App. A]. To rule out the possibility that the observed late-time plateau in the dynamics generated by H 3 originates from a similar effect, in Fig. 19 we compare the time-evolution curves of the entanglement entropy for different system sizes. Notice that the slight shift in the curves by approximately a constant with increasing L, is caused by the slight change in the initial energy density, corresponding to the fixed initial temperature, and the observed shift of the prethermal plateau matches the shift of the initial energy density. Thus, it is caused by a systematic mismatch in the initial energy densities and is not a finite-size effect -a fact corroborated also by the scale on the y-axis.
Although the dynamics of the system at late times is not expected to strongly depend on the initial state, one may want to reason that the observed agreement in the lower panels of Fig. 11(b-c), as well as the corresponding   Fig. 20]. In this setup, we do not find a good agreement of the fitted and ETHpredicted inverse temperatures, using either of the two effective Hamiltonians at small and intermediate driving times [Fig. 20(b-c)]; only at long driving times is the agreement restored when using H (0+1) eff , since the infinite-temperature state is a universal long-time attractor. Thus, using the zeroth-order term, H eff , an agreement of the numerically-fitted and the ETH-predicted temperatures is only reached close to infinite temperature.
Finally, we would also like to emphasize that the improvement brought by higher-order corrections to the effective Hamiltonian depends on the observable of interest, as suggested by Fig. 21. Indeed, Fig. 21 is equivalent to Fig. 11 of the main text, however, this time using H A,(0) eff (i.e. the effective Hamiltonian to leading order on the subsystem) as observable. The observed deviations are small already on the level of the zeroth order effective Hamiltonian. We emphasize that, Fig. 21 contains a feature that indirectly proves that H . This implies that the system is already initially in the correct thermal state w.r.t. the generator of dynamics so that no drive-induced initial quench dynamics occurs.