# Finite-frequency prethermalization in periodically driven ergodic systems

Lorenzo Correale  SISSA — International School for Advanced Studies, via
Bonomea 265, I-34136 Trieste, Italy  INFN — Istituto Nazionale di Fisica
Nucleare, Sezione di Trieste, I-34136 Trieste, Italy  Leticia F. Cugliandolo
Sorbonne Université, Laboratoire de Physique Théorique et Hautes Energies,
CNRS UMR 7589, 4 Place Jussieu, 75252 Paris Cedex 05, France  Institut
Universitaire de France, 1 rue Descartes, 75005 Paris France  Marco Schirò
JEIP, UAR 3573 CNRS, Collège de France, PSL Research University, 11 Place
Marcelin Berthelot, 75321 Paris Cedex 05, France  Alessandro Silva  SISSA —
International School for Advanced Studies, via Bonomea 265, I-34136 Trieste,
Italy

###### Abstract

We investigate the periodically driven dynamics of many-body systems, either
classical or quantum, finite-dimensional or mean-field, displaying an
unbounded phase-space. We find that the inclusion of a smooth periodic drive
atop an otherwise ergodic dynamics leads to a long-lived prethermalization,
even at moderate driving frequencies. In specific asymptotic limits, we
compute the corresponding prethermal Hamiltonian from an analytical
perturbation scheme.

_Introduction_ — Periodically-driven systems are known to host a great deal of
non-trivial phenomena. For example, an inverted pendulum can be stabilized by
an external force oscillating at sufficiently high frequency Kapitza (1951),
while a kicked-rotor exhibits an integrability-to-chaos transition Chirikov
(1971) that is counter-intuitively suppressed by the inclusion of quantum
fluctuations Fishman  _et al._ (1982). Similarly, a driven quantum single
particle system may localize in energy space resulting in bounded heating
Gefen and Thouless (1987). While single-particle systems have been the focus
of research in the past decades, attention has recently shifted towards many-
particle ones, which have proven to be valuable tools to reproduce
experimentally a wide array of novel phases of matter , within the field of
_Floquet engineering_ Oka and Kitamura (2019). However, a major obstacle to
their practical realization has been the fact that periodically driven
isolated many-body systems tend to quickly heat up to a featureless infinite-
temperature state, due to the absence of conservation laws D’Alessio and Rigol
(2014); Lazarides  _et al._ (2014); Ponte  _et al._ (2015); Genske and Rosch
(2015).

Infinite-temperature thermalization can be avoided by many-body localization
Khemani  _et al._ (2016); Else  _et al._ (2016); Yao  _et al._ (2017) and
integrability Russomanno  _et al._ (2012); Citro  _et al._ (2015); Gritsev and
Polkovnikov (2017); Ishii  _et al._ (2018); Prosen (1998a, b, 1999). For
generic short-range interacting systems, heating can be suppressed only in a
finite time window. For example, a prethermalization induced by a high-
frequency drive $\Omega$ has been observed in numerical studies both in
classical Howell  _et al._ (2019); Pizzi  _et al._ (2021); Sadia  _et al._
(2022) and quantum systems Kuwahara  _et al._ (2016); Mori  _et al._ (2016);
Abanin  _et al._ (2017a, b); Peronaci  _et al._ (2018, 2020); Ho  _et al._
(2023); Torre and Dentelski (2021), as well as in some experimental setups
Neyenhuis  _et al._ (2017); Rubio-Abadal  _et al._ (2020); Peng  _et al._
(2021); Shkedrov  _et al._ (2022). In this case, the corresponding prethermal
Hamiltonian is asymptotically close to the average over a period,
$\overline{H(t)}$, of the true periodic Hamiltonian $H(t)$. While most
literature has concentrated on systems with a locally bounded energy spectrum,
high-frequency prethermalization has been conjectured as a generic outcome
also for unbounded spectra Hodson and Jarzynski (2021), as supported by
numerical studies on periodically kicked systems Rajak  _et al._ (2018, 2019);
Sadia  _et al._ (2022) and experiments involving a uniform Fermi gas Shkedrov
_et al._ (2022).

From a broader perspective, prethermalisation is also observed close to an
integrable point, both in quenched systems Mori  _et al._ (2018) as well as
under periodic drive Weidinger and Knap (2017). Beyond these scenarios, the
conventional expectation is that a non-integrable system with intermediate-
coupling interactions and drive frequency should rapidly thermalize. In this
work, we show that this expectation is not met when a system with an unbounded
energy spectrum is driven under a smooth periodic force.

We investigate two systems, namely the 1D lattice $\phi^{4}$ model and the
$p$-spin spherical model (PSM). For the parameters we choose, the dynamics of
both models rapidly relax to thermal equilibrium in absence of the drive.
Instead, the coupling to an external drive leads the systems to prethermalize
with an effective Hamiltonian ${\mathcal{H}}_{F}$ which is not just a time-
average of the unperturbed one. For the $\phi^{4}$ model, we calculate
${\mathcal{H}}_{F}$ in certain asymptotic, but not necessarily high-frequency,
limits. Additionally, in both classical and quantum PSMs, an approximate
classical fluctuation-dissipation relation supports the likelihood of
approaching a prethermal finite-temperature state.

Figure 1: One dimensional lattice $\phi^{4}$ model with $N=100$. Data averaged
over $\mathcal{N}=1000$ configurations. (a)-(b) The energy density
$e(t)=\braket{H_{\phi^{4}}\big{(}\mathbf{x}(0),\dot{\mathbf{x}}(0)\big{)}}_{H_{0}}/N$,
obtained by integrating Eqs. (2), with $B_{0}=g=\lambda=\omega_{0}^{2}=1$. (a)
$\Omega=2.4$ and several $e_{0}$. (b) $e_{0}=2$ and several $\Omega$. Inset:
time $t^{*}$ where $e(t)$ crosses the threshold $e_{1}=8$ for the first time,
for several $\Omega$. (c) Comparison between $e(t)$ and the energy density in
Eq. (8) estimated from the effective dynamics in Eqs. (7), with
$e_{0}=\Omega=\tilde{B}_{0}=\tilde{g}=\tilde{\lambda}=1$,
$\Delta\omega^{2}=3$, $\epsilon=0.05$. The parameters of Eq. (2) are
determined by the rescaling in the main text. Inset: time-averaged difference
$\Delta e^{2}$, between the real and the MSA dynamics, as defined in Eq. (9).
To avoid transient effects, we perform the time-average over a window
$[3/4~{}t_{max},t_{max}]$, with $t_{max}$ the maximum simulation time. (d)
$e_{0}$ dependence of the asymptotic energy density $e_{\infty}$, computed as
the average Eq. (12) over the same late time window of the dynamics in Eqs.
(2). The dependence on $\epsilon$ is introduce through
$\tilde{B}_{0},\tilde{g},\tilde{\lambda}$. The dashed line is a microcanonical
average, Eqs. (11) and (10).

_The classical 1D lattice $\phi^{4}$ model_ — The model is

$H_{\phi^{4}}(\mathbf{x},\mathbf{p})=\sum_{i=1}^{N}\bigg{[}\frac{p_{i}^{2}}{2}+\frac{\omega_{0}^{2}}{2}x_{i}^{2}+\frac{\lambda}{4}x_{i}^{4}+\frac{g}{2}(x_{i+1}-x_{i})^{2}\bigg{]}$
(1)

with $\omega_{0}^{2}>0$ and periodic boundary conditions. The normal modes of
its quadratic part have frequencies
$\omega_{k}=\sqrt{\omega_{0}^{2}+2g[1-\cos(2\pi k/N)]}$, for $k=0\ldots N-1$.
In absence of an external drive and for the parameters we consider here, the
Hamilton dynamics displays ergodicity Parisi (1997); Aarts  _et al._ (2000)
(see also the SM SM ) and chaos Pettini and Cerruti-Sola (1991); Pettini
(1993). Several experimental realizations have been proposed for its quantum
extension Shimshoni  _et al._ (2011); Silvi  _et al._ (2013, 2014). We aim to
study the dynamics driven by the term $H_{1}(\textbf{x},t)=-B_{0}\sin(\Omega
t)\sum_{i}x_{i}$. The classical equations of motion are

$\ddot{x}_{i}+\omega_{0}^{2}x_{i}=B_{0}\sin(\Omega
t)+g(x_{i+1}-2x_{i}+x_{i-1})-\lambda x_{i}^{3}\ ,$ (2)

for $i=1,\dots,N$. In Fig. 1 we plot the energy density,
$e(t)=\braket{H_{\phi^{4}}(\textbf{x}(t),\dot{\textbf{x}}(t))}_{0}/N$,
averaged over $\mathcal{N}$ initial conditions randomly sampled on
$H_{\phi^{4}}\big{(}\mathbf{x}(0),\dot{\mathbf{x}}(0)\big{)}=Ne_{0}$, a
constant energy manifold. For all values of $e_{0}$ (a) and $\Omega$ (b),
$e(t)$ saturates to a finite $e_{\infty}$ at long times, signalling
prethermalization. This result remains robust even when choosing a different
$H_{1}(\textbf{x},t)$ SM . By integrating the dynamics over a longer time
scale and for several values of $\lambda$, we observe deviations of $e(t)$
from $e_{\infty}$ which suggest a finite life-time of the plateau controlled
by $\lambda$ SM . When $\Omega$ is roughly outside the spectrum of normal mode
frequencies of Eq. (1), the growth of $e(t)$ slows down. For large $\Omega$,
we quantitatively estimate the slowdown from the time $t^{*}$ at which $e(t)$
crosses a threshold $e_{1}\in[0,e_{\infty}]$. The inset of Fig. 1(b)
illustrates that $t^{*}\sim\exp(c\Omega)$ for sufficiently large $\Omega$.
This dependence resembles the one of the high-frequency prethermalization time
in classical and quantum spin chains Howell  _et al._ (2019); Ho  _et al._
(2023).

The crucial observation that we make is that $e_{\infty}$ typically differs
from the initial energy $e_{0}$. This result suggests that, assuming the
prethermal state to be Gibbs-like, the corresponding effective Hamiltonian
should be distinct from $H_{\phi^{4}}(\mathbf{x},\mathbf{p})$. To gain some
analytical insight on this phenomenon, we study the limit of small driving
amplitude and interaction, by rescaling the parameters as
$B_{0}=\epsilon\tilde{B}_{0}$, $\lambda=\epsilon\tilde{\lambda}$ and
$g=\epsilon\tilde{g}$, for $\epsilon\ll 1$. We also rescale the driving
frequency to $\omega_{0}=\sqrt{\Omega^{2}-\epsilon\Delta\omega^{2}}$, keeping
$\Omega\sim O(1)$, to enable energy absorption through parametric resonances
Landau and Lifshitz (1976) even when Eqs. (2) are close to the ones of non-
interacting harmonic oscillators at frequency $\omega_{0}$ 111Notice that this
is not the high-frequency limit, which instead would require $\Omega\sim O(1)$
and $\omega_{0}^{2}\sim\epsilon$, instead of the scaling
$\Omega^{2}\sim\omega_{0}^{2}+O(\epsilon)$ considered here.. Standard
perturbative methods fail even at small $\epsilon$, due to the proliferation
of secular terms. We circumvent this problem by employing a multiple-scale
analysis (MSA) Bender and Orszag (1999), a method already used to study
single-particle Floquet systems in their high-frequency limit Rahav  _et al._
(2003). Here we briefly review the method, details are postponed to the
Supplemental Material (SM) SM . We begin by introducing a redundant variable
$\tau=\epsilon t$, representing the ‘slow’ time-scale over which the growth of
$e(t)$ is expected. We propose the asymptotic expansion

$x_{i}(t)=x_{0,i}(t,\tau)+\epsilon\ x_{1,i}(t,\tau)+\dots$ (3)

and treat $t$ and $\tau$ as they were _independent_. Although the exact
solution $\mathbf{x}(t)$ depends only on $t$, this trick is useful to
eliminate secular terms at each perturbative order. We plug Eq. (3), together
with the chain rule $d/dt=\partial_{t}+\epsilon\partial_{\tau}$, into Eq. (2).
After collecting powers of $\epsilon$, the first two perturbative orders read:

