Stages of dynamics in the Fermi-Pasta-Ulam system as probed by the first Toda integral

TL;DR

This paper studies the long-term behavior of FPU system trajectories using the Toda integral, revealing two key timescales and how they scale with energy and system size, providing insights into FPU dynamics and thermalization.

Contribution

It introduces the use of the first Toda integral as a fluctuation-free observable to analyze FPU dynamics and identifies scaling laws for stability and equilibration times.

Findings

Two fundamental timescales: stability and equilibrium.

Exponential scaling of times with inverse energy.

Crossover energy scales as a power law with system size.

Abstract

We investigate the long term evolution of trajectories in the Fermi-Pasta-Ulam (FPU) system, using as a probe the first non--trivial integral $J$ in the hierarchy of integrals of the corresponding Toda lattice model. To this end we perform simulations of FPU--trajectories for various classes of initial conditions produced by the excitation of isolated modes, packets, as well as `generic' (random) initial data. For initial conditions corresponding to localized energy excitations, $J$ exhibits variations yielding `sigmoid' curves similar to observables used in literature, e.g. the `spectral entropy' or various types of `correlation functions'. However, $J(t)$ is free of fluctuations inherent in such observables, hence it constitutes an ideal observable for probing the timescales involved in the stages of FPU dynamics. We observe two fundamental timescales: i) the `time of stability' (in which, roughly, FPU trajectories behave like Toda), and ii) the `time to equilibrium' (beyond which energy equipartition is reached). Below a specific energy crossover, both times are found to scale exponentially as an inverse power of the specific energy. However, this crossover goes to zero with increasing the degrees of freedom $N$ as $\varepsilon _c \sim N^{-b}$, with $b \in [1.5, 2.5]$. For `generic data' initial conditions, instead, $J(t)$ allows to quantify the continuous in time slow diffusion of the FPU trajectories in a direction transverse to the Toda tori.

Equations

$$\displaystyle H_{FPU}=\sum_{n=0}^{N-1}\big{[}{1\over 2}p_{n}^{2}+{1\over 2}(q_{n+1}-q_{n})^{2}+{\alpha\over 3}(q_{n+1}-q_{n})^{3}\big{]}$$

$$\displaystyle H_{T}=\sum_{n=0}^{N-1}\big{[}{1\over 2}p_{n}^{2}+\frac{1}{4\alpha^{2}}(a_{n}^{2}-1)\big{]},$$

$$\displaystyle J={1\over N}\sum_{n=0}^{N-1}\big{[}{\alpha^{2}\over 2}p_{n}^{4}+\frac{1}{2}a_{n}^{2}(p_{n}^{2}+p_{n}p_{n+1}+p_{n+1}^{2})+\frac{a_{n}^{2}}{16\alpha^{2}}(a_{n+1}^{2}+a_{n}^{2}+a_{n-1}^{2})-\frac{3}{16\alpha^{2}}\big{]}$$

$$\tilde{J}(t)=\frac{J(t)-J(0)}{J_{eq}-J(0)}.$$

$$\Delta (t)=\frac{\overline{{\sum }_n{q}_n{q}_{n+1}}}{\overline{{\sum }_n{q}_n^2}}$$

$$\displaystyle J\simeq 2\varepsilon+{1\over{2N}}\sum_{n}\Big{[}p_{n}p_{n+1}+\delta q_{n}\delta q_{n+1}\Big{]}~{}~{}.$$

$$\displaystyle J\simeq 2\varepsilon+{1\over{N}}\sum_{k}\cos\big{(}\frac{\pi k}{N}\big{)}E_{k}~{}~{}.$$

$$\eta =2\sum _{k\ge N/2}{E}_k/E$$

$$f\propto \exp \left(-\frac{\beta }{2}\sum _{k=1}^N(P_k^2+{\omega }_k^2{Q}_k^2)\right).$$

$${\sigma }_J^2(t)=\frac{1}{M}\sum _{l=1}^M{\left(J_l(t)-J_l(0)-{\mu }_J(t)\right)}^2$$

$$t_0(N,\epsilon )=t_0(\epsilon )\simeq c_1\tau$$

$$t_{eq}(N,\epsilon )=t_{eq}(\epsilon )\simeq c_2\tau$$

$$QT=\frac{1}{N}\sum _n\frac{1}{2}(p_n^2+p_n{p}_{n+1}+p_{n+1}^2)+\frac{1}{8}[(q_{n+2}-q_n{)}^2+4(q_{n+1}-q_n{)}^2+(q_{n+1}-q_{n-1}{)}^2]$$

$$QT\approx 2\epsilon +\frac{1}{2N}\sum _n(p_n{p}_{n+1}+\delta q_n\delta q_{n+1})$$

$$QT\approx 2\epsilon +\frac{1}{N}\sum _k\cos (\frac{\pi k}{N})E_k.$$

Full text

1 Introduction

The numerical experiment of Fermi, Pasta and Ulam [20, 18] in 1954 aimed to probe ergodicity in an one–dimensional chain of $N$ weakly nonlinearly coupled oscillators. The discovery of FPU recurrences led to a substantial open question in Statistical Mechanics: Does a system with $N$ degrees of freedom, stemming from a generic perturbation to an integrable model, always tend to an ergodic state as $N$ becomes large? Several numerical studies confirm ergodicity in the FPU problem for large $N$, and for various special types of initial conditions [4, 30, 31]. Nevertheless, despite a vast effort we are still far from reaching a rigorous answer to the above question. One obstruction to rigorous results stems from difficulties in implementing perturbation theory in the ‘thermodynamic limit’, i.e. when the energy per oscillator $\varepsilon$ (specific energy) is small and fixed while $N$ becomes large.

Detailed reviews on the FPU problem can be found in [2, 27]. We quote here some important numerical and theoretical main lines of approach to the FPU problem, which are related to our present study:

i) Departure from the harmonic limit: the lowest order integrable approximation to the FPU is the chain of uncoupled oscillators constituting its linear normal modes. Initial conditions ‘à la Fermi’, i.e., low–frequency normal mode excitations, have been studied extensively [1, 3, 4, 21, 22, 23, 28, 36, 37, 38]. They are well known to lead to exponentially localized energy profiles in the space of normal modes. Such profiles persist for very long times (‘metastable states’), but eventually evolve to states closer to energy equipartition (‘equilibrium states’). In the harmonic limit (small specific energies) this behavior can be interpreted through the stability properties of particular low-dimensional invariant objects of the FPU phase space. Most notably, the ‘$q$–breathers’ are Lyapunov periodic orbits forming the continuation of the linear modes [21, 22, 34]. On the other hand, the ‘$q$–tori’ correspond to the extension, in the nonlinear regime, of quasi-periodic motions pertinent to excitations of packets of normal modes in the linear regime [12, 13, 14]. As far as the $q$–breathers or $q$–tori are stable (e.g. with respect to transverse perturbations), the metastability phenomenon can be understood as a case of ‘stickiness’ around these objects. From this viewpoint, however, it remains puzzling that the metastability phenomenon persists for specific energies higher than the stability thresholds of $q$–breathers or $q$–tori [21, 22, 12, 14], and clearly beyond the harmonic regime (although of course still small, i.e. $\varepsilon<1$). One should note here that slow in time deviation from a quasi-integrable behavior seems to occur also in cases of more general types of initial conditions, where, from the beginning, the energy is distributed in the whole spectrum, instead of isolated packets of modes. In this case, we can compute the evolution of the autocorrelation functions for suitably defined phase space quantities [7, 8, 9, 10]. The equilibrium state is identified as the point of complete decay of the autocorrelation coefficients.

ii) Stochasticity threshold: In 1959 Chirikov published his celebrated works [15, 16, 26] on the existence of a stochasticity threshold in terms of specific energy in the FPU, which separates chaotic from weakly chaotic or regular motions. Chirikov’s threshold $\varepsilon_{s}$ is derived by computing conditions under which the FPU resonances exhibit substantial overlapping. He finds that $\varepsilon_{s}$ vanishes with $N$ like $\varepsilon_{s}\sim 1/N^{4}$. In [40], Shepelyansky extended Chirikov’s results by including the case of low–mode excitations. Thus, the general conclusion from such studies is that one should expect stochasticity to prevail, and metastability phenomena to disappear, in the thermodynamic limit. Similar conclusions are reached when considering the dependence of the stability thresholds on $N$ for the $q$–breathers [21, 22] and a weaker decay of the form $1/N$ for the $q$–tori [14].

However, as noted in [40], any approach based on general criteria of resonance overlap cannot exclude the possibility that the system under investigation is very close to an integrable one, in which case estimates on resonances do not apply. In the words of Shepelyansky: ‘This point is very crucial for the $\alpha$–FPU problem, since at low energy it is very close to the Toda lattice. Due to that generally we should expect that, in contrast to the above estimates and numerical data the dynamics of the $\alpha$–FPU problem will be integrable’. In recent years, there has been increasing attention to Toda as a reference integrable model close to the FPU [5, 6, 14, 25, 37]. Already in 1982 the pioneering work of Ferguson et al. [19] put forward the idea that the higher polynomial order contact between FPU and Toda allows to associate the FPU’s integrable–like behavior (e.g. FPU recurrences) with Toda. In fact, this can be regarded as a ‘discrete’ analogue of the Zabusky–Kruskal approach, which associates the FPU with a different integrable continuous limit, i.e., the Korteweg–de Vries (KdV) equation [32].

In the present paper, we propose a simple method to measure the proximity of (evolving) FPU states, generated by various types of initial conditions, to the dynamics of the nearby integrable Toda model: this is to observe directly the evolution of the functions $J(q,p)$ yielding the Toda integrals, for a single or an ensemble of FPU trajectories. An independent work along the same direction appeared recently in [25]. In the sequel, we focus on the evolution of only the first non-trivial Toda integral (besides the energy), denoted hereafter as $J$. One has a constant value $J(t)=c$ for any trajectory under the exact Toda dynamics, while $J(t)$ varies along the FPU trajectories. As shown below, these variations allow to characterize the FPU stages of dynamics according to their proximity to the Toda dynamics. In particular, one can clearly identify ‘stickiness’ effects, and measure the times up to which the FPU trajectories remain sticky to nearby Toda tori.

Examples of this behavior, along with a corresponding quantitative analysis, are given in cases covering most widespread categories of initial conditions encountered in FPU literature. These roughly cover three classes of initial conditions: i) single mode or packet of modes excitations with coherent or random phases, ii) random data, where the whole energy spectrum has comparable power from the start, and iii) generic data close to energy equipartition. Cases (i) and (ii) are characterized by distinct phases of evolution of the FPU trajectories, which roughly correspond to the stages of dynamics recognized in [1, 37, 5]. We show below how $J(t)$ allows to measure the timescales related with the stages of dynamics as well as the latter’s dependence on the system’s specific energy and number of degrees of freedom. In case (iii), using $J(t)$ we provide direct evidence of diffusion taking place as the FPU trajectories wander in phase space transversally to ‘Toda tori’. This is useful also in interpreting the observed decay of the autocorrelation for the Toda integrals as reported in [5].

The paper is organized as follows: Section 2 deals with definitions and examples related to the utility of $J$ as an observable for featuring FPU–evolution. Section 3 analyses the information obtained by $J$ in numerical experiments covering each of the aforementioned class of initial conditions. Section 4 contains the basic conclusions of the present study.

2 The first Toda integral as an FPU–observable

Fermi, Pasta and Ulam [20] considered the one–dimensional lattice consisting of $N$ weakly nonlinearly coupled oscillators. The dynamics of this model is described by the Hamiltonian:

[TABLE]

where $q_{n}$ is the $n$–th particle’s displacement with respect to equilibrium and $p_{n}$ its canonically conjugate momentum. Fixed boundary conditions are imposed for $q_{0}=q_{N}=p_{0}=p_{N}=0$.

The normal modes of the harmonic limit $\alpha=0$ are defined by the canonical transformation $(Q_{k},P_{k})=\sqrt{2\over N}\sum_{n=1}^{N-1}(q_{n},p_{n})\sin\left({nk\pi\over N}\right)$. The normal mode energies are $E_{k}={1\over 2}(P_{k}^{2}+\omega_{k}^{2}Q_{k}^{2})$ with frequencies $\omega_{k}=2\sin(k\pi/2N)$, $k=1,\ldots,N-1$.

The FPU Hamiltonian has a third–order polynomial contact with the integrable Toda Hamiltonian [41]

[TABLE]

where $a_{n}=e^{\alpha(q_{n+1}-q_{n})}$. Taylor–expanding the exponential $a_{n}$, one obtains $H_{T}=H_{FPU}+{\cal O}([\delta q_{n}]^{4})$, where $\delta q_{n}=q_{n+1}-q_{n}$ are the relative displacements.

The expressions of all the Toda integrals are given in [24, 29]. Besides $H_{T}$, the first non–trivial Toda integral after adding and dividing with some constants, takes the form:

[TABLE]

where fixed boundary conditions are imposed for $a_{0}=a_{-1}=e^{\alpha q_{1}}$, $a_{N}=a_{N-1}=e^{-\alpha q_{N}}$.

We outline the behavior of $J$ along FPU trajectories through the example of Fig.1. Consider the excitation of the first normal mode for the FPU system with $\alpha=1/2$, $\varepsilon=0.01$ and $N=8192$ particles. As will be described below, independently of the initial conditions considered, we find that, after sufficiently long integration, $J(t)$ always tends to stabilize (apart from small fluctuations) to a final value, called hereafter, the ‘equilibrium value’ $J_{eq}$. In practice, $J_{eq}$ is numerically computed by taking averages of $J(t)$ in moving time windows long enough to absorb fluctuations. One notes that eventually the moving average stabilizes to a constant value after $t=t_{eq}$. For $t>t_{eq}$, $J(t)$ only exhibits rapid and irregular fluctuations around the constant average $J_{eq}$. We argue below that this is an indication for a statistical equilibrium state reached by the system. The time $t_{eq}$ is hereafter called the equilibrium time. The normalized $J$ is defined as:

[TABLE]

In Fig.1 $\tilde{J}$ displays a sigmoidal evolution from 0 to 1. The detachment of $\tilde{J}$ from zero, which indicates an essential departure from quasi-integrable behavior, takes place at a time $t_{0}$. During the time interval $[0,t_{0}]$ the behavior of $J(t)$ for FPU and Toda is nearly indistinguishable. The time $t_{0}$ is hereafter called the time of stability. At longer times, $\tilde{J}$ tends sigmoidally towards the value 1, reached at time $t_{eq}$.

It is crucial to clarify in which sense $t_{eq}$ reflects the time in which statistical equilibrium has been reached. We give numerical evidence in Fig.2, and we theoretically justify below, that the stabilization of $J$ to $J_{eq}$ for $t>t_{eq}$ is related to the following: (i) the variables $p_{n}$, $\delta q_{n}$ decorrelate (Fig.2(b),(f)), (ii) the energy spectrum $E_{k}$ reaches equipartition (Fig.2(d)).

Regarding (i), Fig.2(b) displays the evolution of the sums $\sum p_{n}p_{n+1}$, $\sum\delta q_{n}\delta q_{n+1}$ for the same orbit as in Fig.1. Similar behavior is found for all orbits with initial conditions corresponding to localized excitations. For $t<t_{0}$ both sums oscillate, with maxima of the one sum corresponding to minima of the other. Between $t_{0}$ and $t_{eq}$, both sums tend sigmoidally to zero, while, beyond $t_{eq}$ they become uncorrelated and both fluctuate irregularly around zero. Let us note, in respect, that a correlation sum similar to the above was considered by Parisi [33]

[TABLE]

as an accurate and simple observable to be used in specifying the equilibrium time.

We will now show, that property (ii) above is connected with the correlation sums $\sum p_{n}p_{n+1}$, $\sum\delta q_{n}\delta q_{n+1}$ both tending to zero for times beyond $t_{eq}$. To this end, for specific energies much smaller than unity, $J$ is well approximated by expanding the exponentials in (3) in Taylor series in powers of the quantities $p_{n}$, $\delta q_{n}$. Up to quadratic terms we get (see Appendix A):

[TABLE]

Due to Eq.(5), the near–constancy of $J$ for $t<t_{0}$ in Fig.2(a) emerges from the counterbalance of the sums $\sum p_{n}p_{n+1}$ and $\sum\delta q_{n}\delta q_{n+1}$ in Fig.2(b), which oscillate around a non–zero mean with nearly opposite phases. For $t>t_{eq}$, however, they both tend to zero. In fact, as shown in Appendix A, Eq.(5) can be recast as a sum over the energy spectrum $E_{k}$:

[TABLE]

This indicates that $J$ nearly corresponds to a linear combination of the energies $E_{k}$, where more weight is given to the energies $E_{k}$ corresponding to low and high frequency modes, and less weight to modes in the middle of the spectrum ($k\approx N/2$). Most notably, the sigmoid transition of $J$ from the value $J(0)$ to $J_{eq}$ is connected with the evolution of the energy spectrum $E_{k}$ from a localized one (for $t<t_{0}$) to near-equipartition $E_{k}\simeq\varepsilon$ for times $t>t_{eq}$. Since $\sum_{k}\cos(\pi k/N)=0$, one gets $J_{eq}\simeq 2\varepsilon$, implying that $t_{eq}$ can also be interpreted as the ‘equipartition time’. On the other hand, the distribution of the momenta for $t>t_{eq}$ approaches a Gaussian with dispersion $\sigma_{p}^{2}=2\varepsilon$ (cf. Fig.2(e) and (f)). This is not far from the Gibbs measure for an equilibrium state with Hamiltonian (1). Thus, while energy equipartition alone does not necessarily imply an equilibrium state (see also numerical simulations below), in the above example energy equipartition is linked to and appears at the same timescale as the decorrelation between the phase space variables.

It is noteworthy that approximating formulas analogous to Eq.(6) can be derived relating all the Toda integrals with linear combinations of the spectral energies $E_{k}$. Thus, a behavior analogous to the one found above for $J$ is expected to hold for all Toda integrals. This subject is currently under investigation.

Finally, we remark that the evolution of $\tilde{J}$ in Fig.1 is nearly indistinguishable from the evolution of the tail energy 111 A idea similar to the tail energy for generic packet excitations can be traced back to reference [35].

[TABLE]

which has been proposed as an efficient observable for distinguishing stages of dynamics in the FPU problem (see [37, 4]). An additional method based on the times at which the FPU trajectories intersect an ‘equilibrium manifold’ [17] yields results similar to those of [37] found using the tail energy evolution. We find that the normalized Toda integral $\tilde{J}$ and the tail energy practically coincide for trajectories with initial low-mode excitations. However, as discussed below, $\tilde{J}$ is of use also in more generic initial conditions, in which the tail energy cannot be used.

3 Numerical Experiments

In this section we numerically investigate the stability, diffusion, and approach to equilibrium for FPU trajectories with various types of initial conditions, using as a probe the evolution of $J$.

3.1 Classes of initial conditions and the evolution of J

We cluster different types of initial conditions into three main categories referred to below as: (i) with ‘localized energy spectrum’, (ii) with ‘delocalized energy spectrum’, and (iii) ‘close to equipartition’.

(i) Localized energy spectrum: Case (i) contains FPU trajectories with initial conditions corresponding to single–site excitations, and packets of modes with coherent or random phases. The initial conditions are controlled through the initial values assigned to the normal mode variables $(Q_{k},P_{k})$, $k=1,\ldots,N-1$. To an initial energy $E_{k}(0)$ assigned to mode $k$ there corresponds a family of possible initial conditions with $Q_{k}(0)=A_{k}\sin(\phi_{k})$, $P_{k}(0)=\omega_{k}A_{k}\cos(\phi_{k})$ where $A_{k}=(2E_{k}(0)/\omega_{k})^{1/2}$ and $0\leq\phi_{k}<2\pi$. A single–site excitation is obtained by setting $A_{k}=0$ for all $k$ except one $k=k_{0}$. Here we consider the case $k_{0}=1$. A ‘packet of modes’ excitation of percentage $w$ corresponds to setting $E_{k}(0)=N\varepsilon/w$ for all $k$ in $1\leq k\leq[wN]$ with $0<w\leq N$ and $E_{k}(0)=0$ otherwise. In this case, we can have initial phases ‘coherent’ ($\phi_{k}=const$, here $\phi_{k}=0$), or ‘random’ (here $\phi_{k}$ is chosen from a uniform distribution in $[0,2\pi)$).

Figure 3(a) shows the evolution of $J$ along FPU trajectories with initial conditions corresponding to packets of various sizes (values of $w$ as denoted in the figure) and random phases, while the specific energy is kept fixed at $\varepsilon=0.01$. The $12.5\%$ packet (black) gives rise initially to a strongly localized energy spectrum which is nearly flat among the packet modes and develops an exponential tail for the rest of the modes, as shown in Fig.3(b). By progressively increasing the width of the packet (from $12.5\%$ to $87.5\%$), different localization patterns emerge in Fig.3(b). However, the evolution of $J$ in Fig.3(a) shows that equilibrium is reached at nearly equal times $t_{eq}$ for all these trajectories. The dependence of both the ‘stability time’ $t_{0}$ and the equilibrium time $t_{eq}$ on $N$ and $\varepsilon$ is discussed in subsection 3.2.

(ii) Delocalized energy spectrum: case (ii) refers to non–equilibrium initial conditions with substantial power in the whole normal mode energy spectrum $E_{k}$. As a basic example, we consider random initial positions and momenta $(q_{n},p_{n})$ extracted from a uniform distribution in the intervals $-p_{\varepsilon}\leq p_{n}\leq p_{\varepsilon}$, $-q_{\varepsilon}\leq q_{n}\leq q_{\varepsilon}$, with $p_{\varepsilon}=q_{\varepsilon}$ tuned numerically so that the total energy takes a selected value $E=N\varepsilon$. Fig.4(a) shows the initial energy spectrum $E_{k}$ for such a choice of initial conditions, leading to specific energy $\varepsilon=0.02$. At $t=0$ the system appears as not being too far from energy equipartition, however, such initial conditions are not only far from the Gibbs measure but need longer equilibrium times than low–mode excitations. As shown in Fig.4(b) for a single realization, $J(t)$ exhibits a similar sigmoidal evolution to the case (i). The approximations (5) and (6) yield curves nearly indistinguishable from the exact curve $J(t)$. Thus, the sigmoidal evolution is related to the slight redistribution of energies taking place as the system evolves to equilibrium. These are hard to distinguish using measures based directly on the spectrum $E_{k}$, as for example, the spectral entropy $S$; see [38]). However $J(t)$ measures with accuracy the time to equilibrium $t_{eq}\approx 10^{7}$. Remarkably, $t_{eq}$ for random initial positions and momenta turns to be of the same order as for the packet initial conditions, and actually somewhat larger than the times $t_{eq}$ found for any ‘à la Fermi’ type of initial condition. This implies that stickiness to the Toda dynamics is a generic property not restricted to trajectories with localized initial energy spectrum.

(iii) Initial conditions close to equipartition: here we consider two subcases of initial conditions, namely far from or close to the Gibbs distribution $f(\mathbf{q},\mathbf{p})\propto e^{-\beta H_{FPU}(\mathbf{q},\mathbf{p})}$. In the first subcase, we consider initial conditions called hereafter ‘random momenta’, in which we simply set $q_{n}=0$, and $p_{n}$ randomly chosen from a uniform distribution in the interval $-p_{\varepsilon}\leq p_{n}\leq p_{\varepsilon}$, with $p_{\varepsilon}=\sqrt{6\varepsilon}$. One obtains $<P_{k}^{2}>=<p_{n}^{2}>=2\varepsilon$, and hence a spectrum $E_{k}$ with $<E_{k}>=\varepsilon$. As an alternative, we fix the normal mode variables $P_{k}=\omega_{k}A_{k}\cos\phi_{k}$, $Q_{k}=A_{k}\sin\phi_{k}$, with $A_{k}=(2\varepsilon/\omega_{k}^{2})^{1/2}$ and $\phi_{k}$ randomly chosen with uniform distribution in the interval $[0,2\pi)$. Both these types of initial condition lead to spectra close to equipartition, but far from the distribution function associated with statistical equilibrium. Instead, in the second subcase, hereafter called ‘close to equilibrium’, as in the works [5, 6, 7, 8] we approximate a Gibbs distribution function $f$ for the normal mode variables $(Q_{k},P_{k})$ by considering only the quadratic part of the Hamiltonian $H_{FPU}$. The function $f$ becomes a Gaussian

[TABLE]

Setting $\beta=\varepsilon^{-1}$ ensures the specific energy is equal to $\varepsilon$ for any random realization in the above distribution.

One key feature of the numerical experiments with all the above initial conditions is that $J(t)$ (for single trajectories, or averaged over trajectories with different realizations of the initial conditions) exhibits no systematic change of its value (e.g. of sigmoid form) allowing to define a timescale to equilibrium as in the case of the experiments in class (i) and (ii). Fig.5(a) and (b) shows two examples of the evolution of the average $J(t)$ for ten realizations of trajectories with initial conditions belonging to the ‘random phase’ and ‘random momenta’ subclass. Despite the absence of sigmoid evolution, one may still argue that for initial conditions far from equilibrium, a certain time is required for these trajectories to drift substantially in phase space so as to approach a state of equilibrium. Such an approach cannot be detected by the evolution of the energy spectra $E_{k}$, since the latter are set from the start close to energy equipartition. Nevertheless, the form of the curves in Fig.5 suggests the trajectories undergo diffusion in the direction normal to the integral surfaces defined by the Toda integrals, or at least the first one. Such diffusion can be detected by measuring the evolution of the dispersion in the values of $J(t)$ over $M$ trajectories

[TABLE]

where $\mu_{J}(t)={1\over M}\sum_{l=1}^{M}\left(J_{l}(t)-J_{l}(0)\right)$. Subtracting the initial value $J_{l}(0)$ is necessary to absorb the spreading in the initial values of $J_{l}$ due to the random generator of initial conditions. Fig.6a Shows $\sigma^{2}(t)$ as a function of $t$ for a set of 30 trajectories with ‘random initial momenta’. The solid line has logarithmic slope equal to 1, indicating a normal diffusion law $\sigma^{2}\propto t$. Notice that the diffusion stops after a time $t\sim 10^{6}$. The spreading in the values of $J$ after this time has measure $O(\varepsilon^{2})$, indicating that the trajectories have spread over parts of the phase space covering the whole possible range of values of $J$ consistent with a fixed specific energy equal to $\varepsilon=0.01$.

On the other hand, it is really remarkable that a similar spreading is found for initial conditions very close to the Gibbs measure, as exemplified in Fig.6(b). Now, the initial conditions are chosen by the distribution function of Eq.(8), thus they represent a state as close as possible to statistical equilibrium already at the starting point of the simulation. Yet, we observe that the underlying integrable dynamics of the associated Toda model leaves its traces in this case too, since the trajectories diffuse transversally to the integral surface of the integral $J$ with a very slow speed, comparable to the one in initial conditions far from equilibrium. Measuring this speed at various energy levels, as well as for different Toda integrals represents a challenging numerical task, since it requires many trajectory realizations per parameter set considered in order to ensure a good statistics (setting $M=30$ in the above experiments is marginal in this respect). Instead, the speed of approach to equilibrium is measured more easily via the sigmoid curves in experiments of classes (i) and (ii), a computation to which we now turn our attention.

3.2 Intensive quantities and a single timescale

We first investigate for which types of initial conditions the behavior of $J(t)$ and associated timescales (of ‘stability’ or ’approach to equilibrium’, see section 2) exhibit a dependence on $N$. It turns out that only random initial data lead to $N$–independent results. For example, Fig.7(a) shows the sigmoid curves when considering the excitation of the first normal mode (case (i)), with $N$ ranging from 512 to 8192. The curves are distinguished, in particular as regards the ‘stability’ times characterized by the onset of the rising part of the sigmoid evolution, even if the equilibrium times are similar for all $N$. In comparison, Fig.7(b) shows the same computation but for initial conditions extracted by random positions and momenta (case (ii)). Superposing the sigmoid curves for various $N$ shows a near-coincidence at all times.

We note that a similar behavior to Fig.7(b) holds for initial data as those of Fig.3. Thus, packet excitations with random phases exhibit $N$–independence in the timescales associated with the evolution of $\tilde{J}(t)$. This is in agreement with the results first reported in [3], where the authors aim to unify ‘incompatible’ earlier findings in the FPU literature (see [3] and references therein) on the dependence of various scaling laws on the energy $E$, when exciting packets of modes with coherent phases, or the specific energy $\varepsilon$, when the initial phases are random. It is also in agreement with the work [14], where it is found that extensive packets with coherent phases give rise to exponentially localized energy spectra of the form $E_{k}\sim e^{-\sigma k/N}$, with $\sigma\sim(\alpha^{2}E)^{-d}$, $d>0$, i.e., depending on $E$ rather than $\varepsilon$.

Exploiting the property of $N$–independence, we can characterize the timescales involved in the evolution of $J(t)$ for classes of trajectories based on some form of random initial data. Computing $J(t)$ allows to obtain a good estimate of the ‘stability time’, from the size of the initial plateau of the corresponding sigmoid curve, and the time of approach to equilibrium, when the second plateau is reached. Due to the numerical indications for $N$–independence one may reasonably argue that the results hold in the thermodynamic limit ($N\rightarrow\infty$ keeping $\varepsilon$ constant). Focusing, as an example, on trajectories with random initial positions and momenta, the sigmoid curves for the normalized Toda integral $\tilde{J}(t)$ are reproduced in Fig.7(c). With an appropriate time–shift $\tau\sim\varepsilon^{-2.68}$ (inset), those curves fall one onto the other as shown in Fig.7(d). We typically find $\tau\sim\varepsilon^{-a}$ beyond a crossover specific energy (see below) for some $a>0$ depending on the type of initial conditions. Since the shift in the sigmoid curves takes place in a logarithmic time axis, the above facts imply that $\tau$ allows to characterize the scalings with $\varepsilon$ of both the initial ‘stability time’ and the ‘time to equilibrium’. One has

[TABLE]

with $c_{1}<1<c_{2}$. At practical level, the possibility to extrapolate scaling laws from the ‘stability’ to the ‘approach to equilibrium’ phase allows to translate results found for $t_{0}$ to analogous results for $t_{eq}$, even when $t_{eq}$ is much larger than $t_{0}$ and hence hard to compute by direct numerical experiments. One has to compute in respect the initial part of the curve $J(t)$, up to the detachment from the first plateau. We also note that the numerical complexity for computing $J$ is $O(N)$, which is much better than the $O(N^{2})$ complexity of any computation requiring knowledge of the evolution of the energy spectra $E_{k}$.

All the above results hold asymptotically for $N$ sufficiently large. Instead, for $N$ small we may find deviations from the law $\tau\sim\varepsilon^{-a}$ which are $N$–dependent. In particular, we find transient exponential laws $\tau\sim\exp\varepsilon^{-b(N)}$, $b(N)>0$ up to a crossover specific energy $\varepsilon_{c}(N)$. However, we can make use of $J$’s asymptotic independence in order to locate the specific energy crossover $\varepsilon_{c}(N)$. The numerical procedure is the following: the sigmoid curves $\tilde{J}(t)$ or $\tilde{J}(t/\tau)$ (as in Fig.7(d)) serve as exemplary curves which determine the power-law dependence on $\varepsilon$ of the timescales $t_{0}$ and $t_{eq}$. These curves are $N$–independent for random–based initial conditions. Starting with small system sizes $N$ and by lowering the energy, one eventually encounters the crossover: the law describing these timescales changes from power to exponential and a dependence on $N$ emerges. Practically, this means that for $\varepsilon<\varepsilon_{c}$ the $\tilde{J}$–curves do not match with their exemplary counterparts, as in the example of Fig.8.

Two sets of three sigmoid curves are shown in Fig.8(a) (experiments with random initial positions and momenta) for decreasing energies with constant logarithmic step, namely $\log\varepsilon=-2.5$, $-3$, $-3.5$ and $N=8192$ (black set) or $N=128$ (green set). Only the initial parts of the $\tilde{J}(t)$–curves are shown, up to the value $0.1$. At $\log\varepsilon=-2.5$ the curves for the two values of $N$ nearly coincide, while a mismatch appears at $\log\varepsilon_{c}=-2.75$, clearly becoming larger by further lowering the energy. Fig.8(b) explains this mismatching in terms of the stability time $t_{0}$, defined as the times when the curve $\tilde{J}(t)$ first crosses the value $0.025$. Inspecting Fig.8(b), the delineation of $t_{0}$ to the $N$–independent final power-law (dashed line, obtained by fitting the data for $N=8192$) occurs at smaller crossover specific energies $\varepsilon_{c}$ as $N$ increases. Repeating the study for the initial excitation of a $12.5\%$ packet of modes with random phases (Fig.8(c) and (d)) find the same trend of $\varepsilon_{c}$ with $N$, with slightly different exponents. Notice that, for energies smaller than the crossover, one cannot always safely distinguish between an exponential law or a power law with steeper exponent than the exponent of the $N$–independent profile. For example, in Fig.8(d) the $N=128$ case could be fairly well fitted by a steeper power-law, in accordance with the interpretation via the ‘six-wave resonant interactions’ in [39]. Leaving open such questions, we here emphasize the utility of observing $\tilde{J}(t)$ for answering them, as well as for characterizing the asymptotic regime, which indicates that $\varepsilon_{c}\rightarrow 0$ as $N\rightarrow\infty$: the latter fact suffices to exclude exponentially-long lasting deviations from the equilibrium state at the thermodynamic limit ( Fig.8(e)).

Finally, we examine how the scaling laws on equilibrium times are affected by particular choices of initial conditions. To this end we consider six different types of initial conditions for the FPU system with $\alpha=1/2$ and $\varepsilon=0.01$, exciting: (1) the first normal mode, (2) the $12.5\%$ lowermost frequency packet of modes with zero initial phases (as in [12, 13, 14]), (3) the $12.5\%$ lowermost–frequency packet of modes with random initial phases (as in [3]), (4) random initial momenta, (5) random initial positions and finally, (6) random initial positions/momenta. The corresponding sigmoid curves are displayed in Fig.9(a), with solid lines for $N=8192$ and dashed lines for $N=1024$ ($\varepsilon=0.01$ in all cases). We immediately note the $N$–dependence of the behavior of the sigmoid curves for (1) and (2). Also, in (4) the spectrum $E_{k}$ exhibits equipartition already at $t=0$, hence $J(0)=J_{eq}$. On the other hand, in (5) and (6) $J$ starts below $J_{eq}$, since the energy spectrum somewhat favors high modes (see Eq.(6)). Finally, all curves equilibrate at $J_{eq}\approx 0.02015$, i.e. very close to $0.02$, as predicted by Eq.(6) for all $E_{k}$ equal.

Using now the numerical procedure described above, we calculate the equilibrium times $t_{eq}$ for five out of the six types of initial conditions in Fig.9(a) and for $N=8192$ (we exclude (4) in which $J(0)=J_{eq}$). As Fig.9(b) shows, in all cases we find $t_{eq}\sim\varepsilon^{-a}$, with the exponent $a$ varying between 2.5 (for the low–frequency mode excitations), to 3 (random positions and/or momenta). To within numerical uncertainties, one is tempted to conclude that yielding, in the initial data, more power to the high-frequency part of the mode spectrum results in steeper laws $t_{eq}\propto\varepsilon^{-a}$, i.e., trajectories more stable against the approach to equilibrium. We leave open the question of the associated scalings when only high modes are excited, a case not largely considered so far in literature.

4 Conclusion

We examined the role of the first Toda integral $J(q,p)$ as an indicator of the evolution, and crossing of various ‘stages of dynamics’ for FPU trajectories. This is accomplished by computing the time variations of the ‘normalized’ Toda integral $\tilde{J}$ (section 2) when an FPU trajectory $(q(t),p(t))$ is substituted within $J(q,p)$. Our main conclusions are summarized below.

Despite its apparent simplicity, $J(t)$ proves to be a very efficient indicator allowing to clearly distinguish stages of FPU dynamics for a wide class of initial conditions. In particular, it allows to clearly define two times, called the ‘time of stability’ $t_{0}$, and the ‘time to equilibrium’ $t_{eq}$. For times $t<t_{0}$ the FPU trajectories remain extremely close to the original integral surface $J(q_{0},p_{0})=const$. Thus, $t_{0}$ can be interpreted as a time of stickiness to the original Toda integral surface. On the other hand, $t_{eq}$ marks the approach of the system to equilibrium.

For classes of initial conditions based on random initial data (in position or momenta, or the phases for normal mode variables), as long as the initial energy spectrum $E_{k}$ is not in equipartition, the times $t_{0}$ and $t_{eq}$ are connected via a ‘sigmoid’ evolution of the curves $\tilde{J}(t)$. For large $N$, the graph of $\tilde{J}(t)$ becomes $N$–independent, so that we obtain one representative sigmoid curve $J(t/\tau)$ described by one time-scale $\tau\sim\varepsilon^{-a}$ with exponents $a>0$, depending on the class of initial conditions considered. However, by lowering $N$ we cross to an energy regime where the system becomes $N$–depended and its corresponding time-scales extend exponentially like $\tau\sim\exp\varepsilon^{-b(N)}$, $b(N)>0$. Such a phenomenon seems to disappear in the thermodynamic limit for two cases of random initial data with non-equipartitioned energy spectra.

We finally examine the information drawn from computing time fluctuations of $J$ in cases in which the system is initially at energy equipartition, with the initial conditions being distributed either far or close to the Gibbs measure. In such cases, through $J(t)$ we can clearly compute a speed of diffusion transversally to the Toda integral surfaces. This speed is always small, implying that an underlying nearly-integrable dynamics holds even when the initial conditions are close to equipartition, and even close to the final equilibrium state.

As a final comment, we expect similar features, as those encountered for the evolution of the first Toda integral, to appear also for the other Toda integrals, a subject currently under investigation. In fact, treating such integrals analogously as in Eq.(6) allows to conjecture that, at low energies, this behavior will be reflected satisfactorily by their harmonic approximations. A detailed investigation of this subject is proposed for future study.

5 Acknowledgments

We are indebted to G. Benettin and A. Ponno for very fruitful discussions which significantly improved the present work. H.C. was supported by the State Scholarship Foundation (IKY) operational Program: ‘Education and Lifelong Learning–Supporting Postdoctoral Researchers’ 2014-2020, and is co–financed by the European Union and Greek national funds.

Appendix A Analytic estimates

The quadratic terms of (3) derive from the terms $T_{1}=\frac{1}{2}a_{n}^{2}(p_{n}^{2}+p_{n}p_{n+1}+p_{n+1}^{2})$ and $T_{2}=\frac{a_{n}^{2}}{16\alpha^{2}}(a_{n+1}^{2}+a_{n}^{2}+a_{n-1}^{2})$. $T_{1}$ immediately yields the terms $\frac{1}{2}(p_{n}^{2}+p_{n}p_{n+1}+p_{n+1}^{2})$ and $T_{2}$ yields the terms $\frac{1}{8}[(q_{n+2}-q_{n})^{2}+4(q_{n+1}-q_{n})^{2}+(q_{n+1}-q_{n-1})^{2}]$ after Taylor–expanding any of the terms $a_{n}^{2}a_{n+1}^{2}=e^{\alpha(q_{n+2}-q_{n})}$ etc.

Now summing all quadratic terms together,

[TABLE]

their expression easily simplifies by matching terms like $\frac{1}{2}\sum_{n}[p_{n}^{2}+(q_{n+1}-q_{n})^{2}]$ which can be replaced by $\sum_{k}E_{k}\simeq E$. Furthermore, it is $\sum_{n}(q_{n+2}-q_{n})^{2}=\sum_{n}(q_{n+2}-q_{n+1}+q_{n+1}-q_{n})^{2}\approx 2\sum_{n}[\delta q_{n}\delta q_{n+1}+(q_{n+1}-q_{n})^{2}]$ and $\sum_{n}(q_{n+1}-q_{n-1})^{2}\approx 2\sum_{n}[\delta q_{n}\delta q_{n+1}+(q_{n+1}-q_{n})^{2}]$ which approximates the quadratic terms in (10) and gives the expression (4):

[TABLE]

.

Under the Fourier transform the above expression takes the form (5). In particular, it is $\sum_{n}p_{n}p_{n+1}=\sum_{k}(1-\frac{\omega_{k}^{2}}{2})P_{k}^{2}=\sum_{k}2\cos(\frac{\pi k}{N})P_{k}^{2}$ and $\sum_{n}\delta q_{n}\delta q_{n+1}\approx 2\sum_{k}\omega_{k}^{2}\cos(\frac{\pi k}{N})Q_{k}^{2}+2\sum_{l,m}c_{l,m}Q_{l}Q_{m}$, where $c_{l,m}=\sin(\frac{l\pi}{N})\sin(\frac{m\pi}{N})$ and the off–diagonal sum $\sum_{l,m}c_{l,m}Q_{l}Q_{m}$ is approximately zero. Therefore, we get that

[TABLE]