Prethermalization in aperiodically kicked many-body dynamics

Driven many-body systems typically experience heating due to the lack of energy conservation. Heating may be suppressed for time-periodic drives, but little is known for less regular drive protocols. In this work, we investigate the heating dynamics in aperiodically kicked systems, specifically those driven by quasi-periodic Thue-Morse or a family of random sequences with $n$-multipolar temporal correlations. We demonstrate that multiple heating channels can be eliminated even away from the high-frequency regime. The number of eliminated channels increases with multipolar order $n$. We illustrate this in a classical kicked rotor chain where we find a long-lived prethermal regime. When the static Hamiltonian only involves the kinetic energy, the prethermal lifetime $t^*$ can strongly depend on the temporal correlations of the drive, with a power-law dependence on the kick strength $t^*\sim K^{-2n}$, for which we can account using a simple linearization argument.

Introduction.-Time-dependentmany-body systems have attracted sustained interest due to their ubiquity in nature and the potential to realize novel non-equilibrium phases of matter.One typical example is the discrete time crystal which spontaneously breaks discrete time translation symmetry (TTS) [1,2].However, due to the absence of energy conservation, closed driven systems tend to heat up and lose any non-trivial correlations [3,4].Therefore, understanding and controlling the onset of heating in time-dependent systems is key to their stabilization, and to the realization of exotic nonequilibrium phenomena.
It is natural to ask: can heating be efficiently suppressed in many-body systems without TTS, e.g., when drives are quasi-periodic, or even random?This is a notoriously difficult question as breaking TTS generally opens up further deleterious heating channels that can destabilize systems rapidly [44][45][46].For certain piecewise constant and continuous quasi-periodic drives, this is known to be possible in the high-frequency regime [47][48][49][50][51][52][53][54][55][56][57][58].Rigorous bounds on heating rates can also be established by generalizing the Floquet theory [59,60].However, this becomes obscure for kicked systems as the high-frequency limit of kicks in principle allows a divergent rate of energy input into the system.Aperiodically kicked systems have been most limited to few-body settings [61][62][63][64][65] and it remains an outstanding challenge to control heating in the thermodynamic limit.
In this work, instead of focusing on the high-frequency regime, we exploit the self-similarity inherent in the RMD sequence to demonstrate that heating can be parametrically controlled by the kick strength.Through a perturbative expansion, we derive an effective Hamiltonian that governs the initial time evolution.Remarkably, the self-similar multipolar structure leads to exact cancellations of numerous terms in the effective Hamiltonian, thereby eliminating the corresponding heating channels.This mechanism of heating suppression is independent of the specific model and is applicable to both quantum and classical many-body systems.
For numerical efficiency, we demonstrate this effect in a concrete model, namely a kicked chain of classical rotors.Starting from low-temperature initial states, the system exhibits a long-lived prethermal regime before heating up.The lifetime scaling depends on the microscopic details of the kicked system and if the static part only involves the kinetic energy, the lifetime scales as a power law with a tunable exponent 2n.This we account for by analysing the linear stability of the system.In the quasiperiodic TM limit, we also show that the lifetime grows faster than any power law but slower than exponentially.
In the following, we first consider a general kicked system and demonstrate heating suppression in a perturbative expansion of the effective Hamiltonian in the small kick strength.We then present the results on the chain of rotors before a concluding discussion.
Setting.-Consider the time-dependent Hamiltonian H(t) = H + V ∆(t), where H denotes the static part and V defines the kick with ∆(t) = l K l δ(t − lτ ), and the kick strength K l with period τ .We focus on the intermediate frequency regime, i.e. τ is not necessarily small.Suppose the strength K l has binary choices ±K following an n-RMD sequence.For n = 0, there are two possible unitary time evolution operators: We can formally define the time-independent effective Hamiltonian H ± 0 through the relation . For weak kick strength K, we can perturbatively construct the effective Hamiltonian as , describing the dynamics at times t = lτ .Although such an expansion may diverge for many-body systems, we expect that its truncation at low orders in K will approximate the initial time evolution.The lowest order term is simply the static Hamiltonian, Ω ± 0,0 = H.Via the replica resummation of the Baker-Campbell-Hausdorff series, the leading correction can be expressed in a compact form [15] where ad X (Y ) = [X, Y ] is the Lie derivative.It can also be expanded in a power series in τ as where higher order terms only contain nested commutators of the form [H, V ] s := [H, . . ., [H, V ] . . .] with a single kick V but multiple (s) H operators.Generally, Ω 0,1 does not vanish, and initially, the time evolution is dominated by H ± KΩ 0,1 .The term Ω 0,1 occurs randomly with an amplitude linear in K, and one expects it to quickly destabilize the system and induce heating in a short time.
We now use the self-similar structure of RMD protocols to show that many terms of order O(K) in the effective Hamiltonian can be eliminated.Furthermore, if the condition can be satisfied for ∀s ≥ n c with some integer n c , for n−RMD systems with any n ≥ n c , random perturbations start appearing at a higher order O(K 3 ) and hence heating can be significantly suppressed.
To see this, we first observe that higher-order multipolar operators can be recursively obtained using the relation where U ± n generates the time evolution over the duration 2 n τ [51].For n-RMD systems, the time evolution is given by a random sequence of multipolar operators U ± n .Similarly, the effective Hamiltonian H ± n is defined through , governing the stroboscopic time evolution (t = 2 n τ l for integers l).The perturbative expansion is denoted as Notably, the time evolution operators in Eq. 1 possess the special property that U + 0 can be mapped to U − 0 by changing K → −K.Thus, terms in the effective Hamiltonians coincide for even orders in K, while differing by a minus sign for odd orders, given by Similar to the purely random drive (n = 0), the initial stroboscopic time evolution is governed by H ± KΩ n,1 , and the system may still exhibit rapid heating.However, it is noteworthy that the self-similar construction in Eq. 5 and the symmetry property in Eq. 6 lead to an important observation: several terms in Ω n,1 actually vanish, resulting in the remarkable property where the summation starts from s = n, although obtaining the coefficient f n,s can be a cumbersome task.Importantly, Eq. 7 suggests that, to the leading order of O(K), heating can only occur through heating channels in the form of [H, V ] s with s ≥ n, while all other heating channels are strictly forbidden.The derivation of this expression is presented in Sec.SM 1 of the Supplementary Materials (SM).Now, if the condition given by Eq. 4 is satisfied, all terms in Ω n,1 vanish.Consequently, the stroboscopic time evolution of the system is effectively governed by the Hamiltonian H ± n = Hn ± O(K 3 ), where the static part is denoted by Hn = H + K 2 Ω n,2 .Therefore, the RMD kicked systems first relax to a prethermal ensemble determined by Hn before notable heating is induced by random perturbations of order O(K 3 ).
Although we use the perturbative expansion for quantum systems, it is important to note that this mechanism of heating suppression equally applies to classical many-body systems.The Liouville equation, which describes the phase-space distribution of a classical system, exhibits a structural similarity to the Schrödinger equation in quantum systems.Consequently, the effective Hamiltonian for classical systems can be obtained by formally replacing the commutator [. . .]/i in its quantum counterpart with the Poisson bracket {. . .} [18].Due to the computational efficiency of numerical simulations for large classical systems, we proceed to demonstrate this heating suppression in a classical rotor system.
Many-body kicked rotors.-Weconsider many-body rotors with the static kinetic energy (H = H kin ) and the kicked nearest-neighboring interactions (V = V int ), where p j and q j for j = 1, ..., N are the conjugate angular momenta and angles of N rotors, respectively.Periodic boundary conditions are used (q 1 = q N +1 ).The interaction preserves the total angular momentum N j=1 p j .As we will demonstrate below, the condition in Eq. 4 can be approximately satisfied when the angular momentum distribution is narrow.In the prethermal regime, the width of the distribution is determined by the temperature T , which is controllably small and scales with the kick strength as T ∼ K 2 .
To begin, we derive the nested Poisson brackets {H, V } s for kicked rotors, which reduce to j (p j − p j+1 ) s sin(q j −q j+1 ) for odd s and j (p j −p j+1 ) s cos(q j − q j+1 ) for even s.Therefore, for multipolar order n ≥ 1, the expression in Eq. 7 implies that the dominant random perturbation Ω n,1 only contains terms proportional to (p j − p j+1 ) or its higher powers.These terms become negligible when the kinetic energy distribution is sufficiently narrow.The suppression becomes stronger with higher multipolar orders, and Eq. 4 can be more effectively satisfied for larger multipolar order n.
When we start from the initial condition p j = p for all j, Eq. 4 is fulfilled exactly, and the initial time evolution is governed by the static effective Hamiltonian Hn = H kin + K 2 Ω n,2 .As this Hamiltonian is generally nonintegrable, the angular momentum distribution spreads.In the prethermal regime, it approximately reaches the Gibbs distribution N j=1 exp −(p j − p) 2 /2T [43].The width of the distribution is determined by an effective prethermal temperature T .As shown in Sec.SM 3, this temperature can be controlled to be small for weak kicks (T ∼ K 2 ).
The spreading of the angular momentum distribution can be quantified by the kinetic energy density E kin (t) := , making it a suitable measure of heating and temperature increase.The initial conditions are chosen such that the angles q j uniformly distribute between 0 and 2π, and the angular momentum p j = 0.1 is fixed for all rotors, satisfying the condition in Eq. 4. In Fig. 1(a), we depict the time evolution of the averaged kinetic energy E kin , averaged over 350 noise realizations with different initial states, for a fixed kicking strength K = 0.03 and rotor number N = 500 [67].For multipolar order n ≥ 1, the averaged kinetic energy remains almost unchanged for a long timescale t * .However, as the kinetic energy is unbounded, it eventually increases when heating takes over.We observe that for larger n, the timescale t * remarkably extends by several orders of magnitude, reaching its largest value in the quasi-periodic TM limit.In contrast, for the fully random drive (n = 0), unbounded diffusion starts even at very early times, and no prethermal regime can be established.
We quantify the prethermal lifetime t * and its dependence on the kicking strength K. To extract t * , one can fit the averaged kinetic energy up to time t f with a power law t b and monitor the power b for different t f [43].During the prethermal regime, the power b remains close to zero, and t * is determined when b first reaches a threshold.In our numerical simulations, we choose b = 0.05, but our findings are independent of the specific threshold value as long as it is small.Fig. 1(b) illustrates the dependence of t * on the kick strength for different multipolar orders.Using a log-log scale, a linear fit suggests that the prethermal lifetime follows an algebraic dependence on the kick strength, t * ∼ (1/K) α .The scaling exponent α can be determined through numerical fitting.For n = 0, the exponent is close to 2. It remains approximately the same for n = 1, although the prefactor differs by three orders of magnitude, indicating significant suppression of heating due to the dipolar structure.Interestingly, for higher multipolar orders, α notably increases and exhibits a good approximation to the relation α ≈ 2n for n = 1, 2, 3.In the TM limit, the lifetime scaling converts to where the constant C ≈ 0.8 and g ≈ 0.3 as shown in Fig. 1(c).A similar functional form has been reported in Ref. [60] but in the high-frequency regime.We verify that this scaling grows faster than any power law (cf.Sec.SM 5), indicating a significant suppression of heating in a non-perturbative manner.
Linearization.-Although it is expected that higher multipolar orders may further suppress heating, the perturbative expansion of the effective Hamiltonian is insufficient to explain the scaling of the prethermal lifetime.To address this, we develop a simple theory by linearizing the many-body systems.
Assuming small angular differences between neighboring rotors, (q j − q j+1 ) mod 2π 1, we can expand the kicked interaction using the quadratic approximation cos(q j − q j+1 ) ≈ 1 − 1 2 (q j − q j+1 ) 2 in the Hamiltonian Eq. 8. Performing a Fourier transform, we obtain where w := 2πI/N for an integer I, F (w) := 4K sin 2 (w/2), and the Fourier components are defined as √ N .The ± sign follows the RMD sequence.The system now decouples into a set of independent kicked harmonic oscillators labeled by w.For each oscillator, we can integrate its discrete time evolution in a two-dimensional phase space over one period τ : where M ± 0 is the elementary evolution matrix (see derivations in Sec.SM 2 A): where we drop the label w as the following discussion equally applies to all w.Similar to Eq. 5, higher multipolar evolution matrices can be recursively derived as generate stroboscopic time evolution over the duration 2 n τ .Crucially, both M + n and M − n have the property det(M ± n ) = 1, making them area-preserving maps [68].Therefore, when only M + n or M − n is periodically applied, the system, for weak kicking strength, exhibits non-chaotic dynamics confined to a closed elliptical orbit around its fixed point (Q, P ) = (0, 0).However, the random concatenation of two slightly different maps M ± n generally perturbs these stable trajectories, causing them to deviate indefinitely from their fixed points (Fig. 2(a), blue).By quantifying such deviation, one can estimate the heating rate and its relation to the multipolar order.
To analyze this deviation, we define the averaged evolution matrix as Mn := 1  2 (M + n + M − n ) and the difference between the two matrices as ) for non-zero n, indicating that the averaged map Mn does not preserve area in phase space.Instead, the trajectory slowly spirals out with a constant expansion rate scaling as F 2n (magenta in Fig. 2(a)).Additionally, the stochastic term D n possesses eigenvalues that scale as F n (cf.Sec.SM 2 C), and one would expect it to contribute to a diffusive spiral-out process with a rate also scaling as F 2n .
To quantify this process, we introduce the normalized map M = Mn / det Mn , ensuring its area-preserving property with det( M ) = 1, thus generating a closed elliptical orbit (Fig. 2(a), black).The matrix elements of M define the metric of the orbit and determine its conserved area A(Q, P ) [69].The radius of the ellipse, defined as r h = A(Q, P )/π, becomes time-dependent when M ± n is stochastically applied h times.The expansion rate of the radius can be calculated as ∆r h /r h , where ∆r h = r h+1 − r h .By averaging over different random realizations and the polar angle of the ellipse, we find that its leading order contribution scales as F 2n , with a specific expression ∆r h /r h ≈ 3τ 2 F 2 /4 for n = 1 and 6τ 4 F 4 for n = 2, as detailed in Sec.SM 2 B. Consequently, the averaged growth of the radius at early times can be obtained accordingly.
In Fig. 2(b), we present numerical simulations (circles) of the averaged radius for n = 1 (blue) and 2 (orange), which closely match our analytical predictions (dashed lines).As F is proportional to the kicking strength, the expansion rate scales as K 2n , and its inverse corresponds to the observed prethermal lifetime scaling in Fig. 1.
We note that the strong dependence of the multipolar order n in the prethermal lifetime scaling is remarkably robust, even for initial states that deviate significantly from the linearization regime where (q j − q j+1 ) mod 2π 1. Indeed, our numerical results in Fig. 1 are obtained using a random distribution of q j over a wide range [0, 2π].In Sec.SM 4, we confirm that this phenomenon persists as long as the prethermal regime exhibits a low temperature, leading to a narrow distribution of angular momenta.
Discussion.-We have proposed a mechanism to suppress heating in aperiodically kicked systems by intro- ducing self-similar multipolar structures.This mechanism brings about significant changes in the thermalization pathways and effectively blocks a series of heating channels.As a result, it supports the existence of a longlived prethermal regime even in the absence of TTS and away from the high-frequency regime.
To demonstrate this mechanism, we have considered classical many-body rotor systems, where we have discovered a characteristic prethermal lifetime scaling of (1/K) 2n .In the quasi-periodic TM limit, the heating suppression becomes non-perturbative, leading to the scaling given by Eq. 10, which does not follow an exponential or algebraic form.A similar functional form has been rigorously proven in the high-frequency regime [60].However, it remains an interesting open question to justify such a dependence on the kicking strength.
The Hamiltonian Eq. 8 can be experimentally realized, e.g., using an array of bosonic Josephson junctions [42,70,71].This opens up possibilities for the experimental exploration of prethermalization with RMD kicks.
It is important to note that while heating channels can be suppressed by the RMD sequence, the lifetime scaling is not universal and can strongly depend on the microscopic details of the kicked system.The strong dependence of the multipolar order n in the lifetime scaling may not occur if interaction terms are also present in the static Hamiltonian, such as H = H kin + V int .The leading order perturbation Ω n,1 involves terms with more than one V int , e.g., {V int , {V int , H kin }} ∼ j [sin(q j − q j+1 ) − sin(q j−1 − q j )] 2 .These terms are independent of angular momenta and cannot be suppressed even at low prethermal temperatures.Hence, the condition Eq. 4 cannot be satisfied in this case.
We implemented a kicked protocol with a modified kick strength K l = ±K + B such that the additional static interaction can be efficiently simulated at stroboscopic times [72].In Fig. 3(a), we illustrate the results with B = 0.01.The prethermal plateau is still observed, and the corresponding lifetime is shown in panel (b).It is evident that for n > 0, heating can still be significantly suppressed.However, its dependence on the kicking strength now follows t * ∼ (1/K) 2 regardless of the multipolar order.A similar linearization analysis can be performed, and the expansion rate for each decoupled oscillator is K 2 , as detailed in Sec.SM 2 C. Identifying a general mechanism for further suppressing heating with a better scaling remains an intriguing open question.
Finally, we note that the perturbative expansion predicting the suppression of heating also applies to RMD kicked quantum systems.A systematic study of quantum thermalization in kicked systems and its relation to their classical counterparts is an intriguing subject for future study.For the quantum kicked system, we have two different unitary time evolution operators

