Pre-thermalization in a classical phonon field: slow relaxation of the number of phonons

phonons François Huveneers and Jani Lukkarinen Université Paris Dauphine-PSL, CEREMADE, Place du Maréchal de Tassigny, 75016 Paris, France∗ University of Helsinki,Department of Mathematics and Statistics, P.O. Box 68, FI-00014 Helsingin yliopisto, Finland† (Dated: February 26, 2020) We investigate the emergence of an astonishingly long pre-thermal plateau in a classical phonon field, here a harmonic chain with on-site pinning. Integrability is broken by a weak anharmonic on-site potential with strength λ. In the small λ limit, the approach to equilibrium of a translation invariant initial state is described by kinetic theory. However, when the phonon band becomes narrow, we find that the (non-conserved) number of phonons relaxes on much longer time scales than kinetic. We establish rigorous bounds on the relaxation time, and develop a theory that yields exact predictions for the dissipation rate in the limit λ→ 0. We compare the theoretical predictions with data from molecular dynamics simulations and find good agreement. Our work shows how classical systems may exhibit phenomena which at the first glance appear to require quantization.

In this letter, we investigate a classical many-body Hamiltonian H = H 0 + λV in the limit λ → 0, where H 0 is integrable (a harmonic lattice or free phonon field) and V breaks integrability (an on-site anharmonic potential). If the system is started in a translation invariant state, its state evolves swiftly to the generalized Gibbs ensemble (GGE), characterized by all the conserved quantities of H 0 [26][27][28][29]. Next, according to kinetic theory and the Boltzmann-Peierls equation, it approaches equilibrium in a time τ 1 with τ 1 ∼ λ −2 for λ → 0 [30][31][32][33]. However, if the phonon band is sufficiently narrow, the (non-conserved) number of phonons is preserved by kinetic processes [32,34], and only a pre-thermal plateau is reached on kinetic time scales. As our analysis reveals, the proper equilibrium is only reached after a longer time τ 2 , scaling as τ 2 ∼ λ −2p for some p ≥ 1. See Fig. 1 for a summary of the above process.
The presence of an almost conserved quantity (or adiabatic invariant) for the Hamiltonian H studied in this letter, has been realized in [18,19]. Here, we first provide rigorous quantitative bounds on the dissipation of the number of phonons, see Claims 1 and 2 below. In addition, we connect them to the slow dissipation of some quantized fields [13,15], a phenomenon that seemed to require quantization. Second, we provide a theory to compute the dissipation rate of the pseudo-conserved quantity, which we are able to back with numerical results. General predictions for the rate have been derived in [35] for quantum systems, see also [36][37][38]. However, we notice that our system is classical and that an extra time-scale is present since relaxation to the pre-thermal plateau involves kinetic processes. Model -Let the Hamiltonian H be given by with r > 2 an even integer (below we focus on r = 4, 6), and ω 0 > 0 a characteristic frequency of the system. The dynamics is classical:q x = p x anḋ p x = −(1 − 2δ)ω 2 0 q x + δω 2 0 (q x+1 − 2q x + q x ) − λq r−1 x whereẋ denotes the time derivative of x. Stability imposes λ ≥ 0 and 0 ≤ δ ≤ 0.5. The chain is uncoupled for δ = 0 and unpinned for δ = 0.5, and we restrict our attention to 0 < δ < 0.5. If needed, we may obviously restrict the summation in eq. (1) to 1 ≤ x ≤ L for a length L, and consider the limit L → ∞ only afterwards. For λ = 0, the chain is harmonic. For a pseudomomentum k in the Brillouin zone BZ = ]−0.5, 0.5], let the phonon mode a(k) := a − (k) be defined by a ± (k) = 1 √ 2 ω 1/2 (k)q(k) ∓ i 1 ω 1/2 (k)p (k) with the dispersion relation and withf the Fourier transform of f , defined byf (k) = x∈Z f (x)e −2iπkx . We identify n(k) = |a(k)| 2 with the number of phonons with pseudo-momentum k and we denote the total number of phonons by N 0 = BZ dk n(k). From the analyticity of ω(k) in eq. (3), we deduce that N 0 is quasi-local, i.e.
where the kernels K (z) decay exponentially with |z|.
Pseudo-conservation of N 0 -To compute the timescales on which N 0 gets dissipated, let us first assume that the system is quantized, i.e. that a ± are creation/annihilation operators for bososns. Later on, we will see how the conclusions carry over to the classical system. In this case, N 0 has integer spectrum. In the limit λ → 0, only resonant processes, preserving the bare energy H 0 , do effectively destroy the conservation of N 0 . Therefore, in first order in λ, the process of creating two phonons (since r is even, it is not possible to create only one phonon) must satisfy the constraints ω 1 + · · · + ω n/2+1 − ω n/2+2 − · · · − ω n = 0, where n = r (we will later consider n = r when dealing with higher order processes), and where the second constraint is taken modulo 1 and stems from translation invariance, cfr. eq. (4). Since the width of the dispersion relation in (3) scales as 2ω 0 δ for small δ, larger and larger values of r are needed to satisfy the constraints in eq. (5). Thus, for given n ≥ 2 even, there exists δ c (n) such that eq. (5) only has solutions for δ ≥ δ c (n): δ c (2) = δ c (4) = 0.5 (exact), δ c (6) = 0.3 (exact), δ c (8) 0.25 (numerical), and in general δ c (n) ∼ 2/n for n → ∞, see the Supplemental Material (SM). From now on, we will assume that δ is such that δ = δ c (n) for any n, leaving these exceptional cases for further studies.
Higher order processes need obviously to be taken into account in the case δ < δ c (r). The analysis detailed in the SM yields that a process of order λ p involves the creation/annihilation of n = p(r−2)+2 phonons[39], and the above analysis caries over with this new value for n. Given r and δ, we may now determine the smallest p so that a process of order λ p destroys N 0 effectively: p ≥ 1 is the only integer such that Explicit results are gathered on Table I: Table I. Order of the processes (λ p ) destroying effectively the conservation of N0, for r = 4, 6 in eq. (4).
We observe that, even though quantitative statements are obviously model dependent, and in particular the threshold values δ c depend on the specific form of the dispersion relation ω(k) in eq. (3), the conclusion that p ∼ c/δ as δ → 0 is generic (for a polynomial interaction V ), following from the fact that the width of the band in ω(k) decays as δ for δ → 0.
The above analysis may be turned into rigorous results, using the formalism developed in [13] that can be straightforwardly adapted to a classical system through the canonical replacement −i[H, ·] by {H, ·}, see SM and below. The key observation to proceed is that, even though N 0 is no longer quantized, the spectrum of ad N0 = {N 0 , ·} is: Hence, it acts formally in the same way as the quantum super-operator −i[N 0 , ·], and this is eventually what matters. Let us fix r and δ such that p > 1, and let us first express a result in a formulation directly inspired from [13]: There exists a canonical change of variables, bringing H intoH, so that {H, N 0 } = O(λ p ), where the term O(λ p ) is extensive. To formulate a precise claim, we will consider extensive observables ϕ m of the form where m ≥ 2 and whereφ m is analytic in (k 1 , . . . , k m ), ensuring that ϕ m is a sum of quasi-local observables. A function ϕ will simply be said to be a polynomial (of order m) if it is of the form ϕ = m k=2 ϕ k with ϕ k as in eq. (7). The following claim is shown in the SM: Claim 1 Let L be finite and let us assume periodic boundary conditions. For |λ| small enough, there exists a polynomial G = λ p−1 n=1 λ n−1 G n , such thatH := e −{G,·} H is a well defined real analytic function in a neighborhood of the origin in R 2L , and such that {H, N 0 } = λ p ∞ n=p λ n−p J n , where J n are polynomials and where the expansion converges to an analytic function in a neighborhood of the origin in R 2L .
Claim 1 provides a good way to think about the phenomenon but does not yield as such a very powerful result in the thermodynamic limit: The radius of convergence of λ may shrink as L → ∞ (even though we are interested in the regime λ → 0, the limit L → ∞ needs to be taken before the limit λ → 0). As a way out, we may undo the above transformation to obtain a dressed number of phonons, e {G,·} N 0 , and then truncate its expansion at order p − 1. Doing so, we end up with a well defined pseudo-conserved quantity N in the thermodynamic limit, see the SM for details: Claim 2 Let now the chain be defined on the full lattice Z. There exists a quantity N = N 0 + λ p−1 n=1 λ n−1 N n , where N n are polynomials of order nr − 2(n − 1) for 1 ≤ n ≤ p − 1, such that {H, N } = λ p {V, N p−1 }.
These claims furnish upper bounds on the dissipation rate of N , but the determination of this rate requires the knowledge of the instantaneous state of the system. To proceed further, we will invoke additional assumptions, and leave mathematical rigor behind.
Evaluation of the decay rate -For p = 1, N = N 0 . The dissipation rate γ of the density N 0 /L in an instantaneous state ρ and in the infinite volume limit is given by γ = J ρ /δ(0), with a flux J = λJ = λ{V, N 0 } and where "δ(0)" stands for the infinite volume, corresponding to L in a chain of finite length. If the system is prepared in a translation invariant state with zero average, after a short transient time, it evolves towards a GGE characterized by a Wigner function W , i.e. a Gaussian state e − BZ dkn(k)/W (k) /Z. Usual kinetic theory yields the following expression for the rate γ(W ): Given a function ϕ = BZ dkφ(k)n(k), let J ϕ = λJ ϕ = {H, ϕ} = λ{V, ϕ}, then where · W denotes the average over the GGE, and where the dynamics in the time integral is the free dynamics (λ = 0). Expression (8) can thus be evaluated explicitly in leading order. In the SM, we provide a derivation of eq. (8) which is fully consistent with the derivation used in the more general case p > 1. We tested numerically the validity of eq. (8) for r = 6 and δ = 0.35, for an out-ofequilibrium initial state corresponding to W = (β(ω(k)− µ)) −1 with µ = −1 and various values of β. We found excellent agreement with the value of γ(W ) extracted from direct simulation of the dynamics, see Fig. 2. See SM for the numerical protocol. For p > 1, we find it convenient to move to the rotated frame where N 0 is pseudo-conserved quantity ofH, see Claim 1. The GGE is now parametrized by only two extensive quantitiesH and N 0 . Our main assumption (to be discussed later on) is that the system is in a state ρ, λ p -close to the GGE ρ 0 : where f is some correction that will be determined by maximizing local stationarity. The key observation is and the integral of a Poisson bracket vanishes. Hence the instantaneous dissipation rate γ(β, µ) may be written as Since this rate scales as λ 2p , the state ρ itself evolve on that time scale. Stationarity on shorter time scales determines f and an explicit computation yields, see SM: Again, the dynamics in the time integral is the dynamics generated by H 0 and γ(β, µ) can thus be computed explicitly in leading order. We performed two sets of tests in the case p = 2 (accessing larger values of p would require too long simulation times). In all cases, we start from a state of the type ρ 0 with µ = −1 and various values of β. For r = 4 and δ = 0.45, the results reported on Fig. 3 show very good agreement between the prediction from eq. (10) and direct simulation of the dynamics. For r = 6 and δ = 0.28, the observed rate is significantly smaller than the one predicted by eq. (10), but it decreases slower as a function of λ and β −1 , see Fig. 4. Comparing with the discrepancies at the largest values of λ on the upper panel of Fig. 3, makes it plausible that our theory just needs smaller values of λ to be validated. Besides, the fact that the observations are below the theoretical predictions is a second indication that the theory will become accurate for smaller values of λ, since a smaller rate guarantees that our main hypothesis, eq. (9), from which our predictions follow, is more easily satisfied.
Irrespectively of numerical observations, we finally would like to make a consistency check of the main assumption in eq. (9). On the one hand, due to the dissipation of N 0 , the state ρ 0 evolves with time at a rate v 1 ∼ λ 2p , i.e. β, µ evolve at this rate in order to yield correct values for H ρ0 and N 0 ρ0 . On the other hand, the system relaxes towards the instantaneous pseudoequilibrium ρ 0 through kinetic processes. Assuming that the state ρ is at a distance of order λ p from ρ 0 , as required by eq. (9), we conclude that it moves at a rate v 2 ∼ λ 2 × λ p . Consistency of the theory requires that v 1 v 2 , i.e. that the instantaneous fixed point ρ 0 moves slow enough so that the state ρ has the time to relax to it. Clearly, this is wrong for p = 1, marginal for p = 2 and fine for p > 2. We had treated separately the case p = 1, since indeed there is no reason to think that the pre-thermal state should be characterized by the two parameters β, µ only. Unfortunately, the above argument is not conclusive for p = 2, while numerical data are only available in this case.
Conclusions and outlook -Our work provides a new example of long pre-thermal plateau, it shows how a phenomenology initially explored in quantum systems carries over to a classical set-up, and it participates to recent efforts to describe accurately the dissipation of pseudoconserved quantities. The main features of our theory carry over to d > 1 and, for d = 3, we may contemplate the possibility of realizing a pre-thermal Bose-Einstein condensate in this classical system, exploiting the conservation of the number of phonons over a very long period.
We thank W. De Roeck and H. Spohn for helpful discussions, and C. Mendl for providing the original code for numerical simulations. F. H. and J. L. benefited from the support of the project EDNHS ANR-14-CE25-0011, and F. H. from the project LSD ANR-15-CE40-0020-01 of the French National Research Agency (ANR), as well as from the support of the International Centre for Theoretical Sciences (ICTS) during a visit for the program -Thermalization, Many body localization and Hydrodynamics (Code: ICTS/hydrodynamics2019/11). The work has also been supported by the Academy of Finland via the Centre of Excellence in Analysis and Dynamics Research (project 307333) and the Matter and Materials Profi4 university profiling action. * huveneers@ceremade.dauphine.fr † jani.lukkarinen@helsinki.fi  Here we derive bounds on the quantity δ c (n), with n ≥ 2 even, such that eq. (5) only admits solutions for δ ≥ δ c (n). The nearest neighbor dispersion relation is given by eq. (3), i.e. ω(k) := ω 0 1 − 2δ cos (2πk) with ω 0 > 0, 0 < δ ≤ 1 2 . Eq. (5) admits a solution if the number non-preserving collisional manifold is not empty, i.e. if there is k ∈ (BZ) n and σ ∈ {±1} n for which with σ such that ∆N 0 := n =1 σ = 0 (we consider only processes that do not preserve the number of phonons). To clearly make the connection with eq. (5), we notice that by sign-change symmetry, we may focus on the case with ∆N 0 > 0, permute the labels so that all positive signs come before the negative ones, and remark that eq. (11) admits a solution for ∆N 0 > 0 if and only if it admits a solution for ∆N 0 = 2. This last point follows from the fact that Ω(k, σ) ≥ Ω(k, σ ) for any k ∈ (BZ) n , if σ ≥ σ for all 1 ≤ ≤ n, and from the fact that eq. (11) admits no solution if and only if Ω(k, σ) does not change sign, i.e. remains strictly positive/negative, on the set {k ∈ (BZ) n : n =1 k = 0} (BZ) n−1 (and the sign is the sign of ∆N 0 since Ω(0, σ) is proportional to ∆N 0 ). Below, for simplicity, we set ω 0 = 1 since its value will not affect the value of δ c .
Denote the minimum of ω(k) by m − and maximum by m + . The minimum is reached at k = 0 and the maximum at k = 1 2 and thus Denote n ± = #{ : σ = ±1} for which obviously n = n + + n − and ∆N 0 = n + − n − ≥ 2. We then have Since here n ± = (n ± ∆N 0 )/2, it follows that Consider then the following function of δ ∈ ]0, 1 2 ]: Clearly, Here G is strictly decreasing from +∞ to 2, and thus there are unique values δ − (N ) obtained as a solutions of G(δ − (n)) = n.
For n = 8, the above bounds yields 0.235 < δ c (8) < 0.255, and further numerical checks of the values of Ω on the values satisfying the translation invariance constraint show that δ c (8) ≈ 0.25.

