Scaling of the dynamics of a homogeneous one-dimensional anisotropic classical Heisenberg model with long-range interactions
SScaling of the dynamics of a homogeneous one-dimensionalanisotropic classical Heisenberg model with long-rangeinteractions
C. R. Louren¸co
Instituto de F´ısica, Universidade de Bras´ılia, CP 04455,70919-970 - Bras´ılia, Brazil and Instituto Federal Bras´ılia, Rodovia DF-128,km 21, Zona Rural de Planaltina, 73380-900 - Planaltina, Brazil
T. M. Rocha Filho ∗ Instituto de F´ısica and International Center for Condensed Matter Physics,Universidade de Bras´ılia, CP 04455, 70919-970 - Bras´ılia, Brazil
Abstract
The dynamics of quasi-stationary states of long-range interacting systems with N particles canbe described by kinetic equations such as the Balescu-Lenard and Landau equations. In the case ofone-dimensional homogeneous systems, two-body contributions vanish as two-body collisions in onedimension only exchange momentum and thus cannot change the one-particle distribution. Using aKac factor in the interparticle potential implies a scaling of the dynamics proportional to N δ with δ = 1 except for one-dimensional homogeneous systems. For the latter different values for δ werereported for a few models. Recently it was shown by Rocha Filho and collaborators [Phys. Rev. E , 032133 (2014)] for the Hamiltonian mean-field model that δ = 2 provided that N is sufficientlylarge, while small N effects lead to δ ≈ .
7. More recently Gupta and Mukamel [J. Stat. Mech.P03015 (2011)] introduced a classical spin model with an anisotropic interaction with a scaling inthe dynamics proportional to N . for a homogeneous state. We show here that this model reducesto a one-dimensional Hamiltonian system and that the scaling of the dynamics approaches N with increasing N . We also explain from theoretical consideration why usual kinetic theory failsfor small N values, which ultimately is the origin of non-integer exponents in the scaling. ∗ marciano@fis.unb.br a r X i v : . [ c ond - m a t . s t a t - m ec h ] J un . INTRODUCTION Systems with long range interactions may present unusual properties such as as non-ergodicity, anomalous diffusion, aging, non-Gaussian Quasi-Stationary States (QSS), neg-ative microcanonical heat capacity, ensemble inequivalence, and more importantly for thepresent work very long relaxation time to thermodynamic equilibrium of a QSS, divergingwith the number of particles N [1–9]. A pair interaction potential is long-ranged in d spatialdimensions if it decays at long distances as r − α with α ≤ d . The dynamics of such systemscan be decomposed in three stages: a violent relaxation into a QSS in a short time, a slowrelaxation of the QSS and the final thermodynamic equilibrium. In some cases after theviolent relaxation the system may also oscillate for a very long or even infinite time arounda QSS [10]. By introducing a Kac factor proportional to 1 /N in the pair-interaction poten-tial the fluid (Vlasov) limit is well defined and given by N → ∞ [11–14]. The dynamics isexactly described by the Vlasov equation for the one particle distribution function, whilefor finite N it is valid only for short times. Collisional terms must be considered for a moreaccurate description of the dynamics for finite N , leading to kinetic equations such as theLandau or Balescu-Lenard equations [13, 15, 16].As already noted, the dynamics of relaxation to equilibrium depends on the number ofparticles in the system and has been extensively studied in the recent literature [1–5, 14, 17–24]. Its dependence on N can be obtained from collisional corrections to the Vlasov equation,i. e. by determining the relevant kinetic equation. Two-body collisions lead to a collisionalintegral in the kinetic equations of order 1 /N , and thus relaxation occurs in a time scaleproportional to N , an exception being three-dimensional gravity with a relaxation timeof order N/ log N [25, 26]. For one-dimensional homogeneous systems two-body terms inthe kinetic equation vanish identically as collisions between two particles results only inmomentum exchange [27–29]. For instance, the Balescu-Lenard and Boltzmann equationsfor a homogeneous one-dimensional system with a pair interaction potential are respectivelywritten as [16]: ∂∂t f ( p ; t ) = 2 π nN ∂∂p (cid:90) d p (cid:90) d k k ˜ V ( k ) | ε ( k, kp ) | δ ( k ( p − p )) (cid:32) ∂∂p − ∂∂p (cid:33) f ( p ; t ) f ( p ; t ) , (1)2nd ∂∂t f ( p ; t ) = 1 N (cid:90) dp | p − p | [ f ( p (cid:48) ; t ) f ( p (cid:48) ; t ) − f ( p ; t ) f ( p ; t )] , (2)where p is the one-dimensional momentum variable, n the particle density, ε ( k, kp ) thedielectric function and p (cid:48) and p (cid:48) the post-collisional momenta for incoming particles withmomenta p and p . Setting ε ( k, kp ) = 1 is equivalent to neglect collective effects and yieldsthe Landau equation. In both cases the right hand is identically zero due to the Dirac deltafunction in the collisional integral in Eq. (1), while for the Boltzmann equation in Eq. (2)we have p (cid:48) = p and p (cid:48) = p . In both cases we obtain ∂f /∂t = 0 if only two-body collisionsare considered, and the dominant term comes then from three-body or higher order terms.This has been considered recently for the Hamiltonian Mean Field (HMF) model, resultingin a dynamics of the homogeneous states scaling with N [30], at variance with previousresults with scalings proportional to N . and exp( N ) which are due to small N effects [31–33]. The N . was also reported for a classical anisotropic Heisenberg model by Gupta andMukamel in Ref. [21] and the question remains if it is due also to small N effects. In thispaper we investigate this issue for small and large N . We re-obtain a N scaling for large N as predicted from kinetic theory, while non-integer exponents in the scaling are due tofinite N effects, as a result of the failure of basic approximations usually considered for thedetermination of kinetic equation in closed form, as shown below.The structure of the paper is as follows: in Section II we present the model and discusssome of its properties. The scaling of the dynamics of a QSS is determine numerically inSection III. We address the limits of applicability of kinetic theory in Section III and closethe paper with some concluding remarks in Section V. II. THE MODEL
The mean-field classical anisotropic Heisenberg model consists of N classical spins (cid:126)S i =( S ix , S ij , S iz ), i = 1 , , ....N , of unit length globally coupled, with Hamiltonian [21]: H = − J N N (cid:88) i,j =1 (cid:126)S i · (cid:126)S j + D N (cid:88) i =1 S iz . (3)The first term in the right-hand side with J > = 1 and D = 15. Note that the coefficient 1 /N in the coupling term in Hamiltonian is aKac factor that makes the energy extensive. The magnetization of the system is defined by: (cid:126)m = (cid:104) (cid:126)S (cid:105) = 1 N N (cid:88) i =1 (cid:126)S i . (4)Using spherical coordinates the spin components are written as S ix = sin( θ i )cos( φ i ), S iy =sin( θ i )sin( φ i ) and S iz = cos( θ i ), and the equations of motion are given by: d (cid:126)S i dt = { (cid:126)S i , H } , (5)with i = 1 , , ...N and the Poisson bracket: { A, B } = N (cid:88) i =1 (cid:40) ∂A∂φ i ∂B∂S iz , ∂A∂S iz ∂B∂φ i (cid:41) . (6)Thus ˙ S ix = S iy m z − S iz m y − DS iy S iz , ˙ S iy = S iz m x − S ix m z + 2 DS ix S iz , ˙ S iz = S ix m y − S iy m x . (7)These equations of motion admit as first integrals the z -component m z of (cid:126)m , the total energyand the the length of each spin. This allows us to rewrite the equations of motion as˙ θ i = m x sin( φ i ) − m y cos( φ i ) , ˙ φ i = m x cot( θ i ) cos( φ i ) + m y cot( θ i ) sin( φ i ) − m z + 2 D cos( θ i ) . (8)As a first result we observe that these equations are canonical and derive from the Hamil-tonian: H = − N (cid:88) i =1 (cid:20) m x cos( φ i ) (cid:113) − S iz + m y sin( φ i ) (cid:113) − S iz + m z S iz − DS iz (cid:105) . (9)were p i ≡ cos θ i = S iz and q i ≡ φ i are canonically conjugate and correspond to the mo-mentum and position variables respectively. As a consequence the model is effectively one-dimensional and thus explains why a scaling proportional to N δ with δ (cid:54) = 1 is observed. Asthe model is effectively one-dimensional and Hamiltonian the tools of kinetic theory can beused to derive a kinetic equation, as described for instance in Ref. [13]. The first consequenceof this fact is that for a homogeneous state in φ , the collisional integral proportion to 1 /N of the Balescu-Lenard equation vanishes, and one must go to the next order in an expansionin powers of 1 /N (see Ref. [30] and references therein).4 II. SCALING OF THE DYNAMICS
In order to study the dynamics of a homogeneous state, and for comparison purposes, weuse here the same initial condition as in Ref. [21], a waterbag state (uniform distribution)in the intervals φ ∈ [0 , π ) and θ ∈ [ π/ − a, π/ a ], with energy per particle: e = D a, (10)and a chosen such that e = 0 .
24. The state is spatially homogeneous and stable for thisenergy. From Ref. [30] the expected scaling of the dynamics of this QSS is N . On the otherhand Gupta and Mukamel obtained from numerical simulations a different scaling in N . .We argue that, similar to what occurs in the HMF model, the N . scaling only occurs forsufficiently small number of particles, while for larger N the scaling becomes proportionalto N .In a homogeneous stable state the spatial distribution for variable φ is always uniformup to small fluctuations, but the distribution for variable θ slowly varies with time towardsthermodynamic equilibrium [30]. As a consequence, the dynamics can be probed by thestatistical moments M n = (cid:104) ( θ − (cid:104) θ (cid:105) ) n (cid:105) . Odd moments of θ vanish for an even distributionin θ as is the case here. Figure 1 shows the second moment M as a function of time.It varies very slowly for the states considered here (it is almost a constant of motion) sowe consider the time evolution of the fourth moment M which is more responsive to smallchanges in the statistical state of the system. In Ref. [21] Gupta and Mukamel considered theaverage (cid:104) cos θ (cid:105) which is more difficult to characterize the small variations in the distributionfunction of θ (compare for instance Fig. 3 of their paper to our Figures 1 and 2 below).The equations of motion in Eq. (8) are solved using a parallel implementation of a fourthorder Runge-Kutta algorithm in a graphics processing unit using the CUDA extension tothe C language [34, 35]. This allows us to perform simulations with a much greater numberof particles than considered in Ref. [21]. The time-step used is δt = 0 .
01 and ensures amaximum relative error in the energy or order 10 − . Figure 2 shows the time evolution of M for different number of particles up to N = 100 000. Figures 3 and 4 show the sameresults but with 1 /N . and 1 /N time rescalings, respectively. A better data collapse isobtained for the N scaling.In order to compare quantitatively ours with previous results, we performed a series ofsimulations with the same number of particles as in Ref. [21] but also considering values5 M × × × × FIG. 1. Second statistical moment M of variable θ for N = 100 000. M N=10 X -6 × X -6 X -6 X -7 X -7 t FIG. 2. (Color online) Moment (cid:104) M (cid:105) of variable θ i as a function of time for different numbers ofparticles N = 10 000 ,
20 000 ,
40 000 ,
60 000 ,
80 000 ,
100 000. of N up to 60 000. By averaging over many realizations we compare the time evolution of M for a given value of N with the previous smaller number of N in the simulations, andperform a least squares fit for the difference between both time series rescaled by 1 /N δ . Theresults are shown in Table I and corroborate, up to some small deviations, a scaling in N .Figure 5 shows the statistical moment M for the same number of particles as in Table Iwith time scaled as 1 /N with a very good data collapse for N ≥ N = 300 , , , N=10 000N=20 000N=40 000
N=60 000
N=80 000N=100 000 x -6 x -6 x -6 x -7 x x x FIG. 3. (Color online) Same as in Fig. 2 but with a time rescaled by ( N −
10 000) − , . N=10 000N=20 000
N=40 000
N=60 000
N=80 000
N=100 000 M x -6 x -6 x -6 x -7 x x x FIG. 4. (Color online) Same as in Fig. 2 but with a time rescaled by ( N −
10 000) − . N = 3000 and 5000 comes from the fact that considering the magnetization as a relevantvariable yields more imprecise results then when considering the statistical moments of themomenta variables (see also the discussion of this in Ref. [35]). IV. LIMITATIONS OF KINETIC THEORY
Our results are in agreement with what is expected from a kinetic theory derived fromthe BBGKY hierarchy in a series expansion in power of 1 /N . Two and three particle cor-relation functions contribute with terms proportional to 1 /N and 1 /N , respectively. As7 N δ
300 1000 1.7671000 3000 1.7973000 5000 2.0155000 10 000 2.05610 000 20 000 2.07220 000 40 000 2.06640 000 60 000 2.096TABLE I. Best scaling in N δ for the moment M between a pair of simulated data with N and N particles. N=300N=1000N=3000N=5000N=10 000N=20 000N=40 000N=60 000 t/N M x -6 x -6 -6 x x x FIG. 5. (Color online) Moment M of the θ i as a function of time for N =300 , , , ,
10 000 ,
20 000 ,
40 0000 ,
60 000 with a time rescaling ( N − − . Thenumber of realization varies from 300 for N = 300 to 25 for N = 60 000. two-particle contributions to the kinetic equation vanishing in the present case one must re-tain the contributions from three-particle collisions which are proportional to 1 /N . Theseconsiderations are based on the introduction of the Kac factor in the Hamiltonian and thescaling proportional to N − . reported by Gupta and Mukamel is re-obtained here for smallervalues of N . This unusual scaling stems on the failure for small N of the Markovianizationhypothesis used in the determination of the Balescu-Lenard and Landau equations, which8equires that the force auto-correlation function (for homogeneous systems) differs signif-icantly from zero only for very short times if compared to the dynamical time scale overwhich the one-particle distribution function varies significantly. Let us show this explicitlyfor the simpler case of the Landau equation, i. e. for weak coupling, as the same kind ofapproximations are used in the deduction of the Balescu-Lenard equation (see Ref. [13] fora thorough discussion on these assumptions).The N -particle distribution function f N ( r , v , . . . , r N , v N ; t ) is the probability density inthe N -particle phase space for a particle at time t to have position r i and momentum p i .Defining the s -particle distribution function by f s ≡ f s ( r , v , . . . , r s , v s ; t ) = (cid:90) d r s +1 d v s +1 · · · d r N d v N f N ( r , v , . . . , r N , v N ; t ) . (11)where r i and p i are the position and momentum vectors of particle i in d dimensions.Liouville equation then implies that the reduced distribution functions satisfy the BBGKYhierarchy [13, 16]: ∂∂t f s = s (cid:88) j =1 ˆ K j f s + 1 N s (cid:88) j 3) stands for permutationsof particles 1, 2 and 3 and C s is the s -particle correlation function. Let us consider aparameter λ (cid:28) V = O ( λ ). A two-particle correlation requires the interaction of two particles to be created and therefore C is9f order λ . A three-particle correlation requires the interaction between two pairs of particlesand thus C is of order λ , and so on. By considering the case s = 1 in Eq. (12) and usingEq. (14) we have: ∂∂t f ( v ; t ) = N − N (cid:90) d v d r ˆΘ [ f ( v ; t ) f ( v ; t ) + C ( v , v , r − r ; t )] . (16)The two-particle correlation function is the solution of the equation obtained by replacingEq. (15) into Eq. (12) for s = 2 and discarding higher order terms containing three-particlecorrelations: (cid:32) ∂∂t − ˆ K − ˆ K (cid:33) C ( v , v , r − r ; t ) = ˆΘ f ( v ; t ) f ( v ; t ) . (17)Its solution can be written as: C ( v , v , r − r ; t ) = e ( ˆ K + ˆ K ) t C ( v , v , r − r ; 0)+ (cid:90) t dt e ( ˆ K + ˆ K ) τ ˆΘ f ( v ; t − τ ) f ( v ; t − τ ) . (18)The first term in the right-hand side of Eq. (18) is a transient term due to correlation at t = 0 and dies out rapidly [13]. By replacing Eq. (18) into Eq. (16) and noting that themean-field force vanishes in a homogeneous state, we obtain (using N − s → N for large N ): ∂∂t f ( v ; t ) = (cid:90) t idτ (cid:90) d v d r ˆΘ e ( ˆ K + ˆ K ) τ ˆΘ f ( v ; t − τ ) f ( v ; t − τ )= (cid:90) t dτ (cid:90) d v d r ∂ ∇ V ( r ) e ( ˆ K + ˆ K ) τ ∇ V ( r ) ∂ f ( v ; t − τ ) f ( v ; t − τ )= (cid:90) t dτ (cid:90) d v d r ∂ ∇ V ( r ) ∇ V ( r − v τ ) ∂ f ( v ; t − τ ) f ( v ; t − τ )(19)with r ≡ r − r and v ≡ v − v . The force auto-correlation of the F ( r , t ) at position r is defined by C ( t ) ≡ (cid:104) F ( t ) F (0) (cid:105) = (cid:90) d r F ( r , F ( r , t ) = (cid:90) d r ∇ V ( r − v t ) ∇ V ( r ) . (20)Thence we have: ∂∂t f ( v ; t ) = (cid:90) t dτ (cid:90) d v ∂ (cid:104) F ( τ ) F (0) (cid:105) ∂ f ( v ; t − τ ) f ( v ; t − τ ) . (21)This is a master equation which is non-Markovian as it depends on f at previous timesform 0 to t . To obtain a true (Markovian) kinetic equation the usual procedure is to assumethat the dynamic time scale t d over which the one-particle distribution function f varies10ignificantly is much greater than the time scale t c such that the force auto-correlation issufficiently small. In this case, one can replace f ( v ; t − τ ) in the integrand in Eq. (21) by f ( v ; t ), which corresponds to the ballistic approximation (free motion for a homogeneoussystem), and extend the time integration to infinity. We then finally obtain the Landauequation: ∂∂t f ( v ; t ) = (cid:90) ∞ dτ (cid:90) d v ∂ (cid:104) F ( τ ) F (0) (cid:105) ∂ f ( v ; t ) f ( v ; t ) . (22)This form will suffice for the present discussion. The same type of considerations are alsonecessary in the determination of the Balescu-Lenard equation [13, 15]. As discussed above,for a one-dimensional homogeneous system these corrections vanish and one must go oneorder further in the 1 /N expansion. Usually one always considers a Markovianization pro-cedure taking into account the time scales such that t d (cid:29) t c . A failure of this conditionimplies, among other things, that the collisional integral does not vanish exactly for a one-dimensional homogeneous system, and one should expect that the scaling of the dynamicsis therefore affected.In order do address this point, we compute the force auto-correlation function fromnumeric simulations by: C ( t ) = 1 N N (cid:88) i =1 F i ( t ) F i (0) , (23)where F i ( t ) is the force on particle i at time t due to all other particles. Figure 6a shows C ( t ) for different number of particles for the present model and the time evolution of (cid:104) M (cid:105) for variable θ . We observe that the time required for a significant decrease of C ( t ), i. e. thecorrelation time t c , is roughly the same for all values of N , while the dynamical time t d issmaller for smaller N as shown in the right panel of Fig. 6. In this way, the correlationtime can become of the same order of magnitude as the dynamical time, breaking downthe Markovian condition, and therefore the usual derivation of Kinetic equations from theBBGKY hierarchy is no longer valid. Figure 6b show the fourth moment M of variable θ and it becomes evident that Markovianity is not valid for N = 1000 and N = 3000, whileit is approximately valid for N = 5000. For N ≥ 10 000 the system is clearly Markovian,in agreement with the results in Table I. This explains why a different scaling in N δ of thedynamics with δ (cid:54) = 2 is observed for homogeneous one-dimensional systems for small N [30].11 .. 1 100 10000t-20-1001020 N=1000N=10 000N=100 000N=1000 000 M N=1000N=3000N=5000N=7000N=10 000N=100 000N=1000 000 (a) (b) x -7 x -7 x -7 x -7 FIG. 6. (Color online) a) Force auto-correlation C ( t ) as a function of time for different values of N . b) Time evolution for the fourth moment (cid:104) M (cid:105) of variable θ averaged over 1000 realizationsexcept for N = 100 000 and N = 1 000 000 with 300 and 200 realizations, respectively. The initialconditions are the same homogeneous state as in Fig. 2 thermalized up to t = 20 . V. CONCLUDING REMARKS We shown in this paper that the mean-field anisotropic Heisenberg model introducedby Gupta and Mukamel in Ref. [21] is effectively a one-dimensional classical Hamiltoniansystem, and the dynamics of a QSS scales as N for large N while the scaling in N . previously reported is due to small N non-Markovian effects in the dynamics. For large N a kinetic equation for a homogeneous one-dimensional long-range interacting system mustconsider three-body collisions, which are of order 1 /N . This approach is only valid if N is sufficiently large such that the contribution of strictly two-body collisions does vanish,while for small N the arguments leading to the N scaling fail. The small N case can betackled using an approach developed by Ettoumi and Firpo by determining the diffusioncoefficient in terms of action variables who used a mean passage time approach [36] andobtained a N . scaling for the Hamiltonian mean Field model [37]. Based on time evolutionof the auto-correlation of the force for the homogeneous case, one can consider if a similarbehavior occurs for non-homogeneous one-dimensional and for higher dimensional systems,and whether and how it influences the scaling for small N . This will be the subject ofa separate publication. This also raises the question whether similar effects might have arole in astrophysics. Indeed smaller globular clusters can have a number of stars of the12rder of just a few thousand, as opposed to 10 stars in a typical galaxy. Other long-rangeinteracting systems may also have a similar behavior. More recently, Gupta and Mukamelintroduced a different model of classical spins in a sphere described bu a two-dimensionalHamiltonian [38] and also displaying a scaling of the dynamics of a homogeneous QSS in N . . Taking into account the discussion in the present paper and in Ref. [30] this is a strongindication that this model is in fact effectively one-dimensional, which is still to be shown. VI. ACKNOWLEDGMENTS The authors acknowledge partial financing by CAPES (Brazilian government agency).TMRF would like to thank B. Marcos for fruitful discussions. [1] A. Campa, T. Dauxois, D. Fanelli and S. Ruffo, Physics of Long-Range Interacting Systems ,(Oxford Univ. Press, Oxford, 2014).[2] Dynamics and Thermodynamics of Systems with Long-Range Interactions , T. Dauxois,S. Ruffo, E. Arimondo and M. Wilkens Eds. (Springer, Berlin, 2002).[3] Dynamics and Thermodynamics of Systems with Long-Range Interactions: Theory and Ex-periments , A. Campa, A. Giansanti, G. Morigi and F. S. Labini (Eds.), AIP Conf. ProceedingsVol. 970 (2008).[4] Long-Range Interacting Systems, Les Houches 2008, Session XC , T. Dauxois, S. Ruffo andL. F. Cugliandolo Eds. (Oxford Univ. Press, Oxford, 2010).[5] A. Campa, T. Dauxois and S. Ruffo, Phys. Rep. , 57 (2009).[6] T. M. Rocha Filho, A. Figueiredo and M. A. Amato, Phys. Rev. Lett. , 190601 (2005).[7] A. Figueiredo, T. M. Rocha Filho and M. A. Amato, Europhys. Lett. , 30011 (2008).[8] T. M. Rocha Filho, M. A. Amato, and A. Figueiredo, Phys. Rev. E , 062103 (2012).[9] T. M. Rocha Filho, M. A. Amato, B. A. Mello, and A. Figueiredo, Phys. Rev. E , 041121(2011).[10] H. Koyama, T. Konishi and S. Ruffo, Comm. Nonlinear Sc. Num. Simul. , 868 (2008).T. Dauxois, S. Ruffo, E. Arimondo and M. Wilkens Eds. (Springer, Berlin, 2002).[11] M. Kac, G. Uhlenbeck and P. Hemmer, J. Math. Phys , 216 (1963). 12] W. Braun and K. Hepp, Commun. Math. Phys. , 125 (1977).[13] R. Balescu, Statistical Dynamics - Matter out of Equilibrium , Imperial College Press (London,1997).[14] T. M. Rocha Filho, M. A. Amato, A. E. Santana, A. Figueiredo and J. R. Steiner, Phys. Rev.E , 032116 (2014).[15] A. Lenard, Ann. Phys. , 390 (1960).[16] R. L. Liboff, Kinetic Theory - Classical, Quantum, and Relativistic Descriptions , 3rd ed,Springer-Verlag (New York, 2003).[17] F. Baldovin and E. Orlandini, Phys. Rev. Lett. , 100601 (2006).[18] P.-H. Chavanis, Physica A , 102 (2006).[19] P. de Buyl, D. Fanelli, R. Bachelard and G. De Ninno, Phys. Rev. Sep. Top. Acc. Beams ,060704 (2009).[20] A. Gabrielli, M. Joyce and B. Marcos, Phys. Rev. Lett. , 210602 (2010).[21] S. Gupta and D. Mukamel, J. Stat. Mech. P03015 (2011).[22] M. Assllani, D. Fanelli, A. Turchi, T. Carletti and X. Leoncini, Phys. Rev. E , 021148(2012).[23] R. Pakter and Y. Levin, Phys. Rev. Lett. , 140601 (2013).[24] R. Bachelard, C. Chandre, D. Fanelli, X. Leoncini and S. Ruffo, Phys. Rev. Lett. , 260603(2008).[25] S. Chandrasekhar, Astrophys. J. , 47 (1944).[26] P.-H. Chavanis, J. Stat. Mech. P05019 (2010).[27] O. C. Eldridge and M. Feix, Phys. Fluids , 398 (1963).[28] B. B. Kadomtsev and O. P. Pogutse, Phys. Rev. Lett. , 1155 (1970).[29] M. M. Sano, J. Phys. Soc. Japan , 024008 (2012).[30] T. M. Rocha Filho, A. E. Santana, M. A. Amato, and A. Figueiredo, Phys. Rev. E , 032133(2014).[31] D. H. Zanette and M. A. Montemurro, Phys. Rev. E , 031105 (2003).[32] Y. Y. Yamaguchi, J. Barr´e, F. Bouchet, T. Dauxois and S. Ruffo, Physica A , 36 (2004).[33] A. Campa, A. Giansanti and G. Morelli, Phys. Rev. E , 041117 (2007).[34] NVIDIA, CUDA Programming Guide, Ver. 4.0, 2011.[35] T. M. Rocha Filho, Comp. Phys. Comm. , 1364 (2014). 36] W. Ettoumi and M.- C. Firpo, Phys. Rev. E , 030102(R) (2013).[37] M. Antoni and S. Ruffo, Phys. Rev. E , 2361 (1995).[38] S. Gupta and D. Mukamel, Phys. Rev. E , 052137 (2013)., 052137 (2013).