Universal law of thermalization for one-dimensional perturbed Toda lattices

TL;DR

This paper demonstrates that in one-dimensional perturbed Toda lattices, the thermalization time universally scales as the inverse square of the perturbation strength, revealing a fundamental law governing thermalization in weakly nonlinear systems.

Contribution

It establishes a universal law for thermalization time in perturbed Toda lattices, linking it to perturbation strength and extending to weak nonlinear lattices.

Findings

Thermalization time scales as $oxed{T_{eq} \\sim \\epsilon^{-2}}$.

Universal applicability to weak nonlinear lattices.

Thermalization behavior derived in the thermodynamic limit.

Abstract

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.

Equations

$$V(x)=V_T(\alpha ,x)+{\theta }_n{x}^n/n,\text{and}n\ge 3,$$

$$V_T(\alpha ,x)=\frac{e^{2\alpha x}-2\alpha x-1}{4{\alpha }^2}$$

$$V_T(\alpha ,x)=\frac{x^2}{2}+\frac{\alpha x^3}{3}+\sum _{n=4}^{\infty }\frac{{\theta }_n^T{x}^n}{n},$$

$${\theta }_n^T=\frac{(2\alpha {)}^{n-2}}{(n-1)!}.$$

$$V_{\text{gFPU}}(x)=\frac{x^2}{2}+\frac{\alpha x^3}{3}+\sum _{n=4}^{\infty }\frac{{\theta }_n{x}^n}{n},$$

$$H=H_0+H^{\prime },$$

$$H^{\prime }=\frac{{\theta }_n{x}^n}{n}.$$

$${ϵ}_n=|{\theta }_n|{\epsilon }^{(n-2)/2}.$$

$$V(x)=\frac{x^2}{2}+\frac{(\alpha +{\theta }_3)x^3}{3}+\frac{{\alpha }^2{x}^4}{6}+\frac{{\alpha }^3{x}^5}{15}+\cdots .$$

$$V_T(\tilde{\alpha },x)=\frac{x^2}{2}+\frac{\tilde{\alpha }{x}^3}{3}+\frac{{\tilde{\alpha }}^2{x}^4}{6}+\frac{{\tilde{\alpha }}^3{x}^5}{15}+\cdots .$$

$$H^{\prime }=\frac{{ϵ}_4{x}^4}{6}+\frac{{ϵ}_5{x}^5}{15}+\cdots ,$$

$${ϵ}_n=|(\alpha +{\theta }_3{)}^{n-2}-{\alpha }^{n-2}|{\epsilon }^{(n-2)/2},n=4,5,\cdots ,$$

$$H^{\prime }=\sum _{n=4}^{\infty }{ϵ}_n\frac{x^n}{n},$$

$${ϵ}_n=|{\theta }_n-{\theta }_n^T|{\epsilon }^{(n-2)/2},n=4,5,\cdots ,\infty .$$

$$T_{eq}\propto {ϵ}_{n_L}^{-2},$$

$$H=\sum _{j=1}^N\left[\frac{p_j^2}{2}+V(q_j-q_{j-1})\right],$$