Supplementary Material Prethermalization in aperiodically kicked many-body dynamics
One can obtain a time-averaged Hamiltonian H ± ave = H ± K τ V as the effective Hamiltonian to approximate the early time dynamics.However, this is not a suitable expansion if τ is not sufficiently small and one wants to study the perturbative expansion with respect to the kick strength K.For instance, the term [H, V ] s all have an amplitude scaling of O(K) but they are not captured in the averaged Hamiltonian.Instead, it is necessary to perform the replica resummation, whose general expression can be cumbersome to obtain, but there is a systematic approach to achieve it [15], see also examples in [16].For the same reason, the previous heating analysis on RMD systems in the high-frequency regime and the expansion of order O(τ ) will not be applicable here.
To explore the heating effect in RMD kicked systems, we consider the expansion in O(K) as for all n.For n = 0, a systematic method has been established for constructing the expansion ))] and the O(K) term is presented in Eq. 2 [15] .For larger values of n we still begin by considering the leading-order correction with m = 1, and we assume a general structure for n ≥ 1 as From Eq. S.2 we know for any n and s, we have f + n,s = −f − n,s for m = 1.Note, we do not require the specific expression for each coefficient f ± n,s .Instead, it suffices to demonstrate that some of them become exactly zero, thereby prohibiting certain heating channels.A similar expansion can be derived for higher-order multipolar operators where g n,s,l are some coefficients, and importantly, f ± n,s + f ∓ n,s = 0 in the first line cancels.By comparing it with the assumption Eq.S.3 but for n = 1, we have and by matching the coefficients of τ s , we can establish the following relation for the coefficients and obtaining f + n,s for s ≥ n + 1 can be cumbersome.Therefore, the first line in Eq.S.6 implies Eq. 7 in the main text.It suggests that, to the leading order of O(K), heating can only occur via the process [H, V ] s with s ≥ n, while all other heating channels are strictly forbidden.
Higher-order terms with an even order of K do not introduce random perturbations.However, they still contribute to heating in the form of Arnold diffusion, similar to periodically driven systems, but their contribution is exponentially small in the kick strength [43].As a result, the next significant random heating channels emerge at order O(K 3 ).It can also be shown that the self-similarity of the RMD sequence leads to the exact suppression of these heating channels.To see this more easily, we consider a special case where higher-order nested commutators have negligible contributions for ∀s ≥ n c with a certain integer n c .Consequently, Ω n,1 = 0 and the stroboscopic time evolution of the system is effectively governed by the Hamiltonian where for n > k and a certain integer k.By using

