text:
# Lyapunov exponent of many-particle systems: testing the stochastic approach

## I Introduction

In a recent paper  we presented a theory (“the stochastic approach”) that allows to estimate the largest Lyapunov exponent of many-particle systems having a smooth Hamiltonian. This theory is inspired by and complementary to the work of Pettini et al. , and Barnett et al. .

Starting from the linear equations that describe the evolution of tangent vectors, we used in  the techniques of stochastic linear differential equations to obtain an average propagator whose diagonalization provides the desired Lyapunov exponent. Even though the result obtained in this way is formally (almost) exact, approximations must be invoked in practical situations to calculate a concrete Lyapunov exponent.

The aim of this paper is to show that, even if one makes several crude approximations, the theory still keeps the strength to describe quantitatively some weakly chaotic many-particle systems. We do this by comparing the predictions of the theory with the outcomings of numerical simulations of the Hamiltonian many-particle dynamics.

The paper is organized in the following way. For reasons of self-containedness, we begin by presenting a short review of the theory in Sect. II. The systems to be studied (the mean-field XY-Hamiltonians) are introduced in Sect. III, where we also work out the predictions of the theory for these particular systems. Section IV contains the critical comparison of theoretical results with numerical simulations. We close in Sect. V with some remarks.

## II Review of the theory

The theory developed in  can, in principle, be applied to any dynamical system. However, given that the approach is perturbative, its success will in general depend on the choice of the unperturbed system. Here we restrict ourselves to perturbations of ballistic motion, i.e., Hamiltonians of the type

$$=\frac{1}{2I}\underset{i=1}{\overset{N}{}}p_i^2+𝒱(q_1,\mathrm{},q_N),$$
(1)
with $`q_i`$ and $`p_i`$ conjugate position-momentum coordinates, $`I`$ the mass (or moment of inertia) of the particles, and the interaction $`𝒱`$ that we assume to be small.

Let us denote phase space points by $`x=(q_1,\mathrm{},q_N,p_1,\mathrm{},p_N)`$ and tangent vectors by $`\xi =(\delta q_1,\mathrm{},\delta q_N,\delta p_1,\mathrm{},\delta p_N)^T`$. Differentiating the Hamilton equations, one obtains the evolution equations for tangent vectors:

$$\dot{\xi }=𝐀(t)\xi .$$
(2)
For a Hamiltonian of the special form (1), the operator $`𝐀`$ has the simple structure

$$𝐀(t)=\left(\begin{array}{cc}0& 𝟙/𝕀\\ 𝐕(t)& 0\end{array}\right).$$
(3)
Here $`𝐕`$ is the Hessian matrix of the potential $`𝒱`$, namely

$$V_{ij}=\frac{^2𝒱}{q_iq_j}.$$
(4)
Once initial conditions $`x_0`$ and $`\xi _0`$ have been specified, the Lyapunov exponent $`\lambda `$ is found by calculating the limit

$$\lambda =\underset{t\mathrm{}}{lim}\frac{1}{2t}\mathrm{ln}|\xi (t;x_0,\xi _0)|^2.$$
(5)
We assume that for any initial condition $`x_0`$ in phase space, the trajectory $`x(t;x_0)`$ is ergodic on its energy shell. This implies that $`\lambda `$ depends only on energy and other system parameters, but not on $`x_0`$, which can then be chosen randomly according to the microcanonical distribution. There will also be no dependence on initial tangent vectors, because if $`\xi _0`$ is also chosen randomly, it will have a non-zero component along the most expanding direction. Moreover, if the corrections to the exponential law implicit in Eq. (5) are neglected, one can write

$$\lambda \underset{t\mathrm{}}{lim}\frac{1}{2t}\mathrm{ln}|\xi (t;x_0,\xi _0)|^2,$$
(6)
brackets meaning microcanonical averages over $`x_0`$. This approximation will be analyzed in Sect. IV with the help of numerical simulations.

By letting $`x_0`$ be a random variable, $`𝐕(t;x_0)`$ becomes a stochastic process, and Eq. (2) can be treated as a stochastic differential equation. However, as we are interested in the square of the norm of $`\xi `$, we focus, not on Eq. (2) itself, but in the related equation for the evolution of the “density matrix” $`\xi \xi ^T`$:

$$\frac{\mathrm{d}}{\mathrm{d}t}(\xi \xi ^T)=𝐀\xi \xi ^T+\xi \xi ^T𝐀^T\widehat{𝐀}\xi \xi ^T,$$
(7)
the rightmost identity defining the linear superoperator $`\widehat{𝐀}`$. For the purpose of the perturbative approximations to be done, the operator $`\widehat{𝐀}`$ is split into two parts

$$\widehat{𝐀}=\widehat{𝐀}_0+\widehat{𝐀}_1(t),$$
(8)
where $`\widehat{𝐀}_0`$ corresponds to the evolution in the absence of interactions. In our case $`\widehat{𝐀}_0`$ and $`\widehat{𝐀}_1`$ are associated with

