Shun Ogawa

Julien Barré

Hidetoshi Morita

Yoshiyuki Y Yamaguchi

Dynamical pattern formations in two dimensional fluid and Landau pole bifurcation

Keywords: numbers: 47.54.-r, 47.15.ki, 05.45.-a

A phenomenological theory is proposed to analyze the asymptotic dynamics of perturbed inviscid Kolmogorov shear flows in two dimensions. The phase diagram provided by the theory is in qualitative agreement with numerical observations, which include three phases depending on the aspect ratio of the domain and the size of the perturbation: a steady shear flow, a stationary dipole, and four traveling vortices. The theory is based on a precise study of the inviscid damping of the linearized equation and on an analysis of nonlinear effects. In particular, we show that the dominant Landau pole controlling the inviscid damping undergoes a bifurcation, which has important consequences on the asymptotic fate of the perturbation.

I. INTRODUCTION

Patterns in effective two-dimensional (2D) fluids flows are found in nature in various contexts [1,2]: atmospheric [3,4] and oceanic flows [5] are examples. To understand such large scale patterns theoretically, the 2D Euler equation describing perfect fluids flows is a simplified starting point. In this context, the study of nonlinear structures, such as Kelvin's cat's eyes, over a shear flow has a long history [6][START_REF] Lamb | Hydrodynamics[END_REF][START_REF] Stuart | [END_REF], which includes more recent mathematical developments [START_REF] Holm | Nonlinear systems of partial differential equations in applied mathematics[END_REF]10].

Statistical physics has often been invoked to explain the formation of large vortices since Onsager [11,12]. The Miller-Robert-Sommeria (MRS) theory [13][14][15][16], which constructs the microcanonical measure for 2D Euler flows by taking all the invariants into account, is a particularly successful achievement. A difficulty in applying this theory in practice is the fact that the 2D Euler equation has infinitely many invariants. More importantly, the theory assumes the vorticity field on the large scales to be stationary, and therefore cannot describe non-stationary asymptotic behaviors. For these cases, the statistical physics approach should be supplemented by a dynamical understanding of pattern formations.

In this paper, we consider a shear base flow to which a perturbation is added, and propose a phenomenological approach to the pattern formations by analyzing the Landau pole of the linearized equation with its bifurcation, and by taking into account nonlinear effects.

Detailed numerical simulations of this situation are presented in Ref. [17]: starting from a Kolmogorov flow on the doubly periodic domain T 2 = [0, 2π) × [0, 2πΓ), a perturbation of size ǫ is added; depending on the parameters (ǫ, Γ), the perturbation may fully damp, evolve into a stationary dipole, or create four long lived traveling vortices. The regions on the (ǫ, Γ) plane where the dipole or traveling vor- * E-mail: sogawa@amp.i.kyoto-u.ac.jp tices appear are numerically investigated in Ref. [17], but a theoretical understanding is lacking. In this article, we aim at providing such a theoretical explanation, by exploiting the analogy between the 2D Euler equation for the vorticity field and the Vlasov equation for plasmas. Indeed, a similar phenomenon has been investigated in one-dimensional plasmas, described by the Vlasov equation. A naive linear theory predicts that perturbations added to stable stationary states damp exponentially [18]. It is well known, however, that, if the perturbation is large enough, nonlinear effects come into play, prevent complete damping, and may create traveling clusters [19][20][21][22]; see also Ref. [23] for the simpler cases of ferromagnetic and anti-ferromagnetic Hamiltonian mean-field model. To be more precise, according to Refs. [19,23], such a phenomenon occurs when the two following criteria are satisfied:

1. The Landau damping time scale is longer than a nonlinear timescale, called "trapping timescale".

2. If several clusters are formed, they should be so small that they do not overlap and that a nonlinear superposition approximation [24] may work.

Although these criteria are in view of nonlinear dynamics, they can be expressed using the dominant Landau pole computed from the linearized Euler equation. Our goal is to draw the phase diagram in the (ǫ, Γ) plane by using the above two criteria in the context of the Euler equation, that is, by combining the phenomenological nonlinear estimate and the linear Landau damping theory. This article is organized as follows. In Sec. II, we qualitatively analyze criteria 1 and 2, and show that they are related to the imaginary and the real parts of the dominant Landau pole of the initial Kolmogorov flow. We therefore briefly review the linear theory for 2D incompressible and inviscid fluids, and derive the dispersion relation in Sec. III; the computation of the Landau pole requires an analytic continuation procedure similar to the one used in Ref. [25]. Using these computations, we draw a phase diagram in the (ǫ, Γ) plane in Sec. IV, and compare this phenomenological estimate with numerical simulations in Sec. V. Section VI is devoted to the summary and discussions.

II. THE CRITERIA TO OBSERVE A DIPOLE, OR TRAVELING VORTICES

We start from the 2D Euler equation in the domain

T 2 ∂ω ∂t + v • ∇ω = 0, (1)
where the vorticity field ω and the velocity field v are related to a stream function ψ through

ω = ∇ 2 ψ, v = - ∂ψ ∂y , ∂ψ ∂x . ( 2
)
The periods for x and y axes are set as 2π and 2πΓ respectively.

We consider the stationary Kolmogorov flow, called hereafter the "base flow", whose stream function is ψ 0 = -Γ sin(y/Γ) and vorticity and velocity fields are respectively

ω 0 (y) = -U ′ (y), v 0 = (U (y), 0), (3)
with U (y) = cos(y/Γ). We add a large scale perturbation to the base flow, and expand the stream function as

ψ(x, y, t ) = ψ 0 (y) + ψ 1 (x, y, t ), (4)
with

ψ 1 (x, y, t = 0) = ǫ cos x. ( 5
)
In our analytical computations, ǫ is assumed to be small; this restriction obviously does not hold for numerical simulations. Nevertheless, the theory will be qualitatively in good agreement with the simulations. When the base flow is stable (this corresponds to an aspect ratio Γ > 1), the linear theory typically predicts that the perturbation damps and possibly oscillates at complex frequency c = c R + i c I , where c is the root of the dispersion relation yielding the slowest damping, that is, the root closest to the real axis, or, in other words, with the smallest |c I |. The idea is that while it damps, the perturbation will tend to create structures traveling at velocity c R ; if nonlinear effects are strong enough, these structures may persist for long times. If c R = 0, one may expect traveling vortices, while the formation of two stationary vortices (a dipole) should be favored if c R = 0.

A. Criterion 1: damping timescale longer than trapping timescale

The damping time scale τ L is easily estimated as the inverse of the Landau damping rate τ L ≃ 1/|c I |. Actually, in addition to the exponential Landau damping described by c I , there is an algebraic damping coming from the branch points of the dispersion function (see Sec. III); however, the exponential Landau damping is expected to be dominant on the time scale considered in this paper.

Next, we estimate the "trapping time scale", τ T , which is the characteristic time scale concerning nonlinearities [19], as the period of a test point vortex trapped around the edge of a small vortex. The temporal evolution of the position of a test point vortex is governed by the velocity field

ẋ = - ∂ψ ∂y = U (y) - ∂ψ 1 ∂y , ẏ = ∂ψ ∂x = ∂ψ 1 ∂x , (6)
where the dot denotes d/dt . We approximate the perturbation ψ 1 (x, y, t ) phenomenologically. We are interested in the macroscopic behavior corresponding to the k = ±1 mode,

where k is the wave number with respect to x, and we assume that the damping perturbation has created small vortices whose velocity is c R in the x direction; note that c R can be 0. The first order ψ 1 (x, y, t ) can be, therefore, approximated by

ψ 1 (x, y, t ) ≃ ǫ ψ1 (y) cos(x -c R t ), U (y) ∼ c R , (7)
where ψ1 (y) is of O(ǫ 0 ). Substituting Eq. ( 7) into Eq. ( 6), we have the approximate equations of motion

ẋ = U (y) -ǫ ψ′ 1 (y) cos(x -c R t ), ẏ = -ǫ ψ1 (y) sin(x -c R t ), (8)
where the prime represents d/dy. Thanks to the phenomenological approximation of ψ 1 [START_REF] Lamb | Hydrodynamics[END_REF], we can compute the motion of the test point vortex by perturbation techniques. Expanding x and y into series of ǫ as

x(t ) = x 0 (t ) + ǫx 1 (t ) + O(ǫ 2 ), y(t ) = y 0 (t ) + ǫy 1 (t ) + O(ǫ 2 ), (9)
we obtain the solutions x 0 (t ) and y 0 (t ) with initial conditions x 0 (0) = X and y 0 (0) = Y :

x 0 (t ) = X +U (Y )t , y 0 (t ) = Y . ( 10
)
Equations for the first order in ǫ are:

ẋ1 (t ) = U ′ (Y )y 1 (t ) -ψ′ 1 (Y ) cos (X + (U (Y ) -c R )t ) , ẏ1 (t ) = -ψ1 (Y ) sin (X + (U (Y ) -c R )t ). ( 11
)
From the frequency of y 1 in Eq. ( 11), we estimate the trapping time scale τ T as

τ T ≃ 1
|U (Y * + g (ǫ, Γ)) -c R | ≃ 1 g (ǫ, Γ)|U ′ (Y * )| ( 12
)
where Y * satisfies U (Y * ) = c R , and g (ǫ, Γ) > 0 represents the width of a vortex.

The unknown quantity, the width of a vortex g (ǫ, Γ), is determined self-consistently. The solution y 1 (t ) is, from Eq. ( 11),

y 1 (t ) = ψ1 (Y ) U (Y ) -c R [cos(X + (U (Y ) -c R )t ) -cos X ] , ( 13
)
where the initial condition is y 1 (0) = 0. The amplitude in the y direction, ǫ ψ1 (Y )/(U (Y ) -c R ), must be the same as the width g (ǫ, Γ) at the edge of the vortex, Y = Y * + g . Thus, the self-consistent equation for g is

g = ǫ U (Y * + g ) -c R , ( 14
)
where we have introduced another phenomenological approximation by replacing ψ1 (Y ) with 1, since it is of O(ǫ 0 ) and we are looking for an order-of-magnitude estimate. The width g is, therefore, estimated as

g (ǫ, Γ) ≈ ǫ |U ′ (Y * )| ( 15
)
For the Kolmogorov base flow, U (y) = cos(y/Γ), we have

|U ′ (Y * )| = 1 Γ 1 -c 2 R , ( 16
)
and the trapping time scale is

τ T = Γ ǫ (1 -c 2 R ) -1/4 . ( 17
)
Criterion 1 reads τ T < τ L , that is:

Γc 2 I 1 -c 2 R < ǫ. ( 18
)
If c R = 0, this condition simplifies into Γc 2 I < ǫ.

B. Criterion 2: non overlapping vortices

According to Eq. ( 5), we consider the modes k = ±1. We assume that the frequency c R of the perturbation reflects the existence of periodically moving vortices, thanks to the periodic boundary condition for the x axis. The period 2π permits to identify the frequency with the velocity of the moving vortices, and therefore the y positions of the vortices are estimated as the solutions of the equation U (y) = ±c R . We recall the base flow is U (y) = cos(y/Γ).

If c R = 0, one expects that vortices are formed at y = πΓ/2 and y = 3πΓ/2. The estimate of the vortices width (15) shows that the vortices will never overlap for any reasonably small ǫ (say for instance ǫ < 0.5). Hence criterion 2 brings no restriction in this case.

