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

1.   Introduction

The numerical experiment of Fermi, Pasta and Ulam ^[19,21] 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,31,32]. 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 = E/N$ (specific energy) is small and fixed while $N$ becomes large.

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

ⅰ) 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,22,23,24,29,37,38,39]}. 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 ^[22,23,35]. 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 ^{[12,14,22,23]}, and clearly beyond the harmonic regime (although of course still small, i.e., for specific energy $\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.

ⅱ) Stochasticity threshold: In 1959 Chirikov published his celebrated works ^[15,16,27] 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 ^[41], 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 ^[22,23] and a weaker decay of the form $1/N$ for the $q$-tori ^[14].

However, as noted in ^[41], 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,26,38]}. Already in 1982 the pioneering work of Ferguson et al. ^[20] 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 ^[33].

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 ^[26]. 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: ⅰ) single mode or packet of modes excitations with coherent or random phases, ⅱ) random data, where the whole energy spectrum has comparable power from the start, and ⅲ) generic data close to energy equipartition. Cases (ⅰ) and (ⅱ) are characterized by distinct phases of evolution of the FPU trajectories, which roughly correspond to the stages of dynamics recognized in ^[1,5,38]. 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 (ⅲ), 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 ^[21] considered the one-dimensional lattice consisting of $N$ weakly nonlinearly coupled oscillators. The dynamics of this model is described by the Hamiltonian:

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({n k\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 ^[42]

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 ^[25,30]. Besides $H_T$, the first non-trivial Toda integral after rescaling with appropriate constants, takes the form:

As described below, independently of the initial conditions considered, when the system comes close to equilibrium, $J(t)$ stabilizes to a final value, called hereafter, the 'equilibrium value' $J_{eq}$. Then $J(t)$ exhibits only rapid and irregular fluctuations around the value $J_{eq}$ for all $t > t_{eq}$. The time $t_{eq}$ is hereafter called equilibrium time. The normalized $J$ is defined as:

We outline the sigmoidal evolution of $\tilde{J}(t)$ from 0 to 1 along FPU trajectories through the example of Figure 1. Consider the excitation of the first $12.5\%$ of the normal modes for the FPU system with $\alpha = 1/2$, $\varepsilon = 0.01$ and $N = 8192$ particles. The detachment of $\tilde{J}$ from zero, which indicates an essential departure from quasi-integrable behavior, occurs at $t_0$. During the time interval $[0, t_0]$ $J(t)$ is nearly constant, while FPU behaves as Toda. 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}$.

Figure 1.  Evolution of the normalized Toda integral $J$ and of the tail energy $\eta$ when exciting the first $12.5\%$ of the normal modes in the system with $N = 8192$, $\alpha = 1/2$ and $\varepsilon = 0.01$. A 'nearly integrable' behavior (near constancy of $J(t)$) appears for $0\leq t\leq t_0$, a sigmoid increase of both indicators afterwards, leading to a new stabilization after $t\geq t_{eq}$. The times $t_0$ and $t_{eq}$ denote the 'stability time' and 'equilibrium time' respectively.

DownLoad: Full-Size Img PowerPoint

The Toda integral can be compared with another observable used in the literature ^[4,38] for distinguishing the stages of dynamics in the FPU problem, namely the tail energy^*:

^* A idea similar to the tail energy for generic packet excitations can be traced back to reference ^[36].

We find that the normalized Toda integral $\tilde{J}$ and the tail energy practically coincide for trajectories with initial low-mode excitations, as in Figure 1. However, as discussed below, $\tilde{J}$ is of use also in more generic initial conditions, in which the tail energy cannot be used.

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 Figure 2, and we theoretically justify below, that the stabilization of $J$ to $J_{eq}$ for $t > t_{eq}$ is related to the following: (ⅰ) the variables $p_n$, $\delta q_n$ decorrelate (Figure 2(b), (f)), (ⅱ) the energy spectrum $E_k$ reaches equipartition (Figure 2(d)).

Figure 2.  Excitation of the first normal mode. (a) Evolution of the normalized Toda integral $J$. (b) Evolution of the sums $\sum p_np_{n+1}$, $\sum \delta q_n\delta q_{n+1}$ as indicators of correlation between the real space variables. (c) and (d): The energy spectrum $E_k/E$ at $t = t_0$ and $t = t_{eq}$ respectively. (e) and (f): The momenta distributions at $t = t_0$ and $t = t_{eq}$ respectively. The black curve indicates a Gaussian distribution with dispersion $\sigma_p^2 = 2\varepsilon$.

DownLoad: Full-Size Img PowerPoint

Regarding (ⅰ), Figure 2(b) displays the evolution of the sums $\sum p_np_{n+1}$, $\sum \delta q_n\delta q_{n+1}$ for the same orbit as in Figure 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 ^[34]

as an accurate and simple observable to be used in specifying the equilibrium time. Other recent methods to estimate the 'ergodization time' can be found in ^[17], by measuring the times where an FPU trajectory intersects the 'equilibrium manifold', and in ^[18] by following the statistical properties of the orbit's excursions in appropriate action-angle coordinates.

We will now show, that property (ⅱ) above is connected with the correlation sums $\sum p_np_{n+1}$, $\sum \delta q_n\delta q_{n+1}$ both tending to zero for times beyond $t_{eq}$. To this end, the function $J$ in (2.3) can be well approximated by using Taylor expansions up to quadratic terms in terms of the quantities $p_n$, $\delta q_n$, when the specific energies are smaller than unity, and we get (see Appendix A):

Due to Eq.(2.6), the near-constancy of $J$ for $t < t_0$ in Figure 2(a) emerges from the counterbalance of the sums $\sum p_np_{n+1}$ and $\sum \delta q_n\delta q_{n+1}$ in Figure 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.(2.6) can be recast as a sum over the energy spectrum $E_k$:

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. Figure 2(e), (f)). This is not far from the Gibbs measure for an equilibrium state with Hamiltonian (2.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.(2.7) 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.

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: (ⅰ) with 'localized energy spectrum', (ⅱ) with 'delocalized energy spectrum', and (ⅲ) 'close to equipartition'.

(i) Localized energy spectrum: Case (ⅰ) 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_kA_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 is performed when distributing the total energy $E$ by setting $E_k(0) = E w$ among a percentage $1\leq k /N \leq w$ of modes with $0 < w\leq 1$, while setting $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 Figure 3(b). By progressively increasing the width of the packet (from $12.5 \%$ to $87.5 \%$), different localization patterns emerge in Figure 3(b). However, the evolution of $J$ in Figure 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.

Figure 3.  Initial excitation of the modes $0\leq k/N \leq w$ with random initial phases (as indicated in the figure) for $N = 8192$, $\alpha = 1/2$ and $\varepsilon = 0.01$: (a) evolution of $J$, (b) Time–averaged energy spectra during the metastable state (for $t < t_0$).

DownLoad: Full-Size Img PowerPoint

(ii) Delocalized energy spectrum: Case (ⅱ) 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$. Figure 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 Figure 4(b) for a single realization, $J(t)$ exhibits a similar sigmoidal evolution to the case (ⅰ). The approximations (2.6) and (2.7) 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 ^[39]). 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.

