Stalled response near thermal equilibrium in periodically driven systems

The question of how systems respond to perturbations is ubiquitous in physics. Predicting this response for large classes of systems becomes particularly challenging if many degrees of freedom are involved and linear response theory cannot be applied. Here, we consider isolated many-body quantum systems which either start out far from equilibrium and then thermalize, or find themselves near thermal equilibrium from the outset. We show that time-periodic perturbations of moderate strength, in the sense that they do not heat up the system too quickly, give rise to the following phenomenon of stalled response: While the driving usually causes quite considerable reactions as long as the unperturbed system is far from equilibrium, the driving effects are strongly suppressed when the unperturbed system approaches thermal equilibrium. Likewise, for systems prepared near thermal equilibrium, the response to the driving is barely noticeable right from the beginning. Numerical results are complemented by a quantitatively accurate analytical description and by simple qualitative arguments.

In this context, the majority of previous studies focused on the long-time behavior and, in particular, on the properties of the so-called Floquet Hamiltonian.A key aspect of such an approach is that it can only capture the actual behavior of the periodically driven system stroboscopically in time, i.e., at integer multiples of the driving period, whereas the possibly still very rich behavior in between those discrete time points remains inaccessible.For instance, the stroboscopic dynamics may appear nearly stationary even though the full, continuous dynamics still exhibits oscillations with large amplitudes.
We adopt a complementary perspective and explore the continuously time-resolved response on short-tointermediate time scales.Intuitively, one might naturally expect that periodic forcing leads to a clearly noticeable change of the observable properties if its strength and period are of the same order as the main intrinsic energy and time scales of the undriven system.
In this work, we show that such a fairly pronounced response is indeed observed for isolated many-body systems that are far away from thermal equilibrium.Our main discovery, however, is that this intuitively expected response is strongly suppressed near thermal equilibrium, at least as long as heating effects of the driving remain negligible.We dub this phenomenon "stalled response" in view of its two principal manifestations: For a system that is prepared far away from equilibrium, the observ-able response dies out as soon as the corresponding undriven reference system approaches thermal equilibrium.Similarly, when the system already starts out in thermal equilibrium, the driving is barely noticeable right from the beginning.In both cases, it is only at much later times that the driving effects may reappear in the form of very slow heating.Besides numerical evidence from several examples, we support our general prediction of stalled response near thermal equilibrium with simple heuristic arguments and with an analytical theory for large classes of many-body systems.Remarkably, we can also identify the main qualitative signatures of such a stalled response behavior in data from a very recent NMR experiment [8].
driving is monitored by the deviations of ⟨A⟩ ρ(t) from ⟨A⟩ ρ0(t) .

Phenomenology
To illustrate the announced phenomenon of stalled response, we first present a numerical example in Fig. 1.Its specific choice is mainly motivated by the fact that it will admit a direct comparison with our analytical theory (presented below) without any free fit parameter.Further examples will be provided later.
As sketched in Fig. 1, we consider an L × L spin-1 2 lattice with L = 5 and open boundary conditions, where nearest neighbors are coupled by Heisenberg terms in the unperturbed system (solid links in the sketch), The vector σ i,j = (σ x i,j , σ y i,j , σ z i,j ) collects the Pauli matrices acting on site (i, j).The perturbation additionally introduces spin-flip terms in the z direction between next-nearest neighbors (dashed links in the sketch), Since the magnetization S z := i,j σ z i,j commutes with both H 0 and V , we focus on one of the two largest subsectors, namely the one with eigenvalue −1 for S z .To prepare the system out of equilibrium, we fix the spins at sites (2, 2) and (3,3) in the "up" state (red in the sketch in Fig. 1) and orient all other spins randomly.To obtain a well-defined energy, we additionally emulate a macroscopic energy measurement by acting with a Gaussian filter [15][16][17] of a target mean energy E = −12 and standard deviation ∆E = 4 on the so-defined state.Formally, the initial condition can thus be expressed as where |ϕ⟩ is a Haar-random state in the S z = −1 sector.The projector Q := π + 2,2 π + 3,3 with π ± i,j := (1 ± σ z i,j )/2 enforces σ z 2,2 = σ z 3,3 = 1, and this deflection is only weakly reduced by the subsequent Gaussian energy filter (cf.Fig. 1).From a different viewpoint, the situation may also be seen as a small non-equilibrium system in contact with a large thermal bath (red and black vertices, respectively, in the sketch).
Accordingly, an obvious choice for the considered observable is the correlation between the initially disequilibrated sites, A = σ z 2,2 σ z 3,3 .Incidentally, the ground-state energy of H 0 from (3) is approximately −60, whereas the infinite-temperature state has an energy of approximately −1.Hence, our choice of the target energy E = −12 should be reasonably generic and corresponds, as detailed in Supplementary Note 2.2, to an inverse temperature β ≈ 0.08.Further examples for different target energies/temperatures can also be found in Supplementary Note 2.2.
In Fig. 1 we present numerical results, obtained by Suzuki-Trotter propagation, for the unperturbed system H 0 and for a sinusoidally driven system (1) with yielding the solid black and blue lines, respectively.
The key observation is that the driven (blue) and undriven (black) expectation values in Fig. 1 differ quite notably during the initial relaxation of the unperturbed system, but they become (nearly) indistinguishable upon approaching their (almost) steady long-time values.Moreover, both long-time values agree very well with the thermal expectation value A th ≃ −0.026, obtained numerically by evaluating A = σ z 2,2 σ z 3,3 in the microcanonical ensemble of the unperturbed system.In other words, the perturbations by the periodic driving get stalled upon thermalization of the undriven system.
To further highlight this phenomenon, let us also consider the analogous equilibrium initial conditions with Q = 1 in (5).Hence, the initial state populates the same energy window as in the nonequilibrium setting, but the observable expectation values now (approximately) assume the pertinent thermal equilibrium values [15][16][17].The solid green and red lines in Fig. 1 illustrate the so-obtained numerical results for the unperturbed and the driven system.In particular, the initial expectation value is now very close to the thermal equilibrium value A th ≃ −0.026.Moreover, the effects of the driving are indeed barely noticeable, and are even expected to become still smaller for larger system sizes, as detailed in Supplementary Note 2.3.
The bottom line of all these numerical findings is that the same system exhibits a quite significant response to the periodic driving away from thermal equilibrium, but hardly shows any reaction to the same driving as the unperturbed system approaches thermal equilibrium, or if it already started out near thermal equilibrium (stalled response).
Note that the driving amplitudes in Fig. 1 are far outside the linear response regime, as can be inferred, e.g., by comparing the blue curves of Figs.1c and f (see also Supplementary Note 2.1).We also remark that for noncommuting perturbations and observables (as in Fig. 1), linear response theory generically excludes that there is no response at all.The main challenge is to understand why the non-linear response remains so weak at thermal equilibrium.
Likewise, the observable response becomes uninterestingly weak for extremely small or large driving periods T , regardless of the initial conditions and their proximity to thermal equilibrium.Hence, our focus here is on the natural regime of moderate T values that are similar to, or slightly below the relaxation time of the unperturbed system, where the stalling effect is most pronounced and interesting.The interplay of the various time scales is further elaborated in Supplementary Note 1.1.