If c R = 0, in contrast, the vortices can be close one to another, and criterion 2 leads to a restriction. We We assume that each vortex has the same width g . Then the distance between two nearby vortices must be larger than 2g in order to avoid an overlap. Using π/2arccos(c R ) = arcsin(c R ), this condition is expressed as

Γ arcsin(c R ) > g , (|c R | < 1/ 2), Γ arccos(c R ) > g , (|c R | > 1/ 2). ( 19
)
The first inequality comes from the distance in y between vortices A and B (identically C and D), and the second from the distance between B and C (identically A and D). The pair (B,C) moves toward the left while (A,D) toward the right according to the base flow U (y). Moreover, the numerical observations [17] say that the difference in x among each pair stays around π. Thus the distance between B and C (A and D) is large enough to neglect the effect of overlapping; see Fig. 1. Hence we can omit the second inequality of Eq. ( 19), and criterion 2 is finally expressed as:

ǫ < Γ 1 -c 2 R (arcsin(c R )) 2 . ( 20
)
The two conditions ( 18) and ( 20) involve the parameters Γ and ǫ, as well as the dominant Landau pole c = c R + i c I . In the following section, we turn to the computation of this Landau pole.

III. THE LINEAR THEORY AND THE DISPERSION RELATION

This section contains a classical computation for the linearized 2D Euler equation, as well as an analytic continuation in the spirit of Ref. [25], for the self-consistency of the paper. We perform it for a general base flow and perturbation, and then specialize it to Eqs. ( 3) and (5).

A. The dispersion function D(c)

We linearize Eqs. ( 1) and ( 2) around the base flow (3). We add a small perturbation ψ 1 to ψ 0 (y), whose associated vorticity and velocity fields are denoted respectively as ω 1 and v 1 . Substituting ω = ω 0 + ω 1 and v = v 0 + v 1 into the Euler equation ( 1) gives the linearized 2D Euler equation

∂ω 1 ∂t + v 0 • ∇ω 1 + v 1 • ∇ω 0 = 0. ( 21
)
Using the perturbation ψ 1 and Eq. ( 3), we rewrite the linearized Euler equation ( 21) as

∂ω 1 ∂t +U (y) ∂ω 1 ∂x -U ′′ (y) ∂ψ 1 ∂x = 0. ( 22
)
Our goal is now to derive the dispersion relation for the linearized Euler equation; we use a Fourier-Laplace transformation, and follow the route of Ref. [26]. Thanks to the periodicity with respect to x, we expand ω 1 and ψ 1 into Fourier series as

ω 1 (x, y, t ) = k∈Z ωk (y, t )e ik x , ψ 1 (x, y, t ) = k∈Z ψk (y, t )e ik x . ( 23
)
Substituting the Fourier expansions (23) into the linearized Euler equation (22), we obtain the equation for k-th Fourier mode as

∂ ωk ∂t + i kU (y) ωk -i k ψk U ′′ (y) = 0, ∀k ∈ Z. (24)
The Laplace transform of a function ĝ (t ) with respect to t is defined as

g (z) = ∞ 0 ĝ (t )e -zt d t , Re z > 0, ( 25
)
where the condition Re z > 0 is introduced to ensure the convergence of the integral. Performing the Laplace transform of the linearized Euler equation (24), we obtain an algebraic equation for ωk and ψk :

U (y) -c ωk -U ′′ (y) ψk = ωk (y, 0) i k , ( 26
)
where z = -i kc; note that Im(kc) > 0.

The Fourier transform of Eq. ( 2) with respect to

x ωk (y, c) = ∂ 2 ∂y 2 -k 2 ψk (y, c) (27)
gives a closed equation for ψk (y, c), which is the Rayleigh equation

∂ 2 ψk ∂y 2 -k 2 ψk - U ′′ (y) U (y) -c ψk = ωk (y, 0) i k(U (y) -c) . ( 28
)
The stream function ψk (y, t ) is computed through an inverse Laplace transform. Its asymptotic behavior is determined by the singularities of ψk (y, c) with respect to the complex variable c. For simplicity, we introduce the functions q and f as

q(y) ≡ k 2 + U ′′ (y) U (y) -c , f (y) ≡ ωk (y, 0) i k(U (y) -c) . ( 29
)
These functions have no singularity for real y, since c is defined in the region Im (kc) > 0.

Fixing the complex variable c, we analyze the Rayleigh equation of the form

d 2 φ dy 2 -q(y)φ = f (y), (30)
and the corresponding homogeneous Rayleigh equation

d 2 φ dy 2 -q(y)φ = 0. ( 31
)
Let φ 1 and φ 2 be independent solutions to the homogeneous equation (31) with boundary conditions

φ 1 (0) = 1, φ ′ 1 (0) = 0, and φ 2 (0) = 0, φ ′ 2 (0) = 1. (32)
The particular solution φ p to the inhomogeneous equation ( 30) is then given by

φ p (y) = -φ 1 (y) y 0 φ 2 (y ′ ) f (y ′ ) d y ′ + φ 2 (y) y 0 φ 1 (y ′ ) f (y ′ ) d y ′ . ( 33
)
Indeed, the double derivative of φ p (y) is

φ ′′ p (y) = q(y)φ p (y) + W (y) f (y), (34)
where W (y) is the Wronskian

W (y) = φ 1 (y)φ ′ 2 (y) -φ ′ 1 (y)φ 2 (y) (35)
and is constant

W (y) = W (0) = 1 for all y. ( 36
)
Hence, the general solution φ g to Eq. ( 28) is:

φ g = φ p + a 1 φ 1 + a 2 φ 2 , ( 37
)
where a 1 and a 2 are constants determined from the periodic boundary condition

φ g (2πΓ) = φ g (0), φ ′ g (2πΓ) = φ ′ g (0). ( 38
)
From the boundary conditions (32) and φ p (0) = φ ′ p (0) = 0, Eqs. (37) and (38) lead to:

φ 1 (2πΓ) -1 φ 2 (2πΓ) φ ′ 1 (2πΓ) φ ′ 2 (2πΓ) -1 a 1 a 2 = - φ p (2πΓ) φ ′ p (2πΓ) . ( 39
)
Remembering that q(y) and hence φ 1 and φ 2 depend on c, we define the function D(c) as

D(c) = det φ 1 (2πΓ) -1 φ 2 (2πΓ) φ ′ 1 (2πΓ) φ ′ 2 (2πΓ) -1 = 2 -φ 1 (2πΓ) -φ ′ 2 (2πΓ). ( 40
)
To show the last equality in Eq. ( 40), we have used the fact that the Wronskian W (y) is unity. The values of φ 1 (2πΓ) and φ ′ 2 (2πΓ) are computed by integrating the homogeneous Rayleigh equation ( 31) from 0 to 2πΓ.

B. Analytic continuation of D(c)

The general solution (37) has singularities for c satisfying D(c) = 0, and these singularities yield non-trivial modes proportional to e -ikc t by the inverse Laplace transform. This justifies our terminology "dispersion relation" for D(c) = 0, and "dispersion function" for the function D(c). Recall that the dispersion function D(c) is a priori defined only in the region Im(kc) > 0 to ensure the convergence of the Laplace transform. Now, to find roots giving stable modes, we analytically continue D(c) to the whole complex c plane, in a similar manner to what has been performed for the Euler equation on the 2D disc [25,27,28].

Hereafter we set k = 1 without loss of generality, and take the particular form of the base flow: U (y) = cos(y/Γ). Corresponding to a given c, the functions q(y) and f (y) have singularities at which the equation U (y) -c = 0 is satisfied. Such singular points are, together with the integration paths of the homogeneous Rayleigh equation (31), schematically illustrated in Fig. 2 for both the cases Im c > 0 and Im c < 0. The integration paths are determined as follows.

For c on the upper half plane, D(c) is simply defined by integrating Eq. ( 31) along the real y axis, namely along the integration path:

L = y ∈ R | y = y(s) = 2πΓs, s ∈ [0, 1) . ( 41
)
Continuously moving c to the lower half plane allows to define D(c) for Im c < 0. If the path taking c from the upper to the lower half plane crosses the singular interval c ∈ [-1, 1], we have to avoid the singularity that would be created. For this purpose, we deform the integration path L to

L h = y ∈ C | y = y h (s) = 2πΓ(s + i h(s)), s ∈ [0, 1) , ( 42
)
using a C 2 -class, real-valued, and 1-periodic function h(s) satisfying h(0) = h(1) = 0; correspondingly, this amounts to deform the singular line

σ = {c ∈ R | c = U (y(s)), s ∈ [0, 1)}, (43)
which doubly covers the interval [-1, 1], to the curve

σ h = {c ∈ C | c = U (y h (s)), s ∈ [0, 1)}. ( 44
)
We choose the function h(s) so that the deformed curve σ h is below the c for which we want to define D(c). This procedure is similar to the spectral deformation for the Vlasov-Poisson equation in Refs. [29,30].

If the path taking c from the upper to the lower plane crosses the real axis through either one of the half lines (-∞, -1) or (1, ∞), no singularity is crossed, and the analytic continuation does not require any deformation of the integral path; in the lower half plane, D(c) is multi valued unless a Riemann surface associated with the branch points c = ±1 is introduced.

We choose the particular deformation function h(s) in Eq. (42) as

h(s) = a sin 2πs, a ≥ 0. ( 45
)
The extended domains of D(c) are shown in Fig. 3 for two different values of a. We use a large enough so that σ h is below c; the computational cost increases exponentially with increasing a.

By introducing the following functions of s:

ϕ(s) = φ(y h (s)), p(s) = q(y h (s)), (46)
the homogeneous Rayleigh equation ( 31) is rewritten as

ϕ ′′ R - h ′ h ′′ 1 + h ′2 ϕ ′ R + h ′′ 1 + h ′2 ϕ ′ I -(2πΓ) 2 (1 -h ′2 )p R -2h ′ p I ϕ R + (2πΓ) 2 2h ′ p R + (1 -h ′2 )p I ϕ I = 0, ϕ ′′ I - h ′′ 1 + h ′2 ϕ ′ R - h ′ h ′′ 1 + h ′2 ϕ ′ I -(2πΓ) 2 2h ′ p R + (1 -h ′2 )p I ϕ R -(2πΓ) 2 (1 -h ′2 )p R -2h ′ p I ϕ I = 0, ( 47
)
where the subscripts "R" and "I" represent the real and imaginary parts respectively, and the prime denotes the derivative with respect to s. The continued solutions ϕ 1 and ϕ 2 are computed by solving Eq. ( 47) with the boundary con- Corresponding to a given cross point on the c plane, there are two singular points (crosses) of q(y) and f (y) on the y plane in general. As the point on the c plane goes down, the singular points on the y plane pass the real axis and the integration path L must be deformed to L h to avoid the singularities. The mapping of L h , denoted by U (L h ), determines the boundary of the continued domain of the dispersion function D(c). See Fig. 3. A branch cut on the c plane is set as the boundary.

ditions

           ϕ 1,R (0) = 1, ϕ 1,I (0) = 0, ϕ ′ 1,R (0) = 0, ϕ ′ 1,I (0) = 0,
