Lyapunov Exponent and Criticality in the Hamiltonian Mean Field Model

TL;DR

This paper studies how the largest Lyapunov exponent behaves in a Hamiltonian mean field model at equilibrium, revealing critical behavior and chaos transition related to phase changes and resonant particle interactions.

Contribution

It provides numerical confirmation of a critical exponent for the Lyapunov exponent and discusses the limitations of geometrical analytical estimates in the model.

Findings

Existence of a critical exponent for the Lyapunov exponent at phase transition.

Chaos persists in the magnetized phase even in the thermodynamic limit.

Transition from weak to strong chaos coincides with the onset of diffusive center of mass motion.

Abstract

We investigate the dependence of the largest Lyapunov exponent of a $N$-particle self-gravitating ring model at equilibrium with respect to the number of particles and its dependence on energy. This model has a continuous phase-transition from a ferromagnetic to homogeneous phase, and we numerically confirm with large scale simulations the existence of a critical exponent associated to the largest Lyapunov exponent, although at variance with the theoretical estimate. The existence of chaos in the magnetized state evidenced by a positive Lyapunov exponent, even in the thermodynamic limit, is explained by the resonant coupling of individual particle oscillations to the diffusive motion of the center of mass of the system due to the thermal excitation of a classical Goldstone mode. The transition from "weak" to "strong" chaos occurs at the onset of the diffusive motion of the center of mass of the non-homogeneous equilibrium state, as expected. We also discuss thoroughly for the model the validity and limits of a geometrical approach for their analytical estimate.

Equations

$$H=\sum _{i=1}^N\frac{p_i^2}{2m}+\frac{1}{2N}\sum _{i,j=1}^N\left[1-\cos ({\theta }_i-{\theta }_j)\right],$$

$$\mathbf{M}=\left(M_x,M_y\right)=\frac{1}{N}\sum _{i=1}^N(\cos {\theta }_i,\sin {\theta }_i).$$

$${\dot{\theta }}_i$$

$${\dot{p}}_i$$

$$f_{\mathrm{eq}}(p,\theta )=\frac{\sqrt{\beta }}{(2\pi {)}^{3/2}{\mathrm{I}}_0(\beta )}{e}^{-\beta \left(p^2/2-M\cos (\theta )\right)},$$

$$M=\frac{{\mathrm{I}}_1(\beta M)}{{\mathrm{I}}_0(\beta M)},$$

$$\mathbf{x}\equiv ({\mathbf{x}}_1,{\mathbf{x}}_2,....{\mathbf{x}}_{\mathbf{n}}).$$

$$\frac{d\mathbf{x}(\mathbf{t})}{dt}=\mathbf{F}(\mathbf{x}(\mathbf{t})).$$

$$\mathbf{w}=(\delta x_1,\delta x_2,\dots ,\delta x_n).$$

$$\frac{d\mathbf{w}}{dt}=\mathbf{J}(\mathbf{x}(t))\mathbf{w},$$

$$\lambda \equiv \underset{t\to \infty }{\lim }\frac{1}{t}\ln \left(\frac{||\mathbf{w}(t)||}{||\mathbf{w}(0)||}\right).$$