Theory
Our next goal is to establish an analytical theory for reasonably general classes of many-body quantum systems which explains these numerical findings.We start by collecting the basic ingredients and assumptions, then present the main result, and finally sketch the derivation.
First, we focus on initial states ρ(0) with a welldefined macroscopic energy.Denoting by E µ and |µ⟩ the eigenvalues and -vectors of the unperturbed Hamiltonian H 0 , this means that non-negligible level populations ⟨µ|ρ(0)|µ⟩ only occur for energies E µ within a sufficiently small energy interval ∆, such that the density of states can be approximated by a constant D 0 throughout ∆.
Second, within this energy interval ∆, the matrix elements V µν := ⟨µ|V |ν⟩ of the perturbation operator V are assumed to exhibit a well-defined perturbation profile ṽ where [ In passing, we note that at sufficiently high temperatures, v(t) can be approximated by the two-point correlation function ⟨V (t)V ⟩ ρmc /2, where V (t) := e iH0t V e −iH0t and ρ mc is the microcanonical ensemble corresponding to the energy interval ∆; see Supplementary Note 3 for details.Third, the time-dependent perturbations f (t)V in (1) should not become overly strong compared to H 0 , so that establishing a connection between the unperturbed and driven systems remains sensible and the above mentioned heating effects stay reasonably weak.
In terms of the above introduced quantities, our main analytical result is the prediction where A th = tr(ρ mc A) is the thermal expectation value introduced below Eq. (6).The driving effects are encoded in the response function γ τ (t), evaluated at τ = t in (9), which is obtained as the solution of the parametrically τ -dependent family of integro-differential equations with initial condition γ τ (0) = 1 and coefficients where . We emphasize that the theory and Eq. ( 10) in particular are nonlinear, which -in light of the numerically observed response characteristics (see Fig. 1) -is essential to faithfully reproduce the observed behavior.
To derive these results, we combined and advanced three major theoretical methodologies: (i) a Magnus expansion [27] for the propagator U(t) (see below Eq. ( 2)); (ii) a mapping of the time-dependent problem (1) to a parametrically τ -dependent family of time-independent auxiliary systems; (iii) a typicality (or random matrix) framework [28][29][30] to determine the generic behavior (9) for the vast majority of all systems sharing the same H 0 , ṽ(E), and f (t).Details of the derivation are collected in the Methods.
Of the adopted techniques, the Magnus expansion in particular implies that such an approach only covers the transient dynamics up to a certain maximal time, which increases as the driving period T becomes smaller.Since this maximal time has been related to the onset of heating [18,22,31], the result (9) does not capture such heating effects anymore.Yet it may well remain valid over a quite extended time interval since heating is suppressed exponentially for small T [7,8,[23][24][25][26], see also Supplementary Note 1 for a more detailed discussion of the relevant time scales and of the response function γ τ (t).
Due to the employed typicality framework, in turn, the prediction (9) may not reproduce the dynamics accurately in certain setups with strong correlations between the observable A and the perturbation V .
A more in-depth discussion of the expected regime of applicability is provided in the Methods.

Interpretation and further examples
For a quantitative comparison of our theoretical prediction (9) to specific examples, some approximate knowledge of the perturbation profile (7) is clearly indispensable.Qualitatively, however, the theory quite remarkably allows us to make some largely general predictions without any such specific knowledge.
The first and foremost of these predictions is based on the general upper bound |γ t (t)| ≤ 1, whose detailed analytical derivation is provided in Supplementary Note 7 (see also Supplementary Note 1.2).It then immediately follows from (9) that the driving effects are strongly suppressed whenever ⟨A⟩ ρ0(t) ≃ A th , i.e., whenever the unperturbed system is close to thermal equilibrium.The latter in turn is true for all times t if the unperturbed system is at thermal equilibrium from the outset, and for all sufficiently late times t if the unperturbed system starts out far from equilibrium and is know to thermalize in the long run.Altogether, our stalled response phenomenon is thus analytically predicted to occur under very general circumstances.
Next we turn to a more detailed quantitative comparison of the theoretical prediction (9) with concrete nu-merical examples.For the setup considered in Fig. 1, exact diagonalization of a smaller system with L = 4 [32] suggests that the perturbation profile ṽ(E) from (7) can be approximated very well by an exponential decay ṽ(0) e −|E|/∆v .Utilizing Ref. [30], one moreover finds for the L = 5 system in the relevant energy window the numerical estimates ṽ(0)D 0 ≃ 3.6 and ∆ v ≃ 3.0, yielding v(t) via (8).All quantities entering the theoretical prediction ( 9)- (10) are thus either numerically available [⟨A⟩ ρ0(t) , A th ] or otherwise known [v(t), a τ , b τ ], i.e., there remains no free fit parameter.
As can be inferred from the solid blue and dashed purple lines in Fig. 1, the theory indeed describes the nontrivial details of the driven dynamics remarkably well.Notably, it reproduces the pronounced drop compared to the unperturbed curve around t = T /2 and the quite surprising comeback around t = T .Moreover, it indeed also explains the stalled response behavior in Fig. 1 very well, for initial conditions both close to and far from thermal equilibrium.
Within the framework of Floquet theory, a related, but distinct effect is well-known under the name "Floquet prethermalization" [7,8,20,23,24,26,33,34]: The dynamics described by the Floquet Hamiltonian approaches a prethermal plateau value before heating becomes significant and pushes the system towards infinite temperature.However, the dynamics encoded in the Floquet Hamiltonian only agrees with the actual dynamics of the driven system stroboscopically, i.e., only at integer multiples of the driving period.A prethermal plateau of the Floquet-Hamiltonian dynamics therefore still leaves room for strong oscillations of the actual dynamics between the stroboscopic time points where both agree.Accordingly, the salient new insight provided by our present results is that no such strong oscillations are observed if the unperturbed system relaxes to or starts out from a thermal equilibrium state.In other words, our stalled response effect amounts to a highly nontrivial extension of the established Floquet prethermalization phenomenon since it means that the plateau value is assumed not only stroboscopically, but even continuously in t.An extended discussion of the relation between our approach and Floquet theory can be found in Supplementary Note 4.
As a second example, we consider a nonintegrable variant of the transverse-field Ising model in Fig. 2, see the figure caption for details.We particularly emphasize that, for variety and in contrast to Fig. 1, this setup consists of a one-dimensional system and globally out-ofequilibrium initial conditions.
Qualitatively, the numerical results in Fig. 2 once again confirm the main message of our paper, namely the occurrence of stalled response: Initially, the dynamics shows a pronounced response when starting away from equilibrium (solid black vs. blue lines).Stalling of that response appears as the unperturbed system approaches thermal equilibrium, meaning that the oscillations caused by the driving become smaller and smaller.This is highlighted in the insets, in particular.(A special feature and L = 24.The driving operator is a longitudinal magnetic field, V := −J L j=1 σ x j .Time-dependent expectation values ⟨A⟩ ρ(t) of (a) the single-site magnetization A = σ z 1 , (b) the nearest-neighbor correlation A = σ z 1 σ z 2 , and (c) the next-nearest-neighbor correlation A = σ z 1 σ z 3 are shown for the periodically driven system (1), (6), with driving amplitude f0 = 4 and period T = 0.5.Solid black and blue lines: numerical results for nonequilibrium initial conditions (5)  (The corresponding inverse temperature, ground-state energy, and infinite-temperature energy are now approximately 0.2, −18.5, and 0, respectively, see also above Eq.( 6).) Solid green and red lines: same but for equilibrium initial conditions (5), i.e., with a Haar-random state |ϕ⟩.Dashed lines: corresponding theoretical predictions (9), adopting the numerically obtained unperturbed behavior ⟨A⟩ ρ 0 (t) , squared response function |γt(t)| 2 (by numerical integration of ( 10)), and thermal equilibrium values A th ≃ 0.066, 0, 0 in (a), (b), (c), respectively.Insets: Same numerical data, but with rescaled x and y axes to display the long-time behavior.FIG. 3. Imperfect stalling upon breaking a conservation law.Time-dependent expectation values ⟨A⟩ ρ(t) of the single-site magnetization A = σ z 1,1 are shown for a periodically driven 2 × 14 spin double-chain (see inset) with Hamiltonian (1), ( 6), ( 12), ( 13), and driving period T = 0.25.Solid black and blue lines: numerical results for non-equilibrium initial conditions (5) with , for driving amplitudes f0 = 0 (unperturbed, black) and f0 = 3.2 (driven, blue).Solid green and red lines: same but for equilibrium initial conditions (5) with Q = 1.Dashed lines: corresponding theoretical predictions (9), obtained as in Fig. 1 but with A th = 0.
of this example is that already the unperturbed system (black lines) exhibits a relatively complex and longlasting relaxation process.)Likewise, the effects of the driving are barely visible on the scale of the plot when starting directly from a thermal equilibrium state (solid green vs. red lines).
For a quantitative comparison of the numerical results with the theoretical prediction (9), we assume, as in the previous example, an approximately exponential perturbation profile ṽ(E) = ṽ(0)e −|E|/∆v [cf.Eq. ( 7)], and use again the theory from Ref. [30] to estimate ṽ(0)D 0 ≃ 0.46 and ∆ v ≃ 0.6.The resulting theoretical curves in Fig. 2 (dashed lines) describe the numerics reasonably well in the initial regime.In accordance with the discussion below Eq. ( 11), for larger times the theory is no longer quantitatively very accurate (but still correctly predicts the occurrence of stalling per se).For this reason, no dashed lines are shown in the insets.
Yet another interesting general prediction of the theory (9) (see also beginning of this section) is that noticeable effects of the driving (as encoded in |γ t (t)| 2 ) may actually persist even beyond the relaxation time scale of the unperturbed system if its long-time expectation value Ā0 := ⟨A⟩ ρ0(t) (infinite time average) differs from the thermal value A th .This can happen, for example, if the perturbation V breaks a conservation law of H 0 .
To verify this prediction, we consider a third example in Fig. 3.Here the unperturbed system consists of two isolated spin chains of L = 14 sites with periodic bound-ary conditions and Hamiltonian while the perturbation in (1) connects the chains sitewise, see also the sketch in the inset.The initial state is again of the form ( 5) with E = −14 and ∆E = 4, restricted to the sector with vanishing S z := j (σ z 1,j + σ z 2,j ).(The corresponding inverse temperature, ground state energy, and infinite-temperature energy are now approximately 0.12, −50, and −1, respectively see also above Eq.( 6).)However, for the nonequilibrium setup we now fix two spins in the "up" state for the first chain and two in the "down" state for the second chain (red and blue, respectively, in the sketch), i.e., Q : . Since the two chains (i = 1, 2) do not interact in the unperturbed system, their magnetizations S z i := j σ z i,j are conserved individually, and thus maintain their initial expectation values 2 and −2, respectively, under evolution with H 0 .In the driven system, by contrast, only the total S z := S z 1 + S z 2 is conserved.Choosing the single-site magnetization A = σ z 1,1 as our observable, we thus find by symmetry that Ā0 = 2/L is the long-time expectation value of the unperturbed dynamics, whereas the thermal value of the joint system is A th = 0.
The numerics in Fig. 3 (solid blue line) visualizes the aforementioned imperfect stalling upon breaking a conservation law: The suppression of the response is the stronger the closer the unperturbed system is to thermal equilibrium.Crucially, however, the driving effects still remain visible even when the unperturbed dynamics has essentially reached its nonthermal long-time value Ā0 .Altogether, this confirms the prediction of ( 9) that proximity to thermal equilibrium is indeed the decisive condition for stalled response and not, for example, relaxation of the unperturbed system.Furthermore, this example highlights once again that stalled response and Floquet prethermalization are distinct effects: The present system exhibits Floquet prethermalization, meaning that the stroboscopic dynamics approaches a stationary plateau, but no stalled response since ⟨A⟩ ρ(t) continues to oscillate.
For a quantitative comparison with the theory (9), we again adopt the same ansatz as before and estimate ṽ(0)D 0 = 0.98 and ∆ v = 4.2 via [30].The so-obtained prediction (9) (dashed purple) agrees rather well with the numerics for t ≲ 1.At later times, the quantitative deviations between the prediction and the numerics increase.As suggested below (11) and discussed in more detail in the Methods, we can attribute these deviations to the adopted Magnus expansion and its truncation at second order.Yet the above mentioned general qualitative prediction of our theory remains valid nonetheless.

Basic physical mechanisms
Intuitively, the basic physics behind all our above mentioned numerical and analytical findings can also be understood by means of the following simple arguments: As long as heating is insignificant, we may focus on the dynamics within the initially populated energy interval ∆ (see above (7)).Denoting by P the projector onto the eigenstates |µ⟩ with E µ ∈ ∆, the Hamiltonian H(t) from (1) can thus be reasonably well approximated by its projection/restriction H(t) := P H(t)P to ∆.Since the microcanonical ensemble ρ mc := P/ tr{P } commutes with H(t), it is a stationary state with respect to H(t).Within the present approximation, a system in thermal equilibrium is thus completely unaffected by the periodic driving, and analogously the effects remain weak if the system is in a state close to thermal equilibrium.(Incidentally, the relaxation of a non-equilibrium initial state under H(t) can be heuristically understood by similar arguments as in Ref. [19].)On the other hand, subleading effects like small remnant oscillations and slow heating cannot be understood within this simplified picture.Rather, these effects must be attributed to the neglected corrections H(t) − H(t) and, as a consequence, are intimately connected with each other.
A complementary, and even more simplistic argument is based on the well-established fact [35][36][37] that the vast majority of all pure states with energies in ∆ behave akin to ρ mc for sufficiently large many-body systems.This socalled typicality property suggests that once the system has reached (or starts out from) such a state, it remains within this vast majority in the absence as well as in the presence of the periodic driving.
Essentially, our stalled response effect thus seems to be the result of a subtle interplay between the system's many-body character and intriguing peculiarities of thermal equilibrium states.The above intuitive arguments moreover suggest that the indispensable prerequisites for stalled response per se may be substantially weaker than those of our analytical theory (see also Supplementary Note 2.4).

DISCUSSION
Our core message is that the same many-body system may either exhibit a quite significant response when perturbed by a periodic driving, or may not show any notable reaction to the same driving, depending on whether the unperturbed reference system finds itself far from or close to thermal equilibrium.We demonstrated this stalled response effect by numerical examples, and further substantiated it by sophisticated analytical methods and by simple physical arguments.
Complementary to those long-term features for discrete time points, our present focus is on how a manybody system approaches such prethermal regimes continuously in time.Overall, we thus arrive at the following general picture for periodically driven systems with moderate driving periods and amplitudes: Given a thermalizing unperturbed system that is prepared sufficiently far from equilibrium, the periodic perturbations generically lead to quite notable response effects on shortto-intermediate time scales.Subsequently, the expectation values approach a (nearly) time-independent behavior.On even much larger time scales, the system finally heats up to infinite temperature, manifesting itself in a slow drift of the expectation values towards their genuine infinite-time limits.
In principle, our predictions can be readily tested with presently available techniques in, for example, cold-atom [1][2][3][4][5][6] or polarization-echo [6][7][8] experiments.In practice, previous experimental (as well as theoretical) investigations mostly focused on the long-time behavior and stroboscopic dynamics.A notable exception is the NMR experiment from Ref. [8]: In Figs.3(a) and 5(a,b) therein, the NMR signal of the initially out-of-equilibrium system undergoes vigorous oscillations at first (called "transient approach" in [8]).Then, their amplitude gradually decreases as the running mean approaches a quasistationary value (called "prethermal plateau" in [8]).Even later, the only noticeable effect of the driving is a slow drift as the system heats up (called "unconstrained thermalization" in [8]).Unfortunately, the available experimental details are not sufficient to compare the measurements quantitatively with our analytical theory (9).Nevertheless, the observed NMR signal clearly shows the general qualitative features of stalled response as predicted by Eq. ( 9).

METHODS
We first lay out the three main steps in the derivation of ( 9)- (10), and subsequently address the expected validity regime of the employed approximations.

Magnus expansion
The time evolution of the driven quantum system with Hamiltonian H(t) from ( 1) is encoded in the propagator U(t) introduced below Eq. ( 1), which satisfies the Schrödinger-type equation d dt U(t) = −iH(t)U(t).Whereas this equation is formally solved by an (operator-valued) exponential for time-independent Hamiltonians, no such simple solution is available for the driven case.To make progress while keeping the setting as general as possible, we adopt a Magnus expansion [27] of the propagator, writing where the individual terms Ω k (t) in the exponent consist of integrals over k − 1 nested commutators of H(t) at different time points.The virtue of the Magnus series compared to other expansion schemes (e.g., a Dyson series) is that U(t) remains unitary when truncating ( 14) at a finite order.For Hamiltonians of the specific form (1), the first two terms of the general Magnus expansion (see, e.g., Ref. [27]) can be readily rewritten as where [V, H 0 ] := V H 0 − H 0 V (commutator), and F 1,2 (t) are defined below Eq. ( 11).

Mapping to auxiliary systems
Adopting the Magnus expansion (14), the propagator U(t) = e Ω(t) assumes an exponential form similar to the case of time-independent Hamiltonians.However, the time dependence of the exponent is generally still complicated.To proceed, we introduce a one-parameter family of time-independent auxiliary Hamiltonians where τ > 0 is treated as an arbitrary but fixed parameter.Starting from the same initial state ρ(0) as in the actual system of interest, any of these Hamiltonians H (τ ) generates a time evolution with the state at time t given by Since ρ(t) = U(t)ρ(0)U(t) † , the combination of Eqs. ( 14), (16), and (17) implies that the state ρ(t) of the driven system of interest coincides with the time-evolved state of the auxiliary system H (t) at time t, i.e., Hence finding the dynamics of the original driven system is equivalent to determining the behavior of all the auxiliary systems with time-independent Hamiltonians H (τ ) up to time t = τ , respectively.
Restricting ourselves to the second order of the Magnus expansion, we adopt Eqs.(15) in (16) to approximate the auxiliary Hamiltonians as with thereby splitting off the τ -independent reference Hamiltonian H 0 .

Typicality framework
It is empirically well established that the macroscopically observable behavior of systems with many degrees of freedom can be described by a few effective characteristics despite the vastly complicated dynamics of their individual microscopic constituents.Detecting and separating the macroscopically relevant properties of a manybody system from the intractable microscopic details can arguably be considered as the paradigm of statistical mechanics.The final component of our toolbox to describe the driven many-body dynamics aims at adopting such an approach to the observable expectation values ⟨A⟩ ρ(t) .
To this end, we start with the Hamiltonian H(t) = H 0 + f (t)V from (1) and temporarily consider an entire class (or a so-called ensemble) of similar driving operators V .Ideally, we would like to establish that all members of such an ensemble exhibit the same observable dynamics.In practice, what is analytically feasible is a slightly weaker variant of such a statement.Namely, we demonstrate that nearly all members V of the ensemble show in very good approximation the same typical behavior, and that the fraction of exceptional members, leading to noticeable deviations from the typical behavior, is exponentially small in the system's degrees of freedom.
In essence, the defining characteristic of the considered ensembles is the perturbation profile ṽ(E) from (7).Introducing the symbol E[ • • • ] to denote the average over the V ensemble, the matrix elements V µν are treated as independent (apart from the Hermiticity constraint, Hence the property (7) of the true perturbation is built into the ensemble in an ergodic sense, i.e., upon replacing local averages [ • • • ] E (see below Eq. ( 7)) by ensemble averages Due to a generalized central limit theorem (cf.Supplementary Note 6), these first two moments are essentially the only relevant characteristics of the V ensemble, i.e., the precise distribution of the V µν can still take rather general forms.A detailed definition of the admitted ensembles is provided in Supplementary Note 5.
For time-independent Hamiltonians of the form H = H 0 + λV with a constant (time-independent) perturbation, it was demonstrated in Refs.[29,30] that those ensembles can indeed be employed to predict the observed dynamics in a large variety of settings.In the following, we will extend the underlying approach to the auxiliary Hamiltonians H (τ ) of the form (19).The distribution of the V µν thus induces a distribution of the matrix elements V (τ ) µν := ⟨µ|V (τ ) |ν⟩ of V (τ ) from (20).In particular, we obtain E[V (τ ) µν ] = 0 and, together with the definitions ( 7), (11), and (20), As a first step of our typicality argument, we then calculate the ensemble average E[⟨A⟩ ρ(t,τ ) ] of the timeevolved expectation values.Deferring the details to Supplementary Note 6, we eventually obtain the relation Here a Fourier transformation relates the response function (see above (10)) via to the function G(z, τ ), which solves (24) and encodes the ensemble-averaged resolvent of In Supplementary Note 7, we furthermore show that Eqs. ( 23) and (24) imply the relation (10) for γ τ (t).
As a next step, we turn to the deviations ξ(t, τ ) := ⟨A⟩ ρ(t,τ ) − E[⟨A⟩ ρ(t,τ ) ] between the driven dynamics induced by one particular perturbation operator V and the average behavior.More explicitly, we inspect the probability P(|ξ(t, τ )| ≥ x) that a randomly selected perturbation V generates deviations ξ(t, τ ) that are larger than some threshold x.As explained in more detail in Supplementary Note 8, we can find a constant δ = 10 −O(N dof ) (decreasing exponentially with the system's degrees of freedom N dof ) such that where ∆ A is the measurement range of A (difference between its largest and smallest eigenvalues).In other words, observing deviations which exceed some exponentially small threshold value becomes exponentially unlikely as the system size increases, a phenomenon that is also sometimes called "concentration of measure" or "ergodicity" in the literature.Consequently, becomes an excellent approximation for the vast majority of perturbations V in sufficiently large systems.Combining Eqs. ( 18), (22), and (26), we thus finally recover our main result (9).

Limits of applicability
The class of systems whose Hamiltonian can be written in the form (1) is extremely general.However, the methods described above contain three major assumptions or idealizations that restrict the types of admissible setups to some extent.
The first issue arises when adopting the Magnus expansion (14) for the propagator U(t).The question of its convergence is generally a subtle issue and rigorously guaranteed in full generality only up to times t such that the operator norm ∥H(s)∥ satisfies t 0 ds ∥H(s)∥ < π, but can extend to considerably longer times in practice nonetheless [27].Due to the extensive growth of H(t) with the degrees of freedom, guaranteed convergence is thus very limited for typical many-body systems, but the expansion can still remain valuable as an asymptotic series for short-to-intermediate times [23,33].For periodically driven systems in particular, the (Floquet-)Magnus series amounts to a high-frequency expansion and thus works best for small driving periods T [27,43].More generally, the smaller the characteristic time scale of the driving protocol f (t) is, the larger is the time up to which the expansion offers a satisfactory approximation at any fixed order.
Physically, the breakdown of the Magnus expansion has been related to the onset of heating [18,22,31].Generically, many-body systems subject to perpetual driving are expected to absorb energy indefinitely and heat up to a state of infinite temperature [18][19][20][21][22], unless there are mechanisms preventing thermalization such as an extensive number of conserved quantities [38,42] or many-body localization [21,39,44].Nevertheless, under physically reasonable assumptions about the system, such as locality of interactions, it has been shown that the heating rate is exponentially small in the driving frequency [23][24][25][26].For sufficiently fast driving, therefore, energy absorption is essentially suppressed for a long time and the Magnus expansion can provide a good description of the dynamics.A more quantitative discussion of the interdependence of the relevant time scales is provided in Supplementary Note 1.1.
In summary, the Magnus expansion is expected to work as long as the state ρ(t) stays roughly within the initially occupied microcanonical energy window ∆ of the unperturbed reference Hamiltonian introduced above Eq.( 7).Consequently, the stalled-response effect and the applicability of the prediction (9) are generally expected to persist for longer times at larger initial temperatures because the relative influence of heating is smaller in this case.Furthermore, higher temperatures come with a higher density of states, such that finite-size effects are smaller, too.The temperature dependence is discussed in more detail in Supplementary Note 2.2.
A second limitation is our truncation of the Magnus expansion at second order.In general, this will further restrict applicability towards shorter times and/or faster driving, but still leaves room for a broad and in-teresting parameter regime as demonstrated examplarily in Figs.1-3.In principle, including higher-order terms may be possible, even though it leads to severe technical complications in the typicality calculation outlined above (see also Supplementary Note 6), and is thus beyond the scope of our present work.Besides the response function γ t (t), higher-order corrections are also expected to affect the long-time value (A th in Eq. ( 9)): It is well known from Floquet theory that this plateau value of Floquet prethermalization is controlled by the Floquet Hamiltonian [23,24,26,33,34].The latter agrees with H 0 to lowest order, but can yield different long-time behavior in general, even though the corrections are generically expected to be small [23].
A third potentially limiting factor for the applicability of our present approach is the typicality framework, within which we introduce ensembles of matrix representations V µν of the driving operator V in the eigenbasis of the reference Hamiltonian H 0 .Our main result states that the observable dynamics of nearly all members V of such an ensemble is described by Eqs. ( 9) and ( 10) (up to the limitiations discussed earlier).The final point to establish is that the true (non-random) driving operator V of actual interest is one of those typical members of the ensemble, which evidently requires a faithful modeling of the system's most essential properties with regard to the observable dynamics.
The classes of perturbation ensembles considered here are a compromise between what is physically desirable and mathematically feasible.From a physical point of view, we would like to emulate the matrix structure of realistic models as closely as possible.We therefore explicitly incorporate the possibility for sparse (most V µν are strictly zero) and banded (the typical magnitude |V µν | decays with the energy separation |E µ − E ν | of the coupled levels) perturbation matrices.These features indeed commonly arise as a consequence of the local and few-body character of interactions in realistic systems as supported by semiclassical arguments [45,46], analytical studies of lattice systems [47,48], and a large number of numerical examples (e.g.Refs.[49][50][51]).Similar assumptions are also well-established in random matrix theory and in the context of the eigenstate thermalization hypothesis [28,[52][53][54].On the other hand, the geometry of the underlying model and the structure of interactions (for instance their locality) are not explicitly taken into account.Therefore, the existence of macroscopic transport currents as a consequence of macroscopic spatial inhomogeneities can likely invalidate the prediction ( 9)- (10), at least for observables A which are sensitive to such initial spatial imbalances and their equalization in the course of time.This is ultimately related to our idealization of statistically independent matrix elements V µν for µ ≤ ν.In any realistic system, some of the matrix elements will certainly mutually depend on each other.However, it is generally hard to identify (let alone quantify) potential correlations in any given system, so independence may also be understood as unbiasedness in the absence of more detailed information.Moreover, mild correlations will often not have a noticeable impact on the properties relevant for the observable dynamics [55].
A specific case where correlations can become relevant, though, are observables A that are strongly correlated with the perturbation V , most notably if A = V .Since we keep the observable fixed when calculating ensemble averages, most members of the V ensemble will obviously violate such a special relationship.Unfortunately, it is not straightforwardly possible to adapt the method such that the case A = V can be described as well because including A = V in the ensemble averages would also affect the unperturbed reference dynamics ⟨A⟩ ρ0(t) .Numerical explorations and further discussions of this case are provided in Supplementary Note 2.4.Notably, the qualitative predictions of the theory (9) and, in particular, the occurrence of stalled response can still be seen for the observable A = V .
For the rest, we emphasize that it is not necessary for all members V of a certain ensemble to be physically realistic.The decisive question is whether their majority embody the key mechanism underlying the observable dynamics in the same way as the true system of interest.To give an example from textbook statistical mechanics, a large part of states contained in the canonical ensemble (as a mixed density operator) will be unphysical, and yet its suitability to characterize macroscopically observable properties of closed systems in thermal equilibrium is unquestioned provided that the temperature as the pertinent macroscopic parameter is chosen appropriately.
More generally, the probabilistic nature of the result implies that any given system can show deviations even if all prerequisites are formally fulfilled, but the probability for such deviations is exponentially suppressed in the system's degrees of freedom, cf.Eq. ( 25).For generic many-body systems, we therefore cannot but conclude that Eqs. ( 9)-( 10) are expected to hold unless there are specific reasons to the contrary.The explicit example systems from Figs. 1-3 only corroborate this observation, noticeably even though the number of degrees of freedom is still far from being truly macroscopic in those systems.

Response and heating time scales
Generically, a many-body system is expected to absorb energy and thus heat up under periodic driving [18-22, 31, 40, 41, 56].As explained in the main paper, a key prerequisite to observe stalled response is that these heating effects are sufficiently suppressed such that the dynamics remains practically confined to the initially occupied (microcanonical) energy shell for the relevant response and stalling time scales.
As a first general guess, the characteristic time scale t resp of the observable response for a system away from equilibrium will often be on the order of the characteristic driving time scale T and thus decrease as t resp = O(T ) for small T .Within our theory (cf.(9) in the main paper), the response is encoded in |γ t (t)| 2 .Hence, t resp is the typical time scale of the solutions of (10) in the main paper, an example of which is shown in Fig. S1.For concreteness, we take t resp to be the time at which the first minimum of |γ t (t)| 2 is assumed (see also Fig. S1) and illustrate its dependence on T in Fig. S2a, where we indeed observe a linear relationship in the small-T regime.
For a system that is initially out of equilibrium, stalling of the observable response and relaxation to the prethermal plateau are then predicted to occur on the time scale t stall on which the dynamics ⟨A⟩ ρ0(t) of the associated unperturbed system H 0 thermalizes.Since no assumptions about the unperturbed system are made, t stall is largely arbitrary and can vary considerably, similarly to the relaxation times of isolated many-body systems.To observe a noticeable effect of the driving and its subsequent stalling, we need t resp < t stall (hence T ≲ t stall ).Moreover, t stall must be considerably smaller than the heating time scale t heat .
Contrary to t resp , this time scale t heat for heating will usually grow as T is decreased.Specifically, approximate laws and rigorous upper bounds for the energy absorption rate per degree of freedom, Γ ∼ t −1 heat , have been established in various lattice systems, for example: for spins or fermions with local interactions in the linear response regime [25] and beyond [24]; for spins with fewbody interactions based on truncated Floquet-Magnus  expansions [23]; or for hard-core bosons by numerical linked-cluster expansions [20].All those works demonstrate that Γ ≤ e −O(1/T ) asymptotically for small T , opening up a large initial time window where heating effects are insignificant on the typical response and stalling time scales if the driving period is sufficiently small.Intuitively, one expects that the system can no longer follow the driving for T → 0, i.e., the driving effects average out to zero for asymptotically fast driving.Moreover, the leading order corrections for finite T are expected to be invariant under a sign change of T , i.e., they generically should scale quadratically with T .Within the the- ory (cf.( 9) in the main paper), we can assess the magnitude of the response as the amplitude of |γ t (t)| 2 at the first minimum.Recalling that γ τ (0) = 1, we thus inspect the quantity and find that γ indeed scales quadratically with T for small values, as shown in Fig. S3a.To achieve a noticeable observable response, one should thus increase the amplitude f 0 of the driving if T is decreased.As illustrated in Fig. S3b, γ likewise scales quadratically with f 0 for fixed (sufficiently small) T , whereas the time scale t resp is largely unaffected by such variations of f 0 (cf.Fig. S2b).Consequently, a decrease of T can be compensated by a proportional increase of f 0 to maintain a similar magnitude of the observable response; see also Fig. 1 in the main paper for a visualization in a concrete example system.The heating rate Γ ∼ t −1 heat , in turn, also grows quadratically with f 0 for small amplitudes within the linear response regime according to Refs.[20,40,41].However, as mentioned above, the observable response will be very weak if both f 0 and T are small.More precisely, except for finite-size effects, we still expect that the relative difference in response between systems near and far from thermal equilibrium will remain significant in this case, but the stalling effect will be less impressive on an absolute scale.
Hence more interesting is the case of stronger driving beyond the linear response regime.Here, general statements about the heating rate Γ are scarce and require more information about the driving operator, but there is evidence that the dependence of Γ on f 0 is often nonmonotonic, such that Γ may decrease again eventually as f 0 becomes larger [40,41,56].In any case, the dependence is typically subexponential.As T is decreased, and even if f 0 is increased accordingly, the exponential suppression of heating in 1/T will thus eventually dominate.For sufficiently small T and large f 0 , we therefore generically expect a regime where heating is insignificant, with a significant response away from equilibrium, but strongly suppressed response in thermal equilibrium.
In conclusion, the parameter regime for stalled response essentially coincides with the one for the theoretically and experimentally well-established phenomenon of Floquet prethermalization [8,20,23,24,33,34,57,58] (see also the discussion in the fifth paragraph of the section "Interpretation and further examples" in the main paper).

General properties of γ τ (t)
Expanding on the discussion of their time scales and amplitudes in the previous subsection, we collect a few general properties of the solutions γ τ (t) of the nonlinear integro-differential equation (10) in the main paper.Some of these properties are more easily understood from alternative representations of γ τ (t), such as (23) in the main paper, or Eq.(S29), defined below in Supplementary Note 6.The equivalence of these representations and (10) in the main paper will be established in Supplementary Note 7 below.The present section merely serves as a convenient overview.
For all τ ∈ R, γ τ (t) satisfies the initial condition and is bounded, see below Eq. (S55).Furthermore, u(E, τ ) from Eq. (S27) is real-valued and, for the considered perturbation ensembles as defined in Supplementary Note 5, an even function of E. Since γ τ (t) is the Fourier transform of u(E, τ ) according to Eq. (S29), it inherits those properties, i.e., Similarly, taking for granted that u(E, τ ) is sufficiently regular, we can conclude that for any fixed τ .We point out, though, that the longtime behavior of γ τ (t) is of minor interest for our present purposes because, as discussed in the Methods of the main paper, the truncated Magnus expansion adopted to derive the theoretical prediction from (9) in the main paper will eventually invalidate it at long times.
As for the short-term behavior, we can immediately infer from Eq. (S4) (or from the integro-differential equation (10) in the main paper) that γτ (0) = 0 . (S6) More generally, we can expand γ τ (t) into a Taylor series around t = 0. To this end, it is convenient to introduce the abbreviation v τ (t) := a τ v(t) + b τ v(t), see also Eq. (S51).Denoting by γ (n) τ (t) the n-th derivative with respect to t, a straightforward calulation then yields the recurrence relation (S7) for the even derivatives, whereas all odd derivatives vanish (see also Eq. (S4)).Up to third order in t, for example, we thus find + O(t 4 ) , (S8) illustrating how the initial decay of γ τ (t) is controlled by the driving and the decay characteristics of v(t).
In (e-h), only results for nonequilibrium initial conditions are shown; the curves for thermal equilibrium initial conditions are almost identical.Parameter values for β ≈ 0.08, 0.12, 0.16, 0.2: E = −12, −18, −24, −30 (target energy of the Gaussian filter, see also Fig. S6b); A th = −0.026,−0.002, 0.026, 0.057 (see also Fig. S6a); ṽ(0)D0 = 3.6, 2.25, 1.44, 1.3 and ∆v = 3.0., 4.2, 5.0, 4.0 (see also third paragraph in "Interpretation and further examples" of the main paper).dle of the spectrum (highest density of states, smallest relative variations of the density of states), but the range of compatible temperatures is expected to increase with the system size due to the exponential growth of the overall number of levels as well as of the level density.Furthermore, we recall that the stalled-response effect is predicted to occur as long as energy absorption from the periodic driving is negligibile.The relative influence of this heating is typically smaller for higher temperatures, too, in the sense that the departure from the initially occupied energy window is smaller if all other parameters (notably T and f 0 ) are kept fixed (see also Fig. S5).Given the lim-its of our computational resources, we therefore focused on relatively high temperatures in the main paper (see also the subsequent Supplementary Note 2.3).
To further substantiate these general arguments, we discuss the behavior for lower temperatures in Fig. S5, using the same two-dimensional spin-lattice system as in Fig. 1 of the main paper.Likewise, the initial states are again generated according to (5) in the main paper, but now using successively lower target energies E = −12 (similar to Fig. 1), −18, −24, −30.As can be inferred from Fig. S6b, these correspond to inverse temperatures β ≈ 0.08, 0.12, 0.16, 0.2, respectively.This Fig. S6b furthermore demonstrates that the energy in the ground state (β → ∞) and at infinite temperature (β = 0) are indeed approximately −60 and −1, respectively, as stated in the main paper above (6).
Returning to Fig. S5, the top panels (a)-(d) show the time-dependent expectation values of the magnetization correlation A = σ z 2,2 σ z 3,3 as before.(Another example -starting from infinite temperature, β = 0 -can be found in Fig. S9 below, cf.Supplementary Note 3.) To visualize the energy absorption in the driven system, we additionally show in the bottom panels (e)-(h) the timedependent expectation values ⟨H 0 ⟩ ρ(t) of the unperturbed reference Hamiltonian H 0 (thick blue line) along with bands of one "standard deviation" [⟨H 2 0 ⟩ ρ(t) − ⟨H 0 ⟩ 2 ρ(t) ] 1/2 (blue shaded region).For comparison, we also indicate the initially occupied energy window (black).
We observe that the driven system is not strictly confined to this initially occupied energy window in any of the examples from Fig. S5.Instead, all of them show FIG.S7.Same as in Fig. 1 of the main paper, but now for a 4 × 4 lattice (see also Supplementary Note 2.3 for more details).
residual heating effects, which manifest themselves by a drift of the mean energy and broadening of the fluctuations, as expected for finite driving period.Notably, however, the departure from the initial energy window is smaller for states of higher temperature (smaller β).In line with this observation and the discussion in the main paper, the observable dynamics ⟨A⟩ ρ(t) (top panels) shows stronger stalling effects for smaller β, too.We point out, however, that the response at very early times, namely the first minimum of ⟨A⟩ ρ(t) at time t ≈ T /2 = 0.1, is strongly suppressed for all displayed temperatures when comparing the thermal equilibrium initial conditions to the nonequilibrium behavior.
If all other parameters are kept fixed, stalled response is thus typically more pronounced at higher temperatures, and most pronounced at early times, because of relatively weaker heating effects, similarly as in the case of higher driving frequencies (cf.Supplementary Note 1.1).Stronger suppression at lower temperatures, in turn, can be expected upon increasing the driving frequency, or upon increasing the system sizes (see next subsection).

Finite-size effects
We recall that our square lattice model from Fig. 1 of the main paper only exhibits a relatively small extension of L = 5 sites along each of the two spatial directions.Hence, notable finite-size effects may still be expected.
To get an idea of their relevance, we consider a smaller version with L = 4 of the same two-dimensional spinlattice system (cf.Eqs.(1), (3), and (4) in the main paper) as before.Again, we employ a sinusoidal driving protocol ((6) in the main paper) and initial states as in (5) of the main paper with Q = π + 2,2 π + 3,3 , but now choosing E = −8 and ∆E = 2 as the target energy window in order to account for the different absolute energy scale of the L = 4 system, and to obtain the same inverse tem-perature β ≈ 0.08 as in the example with L = 5.The thermal expectation value of our observable A = σ z 2,2 σ z 3,3 in this window now assumes the value A th = −0.040.
As far as the theoretical prediction from Eqs. ( 9)-( 10) of the main paper is concerned, an advantage of the smaller system size is that we can calculate the perturbation profile from (7) of the main paper directly by exact diagonalization.We find that it is well approximated by an exponential decay (see also Ref. [32]) of the form ṽ(E) = 5.08 × 10 −3 e −|E|/8. 4 (S11) for eigenstates of H 0 in the relevant energy window, |E µ − E| ≤ ∆E.Furthermore, the mean density of states in this window is D 0 = 425.The Fourier transform of ṽ(E), cf.(8) of the main paper, is thus from which γ τ (t) can be calculated by integrating (10) of the main paper numerically as before.The so-obtained numerical results along with the corresponding theoretical predictions are shown in Fig. S7.As a technical aside, we note that in order to avoid additional finite-size artifacts from the employed dynamicaltypicality method, we averaged over 100 random states |ϕ⟩ (cf.(5) of the main paper).
Comparison of the these numerical results for L = 4 in Fig. S7 with those for L = 5 in Fig. 1 of the main paper provides quite convincing evidence that the stalled response effect should become more pronounced as the system size is increased.In particular, the small remaining driving effects in case of thermal equilibrium initial conditions (red curves) may be expected to become still smaller upon further increasing L, which, however, is computationally infeasible for us in practice.Furthermore, we observe that the theoretical prediction according to (9) from the main paper agrees better with the FIG.S8.Time-dependent expectation values of s z stag from Eq. (S14) (top) and s x from Eq. (S15) (bottom) for the nonintegrable transverse-field Ising model (cf.Fig. 2 of the main paper) with driving amplitude f0 = 4, driving periods T as indicated in each panel, and two different initial states: Néel ordered in the z direction (left, a-d) and fully polarized in the x direction (right, e-h).All other parameters are as in Fig. 2 of the main paper.Insets: Same numerical data, but with rescaled x and y axes to display the long-time behavior.Often, some of the curves are hidden by others, most notably in panels (e) and (f).
numerics (solid lines) for L = 5 than for L = 4, in accordance with the fact that the derivation of the theory assumes large system sizes (see also "Methods" in the main paper).
Finally, we also note that increasing β (cf.Supplementary Note 2.2) seems to have somewhat similar qualitative effects as decreasing the system size L.
For any given such β, this suggests once again that the agreement between numerics and theory, and thus the manifestation of stalled response, would significantly improve if one increased L further.

Special case A = V
As explained in the "Methods" of the main paper, the typicality framework adopted to derive our main analytical result, (9) in the main paper, is not well suited to describe situations in which the observable A is strongly correlated with the driving operator V .To illustrate this and to clarify whether the effect of stalled response per se still applies (as supported by the heuristic arguments provided in the main paper's "Basic physical mechanisms" section), we investigate the case A = V in the following.Since V is an extensive operator, it is appropriate to consider initial states that are globally out of equilibrium (otherwise the observable A = V is expected to exhibit no difference compared to equilibrium initial states for asymptotically large systems).Let us therefore return to our example of the nonintegrable Ising model from Fig. 2 of the main paper, prepared in a (Gaussian-filtered) Néel state, The subscript 'z' here explicitly indicates that the state is expressed with respect to the spin basis in the z direction.We also compare the resulting dynamics to the one obtained from the corresponding thermal equilibrium conditions, emulated as before by choosing |ϕ⟩ to be a Haar-random state in Eq. (S13).
In both cases, we furthermore employ in Eq. (S13) the same parameters E = −2.4 and ∆E = 1 as in Fig. 2 of the main paper.A natural observable in view of the initial state's Néel order is the staggered magnetization in the z direction, As anticipated from its close connection to the singlesite magnetization σ z 1 shown in Fig. 2a of the main paper, the numerics for s z stag in Fig. S8a-b again exhibits stalled response and good agreement with our analytical prediction, (9) of the main paper, for sufficiently small times.
Next we turn to Fig. S8c-d showing the time evolution of the x magnetization, Up to a trivial factor −1/L, which we henceforth ignore, this observable A = s x thus coincides with the driving operator V in our present setup (see also the caption of Fig. 2 in the main paper).Our first observation is that, numerically, the initial response of s x away from equilibrium is weaker than what is seen in s z stag .(Note that we deliberately chose the same y-axis scaling for the plots in panels (a) through (d) to facilitate their direct comparison.)Second, the numerically observed reaction to the driving is significantly stronger away from equilibrium than near thermal equilibrium, as can be seen by comparing the solid blue and red lines in the insets in particular.In other words, the stalled-response phenomenon also manifest itself for the special observable A = V .
Third, the latter applies despite the fact that the theoretical prediction from (9) of the main paper breaks down for A = V , as anticipated there and above.Indeed, since the thermal expectation value s x th = 0 and the unperturbed dynamics ⟨s x ⟩ ρ0(t) = 0 (the associated black lines in Fig. S8c-d are hidden behind the green ones), the theory from (9) of the main paper predicts ⟨s x ⟩ ρ(t) = 0 for the driven dynamics, too (the dashed purple lines are likewise hidden behind the green solid ones), while the actually observed numerical response quite notably deviates from this prediction.The same theoretical predictions also apply to the thermal equilibrium initial conditions, but here the numerically observed deviations from the theory are less severe due to the occurrence of stalling.
Since the unperturbed dynamics for A = V with Néelordered initial conditions are trivial, we consider a second example where the initial state is of the same form (cf. Eq. (S13)), but with |ϕ⟩ = |↑↑ • • • ⟩ x being fully polarized in the x direction.Hence, we obtain an initial value ⟨s x ⟩ ρ(0) that is far from the thermal equilibrium value s x th = 0 (see Fig. S8g-h).
On the other hand, the unperturbed dynamics of the staggered magnetization is now trivial, ⟨s z stag ⟩ ρ0(t) = (s z stag ) th = 0, for both nonequilibrium and equilibrium initial conditions (cf.Fig. S8e-f; black lines hidden behind green ones).Hence the theory from (9) of the main paper predicts ⟨s z stag ⟩ ρ(t) = 0 in both cases.This time, this theoretically predicted behavior is indeed found in the actual dynamics of s z stag , i.e., one essentially observes no response in Fig. S8e-f for both the nonequilibrium initial state (blue lines hidden behind green ones) and the equilibrium one (red lines).
For A = V in Fig. S8g-h, by contrast, the theory fails to predict the correct dynamics in the nonequilibrium setting (solid blue vs. dashed purple lines), but as emphasized repeatedly, this is understood to result from the adopted typicality framework in the derivation of (9) of the main paper.Notwithstanding, we observe that the difference between the unperturbed and driven dynamics is overall larger when the system is away from equilibrium (black and blue lines) compared to the situation close to thermal equilibrium (green and red lines).Therefore, despite the breakdown of our analytical theory, the stalled-response effect itself remains observable, as sup-ported additionally by the heuristic arguments provided in the main paper's "Basic physical mechanisms" section.3), ( 4), ( 6) of the main paper (see also Fig. 1 there) at infinite temperature (β = 0).Solid lines: numerical results for nonequilibrium initial conditions from Eq. (S21) for driving amplitudes f0 = 0 (unperturbed, black), and for driving periods T and amplitudes f0 as indicated in each panel (driven, blue).Dashed lines: corresponding theoretical prediction from (9) of the main paper, where |γt(t)| 2 is obtained as the solution of (10) of the main paper using the approximation from Eq. (S20) with the numerically obtained two-point correlation function ⟨V (t)V ⟩ρ mc from Eq. (S22) as shown in (g).FIG.S10.Same as in Fig. S9, but now for β ≈ 0.08; see also Supplementary Note 3 for more details.
To test this conjecture, we return to the example setup from Fig. 1 of the main paper.We remove the Gaussian filter in the initial conditions (cf.(5) of the main paper), such that the resulting initial state is effectively at infinite temperature, where |ϕ⟩ is a Haar-random state in the S z = −1 magnetization subsector as before.The solid lines in Fig. S9a-f show the corresponding numerically obtained dynamics for the unperturbed (black) and driven (blue) systems.Furthermore, we numerically calculate the two-point correlation function ⟨V (t)V ⟩ ρmc at infinite temperature (ρ mc = 1/ tr(1)) using dynamical typicality [17,65,66].Concretely, we approximate where |ϕ Q (t)⟩ := e −iH0t |ϕ Q ⟩, |ϕ Q ⟩ := Q|ϕ⟩, and |ϕ⟩ is a Haar-random state as before.The so-obtained correlation function is shown in Fig. S9g.Next we use this result to substitute v(t) = ⟨V (t)V ⟩ ρmc /2 in the integro-differential equation (10) from the main paper for γ τ (t).Integrating the equation numerically as before and utilizing the solution for |γ t (t)| 2 in the prediction from (9) of the main paper, we finally obtain the dashed purple lines in Fig. S9a-f.The agreement between theory and numerics is remarkably good throughout the inspected range of amplitudes and driving periods.
As an illustration of finite-temperature effects in (S20), we repeat the entire procedure for the identical setup as in Fig. 1 from the main paper, meaning that the initial state is now again of the same form as in (5) of the main paper with target energy E = −12 and width ∆E = 4 (such that β ≈ 0.08).The corresponding nu-merical results (solid lines) in Fig. S10a-f are thus identical to those in Fig. 1 of the main paper.To estimate the time-correlation function ⟨V (t)V ⟩ ρmc , we follow the dynamical-typicality approach from Eq. (S22) again, but choose |ϕ Q ⟩ := Qe −(H0−E)/4∆E 2 |ϕ⟩ now, thereby emulating the microcanonical density operator of the energy window around the target energy E = −12.This yields Fig. S10g.
The dashed lines in Fig. S10a-f are then again obtained by employing the approximation from Eq. (S20) in (10) of the main paper.Comparing theory and numerics, we observe stronger deviations than in the infinite-temperature case, but we still recover the qualitative features of the response and even achieve reasonable quantitative proximity.

Supplementary Note 4: Relation to Floquet theory
Numerous insights about the long-time behavior of periodically driven systems, including heating effects, metastable plateau regimes ("Floquet prethermalization"), and their topological properties, have been obtained using so-called Floquet theory; see, for example, Refs.[6-8, 10-14, 18-26, 31, 33, 34, 38-42] and references therein.As explained in the main paper, the focus of our present study is on the complementary regime of shortto-intermediate times.Crucially, our methodological approach is distinct from traditional Floquet theory, too.The following discussion is to clarify their relationship.
Floquet theory is based on the mathematical insight that solutions for the propagator U(t) of the timedependent Schrödinger equation d dt U(t) = −iH(t)U(t) can be decomposed as [11,27,67] U(t) = M(t) e −iHFt (S23) if H(t) = H(t + T ) is time periodic.Here H F is a timeindependent Hermitian operator, the so-called Floquet Hamiltonian.Furthermore, M(t) is a time-dependent unitary operator of the same periodicity as the driving protocol, M(t + T ) = M(t), sometimes called the micromotion operator.Since U(0) = 1, one immediately concludes M(nT ) = 1 for all n ∈ Z and thus U(nT ) = e −iHFnT .In other words, the dynamics generated by the time-independent Floquet Hamiltonian H F agrees with the dynamics of the actual system (generated by the time-dependent Hamiltonian H(t)) at integer multiples of the driving period T .(As an aside, the reference time can be chosen arbitrarily, i.e., by a suitable adaptation of H F , one can instead achieve agreement for all t n = nT + θ with an arbitrary, but fixed θ ∈ [0, T ).Among other things, this implies that the choice of H F is not unique, but these technical details are not important for the ensuing discussion.) Exploiting the stroboscopic agreement between the dynamics generated by H F and H(t), the vast majority of studies investigating the long-time behavior focused on properties of the Floquet Hamiltonian H F and largely ignored the periodic modulations induced by the micromotion operator M(t).For example, "Floquet prethermalization" [7,8,20,23,24,26,33,34] describes the observation that, for sufficiently small driving periods T , the dynamics generated by H F resembles ordinary "prethermalization" [68,69]: The system relaxes from its initial nonequilibrium state to some quasistationary intermediate state, and the true thermal equilibrium state is only approached at much later times.In the Floquet case, the quasistationary intermediate state arises as a result of the strong suppression of heating at fast driving.Crucially, the system can thus spend long times close to this nontrivial intermediate state, before heating eventually takes it towards the featureless infinite-temperature ("thermal equilibrium") state.
The periodic modulations by the micromotion operator M(t) are usually disregarded when discussing this effect.Nevertheless, M(t) can generally still induce a strong time dependence of observable expectation values, even if the stroboscopic dynamics generated by H F relaxes to a plateau value.An example is provided in Fig. 3 of the main paper, where the stroboscopic dynamics becomes stationary, but ⟨A⟩ ρ(t) continues to oscillate.
Our principal observation, the phenomenon of stalled response, implies that those periodic modulations are also suppressed if the accompanying unperturbed system finds itself near thermal equilibrium and heating is negligible.Put differently, the micromotion operator M(t) affects states far away from equilibrium more strongly than states close to thermal equilibrium.Coming back to the example from Fig. 3 of the main paper once again, we recall that the unperturbed system there relaxes to an equilibrium state.Crucially, however, this state is different from the thermal equilibrium state of the full system, and therefore M(t) continues to have an effect even though the H F dynamics has settled down.
From a technical point of view, the reason why we can characterize the response continuously in time rather than stroboscopically is that we do not adopt a decomposition like in Eq. (S23).Instead, we always work with the full propagator U(t), particularly when employing the Magnus expansion according to (14) in the main paper.
On the other hand, the fact that we truncate the Magnus expansion at second order means that our characterization of the "prethermal" plateau state is less accurate than state-of-the-art results obtained from highfrequency expansions of stroboscopic Hamiltonians (such as the Floquet Hamiltonian) [23,24,26,33,34].Within our approximation, and if we tacitly assume that f (t) averages to zero over one driving period, the plateau is essentially determined by the time-averaged Hamiltonian H 0 , see also the "Limits of applicability" subsection in the Methods.This is primarily relevant for the longtime expectation value A th in the prediction from (9) of the main paper.In light of the literature on stroboscopic dynamics and observing that H 0 is the first-order ap-
FIG. S2.Dependence of the response time scale tresp, defined as the time of the first minimum of |γt(t)| 2 , on the driving period T and amplitude f0 of the sinusoidal driving f (t) = f0 sin(2πt/T ) (cf. (6) in the main paper).(a) tresp vs. T for various f0; (b) tresp vs. f0 for various T .The function |γt(t)| 2 was evaluated by numerically integrating (10) from the main paper for an exponential perturbation profile ṽ(E) = e −|E| .Dashed line: tresp ∝ T as a guide to the eye.
FIG. S4.Logarithmically plotted Fourier spectrum |ã(ω)| 2 of the nonequilibrium signal a(t) [see Eqs.(S9)-(S10)] for the magnetization correlation A = σ z 2,2 σ z 3,3 in the two-dimensional spin lattice from Fig. 1 of the main paper for (a) nonequilibrium and (b) equilibrium initial conditions.Solid lines: discrete Fourier transformation of the corresponding numerical data in Fig. 1 as indicated in the legend.The time series extend up to t = tmax = 32 with a step size of ∆t = 0.005, yielding a frequency resolution ∆ω = π 16 ≈ 0.2 in the range |ω| ≤ 200π ≈ 630.Vertical dashed lines: multiples of the driving frequency Ω = 2π/T = 10π.
FIG. S6.Thermal expectation values ofA = σ z2,2 σ z 3,3 and H0 in the two-dimensional spin lattice model from Fig.1of the main paper versus inverse temperature β.Depicted are estimates based on dynamical typicality methods[17,65,66], using imaginary-time propagation of three Haar-random states.Shaded region indicates one standard deviation of the mean (not visible in b).