and

           ϕ 2,R (0) = 0, ϕ 2,I (0) = 0, ϕ ′ 2,R (0) = 2πΓ, ϕ ′ 2,I (0) = 2πΓh ′ (0), (48)
respectively. Solving Eq. ( 47), we obtain ϕ 1 (1) and ϕ ′ 2 (1), and the dispersion function is expressed as

D(c) = 2 -ϕ 1 (1) - ϕ ′ 2 (1) 2πΓ(1 + i h ′ (1)) , (49)
whose real and imaginary parts are respectively

D R (c) = 2 -ϕ 1,R (1) - ϕ 2,R (1) + h ′ (1)ϕ 2,I (1)
2πΓ(1

+ h ′ (1) 2 ) , D I (c) = -ϕ 1,I (1) + h ′ (1)ϕ 2,R (1) -ϕ 2,I (1) 2πΓ(1 + h ′ (1) 2 ) . (50)

C. The main Landau pole and its bifurcation

We now specialize the above computations to the initial perturbation We look for the dominant Landau pole, i.e., the root of D(c) with the largest imaginary part; in a stable situation, it corresponds to the slowest damping. The variation of this pole as a function of the aspect ratio Γ is shown on Fig. 4. For Γ < 1, the flow is unstable. With increasing Γ above 1, the flow becomes stable, with the main Landau pole on the imaginary axis. With further increasing Γ, the dominant Landau pole undergoes a bifurcation at an aspect ratio Γ c ≃ 1.06 and acquires a non-zero real part. We will see that this bifurcation is crucial to understand the appearance of a stationary dipole or of the four traveling vortices.

ω 1 (x, y, 0) = -ǫ cos x, (51) 0 -0.6 0.6 -1 0 1 -1 0 1 0 -1 -2 -3 (a)

IV. THE PHASE DIAGRAM

Using the theoretical analyses above, we construct the phase diagram on the (ǫ, Γ) plane, which consists of zonal, dipole and oscillation phases [17]. The last oscillation phase corresponds to the four traveling vortices.

To obtain the region on the (ǫ, Γ) plane where the appearance of four traveling vortices is expected, we combine criteria 1 and 2, that is Eqs. (18) and (20) respectively, with the result on the dominant Landau pole. Satisfying both criteria then reads:

Γc 2 I 1 -c 2 R < ǫ < 1 -c 2 R (arcsin(c R )) 2 . ( 52
)
The right inequality requires a Landau pole with non-zero c R , which implies Γ > Γ c . The region which satisfies Eq. (52) on the (ǫ, Γ) plane is reported in Fig. 5, which is qualitatively in good agreement with the numerically obtained oscillatory region [17]. We remark that |c R | increases when Γ is increased above Γ c , staying smaller than 1/ 2 in the reported parameter region. This fact validates that we have omitted the lower condition in Eq. ( 19) to derive the second necessary condition.

We expect to find a stationary dipole when criterion 1 is satisfied and c R = 0:

Γc 2 I < ǫ and Γ < Γ c . (53)
The region satisfying these requirements is highlighted in Fig. 5. This theoretically estimated region is again in qualitative agreement with the region where a stationary dipole emerges in numerical simulations [17].

We expect a zonal flow when criterion 1 is not satisfied, that is

Γc 2 I > ǫ, for Γ < Γ c , Γc 2 I 1 -c 2 R > ǫ, for Γ > Γ c . ( 54
)
Physically, it means that nonlinearity is not dominant and hence simple damping of the perturbation is expected. In Fig. 5 the zonal region appears for small ǫ. On Fig. 5, there is a region where the theory makes no prediction; indeed, the vortices that would be created by the phenomenological mechanism considered above would be so close one to another that they would strongly interact; the final fate of the system is then out of scope of the present theory. We note that in this region, zonal flows, dipoles and traveling vortices are observed numerically [17].

V. NUMERICAL TESTS

According to the phenomenological prediction of the preceding sections, the four traveling vortices run in the x direction along the four lines y = y * in T 2 = [0, 2π) × [0, 2πΓ), where y * = Γ arccos(±c R ), and their velocities are U (y * ) = ±c R . The period in the x direction is 2π; hence the phenomenological theory predicts the frequency f pt = |c R |/2π for the vorticity field ω.

We examine our phenomenological prediction by numerically observing the y-positions of the four vortices and the oscillation frequency. The initial condition of the vorticity field is: the numerically computed bump positions.

ω(x, y, 0) = ω 0 (x, y) -ǫ cos x ( 55
Next, we obtain the frequency by computing the power spectra of the (1, 0) Fourier mode of ω(x, y, t ) defined by

Z (t ) = -Re 1 (2π) 2 Γ T 2 ω(x, y, t )e -i x d xd y. ( 57
)
As shown in Fig. 7, the dependence on Γ of f num is in qualitative agreement with the prediction, in the sense that the frequency increases as Γ increases. However, it is not in good agreement quantitatively. One possible explanation for the quantitative discrepancy in frequency is that ǫ = 0.3 is too large to be considered a small perturbation to the base flow. The present theory is based on the linear analysis of the Euler equation. It gives good predictions qualitatively, but quantitatively, nonlinear effects may kick in for rather large ǫ.

VI. SUMMARY AND DISCUSSIONS

Inspired by previous works on the Vlasov equation [19,[22][23][24]31], we have considered i) the balance between Landau damping and nonlinear effects and ii) the nonoverlapping criterion between nonlinear structures, in order to analyze the formation of a stationary dipole and traveling vortices in the 2D Euler equation starting from the perturbed stable Kolmogorov flow. The detailed linear analysis provides information on this nonlinear dynamical phenomenon. We emphasize that the bifurcation of the dominant Landau pole as the aspect ratio is varied plays a crucial role in the theory. Since this kind of bifurcation is probably not a special feature of the Kolmogorov flow, we expect that a similar analysis, predicting the appearance, or not, of a dipole or traveling vortices, would be fruitful for other base flows.

We note that our prediction is in qualitative, though not quantitative, agreement with the numerical results of Ref. [17], not only for the shape of the phase diagram, but also for the vortex frequency. The discrepancy for the phase diagram is rather natural since we have performed phenomenological order-of-magnitude estimates and not taken the prefactors into consideration. The reason for the quantitative disagreement in the frequency is probably that the size of the perturbation ǫ used in the numerical simulation is large enough to trigger nonlinear effects that we do not take into account, while such a large ǫ is necessary to realize the oscillating phase numerically.

  FIG. 1. Schematic picture of four traveling vortices, A, B, C, and D. The arrows indicate traveling directions. Vortices A and B have positive vorticity, while C and D negative. The solid curve represents the velocity of the base flow U (y), suitably rescaled. The two vertical dotted lines correspond to ±c R ; their intersections with the velocity curve, U (y) = ±c R , yield the estimated y positions of the four vortices, which travel along the horizontal dotted lines.

FIG. 2 .

 2 FIG. 2. Schematic picture of complex c and y planes for k = 1. The left two panels are for Im c > 0, and the right two panels for Im c < 0. The integration path L on the y plane, consisting of the solid and the dotted lines, is mapped to the double covering of the interval [-1,1] on the c plane by the mapping U : y → U (y).Corresponding to a given cross point on the c plane, there are two singular points (crosses) of q(y) and f (y) on the y plane in general. As the point on the c plane goes down, the singular points on the y plane pass the real axis and the integration path L must be deformed to L h to avoid the singularities. The mapping of L h , denoted by U (L h ), determines the boundary of the continued domain of the dispersion function D(c). See Fig.3. A branch cut on the c plane is set as the boundary.

FIG. 3 .

 3 FIG. 3. Extension of the domain of D(c). The gray regions represent the domains of D(c) for (a) a = 0.1 and (b) a = 0.3. The curve in the c plane represents U (y h (s)), where y h (s) is the deformed integration path for the homogeneous Rayleigh equation(31). With increasing a, the domain of D(c) is extended to a larger area in the lower half c plane, but the computational cost also increases.

FIG. 4 .

 4 FIG. 4. (color online)The dominant root of the dispersion function (red line) and a favorable region for the existence of traveling vortices (gray region) on the complex c plane. When Γ < 1, one dominant root is obtained on the upper half plane, where the base flow is unstable, and the root is at the origin for Γ = 1. With increasing Γ, the main root goes down to the lower half plane along the imaginary c axis, accompanied by another root going up the axis. When Γ = Γ c ≃ 1.06, these two roots collide at c = -0.28i and they bifurcate to the right and left directions. The gray region is given by the inequality, c 2 I < (1 -c 2 R )(arcsin(c R )) 2 coming from Eq. (52), and c I < 0 required by the stability of the flow. For later convenience, we marked the bifurcation point corresponding Γ = Γ c by the cross, and the point at which the main pole enters the gray region by thestars.

) fixing the small parameter ǫ = 0. 3 .FIG. 5 .

 35 FIG. 5. (color online)Phase diagram on the (ǫ,Γ) plane. The gray, polka-dotted and checkered regions correspond to the oscillatory, zonal and dipole phases, respectively. The upper panel is drawn by the present theory and the lower is a schematic picture of the phase diagram obtained with the numerical simulations by Morita[17]. In the upper panel, the gray region is obtained from the inequality (52). The left green boundary of the gray region is given by ǫ = Γc 2 I / 1 -c 2 R (see Eq. (18)), and the lower red by ǫ = Γ 1 -c 2 R (arcsin c R ) 2 (see Eq. (20)). The left blue boundary of the checkered region is given by ǫ = Γc 2 I , and the upper purple by Γ = Γ c (see Eq. (53)), whose value is determined at the bifurcation point marked by the cross in Fig.4. The star point corresponds to the star point of Fig.4. The broken line represents ǫ = 0.3 on which the theoretical and numerical frequencies are compared in Fig.7. The points marked by the cross and the circle correspond to those in Fig.7, respectively.

4 FIG. 6 .

 46 FIG. 6. (color online) The averaged vorticity field ω(y, t). Four bumps are observed, corresponding to four vortices in the long time regime. The vertical four dotted lines indicate the theoretical predictions for the y-positions of the vortices: y/L y = 0.25 ± 0.097,0.75 ± 0.097. Γ = 1.4 and ǫ = 0.3.

ACKNOWLEDGMENTS

The authors are grateful to F. Bouchet for fruitful discussions and comments on this study. SO thanks for the hospitality of Université de Nice Sophia-Antipolis and acknowledges the support of the JSPS Research Fellowships for Young Scientists (Grant No. 254728). This work was supported by the ANR-09-JCJC-009401 INTERLOP project. YYY acknowledges the support of a Grant-in-Aid for Scientific Research (C) 23560069.


Julien Barré

Alain Olivetti

Yoshiyuki Y Yamaguchi

Landau damping and inhomogeneous reference states

Landau damping is a fundamental phenomenon in plasma physics, which also plays an important role in astrophysics, and sometimes under different names, in fluid dynamics and other fields. Its theoretical discussion in the framework of the Vlasov equation often assumes that the reference stationary state is homogeneous in space. However, Landau damping around an inhomogeneous reference stationary state, a natural setting in astrophysics for instance, induces new mathematical difficulties and physical phenomena. The goal of this article is provide an introduction to these problems and the questions they raise.

Introduction