$$𝐀_0=\left(\begin{array}{cc}0& \\ 0& 0\end{array}\right)\text{and}𝐀_1=\left(\begin{array}{cc}0& 0\\ 𝐕(t)& 0\end{array}\right),$$
(9)
respectively (we have set $`I=1`$, but it can be brought back at any moment by dimensional considerations). Whenever $`𝐀_1(t)`$ is small, it is possible to manipulate Eq. (7) to derive an explicit expression for the evolution of the average of $`\xi \xi ^T`$:

$$\xi \xi ^T(t)=e^{t\widehat{𝚲}}\xi _0\xi _0^T,$$
(10)
where $`\widehat{𝚲}`$ is a time-independent superoperator given by the perturbative expansion:

$`\widehat{𝚲}=\widehat{𝐀}_0+\widehat{𝐀}_1`$ (11)
$`+{\displaystyle _0^{\mathrm{}}}d\tau \delta \widehat{𝐀}_1(t)e^{\tau \widehat{𝐀}_0}\delta \widehat{𝐀}_1(t\tau )e^{\tau \widehat{𝐀}_0}+\mathrm{},`$ (12)
with

$$\delta \widehat{𝐀}_1(t)=\widehat{𝐀}_1(t)\widehat{𝐀}_1.$$
(13)
A clear exposition of this derivation, together with a detailed discussion of its domain of validity has been given by van Kampen .

Let $`L_{\mathrm{max}}`$ be the eigenvalue of $`\widehat{𝚲}`$ which has the largest real part. Taking the trace of Eq. (10), one sees that the largest Lyapunov exponent $`\lambda `$ is related to the real part of $`L_{\mathrm{max}}`$ by

$$\lambda =\frac{1}{2}\text{Re}\left(L_{\mathrm{max}}\right).$$
(14)
In Eq. (11) we give explicitly only the first two cumulants, the dots stand for third cumulants and higher order ones. If all cumulants were summed up, Eq. (10) would be exact in the long-time regime $`t\tau _c`$ . From now on, we restrict our analysis to the propagator $`\widehat{𝚲}`$ truncated at the second order. The perturbative parameter controlling the quality of the truncation is the product of two quantities, the “Kubo number” $`\sigma \tau _c`$. The first factor, $`\sigma `$, characterizes the amplitude of the fluctuations of $`\widehat{𝐀}_1(t)`$. The second, $`\tau _c`$, is the correlation time of $`\widehat{𝐀}_1(t)`$.

To proceed further one needs the matrix of $`\widehat{𝚲}`$ in some basis. The crudest approximation consists in restricting $`\widehat{𝚲}`$ to the subspace spanned by the following three matrices:

$$𝐈_1=\left(\begin{array}{cc}& 0\\ 0& 0\end{array}\right),𝐈_2=\left(\begin{array}{cc}0& 0\\ 0& \end{array}\right),𝐈_3=\left(\begin{array}{cc}0& \\ & 0\end{array}\right).$$
(15)
This choice is equivalent to a mean-field approximation in tangent space . After some lengthy algebra one arrives at the corresponding 3$`\times `$3 matrix for $`\widehat{𝚲}`$:

$$𝚲=\left(\begin{array}{ccc}0& 0& 2\\ 2\sigma ^2\tau _c^{(1)}& 2\sigma ^2\tau _c^{(3)}& 2\mu \\ \mu +2\sigma ^2\tau _c^{(2)}& 1& 2\sigma ^2\tau _c^{(3)}\end{array}\right),$$
(16)
with the definitions

$`\mu `$ $`=`$ $`{\displaystyle \frac{1}{N}}\mathrm{Tr}𝐕,`$ (17)
$`\sigma ^2`$ $`=`$ $`{\displaystyle \frac{1}{N}}\mathrm{Tr}\left(\delta 𝐕\right)^2,`$ (18)
$`\tau _c^{(k+1)}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑\tau \tau ^kf(\tau ),`$ (19)
where we have introduced the normalized correlation function $`f(\tau )`$

$$f(\tau )=\frac{1}{N\sigma ^2}\mathrm{Tr}\delta 𝐕(0)\delta 𝐕(\tau ).$$
(20)
\[A Lyapunov exponent with the correct units (inverse time) is obtained by making the substitutions $`\mu ,\sigma \mu /I,\sigma /I`$ in Eq. (16).\]

In the mean-field approximation the Lyapunov exponent is expressed in terms of the set of four parameters $`\mu `$ and $`\sigma ^2\tau _c^{(k+1)}`$, $`k=0,1,2`$. The parameters $`\mu `$ and $`\sigma `$ are, respectively, the mean and variance of the stochastic process $`𝐕(t)`$, and can in principle be obtained analytically by calculating the corresponding microcanonical averages.

The characteristic time $`\tau _c^{(1)}\tau _c`$ is naturally interpreted as the correlation time of the process $`𝐕(t)`$. For instance, if $`f(\tau )`$ is approximately Gaussian, the expansion of $`𝐕(0)𝐕(\tau )`$ around $`\tau =0`$ gives an explicit formula for the correlation time, namely

