Understanding the FPU state in FPU-like models

Foreword

This paper is dedicated to our colleague and friend Antonio Giorgilli, in occasion of his 70$^{\mathrm{th}}$ birthday. The scientific contribution of Antonio, in different branches of Mathematical Physics, is well known and is invaluable. Concerning in particular FPU, beyond the interest of specific results, Antonio had a special propulsive role in Italy, and influenced very much our own work. "Understanding the FPU state...", which is the title of the present paper, started for us from studying ^[14] and the scaling law (3) there illustrated. Thanks, Antonio!

1.   Introduction

This paper is devoted to the Fermi-Pasta-Ulam (FPU) problem ^[1], namely the problem of energy sharing in weakly nonlinear chains of oscillators. The Hamiltonian has the form

$V$ being some nearest-neighbors potential with single minimum in zero,

In (1) the boundary conditions are still not specified; they are generally either fixed ends, i.e., $q_0 = q_{N+1} = 0$ (like in the original FPU paper), or periodic, $q_{N} = q_0$. The fixed-ends case can be viewed as a submanifold of the periodic case with $N' = 2(N+1)$ degrees of freedom, with symmetry $q_{N'-j} = -q_{j}$.

The value of $\alpha$, if different from zero, is irrelevant, since a trivial rescaling reports it to any prefixed value; the effective parameters determining the dynamics are indeed $|\alpha|\sqrt{\varepsilon }$, $\varepsilon = E/N$ being the specific energy, and then $\beta/\alpha^2$, … The choice of $\alpha$ actually fixes the energy scale, as well as the scale of $\beta$ and of possibly further coefficients in (2).

As is well known, Fermi, Pasta and Ulam aimed to investigate how the system reaches the statistical equilibrium, identified with the equipartition of energy among normal modes, if started very far from equilibrium, the whole energy being initially given to only one or two long-wavelength normal modes. With great surprise they found, within the time scale accessible to their computer, ^* no equilibrium at all: the system apparently reached a stationary state, different both from the initial state and from equipartition, in which only a few normal modes significantly share energy. Such a peculiar state is often called "the FPU state". Besides the lack of equipartition, the FPU state is characterized by long time recurrences, made especially evident a few years later in ^[3].

^*The celebrated MANIAC I of the Los Alamos Laboratories, capable to process approximately $10^4$ instructions per second ^[2]. Compare with modern computers, having hundreds or thousands of CPUs with clock at several GHz.

The literature on FPU, 65 years after the original FPU paper, is so huge, that it is not even conceivable to summarize it through a reasonable number of main references. Two collections of papers ^[4,5], produced in the occasion of the 50th anniversary of FPU, may give an idea of the status of the art around 2004–2005, as well as of the widespread attitudes of different authors, or different schools, in front of such non easy subject which touches the very foundations of classical statistical mechanics.

Two main questions have been extensively studied in the literature, since the very beginning:

ⅰ. Understanding the unexpected FPU state: its characteristic features, its long time recurrences, how is it produced from initial data, how is it possible that, in spite of the nonlinearity which couples normal modes, energy is shared among only a few of them.

ⅱ. Understanding whether, possibly on a much larger time scale, the FPU state further evolves, eventually reaching thermalization and normal statistical behavior, as in the original FPU expectation.

Both questions have been studied rather soon, independently. On the one hand, the recurrences observed in ^[1,3] suggested the presence of a close underlying integrable dynamics, which was initially searched, starting with ^[6] in 1965, among nonlinear wave equations, and later on, starting with ^[8] and ^[9], in the Toda model. On the other hand, almost simultaneously in 1966, the possibility was outlined ^[7] that normal statistical behavior is recovered if the specific energy $\varepsilon = E/N$ exceeds a threshold, which in turn, in ^[7], was conjectured to vanish for $N\to\infty$.^†

^†Ref. ^[7], as well as many later papers investigating the conditions for thermalization in FPU, deal with the so-called $\beta$ model, having $\alpha = 0$. The two models in fact (together with other variants of FPU, like Lennard-Jones) quite often have been considered equivalent and interchangeable, without much criticism. Certainly, they are not for what concerns the FPU state, which is peculiar of the $\alpha$-models. The long time behavior, with thermalization on a suitably long time scale, is instead rather similar; see for a comparison ^[10].