Vlasov equation and Landau damping have been an active subject of research since 70 years, and have been encountered at some point in their career by most physicists. The purpose of this article is to provide an introduction to the specificities of Landau damping close to a stationary state of the Vlasov equation which is not homogeneous in space. This situation is in particular relevant for astrophysical applications of the Vlasov equation, but it is striking that a phenomenology similar to the one we will describe here is encountered in very different physical systems: fluids described by Euler equation [1], Kuramoto model for oscillators synchronization [2], sound propagation in bubbly fluids [3]. . . For self-consistency, we will first recall some basic facts about Vlasov equation and standard Landau damping as it was first introduced in plasma physics, close to an homogeneous in space stationary state.

The Vlasov equation

When one tries to describe the positions and velocities of a large number of interacting particles by a density in phase space, different situations may occur: if the interactions between particles are strong and rare, so that they involve only two particles each time, one expects a Boltzmann-like equation; if on the contrary one particle feels the effect of many others, the individual effect of each one being weak, one enters the realm of Vlasov equation. This scaling limit, sometimes called mean-field limit, appears in many areas of physics, the two main examples being: particles in Coulombian interaction, and particles in Newtonian interaction. In the plasma physics context, it was indeed introduced by Vlasov [4]; in the astrophysical context, the Vlasov equation is usually called "collisionless Boltzmann equation". We write down the equation in Hamiltonian form:

∂ t f + {f, h} = 0, (1)
where f (x, v, t) is the density in phase space, normalized to 1, {f, g} is the

Poisson bracket {f, g} = 1 M (∇ x f • ∇ v g -∇ v f • ∇ x g) (2)
M is the particles' mass and h(x, v, t) the Hamiltonian

h(x, v, t) = M v 2 2 + φ(x, t). (3)
The potential φ(x, t) is created by f . For Coulomb or Newton interaction, φ is related to the spatial density ρ = f dv through the Poisson equation:

∆φ = Cρ (4)
where C = 4πG for Newton with G the gravitational constant, and -4πe 2 /M for Coulomb interaction. More generally, φ may be given by the convolution of the density ρ and a kernel K:

φ(x) = K(x -y)ρ(y)dy. (5)
For instance, K is proportional to 1/|x| for Coulomb (resp. -1/|x| for Newton) interaction in 3D. Physically, Vlasov equation is simply a transport equation, meaning that each particle moves in the field that all particles create. Hereafter we set M = 1 for simplicity. The convergence of the discrete dynamics of N particles towards the continuous Vlasov equation has been rigorously proved [5,6], but under the hypothesis that the kernel K is regular enough. In particular, this excludes the Coulomb and Newton cases!

Stationary states and their stability

The Vlasov equation (1) admits an infinite number of stationary states. Indeed, any function f = ϕ(h(x, v)) that depends on the phase space coordinates only through the Hamiltonian h is a stationary solution. An easy special case is provided by a function f that depends only on v 2 (this implies that f is homogeneous in space), with a vanishing self-consistent potential. Space inhomogeneous stationary solutions, needed for instance to describe self-gravitating systems, are more difficult to construct, since one needs to solve the self-consistent equation (4) for the potential.

Studying the Vlasov dynamics close to a stationary solution proved a difficult endeavor, with surprises even at the linear level. When the linearized Vlasov operator has an eigenvalue with positive real part, the corresponding stationary state is clearly unstable. Futhermore, it is easy to realize that the spectrum of the linearized operator is symmetric with respect to the imaginary axis; hence, if there is an eigenvalue with negative real part, there is also one with positive real part, and the corresponding stationary state is also unstable. Thus, stable stationary states are only marginally stable, the spectrum being included in the imaginary axis. Yet, a perturbation around a stable stationary state may decay in some sense exponentially, with a well defined rate: this is the well-known Landau damping phenomenon, first described in 1946 [7]. In the following, we consider such damping of a perturbation δf around a stable stationary state f 0 .

Landau damping, homogeneous case

We start with a simple setting, 1D and periodic in space, and consider a homogeneous in space stationary solution f 0 (v). The linearized Vlasov equation around f 0 for a small perturbation δf reads:

∂ t δf + v∂ x δf -f ′ 0 ∂ x δφ = 0 (6) δφ(x) = K(x -y)δρ(y)dy (7)
where δρ and δφ are the density and potential perturbations respectively.

To understand the basic mechanism of Landau damping, we start by forgetting the interaction, keeping only the free transport: each particle moves with constant velocity v. Velocities of initially close particles are different; this shear develops finer and finer structures in the (x, v) phase space; any initial spatial organization is thus damped by this phase mixing process: see an illustration Fig. 1. Introducing the interaction term does not modify much this basic physical mechanism, but has an important effect on the type of damping.

The standard mathematical treatment of ( 6)-( 7) is as follows. After a Fourier transform in space and a Laplace transform in time, which is defined for the function u(x, t) as

û(k, ω) = e -ikx ∞ 0 e iωt u(x, t)dt dx, (8)
the linearized Vlasov equation ( 6)-( 7) is reduced to a diagonalized algebraic equation and is solved formally. Simple manipulations give the Fourier-Laplace component as where

δφ(k, ω) = G(k, ω) D(k, ω) , (9)
D(k, ω) = 1 -K k kf ′ 0 (v) kv -ω dv, G(k, ω) = -iK k fini (k, v) kv -ω dv, (10)
and K k and fini are the Fourier transforms of K(x) and the initial distribution respectively. The time evolution of the kth Fourier component of the potential δφ(k, ω) is then given by an inverse Laplace transform. The asymptotic behavior of a time-dependent function is intimately related to the regularity of its Laplace transform: the more singular is the Laplace trasnform, the slower is the decay at large times. Thus to understand the large time behavior of δφ, we have to study the singularities of the D and G functions [START_REF] Weinberg | [END_REF], or rather the singularities of their analytic continuation. Although they are not eigenvalues of the linearized operator ( 6)- (7), roots of the continued D function yield an exponential damping of δφ: this is Landau damping.

We remark that an algebraic damping may happen if the stable stationary state f 0 or perturbation δf is not analytic [8]. In the next section, we show that the algebraic damping is actually unavoidable in the inhomogeneous case even if f 0 and δf are analytic.

Landau damping, inhomogeneous case

In astrophysics, the analog of Landau damping takes place, and is considered a standard mechanism in the relaxation of stellar structures, see [START_REF] Binney | Galactic Dynamics Second Edition[END_REF] for a textbook reference. However, it is remarkable that actual computations of Landau damping rates for self-gravitating structures are very rare; indeed, we are only aware of one instance of such a computation [START_REF] Weinberg | [END_REF]11] 1 . Clearly, the reason lies in the technical difficulty, as will become clear later. Inhomogeneous stationary states are also a common feature in plasma physics; in this context, Landau damping has been studied in the limit of vanishing velocity dispersion [12]. We feel that this situation calls for a more precise study of Landau damping close to inhomogeneous stationary states.

For simplicity, let us stick to a 1D setting with periodic boundary conditions, and consider an inhomogeneous stationary state of the Vlasov equation, of the form f 0 (x, v) = ϕ(h(x, v)).

In the inhomogeneous case, the derivative with respect to v survives in the linearized Vlasov equation, and this prevents us from obtaining an algebraic equation. A solution of this problem is to introduce the angle-action variables (θ, J) associated with the Hamiltonian h, which is always integrable, since we stick to 1D and f 0 (and thus h) is stationary. The angle-action variables eliminate the derivative with respect to J, and the Fourier transform with respect to θ yields an algebraic equation just as in the homogeneous case. However, in the inhomogeneous case, the Fourier bases for x, e ikx , and for θ, e imθ , are not orthogonal, and the algebraic equation is not diagonalized: all modes are mixed. Thus, in general, the algebraic equation is described by using an infinite size matrix D whose (l, m)-element is denoted by D(l, m, ω) and an infinite size vector G. The Fourier-Laplace component of potential is expressed in the form

δφ(l, ω) = m D -1 (l, m, ω)G(m, ω), (11)
where D -1 (l, m, ω) is the (l, m)-element of the inverse matrix D -1 . Each element of D and G has a similar form as in the homogeneous case, it reads as an integral:

F (ω) = ψ(J) Ω(J) -ω/m dJ (12)
where Ω(J) = (dh/dJ)(J). This solution "in matrix form" was introduced in the 70's [13,14]. Let us discuss the singularities of the function F , since we know they control the asymptotic behavior of the potential. Formally, the roots of the determinant of D correspond to poles of δφ, and these poles yield exponential dampings just as in the homogeneous case [15].

Poles are not the only singularities in the inhomogeneous case. Indeed the integration range for the action has a minimum value J v , whereas the integration range for the velocity in [START_REF] Weinberg | [END_REF] was R. As a result, a logarithmic singularity appears on the real axis of ω. A simple example helps to understand the creation mechanism of the new singularity. Suppose m = 1 and ψ(J) = 1, and Ω(J) is expanded around J v as

Ω(J) = Ω(J v ) + a(J -J v ) + • • • . ( 13
)
The F function is then expressed as

F (ω) = ∞ 0 dJ aJ -(ω -Ω(J v )) + • • • , ( 14
)
and it is clear that a logarithmic singularity appears at ω = ω v ≡ Ω(J v ). This singularity lies on the real axis, since the function Ω is real, and gives rise to a branch cut when viewed in the complex plane; hence the analytic continuation of the function F is multi-valued. In the homogeneous case, the integral runs over the whole real line, see [START_REF] Weinberg | [END_REF], hence the singularities on the real axis cancel: integrals on the left and the right of the singularity diverge in the same way, with opposite sign. Formally, the mechanism is that analytic continuation turns the integrals in the D and the G functions into principal values. Existence of the minimum point J v prevents us from applying the cancellation mechanism to the inhomogeneous case. We stress that this logarithmic singularity comes from the inhomogeneous nature of the stable stationary state f 0 . Thus, the singularity is unavoidable, and shows up even when the reference state f 0 and the initial perturbation are perfectly smooth.

Another type of singularity arises when Ω(J) has a local extremum. Generically, detailed computations reveal that the singular part of F is [16,17]

F sing (ω) = C 1 (ω -ω) ν ln |ω -ω v | + C 2 (ω -ω) µ H(ω -ω v ) (15)
where H(x) is the Heaviside step function, and a non-negative integer ν and a positive real value µ + 1 come from the leading powers in the expansion of ψ at J = J v . Known singularities of the Laplace transform δφ are summarized in Table 1.

Singularity Asymptotic Damping Homogeneous Inhomogeneous Pole Exponential Log. branch point Algebraic Heaviside step Algebraic

Table 1: Singularities of the Laplace transform of the potential for an analytic stable stationary state f 0 and perturbation δf with induced asymptotic dampings. The checkmarks signal the presence of the singularities.

The inverse Laplace transform of the logarithmic and the Heaviside step singularities in (15) both yield an oscillatory algebraic damping, with exponent 1+ν and 1+µ and frequency ω v respectively. The analysis shown in this section can be extended to higher dimensional systems [17]. The value of the exponent ν or µ giving the slowest damping can be obtained by a detailed analysis of the leading powers around the singular points on J in the reference state f 0 and perturbation δf .