Figure 4.  FPU with $N = 8192$, $\alpha = 1/2$ and $\varepsilon = 0.01$. (a) Example of the harmonic energy spectrum $E_k$ for an initial condition corresponding to random (non-Gaussian) positions and momenta. (b) Evolution of $J$ (red) and its two quadratic approximations according to Eqs.(2.6) (green) and (2.7) (blue).

DownLoad: Full-Size Img PowerPoint

(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_kA_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

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 (ⅰ) and (ⅱ). Figure 5(a), (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 Figure 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

Figure 5.  Evolution of $J$ for two types of initial conditions in the FPU system with $N = 8192$, $\alpha = 1/2$ and $\varepsilon = 0.01$. (a) Energies starting at equipartition $E_k(0) = \varepsilon$ with randomly chosen initial phases for the normal mode variables, and (b) random positions with uniform distribution (see text). Both panels show the evolution of the mean $J$ from 10 realizations per class of initial conditions. According to the formulas (2.6) and (2.7) in both cases we have $J(0)\simeq J_{eq}$ for each individual trajectory.

DownLoad: Full-Size Img PowerPoint

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. Figure 6(a) 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$.

Figure 6.  The dispersion $\sigma^2_J(t)$ vs. time in logarithmic scale, measuring the diffusion spread transversally to the integral surface defined by the Toda integral $J$. (a) Results from 30 trajectories in the 'random initial momenta' case for $\varepsilon = 0.01$. (b) Same but for the 'close to equilibrium' case.

DownLoad: Full-Size Img PowerPoint

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 Figure 6(b). Now, the initial conditions are chosen by the distribution function of Eq.(3.1), 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 (ⅰ) and (ⅱ), 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, Figure 7(a) shows the sigmoid curves when considering the excitation of the first normal mode (case (ⅰ)), 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, Figure 7(b) shows the same computation but for initial conditions extracted by random positions and momenta (case (ⅱ)). Superposing the sigmoid curves for various $N$ shows a near-coincidence at all times.

Figure 7.  Upper panels: The normalized Toda integral $\tilde{J}$ along the FPU dynamics for $\alpha =1/2$ and $\varepsilon =0.01$ and various $N$ values for: (a) the excitation of the first normal mode, and (b) random initial positions and momenta. Lower panels: (c) The sigmoid evolution of $J$ for random positions and momenta at various energies. (d) Rescaling time by a constant $\tau \sim \varepsilon ^{-2.7}$, all sigmoid curves $J$ cluster in one: $J(t/\tau )$.

DownLoad: Full-Size Img PowerPoint

We note that a similar behavior to Figure 7(b) holds for initial data as those of Figure 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 Figure 7(c). With an appropriate time-shift $\tau \sim \varepsilon ^{-2.68}$ (inset), those curves fall one onto the other as shown in Figure 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

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 Figure 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 Figure 8 and as a result their stability times $t_0$ deviate from $t_0\sim \varepsilon ^{-a}$.

Figure 8.  (a) Random initial positions and momenta: $\tilde{J}$ for the exemplary $N = 8192$ curve (black) and its disagreement with $N = 128$ (green curves), which is evident below $\varepsilon_c = 10^{-2.75}$. (b) Stability times for the case (a): $N = 8192$ (power-law) and the gradual exponential divergence for $N = 128$, $64$ and $32$ particles. (c) and (d): same as (a) and (b) respectively, when the first $12.5 \%$ of modes are excited with random phases. (e) The specific energy crossover of the two cases versus $N$ decays to zero as $N^{-1.5}$ ((a), (b) case with red triangles) and as $N^{-2.5}$ ((c), (d) case with blue polygons).

DownLoad: Full-Size Img PowerPoint

Practically, we evaluate the stability times by computing the time at which the curve $\tilde{J}(t)$, which ranges from 0 to 1, crosses for the first time the value $0.025$, as Figure 8(a), (c) indicates with the red dashed line. Two sets of three sigmoid curves are shown in Figure 8(a) (experiments with random initial positions and momenta) for decreasing energies with constant logarithmic step, namely $\log_{10} \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_{10} \varepsilon = -2.5$ the curves for the two values of $N$ nearly coincide, while a mismatch appears at $\log_{10}\varepsilon_c = -2.75$, clearly becoming larger by further lowering the energy. Figure 8(b) explains this mismatching in terms of the stability time $t_0$, where 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 (Figure 8(c), (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 Figure 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 ^[40]. 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 (Figure 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 Figure 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.(2.7)). Finally, all curves equilibrate at $J_{eq} \approx 0.02015$, i.e., very close to $0.02$, as predicted by Eq.(2.7) for all $E_k$ equal.

Figure 9.  (a) Six types of initial conditions, lead to different sigmoid curves for $J(t)$. Continuous lines: $N = 8192$, dashed line: $N = 1024$. It turns out that the first normal mode and the packets with coherent phases are non-extensive. (b) Equilibrium times versus the specific energy, measured for the FPU system with $N = 8192$, $\alpha = 1/2$ and for 5 classes of initial conditions.

DownLoad: Full-Size Img PowerPoint

We calculate the equilibrium times $t_{eq}$ corresponding to all the cases of Figure 9(a) excluding (4) in which $J(0) = J_{eq}$. Practically, we evaluate $t_{eq}$ by computing the time at which the curve ${J}(t)$ reaches for the first time the value $J_{eq}$ (this is equivalent to $\tilde{J}(t)\approx 1$). As Figure 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.(2.7) 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.

A.   Appendix: Quadratic approximation of the Toda integral

The quadratic terms of (2.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^2a_{n+1}^2 = e^{ \alpha (q_{n+2}-q_n) }$ etc.

Now summing all quadratic terms together,

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 (A.1) and gives the expression (4):

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_lQ_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_lQ_m$ is approximately zero. Therefore, we get that

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.

Conflict of interest

The authors declare no conflict of interest.