The literature on FPU accumulated, along these lines, without much order, and it is not so easy to follow it. Even the basic distinction between formation of the FPU state and possible later evolution to normal statistical behavior, with different phenomena occurring on different time scales, was not always clear at the beginning, and started being emphasized only in 1982 in ^[11,12], where a possible connection between the FPU state and metastability was outlined. Not so many papers are based on this distinction; among them, ^{[10,13,14,15,16,17,18]}.

In this paper we focus the attention on the former problem, namely understanding as much as possible the FPU state. The aim is to revisit the existing literature, tracing possible paths through it, and then to provide a further understanding of the FPU state by looking more closely at the Toda model.

Browsing the literature, we extract two main ideas (see Section 2 for references):

$–$ The FPU state is characterized by a few simple scaling laws; the most relevant one is that the number $M$ of normal modes effectively involved in the energy sharing, when one or a few long-wavelength normal modes are initially excited, scales as

$–$ The reference integrable dynamics underlying FPU, is not the linear one, made of free normal modes, rather the essentially nonlinear Toda dynamics.

Our contribution, in this background, is the following:

$–$ We check in Toda the scaling laws known to hold for FPU, and find they are fully satisfied. Scaling laws appear thus to be proper of Toda, and FPU obeys them as long as its dynamics remains close to the Toda dynamics.

$–$ We exploit the numerical procedure explained in ^[9] to compute numerically the Toda actions; the procedure allows to determine, in any given state $(p, q)$, the number of nonvanishing actions, that is the effective dimensionality of the corresponding Toda torus. We consider states in which only one or a few low frequency modes are excited, and observe that the dimensionality ${\cal M}$ of the torus scales exactly as in (3),

This means: the scaling laws characterizing the FPU state are contained in the nontrivial correspondence between normal modes and Toda actions. This is a more innovative view, in which dynamics (integration of Hamilton equations) does not play a role. By the way, the accuracy of the scaling law (4) is the best we ever met (see in particular Figures 8 and 9 below).

The paper is organized as follows: Section 2, after a few preliminaries, proposes our personal way to look at the literature on the FPU state, and also includes (Sect. 2.4) the new results concerning the scaling laws in the Toda model. Section 3 is instead devoted to the computation of the Toda actions. An Appendix follows, in which the procedure to compute the Toda actions is shortly recalled.

2.   The FPU state: scaling laws from dynamics

2.1.   A first look at the FPU state

Consider, to be definite, the case of fixed ends. Normal modes are then

their energies $E_k$ and frequencies $\omega _k$ are

For given initial data let^‡

^‡Restricting the average to the running window $\frac23 T\le t \le T$ does not change the asymptotics and makes a little quicker the convergence: time averages from the beginning are somehow lazy.

A preliminary qualitative view of the formation of the FPU state, and of its later slow evolution to equipartition, is provided by Figure 1 (a remake of Figure 1 in ^[19]) which refers to $N = 1023$, $\alpha = -1$, $\beta = 2$, specific energy $\varepsilon = 10^{-4}$, and shows $\overline E_k(T)$ vs. $k/n$ at selected times $T$ (marked in the figure) in geometric progression. The energy was initially equidistributed among the lowest $10\%$ of modes (the rectangular profile in the figure). One can see that quite soon, already at $T\simeq10^3$, a well defined profile is formed, in which only some low frequency modes effectively take part to the energy sharing, the energies of the remaining ones decaying exponentially with $k/N$. The energy profile keeps its form nearly unchanged for a rather large time scale, $T\simeq10^5$ or $10^6$, definitely much larger than the time needed to form it. Afterwords, on a much larger time scale, the dynamics slowly evolves towards energy equipartition, the high-frequency modes being progressively involved into the energy-sharing game; in the above conditions, equipartition requires $T\simeq10^{10}$. The initial phases (this will have importance, see later) were chosen randomly; moreover, to improve a little the figure reducing the fluctuations, an average of $E_k(T)$ on 24 different random choices of the phases was performed.