one can again observe the vanishing coefficients
(S.9)

SM 2. LINEARIZATION OF THE MANY-BODY HAMILTONIAN A. Time evolution matrix
Following [73], we can express the Hamiltonian of our model as a collection of decoupled kicked harmonic oscillators in a quadratic approximation: cos(q j − q j+1 ) ≈ 1 − (q j − q j+1 ) 2 /2, provided that the two neighbouring rotor angles are sufficiently close (q j − q j+1 ) mod 2π ≈ 0. Thus, we have where w := 2πI/N , F ± (w) := 4(B ± K) sin 2 (w/2) (the choice of F ± (w) depends on the RMD sequence), q j e −iwj are the Fourier transforms of p j and q j , respectively.Note that Eq. 11 in the main text is a simplified version of the Hamiltonian above with B = 0.Here we use the general expression with non-zero B such that the linear stability of the dynamics in Fig. 3 can also be discussed.
Since p j and q j are real, we have P * w = P −w and Q * w = Q −w (here star denotes complex conjugate).For each w, the classical equations of motion are given by where Ω(t) = +∞ k=−∞ δ(t − kτ ).Consider the evolution of the system over one time period, from t = − to t = τ − with τ .The solution to the above equation is (S.12) During the first part of the period (when t ∈ (− , )), the rotor is kicked and the time evolution is determined by the matrix In the second part of the time period (when t ∈ ( , τ − )), the rotor experiences a free motion, described by the matrix (with → 0) As a result, the phase space evolution of the kicked rotor over one time period is given by the matrix Notice that the matrix M ± w can be reduced to 2 × 2 matrix where the subscript w is dropped from now on.The evolution matrices for higher multipolar orders n can be derived recursively as , which determines the time evolution of duration 2 n τ .For example, when B = 0, n = 1 we have where we have denoted

