Universal route to thermalization in weakly-nonlinear one-dimensional chains

We apply Wave Turbulence theory to describe the dynamics on nonlinear one-dimensional chains. We consider $\alpha$ and $\beta$ Fermi-Pasta-Ulam-Tsingou (FPUT) systems, and the discrete nonlinear Klein-Gordon chain. We demonstrate that resonances are responsible for the irreversible transfer of energy among the Fourier modes. We predict that all the systems thermalize for large times, and that the equipartition time scales as a power-law of the strength of the nonlinearity. Our methodology is not limited to only these systems and can be applied to the case of a finite number of modes, such as in the original FPUT experiment, or to the thermodynamic limit, i.e. when the number of modes approach infinity. In the latter limit, we perform state of the art numerical simulations and show that the results are consistent with theoretical predictions. We suggest that the route to thermalization, based only on the presence of exact resonance, has universal features. Moreover, a by-product of our analysis is the asymptotic integrability, up to four wave interactions, of the discrete nonlinear Klein-Gordon chain.


INTRODUCTION
The numerical experiment of Fermi, Pasta, Ulam and Tsingou (FPUT) in 1955 [15] has been of great influence to physics in many respects. The idea of the experiment was related to the role of chaos on the foundations of statistical mechanics. As chaotic systems are highly irregular, like stochastic ones, they loose memory of initial conditions rapidly, as asked by ergodic hypothesis. In this perspective, the result of Poincaré about the non existence in generic Hamiltonian systems of first integrals of motion, other than the energy, seems positive. In his early theoretical activity, Fermi generalised the result of Poincaré showing that in Hamiltonian systems with N > 2 degrees of freedom no smooth surface can divide the phase space into two regions containing invariant sets. From this correct result, Fermi deduced that non-integrable hamiltonian systems are generically ergodic and, even in absence of a rigorous proof, the ergodic problem was considered basically solved. However, Fermi came back to the problem after the war with this numerical experiment, probably too much absorbed by the development of the just born quantum mechanics in the meanwhile.
With one of the first available computers, they investigated the dynamics of a simple system of springs and masses, where the force is only between nearest neighbors. The novelty consisted in the metodology (it was one of the first numerical studies), but also in the fact that the force was not linear. In the case of linear forces, it is known that the Hamiltonian of the system can be diagonalized, the eigenstates are linear waves, and there is no interaction between them, namely the initial energy distribution among the modes does not change over time. Fermi and collaborators wanted to test whether introducing a small nonlinear term in the equations of motion (they used quadratic, cubic and split-linear terms) would allow energy to spread among the linear modes, and attain equipartition, confirming the earlier idea of Fermi that ergodic hypothesis can be considered true, giving a dynamical justification to statistical mechanics. However, this was not the case and instead quasiperiodic motion was observed or, at most, just a few modes close to the ones initially excited were active for the whole simulation time.
The surprising results of FPUT generated a lot of interest, and sparked a research that resulted in a large body of work, as witnessed by an excellent review [17]. In particular, even though ergodicity and chaos have been clarified to be largely irrelevant for statistical mechanics of macroscopic systems, that is when N 1 [24,27,29], the existence of solitons arise from a possible explanation of the FPUT phenomenology [40]. Moreover, the FPUT experiment was the first showing that numerical simulations are a powerful instrument to get physical insights of complex phenomena, and therefore an unavoidable tool for theoretical progress. In parallel, in the same time it was shown by mathematicians with the KAM theorem [19,22,25]  perturbations. This fundamental result showed that the first guess by Fermi was incorrect.
Despite these outstanding developments and considerable efforts, a clear and exhaustive explanation of the FPUT problem is still lacking. In particular, while there is some consensus concerning the metastability picture [3,16], the precise mechanisms driving to equipartition on very long time-scales is still elusive. One crucial point is that studies have focused on peculiar initial conditions and the specific phenomena that emerge from them in the short or medium timescales. For example, in the original experiments and in most of the successive studies the initial conditions consisted in exciting the low-frequency modes [2], which leads to the formation of coherent nonlinear structures, that is solitons, at least in the short-time. Some attention has been focused also on high-frequency initial conditions and this led to the study of breathers. However, it is difficult to understand all these peculiar direction in one larger vision of the FPUT problem. The problem with generic initial conditions seems to be yet more relevant from a general perspective, as recognised in an important recent contribution [7].
In the present work, we propose a more unifying picture; therefore, we seek universality traits that can explain the dynamics in a way that is only weakly dependent on the particular microscopic processes of the nonlinear chain in consideration. For this reason, we will consider generic initial conditions that are randomly out of equilibrium, involving low and high modes. Moreover, we will use tools and methods of the so-called Wave Turbulence (WT) theory. Originally introduced in fluid dynamics, WT is a statistical mechanics theory of weakly interacting waves, and it has been recently applied in relation to one dimensional anharmonic chains. In some sense, our approach is in line with the recent results that have clearly showed that the α-FPUT chain is a perturbation of Toda lattice, which is an integrable system [1]. Our main motto is that the irreversible transfer of energy in the spectrum in a weakly nonlinear system is achieved by exact resonant wave-wave interactions. Such resonances are the base for the WT theory and are responsible for the phenomenon of thermalization.
In this paper we give an overview on resonances and WT in a comprehensive manner as well as present new results regarding the limit of a large number of modes, providing insight in the fundamental mechanics of anharmonic chains. Our main result is establishment of the power-law scaling of the equipartition time T eq as a function of the nonlinearity strength (to be defined precisely later). In varying the number of elements in the chain, we demonstrate the crossover between two different scaling laws (lower exponent for large systems, that is the thermodynamic limit, steeper for smaller ones).