$\displaystyle\left(\partial_{t}^{2}+\Omega^{2}\right)x_{0,i}$ $\displaystyle\
=\ 0\ ,$ (4) $\displaystyle\left(\partial_{t}^{2}+\Omega^{2}\right)x_{1,i}$
$\displaystyle=-2\partial_{t}\partial_{\tau}x_{0,i}+\Delta\omega^{2}x_{0,i}~{}+$
(5) $\displaystyle\ \ \ -\tilde{\lambda}x_{0,i}^{3}+\tilde{B}_{0}\sin(\Omega
t)~{}.$

The solution of Eq. (4) is given by

$x_{0,i}(t,\tau)=u_{i}(\tau)\cos(\Omega t)+v_{i}(\tau)\sin(\Omega t)~{}.$ (6)

The slowly varying amplitudes $u_{i}(\tau)$ and $v_{i}(\tau)$ are determined
by imposing that secular terms, proportional to $\sin(\Omega t)$ or
$\cos(\Omega t)$, are absent in the right-hand-side of Eq. (5), the next
order, as these would predict a divergence which we do not observe
numerically. After some algebra SM , one finds the Hamiltonian equations

$\frac{du_{i}}{d\tau}=-\frac{\partial\mathcal{H}_{F}}{\partial v_{i}}\
,\qquad\;\;\;\;\frac{dv_{i}}{d\tau}=\frac{\partial\mathcal{H}_{F}}{\partial
u_{i}}\;,$ (7)

generated by the many-body, local effective Hamiltonian

$\displaystyle\\!\\!\\!\\!\mathcal{H}_{F}(\mathbf{u},\mathbf{v})=\frac{1}{2\Omega}\sum_{i}\Big{\\{}\\!-\frac{\Delta\omega^{2}}{2}(u_{i}^{2}+v_{i}^{2})+\frac{\tilde{g}}{2}\Big{[}(u_{i+1}-u_{i})^{2}\;\;$
$\displaystyle\qquad\qquad\quad+(v_{i+1}-v_{i})^{2}\Big{]}+\frac{3\tilde{\lambda}}{16}(u_{i}^{2}+v_{i}^{2})^{2}+\tilde{B}_{0}v_{i}\Big{\\}}\
.$

Plugging the leading order $x_{0,i}(t,\tau)$ into Eq. (1), the estimate of the
energy density is

$\displaystyle e^{\rm MSA}(t)$
$\displaystyle=\frac{1}{N}\braket{H_{\phi^{4}}[\mathbf{x}_{0}(t,\tau),\partial_{t}\mathbf{x}_{0}(t,\tau)]}_{0}\Big{|}_{\tau=\epsilon
t}=$
$\displaystyle=\frac{1}{N}\braket{\mathcal{H}_{0}[\mathbf{u}(\tau),\mathbf{v}(\tau)]}_{0}\Big{|}_{\tau=\epsilon
t}+O(\epsilon)\ .$ (8)

The Hamiltonian
$\mathcal{H}_{0}(\mathbf{u},\mathbf{v})=\Omega^{2}\sum_{i}(u_{i}^{2}+v_{i}^{2})/2$
is obtained omitting the order $\epsilon$ terms in $g$ and $\lambda$ from Eq.
(1). The average of $\mathcal{H}_{0}(\mathbf{u},\mathbf{v})$ is still
performed over microcanonically sampled initial conditions $\mathbf{x}(0)$ and
$\dot{\mathbf{x}}(0)$, corresponding to $\mathbf{u}(0)=\mathbf{x}(0)$ and
$\mathbf{v}(0)=\dot{\mathbf{x}}(0)/\Omega$ from Eq. (6). In Fig. 1(c), we
compare the estimate from Eq. (8) with the results from the full dynamics in
Eqs. (2). The MSA successfully reproduces the dynamics up to a finite time as
well as the saturation to a finite plateau at $e_{\infty}$. Its validity is
further confirmed by an investigation of the time-averaged difference

$\Delta
e^{2}=\lim_{\mathcal{T}\to\infty}\frac{1}{\mathcal{T}}\int_{0}^{\mathcal{T}}dt\;|e(t)-e^{\rm
MSA}(t)|^{2}\ ,$ (9)

between $e(t)$ from the dynamics in Eq. (2), and $e^{\rm MSA}(t)$ from Eq.
(8). In particular, we observe that $\Delta e^{2}$ vanishes for vanishing
$\epsilon$. Our results suggest that the prethermal state that we observe is
indeed generated by an effective Hamiltonian, computed from Eq. (Finite-
frequency prethermalization in periodically driven ergodic systems) at the
leading order, which differs from $H_{\phi^{4}}(\mathbf{x},\mathbf{p})$ 222In
terms of the new variables $\mathbf{u}$ and $\mathbf{v}$, $H_{\phi^{4}}$
coincides with $\mathcal{H}_{0}(\mathbf{u},\mathbf{v})$ at the leading order
in $\epsilon$. even for a small $\epsilon$.

The primary advantage of the MSA is that it quantitatively predicts the value
of $e_{\infty}$, for small $\epsilon$, from equilibrium statistical mechanics.
We first observe that, as $\mathcal{H}_{F}(\mathbf{u},\mathbf{v})$ is
conserved by the dynamics in Eqs. (7), the estimate $e^{\rm MSA}(t)$ in Eq.
(8) is expected to relax to the average

$e^{\rm
MSA}_{\infty}(e_{F})=\braket{\mathcal{H}_{0}(\mathbf{u},\mathbf{v})}_{F}/N\ ,$
(10)

over an _effective_ microcanonical ensemble determined by
$\mathcal{H}_{F}(\mathbf{u},\mathbf{v})=Ne_{F}$. The energy density $e_{F}$ is
in turn determined by the average

$e_{F}(e_{0})=\braket{\mathcal{H}_{F}(\mathbf{u},\mathbf{v})}_{0}/N\ ,$ (11)

computed over a uniform distribution on the surface
$Ne_{0}=\mathcal{H}_{0}(\mathbf{u},\mathbf{v})$, which approximates the
original initial conditions of our protocol, at leading order in $\epsilon$.
Combining Eqs. (10) and (11) we obtain the relationship $e^{\rm
MSA}_{\infty}(e_{0})$. To validate this relation, we compare it to
$e_{\infty}(e_{0})$ obtained from the dynamics in Eq. (2), averaged over a
late time window

$e_{\infty}(e_{0})=\lim_{\mathcal{T}\to\infty}\frac{1}{N\mathcal{T}}\int_{0}^{\mathcal{T}}dt\braket{H_{\phi^{4}}(\mathbf{x}(t),\mathbf{p}(t))}_{0}\
,$ (12)

and initial conditions on the manifold
$H_{\phi^{4}}(\mathbf{x},\mathbf{p})=Ne_{0}$. The comparison, illustrated in
Fig. 1(d), demonstrates that the two relations merge for $\epsilon\to 0$. Our
analysis may be improved by including higher-order terms from Eq. (3).

Figure 2: The classical ($\Gamma=0$) and quantum ($\Gamma=1$) PSM.
$B_{0}=M=J=1$, $p=3$, $t_{b}=40$ and $\Omega=1$. (a) The energy density
$e(t)$, for several values of the bath temperature $T_{0}$. (b) The symmetric
correlation, Eq. (14), evolving from an initial state at temperature
$T_{0}=1.98$. (c) The effective temperature $T_{\text{eff}}$ against $T_{0}$.
The data are obtained through a linear regression of the averaged two-point
functions $d\overline{C}(\tau)/d\tau$ and $\overline{R}(\tau)$ from Eqs. (16)
and (17), using the classical FDT. The corresponding $r^{2}$ parameter is
reported in the inset.

_The classical and quantum fully-connected PSM_ — We can now generalize our
findings to classical and quantum systems of mean-field kind, focusing on the
driven $p$-spin spherical model (PSM):

$\displaystyle\hat{H}_{J}\\!=\\!\sum_{i}\frac{\hat{\Pi}_{i}^{2}}{2M}-\\!\\!\\!\\!\\!\\!\sum_{i_{1}<...<i_{p}}\\!\\!\\!\\!\\!\\!J_{i_{1}...i_{p}}\\!\hat{\sigma}_{i_{1}}...\hat{\sigma}_{i_{p}}\\!-\\!B_{0}\sin(\Omega
t)\sum_{i}\hat{\sigma}_{i}\;$ (13)

The PSM describes a gas of $N$ particles interacting through random, all-to-
all couplings $J_{i_{1},\dots,i_{p}}$, independently sampled from a Gaussian
distribution with zero mean and variance $\overline{J^{2}}=p!J^{2}/2N^{p-1}$.
The particles are globally constrained over an $N$-sphere,
$\sum_{i}\braket{\hat{\sigma}_{i}^{2}}=N$ Kosterlitz  _et al._ (1976). Quantum
fluctuations are implemented by the canonical quantization relations
$[\hat{\sigma}_{i},\hat{\Pi}_{j}]=i\hbar\delta_{ij}$ and are controlled by the
dimensionless parameter $\Gamma=\sqrt{\hbar^{2}/M}$ Cugliandolo  _et al._
(2000). In the following, we will set $M=J=1$. On a critical line
$T_{c}(\Gamma)$, the PSM displays a thermodynamic phase transition between a
paramagnetic and a glassy state, which is either of the first or second order
depending on the strength of $\Gamma$ Cugliandolo  _et al._ (2000, 2001). From
a dynamical point of view, the model is ergodic at high-temperatures and non-
ergodic below a temperature $T_{d}(\Gamma)>T_{c}(\Gamma)$ Cugliandolo and
Lozano (1999). Evidences for dynamical and quantum chaos are retrieved in the
ergodic and non-ergodic phases Bera  _et al._ (2022); Anous and Haehl (2021);
Winer  _et al._ (2022); Correale  _et al._ (2023).

We work at $T>T_{d}(\Gamma)$, where the PSM is ergodic. The symmetric
correlation and linear response function

$\displaystyle C(t,t^{\prime})$
$\displaystyle=\frac{1}{2}\langle\hat{\sigma}(t)\hat{\sigma}(t^{\prime})+\hat{\sigma}(t^{\prime})\hat{\sigma}(t)\rangle\;,$
(14) $\displaystyle R(t,t^{\prime})$
$\displaystyle=\frac{i}{\hbar}\theta(t-t^{\prime})\langle[\hat{\sigma}(t),\hat{\sigma}(t^{\prime})]\rangle\
,$ (15)

obey a closed set of Schwinger-Dyson equations, valid for $N\to\infty$ and
under disorder average Cugliandolo and Lozano (1999); Thomson  _et al._
(2020); Berthier  _et al._ (2001); Busiello  _et al._ (2004) (see the SM SM ).
We initialize the dynamics by evolving the undriven ($B_{0}=0$) system in
contact with a bath of harmonic oscillators at temperature $T_{0}>T_{d}$
during a finite time $0<t<t_{b}$ Thomson  _et al._ (2020). At $t_{b}$ we
switch off the coupling to the bath and we follow the periodically driven
dynamics ($B_{0}>0$) of the isolated system.

In Fig. 2(a) we plot the average energy density $e(t)$
333$\lim_{N\to\infty}e(t)$ can be expressed in terms of $C(t,t^{\prime})$ and
$R(t,t^{\prime})$ (see SM SM ). for $\Omega=1$ and several $T_{0}$. $e(t)$
saturates to a first plateau for $t\lesssim t_{b}$, signalling that the system
reached thermal equilibrium with the bath, and later oscillates around a
second one at $e_{\infty}$. The latter stroboscopic saturation is compatible
with a finite temperature thermalization for both classical and quantum PSM.
Like in the lattice $\phi^{4}$ model, the late saturation of $e(t)$ is robust
against changes in $\Omega$ and the form of the driving force SM .