B. Stability of integrable orbits
We use the method proposed in Ref. [68] to analyse the stability of the elliptical orbits.Let us denote such that M ± n = Mn + ξD n where ξ is a random variable, being either +1 or −1 with the same probability.Hence, its average vanishes ξ = 0 and the variance reads ξ 2 = 1.We note that det( Mn ) = 1 + O((τ F ) 2n ) for non-zero n, implying that the averaged map Mn is not area-preserving.We therefore define a new matrix M := Mn / det Mn so that det( M ) = 1 and this new matrix can be used to define the area of a closed orbit in the linearized system.This orbit is generally a rotated ellipse centered around the fixed point (0, 0).Note, M also depends on the multipolar order n but the following method equally applies for all n.For now we drop it for simplicity.With the matrix elements M ij of M its area can be defined as [74] A which is conserved if M is repeatedly applied.Q w and P −w are generally time-dependent, and in the following, we drop w for simplicity and introduce h to label their time-dependence.We only focus on stroboscopic time evolution and use (Q h , P h ) to represent the trajectory at time h2 n τ .One can use the polar angle φ h to parametrize the points (Q h , P h ) on the rotated ellipse as where R q/p defines length of the major or minor axis, and θ defines the rotation angle with respect to the axis.It can be determined by One can also convert the variables back as The lengths of the major and minor axes are given by with the constant It is worth noting that for n = 1, R q ∼ O(F −1 ), R p ∼ O(F 0 ), so for a weak kick strength, R q can be large.This stretches the ellipse in the Q-direction much more strongly than in the P -direction.
For RMD drives where M ± n is applied stochastically, the area of the closed orbit becomes time-dependent.For a single random realization and at a certain time, this area can either expand or contract.However, if we average over many different random realizations, it generally expands.We can quantify this expansion by first defining the ellipse's radius r h = A(Q h , P h )/π, and calculating the expansion rate ∆r h /r h , where ∆r h := r h+1 − r h .We use the same metric to define the ellipse's area but now the trajectory updates stochastically as The expansion rate now reads We now insert Eq.S.25 into Eq.S.26 to obtain the general expression as a function of the polar angle φ h , the kick strength F and the kick duration τ .Unfortunately, it is usually very complicated and not enlightening.However, if the kick strength F is small, one can perform a Taylor expansion to obtain the most relevant contributions.For our purpose, a Taylor expansion up to the order O(F 2n ) would be sufficient.This process can be done by employing symbolic computation tools such as Wolfram Mathematica.
For n = 1, to the second order of F , we have Terms that are linear in ξ vanish after averaging over different random realizations.Contributions that are quadratic in ξ generally do not vanish, unless for special polar angles, e.g., cos 2 2φ h = 1.Also note that there is a term of order O(F ) which does not depend on ξ.It arises from the fact that the average map Mn is not area-preserving, contributing a constant expansion rate as shown in Fig. 2 (the magenta line).
To further remove the angular dependence in the expansion rate, one can assume that the change in the radial direction is much slower than in the angular direction and integrate the angle φ h over the full range [0, 2π].We further use ξ = 0 and ξ 2 = 1 to obtain the averaged expansion rate that is proportional to F 2 (S.28) This provides an approximation for the mean radius evolution for small h (and small τ F ):

