{'fraction_non_alphanumeric': 0.040573341168655304, 'fraction_numerical': 0.013388180587138148, 'mean_word_length': 4.323721436343852, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stability and pre-thermalization in chains of classical kicked rotors

Periodic drives are a common tool to control physical systems, but have a limited applicability because time-dependent drives generically lead to heating. How to prevent the heating is a fundamental question with important practical implications. We address this question by analyzing a chain of coupled kicked rotors, and find two situations in which the heating rate can be arbitrarily small: (i) linear stability, for initial conditions close to a fixed point, and (ii) marginal localization, for drives with large frequencies and small amplitudes. In both cases, we find that the dynamics shows universal scaling laws that allow us to distinguish localized, diffusive, and sub-diffusive regimes. The marginally localized phase has common traits with recently discovered pre-thermalized phases of many-body quantum-Hamiltonian systems, but does not require quantum coherence.


Introduction
Periodically driven many-body systems are currently being studied intensively, both in theory and in experiments. One common research direction, known as Floquet engineering, uses periodic drives to generate tunable effective couplings. For example, periodic drives can create negative kinetic energies, leading to dynamical localization and kinetic frustration in experiments on ultracold atoms (see Ref. [1] for a recent review). Similarly, periodically oscillating phonons are believed to be responsible for the light-enhanced coherence observed in superconductors [2,3,4,5,6,7]. Periodic drives can additionally be used to create phases of matter that do not have any equilibrium counter-part, such as Floquet topological phases [8,9,10,11,12,13,14] and Floquet time crystals [15,16,17,18,19,20].
The key difference between time-independent and periodically driven Hamiltonian systems is that in the latter, the energy is not conserved. As a consequence, Floquet systems are generically subject to heating, as their average energy increases with time. For many-body closed quantum systems, one can rely on an extension of the eigenstate thermalization hypothesis (see Ref. [21] for an introduction), which is valid at asymptotically long times. Because the energy is not conserved, periodically driven systems are expected to effectively thermalize at an infinite temperature [22,23,24,25,26,27].
In this article, we explore a classical periodically driven many-body system, by considering a chain of coupled kicked rotors. As we will see, this simple system shows many features that were predicted earlier for the more complex quantum case. Our model is a natural extension of the single kicked rotor, whose dynamics is described by the well-known "Chirikov standard map" [55]. This single-particle model offers a paradigmatic example of a transition between regular dynamics and chaos, or equivalently localization and diffusion. Earlier studies revealed that when many rotors are coupled together, the transition is washed out: at exponentially long times, manybody kicked rotors are generically diffusive [56,57].
By focusing on the finite-time dynamics, we find regions in the parameter space that are localized, diffusive, and super-diffusive. These three regimes are characterized by a different behavior of the average kinetic energy of the rotors. In the localized regime, this quantity does not grow with time, indicating that the system does not heat up. In contrast, in the diffusive regime, the system absorbs energy at a constant rate. The transition between the localized and the diffusive regimes depends on the initial conditions. If the system is initialized close to a stable point, the transition corresponds to a sharp resonance, which does not drift with increasing time. If the system is instead initialized with large fluctuations, the critical coupling of the transition depends logarithmically on time. When logarithmic corrections are taken into account, the diffusion parameters show a universal asymptotic scaling, as a function of the kick strength and frequency.
The paper is organized as follows: We first briefly review the physics of a single kicked rotor (Sec. 2) and present its many-body generalization (Sec. 3). Then, we describe the case where the rotors are prepared in vicinity of a stable fixed point, and discuss the system's linear stability and resonances, both analytically (Sec. 4) and numerically (Sec. 5). Next, we address generic initial conditions and study the regime of marginal localization (Sec. 6). Finally, we characterize the different regimes using a Fourier decomposition (Sec. 7), and draw our conclusions (Sec. 8).

Review of a single kicked rotor
The Hamilton of a single kicked rotor is given by Here, the position of the rotor is described by the angle φ, and the associated momentum by p. The parameters κ and τ are the strength and the period of the kicks, respectively. For κ > 0, the static model (Ω(t) → Ω 0 > 0) has a stable fixed point at p = φ = 0, and a unstable one at p = 0, φ = π. Note that this model can be mapped to the κ < 0 model under the transformation φ → φ + π, which flips the role of the stable and unstable fixed points. Thus, without the loss of generality, we focus here on the case κ > 0 only. We now briefly review the linear analysis of Eq. (1) for the sake of completeness, and due to its relevance to the many-body case. We expand cos φ around φ = 0 and apply the quadratic approximation cos φ ≈ 1−φ 2 /2 in Eq. (1). The classical equations of motion are then equivalent to a kicked harmonic oscillator [58], where d dt We consider the evolution of the system over one time period, from t = −ε to t = τ − ε with ε τ . The solution of Eq. (2) is given by where Υ denotes time ordering and M (t) is 2 × 2 matrix of Eq. (2). During the first part of the period, for −ε < t < +ε, the rotor is kicked and the time evolution is determined by the matrix In the second part of the time period, ε < t < τ − ε, the particle 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 The transition between the stable and unstable dynamics of the oscillator can be inferred from the eigenvalues of the matrix M , Eq. (6), which are given by We find that the properties of the eigenvalues are determined by the unitless parameter K = κτ : for K < 4 the eigenvalues are a complex conjugate pair λ ± = exp(±iθ) with cos θ = 1 − K/2, while for K > 4 they are a real reciprocal pair λ ± = exp(±β) with cosh β = 1 − K/2. Thus, the point K = 4 separates the stable (K < 4) and unstable (K > 4) regions, and can be associated with a parametric resonance of the system ‡.
Note that the unitless parameter K = κτ sets the ratio between the driving amplitude and the driving frequency: as expected, the system is stable for weak and fast drives, and unstable for strong and slow drives. The model (1) is furthermore known to undergo a transition between a localized and a diffusive regime at K = K c ≈ 0.9716 [59]. This phenomenon can be explained analyzing the evolution of the phase space as a function of K. For K = 0, the system is integrable and p is a conserved quantity. In this case, the orbits are horizontal lines that slice the phase space in distinct regions. According to the Kolmagorov-Arnold-Moser (KAM) theorem [60], for small values of K, the system's dynamics is governed by stable invariant submanifolds (KAM tori). These KAM tori divide the phase space into many sections, each of which may carry stochastic orbits. Therefore the momentum is bounded between these invariant tori and the kinetic energy does not grow with time (localized regime). As the value of K is increased, the number of KAM surfaces decreases and, at some critical value K = K c , the last KAM tori disappears. For K < K c the momentum of the pendulum is localized between invariant these KAM tori, independently of the initial conditions, while for K > K c the system can show a diffusive behavior (in momentum space) for some initial conditions. Note that this transition occurs in the region where the point φ = 0 is linearly stable. The distinction between linear stability and localization will be a key ingredient for our analysis of the many-body case.
The Hamiltonian in Eq. (1) can be quantized, by assuming p and φ to be quantum operators satisfying canonical conjugation relations. For this system, the quantum evolution operator over one period is the product of free evolution and instantaneous kick operators, which is called Floquet operator. The stroboscopic evolution of the system is completely determined by the eigenstates of the Floquet operator. Using the eigenvalue equation of the Floquet operator, it has been shown before that a quantum kicked rotor can be mapped to a quantum particle moving in a one-dimensional static disordered potential [61,62]. Because of Anderson localization [63], all the eigenfunctions are localized in momentum space and the kicked rotor cannot diffuse. Thus, the quantization of Eq.
(1) leads to a model that is always localized, irrespective of the value of K. This effect ‡ The distinction between the stable and unstable regime is due to the fact that the matrix M is real and has determinant detM = 1. These conditions lead to two distinct options only: either both eigenvalues are real and have reciprocal values (λ and 1/λ), or they are pure phases (e iθ and e −iθ ). The transition between these two situations occurs at an exceptional point, where the two eigenvalues coalesce (i.e. are both equal to 1, or to −1). Physically, this point describes the transition between the stable and unstable regimes of the oscillator. is termed "dynamical localization" and is due to the quantum coherence between the different trajectories in momentum space.
Recently, the effect of quantum interferences in the dynamics of N interacting kicked rotors was investigated by mapping this model to a single quantum particle in a N -dimensional disordered lattice [64]. For N > 2, this model shows a transition between Anderson localized and delocalized phases, but the critical coupling tends to zero as N increases. Thus, in the thermodynamic limit of N → ∞, the effect of quantum coherence is negligible and one would expect the model to be always diffusive [65], as in the classical case [56,57].

Classical many-body kicked rotor
We consider a chain of N coupled kicked rotors, defined by the Hamiltonian This model is a natural extension of the kicked rotor (1), where the coupling between the rotors is induced by the kick §. Without the loss of generality, we assume that κ > 0: Models with positive and negative κ's can be mapped to each other through the transformation φ j → φ j + jπ for all j = 1, 2, 3, · · · . Using the classical Hamilton's equations of motion for the Hamiltonian in Eq. (9), one obtains a relation between the phase space coordinates just before the kick n and the same coordinates after one time period: − sin(φ j−1 (n) − φ j (n)) , φ j (n + 1) = φ j (n) + p j (n + 1)τ.
As mentioned in Sec. 2, one can find an unitless parameter K = κτ such that the dynamics depends on K only. A possible realization of the interacting kicked rotor model is based on the use of superconducting grains or molecules. In the first case, the Hamiltonian, Eq. (9), could be implemented using a chain of voltage-biased superconducting grains, in which the kick term is obtained by suddenly switching on a strong Josephson coupling for short times [66]. In the case of molecules, the Hamiltonian (9) could be implemented by periodic trains of laser pulses [67,68,69]. Alternatively, interacting kicked rotors can be realized using time-modulated potentials with non-commensurate wavelength [70].
The model (9) has been studied before in the context of chaos in almost-integrable models (The model is trivially integrable at K = 0) [71,72,73]. These earlier studies focused on two aspects of the problem: (i) the distribution of the Lyapunov exponents λ; (ii) the scaling of the diffusion coefficient D and of the largest λ max , with K. It was found that λ max ∼ √ K and D ∼ K α with α ≈ 6.5 [73]. This situation was named "fast Arnold diffusion", in contrast to the common Arnold diffusions, where § The model (9) is very convenient for numerical analysis, compared to a model with sinusoidal periodic driving. In this latter case the time evolution needs to be discretized in infinitesimal steps, preventing the calculation of physical observables at very long times. In contrast, for the kicked rotor model we can go up to very long times without facing such problem, thanks to the discrete nature of its time steps. We believe that the physical properties of the system will not change depending on the particular type of driving protocol. [74]. These studies focused on exponentially long times, where the system shows a diffusive behavior. In contrast, we here consider intermediate time scales, where the diffusion coefficient has not reached its long-time limit yet, and the dynamics still depends on the initial conditions.
In this context, our work bears some analogies with the celebrated Fermi-Pasta-Ulam-Tsingou (FPUT) model [75,76]. The FPUT problem originated from one of the first numerical simulations of a physical system on a digital computer and involved the dynamics of a classical many-body system under a non-linear timeindependent Hamiltonian. In contrast to the original expectations, the numerical solution performed by FPUT showed that the system was not ergodic, and did not lead to an equidistribution of the energy among the modes. This effect can be understood within the framework of Arnold diffusion. Because the non-linearities considered by FPUT were relatively small, the time required to observe diffusion were exponentially large and could not been observed within the computational resources available at that time.
In our numerical calculations, we consider an ensemble of initial conditions in which the initial momenta are zero, p j (0) = 0, and the initial positions are given by identical independent variable with a Gaussian distribution f (φ j ) = , with mean φ j = jφ 0 and standard deviation σ. The parameter σ sets the initial fluctuation of the system and will be an important tuning parameter of our theory. In this paper, we will focus on two specific cases of φ 0 , the case φ 0 = 0, where the system is initialized close to the classical ground state of the system (φ j = p j = 0), and the case φ 0 = π, where the nearest neighbor oscillators are initialized in opposite phases.
We characterize the dynamical phases of the model through the diffusion coefficient Here, t i and t f are two different stroboscopic times (t = nτ , where n is an integer), and where · · · represents the averaging over initial conditions. All the averages in this work, are done with 1000 different realizations of initial conditions (see also Appendix B for a case study of the error bars of our results). Because the average kinetic energy of the model is 1 2 p 2 , the coefficient D measures the average heating rate of the system. This coefficient is expected to converge to a finite value in a diffusive phase, while it tends to zero in a localized phase. To distinguish between diffusive and super-diffusive phases, we further calculate the exponent α by fitting the numerically computed p 2 (t) with a generic power-law function At α . The parameter α allows us to identify phases that are diffusive (α = 1), subdiffusive (α < 1), and superdiffusive (α > 1).
Our main aim is to determine the stability diagram of the model (9). To achieve this goal we numerically determine the diffusion parameters D and α as a function of K, for different choices of φ 0 and σ. These quantities are computed for different time ranges, where the initial time is fixed at t i = 1000 and the final time is varied up to t f = 6.5×10 4 . The choice of t i is dictated by the demand that the energy redistributes between the position (φ j ) and momentum (p j ) degrees of freedom (also see discussion in Appendix A). Throughout the paper, we consider systems of N = 1000 rotors. As shown in Appendix B, our numerical results are not sensitive to finite size effects.

Quadratic expansion and beyond
Before presenting our numerical results, let us first consider the linear stability of the model around φ 0 = 0. Using the quadratic approximation cos We can rewrite the Hamiltonian (9) as a collection of decoupled kicked harmonic oscillators where q = 2πk/N , k being an integer, is the wave number, and Here φ j e −iqj are the Fourier transforms of p j and φ j respectively. The Hamiltonian in Eq. (12) is a set of decoupled harmonic oscillators and is therefore integrable for any value of K. This Hamiltonian can nevertheless describe the transition between the localized and delocalized region, as the appearance of dynamical instabilities, where the phase space orbits change from elliptic to hyperbolic.
To determine the dynamical stability of Eq. (12), we proceed in analogy to Sec. 2 and first derive the associated equations of motion. Note that since p j and φ j are real functions, we have P * q = P −q and φ * q = φ −q . As a result, for each q, the classical equations of motion are given by Following a few steps similar to Sec. 2, we obtain the matrix M q for each q, which provides the time evolution of the system over one full period, and is given by Diagonalizing the matrix M q , we get two degenerate eigenvalues for each q: From Eq. (16), the condition of the transition between stable and unstable regions is given by F (q)τ = 4, or equivalently K sin 2 (q/2) = 1, for each q. As a consequence, the system is stable for K < 1: in this regime, all the eigenvalues are pure phases λ ± (q) = e ±iθ(q) with cos θ(q) = 1 − F (q)τ /2. For later reference, we rewrite the eigenvalues as λ ± (q) = e ±iµ(q)τ , where µ(q) = θ(q)/τ is the frequency of q-mode.
Although the system is linearly stable for K < 1, higher-order effects can lead to non-linear resonances. To study these effects, let us consider the 4 th order term in the Taylor expansion of Eq. (9) around φ 0 = 0: Note that the matrix Mq is real and of unit determinant as in the case of the single rotor, (Eq. (2)), although its size is 4 × 4, rather than 2 × 2.
where we rewrote the sum of delta functions (the kicks) as a Fourier series, and neglected the constant m = 0 term. The "unperturbed" Hamiltonian is given by Eq. (12) and describes a set of decoupled kicked harmonic oscillators. In this case, perturbation theory states that the system will absorb energy only if one of the driving frequencies equals the sum of the frequencies of 4 oscillators, or where µ(q) is the frequency of a q mode as defined above. Because µ(q) ≤ µ(π), the right-hand side of Eq. (18) is smaller or equal than 4µ(π) = 4 cos −1 (1 − 2K)/τ . Thus, the resonance condition (18) can be satisfied only if This condition can be satisfied only for m = 1, and predicts that the system is unstable to quartic couplings for K ≥ 0.5. If we consider higher-order terms in the Taylor expansion of the cosine, we can obtain resonances at lower values of K: In general, the k th -order term (∼ φ k ) gives rise to a series of resonances ¶ at K = 1/2(1−cos(2mπ/k)).

Linear stability and resonances
To address the linear stability regime numerically, we consider a system in which all the initial values of φ's are taken from a Gaussian ensemble with φ 0 = 0. First, we consider small initial fluctuations: the standard deviation of the Gaussian ensemble is chosen to be σ = 0.1. The numerically evaluated diffusion coefficient D and exponent α are shown in Fig. 1. In both plots, we find two kinks, respectively at K = 0.5 and K = 1. The positions of these two points do not significantly scale as a function of time.
We can characterize the different regions in terms of the behavior of D and α as a function of the waiting time ∆t. Recall that D and α are calculated by fitting p 2 (t) with linear and power-law function of time, respectively. We thus have the relations p 2 (t) = Dt and p 2 (t) = At α , that lead to D(t) = At α−1 . Therefore it is clear that D does not depend on time for α = 1. As shown in Fig. 1(a), for K > 1 the diffusion coefficient D has a non-zero value and all the curves for different time ranges merge to each other. Accordingly, as shown in Fig. 1(b), the exponent α converges to 1 for all the curves. Both results indicate that in this region the system is diffusive. On the other hand, for K < 0.5, both D and α acquire zero value, and the system is localized. To highlight the behavior of D and α at K = 0.5, we show the same curves around K = 0.5 with logarithmic scale along vertical axis. From Figs. 1(c,d), we observe that both D and α jump from nearly zero values to non-zero one for all the curves at K = 0.5. For 0.5 < K < 1, the exponent α decreases with waiting time with two kinks at the two ends. The kink at the point K = 1 becomes smoother as a function of the waiting time. In contrast, the transition to the localized regime at K = 0.5 is signaled by a kink in α that becomes sharper as a function of time.
The numerically determined stability diagram can be understood within the linear analysis presented in Sec. 4. As shown there, the points K = 1 and K = 0.5 correspond respectively to the stability thresholds of the 2 nd and 4 th order couplings. For K > 1, ¶ As a sanity check of our approach one can verify that second order terms ∼ (φ i − φ i+1 ) 2 give rise to a resonance at K = m, in agreement with the linear stability analysis. the system is linearly unstable and shows a diffusive behavior at all times. In contrast, for K < 0.5, the system is stable to both quadratic and quartic perturbations and is thus localized for all the numerically achievable times. In the intermediate regime (0.5 < K < 1) it is super-diffusive for small waiting times, but becomes diffusive for longer waiting times. As explained in Sec. 4, higher order terms can give rise to resonances at lower values of K. (For example the k = 6 coupling becomes unstable at K = 0.25). To observe these resonances, we repeat the same calculations with larger values of the initial fluctuations σ (at a fixed waiting time). For small σ, the quartic resonance at K = 0.5 is clearly visible as a kink in the diffusion exponent α (Fig. 2(b)) and also a sudden jump in the diffusion coefficient ( Fig. 2(a)). In contrast, for large initial fluctuations (σ > 0.3), the higher order resonances come into play, and α becomes a smooth functions of K. We can find a more clear picture how the jump in D and α at K = 0.5 become smoother with increasing initial fluctuations in Figs. 2(c) and (d), respectively. As we now explain, this situation actually reveals a different type of stability.

Marginal localization
To investigate the continuous behavior of the diffusion parameters, we now focus on σ = 1.0 for different waiting times. In Fig. 3(a), we show that the value of D decreases with time for small K and there is a transition point after which D starts increasing with time. For K 1, all the curves merge to each other. To better characterize the transition, we apply a scaling transformation that leads to a data collapse for small K. As shown in Fig. 3(b), this effect can be achieved by rescaling the x and y axes by log(t f ) and t f respectively. We observe that for K 3.0/ log t f the curves merge to each other, indicating that D is proportional to 1/t f . According to the definition of D in Eq. (11), this finding indicates that the kinetic energy of the system does not grow with time for K 3.0/ log t f , or equivalently that the system is localized (in momentum space) and does not heat up.
The localized region can be alternatively identified by inspecting the graph of α as a function of K for different ranges of time. As shown in Fig. 4(a), these plots show a pronounced peaks, whose position moves to lower K as t f is increased. If we rescale the x-axis by log t f , we observe a collapse of all the curves (see Fig. 4(b)). As expected, the asymptotic behavior of these curves changes at K c ≈ 3.0/ log t f . For K < K c , α decreases as a function of t f , confirming that the system is localized. On the other hand, the system is super-diffusive and diffusive for 3.0 K log t f 6.0 and   K 6.0/ log t f , respectively. Interestingly, at K ≈ 3.0/ log(t f ) the diffusion exponent α shows a steep increase, which sharpens as a function of t f . We conjecture that at very long time, the model transitions directly from the localized regime (α = 0) to the super-diffusive one (α > 1) without showing any sub-diffusive behavior. As a key result, we find that the the localized region survives until t f = exp(3.0/K) which is exponentially long for small K = κT . We therefore refer this case as marginal localization. This finding is in line with earlier studies of quantum systems [41,42,53,45,43], who showed that the thermalization time is exponentially long in the ratio between the driving amplitude and the driving frequency. A similar behavior was observed in time-independent classical systems close to an integrable fixed point as well [77].
The scenario of marginal localization is also achieved by considering initial values of φ's close to the unstable point, φ 0 = π, even when the initial fluctuations are small (σ = 0.1). The corresponding results for D and α are shown in Figs. 5 and 6, respectively. The behavior of both the observables is similar to the previous case of marginal localization, although, the ranges of the different regimes are bit different compared to the previous case. In this case the system is localized for  Fig. 3(b), all the D curves collapse for K 1.5/ log(t f ).  K 1.5/ log t f , whereas it becomes super-diffusive and diffusive for 1.5 K log t f 6 and K 6/ log t f respectively. The transition points differ from the previous case, indicating that the exact location of the phase boundaries of the model depend on the specific choice of the initial conditions, although the scaling behavior is unchanged.

Wave decomposition
The localized and delocalized regimes differ not only in the time evolution of p 2 , but also in the energy distribution among the Fourier components. To show this effect, we compute |P q | 2 as a function of q inside the different phases. Here P q is the wave component of the momentum of the rotors and is formally defined in Sec. 4.
As before, we divide our discussion between initial conditions close to the stable point (φ 0 = 0), and to the unstable one (φ 0 = π). The results of the former case are shown in Fig. 7, where we numerically calculate |P q | 2 as a function of q for K = 0.2, 0.7 and 1.4, which correspond to localized, super-diffusive and diffusive phases, respectively (see Fig. 1). In all the plots the momentum modes close to q = ±π are most excited and the q = 0 mode remains unexcited. The Hamiltonian in Eq. (9) is invariant under the translation transformation j → j + 1 and, therefore, the total momentum of the system is a conserved quantity. Since all the initial momenta of the rotors are chosen to be zero, the momentum at q = 0 remains zero even for later time. In addition, we have also the relation P −q = P * q which makes the spectrum symmetric about q = 0.
Let us first consider the diffusive case (K = 1.4). In this regime, p 2 = q |P 2 q | grows linearly with time. Fig. 7(c) shows that each q mode grows independently according to |P 2 q | ∼ t sin 2 (q/2), where t is the stroboscopic time. Therefore, from the definition of |P q | 2 , we obtain a q-dependent diffusion coefficient D(q) ∝ sin 2 (q/2). We repeated the same analysis with initial angles of rotors chosen from another random distribution, but the behavior of the energy per mode remains intact indicating the form of D(q) is universal in this parameter regime. Note that in a thermal state, P 2 q does not depend on q. Thus, our findings indicate that the system does not thermalize, even at asymptotically long times. This happens due to the absence of upper bound of energy in the spectrum. The kinetic energy of the system increases unboundedly with time, which indicates that it acquires infinite value at infinite temperature ensemble. As a result, the thermalization can not occur at any finite time.
We next consider the localized regime (K = 0.2). As shown in Fig. 7(a), at a given time |P q | 2 randomly oscillates as a function of q. The range of the oscillations does not depend on time. As discussed in Sec. 4, the Hamiltonian in Eq. (9) can be expanded around the stable point φ 0 = 0 and it provides a collection of decoupled harmonic oscillators in the Fourier space. Therefore, for this case, the energy oscillates between the position and momentum of the q mode, φ q and P q , for each q, but the total energy remains constant with time. This phenomenon is reflected in the oscillations of |P 2 q | as a function of q at a fixed time. As expected from the linear model, we do not observe energy exchange among the wave modes.
In the intermediate regime, at K = 0.7, the average energy per mode shows oscillations for shorter times, while the oscillations die out as time increases (see Fig. 7(b)). At long times, the modes around q = π feel the 4 th order resonance described in Sec. 4, and the shape of |P 2 q | starts deviating from the simple sin 2 (q/2) shape. This effect is associated with the super-diffusive behavior described in Sec. 5.
We now discuss a system that shows marginal localization, by considering initial conditions close to the unstable point (φ 0 = π). Following the same strategy as before, we numerically calculate the energy distribution for the localized, super-diffusive and diffusive regimes (see Fig. 8). As shown in Fig. 8(c), the behavior of |P q | 2 in the diffusive regime is very similar to the previous case.
In the marginally localized regime (K = 0.05), the wave modes close to ±π get de-excited, whereas the modes near q = 0 become excited as the number of kicks is increased. This indicates that the excited modes transfer energy to relatively less excited modes, but the system does not absorb energy from the driving. As a result, p 2 does not grow with time in this regime leading to vanishing values of D and α (see Figs. 5 and 6). By extrapolating Fig. 8(a) to asymptotically long times, one finds that |P 2 q | tends to a constant function (with a singularity at q = 0). This finding implies that, although, the system eventually becomes diffusive at exponentially long time (see discussion in Sec. 6), it prethermalizes to a finite effective temperature at intermediate time + This finding is in agreement with recent studies of quantum cases [40,49,50,51,52,53,54], where the authors showed that the driven quantum systems can prethermalize to a thermal steady state at some intermediate time scale.
The intermediate regime (K = 0.25) is shown in Fig. 8(b). Its behavior is also intermediate between the localized and diffusive phases. For short times, the system is in the localized regime and transfers energy from large q to small q. In contrast, for longer times the modes close to q = π become excited, leading to a super-diffusive dynamics.

Summary and conclusions
In this paper, we revisited the problem of stability and thermalization of periodically driven classical many-body system. In quantum systems, the eigenstate thermalization hypothesis (ETH) states that the system should thermalize at long times, but it does not make any prediction about the dynamics of the system at intermediate time scales. For instance, earlier studies showed that the system can reach an intermediate longlived pre-thermalized state before approaching to the infinite temperature ensemble. We addressed the classical analogue of this question by investigating the finitetime dynamics of a classical many-body kicked rotor, whose dynamics is given by a generalization of the "standard map" (see Eq. (10)).
When the initial angles of the rotors are chosen in vicinity of the stable point, we observed two sharp transitions, whose position does not scale with the waiting time. The first transition (K = 1) is observed only at short times, and corresponds to an exceptional point of the linear dynamics, where the fixed point becomes unstable. The second transition (K = 0.5) becomes sharper as a function of time, and is due to quartic terms, which can lead to energy absorption in the linear stability regime. The small initial fluctuations suppress the effects of higher-order terms, and lead to a system that is stable for all numerically achievable times: In this regime the system is nearly integrable and never equilibrates [31,32,33,34,35,36,37,38].
For large initial fluctuations, the linear analysis breaks down and both the diffusion coefficient D and exponent α become continuous functions of K = κτ . A similar situation appears when the initial angles of the rotors are considered around the unstable point. For both the cases, we observe a collapse of the α curves for different waiting times by rescaling the x-axis with log t f . This finding indicates that the localized phase of the system for small K sustains for an exponentially long time and the system is marginally localized. In other words, the rate of approach of the system to the diffusive phase is logarithmically small predicting the existence of a long-lived prethermal regime.
To analyze the phenomenon of prethermalization, and to study the localized and delocalized regimes further, we also considered the Fourier space representation of + The phenomenon of prethermalization is also observed in spinless fermion models with a weak integrability-breaking term after a quantum quench [78]. In this case, the prethermalization time shows a power law variation with interaction term, which is in contrast to our observation of exponentially long localization time for a periodically driven model.

Initial conditions Localized Diffusive Figures
Linear stability Marginal localization the momenta of the rotors. In this connection, we numerically calculated the average fluctuations of the wave components in the different regimes. In the diffusive phase the system absorbs energy from the external driving and the average momentum grows linearly in time. Interestingly, the kinetic energy does not redistribute evenly among the wave-vectors, but remains proportional to sin 2 (q/2) at all times. This observation indicates that although the system is diffusive, it does not thermalize. The reason for this anomaly is found in the un-bound nature of the system's spectrum: Because the kinetic energy is proportional to p 2 , its expectation value should be infinite in an infinite temperature ensemble. Clearly, this cannot happen at any finite time.
In the marginally localized regime, the average energy of the wave modes at small q increases with time, while it decreases for modes close to q = ±π. This indicates that energy is transfered from excited to de-excited modes, but the system does not absorb energy from external driving. The expectation value of P 2 q tends to a constant (with a sharp dip at q = 0), suggesting that the system is flowing towards a prethermalized state, whose energy is set by the fluctuations of the initial state. At exponentially long times, the system escapes from this prethermal state and becomes diffusive. The transition between the marginal and diffusive regimes occurs through an intermediate super-diffusive regime, where energy absorption and thermalization compete.
Our analysis highlights two alternative methods to stabilize periodically driven many-body systems: (i) linear stability, induced by the choice of initial conditions close to a fixed point of the classical dynamics, and (ii) marginal localization, obtained for weak and/or fast drives (see Table I). In both cases, the system can remain stable for an arbitrarily long time: in the former case the instability is suppressed by the closeness to the fixed point, while in the latter it is suppressed (exponentially) by the ratio between the driving frequency and the driving amplitude. Our analysis clarifies the nature of many-body instabilities observed in the earlier literature. In particular, Refs. [28,29,30] considered an initial state with (quasi) long range correlations, which can stabilize the system in the linear stability mechanism. In contrast, the prethermal states observed by Refs. [52,41,45,53,48,38,51] are protected by a large driving frequency, and thus fall into the category of marginal localization. These earlier works considered many-body quantum systems, but the effects that they describe exist in classical systems as well.
discussions. This work was supported by the Israeli Science Foundation Grant No. 1542/14.

Appendix A. Variation of p 2 with time
In this appendix we show typical plots of the time evolution of p 2 , which were used to extract the diffusion parameters D and α. Fig. A1(a) refers to the time evolution in the linearly localized regime. On the other hand, Figs. A1 (b) and (c) exhibit the time evolution for the two scenarios of marginally localized regime (see Sec. 6).
In all cases, the time evolution initially shows fast oscillations, associated with the equipartition of energy between the φ and p degrees of freedom. For Figs. A1 (b) and (c), these oscillations are followed by a smooth evolution in time, which is the focus of the present article. In contrast, the oscillations in Fig. A1(a) sustain for very long time, since the initial conditions are chosen close to the stable point (φ 0 = 0). However, the amplitude of oscillations decreases in time. For this reason, we extract D and α from the time evolution of p 2 at times t > t i = 1000 only. We checked that our results are independent on the specific choice of t i . Appendix B. Effect of finite system size and error due to random initial conditions In this Appendix we investigate whether our numerical results have finite size effects.
In addition, we also analyze the error in our data considering the plots as shown in Figs. 5 and 6. For the first check we calculate both D and α for three different system sizes with one waiting time. As shown in Fig. B1, all the curves for different system sizes perfectly overlap each other exhibiting that our results are free from finite system size effects. However there exists the effect of finite time for which we have provided a clear picture in the main text considering different times and initial conditions. We now consider the scenario of marginal localization as shown in Figs. 5 and 6 as an example to show statistical error in our data due to random initial conditions. Fig. B2 shows the statistical error in diffusion constant D, determined by Eq. (11) with 1000 random initial conditions. On the other hand, the error in α is given by the error in fitting p 2 (t) as a general power-law function At α . (c) and (d) are the same plots as in Fig. 6 with statistical error in the data points indicated by the symbol sizes.