Universal law of thermalization for one-dimensional perturbed Toda lattices

The Toda lattice is a nonlinear but integrable system. Here we study the thermalization problem in one-dimensional, perturbed Toda lattices in the thermodynamic limit. We show that the thermalization time, $T_{eq}$, follows a universal law; i.e., $T_{eq}\sim \epsilon^{-2}$, where the perturbation strength, $\epsilon$, characterizes the nonlinear perturbations added to the Toda potential. This universal law applies generally to weak nonlinear lattices due to their equivalence to perturbed Toda systems.

Many efforts have been made to explain FPU's results  and many theories have been put forward for understanding the relaxation problem, e.g., the soliton theory [9], the Chirikov resonance overlap theory [12,13], the mode-coupling theory [19], the q-breathers theory [20,21], as well as the Kolmogorov-Arnold-Moser theorem [22,23] enunciated by Kolmogorov at nearly the same time of FPU's work. However, the original goal of the FPU problem, that is, to answer whether a simple dynamical system can reach the thermalized state at arbitrarily weak nonlinearity and what properties the equipartition process may have, have not been achieved. To this end, some studies are inconsistent with each other [23]. For instance, for the FPU-β model, while Berchialla et al. showed that the equipartition time T eq depends on the energy density ε in a stretched exponential law in the thermodynamic limit, i.e., T eq ∼ exp(−ε 1/4 ) [32], DeLuca et al. suggested a power-law relationship T eq ∼ ε −3 instead [33]. Benettin et al. indicated however a crossover from the stretched exponential law, T eq ∼ exp(−ε 1/4 ) for the FPU-β model and T eq ∼ exp(−ε 1/8 ) for the FPU-αβ model, respectively, to the power-law T eq ∼ ε −9/4 for both cases [34,35].
Recently, the power-law relationship is affirmed [36][37][38][39] based on the wave turbulence (WT) theory [40][41][42][43][44][45]. It was shown analytically that the exact nontrivial six-wave resonant interactions are responsible for thermalization * yzhang75@xmu.edu.cn of short FPU chains in the weak nonlinear regime and, consequently, lead to T eq ∼ γ −8 for the FPU-α [36] and T eq ∼ γ −4 for the FPU-β [37] model, where γ is the nonlinearity strength defined, respectively, as γ = αε 1/2 and γ = βε for the two models (α and β are the coefficient, respectively, of the cubic and quartic nonlinear term). These results imply that any weak nonlinearity can ensure the system to be thermalized eventually. It was further conjectured (but not verified) that the nontrivial four-wave resonant interactions would dominate the thermalization process in the thermodynamic limit, leading to T eq ∼ γ −4 and T eq ∼ γ −2 for the FPU-α [36] and FPU-β [37] model, respectively.
These conjectures were partially verified in a very recent effort [46] where it was found that in the thermodynamic limit, a universal law, T eq ∼ γ −2 , applies generally to a class of one-dimensional (1D) lattices with interaction potential V (x) = x 2 /2 + λx n /n, where n ≥ 4 is an integer and γ = λε (n−2)/2 is the nonlinearity strength. It also applies to another class of 1D lattices with symmetric interaction potential V (x) = x 2 /2 + λ|x| d /d, where d = m 1 /m 2 > 2 with m 1 and m 2 being two coprime integers and the nonlinearity strength γ = λε (d−2)/2 . The existence of this universal law strongly confirms the assumption that the exact nontrivial wave-wave resonances dominate thermalization. However, it was also found that for a lattice with asymmetric potential interaction, though T eq still depends on γ in a power law, the exponent deviates from −2. In addition, the numerical result T eq ∼ γ −4.6 for the asymmetric FPU-α model [46] deviates from the conjectured T eq ∼ γ −4 [36] seriously as well.
Note that in all these works, the studied nonlinear model was considered to be a perturbed harmonic lattice. The harmonic lattice is integrable and linear. However, for a given nonlinear model, it can also be viewed alternatively as a perturbed Toda lattice [47] that is integrable but nonlinear. Interestingly, for some models, taking the latter viewpoint has been shown to be more consistent. The FPU-α model is a good example, for which supporting evidence from various aspects, e.g., by a normal mode approach [48], by the Lyapunov exponent analysis [49][50][51], and by thermalization process comparison [34,35], has been found. In particular, as inspired by Refs. [34,35], we have revisited the previous studies of thermalization and found that they strongly suggest the FPU-α model (the case of n = 3 in Ref. [46]) be viewed as the perturbed Toda lattice while other models as perturbed harmonic lattices. Given this, we were led to the conclusion that T eq ∼ γ −2 generally applies to the perturbed harmonic lattices in the thermodynamic limit but not to the systems out of this class [46].
In the present work, we study systematically the thermalization rate of 1D perturbed Toda lattices in order to find whether there exists a universal law of T eq for this class as well and if the answer is yes, how it differs from T eq ∼ γ −2 applicable to the 1D perturbed harmonic lattices. In the following, we will first introduce the models in the next section, then provide our theoretical arguments in Sec. III. The numerical approach, as well as simulation results, will be described and presented in Sec. IV, followed by the summary and discussions in Sec. V.