C. Eigenvalues of the matrices Mn and Dn
Instead of rigorously calculating the expansion rates and their dependence on n, one can also estimate them by studying the scaling of eigenvalue properties of the update matrix.We first assume B = 0 and then obtain Mn recursively as (S.32) Since τ 2 F 2 8 the eigenvalues of Mn for each n > 0 form a complex conjugate pair.Furthermore, as det( Mn ) > 1, the deviation from one determines the rate of constant expansion of the dynamics generated by Mn .Specifically, we calculate det( Mn ) for different values of n, which correspond to the norms of the eigenvalues of M n , (S.33) Note that for a weak kick strength, we have (| λn,± | − 1) ∼ (τ F ) 2n for a small value of n, and higher-order terms are negligible.Therefore, we expect the constant expansion rate to scale as F 2n as discussed in the previous section Sec.SM 2 B. However, we also observe that the prefactor for higher powers of F tends to increase for n = 3, suggesting that the perturbative expansion in orders of F may not converge for a large value of n.
We also compute the matrix D n (S. 34) with eigenvalues which scales as |µ n,± | ∼ (τ F ) n for n > 0. As D n appears stochastically in time, we expect its leading-order contribution to vanish.Its second-order effects lead to a diffusive spiral-out process with an expansion rate that scales as F 2n .We perform a similar calculation for non-zero B and by defining F ± ∼ (±K + B), we obtain Importantly, one always finds (| λn,± − 1|) ∼ K 2 , where the scaling exponent does not depend on the multipolar order.This scaling relation corresponds to the observed heating rate scaling of K 2 in Fig. 3 in the main text.

