Classical Heisenberg spins with long-range interactions: Relaxation to equilibrium for finite systems
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b Classical Heisenberg spins with long-range interactions:Relaxation to equilibrium for finite systems
Julien Barr´e and Shamik Gupta Laboratoire J. A. Dieudonn´e, Universit´e de Nice-Sophia Antipolis, UMR CNRS 7351,Parc Valrose, F-06108 Nice Cedex 02, France Laboratoire de Physique Th´eorique et Mod`eles Statistiques, UMR CNRS 8626, Universit´eParis-Sud, Orsay, FranceE-mail: [email protected],[email protected]
Abstract.
Systems with long-range interactions often relax towards statistical equilibriumover timescales that diverge with N , the number of particles. A recent work [S. Gupta and D.Mukamel, J. Stat. Mech.: Theory Exp. P03015 (2011)] analyzed a model system comprising N globally coupled classical Heisenberg spins and evolving under classical spin dynamics. Itwas numerically shown to relax to equilibrium over a time that scales superlinearly with N .Here, we present a detailed study of the Lenard-Balescu operator that accounts at leadingorder for the finite- N effects driving this relaxation. We demonstrate that corrections atthis order are identically zero, so that relaxation occurs over a time longer than of order N ,in agreement with the reported numerical results.PACS numbers: 05.20.Dd, 45.50.-j, 52.25.Dg Keywords: Kinetic theory of gases and liquids, Metastable states ontents1 Introduction 22 The model 43 The Klimontovich equation 54 The Vlasov equation, and a class of stationary solutions 65 The Lenard-Balescu equation 7 D δg φ ∂δf∂φ E . . . . . . . . . . . . . . . . . . . . . . . . . . . 115.3.2 Computing D δg u ∂δf∂u E . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Long-range interacting systems are characterized by an interparticle potential with a rangethat is of the order of the system size. In d dimensions, this corresponds to potentialsdecaying at large separation, r , as 1 /r α , where α lies in the range 0 ≤ α ≤ d [1, 2]. Examplesinclude gravitational systems [3], plasmas [4], two-dimensional hydrodynamics [5], chargedand dipolar systems [6], and many others.Despite obvious differences, these systems often share a common phenomenology (see[1,2]). In particular, the relaxation to equilibrium of an isolated long-range interacting systemof N particles proceeds in two steps: first, a collisionless relaxation, described by a Vlasov-type equation, brings the system close to a nonequilibrium state, called the ”quasistationarystate” (QSS), whose lifetime increases with N ; second, on timescales diverging with N ,discreteness effects due to finite value of N drive the system towards Boltzmann-Gibbsequilibrium. This second step is usually described by a Lenard-Balescu-type equation [7, 8].This scenario is well established for plasmas and self-gravitating systems [9], and has beenstudied in detail in various toy models for long-range interactions [10–14]. The QSS lifetime,which may be regarded to be of the same order of magnitude as the relaxation time, is thusan important quantity: knowing it allows to distinguish between non-relaxed systems, which2hould be described by a QSS, and relaxed ones, for which collisional effects need to betaken into account and an equilibrium description may be relevant. This lifetime dependson the system under consideration. Kinetic theory usually predicts a QSS lifetime of order N , e.g., for 3d plasmas or 1d self-gravitating systems (see [15] for recent numerical tests),but this is not always the case; for example, the QSS lifetime is of order N/ ln N for 3d self-gravitating systems [16]. Furthermore, it is known that the Lenard-Balescu collision termvanishes for 1d systems which do not develop any spatial inhomogeneity (see [17] for the 1dCoulomb case); one thus expects a relaxation time much longer than N in these cases. Thisis indeed numerically observed, see [12], where a time of order N . is reported, and [18],where larger systems sizes are studied and the relaxation time is claimed to be of order N .A similar vanishing of the Lenard-Balescu operator has been found for point vortices in anaxisymmetric configuration [19–21].In a recent work, QSSs have been looked for and found in a dynamical setting differentfrom the ones reviewed above, namely, in an anisotropic Heisenberg model with mean-fieldinteractions [22]. Note that similar spin models with mean field interactions have beensuggested to be relevant to describe some layered spin structures [23]. Specifically, themodel in [22] comprises N globally coupled three-component Heisenberg spins evolving underclassical spin dynamics. An associated Vlasov-type equation is introduced, and QSSs arestationary solutions of this equation. In addition, numerical simulations for axisymmetricQSS suggest a QSS lifetime increasing superlinearly with N . In order to understand thisobservation analytically, we present in this work a detailed study of the Lenard-Balescuoperator that accounts for leading finite- N corrections of order 1 /N to the Vlasov equation.With respect to 1D Hamiltonian systems, the spin dynamics introduces a new term inthe Lenard-Balescu operator; it also complicates the analytical structure of the dispersionrelation for the Vlasov-type equation. Nevertheless, we can still demonstrate that correctionsat order 1 /N are identically zero, so that relaxation occurs over a time longer than of order N , in agreement with the reported numerical results.The paper is structured as follows. In section 2, we describe the model of study and itsequilibrium phase diagram. As a step towards deriving the Vlasov equation to analyze theevolution of the phase space distribution in the limit N → ∞ , in section 3, we first writedown the so-called Klimontovich equation; this leads to a derivation of the Vlasov equationand a discussion of a class of its stationary solutions in section 4. Section 5 contains ourmain results: it is devoted to the derivation of the Lenard-Balescu equation for our modelof study, that is, the leading 1 /N correction to the Vlasov equation; we show that thiscorrection identically vanishes. In section 6, we consider an example of a Vlasov-stationarysolution. In the energy range in which it is Vlasov stable, we demonstrate by performingnumerical simulations of the dynamics that indeed its relaxation to equilibrium occurs overa timescale that does not grow linearly but rather superlinearly with N , in support of ouranalysis. The paper ends with conclusions. 3 . The model The model studied in Ref. [22] comprises N globally coupled classical Heisenberg spins ofunit length, denoted by S i = ( S ix , S iy , S iz ), i = 1 , , . . . , N . In terms of spherical polar angles θ i ∈ [0 , π ] and φ i ∈ [0 , π ], one has S ix = sin θ i cos φ i , S iy = sin θ i sin φ i , S iz = cos θ i . TheHamiltonian of the system is H = − J N N X i,j =1 S i · S j + D N X i =1 S iz . (1)Here, the first term with J >
D >
0, forwhich the energy is lowered by having the magnetization m ≡ N N X i =1 S i (2)pointing in the xy plane. The coupling constant J in equation (1) is scaled by N to makethe energy extensive [24], but the system is non-additive, implying thereby that it cannotbe trivially subdivided into independent macroscopic parts, as is possible with short-rangesystems. In this work, we take unity for J and the Boltzmann constant.In equilibrium, the system (1) shows a continuous phase transition as a function of theenergy density e , from a low-energy magnetized phase in which the system is ordered in the xy plane to a high-energy non-magnetized phase, across a critical threshold given by [22] e c = D (cid:16) − β c (cid:17) , (3)where the inverse temperature β c satisfies2 β c = 1 − β c D + e − β c D √ πβ c D Erf[ √ β c D ] . (4)Here, Erf[ x ] = (2 / √ π ) R x d t e − t is the error function.The microcanonical dynamics of the system (1) is given by the set of coupled first-orderdifferential equationsd S i d t = { S i , H } ; i = 1 , , . . . , N. (5)Here, noting that the canonical variables for a classical spin are φ and u ≡ cos θ, (6)the Poisson bracket { A, B } for two functions of the spins are given by { A, B } ≡ P Ni =1 ( ∂A/∂φ i ∂B/∂u i − ∂A/∂u i ∂B/∂φ i ), which may be rewritten as [25] { A, B } = N X i =1 S i · ∂A∂ S i × ∂B∂ S i . (7)4sing equations (5) and (7), we obtain the equations of motion of the system as˙ S ix = S iy m z − S iz m y − DS iy S iz , (8)˙ S iy = S iz m x − S ix m z + 2 DS ix S iz , (9)˙ S iz = S ix m y − S iy m x , (10)where the dots denote derivative with respect to time. Summing over i in equation (10), wefind that m z is a constant of motion. The dynamics also conserves the total energy and thelength of each spin. Using equations (8), (9), and (10), we obtain the time evolution of thevariables θ i and φ i as˙ θ i = m x sin φ i − m y cos φ i , (11)˙ φ i = m x cot θ i cos φ i + m y cot θ i sin φ i − m z + 2 D cos θ i . (12)
3. The Klimontovich equation
The state of the N -spin system is described by the discrete one-spin time-dependent densityfunction f d ( u, φ, t ) = 1 N N X i =1 δ ( u − u i ( t )) δ ( φ − φ i ( t )) , (13)which is defined such that f d ( u, φ, t )d u d φ counts the number of spins with its canonicalcoordinates in [ u, u + du ] and [ φ, φ + dφ ]. Here, δ is the Dirac delta function, ( u, φ ) arethe Eulerian coordinates of the phase space, while ( u i , φ i ) are the Lagrangian coordinatesof the spins. Note that f d satisfies f d ( u, φ, t ) = f d ( u, φ + 2 π, t ) and the normalization R π d φ R − d u f d ( u, φ, t ) = 1.Differentiating f d with respect to time and using the equations of motion, (11) and (12),we get the Klimontovich equation for the time evolution of f d as ∂f d ( u, φ, t ) ∂t = − g u ∂∂u f d ( u, φ, t ) − g φ ∂∂φ f d ( u, φ, t ) , (14)where g u ≡ g u [ f d ]( u, φ )= √ − u ( m y [ f d ] cos φ − m x [ f d ] sin φ ) , (15) g φ ≡ g φ [ f d ]( u, φ )= m x [ f d ] u √ − u cos φ + m y [ f d ] u √ − u sin φ − m z [ f d ] + 2 Du, (16)( m x , m y , m z )[ f d ]= Z π d φ Z − d u ( √ − u cos φ, √ − u sin φ, u ) f d ( u, φ, t ) . (17)5 . The Vlasov equation, and a class of stationary solutions We now define an averaged one-spin density function f ( u, φ, t ), corresponding to averaging f d ( u, φ, t ) over an ensemble of initial conditions close to the same macroscopic initial state.We write, quite generally, for an initial condition of the ensemble that f d ( u, φ, t ) = f ( u, φ, t ) + 1 √ N δf ( u, φ, t ) , (18)where, denoting by angular brackets the averaging with respect to the initial ensemble, wehave h f d i = f . Here, δf gives the difference between f d , which depends on the given initialcondition, and f , which depends on the average with respect to the ensemble of initialconditions.Using equation (18) in equation (14), we get ∂f ∂t + g u ∂f ∂u + g φ ∂f ∂φ + 1 N h δg u ∂δf∂u + δg φ ∂δf∂φ i + 1 √ N h ∂δf∂t + δg u ∂f ∂u + g u ∂δf∂u + δg φ ∂f ∂φ + g φ ∂δf∂φ i = 0 , (19)where g u = g u [ f ] , g φ = g φ [ f ], and δg u = √ − u ( m y [ δf ] cos φ − m x [ δf ] sin φ ) , (20) δg φ = m x [ δf ] u √ − u cos φ + m y [ δf ] u √ − u sin φ − m z [ δf ] . (21)We now average equation (19) with respect to the ensemble of initial conditions, and notethat h δf i = 0 implies h m x [ δf ] i = h m y [ δf ] i = h m z [ δf ] i = 0. Thus h δg u i = h δg φ i = 0, and weget ∂f ∂t + g u ∂f ∂u + g φ ∂f ∂φ = − N D δg u ∂δf∂u + δg φ ∂δf∂φ E . (22)For finite times and in the limit N → ∞ (or, for times t ≪ N ), we obtain the Vlasovequation satisfied by the averaged one-spin density function f [26] as ∂f ∂t + g u ∂f ∂u + g φ ∂f ∂φ = 0 . (23)Note that the Vlasov equation has been formally obtained after averaging over anensemble of initial conditions. However, if the fluctuations in the initial conditions areweak and do not grow too fast in time, we expect the Vlasov equation to also describethe time evolution of a single initial condition in the limit N → ∞ . This is put on firmmathematical grounds in [27–29], for systems with a standard kinetic energy and a regularenough interaction potential.From equations (15), (16), and (17), it is clear that any distribution that does notdepend on the angle φ (thus, axisymmetric about the z -axis) is a stationary solution of theVlasov equation (23). 6 . The Lenard-Balescu equation We obtain in this section the main result of our paper: for model (1), the Lenard-Balescuoperator, computed for stationary solutions of the Vlasov equation of the form f ( u ),identically vanishes. The Lenard-Balescu equation describes the slow evolution of a stable stationary solution ofthe Vlasov equation under the influence of finite- N corrections, at leading order in 1 /N [26].From equation (22), we get the Lenard-Balescu equation as ∂f ∂t = − N D δg u ∂δf∂u + δg φ ∂δf∂φ E , (24)where, subtracting equation (23) from equation (19), and keeping only the terms of order1 / √ N , we find that δf follows the Vlasov equation linearized around its stable stationarysolution f : ∂δf∂t = √ − u h δm x sin φ − δm y cos φ i ∂f ∂u − (2 Du − m z [ f ]) ∂δf∂φ . (25)Here, we have used g φ = 2 Du − m z [ f ] and δm x,y ≡ m x,y [ δf ]. Now, a natural timescaleseparation hypothesis greatly reduces the complexity of finding the solutions of the coupledsystem of PDEs, equations (24) and (25): The first of the two equations evolves on a slow O (1 /N ) timescale, while the second one evolves on a fast O (1) timescale. Then, we may firstsolve equation (25), and then use its solution to compute the right-hand side of (24) in thelimit t → ∞ .At this point, it is useful to make a comparison of our case with the standard case ofparticles with a kinetic energy moving in a classical potential, for example, a 1d system ofparticles with Coulomb interactions. The equivalent of the axisymmetric stationary solutions(independent of φ ) introduced in section 4 are the homogeneous solutions which depend onlyon velocity in this standard setting, so that the analog of δg φ vanishes, whereas in our caseof the spin dynamics, we have to deal with the extra term D δg φ ∂δf∂φ E . We now solve equation (25) for δf , using Fourier-Laplace transforms δf ( u, φ, t ) = ∞ X k = −∞ Z Γ d ω π f δf k ( u, ω ) e i ( kφ − ωt ) , (26) f δf k ( u, ω ) = Z π d φ π Z ∞ d t δf ( u, φ, t ) e − i ( kφ − ωt ) , (27)7here the Laplace contour Γ is a horizontal line in the complex- ω plane that passes aboveall singularities of f δf k ( u, ω ).We have δm x = Z Γ d ω π e − iωt Z − d u π √ − u (cid:16)f δf − ( u, ω ) + f δf +1 ( u, ω ) (cid:17) , (28) δm y = Z Γ d ω π e − iωt Z − d u πi √ − u (cid:16)f δf − ( u, ω ) − f δf +1 ( u, ω ) (cid:17) , (29) δm z = m z [ δf ] = Z Γ d ω π e − iωt Z − d u πu f δf ( u, ω ) , (30)so that equation (25) gives f δf ± ( u, ω ) = − π √ − u f ′ ( u )2 Du − m z [ f ] ∓ ω Z − d u ′ √ − u ′ f δf ± ( u ′ , ω ) ∓ iδf ± ( u, Du − m z [ f ] ∓ ω , (31)where f ′ ( u ) = ∂f ( u ) /∂u and δf ± ( u,
0) is the Fourier transform of the initial fluctuations δf ( u, φ, √ − u and then integratingover u , we get ǫ ± ( ω ) Z − d u √ − u f δf ± ( u, ω ) = ∓ i Z − d u √ − u δf ± ( u, Du − m z [ f ] ∓ ω , (32)where ǫ ± ( ω ) is the so-called “Plasma response dielectric function” [26]: ǫ ± ( ω ) = 1 + π Z LC d u (1 − u ) f ′ ( u )2 Du − m z [ f ] ∓ ω . (33)To make the dielectric function ǫ +1 (also, ǫ − ) analytic in the vicinity of the real axis(Im( ω ) = 0), which will be needed for later purpose, the above integral has to be performedalong the Landau contour shown in Fig. 1, as discussed in [26]; we have in this case ǫ ± ( ω ) = π R LC d u (1 − u ) f ′ ( u )2 Du − m z [ f ] ∓ ω ; (Im( ω ) > , π P R LC d u (1 − u ) f ′ ( u )2 Du − m z [ f ] ∓ ω ± i π D (1 − u ) f ′ ( u ) | ( m z [ f ] ± ω ) / (2 D ) ; (Im( ω ) = 0) , π R LC d u (1 − u ) f ′ ( u )2 Du − m z [ f ] ∓ ω ± i π D (1 − u ) f ′ ( u ) | ( m z [ f ] ± ω ) / (2 D ) ; (Im( ω ) < , (34)8here P denotes the principal part. If ω / ∈ [ − D − m z [ f ]; 2 D − m z [ f ]] (respectively, ω / ∈ [ − D + m z [ f ]; 2 D + m z [ f ]]), equation (33) already defines an analytic function ǫ +1 (respectively. ǫ − ) in the vicinity of a real ω , without the need to take into account extrapole contributions as in (34). Note that ǫ ± ( ω ) has two branch cut singularities on the realaxis at ω = 2 D − m z [ f ] and ω = − D − m z [ f ], which is the reason why these functionsmay be seen as multi-valued in the lower-half ω -plane. Im( u ) Re( u )Im( u ) Re( u )Im( u ) Re( u )Im( u )000Im( ω ) < , − D − m z [ f ] < Re( ω ) < D − m z [ f ]Im( ω ) = 0 , − D − m z [ f ] < Re( ω ) < D − m z [ f ] − − − ω ) > , − D − m z [ f ] < Re( ω ) < D − m z [ f ] Figure 1.
The Landau contour LC to evaluate ǫ +1 ( ω ), shown as the dashed red line in thecomplex- u plane. Using equation (32) in equation (31) gives f δf ± ( u, ω ) = ± iπ √ − u f ′ ( u ) ǫ ± ( ω )[2 Du − m z [ f ] ∓ ω ] Z − d u √ − u δf ± ( u, φ, Du − m z [ f ] ∓ ω ∓ iδf ± ( u, Du − m z [ f ] ∓ ω . (35)We see from the above expression that the real pole at ω = ± (2 Du − m z [ f ]) is due to thefree part of the evolution that does not involve interaction among the spins, and results inundamped oscillations of the fluctuations δf ( u, φ, t ), see equation (26). The other set of9oles corresponds to the zeros of the dielectric function ǫ ± ( ω ), i.e., values ω p (complex ingeneral) that satisfy ǫ ± ( ω p ) = 0 . (36)Equation (26) implies that these poles determine the growth or decay of the fluctuations δf ( u, φ, t ) in time, depending on their location in the complex- ω plane. For example, whenthe poles lie in the upper-half complex ω -plane, the fluctuations grow in time. On the otherhand, when the poles are either on or below the real- ω axis, the fluctuations do not grow intime, but rather oscillate or decay in time, respectively. Then, the condition ensuring linearstability of a stationary solution of the Vlasov equation reads ǫ ± ( ω p ) = 0 ⇒ Im( ω p ) ≤ . (37)The condition Im( ω p ) = 0 corresponds to marginal stability.To end this section, let us define for later use the quantities δm ± ≡ δm x ± iδm y . (38)On using equation (32), we get the corresponding Laplace transforms as f δm + ( ω ) = 2 iπǫ − ( ω ) Z − d u √ − u δf − ( u, Du − m z [ f ] + ω , (39) f δm − ( ω ) = − iπǫ +1 ( ω ) Z − d u √ − u δf +1 ( u, Du − m z [ f ] − ω . (40)Equation (35) may now be expressed in terms of f δm ± as f δf +1 ( u, ω ) = − √ − u f ′ ( u )2[2 Du − m z [ f ] − ω ] f δm − ( ω ) − iδf +1 ( u, Du − m z [ f ] − ω , (41) f δf − ( u, ω ) = − √ − u f ′ ( u )2[2 Du − m z [ f ] + ω ] f δm + ( ω ) + iδf − ( u, Du − m z [ f ] + ω . (42) We now compute the Lenard-Balescu operator, given by the right hand side of equation (24),in the limit t → ∞ , by using the results of the preceding subsection. We have D δg φ ∂δf∂φ E = X k,l Z Γ Z Γ ′ d ω d ω ′ π e i ( k + l ) φ e − i ( ω + ω ′ ) t il D e δg φ,k ( u, ω ) f δf l ( u, ω ′ ) E , (43) D δg u ∂δf∂u E = X k,l Z Γ Z Γ ′ d ω d ω ′ π e i ( k + l ) φ e − i ( ω + ω ′ ) t D e δg u,k ( u, ω ) ∂ f δf l ( u, ω ′ ) ∂u E . (44)From equations (20) and (21), we have δg φ = (cid:16) δm + + δm − (cid:17) cos φ u √ − u + (cid:16) δm + − δm − i (cid:17) sin φ u √ − u − δm z [ δf ] , δg u = √ − u h(cid:16) δm + − δm − i (cid:17) cos φ − (cid:16) δm + + δm − (cid:17) sin φ i , (46)so that we have e δg φ, ± ( u, ω ) = f δm ∓ ( ω ) u √ − u , e δg φ, ( u, ω ) = f δm z ( ω ) , (47) e δg u, ± ( u, ω ) = ∓ √ − u i f δm ∓ ( ω ) . (48)From equations (39) and (40), we see that f δm + (respectively, f δm − ) depends on δf − ( u, δf +1 ( u, h δf k ( u, δf l ( v, i for the initial fluctuationsat t = 0. Note that t = 0 as a notation is somewhat inappropriate, since actually thiscomputation has to be repeated for any value of the slow time. We may assume that at t = 0, the spins are almost independent (i.e., the two-spin correlation is of order 1 /N ). Inthis case, one gets h δf k ( u, δf l ( v, i = δ k, − l π [ f ( u ) δ ( u − v ) + h ( u, v )] , (49)where δ k, − l is the Kronecker delta function, and h ( u, v ) is a smooth function. The preciseform of this undetermined smooth function will play no role in the computation, as we willshow below. D δg φ ∂δf∂φ E Combining equations (47) and (49), we see that on the righthand side of equation (43), only the terms ( k = +1 , l = −
1) and ( k = − , l = +1) givea non-zero contribution. We detail below the computation for the ( k = +1 , l = −
1) case,the other being similar. Note that D δg φ ∂δf∂φ E is real; thus, we have to compute only the realpart of the ( k = +1 , l = −
1) term, since its imaginary part must cancel with that of the( k = − , l = +1) term. In the following computation, we set h ( u, v ) = 0; we will check atthe end that indeed the contributions containing h vanish.We have h e δg φ, +1 ( u, ω ) f δf − ( u, ω ′ ) i = − uf ′ ( u )4(2 Du − m z [ f ] + ω ′ ) h f δm − ( ω ) f δm + ( ω ′ ) i + i u √ − u (2 Du − m z [ f ] + ω ′ ) h f δm − ( ω ) δf − ( u, i . (50)We thus need h f δm − ( ω ) f δm + ( ω ′ ) i = 2 πǫ +1 ( ω ) ǫ − ( ω ′ ) Z − d v (1 − v ) f ( v )(2 Dv − m z [ f ] − ω )(2 Dv − m z [ f ] + ω ′ ) , (51) h f δm − ( ω ) δf − ( u, i = − i √ − u f ( u ) ǫ +1 ( ω )(2 Du − m z [ f ] − ω ) , (52)11here we have used equation (49). Using these equations in equation (50), we obtain h e δg φ, +1 ( u, ω ) f δf − ( u, ω ′ ) i = − πuf ′ ( u )2(2 Du − m z [ f ] + ω ′ ) ǫ +1 ( ω ) ǫ − ( ω ′ ) × Z − d v (1 − v ) f ( v )(2 Dv − m z [ f ] − ω )(2 Dv − m z [ f ] + ω ′ )+ uf ( u )2(2 Du − m z [ f ] + ω ′ )(2 Du − m z [ f ] − ω ) ǫ +1 ( ω ) ≡ b + b . (53)Now, to compute the contribution to D δg φ ∂δf∂φ E for ( k = +1 , l = − b and b over ω and ω ′ ; we define B ≡ Z Γ Z Γ ′ d ω d ω ′ π ( − i ) e − i ( ω + ω ′ ) t b , B ≡ Z Γ Z Γ ′ d ω d ω ′ π ( − i ) e − i ( ω + ω ′ ) t b . (54)We want to compute B and B in the limit t → ∞ ; thus, we will discard all terms decayingfor large t .Let us start with B . We have B = i πuf ′ ( u )2 Z Γ ′ d ω ′ π e − iω ′ t (2 Du − m z [ f ] + ω ′ )(2 Dv − m z [ f ] + ω ′ ) ǫ − ( ω ′ ) × Z Γ d ω π e − iωt Z − d v (1 − v ) f ( v ) ǫ +1 ( ω )(2 Dv − m z [ f ] − ω ) . (55) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) ω ) Zeros of ǫ +1 ( ω )2 D − m z [ f ] − iδ − D − m z [ f ] − iδ Re( ω )Γ2 Dv − m z [ f ] − iδ ′ Figure 2.
The Γ contour to evaluate the integral in equation (56). Here, δ → δ ′ → t , Z Γ d ω π e − iωt ǫ +1 ( ω )(2 Dv − m z [ f ] − ω ) ∼ i e − i (2 Dv − m z [ f ]) t ǫ +1 (2 Dv − m z [ f ]) , (56)which is obtained by deforming the Γ contour into the Im( ω ) < e − iωt factor, the only contribution to the contour integral in thelimit t → ∞ comes from the pole 2 Dv − m z [ f ] − iδ ′ . Indeed, since f is assumed stable,the zeros of ǫ +1 , if any, have a negative imaginary part; their contributions thus decayexponentially in time. Similarly, the branch-cut singularities, see Fig. 2, contribute termsdecaying algebraically in time. Note also that the important singularity at 2 Dv − m z [ f ] isreal, and is in the range [ − D − m z [ f ]; 2 D − m z [ f ]]. Then, according to the discussion insection 5.2, one needs to use the expression (34) for ǫ +1 . The same remark applies to all thecomputations below, and we will not recall it each time.The integration over ω ′ in equation (55) is a bit more complicated, since there are twopoles on the real axis, at ω ′ = − Du + m z [ f ] and ω ′ = − Dv + m z [ f ]. With the samemethod as for ω , we obtain Z Γ ′ d ω ′ π e − iω ′ t (2 Du − m z [ f ] + ω ′ )(2 Dv − m z [ f ] + ω ′ ) ǫ − ( ω ′ )= − i (cid:20) e i (2 Du − m z [ f ]) t D ( v − u ) ǫ − ( − Du + m z [ f ]) − e i (2 Dv − m z [ f ]) t D ( v − u ) ǫ − ( − Dv + m z [ f ]) (cid:21) . (57)We finally get B = i πuf ′ ( u )2 Z − d v (1 − v ) f ( v ) ǫ +1 (2 Dv − m z [ f ]) × (cid:26) e − i D ( v − u ) t D ( v − u ) ǫ − ( − Du + m z [ f ]) − D ( v − u ) ǫ − ( − Dv + m z [ f ]) (cid:27) . (58)To perform the integral over v , we will use the following lemma. Lemma: lim t →∞ Z − d v ϕ ( v ) v − u (cid:18) e − i D ( v − u ) t ǫ − ( − Du + m ) − ǫ − ( − Dv + m ) (cid:19) = − iπ ϕ ( u ) ǫ − ( − Du + m ) . (59) Proof: Z − d v ϕ ( v ) v − u (cid:18) e − i D ( v − u ) t ǫ − ( − Du + m ) − ǫ − ( − Dv + m ) (cid:19) = Z D (1 − u ) t − D (1+ u ) t d x ϕ ( u + x Dt ) x (cid:18) e − ix ǫ − ( − Du + m ) − ǫ − ( − Du + m − xt ) (cid:19) . (60)Recalling that − < u <
1, and taking the limit t → ∞ , the above expression simplifies to Z ∞−∞ d x ϕ ( u ) x (cid:18) e − ix ǫ − ( − Du + m ) − ǫ − ( − Du + m ) (cid:19) ϕ ( u ) ǫ − ( − Du + m ) Z ∞−∞ d x e − ix − x . (61)The last integral is − iπ , which completes the proof of the lemma.Using the lemma, we conclude that B = π uf ′ ( u )(1 − u ) f ( u )4 Dǫ +1 (2 Du − m z [ f ]) ǫ − ( − Du + m z [ f ]) . (62)We now turn to the computation of B . The integration over ω and ω ′ is performed asabove, deforming the contours in the lower-half ω ′ -plane, and keeping only the contributionsof the poles on the real axis. We obtain B = − i uf ( u )2 1 ǫ +1 (2 Du − m z [ f ])= − i uf ( u )2 ǫ − ( − Du + m z [ f ]) ǫ +1 (2 Du − m z [ f ]) ǫ − ( − Du + m z [ f ]) . (63)As explained above, we need to compute only the real part of B ; thus, we keep only thecontribution coming from the imaginary part of ǫ − ( − Du + m z [ f ]). Using from equation(34) that Im[ ǫ − ( − Du + m z [ f ])] = − π D (1 − u ) f ′ ( u ) , (64)we conclude thatRe( B ) = − π uf ( u )(1 − u ) f ′ ( u )4 Dǫ +1 (2 Du − m z [ f ]) ǫ − ( − Du + m z [ f ]) . (65)From equations (62) and (65), we find that B + Re( B ) = 0. Similar to above, one can showthat the real part of the contribution to D δg φ ∂δf∂φ E from ( k = − , l = +1) also vanishes, while,as discussed above, the imaginary part of the contribution to D δg φ ∂δf∂φ E from ( k = − , l = +1)must cancel that from ( k = +1 , l = − D δg φ ∂δf∂φ E = 0 . (66)We need now to check that the contributions containing the function h introduced in(49) indeed vanish. For example, let us compute its contribution to B . First, its contributionto h f δm − ( ω ) δf − ( u, i is − iǫ +1 ( ω ) Z − d v √ − v Dv − m z [ f ] − ω h ( v, u ) . (67)Thus, its contribution to B is u √ − u Z − d v √ − v h ( v, u ) Z Γ Z Γ ′ d ω d ω ′ π e − i ( ω + ω ′ ) t ǫ +1 ( ω )(2 Dv − m z [ f ] − ω ) × Du − m z [ f ] + ω ′ ) . (68)14he integrals over ω and ω ′ can be performed as before. Since the poles for ω and ω ′ aredifferent, we see that for large t , a factor e − Di ( v − u ) t oscillating rapidly in time remains inthe integral over v ; this leads to this integral vanishing in the limit t → ∞ . A similarphenomenon ensures that all terms containing the function h vanish in the same way. Thus,we set h = 0 in the following, without modifying the results. D δg u ∂δf∂u E As for D δg φ ∂δf∂φ E , we see that on the right hand side of equation(44), only the terms ( k = +1 , l = −
1) and ( k = − , l = +1) are non-zero. We give belowthe computation for the case ( k = +1 , l = −
1) case, the other being similar. Again, we canrestrict the computations to the real part of each term, since D δg u ∂δf∂u E is real.From equations (42) and (48), we get D e δg u, +1 ( u, ω ) ∂ f δf − ( u, ω ′ ) ∂u E = √ − u i ∂∂u h √ − u f ′ ( u )[2 Du − m z [ f ] + ω ′ ] h f δm − ( ω ) f δm + ( ω ′ ) i i − √ − u ∂∂u h Du − m z [ f ] + ω ′ h f δm − ( ω ) δf − ( u, i i = π i √ − u ǫ +1 ( ω ) ǫ − ( ω ′ ) ∂∂u h √ − u f ′ ( u )[2 Du − m z [ f ] + ω ′ ] × Z − d v (1 − v ) f ( v )(2 Dv − m z [ f ] − ω )(2 Dv − m z [ f ] + ω ′ ) i + i √ − u ǫ +1 ( ω ) ∂∂u h √ − u f ( u )(2 Du − m z [ f ] + ω ′ )(2 Du − m z [ f ] − ω ) i ≡ a + a , (69)where, in obtaining the second equality, we have used equations (51) and (52).Now, to compute the contribution to D δg u ∂δf∂u E for ( k = +1 , l = − a and a over ω and ω ′ . Let us define A ≡ Z Γ Z Γ ′ d ω d ω ′ π e − i ( ω + ω ′ ) t a , A ≡ Z Γ Z Γ ′ d ω d ω ′ π e − i ( ω + ω ′ ) t a . (70)We want to compute A and A in the limit t → ∞ , so that we may discard all termsdecaying for large t .Let us first compute A . Integration over ω and ω ′ may be carried out by deformingthe contours Γ and Γ ′ , as done in the preceding subsection. One gets Z Γ d ω π e − iωt ǫ +1 ( ω )(2 Dv − m z [ f ] − ω ) = i e − i (2 Dv − m z [ f ]) t ǫ +1 (2 Dv − m z [ f ]) , (71)and Z Γ ′ d ω ′ π e − iω ′ t (2 Du − m z [ f ] + ω ′ )(2 Dv − m z [ f ] + ω ′ ) ǫ − ( ω ′ )15 − i h e i (2 Du − m z [ f ]) t D ( v − u ) ǫ − ( − Du + m z [ f ]) − e i (2 Dv − m z [ f ]) t D ( v − u ) ǫ − ( − Dv + m z [ f ]) i , (72)so that in the limit t → ∞ , we have A = − i π √ − u ∂∂u h √ − u f ′ ( u ) Z − d v (1 − v ) f ( v )2 D ( v − u ) ǫ +1 (2 Dv − m z [ f ]) × h e − i D ( v − u ) t ǫ − ( − Du + m z [ f ]) − ǫ − ( − Dv + m z [ f ]) ii = − π D √ − u ∂∂u h √ − u (1 − u ) f ( u ) f ′ ( u ) ǫ +1 (2 Du − m z [ f ]) ǫ − ( − Du + m z [ f ]) i , (73)where, in obtaining the second equality, we have used the lemma (59).Now, A may be computed along the same lines as done for B in the precedingsubsection. One getsRe( A ) = i √ − u ∂∂u h √ − u f ( u ) ǫ +1 (2 Du − m z [ f ] i = π D √ − u ∂∂u h √ − u (1 − u ) f ( u ) f ′ ( u ) ǫ +1 (2 Du − m z [ f ]) ǫ − ( − Du + m z [ f ]) i , (74)where, in obtaining the second equality, we have used equation (64). From equations (73)and (74), we see that A + Re( A ) = 0. Similarly, one can show that the real part of thecontribution to D δg u ∂δf∂u E from ( k = − , l = +1) also vanishes, while the imaginary part ofthe contribution to D δg u ∂δf∂u E from ( k = − , l = +1) must cancel that from ( k = +1 , l = − D δg u ∂δf∂u E = 0 . (75)Combining equations (66) and (75), we see that the Lenard-Balescu operator identicallyvanishes for our model, which is the announced result.
6. Example of a Vlasov-stationary state: Relaxation to equilibrium
In this section, let us consider as an example of an axisymmetric Vlasov stationary state f a state prepared by sampling independently for each of the N spins the angle φ uniformlyover [0 , π ] and the angle θ uniformly over an interval of length ( a + b ) asymmetric about θ = π/
2, that is, θ ∈ [ π/ − a : π/ b ]. The corresponding single-spin distribution is f ( u, φ ) = 12 π p ( u ) , (76)with p ( u ), the distribution for u , given by p ( u ) = a +sin b if u ∈ [ − sin b, sin a ] , . (77)16t is easily verified that this state has the energy e = (cid:16) D − (cid:17) (sin a + sin b ) − (cid:16) D − (cid:17) sin a sin b. (78)Since m x [ f ] = m y [ f ] = 0, we have g u = 0, and also, ∂f /∂φ = 0; it then follows that thestate (76) is stationary under the Vlasov dynamics (23). Note that we have m z [ f ] = sin a − sin b . (79)As mentioned after equation (37), the condition Im( ω p ) = 0 will correspond to themarginal stability of the state (76), so that the zeros ω p of the dielectric function lie on thereal- ω axis. Let us denote these zeros as ω ∗ pr . From equation (34), we find that ω ∗ pr satisfies1 + π P Z − d u (1 − u ) f ′ ( u )2 Du − m z [ f ] ∓ ω ∗ pr ± iπ (1 − u ) f ′ ( u ) (cid:12)(cid:12) ( m z [ f ] ± ω ∗ pr ) / (2 D ) = 0 . (80)Equating the real and the imaginary parts to zero, we get1 + π P Z − d u (1 − u ) f ′ ( u )2 Du − m z [ f ] ∓ ω ∗ pr = 0 , (81)(1 − u ) f ′ ( u ) (cid:12)(cid:12) ( m z [ f ] ± ω ∗ pr ) / (2 D ) = 0 . (82)Now, equation (76) gives f ′ ( u ) = 12 π (sin a + sin b ) [ δ ( u + sin a ) − δ ( u − sin b )] , (83)so that we obtain from equation (82) that δ (cid:16) m z [ f ] ± ω ∗ pr D + sin a (cid:17) = δ (cid:16) m z [ f ] ± ω ∗ pr D − sin a (cid:17) , (84)implying that ω ∗ pr = − m z [ f ] . (85)Using equation (85) in equation (81), we get4 D = cos a sin b + cos b sin a sin a sin b (sin a + sin b ) , (86)which when combined with equation (78) gives the energy e = e ∗ = AB ; A = 2(sin a − sin b ) cos a + 2 cot a sin b cos a + 2 cos b cot b sin a + 2 cos b (sin b − sin a ) − a − sin b ) (sin a + sin b ) ,B = 24(sin a + sin b ) , (87)17or which the state (76) is a marginally stable stationary solution of the Vlasov equation(23). For energies e > e ∗ , such a state is linearly stable under the Vlasov dynamics. Onthe basis of our analysis in this paper showing the Lenard-Balescu operator being identicallyzero, we expect that for finite N , the state relaxes to Boltzmann-Gibbs equilibrium on atimescale ∼ N δ , with δ >
1. This was indeed observed in Ref. [22] for the class of initialstates (76) that is non-magnetized, that is, a = b ; in this case, combining equations (86) and(87), we get e ∗ | a = b = D D . (88)For energies e < e ∗ | a = b , the state being linearly unstable under the Vlasov dynamics was seento relax for finite N to the Boltzmann-Gibbs equilibrium state over a timescale ∼ ln N [22]. -9 -8 -7 -6 -5 -4 -3 -2 -1 m ( t ) tN -2 N=100=200=500=1000=2000=4000
Figure 3.
Magnetization m ( t ) as a function of tN − in the energy range in which thestate (76) with a = 0 . , b = 0 . e = 0 . D = 8 . ∼ N δ with δ = 2, but the data do not allow a precise determination of theexponent δ : any value of δ between about 1 . Let us choose a = 0 . , b = 0 .
3. Then, equation (87) gives D ≈ . e ∗ ≈ . D , the state (76) with a = 0 . , b = 0 . e ∗ ≈ . e ∗ < e < e c , where e c can be computed from equation(3) to be e c ≈ . N . For e = 0 .
22, results of numericalsimulations of the dynamics shown in Fig. 3 indeed suggest a relaxation timescale ∼ N δ ,18ith δ >
1; for the range of system sizes explored in this numerical experiment, any valueof δ between about 1 .
7. Conclusions
In this paper, we have shown that the Lenard-Balescu operator identically vanishes for asystem of globally coupled anisotropic Heisenberg spins, in an axially symmetric Vlasov-stable state. This result explains the numerical findings of [22], reporting a relaxationtime for this system that scales superlinearly with N . To our knowledge, it is the firsttime that this kind of results has been obtained for a spin dynamics. This raises furtherquestions, e.g., what are the general conditions to ensure that the Lenard-Balescu operatorvanishes? The classical explanation relies on the structure of resonances between the particletrajectories: in the absence of resonances between particles with different momentum, theLenard-Balescu operator should vanish. This heuristic argument applies to systems ofparticles moving in a 1d position space, thus with a 2d phase space, when the system ishomogeneous [17, 30], implying a relaxation time growing superlinearly with N . This is alsothe case for axisymmetric configurations of point vortices [19–21], where the phase spaceis again two-dimensional. In a similar manner, it can be argued for the model we havestudied that spins with different projections on the z -axis cannot exchange energy becausethey cannot be in resonance. In a sense, our precise computations validate this qualitativepicture. However, recent numerical simulations of a model with a 4d phase space havealso shown a relaxation time that appears superlinear in N over the range of system sizesstudied [14]: one would expect resonances to appear in this case. Thus, understanding thegeneral conditions under which the Lenard-Balescu operator vanishes may still remain apartly open question.One may also wonder how the relaxation occurs when the Lenard-Balescu operatorvanishes. Formally, the Klimontovich expansion suggests that the next leading term is oforder 1 /N . Although writing down this term is possible in principle, its evaluation isdifficult. However, it is not quite clear that the expansion is valid over such long timescales.Finally, let us stress that the standard route to a formal derivation of the Lenard-Balescu equation, as followed in this article, involves an averaging over initial conditions.Just as what happens for the Vlasov equation, one may actually expect that the equationapproximately describes a single initial condition. Putting this on firm mathematical groundsis an outstanding question, on which some preliminary progress has been made recently [31].
8. Acknowledgements
SG acknowledges the support of the Indo-French Centre for the Promotion of AdvancedResearch under Project 4604-3 and the hospitality of Laboratoire J. A. Dieudonn´e, Universit´e19e Nice-Sophia Antipolis. [1] Campa A, Dauxois T and Ruffo S 2009
Phys. Rep.
Physica A
Int. J. Mod. Phys. B Long-Range Interacting Systems , 2010 ed by T Dauxois, S Ruffo and L F Cugliandolo(Oxford University Press, New York)[5] Bouchet F and Venaille A 2012
Phys. Rep.
Long-Range Interacting Systems , 2010 ed by T Dauxois, S Ruffo and L F Cugliandolo(Oxford University Press, New York)[7] Lenard A 1960
Ann. Phys. (N.Y.) Phys. Fluids Phys. Rev. E Phys. Rev. E Physica A
J. Stat. Mech.: Theory Exp.
P11008[14] Gupta S and Mukamel D, e-print:arXiv:1309.0194[15] Joyce M and Worrakitpoonpon T 2010
J. Stat. Mech.
P10012[16] Chandrasekhar S 1944
Astrophys. J. , 47[17] Eldridge O C and Feix M 1963 Phys. Fluids , 398[18] Rocha Filho T M et al. 2013, preprint arXiv:1305.4417[19] Dubin D H E and O’Neil T 1988 Phys. Rev. Lett. , 1286[20] Chavanis P H 2012 J. Stat. Mech.: Theory Exp.
P02019[21] Chavanis P H 2012
Physica A , 3657[22] Gupta S and Mukamel D 2011
J. Stat. Mech.: Theory Exp.
P03015[23] Campa A, Khomeriki R, Mukamel D and Ruffo S 2007
Phys. Rev. B , 064415[24] Kac M, Uhlenbeck G E and Hemmer P C 1963 J. Math. Phys. J. Math. Phys. Introduction to Plasma Physics (Krieger Publishing Company, Florida)[27] Neunzert H and Wick J 1972 Lecture Notes in Math. , Springer, Berlin[28] Braun W and Hepp K 1977
Comm. Math. Phys. Funct. Anal. Appl. Phys. Rev. E J. Stat. Phys.643