THE MODELS
We consider Hamiltonian systems of the form where H lin is the integrable Hamiltonian, corresponding to a linear dispersive dynamics, and H nonlin that contains the anharmonicity of the potential. More in particular we will deal with the α, β-FPUT and the discrete nonlinear Klein-Gordon (DNKG) (also called the φ 4 ) systems, all characterized by the linear Hamiltonian where N is the number of elements in the chain, q j is the j-th coordinate and p j its conjugate momentum and m models the linear part of the site-potential. For the α and β-FPUT systems m = 0. The nonlinear part of the Hamiltonian for the three models is given by: Without loss of generality, we use the same parameter β in front for nonlinearity in the DNKG and the β-FPUT. We will use periodic boundary conditions. It should be noted, given an initial energy, that there is a finite yet very small probability that the α-FPUT chain breaks up, since the nonlinear part of the Hamiltonian is not bounded from below. While being an inherently ill-posed problem, it has been shown that for low enough energies the α-FPUT chain can be stable for long times [6]. For the other models, we restrict ourselves to the case where the quartic terms are positive, that is β > 0, so that no blow-ups can occur. The parameter m is constrained to be greater than zero in the DNKG model, and we can anticipate that it will play a secondary role since we are interested in the effects of the nonlinear interactions. The presence of the quadratic finite difference terms in the linear Hamiltonian is important, because without that the eigenstates of the linear system would not be waves. Our approach is of perturbative nature, hence we need a way to quantify the nonlinearity strength. In general, the linear and nonlinear parts of the Hamiltionan are not conserved separately, but only the total Hamiltonian is. However, in typical situations where the nonlinear interactions are not too strong, operatively, we define the following nondimensional parameter for the α-FPUT and for the β-FPUT and DNKG nlin (j) is the nonlinear part of the Hamiltonian density of the j-th particle in the for α-FPUT. The absolute value in eq. (4) is necessary because h nlin (j) can be negative and cancellations in the sum would cause an incorrect accounting of the nonlinearity strength. The different scalings in equations (4) and (5) are due to the different degrees of nonlinearity between the α-FPUT and the other two models. Definitions (4) and (5) allows us to refine our statement that the nonlinearity must be small, that is, we require that the parameters are much less than one. In general, they fluctuate over time, so an appropriate averaging is needed to have a meaningful estimate of the nonlinearity strength. The parameters are further discussed when presenting the numerical results.