If the initial energy distribution is narrower, then the formation of the FPU state is slower, but the state looks the same, as well at the later evolution; see Figure 2, left panel, where the initial distribution is $0 < k/N < 0.03$. If energy is raised, the width of the FPU state grows, see the right panel of Figure 2 where $\varepsilon = 4\times10^{-4}$.

We shall now try to order the literature on the FPU state, by tracing two possible paths through it.

2.2.   Path one: scaling laws in FPU

The FPU state needs to be characterized quantitatively, looking in particular at its width——how many modes are effectively involved——and at the time scale for its formation, to be thought of as depending on $N$ and $\varepsilon$ as well as on the kind of the initial excitation.

Let us start by numerical evidence. Both in ^[20] and ^[14] a convincing evidence is produced (particularly accurate in ^[14]) that the width of the FPU state, namely $f = M/N$ where $M$ is the (suitably defined) number of modes effectively sharing energy, for large $N$ depends only on the specific energy $\varepsilon$ and not on $N$, and scales as

correspondingly, $M$ obeys (3). One should remark that ^[20] deals with a periodic pure $\alpha$ model, and considers an initial excitation flat on a fraction $f_0$ of modes already proportional to $\varepsilon ^{1/4}$, while ^[14] deals with a fixed ends model with $\beta = 4\alpha^2$, and only the first mode is initially excited. The notion of "number of modes effectively shearing energy", that is the definition of $M$, is also different, see the original papers; it is remarkable that the scaling law (5) is the same——the phenomenon looks, so to speak, robust. Some numerical evidence of (5) was already present in ^[21]. Concerning the time of formation of the FPU state, in ^[20] we find the estimate

holding for large $N$ and initial data as before, that is for $f_0$ obeying (5).^§

^§No explicit estimate of $t$ as function of $\varepsilon$ is provided in ^[14], but an indication is present, if we correctly interpret, $t\sim\varepsilon ^{-1/2}$. We conjecture the authors have seen the correct law $t\sim\varepsilon ^{-3/8}$ valid for fixed $M_0$, see below, which is close to $\varepsilon ^{-1/2}$ and easily confused with it in a non systematic study.

As is remarkable, the scaling laws (5) and (6) have been also justified theoretically, looking at the way normal modes interact. More precisely, (5) is deduced in ^[21] by adapting to the $\alpha$ FPU model the Chirikov method of overlapping of resonances, originally formulated for the $\beta$ model ^[7]; both (5) and (6) are instead deduced in ^[13], on which ^[20] is based, by adapting to FPU some ideas from the theory of weak turbulence. In both cases, the reference uncoupled model is the harmonic chain.

The scaling laws entering the formation of the FPU state have been later studied in ^[22] in a more systematic way. The conventional definition of $M$——actually the same as in ^[20]——is as follows: let

denote the "entropy" of the spectrum ^[23]; then^¶

^¶Let $\overline E_k = E/M$ for $k\le M$, $\overline E_k = 0$ for $k > M$. Then it is $\eta = \log M$.

Consider initial data with the first lowest $M_0$ modes initially excited; think $f$ to be a function of $t$, $N$, $\varepsilon$ and $f_0 = M_0/N$. The results in ^[22] can be summarized in two general scaling laws for $f(t, N, \varepsilon , f_0)$, both holding for large $N$, depending on whether the phases of the initially excited modes are chosen randomly or are coherent.

(a) For random initial phases, $f$ is independent of $N$ and scales as

where ${\cal G}$ is a kind of universal function of sigmoidal shape, reaching a plateau when its argument grows to infinity. According to (7), the width of the FPU state obeys (5), no matter which is $f_0$. The time scale for the formation of the state scales instead as

The main interest, having in mind thermodynamics, is for fixed $f_0$, that is $M_0$ proportional to $N$ (modes excited up to a given frequency). The scaling law (8) however also holds for $f_0\sim\varepsilon ^{1/4}$, as in ^[20], giving in this case (6), as well as down to small fixed $M_0$, as in ^[14], giving in such case times growing with $N$,

