Fermi–Pasta–Ulam–Tsingou problems: Passage from Boltzmann to -statistics
The Fermi-Pasta-Ulam (FPU) one-dimensional Hamiltonian includes a quartic term which guarantees ergodicity of the system in the thermodynamic limit. Consistently, the Boltzmann factor describes its equilibrium distribution of one-body energies, and its velocity distribution is Maxwellian, i.e., . We consider here a generalized system where the quartic coupling constant between sites decays as . Through first-principle molecular dynamics we demonstrate that, for large (above ), i.e., short-range interactions, Boltzmann statistics (based on the additive entropic functional ) is verified. However, for small values of (below ), i.e., long-range interactions, Boltzmann statistics dramatically fails and is replaced by q-statistics (based on the nonadditive entropic functional , with ). Indeed, the one-body energy distribution is q-exponential, with , and its velocity distribution is given by with . Moreover, within small error bars, we verify , which decreases from an extrapolated value q 5/3 to q=1 when increases from zero to , and remains q = 1 thereafter.
keywords:Fermi-Pasta-Ulam, Long-range interactions, q-statistics
Ludwig Boltzmann intensively tried to derive his celebrated weight, generalized by Gibbs into what is currently called the Boltzmann-Gibbs (BG) weight, from Newtonian mechanics and no other hypothesis. This is sometimes referred to as the Boltzmann program. The nonlinearity of the entangled particle dynamics makes the task a formidable one, and Boltzmann did not achieve it. Even today this remains as a very basic unsolved mathematical problem. This in no way means that we do not have a quite neat scenario about the validity of the Boltzmann weight. Indeed, it is clear by now that if the Newtonian dynamics of the system is such that its maximal Lyapunov exponent is positive (strong chaos), then the dynamics is mixing, hence ergodic, and it is on this basis (along the lines of the so-called Stosszahlansatz, molecular chaos hypothesis) that BG statistical mechanics is constructed. The situation is much more complex when the maximal Lyapunov exponent virtually vanishes (weak chaos). This is the discussion that we address here through a paradigmatic system, namely the Fermi-Pasta-Ulam (FPU) Hamiltonian FPU1955 (see details in ZabuskyKruskal1965 ; Gallavotti2008 ; OnoratoVozellaPromentLvov2015 ). This system plays a relevant role in the discussion of Fourier’s law for heat flow (heat conductivity), rapid (or slow) sharing of energy, eventually yielding equipartition of energy and the zeroth law of thermodynamics. It consists of a ring (chain with periodic boundary conditions) of oscillators which, in addition to their kinetic energies, interact through both harmonic and anharmonic terms. We focus on the following Hamiltonian ChristodoulidiTsallisBountis2014 ; BagchiTsallis2016 ; ChristodoulidiBountisTsallisDrossos2016 :
where and are the displacement and momentum of the -th particle with mass (from now on, without loss of generality, we can set , hence the momenta coincide with the velocities, i.e., ); , , and . Here , is the shortest distance between the -th and -th lattice sites (). If () we refer to short-range (long-range) interactions in the sense that the potential energy per particle is integrable (diverges) as ; in particular, the limit corresponds to first-neighbor quartic interactions (i.e., the historical FPU -model), and the value corresponds to a typical mean-field scenario. The Hamiltonian is made (formally) extensive for all values of by adopting the scaling factor AnteneodoTsallis1998 ; CirtoAssisTsallis2013 ; ChristodoulidiTsallisBountis2014 , which depends on . Note that for we have , which recovers the rescaling usually introduced in mean-field approaches, sometimes referred to as the Kac prescription factor. In the thermodynamic limit , remains constant, namely , for , whereas for , and for . It can be verified that the introduction of in Hamiltonian (1) is equivalent to a simple rescaling of time AnteneodoTsallis1998 .
We numerically integrate the Newton’s equations of motion using the symplectic velocity Verlet algorithm vel-Verlet with a small time-step such that the deviation of total energy in an isolated system is of the order of or less, until large times , depending on the system size . The initial conditions for the displacement variables are randomly chosen from a uniform distribution and the momenta from a normal distribution with unit variance, both centered around zero. Starting from a single random initial condition, the system is evolved for to allow the system to fully relax to its stationary (or quasi-stationary) state before starting averaging the steady state quantities for time-steps (here time is measured in units of ).
First we compute the distributions of the one-particle energy and of velocity as the long-range parameter is increased from zero on (). The energy and velocity distributions, and (for a homogeneous system and ), are obtained by collecting and of all the particles in the system and are shown in Fig. 1, a and b respectively.
For large (short-range interactions) the energy distribution does recover the expected Boltzmann distribution. However, for small values of (long-range interactions), the celebrated exponential distribution dramatically fails, and is replaced (within a fairly good numerical precision) by a -exponential one, where the index depends on ; see Figs. 2a and 2b where we show our simulation data fitted to the theoretical curves for two typical values, corresponding to the long-range and the short-range regimes respectively. Likewise, for the single particle velocity distribution, one obtains a Maxwell’s velocity distribution for large whereas for small the velocity histogram can be well approximated by a q-Gaussian function; this is shown in Figs. 2c and d respectively.
The adequacy of the -Gaussian and -exponential forms has been checked with quantities such as the -kurtosis and -ratio respectively. From the velocity histograms we can compute the -kurtosis of the distribution, defined as TsallisPlastinoEstrada2009 ; CirtoAssisTsallis2013 ; ChristodoulidiTsallisBountis2014
the value of is obtained by fitting the velocity histograms obtained from simulation. Using Eq. (2) the -kurtosis of any histogram can be computed. In particular, it can be verified that for any -Gaussian velocity distribution
Note that, for , we recover the well known kurtosis mandated by Gaussian distributions. Analogously, from the energy histograms we can compute the -ratio defined as follows
the value of is obtained by fitting the energy distributions. Using Eq. (4) the -ratio of any histogram can be computed. In particular, it can be verified that for any -exponential energy distribution
Note that, for , we recover the well known ratio mandated by exponential distributions. The numerical data obtained for the kurtosis and the ratio are displayed in Fig. 3,a and b respectively, along with their analytical expressions.
The results for the -indices (as well as for the associated inverse temperatures ’s), presented in Fig. 4a and b are very eloquent, and deserve some comments. The numerical values for the indices and coincide within small error bars, for all the values of that have been computationally run. This fact surely is nice and simple if we take into account the fact that both indices are associated with one-variable marginal distributions coming from the same many-body distribution (in a -dimensional phase space), which is of course numerically inaccessible for the large values of that have been used in the present calculations. Naturally, the precise -exponential and -Gaussian forms fail in the numerical regions of too high one-body energies and too high velocities respectively. The goal is in fact to attain as best as possible the limit, where the analytical forms could possibly be correct for all energies and all velocities. Such extrapolations have been done as shown in Fig. 4a: they exhibit for the limit, and a monotonically decreasing value for when increases up to , and for , with . The arguments based on the positivity of the largest Lyapunov exponent (see, for instance, TirnakliBorges ) rather suggest . However, all the related numerical results available up to now in various classical -dimensional models CirtoAssisTsallis2013 ; ChristodoulidiTsallisBountis2014 ; ChristodoulidiBountisTsallisDrossos2016 ; CirtoRodriguezNobreTsallis2017 systematically and intriguingly indicate . This robust numerical fact remains so even if , time, precision, integrating algorithms and other circumstances, are modified. This unexpected peculiarity has remained irreducible up to now, and it might suggest a distinction between strongly long-range-interacting systems () and weakly long-range-interacting systems (), with roughly in the range ; the strictly short-range-interacting systems would therefore correspond to . This situation is somehow reminiscent of say the quantum Ising ferromagnet with long-range interactions, which is known to present three (and not only two) thermostatistical regimes ThreeRegimes , namely , , and , corresponding therefore to . The physical interpretation of the present most interesting three regimes remains elusive. However, one possibility could be that the emergence of a neat Boltzmann regime (i.e. ) in a microcanonical ensemble demands not only ergodicity over the entire phase space (obviously assured by a positive largest Lyapunov exponent), or over a nonvanishing-Lebesgue-measure part of it, but, in addition to that, an uniform probability distribution over that region of phase space. If so, the thermostatistical regime for would of course ultimately satisfy ergodicity, but not equal probabilities effectively (at least not yet at the largest times that have been computationally attained). The deep understanding of this point would further enlighten the first-principle conditions of validity of the celebrated Boltzmann-Gibbs statistical mechanics.
Acknowledgments: We acknowledge fruitful discussions with L.J.L. Cirto, E.M.F. Curado, F.D. Nobre, A.R. Plastino, P. Rapcan, G. Ruiz, G. Sicuro, A.M.C. Souza, U. Tirnakli and R. Wedemann, as well as partial financial support from CNPq and Faperj (Brazilian agencies), and from the John Templeton Foundation (USA).
- (1) E. Fermi, J. Pasta and S. Ulam, Studies of nonlinear problems, Document Los Alamos-1940 (1955). Although not included as co-author, it was Mary Tsingou who prepared the code and performed the pioneering computational calculations on MANIAC. Following the excellent suggestion by T. Dauxois, Fermi, Pasta, Ulam, and a mysterious lady, Physics Today 6 (1), 55-57 (2008), whenever appropriate we shall from now on refer to the Fermi-Pasta-Ulam-Tsingou problem.
- (2) N.J. Zabusky and M.D. Kruskal, Interactions of solitons in a collisionless plasma and the recurrence of initial states, Phys. Rev. Lett. 15 (6), 240-243 (1965).
- (3) G. Gallavotti, The Fermi-Pasta-Ulam Problem: A Status Report, Lecture Notes in Physics 728 (Springer, 2008).
- (4) M. Onorato, L. Vozella, D. Proment and Y. Lvov, Route to thermalization in the -Fermi-Pasta-Ulam system, Proceedings of the National Academy of Sciences of the United States of America 112 (14), 4208-4213 (2015).
- (5) U. Tirnakli and E. P. Borges, Sci. Rep. 6, 23644 (2016).
- (6) H. Christodoulidi, C. Tsallis and T. Bountis, Fermi-Pasta-Ulam model with long-range interactions: Dynamics and thermostatistics, EPL 108, 40006 (2014).
- (7) D. Bagchi and C. Tsallis, Sensitivity to initial conditions of -dimensional long-range-interacting quartic Fermi-Pasta-Ulam model: Universal scaling, Phys. Rev. E 93, 062213 (2016).
- (8) H. Christodoulidi, T. Bountis, C. Tsallis and L. Drossos, Dynamics and Statistics of the Fermi–Pasta–Ulam –model with different ranges of particle interactions, JSTAT 123206 (2016) (13 pages).
- (9) C. Anteneodo and C. Tsallis, Breakdown of the exponential sensitivity to the initial conditions: Role of the range of the interaction, Phys. Rev. Lett. 80, 5313-5316 (1998).
- (10) L.J.L. Cirto, V.R.V. Assis and C. Tsallis, Influence of the interaction range on the thermostatistics of a classical many-body system, Physica A 393, 286-296 (2014).
- (11) M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1989).
- (12) C. Tsallis, A.R. Plastino and R.F. Alvarez-Estrada, Escort mean values and the characterization of power-law-decaying probability densities, J. Math. Phys. 50, 043303 (2009).
- (13) J.S. Andrade Jr., G.F.T. da Silva, A.A. Moreira, F.D. Nobre and E.M.F. Curado, Thermostatistics of overdamped motion of interacting particles, Phys. Rev. Lett. 105, 260601 (2010).
- (14) L.J.L. Cirto, A. Rodriguez, F.D. Nobre and C. Tsallis (2017), to be published.
- (15) P. Hauke and L. Tagliacozzo, Phys. Rev. Lett. 111, 207202 (2013).