$$\frac{1}{\tau _c}=\left[\frac{2}{\pi \sigma ^2N}\text{Tr}\left(\frac{d𝐕}{dt}\right)^2\right]^{1/2}.$$
(21)
In this case $`\tau _c^{(2)}`$ and $`\tau _c^{(3)}`$ are trivially related to $`\tau _c`$:

$`\tau _c^{(2)}`$ $`=`$ $`{\displaystyle \frac{2}{\pi }}\tau _c^2,`$ (22)
$`\tau _c^{(3)}`$ $`=`$ $`{\displaystyle \frac{2}{\pi }}\tau _c^3.`$ (23)
## III The mean-field XY-Hamiltonian

In this section we start the application of the perturbative/mean-field theory of Sect. II to a specific model. Consider the one-dimensional Hamiltonian

$$H=\frac{1}{2I}\underset{i=1}{\overset{N}{}}L_i^2+\frac{J}{2N}\underset{i,j=1}{\overset{N}{}}\left[1\mathrm{cos}(\theta _i\theta _j)\right].$$
(24)
This is the so-called mean-field XY-Hamiltonian. It represents a lattice of classical spins with infinite-range interactions. Each spin rotates in a plane and is therefore described by an angle $`0\theta _i<2\pi `$, and its conjugate angular momentum $`L_i`$, with $`i=1,\mathrm{},N`$. The constants $`I`$ and $`J`$ are the moment of inertia and the interaction strength, respectively. (Of course, one can also think of point particles of mass $`I`$ moving on a circle.)

The mean-field XY-Hamiltonian (HMF) has been extensively studied in the last few years (see  for a review). The reasons for the interest in this model are various. From a general point of view, the HMF can be considered the simplest prototype for complex systems with long-range interactions like galaxies and plasmas (in fact, the HMF is a descendant of the mass-sheet gravitational model ). But the HMF is also interesting for its anomalies, be them model-specific or not. Especially worth of mention are the long-lived quasi-equilibrium states observed in the ferromagnetic HMF. These states exhibit breakdown of ergodicity, anomalous diffusion, and non-Maxwell velocity distributions . The explanation of these unusual behaviors may require an extension of the standard statistical mechanics, e.g., along the lines proposed by Tsallis . (Interesting anomalies are also present in the anti-ferromagnetic HMF ).

The simplicity of the HMF makes possible a full analysis of its statistical properties either in the canonical  or microcanonical ensemble . If interactions are attractive ($`J>0`$), the system exhibits a ferromagnetic transition at the critical energy $`E_c=3JN/4`$. In the case $`J<0`$ there is no ordered phase with finite magnetization, and, for not too low temperatures, the system behaves like the disordered phase of the $`J>0`$ case. However, at small energies, some kind of order appears, leading to a complex dynamical behavior .

Lyapunov exponents have been studied with detail, both numerically and analytically. For the ferromagnetic case, the simulations  show that in the magnetized phase the Lyapunov exponent remains finite in the thermodynamic limit $`N\mathrm{}`$. In contrast, if $`E>E_c`$, $`\lambda `$ goes to zero when $`N\mathrm{}`$; the same behavior is also observed in the antiferromagnetic case in all the energy range . The first theoretical studies were due to Firpo , who derived analytical expressions using the geometric method  (see also the discussion in Ref. ). Scaling laws in the high-energy regimes were derived as well using a random-matrix approach . At low temperatures the predictions of the geometric method for the HMF’s Lyapunov exponent are not satisfactory. There is solid evidence that the geometric method predicts wrong scaling laws at low temperatures for both ferromagnetic and antiferromagnetic interactions . This means that there are important gaps in the theoretical description of the Lyapunov exponent of the HMF (and of many-particle systems, in general). Further studies are necessary for understanding the precise domain of validity of the existing theories: this knowledge will be used to make the corresponding improvements.

Going back to the stochastic approach, it was proven in  that the mean-field approximation is exact in the HMF (this agrees with the supersymmetric analysis of Tănase-Nicola and Kurchan ). In the disordered phases of the HMF the fluctuations are small (see below), so, it is expected that the second-order perturbative approximation will work well, irrespective of the sign of $`J`$.

In order to test the stochastic approach, in its mean-field second-order-perturbative version, one has to calculate the average, variance, and correlation function of the Hessian $`𝐕(t)`$, i.e., Eqs. (17), (18) and (20). These ingredients are microcanonical averages of the appropriate observables. In the disordered regimes we will consider, microcanonical and canonical averages are equivalent (to leading order in $`N`$), so, we will prefer the simpler canonical averaging. Anyway, we will verify numerically that canonical and time averages coincide.

### A Calculation of $`\mu `$

Before embarking in the calculation of canonical averages, it is convenient to write the Hamiltonian (24) in the simplified form:

$$H=\frac{1}{2I}\underset{i=1}{\overset{N}{}}L_i^2+\frac{JN}{2}\left(1m^2\right),$$
(25)
where we have introduced the magnetization per particle

$$𝐦=\frac{1}{N}\underset{i=1}{\overset{N}{}}\widehat{𝐫}_i,$$
(26)
with

$$\widehat{𝐫}_i=(\mathrm{cos}\theta _i,\mathrm{sin}\theta _i).$$
(27)
In terms of $`𝐦`$ and $`\{\widehat{𝐫}_i\}`$ the elements of the Hessian matrix read