PROOF OF CLAIM 1 AND CLAIM 2
We provide here a rigorous proof of Claim 1 and Claim 2. Let us first deal with Claim 1. Let r and δ be fixed such that p > 1. Let L be the length of the chain, and let us assume that the Hamiltonian in eq. (1) is defined with periodic boundary conditions. Eq. (2-4) still make sense, provided that we define BZ dkφ(k) = 1 L L−1 k=0φ (k/L) witĥ ϕ(k) = L x=1 ϕ(x)e −2iπkx , and δ(k/L) = L for k = 0 and δ(k/L) = 0 otherwise. Let finally the Poisson bracket for two functions f, g on the phase space R 2L be defined as Let us first perform formal computations that we will justify afterwards. Given a function −G = p−1 n=1 λ n G n on the phase space, we can expand the operator e −ad G = n≥0 λ n S n with S 0 = Id and For n ≥ 1, we further decompose S n as S n = ad Gn + T n−1 and we notice that T n only involves coefficients G k with k ≤ n. Hence For m ≥ 2, let us consider functions as in eq. (7), i.e. translation invariant homogeneous polynomials (TIHP) of order n:  (18). TIPs can be decomposed accordingly. The crucial property implied by this decomposition is that ad N0 (ϕ ) = 0. Assuming that our computations involve only TIPs, and we will show below that this assumption is legitimate, we can find the coefficients G n so that {e −ad G H , N 0 } = O(λ p ). For this, we require that G n solve the set of recursive equations Indeed, inserting (19) in (17) yields and thus To show that the above scheme make sense, and establish Claim 1, we need to prove that the equations (19) can be solved, and that the expansion (20) converges for |λ| small enough. Let us start with eq. (19). From the definitions (15) and (2), we derive the canonical commutation rule Together with the rule {f, gh} = {f, g}h + g{f, h}, we can readily evaluate the Poisson bracket between TIHPs. In particular we derive that if ϕ n1 is a TIHP of order n 1 and if ϕ n2 is a TIHP of order n 2 , then {ϕ n1 , ϕ n2 } is a TIHP of order n 1 + n 2 − 2. Moreover, if ϕ n is a TIHP of order n with kernelφ n (k 1 , . . . , k n , σ 1 , . . . , σ n ), then {H 0 , ϕ n } is again a TIHP of order n with kernel Hence, if ϕ n is a TIHP of order n, and if δ < δ c (n), ensuring that eq. (5) has no solution, then the equation admits a solution u given by with the convention 0/0 = 0. If we first do not pay attention to the regularity of the kernels involved, i.e. if we ignore possible singularities stemming from the fact that n j=1 σ j ω j may vanish in eq. (23), we find that G n solving eq. (19) are TIPs of order nr − 2(n − 1). To show next that singularities do not occur and that the kernels are analytic, we use that δ and p satisfy eq. (6). This guarantees in particular that δ < δ c ((p − 1)(r − 2) + 2) ≤ δ c (nr − 2(n − 1)) for all 1 ≤ n ≤ p − 1, and therefore analyticity. Finally, even though this is not needed for the proof, we notice also that we do not expect to be able to find a regular function G p solving eq. (19), since δ > δ c (p(r − 2) + 2) by eq. (6), i.e. we expect to have reached the optimal order p.
We next deal with the convergence of the expansion in (20). Let us consider the Hamiltonian H on the extended space R 2L+1 , so as to explicitly include the dependence of H on λ, and let us consider the Cauchy problem We observe thatH(1, ·) = e −ad G H, if both terms make sense. By the Cauchy-Kowalevski theorem, eq. (24) admits a real analytic solution in the neighborhood of the origin in R × R 2L+1 . Moreover, since G = O(λ), we may assume that the solution is well defined up to τ = 1 by shrinking the neighborhood in λ. This ensures the convergence of eq. (20). Let us finally move to Claim 2. We notice that e ad G {e −ad G H, N 0 } = {H, e ad G N 0 }, where both sides of the equality are well defined and analytic in λ in a neighborhood of the origin, by a similar argument as before. Moreover, {e −ad G H, N 0 } = O(λ p ) by our construction, hence also e ad G {e −ad G H, N 0 } = O(λ p ). Writing e ad G N 0 = n≥0 λ n N n and defining N = p−1 n=0 λ n N n , we conclude that {H, N } = λ p {N p−1 , V }. The quantity N defines an extensive quantity in the thermodynamic limit L → ∞, and the last relation remains true in this limit. This yields thus Claim 2.
Eq. (8): p = 1 -Let a Wigner function W be given, and let us assume that the systems is in a state of the form where λf represents a first order correction. This assumption is the analog of the assumption (9) that will be used for the case p > 1. If ϕ is any observable of the type ϕ = BZ dkφ(k)n(k), its flux J ϕ = {H, ϕ} = λ{V, ϕ} = λJ ϕ vanishes on average in the state ρ 0 : J ϕ ρ0 = 0. Therefore, the occupations of all phonon modes n(k) must evolve on time scales of order λ 2 , and the whole state ρ evolves thus only on these time scales. Expressing this mathematically determines the first order correction f : where ad † H is the adjoint of ad H with respect to the measure ρ 0 . Let us compute the adjoint ad † H . It is defined as the operator such that ρ 0 (ad H u)v = ρ 0 u(ad † H v) for any functions u, v. We compute Combining eq. (25) and eq. (26), we find that f satisfies {H 0 , f } = J 1/W . However, due to resonances, i.e. due to the fact that 0 is in the spectrum of ad H0 , we insert a regularization to solve this equation: Given τ < +∞, we consider instead the equation (ad H0 + 1 τ )f = g and consider the limit τ → ∞. This can be solved as where g(t) is the evolution of g for the free dynamics generated by ad H0 . We can now derive eq. (8). Let J = λJ = λ{V, N 0 }. We compute Eq. (10): p > 1 -We proceed in a very similar way. As explained in the main text, we find it convenient to move to the rotated frame where N 0 is a pseudo-conserved quantity for a dressed HamiltonianH. According to eq. (9), we assume that the system is in a state ρ of the form where λ p f represents the correction at order p. As derived in the main text, the flux of N 0 /L, i.e. J = {H, N 0 } = λ p J vanishes in the state ρ 0 : J ρ0 = 0. Hence, since N 0 is the only quantity that brings the system out of equilibrium, the evolution of the whole state ρ must itself occur on time scales of order λ 2p . This yields in particular a relation analogous to eq. (25): where ad †H is the adjoint with respect to ρ 0 . Again, we compute that this operator acts on a function v as: Thus, combining eq. (29) and eq. (30), we derive that f must satisfy adH f = −βµJ in lowest order in λ, i.e. ad H0 f = −βµJ . Again, this equation needs to be regularized, and we get where, again, J (t) is the evolution of J under the free dynamics generated by H 0 . We come to the conclusion that