SM 3. TEMPERATURE AT THE PRETHERMAL STAGE
Here we analyse the dependence of the prethermal temperature on the kick strength and we demonstrate that it follows T ∼ K 2 .For the initial condition p j = p, the initial kinetic energy density is given by p2 /2.In the prethermal regime with a weak kick, we assume that the distribution for p j and q j decouples [43].The angular momentum distribution approaches the Gibbs distribution Z −1 N j=1 exp −(p j − p) 2 /2T with a normalization factor Z −1 .The corresponding kinetic energy density is E * kin = (T + p2 )/2.We numerically study the dependence of the temperature T on the kick strength K in Fig. S1 (a) and (b) for two different initial conditions: p = 0.1 and p = 0, respectively.The numerical results fit well with a straight line in a log-log scale, with a slope around 2. This suggests that the temperature follows T ∼ K 2 .In contrast, for the modified kicked protocol used in Fig. 3 with a non-vanishing B, the averaged prethermal kinetic energy linearly depends on B but does not notably change with K. We verify this in Fig. S1, panel (c).

SM 4. DISTRIBUTION OF ANGULAR MOMENTA
In the main text, we linearize the many-body Hamiltonian and explain the characteristic scaling of the prethermal lifetime.We note that such scaling can be very stable and persist even away from the linearization regime (q j − q j+1 ) mod 2π 1, as shown in Fig. 1 where a wide initial distribution of q j is used.Interestingly, we notice that this lifetime scaling is not sensitive to the angular dependence, but strongly relies on the angular momentum distribution during the prethermal regime.As discussed in Sec.SM 3, this distribution is governed by the prethermal temperature.We find that, as long as the prethermal regime exhibits a low temperature, or equivalently, a narrow distribution of angular momenta, the prethermal regime can be sufficiently long-lived, and the dependence on n should manifest.The prethermal temperature can be adjusted by the initial condition.For instance, we consider initial an angular momentum distribution following a Gaussian distribution with a zero mean and a standard deviation σ.A larger standard deviation generally increases the prethermal temperature, resulting in a broader angular momentum distribution during the prethermal regime.This is confirmed in Fig. S2, panel (a), where three different values of the initial standard deviation σ are used.Note that we use n = 1 to generate the dynamics but this figure qualitatively represents other multipolar order as well.The kick strength K is chosen such that the prethermal lifetimes are approximately the same for all σ, with values of K set as 0.012, 0.008 and 0.005 for σ = 0.001, 0.01 and 0.1, respectively.
The probability distributions are extracted at t = 1500 just before the system notably heats up.
We now illustrate the dependence of the prethermal lifetime scaling on different initial conditions, focusing on the TM drive.For a fixed kick strength K, it typically determines the longest possible prethermal lifetime for the entire family of n-RMD protocols.Therefore, if we fit t * versus 1/K on a log-log scale, the scaling exponent also sets the upper bound for other n-RMD protocols.As long as the TM drive exhibits a sufficiently large scaling exponent α, n-RMD with any finite n should exhibit the n-dependence in the lifetime scaling.
In Fig. S2, panel (b), we present the prethermal lifetime scaling for different initial conditions.For narrow distributions, such as σ = 0.001 and 0.01, the scaling exponents are still very large, approximately α ≈ 7.9 and 7.2, respectively.We expect that these fitted scaling exponents may increases further if we perform the fit using larger 1/K and longer time windows, similar to Fig. 1(c).However, for a larger standard deviation, such as σ = 0.1, the scaling exponent notably decreases to α ≈ 3.1, and we expect that n-RMD systems with finite n would heat up faster.
Therefore, we expect that as long as the prethermal regime has a low temperature, a long-lived prethermal regime and the n-dependence in the lifetime scaling should emerge.Further systematic investigations of the temperature dependence of the prethermal lifetime scaling will be explored in future work.