$`V_{ii}`$ $`=`$ $`J(𝐦\widehat{𝐫}_i{\displaystyle \frac{1}{N}}),`$ (28)
$`V_{ij}`$ $`=`$ $`J\widehat{𝐫}_i\widehat{𝐫}_j,ij.`$ (29)
Then one has

$$\mu /J=m^2\frac{1}{N}.$$
(30)
Our next objective is the average $`m^2`$. The interacting part of the canonical partition function reads

$$Z(\beta )𝑑\theta \text{e}^{\beta JNm^2/2}.$$
(31)
The integration over the angles can be reduced to an integration over the possible magnetizations:

$$Z(\beta )𝑑𝐦g(𝐦)\text{e}^{\beta JNm^2/2}.$$
(32)
Here $`g(𝐦)`$ represents the density of states in $`\theta `$-space with magnetization $`𝐦`$. The problem of finding $`g(𝐦)`$ has a long history and is known as Pearson’s random-walk problem . There are no simple closed expressions for $`g(𝐦)`$, but as we are interested in the disordered phases, where $`m`$ is of the order of $`1/\sqrt{N}`$, it suffices to use the central-limit approximation:

$$g(𝐦)\text{e}^{Nm^2}.$$
(33)
The relative error of this approximation is

$$1\frac{1}{2N}\left(12Nm^2+\frac{1}{2}N^2m^4\right)+\mathrm{}.$$
(34)
So, we can safely use expression (33) if $`N1`$ and the relevant configurations are such that $`Nm^21`$. The last condition is synonymous with disorder. It is always satisfied in the antiferromagnetic case. In the case $`J>0`$ it is satisfied if the energy is above the critical value $`\epsilon =3J/4`$, but not very close to it.

Putting together the formulas (32) and (33) we get the probability distribution of $`𝐦`$. To leading order in $`1/N`$ we have

$$P(𝐦)\mathrm{exp}\left[N\left(1\frac{\beta J}{2}\right)m^2\right].$$
(35)
And then we arrive at

$$m^2\frac{1}{N\left(1\beta J/2\right)}$$
(36)
(here and in the following, $``$ means “equal to leading order in $`1/N`$”). This expression, together with (25) and the equipartition theorem, $`kT=L_i^2/I`$, allows us to obtain the relationship between the temperature $`T`$ and the energy per particle $`\epsilon =E/N`$:

$$kT2\epsilon J.$$
(37)
As a function of $`\epsilon `$, the average $`\mu `$ reads

$$\mu /J\frac{1}{N\left(4\epsilon /J3\right)}.$$
(38)
Combining Eqs. (30) and (36) one sees that $`\mu `$ is a finite-temperature correction to the $`T=\mathrm{}`$ magnetization. So, for fixed $`N`$, $`\mu `$ goes to zero as the energy is increased. Equation (38) coincides with the result obtained by Firpo using a different technique .

### B Calculation of $`\sigma ^2`$

Given that all degrees of freedom are statistically equivalent, the definition of $`\sigma ^2`$, Eq. (18), can be expressed as

$$\sigma ^2=\left[\left(\delta 𝐕\right)^2\right]_{11}=\underset{j=1}{\overset{N}{}}\left(\delta V_{1j}\right)^2,$$
(39)
the rightmost identity following from the symmetry of $`𝐕`$. The abovementioned equivalence can be used once more, together with the definition of $`V_{ij}`$, to obtain

$`{\displaystyle \frac{\sigma ^2}{J^2}}`$ $`=`$ $`\left(𝐦\widehat{𝐫}_1\right)^2𝐦\widehat{𝐫}_1^2`$ (41)
$`+{\displaystyle \frac{1}{N}}\left[\left(\widehat{𝐫}_1\widehat{𝐫}_2\right)^2\widehat{𝐫}_1\widehat{𝐫}_2^2\right].`$
To proceed with the evaluation of $`\sigma ^2`$ one needs the probability distributions for two and three particles, $`P_2(\theta _1,\theta _2)`$ and $`P_3(\theta _1,\theta _2,\theta _3)`$. Consider first $`P_2`$, which is just

$$P_2(\theta _1,\theta _2)𝑑\theta _3\mathrm{}𝑑\theta _N\text{e}^{\beta JNm^2/2}.$$
(42)
To do the integrations we split the magnetization into two parts: the contribution from particles 1 and 2, and the remainder:

$$𝐦=𝐦_{12}+𝐦_{3N},$$
(43)
where

$`𝐦_{12}`$ $`=`$ $`{\displaystyle \frac{1}{N}}\left(\widehat{𝐫}_1+\widehat{𝐫}_2\right),`$ (44)
$`𝐦_{3N}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=3}{\overset{N}{}}}\widehat{𝐫}_i.`$ (45)
We outline the final steps. Like done before, we introduce a density of configurations $`g(𝐦_{3N})`$:

$$P_2(\theta _1,\theta _2)𝑑𝐦_{3N}g(𝐦_{3N})\mathrm{exp}\left[\frac{\beta JN}{2}(𝐦_{12}+𝐦_{3N})^2\right].$$
(46)
Now invoke the central limit theorem to approximate $`g(𝐦_{3N})`$. Switch to polar coordinates. Integrate over the polar angle, and then over the modulus $`m_{3N}`$ (the upper limit of integration can be extended to infinity). The final result is

$$P_2(\theta _1,\theta _2)A\mathrm{exp}\left[\alpha \left(\widehat{𝐫}_1\widehat{𝐫}_2\right)\right],$$
(47)
where $`A`$ is a normalization constant and

$$\alpha =\frac{2}{N\left(4\epsilon /J3\right)}.$$
(48)
The parameter $`\alpha `$ is the relevant perturbative quantity in this problem. It is small whenever $`N1`$ and the energy exceeds the transition value by a finite, large enough, amount. (Notice that $`\alpha =2\mu /J`$.)

Proceeding in the same way one also finds the three-particle distribution function:

$$P_3(\theta _1,\theta _2,\theta _3)B\mathrm{exp}\left[\alpha \left(\widehat{𝐫}_1\widehat{𝐫}_2+\widehat{𝐫}_2\widehat{𝐫}_3+\widehat{𝐫}_3\widehat{𝐫}_1\right)\right].$$
(49)
With the distribution functions $`P_2`$ and $`P_3`$ in our hands, we go back to the calculation of $`\sigma ^2`$, Eq. (41). The moments of $`\widehat{𝐫}_1\widehat{𝐫}_2`$ are immediate:

$`\widehat{𝐫}_1\widehat{𝐫}_2`$ $``$ $`{\displaystyle \frac{\alpha }{2}},`$ (50)
$`\left(\widehat{𝐫}_1\widehat{𝐫}_2\right)^2`$ $``$ $`{\displaystyle \frac{1}{2}}.`$ (51)
In addition one has

$$𝐦\widehat{𝐫}_1\frac{1}{N}+\widehat{𝐫}_1\widehat{𝐫}_2=𝒪(1/N).$$
(52)
Then it is easy to verify that Eq. (41) becomes

$$\frac{\sigma ^2}{J^2}\left(𝐦\widehat{𝐫}_1\right)^2+\frac{1}{2N}.$$
(53)
Splitting $`𝐦`$ into its $`N`$ parts, and using that particles are statistically equivalent, we arrive at

$$\frac{\sigma ^2}{J^2}\frac{1}{N}+\left(\widehat{𝐫}_1\widehat{𝐫}_2\right)\left(\widehat{𝐫}_1\widehat{𝐫}_3\right).$$
(54)
The three-body average is quickly done (with the help of Mathematica ) by expanding $`P_3`$ in powers of $`\alpha `$:

$$\left(\widehat{𝐫}_1\widehat{𝐫}_2\right)\left(\widehat{𝐫}_1\widehat{𝐫}_3\right)\frac{\alpha }{4}.$$
(55)
So, the final expression for $`\sigma ^2`$ is

$$\sigma ^2\frac{J^2}{N}\left(1+\frac{N\alpha }{4}\right).$$
(56)
### C Correlation function and correlation time

Above the transition, the relative importance of the interactions decreases with increasing $`\epsilon `$ and $`N`$, and the dynamics is dominated by the kinetic part of the Hamiltonian. The picture is that of particles rotating almost freely during times which are long as compared to the mean rotation period. If the system is in equilibrium, the dynamics can be modeled by the free-motion equations

$$\theta _k\theta _k(t=0)+L_kt/I,1kN,$$
(57)
where $`\{\theta _k(t=0)\}`$ and $`\{L_k\}`$ are independent random variables with uniform and Maxwell distributions, respectively. This is a first approximation valid only during short times. Then it is easy to show that the correlation functions of the elements of the Hessian are directly related to the characteristic function of the momentum distribution, i.e.,

$$2\mathrm{cos}\left(\theta _k\theta _j\right)(\tau )\mathrm{cos}\left(\theta _k\theta _j\right)(0)\left|\mathrm{exp}\left(iL_k\tau /I\right)\right|^2.$$
(58)
In equilibrium the characteristic function is Gaussian:

$$\left|\mathrm{exp}\left(iL_k\tau /I\right)\right|^2=\mathrm{exp}\left[\left(\tau /\tau _{}\right)^2\right],$$
(59)
with

$$\tau _{}=\sqrt{\frac{I}{kT}}.$$
(60)
After using the definition (19) we arrive at the simplest estimate of the correlation time in the high-energy regime

$$\tau _c\sqrt{\frac{\pi I}{4kT}}.$$
(61)
The correlation time $`\tau _c`$ is of the order of the mean period of rotation. It is independent of the system size because it is not directly associated with interactions. Of course, interactions are responsible for the Maxwell equilibrium distributions, and for the extinction of the ballistic regime and its substitution by a random-walk one. However, the time scales involved in these processes are much longer than $`\tau _c`$: the first one is the relaxation time for the one-body momentum distribution, the second is the momentum correlation time. Both times grow with $`N`$, and are not related with $`\tau _c`$.