for fixed small $M_0$ the assumption that the initial phases are random is meaningless, and is in fact irrelevant to deduce (9).

(b) If instead many modes are initially excited with coherent initial phases, for example equal phases, or phases alternating in some easy pattern, then the total energy $E$, rather than the specific energy $\varepsilon$, appears to be the the relevant variable, and (7) is replaced by

For large $E$, in particular for fixed $\varepsilon$ and large $N$, the intermediate FPU state practically confuses with equipartition among all modes: unless of course $f_0\sim 1/N$, in which case, with interesting continuity, (9) is recovered.

It is worthwhile to remark that the choice of equal initial phases, although not well motivated, is rather common in the literature: it occurs indeed whenever, for example, either the $q_i$'s or the $p_i$'s initially vanish, and correspondingly energy is purely kinetic or purely potential (the easiest initial data to be handled numerically). With such initial data, unless $f_0\sim N^{-1}$, the FPU state for large $N$ disappears. The lack of attention to initial phases, in our opinion, is among the reasons why it not so easy, sometimes, to compare papers.

In ^[22] the general scaling laws (7) and (10) are also motivated theoretically, but in a way substantially different from ^[21] and ^[13], namely using, as a normal form for FPU, the KdV nonlinear wave equation. Considering FPU as a perturbation of a nonlinear integrable system is indeed an old idea, that we shall shortly revisit in the next subsection.

2.3.   Path two: FPU state and integrability

As already remarked in the Introduction, to explain the presence of long time recurrences in FPU it was soon suggested that the FPU dynamics is close to the dynamics of a nonlinear but nevertheless integrable system, and this was first searched among nonlinear wave equations. The first attempt in this direction, going back to 1965, is ^[6], a fundamental paper which is at the basis of modern theory of nonlinear integrable wave equations, in which FPU appears as a main motivation to study again, after years, the KdV equation,

The progress in the field is then rapid: a couple of years later the method of inverse scattering ^[24] and Lax pairs ^[25] are introduced, and the presence of infinitely many constants of motion is established ^[26]. Finally, in 1971, KdV is shown to be an infinite dimensional completely integrable Hamiltonian system ^[27].

In fact, the nonlinear wave equation which is more immediately related to FPU, if one searches for a continuum limit in which the first nonlinear and the first dispersive terms beyond the wave equation are kept, is the Boussinesq equation; a possible form is

From Boussinesq it is possible to deduce, in a convenient nontrivial limit, KdV, but the Boussinesq equation itself was soon proved, in 1973, to be integrable ^[28]. In ^[28] the connection between FPU and Boussinesq is quite emphasized, the lack of integrability of FPU being explained as due to the discretization.

In the meanwhile, statistical physicists become interested in the Toda model ^{[29,30,31,32]}. As is well known, this is a Hamiltonian system with the same form as (1), $V$ being the Toda exponential potential

In 1974 the Toda model was proved to be completely integrable, remarkably in three independent papers ^[8,33,34]. Reference ^[8] is particularly important for FPU, because the connection with FPU is there emphasized. Indeed, for $\lambda = 2\alpha$ it is

with

so the model is tangent to FPU up to the third order included and provides an integrable approximation to FPU better than the harmonic chain. In ^[8] the slow stochastization of FPU is no more attributed to the discretization with respect to an integrable wave equation, but to the small difference, with dominating term $\frac14(\beta-\beta_T)r^4$ in the potential, between FPU and Toda. The reference to Toda as the best integrable approximation for FPU, is a considerable change of paradigm in the FPU problem.

The connection between FPU and Toda was proposed again, and emphasized, in 1982, in ^[9]. In this paper, on the one hand a striking evidence is provided that the dynamics of FPU, in the FPU state before the slow drift towards equipartition, is hardly distinguishable from the Toda dynamics; on the other hand, the integrability of Toda is studied very constructively, and an algorithm is proposed to compute numerically the Toda actions. We shall come back to this point in Section 3.

2.4.   Scaling laws in the Toda model

The obvious conjecture, at this point, is that the scaling laws (7) and (10), characterizing the FPU state and implying in particular (5) and (6), are not proper of FPU but come from the underlying Toda integrable dynamics. This is what we shall see in this section, where we revisit and remake for Toda some numerical work originally done in ^[22] for FPU.