These results are consistent with the same prethermalization which we also
observed in the $\phi^{4}$. However, we cannot observe deviations of $e(t)$
from $e_{\infty}$ within our simulation time scales. We also observe that,
unlike for the $\phi^{4}$ model, we cannot perform a multi-scale analysis, due
to difficulties posed by the spherical constraint. Instead, as we have direct
access to the two-point functions, we search for other indications of
relaxation to a Gibbs-like state from the Fluctuation-Dissipation Theorem
(FDT) Nielsen and Dyre (1996). Its classical ($\Gamma=0$) version states that
whenever a system is in equilibrium at temperature $T$,
$dC(\tau)/d\tau=-T~{}R(\tau)$ with $\tau=t-t^{\prime}>0$. Although the
dynamics are non-stationary due to the driving, $C$ and $R$ exhibit an
_approximate_ discrete stationarity in the late-time regime, corresponding to
the invariance under simultaneous translations of $t$ and $t^{\prime}$ by the
period $2\pi/\Omega$. Motivated by this observation, illustrated in Fig. 2(b),
we test the FDT for the averaged two-point functions:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\overline{C}(\tau)=\frac{\Omega}{4\pi}\int_{t_{0}}^{t_{0}+\tfrac{2\pi}{\Omega}}\\!\\!\\!\\!ds\left[C(s+\tau,s)-m(s+\tau)m(s)\right]$
(16)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\overline{R}(\tau)=\frac{\Omega}{4\pi}\int_{t_{0}}^{t_{0}+\tfrac{2\pi}{\Omega}}\\!\\!\\!\\!ds\;R(s+\tau,s)\;,$
(17)

with $m(t)=\int_{0}^{t}ds\;R(t,s)B(s)$ the oscillating magnetization SM . From
the slope of a linear regression of $\overline{R}(\tau)$ against
$d\overline{C}(\tau)/d\tau$, we obtain an _effective temperature_
$T_{\text{eff}}$ at late times. The $T_{0}$ dependence of $T_{\text{eff}}$ is
displayed, along with the corresponding $r^{2}$ parameter from the regression,
in Fig. 2(c). The profile of $T_{\text{eff}}$ is similar for the classical and
quantum cases and resembles the behavior observed in $e_{\infty}(e_{0})$ for
the $\phi^{4}$ model. Notably, we find that $r^{2}\simeq 1$ for both classical
and quantum models, indicating that the classical FDT effectively
characterizes the relation between $\overline{C}$ and $\overline{R}$, even in
the latter case.

_Discussion and conclusion_ — In conclusion, our investigation of the dynamics
of ergodic systems, with unbounded phase space and driven by a smooth periodic
force, revealed clear signatures of prethermalization to a Gibbs-like state.
This is indicated either by a saturation of the energy density to a finite
value or by the emergence of an approximate fluctuation-dissipation relation.
For the $\phi^{4}$ model, we also proposed a perturbative scheme which allowed
us to determine the corresponding thermal asymptotic state. Differently from
previous studies, prethermalization occurs at all frequencies and the
corresponding effective Hamiltonian ${\mathcal{H}}_{F}$ is generically
distinct from the time-averaged one, $\overline{H(t)}$. A deeper understanding
of the connection between ${\mathcal{H}}_{F}$ and $\overline{H(t)}$ could be
crucial for advancing the experimental realization of long-lived many-body
Hamiltonians even at moderate frequencies. Joanes and Gill (1998)

_Acknowledgements —_ We thank D. Abanin, A. Delmonte, A. Polkovnikov and A.
Russomanno for discussions. L.C. acknowledges hospitality from the Institute
of Physics of Collège de France. L.F.C. acknowledges support from the grant
ANR-19-CE30-0014. A.S. acknowledges funding from the grants PNRR MUR project
PE0000023-NQSTI and Quantera project SuperLink. M.S. acknowledges funding from
the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (Grant agreement No. 101002955 – CONQUER).

## References

  * Kapitza (1951) P. Kapitza, Uspekhi Fizicheskih Nauk 44, 7 (1951).
  * Chirikov (1971) B. V. Chirikov,  _Research concerning the theory of non-linear resonance and stochasticity_ (CERN, Geneva, 1971) translated at CERN from the Russian (IYAF-267-TRANS-E).
  * Fishman  _et al._ (1982) S. Fishman, D. R. Grempel, and R. E. Prange, Phys. Rev. Lett. 49, 509 (1982).
  * Gefen and Thouless (1987) Y. Gefen and D. J. Thouless, Phys. Rev. Lett. 59, 1752 (1987).
  * Oka and Kitamura (2019) T. Oka and S. Kitamura, Annu. Rev. Condens. Matter Phys. 10, 387 (2019).
  * D’Alessio and Rigol (2014) L. D’Alessio and M. Rigol, Phys. Rev. X 4, 041048 (2014).
  * Lazarides  _et al._ (2014) A. Lazarides, A. Das, and R. Moessner, Phys. Rev. E 90, 030402 (2014).
  * Ponte  _et al._ (2015) P. Ponte, A. Chandran, Z. Papić, and D. A. Abanin, Ann. Phys. 353, 196 (2015).
  * Genske and Rosch (2015) M. Genske and A. Rosch, Phys. Rev. A 92, 062108 (2015).
  * Khemani  _et al._ (2016) V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 116, 250401 (2016).
  * Else  _et al._ (2016) D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. 117, 090402 (2016).
  * Yao  _et al._ (2017) N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Phys. Rev. Lett. 118, 030401 (2017).
  * Russomanno  _et al._ (2012) A. Russomanno, A. Silva, and G. E. Santoro, Phys. Rev. Lett. 109, 257201 (2012).
  * Citro  _et al._ (2015) R. Citro, E. G. Dalla Torre, L. D’Alessio, A. Polkovnikov, M. Babadi, T. Oka, and E. Demler, Ann. Phys. 360, 694 (2015).
  * Gritsev and Polkovnikov (2017) V. Gritsev and A. Polkovnikov, SciPost Phys. 2, 021 (2017).
  * Ishii  _et al._ (2018) T. Ishii, T. Kuwahara, T. Mori, and N. Hatano, Phys. Rev. Lett. 120, 220602 (2018).
  * Prosen (1998a) T. Prosen, Phys. Rev. Lett. 80, 1808 (1998a).
  * Prosen (1998b) T. Prosen, J. Phys. A: Math. Gen. 31, L645 (1998b).
  * Prosen (1999) T. Prosen, Phys. Rev. E 60, 3949 (1999).
  * Howell  _et al._ (2019) O. Howell, P. Weinberg, D. Sels, A. Polkovnikov, and M. Bukov, Phys. Rev. Lett. 122, 010602 (2019).
  * Pizzi  _et al._ (2021) A. Pizzi, A. Nunnenkamp, and J. Knolle, Phys. Rev. Lett. 127, 140602 (2021).
  * Sadia  _et al._ (2022) Y. Sadia, E. G. Dalla Torre, and A. Rajak, Phys. Rev. B 105, 184302 (2022).
  * Kuwahara  _et al._ (2016) T. Kuwahara, T. Mori, and K. Saito, Ann. Phys. 367, 96 (2016).
  * Mori  _et al._ (2016) T. Mori, T. Kuwahara, and K. Saito, Phys. Rev. Lett. 116, 120401 (2016).
  * Abanin  _et al._ (2017a) D. Abanin, W. D. Roeck, W. W. Ho, and F. Huveneers, Commun. Math. Phys. 354, 809 (2017a).
  * Abanin  _et al._ (2017b) D. A. Abanin, W. D. Roeck, W. W. Ho, and F. Huveneers, Phys. Rev. B 95, 014112 (2017b).
  * Peronaci  _et al._ (2018) F. Peronaci, M. Schiró, and O. Parcollet, Phys. Rev. Lett. 120, 197601 (2018).
  * Peronaci  _et al._ (2020) F. Peronaci, O. Parcollet, and M. Schiró, Phys. Rev. B 101, 161101 (2020).
  * Ho  _et al._ (2023) W. W. Ho, T. Mori, D. A. Abanin, and E. G. D. Torre, Ann. Phys. 454, 169297 (2023).
  * Torre and Dentelski (2021) E. G. D. Torre and D. Dentelski, SciPost Phys. 11, 040 (2021).
  * Neyenhuis  _et al._ (2017) B. Neyenhuis, J. Zhang, P. W. Hess, J. Smith, A. C. Lee, P. Richerme, Z.-X. Gong, A. V. Gorshkov, and C. Monroe, Sci. Adv. 3, e1700672 (2017).
  * Rubio-Abadal  _et al._ (2020) A. Rubio-Abadal, M. Ippoliti, S. Hollerith, D. Wei, J. Rui, S. L. Sondhi, V. Khemani, C. Gross, and I. Bloch, Phys. Rev. X 10, 021044 (2020).
  * Peng  _et al._ (2021) P. Peng, C. Yin, X. Huang, C. Ramanathan, and P. Cappellaro, Nat. Phys. 17, 444 (2021).
  * Shkedrov  _et al._ (2022) C. Shkedrov, M. Menashes, G. Ness, A. Vainbaum, E. Altman, and Y. Sagi, Phys. Rev. X 12, 011041 (2022).
  * Hodson and Jarzynski (2021) W. Hodson and C. Jarzynski, Phys. Rev. Res. 3, 013219 (2021).
  * Rajak  _et al._ (2018) A. Rajak, R. Citro, and E. G. D. Torre, J. Phys. A: Math. Theor. 51, 465001 (2018).
  * Rajak  _et al._ (2019) A. Rajak, I. Dana, and E. G. D. Torre, Phys. Rev. B 100, 100302(R) (2019).
  * Mori  _et al._ (2018) T. Mori, T. N. Ikeda, E. Kaminishi, and M. Ueda, J. Phys. B: At. Mol. Opt. 51, 112001 (2018).
  * Weidinger and Knap (2017) S. A. Weidinger and M. Knap, Sci. Rep. 7, 45382 (2017).
  * Parisi (1997) G. Parisi, Europhysics Letters (EPL) 40, 357 (1997).
  * Aarts  _et al._ (2000) G. Aarts, G. F. Bonini, and C. Wetterich, Nucl. Phys. B 587, 403 (2000).
  * (42) See Supplemental Material.
  * Pettini and Cerruti-Sola (1991) M. Pettini and M. Cerruti-Sola, Phys. Rev. A 44, 975 (1991).
  * Pettini (1993) M. Pettini, Phys. Rev. E 47, 828 (1993).
  * Shimshoni  _et al._ (2011) E. Shimshoni, G. Morigi, and S. Fishman, Phys. Rev. Lett. 106, 010401 (2011).
  * Silvi  _et al._ (2013) P. Silvi, G. D. Chiara, T. Calarco, G. Morigi, and S. Montangero, Ann. Phys. 525, 827 (2013).
  * Silvi  _et al._ (2014) P. Silvi, T. Calarco, G. Morigi, and S. Montangero, Phys. Rev. B 89, 094103 (2014).
  * Landau and Lifshitz (1976) L. Landau and E. Lifshitz,  _Mechanics_ (Elsevier, 1976).
  * Note (1) Notice that this is not the high-frequency limit, which instead would require $\Omega\sim O(1)$ and $\omega_{0}^{2}\sim\epsilon$, instead of the scaling $\Omega^{2}\sim\omega_{0}^{2}+O(\epsilon)$ considered here.
  * Bender and Orszag (1999) C. M. Bender and S. A. Orszag,  _Advanced Mathematical Methods for Scientists and Engineers I_ (Springer New York, 1999).
  * Rahav  _et al._ (2003) S. Rahav, I. Gilary, and S. Fishman, Phys. Rev. A 68, 013820 (2003).
  * Note (2) In terms of the new variables $\mathbf{u}$ and $\mathbf{v}$, $H_{\phi^{4}}$ coincides with $\mathcal{H}_{0}(\mathbf{u},\mathbf{v})$ at the leading order in $\epsilon$.
  * Kosterlitz  _et al._ (1976) J. M. Kosterlitz, D. J. Thouless, and R. C. Jones, Phys. Rev. Lett. 36, 1217 (1976).
  * Cugliandolo  _et al._ (2000) L. F. Cugliandolo, D. R. Grempel, and C. A. da Silva Santos, Phys. Rev. Lett. 85, 2589 (2000).
  * Cugliandolo  _et al._ (2001) L. F. Cugliandolo, D. R. Grempel, and C. A. da Silva Santos, Phys. Rev. B 64, 014403 (2001).
  * Cugliandolo and Lozano (1999) L. F. Cugliandolo and G. Lozano, Phys. Rev. B 59, 915 (1999).
  * Bera  _et al._ (2022) S. Bera, K. Y. Venkata Lokesh, and S. Banerjee, Phys. Rev. Lett. 128, 115302 (2022).
  * Anous and Haehl (2021) T. Anous and F. M. Haehl, J. Stat. Mech. 2021, 113101 (2021).
  * Winer  _et al._ (2022) M. Winer, R. Barney, C. L. Baldwin, V. Galitski, and B. Swingle, J. High Energ. Phys. 2022, 32 (2022).
  * Correale  _et al._ (2023) L. Correale, A. Polkovnikov, M. Schirò, and A. Silva, SciPost Phys. 15, 190 (2023).
  * Thomson  _et al._ (2020) S. Thomson, P. Urbani, and M. Schiró, Phys. Rev. Lett. 125, 120602 (2020).
  * Berthier  _et al._ (2001) L. Berthier, L. F. Cugliandolo, and J. L. Iguain, Phys. Rev. E 63, 051302 (2001).
  * Busiello  _et al._ (2004) G. Busiello, R. V. Saburova, V. G. Sushkova, and G. P. Chugunova, Phys. Solid State 46, 316 (2004).
  * Note (3) $\lim_{N\to\infty}e(t)$ can be expressed in terms of $C(t,t^{\prime})$ and $R(t,t^{\prime})$ (see SM SM ).
  * Nielsen and Dyre (1996) J. K. Nielsen and J. C. Dyre, Phys. Rev. B 54, 15754 (1996).
  * Joanes and Gill (1998) D. N. Joanes and C. A. Gill, Journal of the Royal Statistical Society: Series D (The Statistician) 47, 183–189 (1998).
  * Cugliandolo  _et al._ (2017) L. F. Cugliandolo, G. S. Lozano, and E. N. Nessi, J. Stat. Mech. 2017, 083301 (2017).
  * Kamenev (2011) A. Kamenev,  _Field Theory of Non-Equilibrium Systems_ (Cambridge University Press, 2011).
  * Haldar  _et al._ (2020) A. Haldar, P. Haldar, S. Bera, I. Mandal, and S. Banerjee, Phys. Rev. Res. 2, 013307 (2020).