From the discussion of this section, it is clear that the slowest decay is algebraic, and that it should rule the asymptotic behavior of a perturbation (assuming of course that the linearization remains valid at large times). Does it mean that exponential "Landau damping", coming from the roots of the determinant of the matrix D, is irrelevant? Of course not: it turns out that it often plays an important role at intermediate time scales, before algebraic damping dominates.

Some explicit computations on a simple model

The very small number of explicit computations of Landau damping rates for self gravitating structures is due to the inherent analytical and numerical difficulties of the needed analytic continuation; in [START_REF] Weinberg | [END_REF]11] this part is performed entirely numerically using Padé approximants, a procedure which is probably stable only close to the real axis and far from the branch cut singularities discussed above. If one accepts to give up the realism of 3D 1/|x| potential and turn to a toy model for gravitation, it becomes possible to perform explicit computations of the Landau poles. For this purpose, we will use in this section the Hamiltonian mean-field model [18]. The spatial domain is [0, 2π), and the interaction kernel is a cosine:

φ(x) = 2π 0 [1 -cos(x -y)]ρ(y)dy. (16)
Hence the self-consistent potential itself is of the form -M cos(x -ϑ). In a stationary state, M and ϑ are constant, thus the motion of each particle is that of a simple pendulum. It is well known that all trajectories can then be expressed in terms of elliptic functions and integrals [15], which can be readily analytically continued. Analytic continuation of the coefficient appearing in the matrix method is then relatively easy. Furthermore, thanks to the simple structure of the interaction kernel (only one Fourier mode is present), the determinant one has to compute is actually 2 × 2 and diagonal, rather than infinite in the general case.

Thanks to these simplifications, it becomes possible to explore all the "Landau poles" and follow their bifurcations when varying the stationary state. It is still a delicate task, as the landscape for δφ presents the full complexity of inhomogeneous stationary states: one should in principle explore an infinity of Riemann sheets. An example of the poles and branch cut singularities is given on Fig. 2.

Direct numerical simulations confirm that the dominant pole of the dispersion function (closest to the real axis) indeed dominate the dynamics over an intermediate time scale [15], and the asymptotic algebraic damping as shown on Fig. 3.

Conclusion

Although Landau damping around inhomogeneous stationary solutions is often considered as a direct analog of the much easier homogeneous case, we hope to have illustrated that spatial inhomogeneity requires a somewhat different mathematical structure, which in turn induces important physical differences. In the present short article, we have restricted to the linearized Vlasov equation, which is already very rich. The weakly non linear behavior of the Vlasov equation has also been the subject of many studies, old and new (for a few references [20,21,22,23,24,25]); a similar systematic study close to non homogeneous stationary state largely remains to be done, and the linear theory will be a necessary building block. In view of the differences showing up already at the linear level, surprises may be expected. 1) close to an inhomogeneous stationary state (see [16] for a precise definition of M (1) ). At small times, the amplitude oscillates and decreases exponentially (not easily seen in log-log scale): the dynamics is governed by the dominant Landau pole. At large times t 100, the algebraic asymptotic damping becomes clear. The simulations use a semi-Lagrangian method, (see [19] for details on the algorithm); the results are consistent with the N -body simulations of [16].

Figure 1 :

 1 Figure 1: Damping of an initial spatial structure by free transport (periodic spatial boundary condition are used). Left = initial condition; Center and Right: the shear in phase space creates finer and finer structures. After integration over the v variable, the spatial structure damps.

Figure 2 :

 2 Figure 2: Left panel: Contour plot of the modulus of the dispersion function Det(D) in the ω complex plane; the horizontal (vertical) axis is the real (imaginary) part of ω. The vertical dashed lines are branch cuts, separating different Riemann sheets; the contour levels are chosen so that the roots are clearly visible. Right panel: the reference stationary state ϕ around which the dispersion function shown on the left panel is computed. It is represented here as a function of the variable √ 2h, so that the statistical equilibrium distribution would be a gaussian.

Figure 3 :

 3 Figure3: Evolution of the perturbed field amplitude M(1) close to an inhomogeneous stationary state (see[16] for a precise definition of M(1) ). At small times, the amplitude oscillates and decreases exponentially (not easily seen in log-log scale): the dynamics is governed by the dominant Landau pole. At large times t 100, the algebraic asymptotic damping becomes clear. The simulations use a semi-Lagrangian method, (see[19] for details on the algorithm); the results are consistent with the N -body simulations of[16].

Computing instability rates is easier, since it does not require any analytic continuation

Acknowledgement

YYY acknowledges the support of a Grant-in-Aid for Scientific Research (C) 23560069.


Julien Barré
email: julien.barre@unice.fr

Shamik Gupta
email: shamikg1@gmail.com

Classical Heisenberg spins with long-range interactions: Relaxation to equilibrium for finite systems

Keywords: numbers: 05.20.Dd, 45.50.-j, 52.25.Dg Kinetic theory of gases and liquids, Metastable states

)] analyzed a model system comprising N globally coupled classical Heisenberg spins and evolving under classical spin dynamics. It was 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 leading order for the finite-N effects driving this relaxation. We demonstrate that corrections at this order are identically zero, so that relaxation occurs over a time longer than of order N , in agreement with the reported numerical results.

Introduction

Long-range interacting systems are characterized by an interparticle potential with a range that is of the order of the system size. In d dimensions, this corresponds to potentials decaying at large separation, r, as 1/r α , where α lies in the range 0 ≤ α ≤ d [1,2]. Examples include gravitational systems [3], plasmas [START_REF] Escande | Long-Range Interacting Systems[END_REF], two-dimensional hydrodynamics [START_REF] Bouchet | [END_REF], charged and dipolar systems [START_REF] Bramwell | Long-Range Interacting Systems[END_REF], 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 system of N particles proceeds in two steps: first, a collisionless relaxation, described by a Vlasovtype equation, brings the system close to a nonequilibrium state, called the "quasistationary state" (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-Gibbs equilibrium. This second step is usually described by a Lenard-Balescu-type equation [START_REF][END_REF]8]. This scenario is well established for plasmas and self-gravitating systems [START_REF] Binney | Galactic Dynamics: Second Edition[END_REF], and has been studied in detail in various toy models for long-range interactions [START_REF] Ruffo | [END_REF][11][12][13][14]. The QSS lifetime, which may be regarded to be of the same order of magnitude as the relaxation time, is thus an important quantity: knowing it allows to distinguish between non-relaxed systems, which should be described by a QSS, and relaxed ones, for which collisional effects need to be taken into account and an equilibrium description may be relevant. This lifetime depends on 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 selfgravitating systems [16]. Furthermore, it is known that the Lenard-Balescu collision term vanishes for 1d systems which do not develop any spatial inhomogeneity (see [17] for the 1d Coulomb case); one thus expects a relaxation time much longer than N in these cases. This is indeed numerically observed, see [12], where a time of order N 1.7 is reported, and [18], where larger systems sizes are studied and the relaxation time is claimed to be of order N 2 . A similar vanishing of the Lenard-Balescu operator has been found for point vortices in an axisymmetric configuration [19][20][21].

In a recent work, QSSs have been looked for and found in a dynamical setting different from the ones reviewed above, namely, in an anisotropic Heisenberg model with mean-field interactions [22]. Note that similar spin models with mean field interactions have been suggested to be relevant to describe some layered spin structures [23]. Specifically, the model in [22] comprises N globally coupled three-component Heisenberg spins evolving under classical spin dynamics. An associated Vlasov-type equation is introduced, and QSSs are stationary solutions of this equation. In addition, numerical simulations for axisymmetric QSS suggest a QSS lifetime increasing superlinearly with N. In order to understand this observation analytically, we present in this work a detailed study of the Lenard-Balescu operator 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 in the Lenard-Balescu operator; it also complicates the analytical structure of the dispersion relation for the Vlasov-type equation. Nevertheless, we can still demonstrate that corrections at 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 its equilibrium phase diagram. As a step towards deriving the Vlasov equation to analyze the evolution of the phase space distribution in the limit N → ∞, in section 3, we first write down the so-called Klimontovich equation; this leads to a derivation of the Vlasov equation and a discussion of a class of its stationary solutions in section 4. Section 5 contains our main results: it is devoted to the derivation of the Lenard-Balescu equation for our model of study, that is, the leading 1/N correction to the Vlasov equation; we show that this correction identically vanishes. In section 6, we consider an example of a Vlasov-stationary solution. In the energy range in which it is Vlasov stable, we demonstrate by performing numerical simulations of the dynamics that indeed its relaxation to equilibrium occurs over a timescale that does not grow linearly but rather superlinearly with N, in support of our analysis. The paper ends with conclusions.

The model

The model studied in Ref. [22] comprises N globally coupled classical Heisenberg spins of unit length, denoted by S i = (S ix , S iy , S iz ), i = 1, 2, . . . , N. In terms of spherical polar angles θ i ∈ [0, π] and φ i ∈ [0, 2π], one has S ix = sin θ i cos φ i , S iy = sin θ i sin φ i , S iz = cos θ i . The Hamiltonian of the system is

H = - J 2N N i,j=1 S i • S j + D N i=1 S 2 iz . (1)
Here, the first term with J > 0 describes a ferromagnetic mean-field coupling between the spins, while the second term is the energy due to a local anisotropy. We consider D > 0, for which the energy is lowered by having the magnetization

m ≡ 1 N N i=1 S i (2)
pointing in the xy plane. The coupling constant J in equation ( 1) is scaled by N to make the energy extensive [24], but the system is non-additive, implying thereby that it cannot be trivially subdivided into independent macroscopic parts, as is possible with short-range systems. 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 the energy 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 1 - 2 β c , (3)
where the inverse temperature β c satisfies

2 β c = 1 - 1 2β c D + e -βcD √ πβ c DErf[ √ β c D] . (4)
Here, Erf[x] = (2/ √ π)

x 0 dt e -t 2 is the error function. The microcanonical dynamics of the system (1) is given by the set of coupled first-order differential equations dS i dt = {S i , H}; i = 1, 2, . . . , N.

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} ≡ N i=1 (∂A/∂φ i ∂B/∂u i -∂A/∂u i ∂B/∂φ i ), which may be rewritten as [25] {A, B} =

N i=1 S i • ∂A ∂S i × ∂B ∂S i . (7)
Using equations ( 5) and [START_REF][END_REF], we obtain the equations of motion of the system as

Ṡix = S iy m z -S iz m y -2DS iy S iz , (8) Ṡiy
= S iz m x -S ix m z + 2DS ix S iz , (9) Ṡiz = S ix m y -S iy m x , (10)
where the dots denote derivative with respect to time. Summing over i in equation [START_REF] Ruffo | [END_REF], we find that m z is a constant of motion. The dynamics also conserves the total energy and the length of each spin. Using equations ( 8), [START_REF] Binney | Galactic Dynamics: Second Edition[END_REF], and [START_REF] Ruffo | [END_REF], we obtain the time evolution of the variables θ 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 + 2D cos θ i . ( 12
)

The Klimontovich equation

The state of the N-spin system is described by the discrete one-spin time-dependent density function