The scaling laws (7) and (10), reducing the variables from four to two, are obtained in ^[22] as combinations of three independent homogeneity relations characterizing the function $f(t, N, \varepsilon , f_0)$ for large $N$. Two of them do not depend on the choice of the phases of the excited modes and are

We cannot enter here the way (11) and (12), which might appear obscure, are (heuristically) deduced from KdV, demanding for this to ^[22].

Figure 3 provides a test of the scaling law (11). The left panel shows $f$ vs. $t$ for eight different choices of $N$, $\varepsilon$ and $f_0$ obeying (11), namely $N = 2896-9742$ with step $2^{-1/4}$, $\varepsilon = 10^{-3}-7.8\times10^{-6}$ with step $1/2$, and $f_0 = M_0/N$ with constant $M_0 = 102$. The right panel shows the same curves rescaled according to (11), namely rescaled horizontally by a factor $(\varepsilon /\varepsilon _*)^{-3/4}$ and vertically by a factor $(\varepsilon /\varepsilon _*)^{1/4}$, with $\varepsilon _* = 2.5\times10^{-4}$ (the value of $\varepsilon$ of the third curve from above, which does not move). One can see the rescaled curves superimpose almost exactly. The phases were chosen randomly, and each curve is in fact the average over 24 different random choices of the phases.

Figure 4 provides similarly a test of the scaling law (12). Here $\varepsilon = 2.4\times10^{-4}$ is kept fixed, while $N$ ranges between $1023$ and $16383$; $f_0 = M_0/N$ with fixed $M_0 = 26$. The left panel shows the curves as they are, while in the right panel $t$ is scaled by a factor $(N/N_*)^{-3/2}$, with $N_* = 4095$ (corresponding to the central curve, which stays unchanged). As before, each curve is the average of 24 different random choices of the initial phases.

Putting together (11) and (12) one gets

with a convenient ${\cal F}$ depending only on two variables. Such a relation can be further tested by considering combinations of $N$, $\varepsilon$ and $f_0$ different from the previous ones, for example one can keep $N$ constant and move $\varepsilon$; the corresponding scaling factors are in this case $\varepsilon ^{-3/8}$ on the horizontal axis, $\varepsilon ^{1/4}$ on the vertical axis. Figure 5 shows the result for $N = 1023$, seven values of $\varepsilon$ in the range $10^{-3}-1.56\times 10^{-5}$, $f_0 = M_0/N$ with small $M_0 = 3$. As before the left panel shows $f$ vs. $t$ as it is, while the right panel reports the rescaled curves.

The third relation which completes (11) and (12) depends instead in an essential way on the choice of the phases. For random phases, it simply tells that, for large $N$ and given $f_0$, $f$ is independent of $N$. This means ${\cal F}(x, y) = {\cal G}(xy^{3/2})$, and (7) is deduced. A test of the independence of $f$ on $N$ for large $N$ at fixed $f_0$ is provided in Figure 6, where $\varepsilon$ is either $10^{-3}$ or $2.5\times10^{-4}$, while $N$ assumes five different values in the range $1023-16383$; $f_0 = 0.025$ is fixed. The left panel displays $f$ vs. $t$; quite clearly, $f$ depends on $\varepsilon$ but not on $N$. The further rescaling of the vertical axis by a factor $\varepsilon ^{1/4}$, shown in the right panel, makes all curves superimpose.

For coherent phases instead the additional scaling law is

that is, $f$ is independent of $N$ at fixed $f_0$ only if the total energy $E = N\varepsilon$ is kept fixed. From (14) and (9), (10) is deduced. A test of (14) is exhibited in Figure 7, left panel, where $E = 0.265$ and $f_0 = 0.1$ are kept fixed, while $N$ assumes six values in the range $511-16383$. The right panel shows that if $\varepsilon$, instead of $E$, is fixed, then the independence of $N$ is lost and for large $N$ the FPU state does not exist any more. The phases are here constant, namely $\phi_k = \pi/2$ for all excited modes (purely potential initial energy), but a qualitative similar picture is obtained by choosing for example $\phi_k = 0$, or $\phi_k = k\pi/2$, or similar regular patterns.