##  Appendix A Supplemental Material: Finite-frequency prethermalization
in periodically driven ergodic systems

In Section 1, we show that the driven dynamics of the classical one-
dimensional lattice $\phi^{4}$ model deviates from its prethermal plateau when
integrated at longer times and discuss the dependence of prethermalization on
the degree of non-linearity in the system. In Section 2, we provide a
comprehensive exposition of the multiple-scale analysis, applied to the
lattice $\phi^{4}$ model. In Section 3, we write down and derive explicitly
the mode coupling equations, which govern the dynamics of the correlation and
response functions, in the thermodynamic limit of the $p$-spin spherical
model. In Section 4, we discuss the ergodic properties of both the two models,
in absence of a time-dependent drive. Finally, in Section 5, we study the
Floquet dynamics of both models, driven by smooth periodic forces different
from the ones considered in the Letter.

##  Appendix B 1\. Deviations from the prethermal plateau in the classical
lattice $\phi^{4}$ model

In this section, we extend the integration of Eqs. (19) over a longer time
scale to investigate the behavior of the classical lattice $\phi^{4}$ model.
We compute the resulting average energy density
$e(t)=\braket{H_{\phi^{4}}(\textbf{x}(t),\dot{\textbf{x}}(t))}_{0}/N$, where
$H_{\phi^{4}}(\textbf{x}(t),\textbf{p}(t))$ is defined in Eq. (1) of the
Letter. As in the Letter, the initial condition is chosen as ensemble of
microcanonical configurations, sampled from the surface
$H_{\phi^{4}}\big{(}\mathbf{x}(0),\dot{\mathbf{x}}(0)\big{)}=Ne_{0}$ for a
fixed initial energy density $e_{0}$. To integrate up to a larger maximum
time, we fix $\mathcal{N}$ to a smaller value then in the Letter. We integrate
the dynamics for different values of $\lambda$, the parameter controlling the
non-linearity of the system. From Fig. 3-(a), it is already evident that the
energy density $e(t)$ displays deviations from its finite plateau at
$e_{\infty}$ at sufficiently long-time, confirming the prethermal nature of
the plateau. We also observe that the height of the plateau increases by
decreasing $\lambda$.

To qualitatively assess the influence of the non-linearity on the prethermal
plateau, we introduce the normalized quantity

$Q(t)=\frac{e(t)-e_{0}}{e_{\infty}-e_{0}}\ ,$ (18)

which saturates to $1$ when $e(t)$ saturates to $e_{\infty}$ and such that
$Q(0)=0$. Here, we estimate $e_{\infty}$ as $e(t=t_{0})$ for a fixed
intermediate time $t_{0}$. As depicted in Fig. 3-(b), deviations of $Q(t)$
from its plateau appear at earlier times as $\lambda$ increases, suggesting
that a stronger non-linearity suppresses prethermalization.

Figure 3: Average energy densities for the one dimensional lattice $\phi^{4}$
model, evolving from Eqs. (87), with $N=100$. Data are averaged over
$\mathcal{N}=50$ configurations. Other parameters are set as
$\omega_{0}^{2}=B_{0}=g=1$, $e_{0}=1$ and $\Omega=2$. (a)
$e(t)=\braket{H_{\phi^{4}}(\textbf{x}(t),\textbf{p}(t))}_{0}/N$. (b) $Q(t)$
from Eq. (18). We estimate $e_{\infty}$ as $e(t_{0}=100)$.

##  Appendix C 2\. Detailed calculations on multiple-scale analysis

In this section, we aim to provide a review on the multiple-scale analysis
(MSA) that we employed in the Letter. We will focus again on the Hamilton
dynamics of the classical lattice $\phi^{4}$ model, governed by the equations

$\ddot{x}_{i}+\omega_{0}^{2}x_{i}+\lambda x_{i}^{3}=B_{0}\sin(\Omega t)+g\
(x_{i+1}-2x_{i}+x_{i-1})\ ,$ (19)

for $i=1\ldots N$. We rescale the model parameters as
$B_{0}=\epsilon\tilde{B}_{0}$, $\lambda=\epsilon\tilde{\lambda}$ and
$g=\epsilon\tilde{g}$. We also rescale the harmonic coupling as
$\omega_{0}^{2}=\Omega^{2}-\epsilon\Delta\omega^{2}$. We will work in the
limit of small $\epsilon$, while keeping all other parameters of order one. As
outlined in the Letter, the starting point of the MSA consists in introducing
a redundant time variable $\tau=\epsilon t$. $\tau$ defines a long time scale,
since it becomes significant only when $t$ is of order $\epsilon^{-1}$ or
larger. Although we expect that the actual solution $\mathbf{x}(t)$ to Eq.
(19) is a function of $t$ alone, MSA seeks solutions which are functions of
both variables $t$ and $\tau$, treated as independent from each other. Such
expression of $\mathbf{x}(t)$ as a function of two variables is an artifice to
remove secular effects; the actual solution has $t$ and $\tau$ related by
$\tau=\epsilon t$. The formal procedure consists in assuming a perturbative
expansion in powers of $\epsilon$:

$x_{i}(t)=x_{0,i}(t,\tau)+\epsilon\ x_{1,i}(t,\tau)+\dots\ ,$ (20)

for small $\epsilon$. We use the chain rules
$d/dt=\partial_{t}+\epsilon\partial_{\tau}$ for partial differentiation to
compute the derivatives of the variables $x_{i}(t)$, obtaining

$\frac{dx_{i}}{dt}=\frac{\partial x_{0,i}}{\partial
t}+\epsilon\Big{(}\frac{\partial x_{1,i}}{\partial t}+\frac{\partial
x_{0,i}}{\partial\tau}\Big{)}+O(\epsilon^{2})$ (21)

and

$\frac{d^{2}x_{i}}{dt^{2}}=\frac{\partial^{2}x_{0,i}}{\partial
t^{2}}+\epsilon\Big{(}\frac{\partial^{2}x_{1,i}}{\partial t^{2}}+2\
\frac{\partial^{2}x_{0,i}}{\partial t\ \partial\tau}\Big{)}+O(\epsilon^{2})\
.$ (22)

By substituting Eq. (22) in Eq. (19) and collecting powers of $\epsilon$, we
deduce

$\displaystyle\frac{\partial^{2}x_{0,i}}{\partial t^{2}}+\Omega^{2}x_{0,i}=0\
,$ (23) $\displaystyle\frac{\partial^{2}x_{1,i}}{\partial
t^{2}}+\Omega^{2}x_{1,i}=-2\frac{\partial^{2}x_{0,i}}{\partial t\
\partial\tau}+\Delta\omega^{2}x_{0,i}+\tilde{g}(x_{0,i+1}-2x_{0,i}+x_{0,i-1})-\tilde{\lambda}x_{0,i}^{3}+\tilde{B}_{0}\sin(\Omega
t)\ .$ (24)

The general solution to Eq. (23) is

$x_{0,i}(t,\tau)=A_{i}(\tau)e^{i\Omega t}+A_{i}^{*}(\tau)e^{-i\Omega t}\ .$
(25)

By plugging this expression into Eq. (24), we derive

$\frac{\partial^{2}x_{1,i}}{\partial
t^{2}}+\Omega^{2}x_{1,i}=2~{}\text{Re}\Big{[}e^{i\Omega
t}\big{(}-2i\Omega\frac{dA_{i}}{d\tau}+\Delta\omega^{2}A_{i}+\tilde{g}(A_{i+1}-2A_{i}+A_{i-1})-3\tilde{\lambda}A_{i}^{2}A_{i}^{*}-\frac{i}{2}\tilde{B}_{0}\big{)}-\tilde{\lambda}e^{i3\Omega
t}A_{i}^{3}\Big{]}\ .$ (26)

The slowly varying amplitude $A_{i}(\tau)$ and its complex conjugate
$A_{i}^{*}(\tau)$ are determined by imposing that secular terms, proportional
either to $e^{i\Omega t}$ or $e^{-i\Omega t}$, do not appear in the next
order, that is on the right-hand-side of Eq. (24), as these would predict a
divergence which we do not observe in the complete numerical solution of this
problem. The vanishing of such secular terms thus leads to the following
equation of motion for $A_{i}(\tau)$:

$\frac{dA_{i}}{d\tau}=\frac{1}{2i\Omega}\left[\Delta\omega^{2}A_{i}-3\tilde{\lambda}A_{i}^{2}A_{i}^{*}-\frac{i}{2}\tilde{B}_{0}+\tilde{g}(A_{i+1}-2A_{i}+A_{i-1})\right]\
.$ (27)

The equation of motion for $A_{i}^{*}(\tau)$ is just the complex conjugate of
Eq. (27). We rewrite Eq. (27) in terms of the coordinates
$u_{i}(\tau)=A_{i}(\tau)+A_{i}^{*}(\tau)$ and
$v_{i}(\tau)=[A_{i}(\tau)-A_{i}^{*}(\tau)]/i$, which consequently obey the
equations:

$\left\\{\begin{split}\frac{du_{i}}{d\tau}&=\frac{1}{2\Omega}\left[\Delta\omega^{2}v_{i}+\tilde{g}(v_{i+1}-2v_{i}+v_{i-1})-\frac{3}{4}\tilde{\lambda}(u_{i}^{2}+v_{i}^{2})v_{i}-\tilde{B}_{0}\right]\;,\\\
\frac{dv_{i}}{d\tau}&=\frac{1}{2\Omega}\left[-\Delta\omega^{2}u_{i}-\tilde{g}(u_{i+1}-2u_{i}+u_{i-1})+\frac{3}{4}\tilde{\lambda}(u_{i}^{2}+v_{i}^{2})u_{i}\right]\
.\end{split}\right.$ (28)

In terms of $u_{i}(\tau)$ and $v_{i}(\tau)$, the leading order solution, which
satisfies Eq. (23), reads $x_{0,i}(t)=u_{i}(\tau)\cos(\Omega
t)+v_{i}(\tau)\sin(\Omega t)$. The system of Eqs. (28) is Hamiltonian, as it
can be rewritten as

$\left\\{\begin{split}\frac{du_{i}}{d\tau}&=-\frac{\partial\mathcal{H}_{F}}{\partial
v_{i}}\;,\\\ \frac{dv_{i}}{d\tau}&=\frac{\partial\mathcal{H}_{F}}{\partial
u_{i}}\;,\end{split}\right.$ (29)

with the many-body effective Hamiltonian:

$\mathcal{H}_{F}(\mathbf{u},\mathbf{v})=\frac{1}{2\Omega}\sum_{i}\Big{\\{}-\frac{\Delta\omega^{2}}{2}(u_{i}^{2}+v_{i}^{2})+\frac{\tilde{g}}{2}\left[(u_{i+1}-u_{i})^{2}+(v_{i+1}-v_{i})^{2}\right]+\frac{3\tilde{\lambda}}{16}(u_{i}^{2}+v_{i}^{2})^{2}+\tilde{B}_{0}v_{i}\Big{\\}}\
.$

##  Appendix D 3\. Dynamical equations for the correlation and response
function in the p-spin spherical model

The aim of this section is to derive and discuss the closed set of equations
of motions governing the dynamics of the correlation and response functions,
$C(t,t^{\prime})$ and $R(t,t^{\prime})$, of the $p$-spin spherical model (PSM)
investigated in the Letter. We anticipate that, in analogy to the results
obtained in Refs. Berthier  _et al._ (2001); Busiello  _et al._ (2004), which
studied the periodically driven dynamics of the PSM coupled to a thermal bath,
the equations describing its Floquet dynamics have Schwinwger-Dyson structure
and are given by

$\displaystyle[M\partial_{t}^{2}+z(t)]R(t,t^{\prime})$
$\displaystyle=\delta(t-t^{\prime})+\int_{0}^{t}dt^{\prime\prime}\;\Sigma(t,t^{\prime\prime})R(t^{\prime\prime},t^{\prime})\;,$
$\displaystyle[M\partial_{t}^{2}+z(t)]C(t,t^{\prime})$
$\displaystyle=\int_{0}^{t}dt^{\prime\prime}\;\Sigma(t,t^{\prime\prime})C(t^{\prime\prime},t^{\prime})+\int_{0}^{t^{\prime}}dt^{\prime\prime}\;D(t,t^{\prime\prime})R(t^{\prime},t^{\prime\prime})\
+$ (30)
$\displaystyle+B(t)\int_{0}^{t^{\prime}}dt^{\prime\prime}\;R(t^{\prime},t^{\prime\prime})B(t^{\prime\prime})\;.$
(31)

We will define the self-energy $\Sigma(t,t^{\prime})$ and the vertex
$D(t,t^{\prime})$ later on in this section. The time-dependent Lagrange
multiplier $z(t)$ enforces the spherical constraint and has an expression in
terms of $C$ and $R$ as well. In the Letter, we integrate numerically Eqs. (D)
and (31) to compute the correlation and response function, evolving from the
equal-time conditions:

$\displaystyle C(t,t)$ $\displaystyle=1\
,\quad\quad\quad\quad\partial_{t}C(t,t^{\prime})|_{t^{\prime}\to
t^{-}}=\partial_{t}C(t,t^{\prime})|_{t^{\prime}\to t^{+}}=0\ ,$ (32)
$\displaystyle R(t,t)$ $\displaystyle=0\
,\quad\quad\quad\quad\partial_{t}R(t,t^{\prime})|_{t^{\prime}\to
t^{-}}=\frac{1}{M}\
,\quad\quad\quad\quad\partial_{t}R(t,t^{\prime})|_{t^{\prime}\to t^{+}}=0\ ,$
(33)

resulting from the spherical constraint and the causal structure of the
response function. Then, following Refs. Thomson  _et al._ (2020); Cugliandolo
_et al._ (2017), we compute the energy density of the system as

$e(t)=-\frac{M}{2}\partial_{t}^{2}C(t,t^{\prime})\big{|}_{t^{\prime}\to
t^{-}}-\frac{1}{p}\int_{0}^{t}dt^{\prime\prime}\;\left[\Sigma(t,t^{\prime\prime})C(t,t^{\prime\prime})+D(t,t^{\prime\prime})R(t,t^{\prime\prime})\right]\
.$ (34)

This section is organized as follows. In Section E, we derive Eqs. (D) and
(31) using the Schwinger-Keldysh formalism in the large $N$ limit. In Section
F, we show how to determine the Lagrange multiplier $z(t)$ explicitly, in
terms of the correlation and response functions. Finally, in Section G, we
discuss the predictor-corrector algorithm, that we use to integrate Eqs. (D)
and (31) numerically.

##  Appendix E 3.1 Derivation of the Schwinger-Dyson equations

We begin by reminding that the Hamiltonian of the driven quantum PSM is
defined as

$\hat{H}_{J}=\frac{1}{2M}\sum_{i=1}^{N}\
\hat{\Pi}_{i}^{2}-\mathcal{J}(t)\sum_{i_{1}<\dots<i_{p}}J_{i_{1},\dots,i_{p}}\hat{\sigma}_{i_{1}}\cdots\hat{\sigma}_{i_{p}}-B(t)\sum_{i}\hat{\sigma}_{i}$
(35)

with all-to-all couplings $J_{i_{1},\dots,i_{p}}$, independently sampled from
a Gaussian distribution with zero mean and variance
$\overline{J^{2}}=2J^{2}p!/N^{p-1}$. The spins $\hat{\sigma}_{i}$ obey the
usual spherical constraint $\sum_{i}\braket{\hat{\sigma}_{i}^{2}}=N$ (on
average). The canonical commutation relations
$[\hat{\sigma}_{j},\hat{\Pi}_{k}]=i\hbar\delta_{jk}$ hold. We added a time-
dependent factor ${\mathcal{J}}(t)$ to allow for other kinds of drives.

The unitary dynamics generated by the Hamiltonian in Eq. (35) can be expressed
through a path integral on the Schwinger-Keldysh contour Kamenev (2011), whose
generating functional reads

$\displaystyle\mathcal{Z}[J]$
$\displaystyle=\int\mathcal{D}\boldsymbol{\sigma}^{+}\mathcal{D}\boldsymbol{\sigma}^{-}\exp\left[i\left(S[\boldsymbol{\sigma}^{+}]-S[\boldsymbol{\sigma}^{-}]\right)/\hbar\right]\braket{\boldsymbol{\sigma}^{+}}{\hat{\rho}(0)}{\boldsymbol{\sigma}^{-}}\
.$ (36)

Here, $\hat{\rho}(0)$ represents the element of the initial density matrix at
$t=0$ and is chosen to be a random infinite-temperature initial state, which
is uncorrelated with the disorder. The last factor in the path integral is
then just an irrelevant factor. The action $S$ can be defined in terms of a
quadratic term $S_{0}$ and a disordered interaction term
$V_{J}(\boldsymbol{\sigma})$ as follows:

$\displaystyle S[\boldsymbol{\sigma},J]$
$\displaystyle=S_{0}[\boldsymbol{\sigma}]-\int_{0}^{\infty}dt\;\phantom{.}\mathcal{J}(t)V_{J}\big{(}\boldsymbol{\sigma}(t)\big{)},$
(37) $\displaystyle S_{0}[\boldsymbol{\sigma}]$
$\displaystyle=\int_{0}^{\infty}dt\;\left[\frac{M}{2}\dot{\boldsymbol{\sigma}}^{2}-\frac{z(t)}{2}(\boldsymbol{\sigma}^{2}-N)\right]-\sum_{i=1}^{N}\int_{0}^{\infty}dt\;B(t)\sigma_{i}(t),$
(38) $\displaystyle V_{J}(\boldsymbol{\sigma})$
$\displaystyle=-\sum_{i_{1}<...<i_{p}}^{N}J_{i_{1}...i_{p}}\sigma_{i_{1}}...\sigma_{i_{p}}.$
(39)

The quadratic part $S_{0}$ contains a kinetic term, chosen such that the
eventual dynamical equations are written in terms of second derivatives with
respect to time. The time-dependent Lagrange multiplier $z(t)$ enforces the
constraint $\sum_{i}^{N}\sigma_{i}^{2}=N$. We also included in $S_{0}$ the
coupling to the time-dependent field $B(t)$. This action can be split into
components $\boldsymbol{\sigma}^{+}$ and $\boldsymbol{\sigma}^{-}$ residing on
the forward and backwards Keldysh contours respectively, and recast as

$\displaystyle S[\boldsymbol{\sigma}^{+},\boldsymbol{\sigma}^{-},J]$
$\displaystyle=S_{0}[\boldsymbol{\sigma}^{+}]-S_{0}[\boldsymbol{\sigma}^{-}]-\int_{0}^{\infty}dt\;\mathcal{J}(t)\left[V_{J}\big{(}\boldsymbol{\sigma}^{+}(t)\big{)}-V_{J}\big{(}\boldsymbol{\sigma}^{-}(t)\big{)}\right]$
(40)

where the relative minus sign comes from reversing the integration limits on
the reverse contour.

###  E.1 System-Bath Coupling

The coupling between the system and bath can be treated exactly as in Ref.
Cugliandolo and Lozano (1999). We couple the quantum $p$-spin Hamiltonian
linearly to a bath of harmonic oscilllators, assumed to be in thermal
equilibrium, for a time window $[0,t_{b}]$. This coupling can be described by
the functional Kamenev (2011):

$\displaystyle
S_{bath}=\int_{0}^{t_{b}}dt\int_{0}^{t_{b}}dt^{\prime}\;\big{\\{}-[\boldsymbol{\sigma}^{+}(t)-\boldsymbol{\sigma}^{-}(t)]\eta(t-t^{\prime})[\boldsymbol{\sigma}^{+}(t^{\prime})+\boldsymbol{\sigma}^{-}(t^{\prime})]$
(41)
$\displaystyle+i~{}[\boldsymbol{\sigma}^{+}(t)-\boldsymbol{\sigma}^{-}(t)]\nu(t-t^{\prime})[\boldsymbol{\sigma}^{+}(t^{\prime})-\boldsymbol{\sigma}^{-}(t^{\prime})]~{}\big{\\}}$

where $\eta$ and $\nu$ are the correlation and response functions of the bath,
and are time-translation invariant due to the bath being in equilibrium. They
are given by

$\displaystyle\eta(t-t^{\prime})$
$\displaystyle=-\Theta(t-t^{\prime})\int_{0}^{\infty}d\omega
I(\omega)\sin[\omega(t-t^{\prime})]\;,$ (42) $\displaystyle\nu(t-t^{\prime})$
$\displaystyle=\int_{0}^{\infty}I(\omega)\coth\left(\frac{1}{2}\beta\hbar\omega\right)\cos[\omega(t-t^{\prime})]\;,$
(43)

where $I(\omega)$ is the spectral function of the bath. We choose an Ohmic
bath with $I(\omega)=\frac{1}{\pi}\exp\left(-|\omega|/\Lambda\right)$ and set
the integration cutoff to be $\Lambda=5$.

###  E.2 Disorder Averaging

For an initial condition $\hat{\rho}(0)$ uncorrelated with the disorder, we
can perform the disorder average explicitly. As in the Letter, we assume that
the distribution of the disorder variable $J_{i_{1}...i_{p}}$ is given by the
Gaussian distribution

$\displaystyle P[J]$ $\displaystyle=\sqrt{\frac{N^{p-1}}{\pi
p!J^{2}}}\exp\left(-\frac{N^{p-1}}{p!J^{2}}\sum_{i_{1}\neq...\neq
i_{p}}(J_{i_{1}...i_{p}})^{2}\right)$ (44)

with zero average and variance
$\overline{(J_{i_{1}...i_{p}})^{2}}=(p!J^{2})/(2N^{p-1})$. The disorder
average reduces to just averaging over the terms including the
$V_{J}(\boldsymbol{\sigma})$, which are the only ones depending on the
disorder. By an straightforward calculation, we obtain the disorder-averaged
generating functional:

$\displaystyle\overline{\mathcal{Z}[J]}=\int\mathcal{D}{\sigma}^{-}\mathcal{D}\sigma^{+}\exp\left[i\left(S_{eff}[\boldsymbol{\sigma^{+}},\boldsymbol{\sigma^{-}}]\right)/\hbar\right]\
.$ (45)

The effective action at the exponent is given by:

$\displaystyle S_{\rm eff}[\boldsymbol{\sigma^{+}},\boldsymbol{\sigma^{-}}]$
$\displaystyle=S_{0}[\boldsymbol{\sigma}^{+}]-S_{0}[\boldsymbol{\sigma}^{-}]-V_{D}[\boldsymbol{\sigma}^{+},\boldsymbol{\sigma}^{-}]+S_{bath}[\boldsymbol{\sigma^{+}},\boldsymbol{\sigma^{-}}],$
(46) $\displaystyle V_{D}[\boldsymbol{\sigma}^{+},\boldsymbol{\sigma}^{-}]$
$\displaystyle=\frac{iN}{4}\int
dtdt^{\prime}\;{\mathcal{J}(t)\mathcal{J}(t^{\prime})J^{2}}\sum_{\alpha,\beta=\pm}\sum_{i=1}^{N}\alpha\beta\left(\frac{1}{N}\sigma_{i}^{\alpha}(t)\sigma_{i}^{\beta}(t^{\prime})\right)^{p}$
(47)

where $\alpha,\beta=\pm$ are the Schwinger-Keldysh contours.

###  E.3 Transformed Order Parameters

The contribution to the action containing at most quadratic terms in the spins
$\sigma_{i}^{\alpha}(t)$ can be written down as:

$\displaystyle\frac{1}{\hbar}S_{\rm
eff}^{(2)}[\boldsymbol{\sigma}^{+},\boldsymbol{\sigma}^{-}]$
$\displaystyle=-\frac{1}{2}\sum_{i,\alpha,\beta}\int
dtdt^{\prime}\;\sigma_{i}^{\alpha}(t)O_{p}^{\alpha\beta}(t,t^{\prime})\sigma_{i}^{\beta}(t^{\prime})-\sum_{i,\alpha}\int\frac{dt}{\hbar}\
\alpha B(t)\sigma_{i}^{\alpha}(t)\ ,$ (48)

where again $\alpha,\beta=\pm$. The matrix elements
$O^{\alpha\beta}_{p}(t,t^{\prime})$, appearing in Eq. (48), are explicitly
given by

$\displaystyle O_{p}^{++}(t,t^{\prime})$
$\displaystyle=\frac{1}{\hbar}[M\partial_{t}^{2}+z^{+}(t)]\delta(t-t^{\prime})-\frac{2}{\hbar}(i\nu(t-t^{\prime})-\eta(t-t^{\prime}))\Theta(t-t_{b})\Theta(t^{\prime}-t_{b}),$
(49) $\displaystyle O_{p}^{+-}(t,t^{\prime})$
$\displaystyle=\frac{1}{\hbar}(2\eta(t-t^{\prime})+2i\nu(t-t^{\prime}))\Theta(t-t_{b})\Theta(t^{\prime}-t_{b}),$
(50) $\displaystyle O_{p}^{-+}(t,t^{\prime})$
$\displaystyle=\frac{1}{\hbar}(-2\eta(t-t^{\prime})+2i\nu(t-t^{\prime}))\Theta(t-t_{b})\Theta(t^{\prime}-t_{b}),$
(51) $\displaystyle O_{p}^{--}(t,t^{\prime})$
$\displaystyle=-\frac{1}{\hbar}[M\partial_{t}^{2}+z^{-}(t)]\delta(t-t^{\prime})-\frac{2}{\hbar}(i\nu(t-t^{\prime})+\eta(t-t^{\prime}))\Theta(t-t_{b})\Theta(t^{\prime}-t_{b}).$
(52)

###  E.4 Saddle-Point Equations

To decouple the $p$-body interaction term, we introduce a new set of variables
$Q^{\alpha\beta}(t,t^{\prime})$, by using the following trivial representation
of the number $1$:

$\displaystyle 1$
$\displaystyle=\int\prod_{\alpha\beta}\mathcal{D}Q^{\alpha\beta}\delta\left(\frac{1}{N}\boldsymbol{\sigma}^{\alpha}(t)\boldsymbol{\sigma}^{\beta}(t^{\prime})-Q^{\alpha\beta}(t,t^{\prime})\right),$
(53)
$\displaystyle\propto\int\prod_{\alpha\beta}\mathcal{D}Q^{\alpha\beta}\mathcal{D}\lambda^{\alpha\beta}\exp\left(-\frac{i}{2}\lambda^{\alpha\beta}\left(\boldsymbol{\sigma}^{\alpha}(t)\boldsymbol{\sigma}^{\beta}(t^{\prime})-NQ^{\alpha\beta}(t,t^{\prime})\right)\right).$
(54)

Then, using a compact notation, the final form of the generating functional is

$\overline{\mathcal{Z}[J]}=\int\mathcal{D}Q\mathcal{D}\lambda\;\exp\\{NS[Q,\lambda]\\}$
(55)

where

$S[Q,\lambda]=\sum_{\alpha\beta}\int
dtdt^{\prime}\Big{\\{}\frac{i}{2}\lambda^{\alpha\beta}(t,t^{\prime})Q^{\alpha\beta}(t,t^{\prime})-\frac{1}{4}Q_{\alpha\beta}(t,t^{\prime})\Big{\\}}+\log
Z[Q,\lambda]$ (56)

and

$\begin{split}Z[Q,\lambda]=\int\mathcal{D}{\sigma}^{-}\mathcal{D}\sigma^{+}\exp\Big{\\{}&-\frac{1}{2}\sum_{\alpha\beta}\int
dtdt^{\prime}\sigma^{\alpha}(t)[iO_{p}^{\alpha\beta}(t,t^{\prime})+i\lambda^{\alpha\beta}(t,t^{\prime})]\sigma^{\beta}(t^{\prime})+\\\
&-i\sum_{\alpha}\int\frac{dt}{\hbar}\alpha
B(t)\sigma_{i}^{\alpha}(t)\Big{\\}}\end{split}$ (57)

is the action of a single effective spin. Then, defining the matrix
$M^{\alpha\beta}(t,t^{\prime})=iO_{p}^{\alpha\beta}(t,t^{\prime})+i\lambda^{\alpha\beta}(t,t^{\prime})$,
the saddle point equations for the $N\to\infty$ limit read

$\displaystyle
Q^{\alpha\beta}(t,t^{\prime})=\braket{\sigma^{\alpha}(t)\sigma^{\beta}(t^{\prime})}=(M^{-1})^{\alpha\beta}(t,t^{\prime})+m(t)m(t^{\prime})$
(58) $\displaystyle
i\lambda^{\alpha\beta}(t,t^{\prime})=\frac{p}{2}F[Q]^{\alpha\beta}(t,t^{\prime})$
(59)

The average $\braket{\dots}$ in Eq. (58) is performed over the partition
function $Z[Q,\lambda]$ and the magnetization,

$m(t)=\braket{\sigma(t)}=\sum_{\beta}\int
dt^{\prime}\;\frac{i}{\hbar}(M^{-1})^{\alpha\beta}(t,t^{\prime})\beta
B(t^{\prime})\ ,$ (60)

does not actually depend on the Keldysh index $\alpha$, like every other one-
time quantity. The matrices used in Eq. (58) are defined as

$\displaystyle\mathcal{Q}(t,t^{\prime})$
$\displaystyle=\begin{bmatrix}Q^{++}(t,t^{\prime})&Q^{+-}(t,t^{\prime})\\\
Q^{-+}(t,t^{\prime})&Q^{--}(t,t^{\prime})\end{bmatrix},\quad\mathcal{M}(t,t^{\prime})=\begin{bmatrix}M^{++}(t,t^{\prime})&M^{+-}(t,t^{\prime})\\\
M^{-+}(t,t^{\prime})&M^{--}(t,t^{\prime})\end{bmatrix}$ (61) $\displaystyle
F[Q](t,t^{\prime})$
$\displaystyle=\begin{bmatrix}[Q^{++}(t,t^{\prime})]^{p-1}&-[Q^{+-}(t,t^{\prime})]^{p-1}\\\
-[Q^{-+}(t,t^{\prime})]^{p-1}&[Q^{--}(t,t^{\prime})]^{p-1}\end{bmatrix}\frac{\mathcal{J}(t)\mathcal{J}(t^{\prime})J^{2}}{\hbar^{2}}$
(62)

We used the same notation of Ref. Cugliandolo and Lozano (1999). Manipulating
Eq. (58), we can rewrite the response function as

$R(t,t^{\prime})=\frac{i}{\hbar}[Q^{++}(t,t^{\prime})-Q^{+-}(t,t^{\prime})]=\frac{i}{\hbar}[(M^{-1})^{++}(t,t^{\prime})-(M^{-1})^{+-}(t,t^{\prime})]\
.$ (63)

Then, by fixing $\alpha=+1$ in Eq. (60), the magnetization becomes

$m(t)=\int dt^{\prime}\;R(t,t^{\prime})B(t^{\prime})\ .$ (64)

Using Eq. (59) and applying the matrix $\mathcal{M}$ on both sides of Eq. Eq.
(58), we obtain the dynamical equations for all the two-point correlators on
the Keldysh contour, represented in the following compact form:

$\displaystyle
iO_{p}\otimes\mathcal{Q}(t,t^{\prime})=\mathcal{I}-\frac{p}{2}F[Q]\otimes\mathcal{Q}(t,t^{\prime})+\mathcal{S}\
B(t)\braket{\sigma(t^{\prime})}$ (65)

where the matrix elements $\mathcal{S}^{\alpha\beta}=\alpha$ do not depend on
the time indices.

###  E.5 Dynamical Equations

Following the same prescription described in Ref. Cugliandolo and Lozano
(1999), it is straightforward to obtain a closed set of dynamical equations
for the response function $R(t,t^{\prime})$ and the (symmetric) correlation
function

$\displaystyle
C(t,t^{\prime})=\frac{1}{2}[Q^{+-}(t,t^{\prime})+Q^{-+}(t,t^{\prime})]\ .$
(66)

In particular, the equations of motion for the response function
$R(t,t^{\prime})$ are obtained by taking the difference of the $++$ and $+-$
components of Eq. (65), while the ones for $C(t,t^{\prime})$ are obtained by
taking the addition of the $+-$ and $-+$ components. After some algebra, the
result is given by:

$\displaystyle[M\partial_{t}^{2}+z(t)]R(t,t^{\prime})$
$\displaystyle=\delta(t-t^{\prime})+\int_{0}^{t}dt^{\prime\prime}\;\Sigma(t,t^{\prime\prime})R(t^{\prime\prime},t^{\prime})\
,$ (67) $\displaystyle[M\partial_{t}^{2}+z(t)]C(t,t^{\prime})$
$\displaystyle=\int_{0}^{t}dt^{\prime\prime}\;\Sigma(t,t^{\prime\prime})C(t^{\prime\prime},t^{\prime})+\int_{0}^{t^{\prime}}dt^{\prime\prime}\;D(t,t^{\prime\prime})R(t^{\prime},t^{\prime\prime})+\
,$ (68)
$\displaystyle+B(t)\int_{0}^{t^{\prime}}dt^{\prime\prime}\;R(t^{\prime},t^{\prime\prime})B(t^{\prime\prime})\
.$

Here, we have defined the self-energy $\Sigma(t,t^{\prime})$ and the vertex
$D(t,t^{\prime})$ as

$\displaystyle\Sigma(t,t^{\prime})=$
$\displaystyle-4\eta(t-t^{\prime})\Theta(t_{b}-t)\Theta(t_{b}-t^{\prime})\ +$
(69)
$\displaystyle-\frac{p\mathcal{J}(t)\mathcal{J}(t^{\prime})J^{2}}{\hbar}\textrm{Im}\left[C(t,t^{\prime})-\frac{i\hbar}{2}R(t,t^{\prime})\right]^{p-1},$
$\displaystyle D(t,t^{\prime})=$ $\displaystyle\
2\hbar\nu(t-t^{\prime})\Theta(t_{b}-t)\Theta(t_{b}-t^{\prime})\ +$ (70)
$\displaystyle+\frac{p\mathcal{J}(t)\mathcal{J}(t^{\prime})J^{2}}{2}\textrm{Re}\left[C(t,t^{\prime})-\frac{i}{2}(\hbar
R(t,t^{\prime})+\hbar R(t^{\prime},t))\right]^{p-1}.$

This result is exactly the one that we anticipated at the beginning of Section
3, now with explicit functional forms for the two kernels.

The equations for the classical PSM can be readily obtained from the ones
above by simply taking the $\hbar\to 0$ limit. The global structure of the
equations is the same, only the $\Sigma$ and $D$ kernels are affected and
considerably simplified.

##  Appendix F 3.2 The evolution equation for the Lagrange multiplier

In this section we discuss how to self-consistently determine the Lagrange
multiplier $z(t)$. The approach typically used in the literature Cugliandolo
and Lozano (1999) consists in evaluating Eqs. (D) and (31) at equal times
$t^{\prime}=t$:

$z(t)=\int_{0}^{t}dt^{\prime\prime}\;\Big{[}\Sigma(t,t^{\prime\prime})C(t^{\prime\prime},t)+D(t,t^{\prime\prime})R(t,t^{\prime\prime})\Big{]}-M\partial_{t}^{2}C(t,t^{\prime})\big{|}_{t^{\prime}\to
t}\ .$ (71)

We observe that Eq. (71) leads to ambiguities, as it tautologically depends on
the second derivative of $C(t,t^{\prime})$ at equal times. This issue is
irrelevant while the system is coupled to a thermal bath Cugliandolo and
Lozano (1999), which naturally leads $z(t)$ to thermalize, or when the total
energy is conserved, as the second derivative can be replaced by a causal
expression containing the conserved energy density Thomson  _et al._ (2020);
Cugliandolo  _et al._ (2017). However, for $t>t_{b}$ none of these two
condition is met in our protocol, as the system evolves under a unitary
dynamics generated by a time-dependent Hamiltonian. In this case, Eq. (71)
does not determine $z(t)$. Here we solve this issue by determining $z(t)$ from
the spherical constraint, according to the following procedure. We take the
total derivative of the constraint equation $C(t,t)=1$ multiple times and
obtain:

$\begin{split}0&=\frac{d}{dt}C(t,t)\ =\lim_{t^{\prime}\to
t}\partial_{t}C(t,t^{\prime})\;,\\\
0&=\frac{d^{2}}{dt^{2}}C(t,t)=\lim_{t^{\prime}\to t}\
[\partial^{2}_{t}C(t,t^{\prime})+\partial_{t}\partial_{t^{\prime}}C(t,t^{\prime})]\;,\\\
0&=\frac{d^{3}}{dt^{3}}C(t,t)=\lim_{t^{\prime}\to t}\
[\partial^{3}_{t}C(t,t^{\prime})+3\partial^{2}_{t}\partial_{t^{\prime}}C(t,t^{\prime})]\;.\\\
\end{split}$ (72)

We used a compact notation $\partial_{t}=d/dt$ and we also used the symmetry
relation $C(t,t^{\prime})=C(t^{\prime},t)$ in the last equality of each line
of Eqs. (72). It is easy to realize that the first line of Eqs. (72)
corresponds to some of the equal-time conditions described in Eqs. (33) and
that the second line is useless, as
$\partial_{t}\partial_{t^{\prime}}C(t,t^{\prime})$ cannot be computed using
the dynamical equations. The third line is instead the one that we use to
determine $z(t)$. Specifically, we observe that the right-hand-side of its
second equality can be rewritten in terms of partial derivatives, performed
with respect to $t$ and $t^{\prime}$ respectively, of both sides of Eq. (31).
After some algebra, we obtain

$\displaystyle\frac{dz}{dt}$
$\displaystyle=\int_{0}^{t}dt^{\prime\prime}\;\Big{[}\partial_{t}\Sigma(t,t^{\prime\prime})C(t,t^{\prime\prime})+\partial_{t}D(t,t^{\prime\prime})R(t,t^{\prime\prime})+3\Sigma(t,t^{\prime\prime})\partial_{t}C(t,t^{\prime\prime})+3D(t,t^{\prime\prime})\partial_{t}R(t,t^{\prime\prime})\Big{]}\
+$
$\displaystyle+\partial_{t}B(t)\int_{0}^{t}dt^{\prime\prime}\;R(t,t^{\prime\prime})B(t^{\prime\prime})+3B(t)\int_{0}^{t}dt^{\prime\prime}\;\partial_{t}R(t,t^{\prime\prime})B(t^{\prime\prime})$
(73)

From Eq. (F) it is easy to determine the Lagrange multiplier in a causal form,
as

$\displaystyle z(t)$
$\displaystyle=\int_{t_{0}}^{t}dt^{\prime}\int_{0}^{t^{\prime}}dt^{\prime\prime}\;\Big{[}\partial_{t^{\prime}}\Sigma(t^{\prime},t^{\prime\prime})C(t,t^{\prime\prime})+\partial_{t^{\prime}}D(t^{\prime},t^{\prime\prime})R(t^{\prime},t^{\prime\prime})+3\Sigma(t^{\prime},t^{\prime\prime})\partial_{t^{\prime}}C(t^{\prime},t^{\prime\prime})+3D(t^{\prime},t^{\prime\prime})\partial_{t^{\prime}}R(t^{\prime},t^{\prime\prime})\Big{]}+$
$\displaystyle\
+\int_{t_{0}}^{t}dt^{\prime}\int_{0}^{t^{\prime}}dt^{\prime\prime}\;\Big{[}\partial_{t^{\prime}}B(t^{\prime})R(t,t^{\prime\prime})B(t^{\prime\prime})+3B(s)\partial_{t^{\prime}}R(t^{\prime},t^{\prime\prime})B(t^{\prime\prime})\Big{]}+z(t_{0})$
(74)

provided that we can access the value $z(t_{0})$, for some $t_{0}$. Thus we
first solve numerically the dissipative dynamics generated from Eqs. (D) and
(31), for $0<t<t_{b}$, using the standard expression from Eq. (71).
Subsequently, we solve the “closed dynamics” for $t>t_{b}$ using the
expression from Eq. (F) for $t_{0}=t_{b}$ and by using the knowledge of
$z(t_{b})$ that we obtained by integrating the dissipative dynamics.

##  Appendix G 3.3 Predictor-corrector scheme for the Mode-Coupling equations

To solve Eqs. (D) and (31) in the Letter, we first introduce a discrete time
step $\Delta t$ and we discretize the times as $t=n\Delta t$,
$t^{\prime}=m\Delta t$, etc. for $n$ and $m$ non-negative integers. We use a
“forward” discretization scheme for the time derivatives, that is

$\partial_{t}f(t,t^{\prime})\approx(f_{n+1,m}-f_{n,m})/\Delta t\ ,$ (75)

for every two-point function $f(t,t^{\prime})$. With this notation, the
discretized version of Eqs. (D) and (31) equations is

$\left\\{\begin{split}R_{n+1,m}&=R_{n,m}+\Pi^{(R)}_{n,m}\Delta t/M\;,\\\
C_{n+1,m}&=C_{n,m}+\Pi^{(C)}_{n,m}\Delta t/M\;,\\\
\Pi^{(R)}_{n+1,m}&=\Pi^{(R)}_{n,m}-z_{n}\Delta tR_{n,m}+\delta_{n,m}+\Delta
tF^{(1)}_{n,m}\;,\\\ \Pi^{(C)}_{n+1,m}&=\Pi^{(C)}_{n,m}-z_{n}\Delta
tC_{n,m}+\Delta tF^{(2)}_{n,m}\;,\end{split}\right.$

where

$\left\\{\begin{split}F^{(1)}_{n,m}&\ =\ \Delta
t\sum_{j=m}^{n}\Sigma_{nj}R_{jm}\;,\\\ F^{(2)}_{n,m}&\ =\ \Delta t\
\Big{[}\sum_{j=0}^{m}D_{nj}R_{jm}+\sum_{j=0}^{n}\Sigma_{nj}C_{jm}\Big{]}\;,\\\
\Sigma_{n,m}&\ =\ -\frac{p\mathcal{J}(n\Delta t)\mathcal{J}(m\Delta
t)J^{2}}{\hbar}\
\text{Im}\Big{[}C_{n,m}-\frac{i\hbar}{2}R_{n,m}\Big{]}^{p-1}~{}-\
4\eta_{n,m}\Theta(t_{b}-n\Delta t)\Theta(t_{b}-m\Delta t)\;,\\\ D_{n,m}&\ =\
\frac{p\mathcal{J}(n\Delta t)\mathcal{J}(m\Delta t)J^{2}}{\hbar}\
\text{Re}\Big{[}C_{n,m}-\frac{i\hbar}{2}\big{(}R_{n,m}+R_{m,n})\Big{]}^{p-1}~{}+\
2\hbar\nu_{n,m}\Theta(t_{b}-n\Delta t)\Theta(t_{b}-m\Delta
t)\;.\end{split}\right.$ (76)

The equal time conditions, from Eqs. (33), here assume the simple form

$C_{n,n}=0\ ,\qquad\qquad\quad R_{n,n}=\Pi^{(R)}_{n,n}=\Pi^{(C)}_{n,n}=0\ .$
(77)

It is easy to note that Eqs. (G) are causal, implying that we can iteratively
compute each minor matrix $C_{0\leq i\leq n,0\leq j\leq n}$ from the knowledge
of $C_{0\leq i\leq n-1,0\leq j\leq n-1}$. To understand in detail why, we
notice that the row $C_{n,0\leq j\leq n-1}$ can be straightforwardly computed
from the minor $C_{0\leq i\leq n-1,0\leq j\leq n-1}$ using Eqs. (G), then the
column $C_{0\leq i\leq n-1,n}$ is immediately obtained from the symmetry
relation $C_{n,m}=C_{m,n}$ and the diagonal element $C_{n,n}=1$ is given by
the spherical constraint. The same reasoning holds for the response function,
with the only difference that $R_{n,m}$ is not symmetric and $R_{0\leq i\leq
n-1,n}=0$ due to causality.

We aim to integrate Eqs. (D) and (31) up to some finite time
$t_{max}=n_{steps}\Delta t$. With the discretization we chose, the error we
make is of order $n_{steps}\Delta t^{2}$. In order to increase the precision
of our results, we improve our method by employing a predictor-corrector
algorithm, already used in a similar context in Ref. Haldar  _et al._ (2020).
In a nutshell, we first predict the $(n+1)$-th row of the correlation and
response function using Eq. (G), and we then correct our result by inserting
the result of the prediction in the right-hand sides of the following
equations of motion:

$\left\\{\begin{split}\Pi^{(R)}_{n+1,m}&=\Pi^{(R)}_{n,m}-\Delta
t\frac{z_{n}R_{n,m}+z_{n+1}R_{n+1,m}}{2}+\delta_{n,m}+\Delta
t\frac{F^{(1)}_{n,m}+F^{(1)}_{n+1,m}}{2}\;,\\\
\Pi^{(C)}_{n+1,m}&=\Pi^{(C)}_{n,m}-\Delta
t\frac{z_{n}C_{n,m}+z_{n+1}C_{n+1,m}}{2}+\delta_{n,m}+\Delta
t\frac{F^{(2)}_{n,m}+F^{(2)}_{n+1,m}}{2}\;,\\\ R_{n+1,m}&=R_{n,m}+\Delta
t\frac{\Pi^{(R)}_{n,m}+\Pi^{(R)}_{n+1,m}}{2m}\;,\\\ C_{n+1,m}&=C_{n,m}+\Delta
t\frac{\Pi^{(C)}_{n,m}+\Pi^{(C)}_{n+1,m}}{2m}\;.\end{split}\right.$ (78)

For each $n$-th step, we take a loop over the predictor-corrector scheme
$\mathcal{N}_{L}$ times. In this way, the error we make is of order
$n_{steps}\Delta t^{(1+\mathcal{N}_{L})}$. For the results presented in the
Letter, we have always fixed $\mathcal{N}_{L}=2$.

At this point, we observe that Eqs. (78) alone still do not determine the
Lagrange multiplier $z_{n}$, which is in principle determined by the non-
causal equation (equivalent to Eq. (71) of the Letter):

$z_{n}=-\frac{\Pi^{(C)}_{n+1,n}-\Pi^{(C)}_{n,n}}{\Delta t}+\Delta
tF^{(2)}_{n,m}\ .$ (79)

As long as the system is coupled to a thermal bath, we solve this issue by
making the physical assumption that, due to dissipation, $z_{n}$ converges at
large times to a stationary value, as also observed in Ref. Cugliandolo and
Lozano (1999). Due to this asymptotic convergence, we can safely replace the
difference $\Pi^{(C)}_{n+1,n}-\Pi^{(C)}_{n,n}$ with the one evaluated at the
previous time step, $\Pi^{(C)}_{n,n}-\Pi^{(C)}_{n,n-1}$, in Eq. (79). However,
this substitution is not valid for $t>t_{b}$, where the dynamics is isolated
and periodically driven. In the latter scenario, we proceed as discussed in
Section F and include the discretized version of Eq. (F) in the system of Eqs.
(G), so that also $z_{n}$ can be computed using the predictor-corrector
algorithm.

Figure 4: (Color online). (a-c) Dynamics of the time-independent 1D lattice
$\phi^{4}$ model. We fix $N=100$ and $\omega_{0}^{2}=g=\lambda=1$. (a) Initial
energy density as a function of the inverse temperature, from Eq. (82). (b-c)
Kurtosis of the position and momentum degrees of freedom, from Eq. (83). Data
are averaged over $\mathcal{N}=2000$ configurations. (d-f) Dynamics of the
time-independent $p$-spin spherical model. (d) Equilibrium dependence of
temperature from the energy density $e$, from Eq. (84). (e-f) Kinetic energy
density from Eq. (85), for the classical ($\Gamma=0$) and quantum ($\Gamma=1$)
cases. Dashed lines corresponds to the corresponding classical equilibrium
expectation values from Eq. (86). We fix $M=J=1$ and $p=3$.

##  Appendix H 4\. Ergodic properties of lattice $\phi^{4}$ model and the
$p$-spin spherical model in absence of a drive

In this section, we discuss the ergodic properties of both the 1D lattice
$\phi^{4}$ model and of the $p$-spin spherical model, in absence of a periodic
drive.

We begin by focusing on the $\phi^{4}$ model, whose dynamics is defined by Eq.
(19) at $B_{0}=0$. We initialize the system in a Gaussian Gibbs distribution,
corresponding to the $\phi^{4}$ Hamiltonian at $\lambda=0$ and given by

$\rho_{0}(\mathbf{x},\mathbf{p})=\frac{1}{\mathcal{Z}_{0}}\exp\Big{\\{}-\frac{\beta}{2}\sum_{k=0}^{N-1}(\tilde{p}_{k}\tilde{p}_{-k}+\omega_{k}^{2}\tilde{x}_{k}\tilde{x}_{-k})\Big{\\}}$
(80)

Here, the variables

$\tilde{x}_{k}=\frac{1}{\sqrt{N}}\sum_{j}e^{i\frac{2\pi
k}{N}j}x_{j}~{},\qquad\tilde{p}_{k}=\frac{1}{\sqrt{N}}\sum_{j}e^{i\frac{2\pi
k}{N}j}p_{j}~{},$ (81)

are the Fourier transform of the phase space variables, $\mathcal{Z}_{0}$ is a
normalization factor and $\omega_{k}=\sqrt{\omega_{0}^{2}+2g[1-\cos(2\pi
k/N)]}$ are the normal mode frequencies of the chain. We study the non-
equilibrium dynamics evolving from $\rho_{0}(\mathbf{x},\mathbf{p})$, realized
through a quench in the coupling $\lambda$. It is straightforward to show that
the initial energy density of the system is given by

$e_{0}(\beta)=\frac{1}{N}\braket{H_{\phi^{4}}(\mathbf{x},\mathbf{p})}=\frac{1}{\beta}+\frac{3\lambda}{4\beta^{2}}\left(\frac{1}{N}\sum_{k=0}^{N-1}\frac{1}{\omega_{k}^{2}}\right)^{2}\
.$ (82)

The two terms on the right-hand side of Eq. (82) correspond to the averages of
the quadratic and the quartic part of the Hamiltonian, respectively. We plot
the function $e_{0}(\beta)$ in Fig. 4-(a), for a range of $\beta$ roughly
corresponding to values of $e_{0}$ used in the Letter. We integrate the
dynamics for some values of $\beta$ from the same range and study the time-
evolution of the Kurtosis parameters of the phase space degrees of freedom,
defined as:

$K_{x}(t)=\frac{\braket{\sum_{i}x_{i}(t)^{4}}}{3\braket{\sum_{i}x_{i}(t)^{2}}}-1\
,\qquad
K_{p}(t)=\frac{\braket{\sum_{i}p_{i}(t)^{4}}}{3\braket{\sum_{i}p_{i}(t)^{2}}}-1\
.$ (83)

Both $K_{x}(t)$ and $K_{p}(t)$ vanish for a Gaussian distribution Joanes and
Gill (1998), like the one in Eq. (80), so that $K_{x}(0)=K_{p}(0)=0$. We plot
both $K_{x}(t)$ and $K_{p}(t)$ in Fig. 4-(b), for several values of $\beta$.
We observe that $K_{x}(t)$ and $K_{p}(t)$ approach an asymptotic value at late
times, although the dynamics moderately slows down at low values of $\beta$.
In particular, $K_{p}(t)$ asymptotically vanishes, while $K_{x}(t)$ does not.
This result is compatible with the approach to a thermal state of the
$\phi^{4}$ model with $\lambda>0$, which is Gaussian in the momentum
variables, but not in the position ones.

Next, we explore the isolated dynamics of the classical and quantum PSM,
defined in Eqs. (D) and (31). We fix $B(t)=0$ and $\mathcal{J}(t)=1$, to
eliminate the periodic drive, and set $t_{b}=0$ in Eqs. (69) and (70), to
isolate the system from the external bath. We initialize the the Lagrange
multiplier $z(t)$, from Eq. (71), to $z(0)=2e$. This choice corresponds to a
non-equilibrium initial state, uncorrelated with the disorder and
characterized by a vanishing potential energy and kinetic energy density equal
to $e$ Cugliandolo  _et al._ (2017). As in the classical case,
$\Gamma\equiv\sqrt{\hbar^{2}/M}=0$, the equilibrium energy density $e$ is
related to the temperature by Cugliandolo  _et al._ (2017)

$T_{cl}(e)=e+\sqrt{e^{2}+J^{2}}\ ,$ (84)

we integrate the dynamics for a set of $e$ roughly corresponding to the
initial temperatures $T_{0}$ used in the Letter, as it can be inferred from
Fig. 4-(d). In Fig. 4-(e), we plot the evolution of the kinetic energy
density,

$e_{k}(t)=\frac{1}{N}\sum_{i=1}^{N}\braket{\dot{\sigma}_{i}^{2}(t)}\ ,$ (85)

for the classical case $\Gamma=0$. We observe that $e_{k}(t)$ quickly
saturates to a late time plateau, which is compatible with the equilibrium
expectation value (dashed lines)

$e_{k}^{(mc)}=\frac{e+\sqrt{e^{2}+J^{2}}}{2}.$ (86)

Equation (86) obtained from Eq. (84) and equipartition theorem. This result
proves the fast ergodicity of the system, for the range of temperatures
considered in the Letter. A similar result is obtained in the quantum case
$\Gamma=1$, as shown in Fig. 4-(f), where the kinetic energy density relaxes
roughly to the same plateaus. The similarity between the classical and quantum
case suggest that quantum effects are negligible for the parameter we
considered.

##  Appendix I 5\. Floquet dynamics with different driving terms

In this section, we investigate the periodically driven dynamics of the
classical lattice $\phi^{4}$ model and the $p$-spin spherical model (PSM). We
demonstrate that the energy density profiles highlighted in the Letter remain
qualitatively unchanged even when altering the form of the periodic drives.

Figure 5: Average energy density
$e(t)=\braket{H_{\phi^{4}}(\textbf{x}(t),\textbf{p}(t))}_{0}/N$, for the one
dimensional lattice $\phi^{4}$ model, evolving from Eqs. (87). We fix $N=100$,
$\Omega=2.3$ and $a=\omega_{0}^{2}=g=1$. The average displayed here is
performed over $\mathcal{N}$ configurations, sampled from a microcanonical
manifold at various energy density $e_{0}$, for the values of $e_{0}$ listed
in the legends. We fix $\mathcal{N}=1000$ and $\lambda=1$ in panels (a) and
(b), while $\mathcal{N}=500$ and $\lambda=0.1$ in panel (c). Each panel
corresponds to a different form of drive $F_{i}(\mathbf{x})$. (a)
$F_{i}(\mathbf{x})=\cos(x_{i})$. (b) $F_{i}(\mathbf{x})=x_{i}^{2}$. (c)
$F_{i}(\mathbf{x})=x_{i}$.

We begin with the lattice $\phi^{4}$ model, the Hamilton equations of motion
of which are given by

$\ddot{x}_{i}+\omega_{0}^{2}x_{i}+\lambda x_{i}^{3}=a\sin(\Omega
t)F_{i}(\mathbf{x})+g(x_{i+1}-2x_{i}+x_{i-1})\ ,$ (87)

where $i=1\ldots N$. If we set $F_{i}(\mathbf{x})=1$, we recover the same
Floquet dynamics examined in the Letter. Instead, here we focus on driving
terms corresponding to one of the following possibilities:
$F_{i}(\mathbf{x})=x_{i}^{2}$ or $F_{i}(\mathbf{x})=\cos(x_{i})$. We
numerically integrate Eqs. (87) over an ensemble of $\mathcal{N}$ initial
configurations $\\{\mathbf{x}(0),\dot{\mathbf{x}}(0)\\}$, randomly sampled on
the manifold
$H_{\phi^{4}}\big{(}\mathbf{x}(0),\dot{\mathbf{x}}(0)\big{)}=Ne_{0}$, for
several values of $e_{0}$. The resulting average energy density,
$e(t)=\braket{H_{\phi^{4}}(\textbf{x}(t),\dot{\textbf{x}}(t))}_{0}/N$, is
depicted in Fig. 5-(a) and (b). For both driving terms, we consistently
observe that $e(t)$ saturates at long times to a finite value $e_{\infty}$,
underscoring the robustness of the results presented in the Letter across
different driving forms. In Fig. 5-(c) we examine a third, more subtle case
involving a periodic modulation of $\omega_{0}^{2}$, achieved by setting
$F_{i}(\mathbf{x})=x_{i}$. The key observation here is that, for a low initial
energy $e_{0}$, the average energy density $e(t)$ exhibits slow but persistent
growth within our simulation time scales. This profile, also observed in the
closely related $O(N)$ model under the same driving protocol Weidinger and
Knap (2017), does not provide a definitive indication of either unbounded
growth or saturation to a finite value $e_{\infty}$. However, we note that the
heating diminishes with an increase in the initial energy and eventually
disappears for a sufficiently high $e_{0}$. This observation may indicate the
existence of an energy threshold at which heating is expected to cease,
leading to the saturation of $e(t)$ to a finite value even when evolving from
a low initial energy. To validate this conclusion, further investigation is
required.

Figure 6: Average energy density $e(t)$ for the $p$-spin spherical model,
defined in Eq. (34) and evolved using Eq. (D) and Eq. (31). We drive the
system by setting $B(t)=0$ and $\mathcal{J}(t)=1+\Delta\sin(\Omega t)$, for
$\Delta=0.5$ and $\Omega=5$. We also fix $p=3$, $M=J=1$, $t_{b}=50$. Each
panel corresponds to a different value of the dimensionless parameter
$\Gamma=\sqrt{\hbar^{2}/M}$ which quantifies the strength of quantum
fluctuations. The different curves in each panel correspond to different
initial conditions characterized by their temperature $T_{0}$ specified in the
keys.

We conduct a similar analysis for the PSM, whose dynamics is governed by the
equations of motion (D) and (31). In the Letter, we investigated Floquet
dynamics driven by a homogeneous field, achieved by setting
$B(t)=B_{0}\sin(\Omega t)$ and $\mathcal{J}(t)=1$ in the Hamiltonian in Eq.
(35). Here, our focus shifts to a modulation in the interaction term, achieved
by setting $B(t)=0$ and $\mathcal{J}(t)=1+\Delta\sin(\Omega t)$. As
illustrated in Fig. 6, the corresponding energy density, Eq. (34), saturates
at long times to a finite value, consistently with the results presented in
the Letter for a different kind of drive.