f d (u, φ, t) = 1 N N i=1 δ(u -u i (t))δ(φ -φ i (t)), (13)
which is defined such that f d (u, φ, t)dudφ counts the number of spins with its canonical coordinates in [u, u + du] and [φ, φ + dφ]. Here, δ is the Dirac delta function, (u, φ) are the Eulerian coordinates of the phase space, while (u i , φ i ) are the Lagrangian coordinates of the spins. Note that f d satisfies f d (u, φ, t) = f d (u, φ + 2π, t) and the normalization 2π 0 dφ 1 -1 du 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, φ) = √ 1 -u 2 (m y [f d ] cos φ -m x [f d ] sin φ), (15)
g φ ≡ g φ [f d ](u, φ) = m x [f d ] u √ 1 -u 2 cos φ + m y [f d ] u √ 1 -u 2 sin φ -m z [f d ] + 2Du, ( 16
) (m x , m y , m z )[f d ] = 2π 0 dφ 1 -1 du ( √ 1 -u 2 cos φ, √ 1 -u 2 sin φ, u)f d (u, φ, t). ( 17
)

The Vlasov equation, and a class of stationary solutions

We now define an averaged one-spin density function f 0 (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 0 (u, φ, t) + 1 √ N δf (u, φ, t), (18)
where, denoting by angular brackets the averaging with respect to the initial ensemble, we have f d = f 0 . Here, δf gives the difference between f d , which depends on the given initial condition, and f 0 , which depends on the average with respect to the ensemble of initial conditions.

Using equation (18) in equation ( 14), we get

∂f 0 ∂t + g 0 u ∂f 0 ∂u + g 0 φ ∂f 0 ∂φ + 1 N δg u ∂δf ∂u + δg φ ∂δf ∂φ + 1 √ N ∂δf ∂t + δg u ∂f 0 ∂u + g 0 u ∂δf ∂u + δg φ ∂f 0 ∂φ + g 0 φ ∂δf ∂φ = 0, ( 19
)
where

g 0 u = g u [f 0 ], g 0 φ = g φ [f 0 ], and δg u = √ 1 -u 2 (m y [δf ] cos φ -m x [δf ] sin φ), (20)
δg φ = m x [δf ] u √ 1 -u 2 cos φ + m y [δf ] u √ 1 -u 2 sin φ -m z [δf ]. (21)
We now average equation (19) with respect to the ensemble of initial conditions, and note that δf = 0 implies m x [δf ] = m y [δf ] = m z [δf ] = 0. Thus δg u = δg φ = 0, and we get

∂f 0 ∂t + g 0 u ∂f 0 ∂u + g 0 φ ∂f 0 ∂φ = - 1 N δg u ∂δf ∂u + δg φ ∂δf ∂φ . (22)
For finite times and in the limit N → ∞ (or, for times t ≪ N), we obtain the Vlasov equation satisfied by the averaged one-spin density function f 0 [START_REF] Nicholson | Introduction to Plasma Physics[END_REF] as

∂f 0 ∂t + g 0 u ∂f 0 ∂u + g 0 φ ∂f 0 ∂φ = 0. ( 23
)
Note that the Vlasov equation has been formally obtained after averaging over an ensemble of initial conditions. However, if the fluctuations in the initial conditions are weak and do not grow too fast in time, we expect the Vlasov equation to also describe the time evolution of a single initial condition in the limit N → ∞. This is put on firm mathematical grounds in [START_REF] Neunzert | [END_REF][28][29], for systems with a standard kinetic energy and a regular enough interaction potential.

From equations ( 15), (16), and (17), it is clear that any distribution that does not depend on the angle φ (thus, axisymmetric about the z-axis) is a stationary solution of the Vlasov equation (23).

The Lenard-Balescu equation

We obtain in this section the main result of our paper: for model (1), the Lenard-Balescu operator, computed for stationary solutions of the Vlasov equation of the form f 0 (u), identically vanishes.

Formal derivation

The Lenard-Balescu equation describes the slow evolution of a stable stationary solution of the Vlasov equation under the influence of finite-N corrections, at leading order in 1/N [START_REF] Nicholson | Introduction to Plasma Physics[END_REF]. From equation ( 22), we get the Lenard-Balescu equation as

∂f 0 ∂t = - 1 N δg u ∂δf ∂u + δg φ ∂δf ∂φ , (24)
where, subtracting equation ( 23) from equation ( 19), and keeping only the terms of order 1/ √ N, we find that δf follows the Vlasov equation linearized around its stable stationary solution f 0 :

∂δf ∂t = √ 1 -u 2 δm x sin φ -δm y cos φ ∂f 0 ∂u -(2Du -m z [f 0 ]) ∂δf ∂φ . (25)
Here, we have used

g 0 φ = 2Du -m z [f 0 ] and δm x,y ≡ m x,y [δf ]
. Now, a natural timescale separation hypothesis greatly reduces the complexity of finding the solutions of the coupled system 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 first solve equation (25), and then use its solution to compute the right-hand side of (24) in the limit t → ∞.

At this point, it is useful to make a comparison of our case with the standard case of particles with a kinetic energy moving in a classical potential, for example, a 1d system of particles with Coulomb interactions. The equivalent of the axisymmetric stationary solutions (independent of φ) introduced in section 4 are the homogeneous solutions which depend only on velocity in this standard setting, so that the analog of δg φ vanishes, whereas in our case of the spin dynamics, we have to deal with the extra term δg φ ∂δf ∂φ .

Solution of the linearized Vlasov equation

We now solve equation (25) for δf , using Fourier-Laplace transforms

δf (u, φ, t) = ∞ k=-∞ Γ dω 2π δf k (u, ω)e i(kφ-ωt) , (26)
δf k (u, ω) = 2π 0 dφ 2π ∞ 0 dt δf (u, φ, t)e -i(kφ-ωt) , (27)
where the Laplace contour Γ is a horizontal line in the complex-ω plane that passes above all singularities of δf k (u, ω).

We have

δm x = Γ dω 2π e -iωt 1 -1 du π √ 1 -u 2 δf -1 (u, ω) + δf +1 (u, ω) , (28)
δm y = Γ dω 2π e -iωt 1 -1 du π i √ 1 -u 2 δf -1 (u, ω) -δf +1 (u, ω) , (29)
δm z = m z [δf ] = Γ dω 2π e -iωt 1 -1 du πu δf 0 (u, ω), (30)
so that equation (25) gives

δf ±1 (u, ω) = - π √ 1 -u 2 f ′ 0 (u) 2Du -m z [f 0 ] ∓ ω 1 -1 du ′ √ 1 -u ′2 δf ±1 (u ′ , ω) ∓ iδf ±1 (u, 0) 2Du -m z [f 0 ] ∓ ω , (31)
where f ′ 0 (u) = ∂f 0 (u)/∂u and δf ±1 (u, 0) is the Fourier transform of the initial fluctuations δf (u, φ, 0). Multiplying both sides of the above equation by √ 1 -u 2 and then integrating over u, we get

ǫ ±1 (ω) 1 -1 du √ 1 -u 2 δf ±1 (u, ω) = ∓i 1 -1 du √ 1 -u 2 δf ±1 (u, 0) 2Du -m z [f 0 ] ∓ ω , (32)
where ǫ ±1 (ω) is the so-called "Plasma response dielectric function" [START_REF] Nicholson | Introduction to Plasma Physics[END_REF]:

ǫ ±1 (ω) = 1 + π LC du (1 -u 2 )f ′ 0 (u) 2Du -m z [f 0 ] ∓ ω . ( 33
)
To make the dielectric function ǫ +1 (also, ǫ -1 ) analytic in the vicinity of the real axis (Im(ω) = 0), which will be needed for later purpose, the above integral has to be performed along the Landau contour shown in Fig. 1, as discussed in [START_REF] Nicholson | Introduction to Plasma Physics[END_REF]; we have in this case

ǫ ±1 (ω) =                                  1 + π LC du (1-u 2 )f ′ 0 (u) 2Du-mz[f 0 ]∓ω ; (Im(ω) > 0), 1 + πP LC du (1-u 2 )f ′ 0 (u) 2Du-mz[f 0 ]∓ω ±i π 2 2D (1 -u 2 )f ′ 0 (u)| (mz [f 0 ]±ω)/(2D) ; (Im(ω) = 0), 1 + π LC du (1-u 2 )f ′ 0 (u) 2Du-mz[f 0 ]∓ω ±i 2π 2 2D (1 -u 2 )f ′ 0 (u)| (mz [f 0 ]±ω)/(2D) ; (Im(ω) < 0), ( 34
)
where P denotes the principal part. 33) already defines an analytic function ǫ +1 (respectively. ǫ -1 ) in the vicinity of a real ω, without the need to take into account extra pole contributions as in (34). Note that ǫ ±1 (ω) has two branch cut singularities on the real axis at ω = 2D -m z [f 0 ] and ω = -2D -m z [f 0 ], which is the reason why these functions may be seen as multi-valued in the lower-half ω-plane.