EXPLICIT EVALUATION OF THE DISSIPATION RATE IN SPECIFIC CASES
We compute explicitly the dissipation rate γ in the leading order in λ for the three cases where we want to compare our predictions with numerical data. Our starting point is always the expression (32) (even for p = 1, since if we take W = (β(ω(k) − µ) −1 , the expressions (28) and (32) coincide). Eq. (32) still contains some hidden dependence in λ through J and ρ 0 . To obtain the leading order, we replace · ρ0 by the average · over a Gaussian measure with density e −β(H0−µN0) /Z 0 . Omitting O(λ 2p+1 ) terms in our formulas for simplicity, and writing ε = τ −1 , we get r = 6 and δ > 0.3: In this case p = 1. Our aim is to show that where γ 0 is a number that depends only on the value of δ and that can be evaluated explicitly.

NUMERICAL PROCEDURE
All data points are generated by directly simulating the dynamics for the Hamiltonian H in eq. (1) with L = 1024 and periodic boundary conditions. The numerical scheme is a standard Strömer-Verlet algorithm with a time step ∆t = 0.1. For large λ, where one does not need to follow the dynamics on very long time scales, we have checked that changing L and ∆t only produces marginal differences.
A similar kind of initial state (with different choices for W ) is used in [32]. The data are averaged over 250 − 4000 initial configurations, corresponding to different realizations of ϕ k . Starting from the initial state (43), we expect that the pre-thermal state is reached on very short times, and this seems to be indeed the case, see the left panel on Fig. 5. Next, to measure the rate γ, we measure how N 0 /L evolves with time, and we observe that the evolution is first approximately linear, see the middle panel on Fig. 5. We identify the slope of this linear piece with γ(β, λ, δ). This should become exact in the limit λ → 0 that we investigate. For large λ, there is some arbitrariness in determining the time interval where the evolution is approximately linear. However, for smaller values of λ, this interval simply corresponds to the longest time on which one is reasonably able to run the simulations and perform sufficient averaging (t = 10 8 ). See Fig. 6. Finally, the value of N 0 /L reaches its equilibrium value on longer time scales, see the right panel on Fig. 5. Figure 5. N0/L as a function of time for r = 6, δ = 0.35 and λ = 10 −3 on three different time scales. Average over more than 2000 initial configurations. Left panel: 0 ≤ t ≤ 4 × 10 3 . We observe a leap between t = 0 and t = 40, corresponding presumably to the stage where the system moves to the pre-thermal state. Middle panel: 0 ≤ t ≤ 4 × 10 4 . The value of N0/L increases approximately linearly. The slope is taken as the value for γ to obtain the corresponding point on Fig. 2 in the main text. Right panel: 0 ≤ t ≤ 10 6 . N0/L reaches eventually its thermal value (N0/L) th 0.54 (computed at λ = 0) (orange). Figure 6. N0/L as a function of time for r = 4, δ = 0.45 and λ = 10 −3 (left panel) and λ = 5 × 10 −4 (right panel). The rate γ is determined by a mean square fit (orange). Average over more than 250 initial configurations.