A more precise estimate for $`\tau _c`$ can be derived if the correlation function is indeed Gaussian. Then the correlation time is that given by Eq. (21). We will not detail the calculation of the canonical averages of Eq. (21) because they are very similar to those of Sect. III B. We just show the result:

$$\tau _c\sqrt{\frac{\pi I}{4kT}\left(\frac{1+N\alpha /4}{1+N\alpha /8}\right)},$$
(62)
very close to the simple estimate of Eq. (61).

### D The Lyapunov exponent

Gathering the results of previous sections, one can construct the $`3\times 3`$ matrix $`𝚲`$ associated to the average propagator for the HMF. The Lyapunov exponent is extracted from the eigenvalue of $`𝚲`$ with the largest real part. Though the general expression for $`\lambda `$ might be written down explicitly, its content would not justify its extension. Notwithstanding, it is worth exhibiting the leading term in $`1/N`$. Notice that both $`\mu `$ and $`\sigma ^2`$ are of order $`1/N`$, and that $`\tau _c^{(k)}=𝒪(1)`$ (once the correlation function is assumed Gaussian). Then one gets

$$\lambda \left(\frac{\sigma ^2\tau _c}{2}\right)^{1/3}+𝒪\left(N^{2/3}\right)=𝒪\left(N^{1/3}\right),$$
(63)
with $`\sigma ^2`$ and $`\tau _c`$ given by Eqs. (56) and (62). The absence of $`\mu `$ in the leading-order expression for $`\lambda `$ is a reflection of the fact that in the disordered phases of the HMF, fluctuations are much larger than the average, i.e. $`\sigma \mu `$, and dominate the tangent dynamics.

## IV Numerical studies

For comparing the theory with simulations, we will use the data existing in the literature  as well as some additional data generated by us. We give a succinct description of how our simulations were made. Hamilton’s equations of the $`N`$-particle system were evolved using the Neri-Yoshida fourth-order symplectic algorithm . The time step was fixed through all simulations to $`\mathrm{\Delta }t=0.1`$ (units are such that $`I=1`$ and $`|J|=1`$). Initial conditions were chosen randomly: angles with a uniform distribution in $`[0,2\pi ]`$, and angular momenta Gaussian distributed. Then all momenta were shifted by a fixed amount to set to zero the total momentum (a constant of motion). Finally, velocities were multiplied by an appropriate factor to fix the total energy at a chosen value $`N\epsilon `$. The initial one-body distributions generated in this way are close to their equilibrium values in the disordered regimes defined by $`\epsilon >3J/4`$ and $`J>0`$, or $`\epsilon >0`$ and $`J<0`$. Anyway, before doing any “measurements” the systems were allowed to relax during a time we call $`t_{\mathrm{eq}}`$, typically $`t_{\mathrm{eq}}[1000,10000]`$. Concerning Lyapunov exponents, we used Benettin’s standard algorithm . The initial Euclidian distance between a trajectory and its companion were set to $`d_0=10^6`$. We used two alternative renormalization procedures: (a) The distance vector was compressed each time the distance exceeded $`d_0f`$; usually with $`f=10`$ or $`f=20`$; (b) Renormalizations were made at equally spaced instants, the interval between successive renormalizations corresponding, in average, to expansion factors $`f10,20`$.

Numerically, the Lyapunov exponent is estimated by averaging over initial conditions which are propagated during a finite time:

$$\lambda (t)\frac{1}{2t}\mathrm{ln}|\xi (t;x_0,\xi _0)|^2.$$
(64)
We usually considered 10 pairs of randomly-chosen initial conditions. Each pair was propagated until a time $`t=t_{\mathrm{prop}}`$, when we judged that $`(1/t)\mathrm{ln}|\xi |`$ had converged to a limiting value. Typically, $`t_{\mathrm{prop}}[1000,5000]`$. Remember, however, that our theoretical scheme commutes logarithm and average. In order to test this approximation, we also computed the average

$$\lambda ^{}(t)\frac{1}{2t}\mathrm{ln}|\xi (t;x_0,\xi _0)|^2.$$
(65)
In many cases we run tests to verify that our numerical results are robust against suitable changes of initial conditions, renormalization procedure, or the set of parameters $`\{t_{\mathrm{eq}},t_{\mathrm{prop}},\mathrm{\Delta }t,d_0,f\}`$.

### A Ferromagnetic case

Figure 1 shows the largest Lyapunov exponent of the ferromagnetic HMF as a function of system size $`N`$, for some selected energies $`\epsilon `$. The symbols correspond to the simulations of Refs.  ($`\epsilon =5.0`$) and  ($`\epsilon =1.0,10.0,50.0`$). We have considered large particle numbers $`(N100)`$ and energies well above the transition ($`\epsilon >1.0`$) to ensure (i) the validity of the approximations invoked in calculating the averages of Sect. III, and also (ii) to guarantee that we are in a disordered, quasi-ballistic regime. Full lines correspond to the theoretical Lyapunov exponent obtained by diagonalizing the $`3\times 3`$ matrix of Sect. (II), but we could as well have used the asymptotic expression of Eq. (63) (minor differences would only be visible in the case $`\epsilon =1`$). Dotted lines correspond to the geometric prediction, as calculated by Firpo, i.e., Eqs. (1-3) and (21-22) of Ref. . An inspection of Fig. 1 allows us to verify that there is a satisfactory agreement between theory and simulation, especially for $`N500`$. For larger systems numerical data deviate upward from theory. (Nobre and Tsallis also observed deviations from the $`N^{1/3}`$ law in the three-dimensional version of the XY-Hamiltonian .) This is the opposite to natural expectation “the larger the system, the better the approximations”. We believe that the reason for this deviation may be that the system has not reached microcanonical equilibrium – as far as the Lyapunov exponent is concerned.