3.   More inside Toda: scaling laws from Toda actions

In this section we shall investigate the basic scaling law (5) following a substantially different road, namely looking at the correspondence between linear normal modes and action variables in the Toda model. The aim is understanding the scaling law (5) governing the FPU state in Toda by looking only at the initial state——how many actions are excited, which is the effective dimensionality of the torus supporting the dynamics——with no need to integrate the equations of motion.

Let us remark that computing the Toda actions for large $N$ is a rather nontrivial work, facing delicate questions of numerical precision. We present here only a few preliminary results, without entering in details, with the only aim to show this is a promising way (see also ^[35]). We shall come back soon to the subject, in a wider devoted paper.

The appropriate definition of Toda actions, as well as the numerical procedure to compute them, are provided in ^[9]; for the convenience of the unfamiliar reader, however, the essentials are reported in the Appendix. The natural model to look at, in this connection, is a periodic chain with $N$ particles, that is, eliminating the barycenter, with $N-1$ degrees of freedom. Let us recall that for the corresponding harmonic periodic chain, normal modes are not stationary waves but traveling waves; a convenient choice is

where

(similarly for $p_i(t)$). The harmonic energies and harmonic actions are respectively

Let ${\cal I}_1, \ldots, {\cal I}_{N-1}$ denote the Toda actions, as defined in ^[9]. They turn out to be perturbations of the harmonic actions $I_1, \ldots, I_n$, and in particular, if a single harmonic action $I_k$ is different from zero, then in the limit of vanishing energy $E$, it is

The limit however, see ^[36,37,38], is highly non uniform in $N$, and for small $k$, to get convergence, $EN^3$, equivalently $\varepsilon N^4$, needs to be small. The perturbative approach instead fails, and the correspondence between normal modes and actions gets out of theoretical control, for large $N$ at fixed $\varepsilon$, that is in the thermodynamic limit.

Having in mind the FPU state and the scaling law (5) it satisfies, we naturally expect that for large $N$, if one or a small number $M_0$ of waves of low $k$ are present, then the number of Toda actions that are substantially different from zero scales as $N\varepsilon ^{1/4}$. This is precisely what we found in our numerical computations.

To define and compute the effective number ${\cal M}$ of non vanishing actions in any given state of the system, we found convenient to resort to the analog of spectral entropy, namely

then we let, in complete analogy, ${\cal M} = e^{ {\cal H}}$.

Figure 8 reports ${\cal M}$ vs. $\varepsilon$ in the wide range $10^{-15}\le\varepsilon \le1$, for several values of $N$ between $32$ and $8192$; the state is a single traveling wave with $k = 1$. The solid lines drawn in the figure are obtained by a least square fit, and all slopes, denoted $a$ and reported in the figure aside $N$, turn out to be very close to $1/4$. So, the behavior is as expected: in agreement with (15), all curves converge to ${\cal M} = 1$ for $\varepsilon \to0$, while for larger $\varepsilon$, with crossover depending on $N$, ${\cal M}$ grows proportionally to $\varepsilon ^{1/4}$.

Figure 9, left panel, is the same as Figure 8, but the fraction $f = {\cal M}/(N-1)$ of nonvanishing actions, instead of ${\cal M}$, is there reported vs. $\varepsilon$. Quite clearly all curves flatten, by growing $\varepsilon$, on the same line. The solid line drawn in the figure interpolates the data for $N = 8192$ and has slope $a = 0.25$. An overall check that the scaling law is

is provided by the right panel of Figure 9, where ${\cal M}$ is reported vs. $(N-1)\varepsilon ^{1/4}$. Here too, the line interpolates data for $N = 8192$; its computed slope is $0.998$. In conclusion, in the state of single traveling wave with $k = 1$, the expected scaling law is satisfied with quite great accuracy.