II. THE MODELS
We study the perturbed Toda models with potential where is the Toda potential [47] with α being a free parameter and θ n is the coefficient of the perturbation. The Toda potential can be expanded as Taylor's series: with the coefficients From Eq. (3), we can see that the harmonic model is a special case of the Toda model when α = 0. For this reason, any nonlinear model can be regarded as a perturbed Toda model as well.
To make our analysis more generalisable, we also study the generalized FPU model [34] with potential where α and θ n are free parameters. In principle, any smooth nonlinear potential can be written in this form. Comparing with the perturbed Toda model given by Eq. (1) whose potential is perturbed only in a single high-order term, the potential of this model can be perturbed in multiple high-order terms.

III. DEFINITION OF PERTURBATION STRENGTH AND THEORETICAL ANALYSIS
The Hamiltonian of our systems can be written as where H 0 and H ′ denote, respectively, the integrable part and the perturbation. Intuitively, the larger the perturbation is, the easier the system will be thermalized. A conventional practice is to take the Hamiltonian of the harmonic lattice as H 0 and defines the rest nonlinear part as the perturbation. However, the nonlinearity may not always be a good indicator to characterize the equipartition time. For example, the Toda model can own a very strong nonlinearity but will never be thermalized due to its integrability. Therefore, it is more reasonable to define perturbation strength as the nonintegrability, instead. To this end, it would be superior to adopt the Hamiltonian of the Toda model as H 0 . This scenario is general; it also covers the conventional one where H 0 is the Hamiltonian of the harmonic lattice when α = 0. The definition of perturbation strength for different cases is given below: The perturbed Toda model with n ≥ 4.-We get the perturbation by comparing Eq. (1) and Eq. (3) as By normalizing the Hamiltonian, i.e., by rescaling the relative displacement with the energy density so that x ′ = ε 1/2 x, we can obtain the dimensionless perturbation strength as The perturbed Toda model with n = 3.-For this case, Eq. (1) can be rewritten in Taylor's series as To setα = α + θ 3 and with the help of Eq. (4), we can get another Toda potential much closer to the perturbed system than the original one: Comparing Eq. (9) with Eq. (10), we get perturbation where denotes the nth-order perturbation strength. Note that the leading perturbation is still of the 4th-order for n = 3.
Comparing with ǫ n ∼ |θ n | for the case of n ≥ 4 at a fixed ε, here ǫ n has a more complicated relationship with θ 3 . The generalized FPU model. -Similarly, comparing Eq. (5) with Eq. (3), the perturbation can be identified to be where ǫ n is the dimensionless strength of the nth-order perturbation, given by Again, for this case the 4th-order perturbation is the lowest order one. For the FPU-α model, the leading perturbation strength is ǫ 4 = 2 3 α 2 , which is very different from the nonlinear strength α with respect to the linear integrable (harmonic) model. Now let us evaluate the equipartition time. Based on to the WT theory, it has been proved that either there are no resonances, or all of the scattering matrices are zero at all orders on the resonant manifold for integrable systems [52] so that they are characterized by trivial scattering processes and are never thermalized [36]. All the scattering matrix being zero is broken when the integrable system is perturbed. We assume that the exact nontrivial n-wave scattering processes caused by the nth-order perturbation, in the thermodynamic limit, dominate the thermalization process of the perturbed system. Then based on the theoretical results derived from the WT theory [36][37][38]46], the time scale of equipartition is where ǫ nL and n L denote, respectively, the strength and the order of the leading perturbation. Hereafter we will show that the above assumption is supported by extensive numerical simulations.