It is well known that the HMF relaxes very slowly, e.g., Latora, Rapisarda, and Ruffo  noticed that for $`\epsilon =1`$, the system may get trapped in quasi-equilibrium states similar to those observed below the transition. A very slow relaxation led the same authors to conclude that dynamics is ballistic for $`\epsilon =5.0`$, $`N=1000`$, $`t<10^5`$ . If slow relaxation is ruling, the numerical simulations are likely to provide “quasi-equilibrium” Lyapunov exponents, rather than the microcanonical ones. Actually, as energies and/or system sizes are increased, the system becomes progressively more integrable and deviations from microcanonical predictions will be necessarily observed.

At this point it is desirable to verify if the canonical calculations of $`\mu `$, $`\sigma ^2`$, and $`\tau ^{(k)}`$, reproduce the corresponding dynamical averages. In Fig. 2 we show numerical values for $`\mu `$ and $`\sigma ^2`$, obtained through averaging over 1000 initial conditions, and in the time window $`[5000,10000]`$. One sees that $`\sigma ^2`$ is very close to the theoretical value, within a few percent, but the results for $`\mu `$ show larger deviations. We can understand why $`\sigma ^2`$ is so close to its expectation value and $`\mu `$ is not: The initial value of $`\sigma ^2`$ is already very close to the canonical value, within a few percent for the cases considered. On the contrary, the initial value for $`\mu `$ is zero, far from its equilibrium value; the slow growth of $`\mu `$ is one of the aspects of the slow relaxation to equilibrium. In any case we verified (graphics not shown) that the errors in $`\mu `$ and $`\sigma ^2`$ do not account for the deviations of Lyapunov exponents, i.e., feeding the theory with the numerical $`\mu `$ and $`\sigma ^2`$ does not improve the agreement between theoretical and simulated $`\lambda `$’s.

In order to discuss the correlation times, let us now look at the correlation functions of Fig. 3. Theory and simulations agree almost perfectly in the central part of the distributions. There are long tails, which cannot be appreciated in the figure but might be responsible for larger-than-predicted correlation times $`\tau ^{(k)}`$. We recall that the quantities $`\tau ^{(k)}`$ are moments of the correlation function $`f(\tau )`$, and as such, very sensitive to the precise shape of the tails We verified that a pure numerical evaluation of $`\tau ^{(k)}`$ fails due to poor convergence of the corresponding integrals – this effect is stronger for $`k=1,2`$. So, at present we cannot say with certainty if the theoretical estimates for $`\tau ^{(k)}`$, which assume a Gaussian correlation function, agree with the corresponding dynamical ones. To resolve this issue one should make a very careful study of the tails of the correlation functions. Though desirable such analysis exceeds the scope of this paper. (Interesting information about correlation functions, and, in general, about the geometric method, can be found in Refs. .)

In passing, let us comment that, for the theory to work, $`f(\tau )`$ must decay fast enough, e.g., a correlation function like $`\mathrm{\Gamma }(\tau )=\mathrm{sin}\omega t/\omega t`$  would explode the second cumulant.

To conclude the analysis of the ferromagnetic case, a possible difference between $`\lambda `$ and $`\lambda ^{}`$ \[Eqs. (64) and (65)\] cannot be the explanation for the observed deviations between theory and simulations because it is expected that $`\lambda <\lambda ^{}`$ (see below, and Ref. ).

### B Antiferromagnetic case

Figure 4 shows the results of numerical simulations for $`\lambda `$ and $`\lambda ^{}`$. (The numerical results for $`\lambda `$, in the case $`N=100`$, are consistent with those obtained in Ref. .) Though not shown, as temperature goes to zero ($`\epsilon 1/2`$), the numerical $`\lambda `$’s go to zero too, the empirical scaling law being $`\lambda T^{1/2}`$ . Notice that the relative differences between $`\lambda `$ and $`\lambda ^{}`$ are indeed small, and decrease with increasing $`N`$. We have not analyzed the behavior of $`\lambda ^{}`$ at very low temperatures; in particular, we cannot say whether the scaling law $`\lambda ^{}T^{1/2}`$ is verified or not.