Such an accuracy was indeed a nice surprise, and was not really expected, the numerical procedure to compute the actions being new for us and somehow problematic (see the Appendix, Section C). It is hard to say how reliable are the values of the slopes reported in Figure 8. The only easy analysis is the uncertainty of the least-squares fit giving $a$. Figure 10 reports the computed slopes vs. $N$, the error bars corresponding, as is typical, to $\pm3$ std. deviations; the numerical results are clearly compatible with $a = 0.25$. Such an analysis of the accuracy does not take into consideration other possible sources of error: At the moment we are not able to draw error bars around the symbols of Figure 8, and the choice of the intervals of $\varepsilon$ where performing the best-fit analysis, after the crossover with the horizontal line $N = 1$, is also somehow arbitrary. Anyhow, the overall precision looks remarkable.

Nothing essentially changes if, in place of a single traveling wave with $k = 1$, one considers, for example, three waves with $k = 1, 2, 3$, having the same ${\cal E}_k$. Figure 11 shows, in the left panel, the analog of Figure 8, that is ${\cal M}$ vs. $\varepsilon$, and in the right panel the same as Figure 9, right panel, that is ${\cal M}$ vs. $N\varepsilon ^{1/4}$. Now the limit for $N\varepsilon ^{1/4}\to 0$ is of course ${\cal M} = 3$. The figure shows a curious phenomenon, that we did not try to investigate: before growing, here too as $\varepsilon ^{1/4}$, ${\cal M}$ decreases from $3$ to about $2$, and this for any $N$. The accuracy of the scaling law is not as good as for a single wave, but remains satisfactory.

Considering long waves with small $k$ is essential to produce a FPU state obeying the scaling law (16). If instead we consider a single short wave excitation with, for example, $k = 0.4N$ (recall the maximal frequency is at $k = N/2$), then the limit ${\cal M}\to 1$ for $\varepsilon \to0$ looks uniform in $N$, and up to rather large $\varepsilon \simeq0.1$, see Figure 12, only one action is appreciably different from zero. Note that the vertical scale is quite enlarged with respect to Figure 8: passing from ${\cal M} = 1$ to $2$ requires $\varepsilon \simeq1$.

Acknowledgements

We feel indebted to Martino Cecchetto (Padova) for helping us in understanding the computation of the Toda actions.

Conflict of interest

The authors declare no conflict of interest.

The Appendix: the procedure to compute the actions

We recall here the procedure to compute the actions proposed in ^[9]. Let us stress that the procedure is very clear in principle, but not so easy to implement in practice for large $N$. To understand, ref. ^[35] was very helpful to us.

A. The theoretical frame. Consider a periodic Toda chain with $N$ particles, assume the barycenter is at rest. Denote

so that

and let

$L^+$ is the Lax matrix associated to the system, used to prove integrability in ^[8,34], $L^-$ is an auxiliary matrix. Following ^[9], consider the characteristic polynomials

of $L^\pm$, which turn out to satisfy, for any $x$,

Though the matrices change during the evolution, the characteristic polynomials do not (this is indeed the essence of Lax method to prove integrability). In particular, the eigenvalues $\lambda^\pm_1, \ldots, \lambda^\pm_N$ of $L^\pm$ stay constant, and for either choice of the sign they provide an independent set of $N$ constants of motion. The eigenvalues will be supposed to be ordered decreasingly.

To define the actions it is convenient to introduce the discriminant $\Delta$, defined as

which then satisfies

We recall here some basic properties of $\Delta$, for details and discussion see ^[9]. To be definite we assume $N$ is even.

(a) At mechanical equilibrium ($a_i = 1$, $b_i = 0$), the shape of $\Delta$ is as in Figure 13, left panel, which refers to $N = 16$; $\Delta$ oscillates exactly between $\pm 2$ and all eigenvalues $\lambda^\pm_i$, but the extremal ones $\lambda^+_1$ and $\lambda^+_N$, are double. All eigenvalues stay between^|| $\pm2$.

^||Between $\pm1$ in ^[9]; for a certain convenience, our notation differs from ^[9] by a few sparse factors $2$.