SM 5. SCALING OF THE THUE-MORSE PRETHERMAL LIFETIME
In the main text we show that for the TM drive, the lifetime scaling becomes t * ∼ exp(C[ln(K −1 /g)] 2 ) with constants C and g.Here, we compare this result with other fitting methods.For instance, in Fig. S3(a) we use a log-log scale and clearly the numerical data tends to curve up.In contrast, panel (b) depicts the same data but in log scale and the numerical result bends down.Therefore, this scaling grows faster than any power-law but slower than exponentially.

A. Finite size effects
In Fig. S4 we compare the dynamics using the TM drive for different system sizes.The simulations converge as the system size increases.In the main text, we use N = 500 to generate the data, which is already sufficient to mimic the heating behaviour in thermodynamically large systems.

FIG. 1 :
FIG. 1: (a) Time evolution of the averaged kinetic energy for n−RMD and Thue-Morse (TM) drive with K = 0.03 in a log-log scale.(b) Dependence of the prethermal lifetime t * on 1/K in a log-log scale.Dashed lines (K −2n for n > 0 and K −2 for n = 0) are plotted to guide the eyes.(c) Prethermal lifetime t * scaling for TM drives.

FIG. 2 :
FIG.2:(a) Trajectories in phase space.The black orbit is obtained by the area-preserving map M .The magenta curve with a constant expansion rate is generated via M1 .The blue curve is a single realization obtained by stochastically applying the matrix M ± 1 .(b) The averaged radius r h matches well with the theoretical prediction (dashed lineas).F = 0.08 is used in both panels.