For high energies, i.e., $`\epsilon >1`$, our theoretical predictions agree reasonably well with simulations. However, when the energy tends to its ground-state value, $`\epsilon 1/2`$, our theoretical Lyapunov exponent diverges dramatically from experiments. These behaviors can be easily understood by looking at the energy dependence of the perturbative parameter $`\sigma \tau `$. Figure (5) shows that, for a fixed system size, $`\sigma ^2`$ remains bounded for all energies, growing by a factor of two as $`T`$ goes from zero to infinity. On the other side, $`\tau `$ scales with temperature like $`1/\sqrt{T}`$, and then, $`\sigma \tau 1/\sqrt{T}`$. This means that the theory must fail as $`T0`$. The divergence of the Kubo number at $`T=0`$ is a manifestation of the complexity of the dynamics close to the antiferromagnetic ground state . Notwithstanding the overall agreement at high temperatures, we must point out a non-resolved source of annoyance. At a fixed energy the Kubo parameter scales with the system size like $`\sigma \tau 1/\sqrt{N}`$. Consequently the relative error of the theoretical prediction should decrease with increasing $`N`$, like $`1/\sqrt{N}`$. However, e.g., if we take $`\epsilon =7.5`$, and compare theory versus $`\lambda ^{}`$, we see that the relative error is about 11%, insensitive to variations of $`N`$ from 100 to 2500.

For the sake of completeness, we exhibit in Figs. (5) and (6) the comparison of canonical and dynamical averages for $`\mu `$, $`\sigma `$, and $`f(\tau )`$: No significant differences can be seen. (Again, the tails of the correlation functions may deserve a deeper analysis.)

Figure (4) also displays the predictions of the geometric method . In the high-temperature regime, the relative error of the geometric prediction is somewhat larger than ours. At low temperatures, the geometric method predicts a qualitatively correct behavior, $`\lambda 0`$, but with the wrong scaling law: $`\lambda T^2`$  instead of $`\lambda T^{1/2}`$ .

## V Summary and concluding remarks

In a previous paper  we proposed a theoretical approach to the determination of the largest Lyapunov exponent of many-particle Hamiltonian systems. Despite several crude approximations that had to be made, it was not completely clear whether the resulting perturbative/mean-field script could be carried out for a specific system, one of the obstacles being the analytical determination of the correlation functions.

Now we have verified that the stochastic recipe can indeed be executed and works satisfactorily in the quasi-ballistic regimes of the mean-field XY-Hamiltonians. These systems are especially useful for testing the theory because they make the “mean-field” diagonalization exact. However, three additional sources of error have to be considered. First of all, we made the simplification of taking the logarithm out of the average \[Eq. (6)\]. That is, using a spin-glass analogy, we are estimating an “annealed” Lyapunov exponent $`\lambda ^{}`$ instead of the usual “quenched” $`\lambda `$ . Quenched averages can be calculated, as shown by Tănase-Nicola and Kurchan , but they require more sophisticated tools, like the replica trick or the supersymmetric approach. In the cases considered in this paper we verified that the difference between both Lyapunov exponents is very small. Then it was not necessary to go beyond the annealed averaging.

The second approximation, and most important one, is the truncation of the cumulant expansion at the second order. By doing so, we introduce a relative error of the order of the Kubo number $`\sigma \tau _c`$. This explains the failure of the theory in the low-temperature regime of the antiferromagnetic XY-Hamiltonian, given that the amplitude of the fluctuations $`\sigma `$ remains bounded but the thermal time $`\tau _c`$ diverges. It is not clear at present if the theory can be improved to account for this regime. Further studies are necessary. In particular we have to understand what kind of dynamical correlations develop as the system approaches the ground state. Another problem that requires an answer: We have not observed that the agreement between theory and simulations betters when the size of the system is increased, for fixed energy. This is in contradiction with the Kubo number decreasing like $`1/\sqrt{N}`$.

The difficulties of the preceding paragraph may be related to the third approximation, which concerns the interaction representation we have used. In order to isolate the fluctuations, one should work in the representation associated to the average Hessian. Instead, we have chosen the free-rotator representation. That is, in constructing the “free propagator”, we neglected the average Hessian of the interaction potential. As long as the motion is quasi-ballistic, this approximation seems justified. It will be certainly wrong, e.g., in the low-temperature phase of the ferromagnetic XY, where the quadratic part of the interaction is really important. We are currently working on the implementation of the best interaction representation for the infinite-range XY models. The results will be presented in a forthcoming paper.

The mean-field diagonalization is exact for the systems considered here. In other cases it represents just a truncation of the basis for diagonalizing $`\widehat{\mathrm{\Lambda }}`$ \[Eq. (11)\]. Notice, however, that we are dealing with a hermitian problem, i.e., finding the largest eigenvalue of $`\xi \xi ^T(t)`$ \[Eq. (10)\]. Truncation of the basis will produce a lower bound to the largest Lyapunov exponent (provided that $`\widehat{\mathrm{\Lambda }}`$ is calculated accurately). So, this problem is analogous to finding the ground-state energy of a quantum Hamiltonian. Any small subspace can be considered, the choice being guided by the special characteristics of the system under study.

Like its quantum analogue, the Lyapunov problem seems reluctant to admit simple general solutions. Each class of systems may require special consideration.

## Acknowledgements

We are grateful to F. V. Roig and C. Tsallis for fruitful discussions. We acknowledge Brazilian Agencies CNPq, FAPERJ and PRONEX for financial support.