$${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_N.$$

$${\lambda }_i=-{\lambda }_{2N-i+1},{\lambda }_{N+1}={\lambda }_N=0.$$

$$\frac{d\mathbf{w}}{dt}=\mathbf{J}{|}_{x=x^{\ast }}\mathbf{w},$$

$${\lambda }_1=\frac{1}{KT}\sum _{k=1}^K\ln \frac{||{\mathbf{\delta }}_k||}{ϵ},$$

$$\dot{\delta {\theta }_i}$$

$$\dot{\delta p_i}$$

$$\delta M_x(t)\equiv -\frac{1}{N}\sum _{j=1}^N\delta {\theta }_j(t)\sin {\theta }_j^{\ast }(t),$$

$$\delta M_y(t)\equiv \frac{1}{N}\sum _{j=1}^N\delta {\theta }_j(t)\cos {\theta }_j^{\ast }(t).$$

$$H=\sum _{i=1}^N\frac{p_i^2}{2}+V({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N),$$

$$\frac{d^2{\mathbf{\xi }}_i(t)}{d{t}^2}+\sum _{j=1}^N\frac{{\partial }^{2}V}{\partial {\mathbf{r}}_i\partial {\mathbf{r}}_j}\cdot {\mathbf{\xi }}_j(t)=0.$$

$$\lambda =\underset{t\to \infty }{\lim }\frac{1}{2t}\ln \left(\frac{||\mathbf{\xi }(t)|{|}^2}{||\mathbf{\xi }(0)||}\right),$$

$$||\mathbf{\xi }(t)||\equiv \sqrt{{\mathbf{\xi }}_1(t{)}^2+\cdots +{\mathbf{\xi }}_N(t)}.$$

$$\frac{d^2\psi (t)}{d{t}^2}+k(t)\psi (t),$$

$$\lambda =\underset{t\to \infty }{\lim }\frac{1}{t}\frac{⟨\psi (t)⟩}{⟨\psi (0)⟩}.$$

$${\lambda }_1=\frac{\Lambda }{2}-\frac{2{\kappa }_0}{3\Lambda },$$

$$\Lambda ={\left(2{\sigma }_k^2\tau +\sqrt{\frac{64}{27}{\kappa }_0^3+4{\sigma }_k^4{\tau }^2}\right)}^{1/3},$$

$$\tau =\frac{\pi \sqrt{{\kappa }_0}}{2\sqrt{{\kappa }_0}\sqrt{{\kappa }_0+{\sigma }_k}+\pi {\sigma }_k},$$

$${\sigma }_k^2=⟨{\delta }^2{k}_R{⟩}_c{\left(1+\frac{{\beta }^2}{2}⟨{\delta }^2{k}_R{⟩}_c\right)}^{-1},$$

$$⟨{\delta }^2{k}_R{⟩}_c=4M\frac{\partial M}{\partial \beta },$$

$${\lambda }_1\propto (e_c-e{)}^{\xi },$$

$${\rho }_{eq}(\theta )=C_N\exp \left(\beta M\cos (\theta )\right),$$

Full text

I Introduction

Many body systems with long range interactions are known to have several properties that set them apart from more “usual” systems with short range interactions, such as ensemble inequivalence, negative heat capacity (with no second law violation), anomalous diffusion and non-Gaussian (quasi-) stationary states booklri . An interparticle interaction potential is said to be long ranged if it decays at large distances as $r^{-\alpha}$ with $\alpha\leq d$, $d$ the spatial dimension, with a consequence that the total potential energy increasing superlinearly with volume campa ; dauxois . Some important physical system with long range interactions are non-neutral plasmas rizzato , self-gravitating systems tanu , vortices in two-dimensional turbulent hydrodynamics venaille and free electron laser anton1 . Simplified models were also largely considered in the literature and allowed a better understanding of the statistical mechanics of equilibrium and non-equilibrium of systems with long range interactions, such as one and two-dimensional self-gravitating systems mila ; teles , the Hamiltonian Mean Field (HMF) and self-gravitating ring models ruffoHMF ; sota .

Much progress in the understanding of the relaxation properties in many-particle systems with long range forces came from numerical simulations of model systems steiner ; rl1 ; rl2 ; rl3 ; rl4 ; rl5 ; rl6 ; rl7 ; rl8 . Although the scaling of the relaxation time to equilibrium with $N$ depends on the type of system and spatial dimension scaling ; chris , as a common feature it diverges with $N$, and as a consequences it never attains thermodynamic equilibrium for $N\rightarrow\infty$. In many cases this relaxation time is sufficiently large that even for finite $N$ it can be considered infinite for practical purposes. If the equilibrium is reached, then its properties can be studied using the usual techniques of equilibrium statistical mechanics booklri ; campa ; dauxois ; levin .

Simplified models have been important in the study of the intricate interplay between chaotic dynamics, ergodic properties and statistical mechanics of systems with long range interactions, while Lyapunov exponents has proven to be a useful tool in the study of chaos in dynamical systems ott and particularly also for long range systems vallejos1 ; firpo1 ; anteneodo1 . The precise determination of Lyapunov exponents is an intricate task and usually requires a great numerical effort with very long integration times, that can become prohibitive for a system with a very large number of particles. A prescription for their analytic estimation is therefore of great relevance. Casetti, Pettini and collaborators developed an analytic method to obtain the scaling behavior of the Largest Lyapunov Exponent (LLE) pettini1 ; pettini2 ; pettini3 , and applied to the HMF model by Firpo in Ref. firpo2 .

The HMF model has been widely studied in the literature as a prototype for observations of some dynamical features of long range interacting systems booklri . It consists of $N$ classical particles moving on a unit circle and globally coupled with Hamiltonian ruffoHMF :

[TABLE]

where $\theta_{i}$ and $p_{i}$ are the position angle on the circle and $p_{i}$ its conjugate momentum. The $1/N$ factor in the potential energy is the Kac factor introduced such that the total energy is extensive, and can be obtained from a change in the time unit. We can define by analogy the total magnetization and its components by:

[TABLE]

The equations of motion are then

[TABLE]

As a thermodynamic system this system is exactly solvable, i. e. its equilibrium partition function is obtained in closed form, and its equilibrium distribution function is given by ruffoHMF :

[TABLE]

where ${\rm I}_{k}$ is the modified Bessel function of the first kind with index $k$, and the origin for angles is chosen such that $M_{y}=0$ and $M_{x}=M$. For a given inverse temperature $\beta$, the magnetization is obtained from the equation:

[TABLE]

with $M\equiv||{\bf M}||$. The dependence of $M$ on temperature is thus obtained by solving Eq. (5). A second order phase transition occurs at the critical energy per particle $e_{c}=E_{c}/N=3/4$ and $T=1/\beta=0.5$ from a lower energy ferromagnetic phase to a higher energy phase with zero magnetization. Canonical or microcanonical ensembles are fully equivalent for the HMF model. Out of equilibrium phase transitions for this model were studied in some detail in amato3 .

In the present work we investigate the applicability of the geometrical approach of Refs pettini1 ; pettini2 by directly testing its underlying assumptions. We also consider the scaling with $N$ of the LLE for the HMF model at different energy ranges, and compare our numerical results to theoretical results and other similar numerical investigations, for larger values of $N$ than in previous studies. Particularly we confirm the existence of a new critical exponent corresponding for the LLE theoretically predicted in firpo2 although with a small deviation from the predicted value of the exponent.

This paper is structured as follows: in section II we briefly recall the theory of Lyapunov Exponents and the numerical determination of the LLE. In section III we present and discuss our results for the HMF model and we close in section IV with some concluding remarks.

II Lyapunov Exponents

A Lyapunov Exponent (LE) quantify how the dynamics of the system is sensible to small differences in the initial conditions. With this aim, let us define the vector formed by coordinates in a $n$-dimensional phase space:

[TABLE]

which we suppose satisfy a set of $n$ autonomous first-order differential equations:

[TABLE]

Equation (7) generates a flows in the phase space, and $F(\bf{x}(t))$ is the velocity field of the flow. In order to measure contraction or stretching in the neighborhood of $\bf{x}(t)$, we consider two different solutions of Eq. (7) $\bf{x}^{(1)}(t)$ and $\bf{x}^{(2)}(t)$ and the difference vector $\mathbf{w}\equiv{\bf{x}}^{(2)}(t)-{\bf{x}}^{(1)}(t)$:

[TABLE]

The evolution equation for $\mathbf{w}$ is then:

[TABLE]

with $\bf{J}$ the $N\times N$ Jacobian matrix of the flow. Assuming that the elements of $\bf{J}$ are continuous bounded functions of $t$ for $t\rightarrow\infty$, then the solutions of (9) grow no faster than $\exp(\lambda t)$, for some constant $\lambda$.

The Lyapunov Exponent for a given initial condition $\mathbf{w}(0)$ is defined by

[TABLE]

In a $n$-dimensional problem we have $n$ Lyapunov exponents, each one referring to the divergence degree of specific directions of the system. All of them form a set called Lyapunov Spectrum (LS), which usually are organized as:

[TABLE]

If the LLE $\lambda_{1}$ is positive then neighbor trajectories tend to diverge exponentially which implies a chaotic regime. Due to the Liouville theorem, the Lyapunov spectrum of a Hamiltonian system is, as those considered below, satisfies the relations (Pesin’s theorem):

[TABLE]

Another important result relates the Kolmogorov-Sinai entropy and the Lyapunov spectrum. The former measures exponential rate of information production in a dynamical system ott and according to Pesin’s theorem can be obtained as the sum of all positive Lyapunov exponents eckmann .

II.1 Numeric determination of the LLE

In the tangent map method, one considers the linearized form of Eq. (7) around the point $x=x^{*}$:

[TABLE]

where $\mathbf{J}$ is the Jacobian matrix of the vector function ${\mathbf{F}}({\bf x})$. One then solves the original nonlinear system in Eq. (7) and the linearized equations (13). The steps for determining the Lyapunov exponent are strelcyn ; parker :

For the nonlinear system (7) impose an initial condition $\mathbf{x}_{0}$, and an initial condition ${\mathbf{w}}_{0}=\boldsymbol{\delta}_{0}$ for the linearized equations, with $||\boldsymbol{\delta}||=\epsilon$ and $\epsilon\ll 1$. 2. 2.

Both differential equations are integrated for a time interval $T$. This results in $\mathbf{x}_{0}\rightarrow\mathbf{x}(T)$ and $\boldsymbol{\delta}_{0}\rightarrow\boldsymbol{\delta}_{1}\equiv{\mathbf{w}}(T)$; 3. 3.

After each integration interval $T$, normalize the corresponding difference vector $\boldsymbol{\delta}_{k}$ to $\epsilon$ and use the resulting vector as a new initial condition for solving the linearized equations; 4. 4.

The LLE is obtained from the average:

[TABLE]

where $K$ is chose in order to achieve convergence in the value of $\lambda_{1}$.

By considering a solution $(\theta_{i}^{*}(t),p_{i}^{*}(t))$ of the equations of motion of the HMF model, the linearized equations are obtained by plugging $\theta_{i}(t)=\theta_{i}^{*}(t)+\delta\theta_{i}(t)$ and $p_{i}(t)=p_{i}^{*}(t)+\delta p_{i}(t)$, with small $\delta\theta_{i}(t)$ and $\delta p_{i}(t)$, into Eq. (3):

[TABLE]

where the components of the magnetization are computed at the angles $\theta_{i}^{*}$ and are denoted $M_{x}^{*}$ and $M_{y}^{*}$ and

[TABLE]

Both sets of equation in Eq. (3) and Eq. (15) must be solved simultaneously.

In order to compute the LLE for very large values of $N$ the tangent map method was implemented in a parallel code in graphic processing units eu2 using a fourth-order symplectic integrator for both system yoshida . Figure 1 shows the results for the computation of the LLE for some different values of $N$ and energy per particle $e=0.5$. The error bars decrease rapidly with $N$, as expected. The right-panel of the same figure shows the that convergence is achieved for a total simulation time $t_{f}=10^{5}$. In all the results below we thus chose to use the same parameter values and twice larger a value for $t_{f}$ to ensure proper convergence in all cases.

II.2 An analytic estimate for the LLE

The LLE for the HMF model was investigated numerically by Yamaguchi yamaguchi , and then latter estimated by Latora, Rapisarda and Ruffo latora from a random matrix approach of Parisi and Vulpiani parisi , and by Firpo firpo2 using the differential geometry approach by Pettini and collaborators casetti ; caiani ; pettini2 . In the latter approach, the dynamics of the $N$ particle system is reformulated in the framework of Riemannian geometry, where the trajectories correspond to geodesics of an underlying metric. Chaos then comes from the instability of the geodesic flow, that at its turn depends on the properties of the curvature of the Riemannian manifold. Chaos can also result from a parametric instability of the fluctuation of the curvature along the system trajectory as represented in the manifold.

In order for the present paper to be self contained, we succinctly present here the the main results of the geometrical approach to the computation of Lyapunov exponents and chaos from Refs. pettini2 and casetti (where the reader can find more details). We consider a system of $N$ identical particles with unit mass and Hamiltonian:

[TABLE]

where ${\bf r}_{i}$ is the position vector of particle $i$, ${\bf p}_{i}$ its canonically conjugate momentum and $V$ the potential energy. Considering two solutions $\overline{\bf r}_{i}(t)$ and ${\bf r}_{i}(t)$, $i=1,\ldots,N$, initially close to one another: $\overline{\bf r}_{i}(t)={\bf r}_{i}(t)+\boldsymbol{\xi}_{i}(t)$. The linearized equations of motion for the small variations $\xi_{1}(t)$ are given by:

[TABLE]

The maximal Lyapunov exponent is then obtained from:

[TABLE]

where

[TABLE]

The norm $\psi=||\boldsymbol{\xi}(t)||$ satisfies the equation:

[TABLE]

where $k(t)$ is a stochastic process describing the time evolution of the curvature along a trajectory in phase space. Since the solutions of Eq. (21) are given in term of the averages $\langle\cdots\rangle$ over many realizations of the stochastic process, Eq. (19) assumes the form:

[TABLE]

The solution for the average is obtained from a perturbative expansion in the amplitude of the fluctuations of the stochastic process, which are small for large $N$. By introducing a smallness multiplicative parameter in the stochastic process as $k(t)\rightarrow\alpha k(t)$, where $\alpha\ll 1$, an supposing that fluctuations are delta correlated, a perturbative solution of Eq. (21) in powers of $\alpha$ can be obtained pettini2 ; kampen . In the case that $k(t)$ is a Gaussian process, this solution becomes exact.

An estimate of the LLE $\lambda_{1}$, with the extra assumption that the curvature along a trajectory is well represented by a Gaussian process, was obtained in Ref. casetti , with very good agreement with numerical results for the Fermi-Pasta-Ulam model and the 1D $XY$ model pettini2 . The LLE is given in then given by casetti :

[TABLE]

where

[TABLE]

[TABLE]

with $\kappa_{0}$ and $\sigma_{k}^{2}$ the average curvature and the fluctuations around its mean value, respectively, $\tau$ being a characteristic time for the stochastic process. For the HMF model, Firpo obtained a closed form expression for the quantities $\kappa_{0}$ and $\sigma_{k}$ firpo2 , such that the (Ricci) scalar curvature in the Riemannian manifold is given by $\kappa_{R}=M^{2}$. The next step consists to take $\kappa_{0}=\langle M^{2}\rangle_{\mu}$, i. e. the microcanonical average of $M^{2}$. The variance of the curvature fluctuations in the microcanonical ensemble was obtained in firpo2 as:

[TABLE]

where $\langle\delta^{2}k_{R}\rangle_{c}$ is the variance of the fluctuations in the canonical ensemble:

[TABLE]

with $M$ given by the solution of Eq. (4). It is worth noting that even if $k(t)$ is not Gaussian, the expressions above remain valid up to first order in $\alpha$.

Previous results showed that in the non-magnetized phase the Lyapunov exponent tends to zero as $N^{-1/3}$ obtained in numerical simulations in Ref. latora and predicted theoretically in firpo2 . In the ferromagnetic phase, a more complicate picture emerges. Manos and Ruffo observed numerically a transition from a weak to a strong chaoticity regimes at low energy manos , and related it to the time dependence of the phase of the magnetization vector, which becomes strongly time dependent around the same energy (a more detailed explanation of this point is given in Ref. rochamarcos ). A critical exponent for the LLE was predicted by Firpo firpo2 in the vicinity of the second order phase transition for $e<e_{c}$ in the form

[TABLE]

with an exponent $\xi=1/6$. Ginelli and collaborators obtained a different value $\xi=1/2$ from numerical results, the same critical behavior as the magnetization ginelli . Below we obtain a value of $\xi$ close to the theoretical value by considering much higher values of $N$.

III Results

The first point to consider is whether the fluctuations of the curvature, i. e. of $M^{2}$ for the HMF model, can be modeled by an uncorrelated Gaussian process, as considered in Refs firpo2 ; casetti . Figure 2 shows the distributions of the fluctuations of the curvature $\kappa_{R}=M^{2}$ for a few values of energy, and the correlation function for the fluctuations $\langle M(t_{0}-\tau)^{2}M(t_{0})^{2}\rangle-\langle M^{2}\rangle^{2}$ for a few energy values. In the ferromagnetic state the fluctuations are well described by a Gaussian distribution, but the correlation time, i. e. the time for correlations to be negligible, can be very large. The correlation time is small only at higher energies. In the homogeneous phase, correlations of the fluctuations of the curvature are also non-negligible, and their distribution is non-Gaussian quite close to an exponential function. In fact in this case it is more natural to expect that the fluctuations of the magnetization components are Gaussian rather than those of $M^{2}$, thus explaining the form of the distributions in Fig. 2e and 2g. The distributions for the values of $M$ are given in Fig. 3. Below the critical energy the distribution is Gaussian, while above the phase transition it is well described by a function of the form $bM\exp(-aM^{2})$, with $a$ and $b$ constants. In obtaining Eqs. (23–25) the central assumption was that fluctuations are delta correlated. This is clearly valid only for higher energies, where as shown below the predicted $N^{-1/3}$ scaling of the LLE is observed. Deviations from the theoretical predictions are thus expected for lower energies due to strong correlations in the fluctuations of the curvature.

Figure 4 shows the plot of the LLE $\lambda_{1}$ as a function of energy for some values of $N$, alongside the theoretical prediction of Ref. firpo2 . The parameters used in the numeric integration are $T=10.0$ and $\delta t=0.05$ which are used in all simulations below unless explicitly stated. In the ferromagnetic phase the theoretical results agree only qualitatively with numerical results, predicting a maximum of the LLE for an energy below the critical energy $e_{c}$, but not its position, and also that $\lambda_{1}$ goes to zero at the phase transition. The left panel in Fig. 5 shows a reasonable data collapse if the exponent are rescaled by $N^{-1/3}$, that nevertheless becomes not so good for energies closer to the phase transition as seen on the left panel of Fig. 5.

Figure 6 shows $\lambda_{1}$ as a function of $N$ for some energy energy values in the non-magnetized state. The predicted $N^{-1/3}$ scale is observed far from the phase transition. Nevertheless for higher values of $N$ we slowly approach the $N^{1/3}$ scaling as shown in Fig. 7 for the energy $e=0.8$.

This can be explained by the fact that the fluctuations of the Riemannian curvature are not delta correlated close to the phase transition as seen from the correlation functions in Fig. 8. It is also important to note that for non-Gaussian fluctuations the solution of Eq. (21) is only valid at order $\alpha^{2}$, and therefore is more accurate for smaller $\alpha$ and equivalently greater $N$.

For the ferromagnetic phase, Figure 9 shows the LLE $\lambda_{1}$ as a function of $N$ for a few energy values. The scaling of the LLE with $N$ is close to $N^{-1/3}$ for very low energies while it is much slower for energies above $e_{w}\approx 0.15$ and below the critical energy.

Manos and Ruffo studying the same system observed a transition from weak to strong chaos at the same energy $e_{w}$, such that below it the LLE is much smaller and scales as $N^{-1/3}$, while no results for the scaling of the LLE with $N$ were obtained for $e>e_{w}$ manos . The same authors using the generalized alignment indices method skokos showed that at this energy the fraction of chaotic orbits of the system increases rapidly from a very low (less than 1%) to a very large value (close to 100%). As a consequence, the convergence of the LLE to zero in the mean-field limit is non-uniform, which characterizes two distinct energy intervals. For $e<e_{w}$ (weak chaos) the LLE rapidly tends to zero, while having a significant positive value for $e_{w}<e<e_{c}$ (strong chaos) up to relatively high values of N.

The transition from weak to strong chaos can be explained from the equilibrium properties of the system. The equilibrium spatial distribution obtained by integrating $f_{eq}$ in Eq. (4) over the momentum, is given by

[TABLE]

with $C_{N}$ a normalization constant. In Eq. (29) the maximum of $\rho_{eq}(\theta)$ occurs at $\theta=0$ by a choice of the origin for the angles. The values of the spatial distribution at $\theta=\pi$ as a function of energy are shown in Fig. 11.

As already pointed out by Manos and Ruffo manos , the transition from weak to strong chaos occurs at the energy value $e_{w}$ when $\rho_{eq}(\pi)$ attains a significant value and particles start to cross at the border $\theta=\pi$, causing a time variation of the phase of the magnetization due to asymmetries in the fluctuations of the distribution in Eq. (29) which is valid for $N\rightarrow\infty$. Indeed, the equations of motion in Eq. (3) for any particle in the system can be written as the equation of a pendulum:

[TABLE]

where $M=\sqrt{M_{x}^{2}+M_{y}^{2}}$ and $\phi=\arctan{M_{y}/M_{x}}$. If the phase $\phi$ is time independent the solutions of Eq. (30) are non-chaotic, while having a positive Lyapunov exponent for a time varying phase, which occurs significantly in the strong chaos energy interval. This point is explored in more detail in Ref. rochamarcos .

As a last result, we investigate the possible critical behavior of the LLE for energies close to $e_{c}$ from below as the theoretical prediction in Eq. (28), by numerically determining $\lambda_{1}$. Although theoretically predicted no numerical verification has been obtained previous to the present work. The results are shown in Fig. 12 for $N=100\>000$ and $N=1000\>000$ and some energy values. The fitting of Eq. (28) is very good for both values of $N$, with an exponent close to the theoretical value $1/6$. Small deviations are possibly due to important correlations in the fluctuations of $\kappa_{R}=M^{2}$, which become more important close to the phase transition, as discussed above. The difference with respect to the exponent $\xi\approx 1/2$ obtained in Ref. ginelli can be explained by our longer simulation times and higher values of $N$ which were made feasible by a massively parallel implementation of our numeric code.

IV Concluding Remarks

This paper addressed the study of chaoticity in the HMF model from the determination of LLE. This paradigmatic model has been widely used in the literature to understand the behavior and some properties of long range interacting systems. Our numerical implementation CUDA allowed to investigate the LLE for a wide range of energies, and values of $N$ as large as $2\times 10^{7}$. The size of the system has been essential to describe the main characteristics features of the exponents. For the homogeneous phase ($e\geq 0.75$,) at all energies, it was shown clearly that the exponents scales with the system size as $N^{\beta}$ with $\beta$ approaching $-1/3$, the theoretical predicted value. Close to the phase transition we must go to higher values of $N$ in order to observe the expected scaling. This comes from non negligible self-correlations in time of the fluctuations of the scalar curvature used in the geometric approach for the theoretical determination of the LLE. For energies below the transition energy we observe two different scaling for the LLE: for energies below $e_{w}\approx 0.15$ the LLE scales approximately with $1/N^{1/3}$, while for $e_{w}<e<e_{c}$ the exponent of the scaling is much smaller than $1/3$. This is explained first by non-negligible correlations in the fluctuations of the curvature of the underlying Riemannian manifold and second by the coupling of the motion of individual particles to a time varying phase of the magnetization.

We also confirmed numerically the existence of a critical exponent associated to the Lyapunov exponent as defined in Eq. (28). The value we have obtained for this exponent is $\xi\approx 0.138$ with is reasonably close to the predicted theoretical value of $1/6$, and far from the value of $1/2$ obtained in Ref. ginelli . With respect to the former, this difference is explainable by the fact that the stochastic process representing the Riemannian curvature on the underlying manifold in not delta correlated, as shown in Fig. 2d. Our parallel implementation of the algorithm for computing the LLE allowed a significant improvement in the accuracy of the numerical results, which possibly explains the variance with the result in ginelli .

Whether such a critical exponent also occurs for other long range interacting systems is an open question that requires to be investigated. A similar but much computationally demanding study for the self-gravitating ring model sota is the subject of ongoing work.

Finally we close this section by pointing out that the theoretical results of Firpo firpo2 , although based on some necessary simplifying assumptions with respect to the geometrical approach of Pettini and collaborators yields results quite often close to our numerical findings. The discrepancies are then explained when those assumptions are not valid, as for instance when the fluctuations of the curvature are non-negligible.

V Acknowledgments

MAA and TMRF would like to thank M.-C. Firpo for fruitful discussions. TMRF was partially financed by CNPq (Brazil) and LAMF was financed by CAPES (Brazil).