FIG. 3 :
FIG. 3: (a) Time evolution of the averaged kinetic energy density for different multipolar order n with K = 0.005 in a log-log scale.(b) Prethermal lifetime t * as a function of 1/K in a log-log scale.The system starts from the initial angular momentum p j (0) = 0 with B = 0.01.Other numerical parameters are same as in Fig. 1.Dashed lines correspond to the scaling K −2 .

CONTENTS SM 1 . 8 SM 2 . 12 SM 3 . 13 SM 4 . Distribution of angular momenta 14 SM 5 .
Many-body effective Hamiltonian for kicked systems Linearization of the many-body Hamiltonian 9 A. Time evolution matrix 9 B. Stability of integrable orbits 10 C. Eigenvalues of the matrices Mn and D n Temperature at the prethermal stage Scaling of the Thue-Morse prethermal lifetime 15 A. Finite size effects 15 SM 1. MANY-BODY EFFECTIVE HAMILTONIAN FOR KICKED SYSTEMS FIG. S1: The averaged prethermal kinetic energy density E * kin for n = 1 RMD for (a) B = 0 with p = 0.1, (b) B = 0 with p = 0 and (c) B = 0.01 with p = 0.The prethermal temperature scales as T ∼ K 2 for (a) and (b) when the kick is weak.In the last case it is fixed by B, so all the plateau parts of energy curves are at the same value E * kin = 0.0048.The rest numerical settings are the same as in Fig.1.All plots use a log-log scale.
FIG. S2: (a) Momentum distribution at the end of the prethermal stage for a Gaussian initial momentum distribution with a zero mean and a standard deviation σ indicated in different colours, and (b) corresponding prethermal lifetime t * as a function of 1/K for n = 20 in log-log.The rest numerical settings are the same as in Fig.1.

8 FIG
FIG.S3: (a) and (b) depict the power-law and exponential fitting of the prethermal lifetime t * (1/K) for the TM drive.
FIG. S4: Time evolution of the averaged kinetic energy for the TM drive.The simulation results converge for large systems and N = 500 is already sufficient to produce thermodynamically large systems.Here we use the kick strength K = 0.07 and the results are averaged over 200 random realizations.Initial states are the same as in Fig. 1.