If ω / ∈ [-2D -m z [f 0 ]; 2D -m z [f 0 ]] (respectively, ω / ∈ [-2D + m z [f 0 ]; 2D + m z [f 0 ]]), equation (
Im(u) Re(u) Im(u) Re(u) Im(u) Re(u) Im(u) 0 0 0 Im(ω) < 0, -2D -m z [f 0 ] < Re(ω) < 2D -m z [f 0 ] Im(ω) = 0, -2D -m z [f 0 ] < Re(ω) < 2D -m z [f 0 ] -1 -1 -1 1 1 1 Im(ω) > 0, -2D -m z [f 0 ] < Re(ω) < 2D -m z [f 0 ]
Figure 1. The Landau contour LC to evaluate ǫ +1 (ω), shown as the dashed red line in the complex-u plane.

Using equation (32) in equation (31) gives

δf ±1 (u, ω) = ± iπ √ 1 -u 2 f ′ 0 (u) ǫ ± (ω)[2Du -m z [f 0 ] ∓ ω] 1 -1 du √ 1 -u 2 δf ±1 (u, φ, 0) 2Du -m z [f 0 ] ∓ ω ∓ iδf ±1 (u, 0) 2Du -m z [f 0 ] ∓ ω . ( 35
)
We see from the above expression that the real pole at ω = ±(2Du -m z [f 0 ]) is due to the free part of the evolution that does not involve interaction among the spins, and results in undamped oscillations of the fluctuations δf (u, φ, t), see equation [START_REF] Nicholson | Introduction to Plasma Physics[END_REF]. The other set of poles corresponds to the zeros of the dielectric function ǫ ± (ω), i.e., values ω p (complex in general) that satisfy

ǫ ±1 (ω p ) = 0. ( 36
)
Equation [START_REF] Nicholson | Introduction to Plasma Physics[END_REF] 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, when the poles lie in the upper-half complex ω-plane, the fluctuations grow in time. On the other hand, when the poles are either on or below the real-ω axis, the fluctuations do not grow in time, but rather oscillate or decay in time, respectively. Then, the condition ensuring linear stability of a stationary solution of the Vlasov equation reads

ǫ ±1 (ω p ) = 0 ⇒ Im(ω p ) ≤ 0. ( 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

δm + (ω) = 2iπ ǫ -1 (ω) 1 -1 du √ 1 -u 2 δf -1 (u, 0) 2Du -m z [f 0 ] + ω , (39)
δm

-(ω) = - 2iπ ǫ +1 (ω) 1 -1 du √ 1 -u 2 δf +1 (u, 0) 2Du -m z [f 0 ] -ω . (40)
Equation ( 35) may now be expressed in terms of δm ± as

δf +1 (u, ω) = - √ 1 -u 2 f ′ 0 (u) 2[2Du -m z [f 0 ] -ω] δm -(ω) - iδf +1 (u, 0) 2Du -m z [f 0 ] -ω , (41)
δf -1 (u, ω) = - √ 1 -u 2 f ′ 0 (u) 2[2Du -m z [f 0 ] + ω] δm + (ω) + iδf -1 (u, 0) 2Du -m z [f 0 ] + ω . ( 42
)

Computing the Lenard-Balescu operator

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

δg φ ∂δf ∂φ = k,l Γ Γ ′ dωdω ′ 4π 2 e i(k+l)φ e -i(ω+ω ′ )t il δg φ,k (u, ω) δf l (u, ω ′ ) , ( 43
)
δg u ∂δf ∂u = k,l Γ Γ ′ dωdω ′ 4π 2 e i(k+l)φ e -i(ω+ω ′ )t δg u,k (u, ω) ∂ δf l (u, ω ′ ) ∂u . (44)
From equations ( 20) and ( 21), we have

δg φ = δm + + δm - 2 cos φ u √ 1 -u 2 + δm + -δm - 2i sin φ u √ 1 -u 2 -δm z [δf ], (45)
δg u = √ 1 -u 2 δm + -δm - 2i cos φ - δm + + δm - 2 sin φ , (46)
so that we have

δg φ,±1 (u, ω) = δm ∓ (ω) u 2 √ 1 -u 2 , δg φ,0 (u, ω) = δm z (ω), ( 47
)
δg u,±1 (u, ω) = ∓ √ 1 -u 2 2i δm ∓ (ω). ( 48
)
From equations ( 39) and (40), we see that δm + (respectively, δm -) depends on δf -1 (u, 0) (respectively, δf +1 (u, 0)). Then, to compute the right hand hand side of equations ( 43) and (44), we need to evaluate averages of the type δf k (u, 0)δf l (v, 0) for the initial fluctuations at t = 0. Note that t = 0 as a notation is somewhat inappropriate, since actually this computation 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). In this case, one gets

δf k (u, 0)δf l (v, 0) = δ k,-l 2π [f 0 (u)δ(u -v) + h(u, v)] , (49)
where δ k,-l is the Kronecker delta function, and h(u, v) is a smooth function. The precise form of this undetermined smooth function will play no role in the computation, as we will show below.

Computing δg φ ∂δf ∂φ

Combining equations ( 47) and (49), we see that on the right hand side of equation (43), only the terms (k = +1, l = -1) and (k = -1, l = +1) give a non-zero contribution. We detail below the computation for the (k = +1, l = -1) case, the other being similar. Note that δg φ ∂δf ∂φ is real; thus, we have to compute only the real part of the (k = +1, l = -1) term, since its imaginary part must cancel with that of the (k = -1, l = +1) term. In the following computation, we set h(u, v) = 0; we will check at the end that indeed the contributions containing h vanish.

We have

δg φ,+1 (u, ω) δf -1 (u, ω ′ ) = - uf ′ 0 (u) 4(2Du -m z [f 0 ] + ω ′ ) δm -(ω) δm + (ω ′ ) + i u 2 √ 1 -u 2 (2Du -m z [f 0 ] + ω ′ ) δm -(ω)δf -1 (u, 0) . ( 50
)
We thus need

δm -(ω) δm + (ω ′ ) = 2π ǫ +1 (ω)ǫ -1 (ω ′ ) 1 -1 dv (1 -v 2 )f 0 (v) (2Dv -m z [f 0 ] -ω)(2Dv -m z [f 0 ] + ω ′ ) , (51)
δm -(ω)δf -1 (u, 0) = -i √ 1 -u 2 f 0 (u) ǫ +1 (ω)(2Du -m z [f 0 ] -ω) , (52)
where we have used equation (49). Using these equations in equation (50), we obtain

δg φ,+1 (u, ω) δf -1 (u, ω ′ ) = - πuf ′ 0 (u) 2(2Du -m z [f 0 ] + ω ′ )ǫ +1 (ω)ǫ -1 (ω ′ ) × 1 -1 dv (1 -v 2 )f 0 (v) (2Dv -m z [f 0 ] -ω)(2Dv -m z [f 0 ] + ω ′ ) + uf 0 (u) 2(2Du -m z [f 0 ] + ω ′ )(2Du -m z [f 0 ] -ω)ǫ +1 (ω) ≡ b 1 + b 2 . ( 53
)
Now, to compute the contribution to δg φ ∂δf ∂φ

for (k = +1, l = -1), we have to integrate b 1 and b 2 over ω and ω ′ ; we define

B 1 ≡ Γ Γ ′ dωdω ′ 4π 2 (-i)e -i(ω+ω ′ )t b 1 , B 2 ≡ Γ Γ ′ dωdω ′ 4π 2 (-i)e -i(ω+ω ′ )t b 2 . ( 54
)
We want to compute B 1 and B 2 in the limit t → ∞; thus, we will discard all terms decaying for large t.

Let us start with B 1 . We have First, note that for large t,

B 1 = i πuf ′ 0 (u) 2 Γ ′ dω ′ 2π e -iω ′ t (2Du -m z [f 0 ] + ω ′ )(2Dv -m z [f 0 ] + ω ′ )ǫ -1 (ω ′ ) × Γ dω 2π e -iωt 1 -1 dv (1 -v 2 )f 0 (v) ǫ +1 (ω)(2Dv -m z [f 0 ] -ω) . ( 55
) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 Im(ω) Zeros of ǫ +1 (ω) 2D -m z [f 0 ] -iδ -2D -m z [f 0 ] -iδ Re(ω) Γ 2Dv -m z [f 0 ] -iδ ′
Γ dω 2π e -iωt ǫ +1 (ω)(2Dv -m z [f 0 ] -ω) ∼ i e -i(2Dv-mz [f 0 ])t ǫ +1 (2Dv -m z [f 0 ]) , (56)
which is obtained by deforming the Γ contour into the Im(ω) < 0 half-plane, see Fig. 2, and noting that due to the e -iωt factor, the only contribution to the contour integral in the limit t → ∞ comes from the pole 2Dv -m z [f 0 ] -iδ ′ . Indeed, since f 0 is assumed stable, the zeros of ǫ +1 , if any, have a negative imaginary part; their contributions thus decay exponentially in time. Similarly, the branch-cut singularities, see Fig. 2, contribute terms decaying algebraically in time. Note also that the important singularity at 2Dv

-m z [f 0 ] is real, and is in the range [-2D -m z [f 0 ]; 2D -m z [f 0 ]].
Then, according to the discussion in section 5.2, one needs to use the expression (34) for ǫ +1 . The same remark applies to all the computations below, and we will not recall it each time.

The integration over ω ′ in equation ( 55) is a bit more complicated, since there are two poles on the real axis, at

ω ′ = -2Du + m z [f 0 ] and ω ′ = -2Dv + m z [f 0 ]
. With the same method as for ω, we obtain

Γ ′ dω ′ 2π e -iω ′ t (2Du -m z [f 0 ] + ω ′ )(2Dv -m z [f 0 ] + ω ′ )ǫ -1 (ω ′ ) = -i e i(2Du-mz [f 0 ])t 2D(v -u)ǫ -1 (-2Du + m z [f 0 ]) - e i(2Dv-mz [f 0 ])t 2D(v -u)ǫ -1 (-2Dv + m z [f 0 ]) . ( 57
)
We finally get

B 1 = i πuf ′ 0 (u) 2 1 -1 dv (1 -v 2 )f 0 (v) ǫ +1 (2Dv -m z [f 0 ]) × e -i2D(v-u)t 2D(v -u)ǫ -1 (-2Du + m z [f 0 ]) - 1 2D(v -u)ǫ -1 (-2Dv + m z [f 0 ]) . ( 58
)
To perform the integral over v, we will use the following lemma.

Lemma:

lim t→∞ 1 -1 dv ϕ(v) v -u e -i2D(v-u)t ǫ -1 (-2Du + m 0 ) - 1 ǫ -1 (-2Dv + m 0 ) = -iπ ϕ(u) ǫ -1 (-2Du + m 0 )
.( 59)

Proof:

1 -1 dv ϕ(v) v -u e -i2D(v-u)t ǫ -1 (-2Du + m 0 ) - 1 ǫ -1 (-2Dv + m 0 ) = 2D(1-u)t -2D(1+u)t dx ϕ(u + x 2Dt ) x e -ix ǫ -1 (-2Du + m 0 ) - 1 ǫ -1 (-2Du + m 0 -x t )
.( 60)

Recalling that -1 < u < 1, and taking the limit t → ∞, the above expression simplifies to

∞ -∞ dx ϕ(u) x e -ix ǫ -1 (-2Du + m 0 ) - 1 ǫ -1 (-2Du + m 0 ) = ϕ(u) ǫ -1 (-2Du + m 0 ) ∞ -∞ dx e -ix -1 x . ( 61
)
The last integral is -iπ, which completes the proof of the lemma.

Using the lemma, we conclude that

B 1 = π 2 uf ′ 0 (u)(1 -u 2 )f 0 (u) 4Dǫ +1 (2Du -m z [f 0 ])ǫ -1 (-2Du + m z [f 0 ]) . ( 62
)
We now turn to the computation of B 2 . The integration over ω and ω ′ is performed as above, deforming the contours in the lower-half ω ′ -plane, and keeping only the contributions of the poles on the real axis. We obtain

B 2 = -i uf 0 (u) 2 1 ǫ +1 (2Du -m z [f 0 ]) = -i uf 0 (u) 2 ǫ -1 (-2Du + m z [f 0 ]) ǫ +1 (2Du -m z [f 0 ])ǫ -1 (-2Du + m z [f 0 ]) . (63)
As explained above, we need to compute only the real part of B 2 ; thus, we keep only the contribution coming from the imaginary part of ǫ -1

(-2Du + m z [f 0 ]). Using from equation (34) that Im[ǫ -1 (-2Du + m z [f 0 ])] = - π 2 2D (1 -u 2 )f ′ 0 (u), (64)
we conclude that

Re(B 2 ) = - π 2 uf 0 (u)(1 -u 2 )f ′ 0 (u) 4Dǫ +1 (2Du -m z [f 0 ])ǫ -1 (-2Du + m z [f 0 ]) . (65)
From equations (62) and (65), we find that B 1 + Re(B 2 ) = 0. Similar to above, one can show that the real part of the contribution to δg φ ∂δf ∂φ from (k = -1, l = +1) also vanishes, while, as discussed above, the imaginary part of the contribution to δg φ ∂δf ∂φ from (k = -1, l = +1) must cancel that from (k = +1, l = -1). So, we conclude that for our model,

δg φ ∂δf ∂φ = 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 2 . First, its contribution to δm -(ω)δf -1 (u, 0) is

- i ǫ +1 (ω) 1 -1 dv √ 1 -v 2 2Dv -m z [f 0 ] -ω h(v, u). (67)
Thus, its contribution to B 2 is

u 2 √ 1 -u 2 1 -1 dv √ 1 -v 2 h(v, u) Γ Γ ′ dωdω ′ 4π 2 e -i(ω+ω ′ )t 1 ǫ +1 (ω)(2Dv -m z [f 0 ] -ω) × 1 (2Du -m z [f 0 ] + ω ′ ) . (68)
The integrals over ω and ω ′ can be performed as before. Since the poles for ω and ω ′ are different, we see that for large t, a factor e -2Di(v-u)t oscillating rapidly in time remains in the integral over v; this leads to this integral vanishing in the limit t → ∞. A similar phenomenon 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.

Computing δg u ∂δf ∂u

As for δg φ ∂δf ∂φ , we see that on the right hand side of equation (44), only the terms (k = +1, l = -1) and (k = -1, l = +1) are non-zero. We give below the computation for the case (k = +1, l = -1) case, the other being similar. Again, we can restrict the computations to the real part of each term, since δg u ∂δf ∂u is real. From equations ( 42) and (48), we get

δg u,+1 (u, ω) ∂ δf -1 (u, ω ′ ) ∂u = √ 1 -u 2 4i ∂ ∂u √ 1 -u 2 f ′ 0 (u) [2Du -m z [f 0 ] + ω ′ ] δm -(ω) δm + (ω ′ ) - 1 2 √ 1 -u 2 ∂ ∂u 1 2Du -m z [f 0 ] + ω ′ δm -(ω)δf -1 (u, 0) = π 2i √ 1 -u 2 ǫ +1 (ω)ǫ -1 (ω ′ ) ∂ ∂u √ 1 -u 2 f ′ 0 (u) [2Du -m z [f 0 ] + ω ′ ] × 1 -1 dv (1 -v 2 )f 0 (v) (2Dv -m z [f 0 ] -ω)(2Dv -m z [f 0 ] + ω ′ ) + i 2 √ 1 -u 2 ǫ +1 (ω) ∂ ∂u √ 1 -u 2 f 0 (u) (2Du -m z [f 0 ] + ω ′ )(2Du -m z [f 0 ] -ω) ≡ a 1 + a 2 , (69)
where, in obtaining the second equality, we have used equations (51) and (52). Now, to compute the contribution to δg u ∂δf ∂u

for (k = +1, l = -1), we have to integrate a 1 and a 2 over ω and ω ′ . Let us define

A 1 ≡ Γ Γ ′ dωdω ′ 4π 2 e -i(ω+ω ′ )t a 1 , A 2 ≡ Γ Γ ′ dωdω ′ 4π 2 e -i(ω+ω ′ )t a 2 . ( 70
)
We want to compute A 1 and A 2 in the limit t → ∞, so that we may discard all terms decaying for large t.

Let us first compute A 1 . Integration over ω and ω ′ may be carried out by deforming the contours Γ and Γ ′ , as done in the preceding subsection. One gets

Γ dω 2π e -iωt ǫ +1 (ω)(2Dv -m z [f 0 ] -ω) = i e -i(2Dv-mz [f 0 ])t ǫ +1 (2Dv -m z [f 0 ]) , (71)
and

Γ ′ dω ′ 2π e -iω ′ t (2Du -m z [f 0 ] + ω ′ )(2Dv -m z [f 0 ] + ω ′ )ǫ -1 (ω ′ ) = -i e i(2Du-mz[f 0 ])t 2D(v -u)ǫ -1 (-2Du + m z [f 0 ]) - e i(2Dv-mz [f 0 ])t 2D(v -u)ǫ -1 (-2Dv + m z [f 0 ]) , (72)
so that in the limit t → ∞, we have

A 1 = -i π 2 √ 1 -u 2 ∂ ∂u √ 1 -u 2 f ′ 0 (u) 1 -1 dv (1 -v 2 )f 0 (v) 2D(v -u)ǫ +1 (2Dv -m z [f 0 ]) × e -i2D(v-u)t ǫ -1 (-2Du + m z [f 0 ]) - 1 ǫ -1 (-2Dv + m z [f 0 ]) = - π 2 4D √ 1 -u 2 ∂ ∂u √ 1 -u 2 (1 -u 2 )f 0 (u)f ′ 0 (u) ǫ +1 (2Du -m z [f 0 ])ǫ -1 (-2Du + m z [f 0 ]) , (73)
where, in obtaining the second equality, we have used the lemma (59). Now, A 2 may be computed along the same lines as done for B 2 in the preceding subsection. One gets

Re(A 2 ) = i 2 √ 1 -u 2 ∂ ∂u √ 1 -u 2 f 0 (u) ǫ +1 (2Du -m z [f 0 ] = π 2 4D √ 1 -u 2 ∂ ∂u √ 1 -u 2 (1 -u 2 )f 0 (u)f ′ 0 (u) ǫ +1 (2Du -m z [f 0 ])ǫ -1 (-2Du + m z [f 0 ]) , (74)
where, in obtaining the second equality, we have used equation (64). From equations ( 73) and (74), we see that A 1 + Re(A 2 ) = 0. Similarly, one can show that the real part of the contribution to δg u ∂δf ∂u from (k = -1, l = +1) also vanishes, while the imaginary part of the contribution to δg u ∂δf ∂u from (k = -1, l = +1) must cancel that from (k = +1, l = -1). So, we conclude that for our model,

δg u ∂δf ∂u = 0. ( 75
)
Combining equations (66) and (75), we see that the Lenard-Balescu operator identically vanishes for our model, which is the announced result.

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 0 a state prepared by sampling independently for each of the N spins the angle φ uniformly over [0, 2π] and the angle θ uniformly over an interval of length (a + b) asymmetric about

θ = π/2, that is, θ ∈ [π/2 -a : π/2 + b]. The corresponding single-spin distribution is f 0 (u, φ) = 1 2π p(u), (76)
with p(u), the distribution for u, given by

p(u) =     
Since m x [f 0 ] = m y [f 0 ] = 0, we have g 0 u = 0, and also, ∂f 0 /∂φ = 0; it then follows that the state (76) is stationary under the Vlasov dynamics (23). Note that we have

m z [f 0 ] = sin a -sin b 2 . ( 79
)
As mentioned after equation (37), the condition Im(ω p ) = 0 will correspond to the marginal stability of the state (76), so that the zeros ω p of the dielectric function lie on the real-ω axis. Let us denote these zeros as ω * pr . From equation (34), we find that ω * pr satisfies

1 + πP 1 -1 du (1 -u 2 )f ′ 0 (u) 2Du -m z [f 0 ] ∓ ω * pr ± iπ 2 (1 -u 2 )f ′ 0 (u) (mz [f 0 ]±ω * pr )/(2D) = 0. ( 80
)
Equating the real and the imaginary parts to zero, we get

1 + πP 1 -1 du (1 -u 2 )f ′ 0 (u) 2Du -m z [f 0 ] ∓ ω * pr = 0, (81)
(1 -u 2 )f ′ 0 (u) (mz[f 0 ]±ω * pr )/(2D) = 0. (82)
Now, equation (76) gives

f ′ 0 (u) = 1 2π(sin a + sin b) [δ(u + sin a) -δ(u -sin b)], (83)
so that we obtain from equation (82) that

δ m z [f 0 ] ± ω * pr 2D + sin a = δ m z [f 0 ] ± ω * pr 2D -sin a , (84)
implying that

ω * pr = -m z [f 0 ]. ( 85
for 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. On the basis of our analysis in this paper showing the Lenard-Balescu operator being identically zero, we expect that for finite N, the state relaxes to Boltzmann-Gibbs equilibrium on a timescale ∼ N δ , with δ > 1. This was indeed observed in Ref. [22] for the class of initial states (76) that is non-magnetized, that is, a = b; in this case, combining equations ( 86) and (87), we get

e * | a=b = D 3 + 12D . ( 88
)
For energies e < e * | a=b , the state being linearly unstable under the Vlasov dynamics was seen to relax for finite N to the Boltzmann-Gibbs equilibrium state over a timescale ∼ ln N [22]. , where e c can be computed from equation (3) to be e c ≈ 0.258. We expect on the basis of the analysis presented in this paper that in this energy range, when the state (76) is Vlasov-stable, the relaxation to equilibrium should occur over a timescale that scales superlinearly with N. For e = 0.22, results of numerical simulations of the dynamics shown in Fig. 3 indeed suggest a relaxation timescale ∼ N δ , with δ > 1; for the range of system sizes explored in this numerical experiment, any value of δ between about 1.7 and 2 is compatible with the data.

Conclusions

In this paper, we have shown that the Lenard-Balescu operator identically vanishes for a system of globally coupled anisotropic Heisenberg spins, in an axially symmetric Vlasovstable state. This result explains the numerical findings of [22], reporting a relaxation time for this system that scales superlinearly with N. To our knowledge, it is the first time that this kind of results has been obtained for a spin dynamics. This raises further questions, e.g., what are the general conditions to ensure that the Lenard-Balescu operator vanishes? The classical explanation relies on the structure of resonances between the particle trajectories: in the absence of resonances between particles with different momentum, the Lenard-Balescu operator should vanish. This heuristic argument applies to systems of particles moving in a 1d position space, thus with a 2d phase space, when the system is homogeneous [17,30], implying a relaxation time growing superlinearly with N. This is also the case for axisymmetric configurations of point vortices [19][20][21], where the phase space is again two-dimensional. In a similar manner, it can be argued for the model we have studied that spins with different projections on the z-axis cannot exchange energy because they cannot be in resonance. In a sense, our precise computations validate this qualitative picture. However, recent numerical simulations of a model with a 4d phase space have also shown a relaxation time that appears superlinear in N over the range of system sizes studied [14]: one would expect resonances to appear in this case. Thus, understanding the general conditions under which the Lenard-Balescu operator vanishes may still remain a partly open question.

One may also wonder how the relaxation occurs when the Lenard-Balescu operator vanishes. Formally, the Klimontovich expansion suggests that the next leading term is of order 1/N 2 . Although writing down this term is possible in principle, its evaluation is difficult. 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 equation approximately describes a single initial condition. Putting this on firm mathematical grounds is an outstanding question, on which some preliminary progress has been made recently [31].

Figure 2 .

 2 Figure 2. The Γ contour to evaluate the integral in equation (56). Here, δ → 0 and δ ′ → 0 are infinitesimal positive numbers. The thick vertical lines denote branch-cut singularities.

1 sin

 1 a+sin b if u ∈ [-sin b, sin a],

)A = 2 (

 2 Using equation (85) in equation (81), we get 4D = cos 2 a sin b + cos 2 b sin a sin a sin b (sin a + sin b) , (86) which when combined with equation (78) gives the energy e = e * = A B ; sin a -sin b) cos 2 a + 2 cot a sin 2 b cos a + 2 cos b cot b sin 2 a + 2 cos 2 b (sin b -sin a) -3(sin a -sin b) 2 (sin a + sin b), B = 24(sin a + sin b),

Figure 3 .

 3 Figure 3. Magnetization m(t) as a function of tN -2 in the energy range in which the state (76) with a = 0.1, b = 0.3 is linearly stable under the Vlasov dynamics. Here, e = 0.22, D = 8.224, while the system sizes are marked in the figure. The figure suggests a relaxation timescale ∼ N δ with δ = 2, but the data do not allow a precise determination of the exponent δ: any value of δ between about 1.7 and 2 is compatible with the data.

Acknowledgements

SG acknowledges the support of the Indo-French Centre for the Promotion of Advanced Research under Project 4604-3 and the hospitality of Laboratoire J. A. Dieudonné, Université