IV. NUMERICAL METHOD AND RESULTS
For a homogeneous lattice we consider here that consists of N + 1 particles of unit mass, labelled 0, 1, 2, · · · , N from the left to the right, its Hamiltonian is where p j and q j are, respectively, the momentum and the displacement from the equilibrium position of the jth particle, and V is the nearest-neighboring interaction potential.
For the fixed boundary conditions, i.e., q 0 = p 0 = q N = p N = 0, the normal modes are defined as To each mode k one can associate a harmonic energy and a phase ϕ k defined via Following the definition of equipartition, one expects where ε = E/(N − 1) is the energy density (E denotes the total energy of the system) andĒ k (T ) represents the time average of E k up to time T ; i.e., Here µ ∈ [0, 1) controls the size of time average window. In our numerical simulations, µ = 2/3 is fixed, which not only can speed up the calculations, but also has the advantage of a quicker loss of the memory of the very special initial state as proposed in Ref. [34]. Based on the definedĒ k (T ), we need introduce a parameter to measure how close the system is to equipartition. A frequently used parameter is the effective relative number of degrees of freedom [53,54]. Here we employ the quantity ξ(t) as in Ref. [34], i.e., where is the spectral entropy and When equipartition is approached, ξ will saturate at 1.
To integrate the motion equations numerically, we take the eighth-order Yoshida method [55]. The typical time step is ∆t = 0.1; the corresponding relative error in energy conservation, when all modes are excited and do contribute to the total energy, is around 10 −5 . A further decrease of the time step by one order of magnitude, i.e., ∆t = 0.01, does not change the results. To suppress fluctuations, the average is done over 24 phases uniformly distributed in [0, 2π], and we use · to denote the ensemble average results. Initially the lowest 10% of frequency modes are excited, α = −1 and N = 2048 are kept fixed throughout for all the numerical results presented. We have checked and verified that no qualitative difference will be resulted in neither when the percentage of the excited modes is changed nor when the system size is increased further.
In Fig. 1(a), the results of E k (t)/ε versus k/N for the Toda model are presented. It can be seen that only a small portion of the energy spread quickly from the initial excited low-frequency modes to the high-frequency modes, then the energy profile keeps a stable localized form with an exponential decaying tail. It suggests that for the Toda model, the thermalized state can never be reached. In Fig. 1(b), the results for the perturbed Toda model with θ 3 = −0.05, ǫ 4 = 1.025 × 10 −4 (solid lines) and θ 3 = 0.05, ǫ 4 = 0.975 × 10 −4 (dotted lines) are plotted. It can be seen that thermalization is faster approached in the former case as a consequence of the tiny difference of ǫ 4 , suggesting that the thermalization rate depends on the perturbation strength sensitively. It can also be seen that the energy of low-frequency modes remains for a long time and the energy of high-frequency modes increases very slowly, known as a signature of the metastable state [34,35,[56][57][58]. Nevertheless, the system will be thermalized eventually. As a comparison, Fig. 1(c) shows the results of the perturbed Toda model with θ 4 = −3, ǫ 4 = 3 × 10 −3 (solid lines) and θ 4 = 3, ǫ 4 = 3 × 10 −3 (dotted lines). Note that the thermalization rates of the two cases keep the same as the perturbation strengthes are identical. In addition, the system is fully thermalized at time T ∼ 10 7 , when E k /ε = 1.
To obtain the equipartition time, we study the properties of ξ(t) defined by Eq. (22). Figure 2(a) shows the results for the perturbed Toda model with θ 3 = 0.5. By varying energy density, ε, ǫ n is changed. Note that on a sufficiently large time scale, all values of ξ(t) increase from 0 to 1 with very similar sigmoidal profiles. It suggests that energy equipartition is finally achieved. Meanwhile, when the energy density decreases, the time required to reach the thermalized state increases. Now we adopt the definition of the equipartition time, T eq , as that when ξ(t) reaches the threshold value 0.5 as in Refs. [34,59]. Though assuming the threshold value 0.5 is artificial, it does not influence the scaling law of T eq [23]. This can be seen from Fig. 2(b), where the sigmoidal profiles in Fig. 2(a) can overlap with each other upon suitable shifts, which suggests that the concrete threshold value does not affect the scaling exponent of T eq . With these preparations, we are ready to present the results of T eq as a function of ǫ n . In Fig. 3(a), the numerical results of T eq as a function of θ 3 , θ 4 , θ 5 , θ 6 , and θ 7 are shown in semi-log scale for the perturbed Toda model. As for a given energy density, T eq for different n could be remarkably distinct, here we adopt different energy density for different n in order to present all the data in a single picture. We can see that all the numerical points can be well fitted with a Λ-shape curve of form T eq ∼ |θ n | −2 for n ≥ 4 and T eq ∼ |(θ 3 − 1) 2 − 1| −2 for n = 3, respectively, which is exactly what predicted by Eq. (15) when Eq. (12) and Eq. (8) for the two cases, respectively, are substituted into. The data of T eq are also plotted versus ǫ n in Fig. 3(b) in loglog scale. Note that all the points fall on the lines with a slope of −2, suggesting T eq ∼ ǫ −2 n holds for all the cases. As shown by Eq. (12), the relationship between ǫ n and θ 3 is very complicated. Next, we will carefully check if this relationship is true in a wider parameter range.  Figure 4(a) shows T eq versus θ 3 in a large parameter range, where three peaks of T eq can be clearly recognized. In order to understand the underlying mechanism for the formation of these peaks, we plot ǫ n as a function of θ 3 [see Eq. (12)] in semi-log scale in Fig. 4(b). The first peak resides at θ 3 = 0, corresponding to the integrable point of the Toda lattice, i.e., ǫ n = 0 for all n. The theoretical prediction gives T eq ∼ ǫ −2 4 ∼ |(θ 3 − 1) 2 − 1| −2 near this point, which agrees very well with the numerical results [see Fig. 4(a)]. It is amazing to realize that this theoretical result also predicts the location and the height of the third peak near θ 3 = 2 = −2α, which corresponds to the special point of ǫ n = 0 for all even n [see Fig. 4(b)]. Near this point the perturbed Toda model is very close to its 'mirror image' that adopts a minus α in the Toda potential [see Eq. (2)]. In particular, for θ 3 = 2, the perturbed system can be regarded to have a mirrored Toda potential with additional the fifth and higher odd order perturbations [see blue dashed line and green dashed dots line in Fig. 4(b)]. Note that the third peak is not located at θ 3 = 2 exactly due to the influence of high order resonances becoming non-negligible near this point. The middle peak is near θ 3 = 1 = −α, where the cubic coefficient is approximately zero. Namely, the perturbed system is close to the linear integrable point (harmonic lattice) and hence has an approximately symmetric potential, such that the additional energy mixing channel introduced by the asymmetry of interaction potential [46] is closed and thus gives rise to the middle peak.
So far our investigation has suggested that thermalization of weakly perturbed Toda lattices follows a common feature, i.e., T eq ∝ ǫ −2 . This general behavior coincides completely with that for the perturbed harmonic lattices [46]. The validity of our theoretical analysis can be further tested in the generalized FPU models, and meanwhile we find that it also gives satisfactory explanations to the numerical results reported previously [34]. Figure 5 summarizes the numerical results of the generalized FPU model. Figures 5(a)-(d) show T eq as a function of θ 4 , θ 5 , θ 6 , and θ 7 in semi-log scale, at two different energy densities, respectively. In Fig. 5(a) and (e), the results for repeating the previous study in Ref. [34] are presented, and extended study results are presented in Figs. 5(b)-(d) and (f)-(h). From (a) to (d), each figure shows a very marked peak near θ T n , and all the numerical points are well fitted with T eq ∼ |θ n − θ 0 n | −2 . But the center of the peak is not exactly at θ T n , which is mainly because the effect of higher order resonances becomes increasingly significant as θ n tends to θ T n . However, with the increase of n the difference between θ 0 n and θ T n decreases due to the fact that the higher the order, the weaker its effect [see the values indicated in Figs. 5(a)-(d)]. In Figs. 5(e)-(h), T eq is redrawn as a function of ǫ n that is defined by Eq. (14) with θ T n being replaced by θ 0 n , considering the higher order correction. Note that all the data points fall onto the lines with slope −2, suggesting that again, T eq ∼ ǫ −2 n is confirmed convincingly. Notice that expression (15) can converge to the theoretical results of Ref. [46], i.e., T eq ∝ λ −2 ε −(n−2) for n ≥ 4, which can be regarded as a perturbed linear integrable systems, i.e., a special case of α = 0 in our study here. What is interesting is that the FPU-α model with α = 0 and θ n = 0 is covered by expression (15) automatically. In this case, from Eq. (15) we have , which is the same as that given in Refs. [36,46].

V. SUMMARY AND DISCUSSIONS
In this work, we have shown that thermalization of a 1D weak nonlinear lattice exhibits a universal feature in the thermodynamic limit, i.e., T eq ∝ ǫ −2 , where ǫ is the perturbation strength defined as the difference in the potential of the system from the Toda potential. This universal behavior supports the assumption within the WT framework that the exact nontrivial wave-wave resonances dominate the thermalization process of a weak FIG. 5. The thermalization time Teq as function of θ4 with θ5, · · · , θ∞ = 0 (a); of θ5 with θ4 = θ T 4 and θ6, · · · , θ∞ = 0 (b); of θ6 with θn = θ T n (n = 4, 5) and θ7, · · · , θ∞ = 0 (c); and of θ7 with θn = θ T n (n = 4, 5, 6) and θ8, · · · , θ∞ = 0 (d), in semi-log scale, respectively. The vertical dashed (solid) lines are for θ = θ T n (θ = θ 0 n ), and the Λ-shape solid lines are for Teq ∼ |θn − θ 0 n | −2 , which are plotted for reference. Panels (e)-(h): the same as in (a)-(d) but shown as a function of ǫn in log-log scale instead. The solid lines with slope −2 are drawn for reference. The letters L and R in the legend indicate the points to the left and right of the peak in (a)-(d), respectively. nonlinear lattice. The key to identify the universal exponent −2 is to select the Toda lattice as the reference integrable system. In doing so, the third order nonlinearity has been found to be so crucial that it governs how we should assign the reference integrable system consistently. In particular, the system with (without) the cubic term of interactions should be regarded to be the perturbed Toda (harmonic) model. Comparing with previous studies, the resultant thermalization law (T eq ∝ ǫ −2 ) provides a unified and consistent picture for thermaliza-tion of one-dimensional nonlinear chains. ACKNOWLEDGMENT We are indebted to Prof. Jiao Wang for his kind help in preparing the manuscript. This work is supported by NSFC (Grant No. 11335006).