$$\{\begin{array}{ll}{Q}_k& =\sqrt{\frac{2}{N}}{\sum }_{n=1}^N{q}_n\sin \left(\frac{nk\pi }{N}\right),\\ P_k& =\sqrt{\frac{2}{N}}{\sum }_{n=1}^N{p}_n\sin \left(\frac{nk\pi }{N}\right).\end{array}$$

$$E_k=\frac{1}{2}\left(P_k^2+{\omega }_k^2{Q}_k^2\right)$$

$$Q_k=\sqrt{2E_k/{\omega }_k^2}\sin \left({\phi }_k\right),P_k=\sqrt{2E_k}\cos \left({\phi }_k\right).$$

$$\underset{T\to \infty }{\lim }{\bar{E}}_k(T)\simeq \epsilon ,k=1,\cdots ,N,$$

$${\bar{E}}_k(T)=\frac{1}{(1-\mu )T}{\int }_{\mu T}^T{E}_k(P(t),Q(t))dt.$$

$$\xi (t)=\tilde{\xi }(t)\frac{e^{\eta (t)}}{N/2},$$

$$\eta (t)=-\sum _{k=N/2}^N{w}_k(t)\log [w_k(t)]$$

$$\tilde{\xi }(t)=\frac{{\sum }_{k=N/2}^N{\bar{E}}_k(t)}{\frac{1}{2}{\sum }_{1\le k\le N}{\bar{E}}_k(t)},w_k(t)=\frac{{\bar{E}}_k(t)}{{\sum }_{j=N/2}^N{\bar{E}}_j(t)}.$$

Full text

I Introduction

The problem of thermalization in weak nonlinear systems has a long studying history but is still far from being resolved. The first numerical experiment, aiming to observe the rates of mixing and thermalization in a system of reversible microscopic dynamics, was accomplished in the 1950’s by Fermi, Pasta, Ulam (FPU) Fermi et al. (1955), and Tsingou Fermi et al. (1955); Dauxois (2008). Their numerical results showed very little tendency toward equipartition of energy among the degrees of freedom, which is known as the FPU recurrence Fermi et al. (1955); Dauxois (2008); Wu and Patton (2007); Mussot et al. (2014); Bao et al. (2016); Guasoni et al. (2017); Mussot et al. (2018); Pierangeli et al. (2018); Zabusky and Kruskal (1965). This seminal work failed to observe the expected picture but opened up two entirely new research fields: nonlinear science Zabusky and Kruskal (1965); Dauxois and Peyrard (2006); Gardner et al. (1967); Izrailev and Chirikov (1966); Chirikov (1979); Campbell et al. (2005); Berman and Izrailev (2005); Pettini et al. (2005) and computational science Zabusky (2005); Porter et al. (2009).

Many efforts have been made to explain FPU’s results Fermi et al. (1955); Dauxois (2008); Wu and Patton (2007); Mussot et al. (2014); Bao et al. (2016); Guasoni et al. (2017); Mussot et al. (2018); Pierangeli et al. (2018); Zabusky and Kruskal (1965); Dauxois and Peyrard (2006); Gardner et al. (1967); Izrailev and Chirikov (1966); Chirikov (1979); Campbell et al. (2005); Berman and Izrailev (2005); Pettini et al. (2005); Zabusky (2005); Porter et al. (2009); Lepri (1998); Flach et al. (2005, 2006); Rangarajan (1998); Gallavotti (2008); Ford and Lunsford (1970); Tuck and Menzel (1972); Bivins et al. (1973); Sholl (1990); Ford (1992); Zaslavsky (2005); Carati et al. (2005); Danieli et al. (2017) and many theories have been put forward for understanding the relaxation problem, e.g., the soliton theory Zabusky and Kruskal (1965), the Chirikov resonance overlap theory Izrailev and Chirikov (1966); Chirikov (1979), the mode-coupling theory Lepri (1998), the $q$-breathers theory Flach et al. (2005, 2006), as well as the Kolmogorov-Arnold-Moser theorem Rangarajan (1998); Gallavotti (2008) 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 Gallavotti (2008). For instance, for the FPU-$\beta$ model, while Berchialla et al. showed that the equipartition time $T_{eq}$ depends on the energy density $\varepsilon$ in a stretched exponential law in the thermodynamic limit, i.e., $T_{eq}\sim\exp(-\varepsilon^{1/4})$ Berchialla et al. (2004), DeLuca et al. suggested a power-law relationship $T_{eq}\sim\varepsilon^{-3}$ instead DeLuca et al. (1995); *PhysRevE.60.3781. Benettin et al. indicated however a crossover from the stretched exponential law, $T_{eq}\sim\exp(-\varepsilon^{1/4})$ for the FPU-$\beta$ model and $T_{eq}\sim\exp(-\varepsilon^{1/8})$ for the FPU-$\alpha\beta$ model, respectively, to the power-law $T_{eq}\sim\varepsilon^{-9/4}$ for both cases Benettin and Ponno (2011); Benettin et al. (2013).

Recently, the power-law relationship is affirmed Onorato et al. (2015); Lvov and Onorato (2018); Pistone et al. (2018); Bustamante et al. (2018) based on the wave turbulence (WT) theory Zakharov et al. (1992); Majda et al. (1997); Zakharov et al. (2001, 2004); Nazarenko (2011); Sagaut and Cambon (2018). It was shown analytically that the exact nontrivial six-wave resonant interactions are responsible for thermalization of short FPU chains in the weak nonlinear regime and, consequently, lead to $T_{eq}\sim\gamma^{-8}$ for the FPU-$\alpha$ Onorato et al. (2015) and $T_{eq}\sim\gamma^{-4}$ for the FPU-$\beta$ Lvov and Onorato (2018) model, where $\gamma$ is the nonlinearity strength defined, respectively, as $\gamma=\alpha\varepsilon^{1/2}$ and $\gamma=\beta\varepsilon$ for the two models ($\alpha$ and $\beta$ 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}\sim\gamma^{-4}$ and $T_{eq}\sim\gamma^{-2}$ for the FPU-$\alpha$ Onorato et al. (2015) and FPU-$\beta$ Lvov and Onorato (2018) model, respectively.

These conjectures were partially verified in a very recent effort Fu et al. (2018) where it was found that in the thermodynamic limit, a universal law, $T_{eq}\sim\gamma^{-2}$, applies generally to a class of one-dimensional (1D) lattices with interaction potential $V(x)=x^{2}/2+\lambda x^{n}/n$, where $n\geq 4$ is an integer and $\gamma=\lambda\varepsilon^{(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+\lambda|x|^{d}/d$, where $d=m_{1}/m_{2}>2$ with $m_{1}$ and $m_{2}$ being two coprime integers and the nonlinearity strength $\gamma=\lambda\varepsilon^{(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 $\gamma$ in a power law, the exponent deviates from $-2$. In addition, the numerical result $T_{eq}\sim\gamma^{-4.6}$ for the asymmetric FPU-$\alpha$ model Fu et al. (2018) deviates from the conjectured $T_{eq}\sim\gamma^{-4}$ Onorato et al. (2015) 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 Toda (1967) that is integrable but nonlinear. Interestingly, for some models, taking the latter viewpoint has been shown to be more consistent. The FPU-$\alpha$ model is a good example, for which supporting evidence from various aspects, e.g., by a normal mode approach Ferguson et al. (1982), by the Lyapunov exponent analysis Casetti et al. (1997); Cerruti-Sola et al. (2000); Benettin et al. (2018), and by thermalization process comparison Benettin and Ponno (2011); Benettin et al. (2013), has been found. In particular, as inspired by Refs. Benettin and Ponno (2011); Benettin et al. (2013), we have revisited the previous studies of thermalization and found that they strongly suggest the FPU-$\alpha$ model (the case of $n=3$ in Ref. Fu et al. (2018)) 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}\sim\gamma^{-2}$ generally applies to the perturbed harmonic lattices in the thermodynamic limit but not to the systems out of this class Fu et al. (2018).

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}\sim\gamma^{-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

[TABLE]

where

[TABLE]

is the Toda potential Toda (1967) with $\alpha$ being a free parameter and $\theta_{n}$ is the coefficient of the perturbation. The Toda potential can be expanded as Taylor’s series:

[TABLE]

with the coefficients

[TABLE]

From Eq. (3), we can see that the harmonic model is a special case of the Toda model when $\alpha=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 Benettin and Ponno (2011) with potential

[TABLE]

where $\alpha$ and $\theta_{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

[TABLE]

where $H_{0}$ and $H^{\prime}$ 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 $\alpha=0$. The definition of perturbation strength for different cases is given below:

The perturbed Toda model with $n\geq 4$.—We get the perturbation by comparing Eq. (1) and Eq. (3) as

[TABLE]

By normalizing the Hamiltonian, i.e., by rescaling the relative displacement with the energy density so that $x^{\prime}=\varepsilon^{1/2}x$, we can obtain the dimensionless perturbation strength as

[TABLE]

The perturbed Toda model with $n=3$.—For this case, Eq. (1) can be rewritten in Taylor’s series as

[TABLE]

To set $\tilde{\alpha}=\alpha+\theta_{3}$ and with the help of Eq. (4), we can get another Toda potential much closer to the perturbed system than the original one:

[TABLE]

Comparing Eq. (9) with Eq. (10), we get perturbation

[TABLE]

where

[TABLE]

denotes the $n$th-order perturbation strength. Note that the leading perturbation is still of the 4th-order for $n=3$. Comparing with $\epsilon_{n}\sim|\theta_{n}|$ for the case of $n\geq 4$ at a fixed $\varepsilon$, here $\epsilon_{n}$ has a more complicated relationship with $\theta_{3}$.

The generalized FPU model. — Similarly, comparing Eq. (5) with Eq. (3), the perturbation can be identified to be

[TABLE]

where $\epsilon_{n}$ is the dimensionless strength of the $n$th-order perturbation, given by

[TABLE]

Again, for this case the 4th-order perturbation is the lowest order one. For the FPU-$\alpha$ model, the leading perturbation strength is $\epsilon_{4}=\frac{2}{3}\alpha^{2}$, which is very different from the nonlinear strength $\alpha$ 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 Zakharov and Schulman (1988) so that they are characterized by trivial scattering processes and are never thermalized Onorato et al. (2015). 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 $n$th-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 Onorato et al. (2015); Lvov and Onorato (2018); Pistone et al. (2018); Fu et al. (2018), the time scale of equipartition is

[TABLE]

where $\epsilon_{n_{L}}$ 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 [math], $1$, $2$, $\cdots$, $N$ from the left to the right, its Hamiltonian is

[TABLE]

where $p_{j}$ and $q_{j}$ are, respectively, the momentum and the displacement from the equilibrium position of the $j$th 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

[TABLE]

To each mode $k$ one can associate a harmonic energy

[TABLE]

and a phase $\varphi_{k}$ defined via

[TABLE]

Following the definition of equipartition, one expects

[TABLE]

where $\varepsilon=E/(N-1)$ is the energy density ($E$ denotes the total energy of the system) and $\bar{E}_{k}(T)$ represents the time average of $E_{k}$ up to time $T$; i.e.,

[TABLE]

Here $\mu\in[0,1)$ controls the size of time average window. In our numerical simulations, $\mu=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. Benettin and Ponno (2011).

Based on the defined $\bar{E}_{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 Livi et al. (1985); Goedde et al. (1992). Here we employ the quantity $\xi(t)$ as in Ref. Benettin and Ponno (2011), i.e.,

[TABLE]

where

[TABLE]

is the spectral entropy and

[TABLE]

When equipartition is approached, $\xi$ will saturate at $1$.

To integrate the motion equations numerically, we take the eighth-order Yoshida method Yoshida (1990). The typical time step is $\Delta 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., $\Delta t=0.01$, does not change the results. To suppress fluctuations, the average is done over $24$ phases uniformly distributed in $[0,2\pi]$, and we use $\langle\cdot\rangle$ to denote the ensemble average results. Initially the lowest $10\%$ of frequency modes are excited, $\alpha=-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 $\langle E_{k}(t)/\varepsilon\rangle$ 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 $\theta_{3}=-0.05$, $\epsilon_{4}=1.025\times 10^{-4}$ (solid lines) and $\theta_{3}=0.05$, $\epsilon_{4}=0.975\times 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 $\epsilon_{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 Benettin et al. (2009); Benettin and Ponno (2011); Ponno et al. (2011); Benettin et al. (2013); Goldfriend and Kurchan (2018). Nevertheless, the system will be thermalized eventually. As a comparison, Fig. 1(c) shows the results of the perturbed Toda model with $\theta_{4}=-3$, $\epsilon_{4}=3\times 10^{-3}$ (solid lines) and $\theta_{4}=3$, $\epsilon_{4}=3\times 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\sim 10^{7}$, when $\langle E_{k}/\varepsilon\rangle=1$.

To obtain the equipartition time, we study the properties of $\langle\xi(t)\rangle$ defined by Eq. (22). Figure 2(a) shows the results for the perturbed Toda model with $\theta_{3}=0.5$. By varying energy density, $\varepsilon$, $\epsilon_{n}$ is changed. Note that on a sufficiently large time scale, all values of $\langle\xi(t)\rangle$ increase from [math] 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 $\langle\xi(t)\rangle$ reaches the threshold value $0.5$ as in Refs. Benettin and Ponno (2011); Zhang et al. (2016). Though assuming the threshold value $0.5$ is artificial, it does not influence the scaling law of $T_{eq}$ Gallavotti (2008). 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 $\epsilon_{n}$.

In Fig. 3(a), the numerical results of $T_{eq}$ as a function of $\theta_{3}$, $\theta_{4}$, $\theta_{5}$, $\theta_{6}$, and $\theta_{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 $\Lambda$-shape curve of form $T_{eq}\sim|\theta_{n}|^{-2}$ for $n\geq 4$ and $T_{eq}\sim|(\theta_{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 $\epsilon_{n}$ in Fig. 3(b) in log-log scale. Note that all the points fall on the lines with a slope of $-2$, suggesting $T_{eq}\sim\epsilon_{n}^{-2}$ holds for all the cases. As shown by Eq. (12), the relationship between $\epsilon_{n}$ and $\theta_{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 $\theta_{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 $\epsilon_{n}$ as a function of $\theta_{3}$ [see Eq. (12)] in semi-log scale in Fig. 4(b). The first peak resides at $\theta_{3}=0$, corresponding to the integrable point of the Toda lattice, i.e., $\epsilon_{n}=0$ for all $n$. The theoretical prediction gives $T_{eq}\sim\epsilon_{4}^{-2}\sim|(\theta_{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 $\theta_{3}=2=-2\alpha$, which corresponds to the special point of $\epsilon_{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 $\alpha$ in the Toda potential [see Eq. (2)]. In particular, for $\theta_{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 $\theta_{3}=2$ exactly due to the influence of high order resonances becoming non-negligible near this point. The middle peak is near $\theta_{3}=1=-\alpha$, 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 Fu et al. (2018) 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}\propto\epsilon^{-2}$. This general behavior coincides completely with that for the perturbed harmonic lattices Fu et al. (2018). 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 Benettin and Ponno (2011).

Figure 5 summarizes the numerical results of the generalized FPU model. Figures 5(a)-(d) show $T_{eq}$ as a function of $\theta_{4}$, $\theta_{5}$, $\theta_{6}$, and $\theta_{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. Benettin and Ponno (2011) 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 $\theta_{n}^{T}$, and all the numerical points are well fitted with $T_{eq}\sim|\theta_{n}-\theta_{n}^{0}|^{-2}$. But the center of the peak is not exactly at $\theta_{n}^{T}$, which is mainly because the effect of higher order resonances becomes increasingly significant as $\theta_{n}$ tends to $\theta_{n}^{T}$. However, with the increase of $n$ the difference between $\theta_{n}^{0}$ and $\theta_{n}^{T}$ 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 $\epsilon_{n}$ that is defined by Eq. (14) with $\theta_{n}^{T}$ being replaced by $\theta_{n}^{0}$, considering the higher order correction. Note that all the data points fall onto the lines with slope $-2$, suggesting that again, $T_{eq}\sim\epsilon_{n}^{-2}$ is confirmed convincingly.

Notice that expression (15) can converge to the theoretical results of Ref. Fu et al. (2018), i.e., $T_{eq}\propto\lambda^{-2}\varepsilon^{-(n-2)}$ for $n\geq 4$, which can be regarded as a perturbed linear integrable systems, i.e., a special case of $\alpha=0$ in our study here. What is interesting is that the FPU-$\alpha$ model with $\alpha\neq 0$ and $\theta_{n}=0$ is covered by expression (15) automatically. In this case, from Eq. (15) we have $T_{eq}\propto\epsilon_{4}^{-2}=|0-\theta_{4}^{T}|^{-2}\varepsilon^{-2}=4/9\alpha^{-4}\varepsilon^{-2}$, which is the same as that given in Refs. Onorato et al. (2015); Fu et al. (2018).

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}\propto\epsilon^{-2}$, where $\epsilon$ 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 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}\propto\epsilon^{-2}$) provides a unified and consistent picture for thermalization 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).