$(b)$ For nonvanishing energy, one or more pairs of degenerate eigenvalues $\lambda^\pm_k, \lambda^\pm_{k+1}$ open, producing a "gap", and correspondingly there are maxima of $\Delta$ exceeding $2$ or minima exceeding $-2$. The shape of $\Delta$ looks as in the right panel of Figure 13, blue curve. Denote by $J_1, \ldots, J_{N-1}$, right to left, the intervals between nearby eigenvalues:

To each nonvanishing interval $J_k$ corresponds a positive action ${\cal I}_k$, defined by

Following ^[9], one can trace inside each open gap an arc $\rho'(x) = \pm (2+\rho(x))$, see the red curves of Figure 13. The actions are precisely, up to a trivial factor, the areas between $\rho'(x)$ and the lines $\pm2$. Notice that if $|\Delta|$ is large, then $\rho\simeq\log|\Delta|$.

B. An easy example. To exemplify the computation in an easy case, consider a state with a single sinusoidal wave with not small $k$. Figure 14 refers to $N = 32$, $k = 12\simeq0.4N$, $\varepsilon$ in the range $10^{-4}-10^{-1}$. Quite evidently, the action ${\cal I}_{12}$ is dominant, although it is clear that by increasing $\varepsilon$ other actions start to be appreciable. The peak of $\Delta$, already for $\varepsilon = 0.01$, is large, approximately $130$, and grows to approximately $4\times 10^6$ for $\varepsilon = 0.1$.

By increasing $N$ the peaks of $\Delta$ grow explosively, and this is the main problem in the numerical computation of the actions. Table 1 shows the peaks for some values of $N$ up to $1024$, $\varepsilon$ between $10^{-4}$ and $10^{-1}$. High peaks should not surprise: in a state with a single traveling wave it is $I_k = E_k/\omega _k = \varepsilon N/\omega _k$; should it be ${\cal I}_k\simeq I_k$, then (assuming the gap, i.e., the base of the lobe formed by $\rho'(x)$, is large) $\rho$ is of order $\varepsilon N$ and $\Delta$ is its exponential.^**

^**In ^[9], when peaks as large as $10^{15}$ are met, the comment is: The value of the discriminant at its maximum $(10^{15})$ is entirely typical, and perhaps surprisingly large for someone who has done inverse spectral theory with paper and pencil only.

C. Technical questions. Problems of theoretical precision, in presence of such large numbers, are serious. The results reported in Section 3 should be consider as highly tentative, subject to criticism and justified mainly by their overall coherence.

All computations have been made in the so-called quadruple precision, that is representing numbers with approximately $33$ decimal digits. This is fairly enough for small $N$, but questions arise for, say, $N = 1024$ or larger, unless $\varepsilon$ is very small.

To compute the actions one needs to compute with the best possible precision the eigenvalues $\lambda^\pm_i$. Since this is done only once per each state, computer time can be invested. The procedure we found convenient is:

$–$ Using standard routines, namely routines (in double precision, $16$ decimal digits) designed for sparse symmetric matrices, to compute the eigenvalues of $L^\pm$. This provides a first approximation of $\lambda^\pm_i$, with errors of the order, say, $10^{-14}$.

$–$ Improving the precision by searching the eigenvalues as zeros of $P^{\pm}(x)$, with the Newton method. Easy recursive expressions are available for the determinant of three-diagonal periodic matrices, that is for $P^{\pm}(x)$, as well as for the derivative. A little problem arises when the zeros of $P^{\pm}(x)$ are tangent or nearly tangent, but this is solved by approximating $P^{\pm}(x)$ with a parabola, and moreover, the corresponding gaps $J_k$ and actions ${\cal I}_k$ in this case are small, so one should not worry too much. The serious numerical problem instead is the usual loss of precision, when a zero results from a huge sum of large numbers. Here the only way out seems to us to be the use of multiprecision arithmetics, as we plan to do in a nearby future.

Once the $\lambda^{\pm}_i$ are obtained, the characteristic polynomials can be computed in the standard way,

$\Delta(x)$, $\rho(x)$ and the actions are immediately obtained. When $|P^\pm(x)|$ starts to be large, one can directly compute $\rho(x)\simeq\log|\Delta(x)|\simeq\log|P^\pm(x)|$,

Everything was somehow unusual in our experience, and a little painful.