Fourier space
The eigenstates of the linear Hamiltonian are waves, hence it is useful to work in Fourier space. We define the direct and inverse discrete Fourier transform of the q j variables, and similar definitions forp k . We will use the convention that 0 ≤ k < N , and we noteq * k =q N −k andp * k =p N −k . After this change of variables, and using that is the Kronecker's delta modulo N , the linear part of the Hamiltonians can be written as where we have defined the linear dispersion relation (we recall that m = 0 for the FPUT models). The nonlinear part of the Hamiltonians couple the modes and can be interpreted as n-wave collision terms. For the sake of brevity, we denoteq j =q kj and δ (N ) 1+2+... = δ(k 1 + k 2 + ... mod N ). For the FPUT models we obtain where the collision matrices A 1,2,3 and B 1,2,3 are (after symmetrization) The complex exponential is relevant only when crossing Brillouin zones in the sum of wave numbers. For the DNKG model, it turns out that the interaction matrix is equal to unity, hence the nonlinear part of the Hamiltonian is simply From eq. (7), we can further simplify the problem using the normal modes representation, as the linear part of the Hamiltonians becomes simply whereas the nonlinear terms in the FPUT models and in the DNKG become (after renaming N − k i = k i where needed and dropping the prime) with the interaction matrices given by The nonlinear terms show that a nonlinearity in physical space turns into n-mode collision terms such as a 1 a 2 a 3 δ (N ) 1+2−3 . From the Hamilton equations, we can obtain the dynamical equationsȧ k = −i(1/N )δH/δa * k , which read for the α-FPUT and for the β-FPUT and DNKG These types of Hamiltonians are the canonical form of Hamiltonians for system composed of interacting particles or waves. If nonlinearity is small in the above mentioned sense, then there is a well developed theory called Wave Turbulence theory, which is described in details in [14,32] and briefly in the subsequent section.
The number of modes in the nonlinear terms defines the collision order, that is the α-FPUT contains 3-wave collisions, while β-FPUT DNKG 4-wave collisions (we will see that the effective collision order in the α-FPUT is also of the 4-wave type). The number of non-conjugate and conjugate variables (which matches the number of positive vs. negative terms in the Kronecker's deltas) defines different collision processes, e.g. a 1 a 2 a 3 δ 1+2−3−4 is a 2 → 2 collision process (scattering). The WT approach starts from the equations of motion in the canonical variables a k . Note that so far we have not taken any approximation; we only made the choice of using periodic boundary conditions.

WAVE-TURBULENCE THEORY IN A NUTSHELL
When the number of modes is large, the microscopic dynamics given by equations (18) and (19) do not provide much analytical insight and a statistical approach is needed. WT is the general statistical theory valid in the weak-nonlinear regime. While the theory has been recently developed for higher-order statistical observables [9,13,32] , the core of the WT is the kinetic equation developed for the prediction of the energy spectrum. The basic statistical observable is the two point correlator a 1 a * 2 which, under the hypothesis of homogeneity, it is given by with n k = n(k, t) the wave actions spectral density, i.e. the Fourier transform of the autocorrelation function of a(x, t). The average is over an ensemble of realizations of the system with the same initial linear energy H lin . Then, the energy spectrum that gives the mean energy per mode is The evolution equation for the wave action spectral density can be obtained by various techniques [9,13,14,32,34]. The main issue is due to the fact that calculating the evolution of n k from the equation of motion (18) or (19), one encounters the well know problem of the BBGKY hierarchy [11], that is, the evolution of lower order correlators depends on higher order correlators and a closure problem arises. Originally, a closure was obtained with a quasi-gaussian approximation which allowed to expresses higher cumulants in terms of lower ones [14]. It turns out that less stringent hypotheses are necessary. It is sufficient to make a random-phase approximation and to consider initial conditions which have also random amplitudes [9,13,32]. The random phase assumption is generally deemed to be solid because the linear waves decorrelate quickly even in the linear regime, hence it is expected that such property still holds true in the weakly nonlinear regime, as corroborated by numerical experiments [8] In all the following treatment, we assume that this property holds, and that the phases are independent random variables uniformly distributed between [0, 2π). The resulting kinetic equation can be obtained in the thermodynamic limit N → ∞. Generally also the continuum limit is taken at the same time, and both physical x and momentum space k are continuous.
The main concept underneath the kinetic equation is the existence of conservation laws associated to the wave scattering processes. Indeed, each nonlinear term in the equations (18) and (19) contains an appropriate Kroneker δ function over wave numbers. Depending on the shape of the dispersion relation, it may be possible to associate also a related conservation of energy. To be more specific, let us consider for example the scattering processes, contained in equation (18), then the main question is to verify if a similar relation holds for frequencies, i.e.: Equations (22) and (23) define a resonant manifold in wave number space, whose existence, as anticipated, depends only the analytical form of the dispersion curve ω k . If equation (22) and (23) are satisfied for the same wave numbers, then we are dealing with a resonant process which, according to the kinetic equation, will lead to an irreversible transfer of energy. There are two main types of wave systems that are typically considered in the framework of WT: the one dominated by three wave resonances and by four wave resonances (higher order are also possible). Examples of the systems dominated by three wave resonances include capillary waves on a surface of fluid and internal waves in the ocean. Examples of four wave resonant systems include gravity waves on deep water, as well as the celebrated 2D-Nonlinear Schrödinger equation.
For capillary waves on a surface three wave resonant interactions are possible and a kinetic equation can be derived, [14]: where k is a two dimensional wave vector, dk 23 = dk 2 dk 3 and {(k 1 ↔ k 2 )} implies the same integral with k 1 and k 2 exchanged. A 123 is a matrix that can be derived directly from the dynamical equations for capillary waves.
For surface gravity waves in deep water, it can be show that the leading order resonant process is the four-wave one, characterized by the following resonant process: The appropriate kinetic equation describing such process is Note that in both kinetic equation the interaction matrix is squared. The latter equation, with a transport term, a forcing term due to the wind and term mimicking the dissipation, is integrated daily for operationally ocean wave forecasting pourposes, see [21]. The kinetic equations (24) and (26) have a number of important properties. First of all, it can be shown that both equations conserves the total energy and momentum: The kinetic equation for four-wave interactions preserves also the number of waves (particles) in the scattering: The kinetic equation (26) is reminiscent of the Boltzmann equation for hard spheres [38] and there exists an entropy function such that dS/dt ≥ 0 (the same entropy can be defined for the three wave kinetic equation). The isotropic equilibrium is reached when dS/dt = 0, i.e. at the Rayleigh-Jeans distribution: where T is some normalization constant linked to the total energy and µ is a chemical potential linked to the conservation of the total number of particles (28) and it is zero in the case of a three wave process. We see that when µ = 0 then WT predicts a relaxation towards equipartition of energy. In this framework, exact resonances bring the system to thermalization with a time-scale which is proportional to the coupling coefficient in front of the collision integral. Some remarks are in order.
• Given a physical dispersive waves system, it has to be checked if exact resonances are allowed. In some cases, only trivial collisions or processes which are not able to transfer energy are found. In those cases, it is necessary to pursue the perturbative procedure and derive a kinetic equation with higher-order processes to check if at some point resonances are found. This issue is related to the dispersion relation, and in physical phenomena 3−, 4−, 5− and 6−wave processes have been recognised [32]. In other cases, it can be shown that no exact resonances exist and therefore no exchange of energy is possible, as rigorously shown to be the case for integrable systems [41]. Other systems may show some resonances which are not yet able to transfer energy among modes because resonating waves belong to isolated clusters [5].
• The WT framework is an asymptotic theory valid in the weakly nonlinear régime. If the nonlinearity is not small, different effects can play a role. In particular, the broadening of the frequencies due to the finite value of the nonlinear coupling, may trigger some resonances which were not present in the weak asymptotic limit, therefore allowing a more efficient transfer of energy.
In the present paper, we estimate the time scale of equipartition; therefore, it is of major importance to be able to derive a kinetic equation for our chain models. One should note that chains are intrinsically discrete in physical space and, for a finite number of masses, also the Fourier space is discrete. This issue makes the derivation of the kinetic equation not straightforward and, only in the thermodynamic limit, a formal derivation of the kinetic equation is possible. We shall deal with these points in the following section in some detail.

RESONANCE AND EFFECTIVE HAMILTONIANS IN THE THERMODYNAMIC LIMIT
We consider the thermodynamic limit which consists in taking the number of particles going to infinity, N → ∞, keeping constant the mass linear density. The general linear dispersion relation becomes with k real in the interval 0 ≤ k < 2π. We consider interaction processes from X to Y waves for which the following equations are satisfied: Our goal is to find resonances. Once the leading order resonance (lowest value of X and Y ) has been found, using tools from Hamiltonian mechanics, an effective Hamiltonian can be established and a kinetic equation can be built. The first trivial observation is that resonances of the type X to 0 are excluded for all three models; this is because ω k ≥ 0, and it is equal to zero only for m = 0 for mode k = 0 which does not play any role in the dynamics (for m = 0).

Three-wave resonant interactions in the α-FPUT
The leading order nonlinearity in the evolution equation for the α-FPUT is cubic; this implies three-wave interactions of the form 2 → 1, 1 → 2 or 3 → 0 (the latter have been already excluded). It is known, see for example [4,31], that the function f k = 2| sin(k/2)| is a subadditive function, i.e.
The equality holds only for k 1 or k 2 = l2π, with l ∈ Z; those wave numbers do not enter into the dynamics. The subadditivity implies the non existence of three-wave resonant interactions for m = 0.
Four-wave resonant interactions in the α, β-FPUT and DNKG models The β-FPUT and DNKG models include processes involving four waves, while, at first sight, the α-FPUT model does not include them; however, as it will be shown in the following, the absence of three-wave resonances allows for a perturbative approach on the α-FPUT model that leads to four and higher order wave-wave interaction processes. In [4] it has been shown that resonant processes 3 → 1 or 1 → 3 are excluded for the FPUT models, i.e. in the case m=0. In the Appendix we show that this implies that also for the case of m = 0 there are no such processes. It turns out that in the continuous limit there exist, apart from trivial resonances k 1 = k 2 = k 3 = k 4 or k 1 = k 3 , k 2 = k 4 , a nontirvial resonant manifold for interactions of the type 2 → 2 (see figure 1 for a visualization of the case m = 0, similar result apply for m > 0). The result has been obtained using the software Mathematica Wolfram.
The effective Hamiltonian and the kinetic equation for α, β-FPUT and DNKG models The effective Hamiltonian is obtained by removing all non resonant interactions from the original one. As standard in Hamiltonian systems, this can be carried out through a canonical change of variables (quasi-identity transformation), from the normal modes a k to some other variables b k , such that in the new variables the nonresonant terms are removed from the Hamiltonian, generating higher order nonlinearities. It is worth underlying that one has to ensure the canonicity of the transformation, at least up to the order of the new interaction terms that appear in the new variables. For example, to remove all the three-wave terms in the α-FPUT model one uses where the integral corresponds to three integrals from 0 to 2π and not from −∞ to +∞ and δ 2π account for periodicity of the Fourier space; this is because of the discreteness of our system in physical space. In eq. (34) we recognize terms equivalent to the nonlinear interactions in equation (18). The matrices X (i) 1,2,3 are determined by inserting eq. (34) into the α-FPUT Hamiltonian and by grouping the terms corresponding to different wave processes. The transformation generates four-wave interactions (and higher) and again one has to look for resonances and remove the non resonant ones. The calculation leads to the following effective Hamiltonian for the α-FPUT For the β-FPUT and for the DNKG a transformation is used to remove the 3→1, 1→3 and 4→0 terms, so that the effective Hamiltonian takes the form: For an interested reader, a comprehensive procedure for removing three and four-wave interactions, is presented in [12,26,28]. All three models are dynamically described by the Zakharov equation: where the coefficient W 1,2,3,4 takes different form, depending on the particular system under consideration. Its actual analytical form is irrelevant in our discussion (except in those cases for which the coefficient is zero on the resonant manifold [42,43]); it is just sufficient to mention that it is proportional to α 2 (A 1,2,3 ) 2 for the α-FPUT model and to βB (β,KG) 1,2,3,4 for the β-FPUT and DNKG models. The Zakharov equation is the starting point for developing the statistical theory, i.e. the wave kinetic equation, which takes the following form (see for details [9,14,32]): Note that, because of the δ (2π) , instead of δ, the momentum is not conserved. The time scales, T eq , of four wave resonant interaction, which lead to a thermalized state, are: for the α-FPUT model and for the DNKG and β-FPUT models. We we stress that the theory developed in this section is rigorously valid when the thermodynamic limit and the weak non-linear limit are taken, in the given order. It is important to wonder if the scaling laws obtained can be observed in actual finite-size systems and in particular in numerical experiments. We expect a positive answer, at least for sufficient large N and small nonlinearity. Indeed, it is well known that a small nonlinearity cause a frequency shift of the linear modes, and also a stocasticization of the frequencies, or in other words, a broadening [10,20]. It is reasonable to assume that resonances do not need to be satisfied exactly in practical applications [23]. It is sufficient that broadening of frequencies becomes comparable with the spacing of the frequencies, which decreases with the number of modes, so that ω k becomes continuous in the thermodynamic limit. In this sense, the discrete representation should converge to the continuous theory, and the higher the number of modes the smaller the nonlinearity required to be in agreement with the theory. We shall check this argument with extensive numerical simulations in the following.

DISCRETE EXACT RESONANCES
In many interesting cases, the number of modes N is not too large, like for instance in the original FPU numerical experiments. In this case, it is not possible to consider the system a good approximation of the continuous one, and it has to be studied in its discrete form. As mentioned, one should recall that in discrete systems the condition in the Kroneker δ over wave numbers should be intended mod N .

Four-wave resonant interactions, effective Hamiltonian and itegrability
The first problem is to find out if discrete exact resonances exist and at which order. The previous discussion on three-wave resonance still applies here, hence three-wave resonances are always excluded from the α-FPUT model. Concerning 4-wave resonances, the only process that is potentially active for the dispersion relation (8) is the scattering process 2 → 2. These resonances have been considered extensively in [4,35,36], and we will recap here the results. Obviously, even discrete processes of the type X → X admit trivial resonances of the type k 1 + k 2 = k 1 + k 2 mod N, ω 1 + ω 2 = ω 1 + ω 2 (41) or k 1 = k 2 = k 3 = k 4 . These processes, however, do not result into an exchange of energy, but rather cause a frequency shift of the linear modes [32]. Because the system is periodic and the dispersion relation is symmetric, it turns out that nontrivial resonances of the type 2 → 2 are possible with the crossing of Brillouin zones (Umklapp scatterings). For N even, these resonances take the form and are known as pairing-off resonances [4]. However, one can easily check that these resonances are all disconnected, and they actually give rise to integrable dynamics [18,37]. In [36] other special resonances have been considered for some very specific values of m = 0, but they are limited in number and they do not appear in general to be able to cover the whole Fourier space (a detailed study should be performed for such specific cases). The effective integrable Hamiltonian, up to four-wave interactions for the α, β and DNKG equations (with m different for those special values for which other resonant quartets exists), can be recast as follows: where c.c. denotes complex conjugate. The first three terms depend on the moduli of the amplitudes only and underline an integrable dynamics (Birkhoff normal form); interestingly, the last term does not break integrability (resonant Birkhoff normal form) [18,37], due to the relations (42) and the following symmetries: valid for all admissible values of k 1 , k 2 . It can be proven that N − 1 functionally independent invariants exist and are in involution (for α and β-FPUT see the proof of integrability in [18,37]). We present these invariants in a way that highlights the fact that the resonant quartets are disconnected: For each k = 1, . . . , N/4 , define an irreducible quartet as the set There are four different modes in Q k , except for the degenerate case k = N/4, valid when N/4 is integer, where Q N/4 has two different modes. In the non-degenerate case, the following four invariants depend on the four modes in Q k : • 3 quadratic invariants: • 1 quartic invariant: In the degenerate case, valid when N/4 is integer, we have and the following two invariants depend on the two modes in Q N/4 : • 1 quartic invariant: Thus, in terms of counting, any irreducible quartet (degenerate or not) contributes with a number of invariants that is equal to the number of modes in the quartet. Finally, notice that the mode b N 2 is not in any irreducible quartet. In fact, the following is an invariant: |b N 2 | 2 . It is thus easy to show by simple counting of the modes in the irreducible quartets that the system has a total of N − 1 functionally independent invariants: The result proofs the asymptotic integrability of the DNKG model up to four wave interactions and justifies the metastable states observed in numerical simulations previously performed, [16].

Higher order resonances and break-down of integrability
For the discussion above, we see that the 3-and 4-wave collision terms cannot be effective in bringing the system to equipartition. The isolated resonant quartets of the 4-wave integrable Hamiltonian do not bring the system in a thermalized state and resonances at higher order are to be investigated. For all discrete models considered here we have to perform an extra canonical change of variables, from the normal modes b k to some other variables c k , such that in the new variables the nonresonant terms are removed from the Hamiltonian. The question is now what is the leading order resonant wave interaction. In [35], only power-of-two values of N were investigated (akin to the original paper on the α-FPUT problem) and five-wave interactions were excluded on numerical grounds. From a more recent investigation [4], it turned out that five-wave resonant interactions exist only if N is divisible by 3 when m = 0. When N = 2 a 3 b with a, b > 1, the resulting five-wave clusters are connected and, in principle, a thermalized state can be reached, but this requires further study, to be discussed in a subsequent paper.
Excluding such specific values of N , six-wave resonant interactions are the leading order processes: it is always possible to find resonant six-wave tuples that are all interconnected and cover the whole Fourier space. These resonances are due to the same symmetries of ω k and are of the form: Such resonances are valid for even and odd values of N . Additional resonances can be found both in pairing-off form for even N , and recently also other resonances not in pairing-off form have been found [4], but they are not necessarily to cover the whole Fourier space and we will not discuss them. Another six-wave processes is theoretically possible, that is 4 → 2. Explicit formulas for 4 → 2 resonances have been found recently in [4], at least for N divisible by 3.
Excluding those specific values of N for which five wave interactions and 4 → 2 (or 4 → 2) exist, the six wave interaction Hamiltonian for the α, β-FPUT and DNKG equation can be written as where δ N is the Kroneker modulo N . The new matrix interaction can be calculated [28], but its precise form is not relevant to our scope.

The time scale of thermalization
The next step to be taken in order to evaluate the thermalisation time-scale should be the derivation of the corresponding kinetic equation, as done in the continuous case. Unfortunately, the kinetic equation can be rigorously derived only in the thermodynamic limit, that is N → ∞. Deriving a discrete version (N finite) of the kinetic equation poses significant mathematical problems [13].
Here, we make the conjecture that the time scale for thermalization in the discrete dynamics corresponding to the 6-wave interaction given by eq. (47) is at the leading order equivalent to the corresponding continuous 6-wave process. We also assume that the integrable part of the dynamics does not lead to any irreversible transfer of energy and the irreversible dynamics is fully contained in the six-wave interaction term. Our conjecture is therefore tantamount to saying that the thermalisation time-scale in the discrete case will be always proportional to the square of the coupling coefficient Z 1,2,3,4,5,6 in the six by Eq. (47), the latter being proportional to α 4 in the α-FPUT model and β 2 in the β-FPUT and DNKG system. We get the following estimates for the thermalisation time-scale for the α-FPUT model and for the β-FPUT, DNKG models. The argument used in the discrete case is always made for the limit of vanishing nonlinearity, where only exact resonances are important. However, the presence of a small but finite nonlinearity is expected to broaden the frequencies. In the case of a small number of modes N , that is in the present discrete case, an important difference with respect to the continuous case is expected. Since the wave-space is discrete, a large enough nonlinearity may cause a sufficient broadening to trigger resonances which are not in the exact-resonance manifolds and are found at a lower-order than the exact ones. Therefore, we can expect that at small N and sufficient nonlinearity the transfer of energy can be made through lower-order processes on smaller time-scales. A definite answer on the role of quasi-resonance is beyond the scope of the present manuscript.

Final remarks on the Wave Turbulent approach
Concluding the theoretical analysis, the WT formalism applied to the different anharmonic chains leads to consider the problem of transfer of energy as related to the collisions between waves in the system. In this framework, the energy is redistributed among modes if resonances exist which connect the whole wave-space. In the thermodynamic limit, the theory is continuous and is almost rigorous. It predicts that 4-wave processes are the leading order process in all chains and this fact allows a transparent estimate of the equipartition time-scale. In the case of a finite small N , we have considered the system intrinsically discrete and we have searched for the exact discrete resonances. We have seen that in this case, 4-waves processes are not able to transfer energy among modes, nor 5-wave collisions. The leading order is found to be 6-waves. Although we are not able to write the corresponding discrete kinetic equation, we have made the hypothesis that the same reasoning as in the continuous case can be followed, at least when looking at statistical observables like equipartion time-scale. In this way, an estimate of this time-scale is obtained always proportional to the square of the nonlinear coupling coefficient. Moreover, also when N is small, nonlinear effects can be important. We expect that the broadening of the frequencies due to nonlinearity can trigger lower-order process to transfer energy among modes when it becomes of the same order of the spacing in the frequency space, which is inversely proportional to N . The discrete case stands therefore on a less rigorous ground and the conjectures proposed can be only verified numerically a posteriori.

NUMERICAL SIMULATIONS
We now present the result of numerical simulations in support of the previous discussion. The equations (18) and (19) have been implemented in physical space, with a symplectic algorithm of the sixth order [39]. It is important to use a high-order integration scheme that preserves accurately the energy of the system, because simulation times are very long compared to the typical wave periods. The initial random energy per mode, e k , were drawn from an uniform distribution, and then scaled in order to obtain the desired total linear energy E and an initial number of particles N compatible with a final relaxation distribution with µ = 0 in eq. (30), in order to better observe equipartition. As an example, we show the initial distribution of e k for the simulations of the DNKG system with N = 64 in figure 2, together with its final thermalized state. The linear energy E and the wave action N change only for one part in 10 4 across the whole simulation time, being quasi-conserved quantities. Each realization in the ensemble shares the same initial energy per mode, but phases are randomized. This scheme of initialization ensures that all the realization share the same initial linear energy. The time step was set to 0.1 in all simulations, and it was checked that from beginning to the end of the simulation the total energy is conserved up to one part in 10 7 . The computation is very easy to parallelize, since ensemble averages need to be computed. For this reason we implemented the code on GPU hardware, which is powerful when the code can be parallelized.
The parameter α and β are chosen in such a way that the nonlinear energy is small compared to the linear one and the parameters α and β,KG are small. We will present the results of T eq as a function of α and β,KG . The approach to equipartition is often monitored with the information entropy, where In general, the information entropy is not guaranteed to have entropy-like properties, that is monotonic increase or decrease over time. Nevertheless, it has the useful property that at perfect equipartition S inf = 0 and it is greater than zero in any other state. We argue that when µ = 0 the information entropy is essentially equivalent to the WT entropy (29). The WT entropy has the remarkable property that one can replace n k in eq. (29) with f (k)n k , with f (k) being any function of k constant in time and still it maintains the property dS/dt ≥ 0. We exploit this and use f (k) = N/Eω k , so that the WT entropy now is computed as This expression preserves the monotonicity of the WT entropy with an inverted sign (dS WT /dt ≤ 0, S WT = 0 at equipartition). Since log(e k ) 0 close to equipartition, we see that in practice the two entropies (50) and (52) are equivalent for the purpose of determining the equipartition of the system. In numerical simulations, the entropy shows some fluctuations and never reaches exactly zero, because the simulation time is finite, but most importantly because the ensemble size is finite. In the following, we will determine T eq looking at the time when the entropy crosses a threshold value from above. The selection of the threshold depends on the ensemble size and the number of particles, but it does not depend on the degree of nonlinearity. Generally we chose a value of S WT two orders of magnitude smaller than that of the initial conditions.

Small number of particles, discrete resonances
We first present the results on the case of a limited number of particles, N = 64, so that discrete resonances need to be considered. To recap, we expect for the α-FPUT system the scaling and for β-FPUT and DNKG systems the scaling We show the results of the simulations for N = 64 in figures 3, 4, 5. We see a good alignment with the expected power-laws for all the three models. We chose m = 1 in the DNKG model arbitrarily: in [36] it is argued that since the resonances in the discrete case do not depend on m, then in the scope of the thermalization dynamics the value of m is not relevant. However, we should add here a comment for the limit of m very small. In such cases, the contribute to H nlin increases for the low frequency modes. On the other hand, for the FPUT systems the nonlinearity depends on a power of the finite differences, which are approximately proportional to k, hence the contribute to the nonlinear part of the Hamiltonians for modes with small k does not grow. As a consequence, in the case of the DNKG system with m = 0, a complete thermalization of the lower modes would invalid the weak-nonlinear assumption, because the nonlinearity would be too high. What is observed actually (see figure 6) is a partial thermalization of the system, with the lower modes stationary at a lower energy than equipartition.

Large number of particles, continuous resonances
We now present the results of numerical simulations for the continuous resonance manifold. The expected results are for the α-FPUT system, and   are available in [30,32], and they can be compared to the typical frequency spacing for a fixed number of particles. We point out that in general a higher value for m in eq. (8) makes the frequency gaps smaller among modes, hence the DNKG model in general requires a small number of particles to observe the continuous resonances. In figure 7,8, 9 we show the results for simulations with N = 1024, again with m = 1 for the DNKG system. The figure highlight a reasonable agreement between numerics and the theoretical predictions (the exponent in the α-FPUT seems to be slightly steeper than predicted, probably because the thermodynamic limit has not been reached numerically yet). Therefore, from the simulation data we can deduce that the WT picture works reasonably well in all cases, at small and large N .

CONCLUSIONS
In this work, we have presented recent results on the dynamics of nonlinear chains, focusing on the thermalisation time. We have proposed to encompass the behaviour of all systems of this kind within the Wave Turbulence framework. In this context, the universal mechanism invoked to lead the systems to equipartition is the resonance among different modes, at least when generic initial conditions are considered. The scaling of the thermalization time T eq as a function of the nonlinearity level of the system is obtained from the theory. The result is a power-law for all the cases considered, with an exponent dependent on the order of the active resonances in the system. Although the point if view is rather different, present results are not in contradiction with the most recent and accurate estimates obtained also for particular low-wavenumber initial conditions [3]. Since the resonant manifold is significantly different in the discrete and continuous N → ∞ limit, we have considered the two cases separately, and they lead to two different scalings. Notably, large-size chains reach equipartition faster. Our aim in this paper has been to present the results anticipated in works [30,35,36] in an unified manner to underline how our approach can be systematic, rather than dependent on the particular features of the systems. In the continuous limit, we have shown that there exists a resonant manifold of 2 → 2 waves. We have analysed numerically the continuous limit N → ∞ for the α and β-FPUT systems, as done previously for the DNKG chain. In the discrete case, for N a power of two (which excludes the existence of five-wave resonant interactions [4]), we have extensively verified that, for all the considered systems, the leading order interaction consists of six wave interactions. All the results seem to indicate a universal route to thermalization predicted. We have also shown that, besides the α and β-FPUT systems, also the DNKG equation is asymptotically integrable, if the expansion is truncated at four wave interactions.
One should note that the WT theory predicts thermalization for arbitrary small nonlinearity, with always the same power law scaling. It would be interesting to understand how this is related to the result of Nekhoroshev [33] who demonstrated that solutions stay close to their integrable counterparts for an exponentially long time.
The function f k is m-independent, while ω k with m = 0 is proportional to m. The terms of the type −ω 2 ki all cancel out. It is easy to observe that F (k 1 , k 2 , .., k X ; m) is an increasing monotnic function of m. For m = 0, because of the inequality in (57), the F (k 1 , k 2 , .., k X ; 0) ≥ 0; then, for its monotonicity, for any m > 0, F (k 1 , k 2 , .., k X ; 0) > 0. This proves that there are no resonances in the DNKG model of the type X → 1.