W Ettoumi
email: wahb.ettoumi@ens-cachan.fr

M.-C Firpo

Linear theory and violent relaxation in long-range systems: a test case

Keywords:

In this article, several aspects of the dynamics of a toy model for longrange Hamiltonian systems are tackled focusing on linearly unstable unmagnetized (i.e. force-free) cold equilibria states of the Hamiltonian Mean Field (HMF). For special cases, exact finite-N linear growth rates have been exhibited, including, in some spatially inhomogeneous case, finite-N corrections. A random matrix approach is then proposed to estimate the finite-N growth rate for some random initial states. Within the continuous, N → ∞, approach, the growth rates are finally derived without restricting to spatially homogeneous cases. Then, these linear results are used to discuss the large-time nonlinear evolution. A simple criterion is proposed to measure the ability of the system to undergo a violent relaxation that transports the mean field modulus in the vicinity of its equilibrium value within some linear e-folding times.

Introduction

Systems of particles interacting via two-body long-range forces are well-known to have peculiar equilibrium and non-equilibrium statistical mechanics (see e.g. Ref. [1] and references therein). As far as their relaxation properties are concerned, much progress originated from numerical simulations of the one-dimensional gravitational system. In 1982, Wright, Miller and Stein [2] observed its reluctance to thermalize due to the existence of quasi-stationary states (QSSs). These observations were later refined by various authors showing that the relaxation of the one-dimensional gravitational system usually proceeds through a rapid approach to a QSS, referred to as violent relaxation, followed by a very slow drift toward equilibrium [3,4]. Such studies initiated a still very active line of research (see e.g. Refs. [5,6,7]) on the intricate interplay between dynamics, ergodic properties and statistical mechanics in self-gravitating Hamiltonian systems.

Moreover, in the case where space dimension is larger than one, the thoroughly investigated gravitational system, as well as the Coulomb system, combine the difficulties of long-range interaction with a short-range divergence. This was a motivation to introduce models in which the potential was truncated to retain only its long-range components. In addition, the periodic boundary conditions considered in such models amount to work with a compact space which is numerically convenient. Various numerical simulations and theoretical arguments [1,8,9,10] gave indications that the corresponding toy models obtained in this way were sharing purely long-range relaxation characteristics similar to the original systems.

The Hamiltonian Mean Field (HMF) model [9] derives from such a truncation procedure as it amounts, in its attractive ferromagnetic-like form, to the one-dimensional gravitational system with periodic boundary conditions where only the lowest Fourier mode is retained. It has become a well-known toy model to address the intricate relationships between dynamics and statistical mechanics of long-range interacting systems. It is defined by the following Hamiltonian

H = N i=1 p i 2 2 + 1 2N N i=1 N j=1 [1 -cos (θ i -θ j )] , (1)
where N is the number of particles, and θ i and p i denote respectively the position and momentum of the i th particle. A useful collective quantity to introduce is the so-called magnetization vector (M x , M y ) with

M x = 1 N N i=1 cos θ i and M y = 1 N N i=1 sin θ i (2)
The average energy per particle U = H/N reads then

U = N i=1 p i 2 2N + 1 2 1 -M 2 , (3)
where M ≡ M x 2 + M y 2 denotes the modulus of the magnetization vector.

Recently, much interest has been devoted to the QSSs which are known to be responsible for the very slow convergence towards the statistical mechanics equilibrium predictions. Far from being difficult to generate, these QSSs naturally emerge in the HMF model from waterbag initial distributions (see e.g. [11,12,14,15,16] and the recent review [17]). It is also known that initial waterbag conditions in momenta, associated to zero or almost zero initial magnetization, induce the longest lasting QSSs. However, when lowering towards zero the initial temperature of the particles, it is possible to exhibit waterbag momenta configurations with vanishing magnetization in which the magnetization eventually converges exponentially towards its Boltzmann-Gibbs equilibrium value. This calls for a linear theory approach.

Linear stability of the HMF model about unmagnetized equilibrium states has been up to now only studied within the Vlasov framework [9,12,28], which assumes in particular an infinite number of particles. Moreover, let alone some very recent publications [28], the linear stability of spatially inhomogeneous, unmagnetized, equilibria has never been considered yet.

The motivation of the present study is then twofold: Firstly and mostly, one wishes to tackle the linear study of the unmagnetized cold HMF equilibria, within a finite-N , therefore exact, framework; secondly, the ensuing nonlinear dynamics is briefly addressed to show that the thermalization of cold unmagnetized HMF systems finely illustrates Lynden-Bell's concept of violent relaxation for long-range systems.

In Section 2, we shall establish the finite-N framework used for the linear stability derivation. In Section 3, we shall calculate the exact linear growth rates for two finite-N equilibria, both of zero temperature and zero magnetization, and compare them to numerical simulations. Section 4 is dedicated to a random matrix approach for the calculation of symmetric non-deterministic initial states growth rates. In Section 5, we eventually derive the linear theory in the N → ∞ limit using a fluid approach derived from the Vlasov equation for a vanishing temperature. Section 6 ends this study by discussing the connections between the linear features just derived and the HMF thermalization properties. The dynamics of the cold unmagnetized HMF model is proposed as a paradigm of violent relaxation.

Linear dynamics about cold unmagnetized finite-N equilibria

The equations of motion can straightforwardly be written from Equation (1) as

∀ k ∈ {1, ..., N } ,    θk = p k ṗk = 1 N N i=1 sin (θ i -θ k ) ≡ F k (4)
Using Equation ( 2), the force acting on the particle k may be written as

F k = M y cos θ k -M x sin θ k . (5)
Let us consider unmagnetized finite-N equilibria, namely stationary states of the equations of motion (4), with M x = M y = 0. This amounts to have p * k = 0 and angles θ * k distributed in such a way that M x = M y = 0. Let us remark that, since the total momentum P = N i=1 p i is a constant of motion, the zero momentum equilibrium case considered here is just the cold (i.e. monokinetic) case for the special choice P = 0.

Let us perform the linear stability of this system. In that purpose, we write θ k = θ * k + δθ k and p k = δp k where the asterisk denotes the unperturbed solution. At first order in δθ k , the force δF k felt by the particle k verifies

δF k = 1 N cos θ * k N i=1 cos θ * i δθ i + sin θ * k N i=1 sin θ * i δθ i . (6)
This yields the following linear system

  δ θ k δ p k   = 0 N I N A 0 N {δθ k } {δp k } (7)
where I N is the N × N identity matrix and A a N × N matrix defined by

A i,j = 1 N cos(θ * i -θ * j ). (8)
The stability depends on the eigenvalues {λ k } of the Jacobian matrix in Equation (7).

Let us name it J. The eigenvalue problem can now be reduced to the unique study of A using the following transformation

J -λ I 2N = -λI N 0 N A I N • I N -1 λ I N 0 N -λI N + 1 λ A (9)
provided λ = 0. Hence,

det(J-λI 2N ) = det(-λI N ) det(-λI N + 1 λ A) = (-1) N det A -λ 2 I N .( 10
)
In other words, writing χ M the characteristic polynomial of M , one has

χ J (λ) = (-1) N χ A (λ 2 ) (11)
This factorization allows us to focus only on the matrix A, and deduce the eigenvalues of the higher-order matrix J by taking the square root of A ones. We shall now study the case of two particular finite-N equilibria that advantageously simplify A.

3. Exact finite-N treatment for two special cold force-free equilibria

The quiet start case

Having in mind the computational plasma terminology, we define the so-called "quiet start" configuration as the equilibrium characterized by an equipartition of the particles on the circle. It is here formally described by ∀k, θ * k = 2πk/N , up to some constant phase, and p * k = 0. Using this definition in Equation ( 8), the matrix A immediately reduces to

A i,j = 1 N cos 2π N (i -j) . (12)
We notice that we can rewrite the coefficients as

A i,j = A (i-j)modN , (13)
which makes evident that A is a circulant matrix. Using Equation (11), it can be easily shown that its eigenvalues {λ 2 k } are expressed as

λ 2 k = N j=1
A j e 2ijkπ/N . (

Moreover, since A is real and symmetric, its eigenvalues are real, and one can identify the previous equation with its real part yielding

λ 2 k = 1 N N j=1 cos 2jπ N cos 2jkπ N . (15)
This is just

λ 2 k = 1 2 (δ k,1 + δ k,N -1 ) , (16)
so that A has only one double non-zero eigenvalue equal to 1/2. Using Equation (11) finally yields the expected growth rate γ QS coming from the spatially homogeneous Vlasov linear theory [9] for the cold waterbag as

γ QS = λ 2 1 = 1 √ 2 , (17)
with no finite-N correction. In order to test the validity of this linear study, we performed numerical simulations based on a fourth-order symplectic integrator [18]. Starting from a quiet start configuration, every particle is moved by a uniformly randomized quantity ǫ ≪ 2π/N . As shown in Figure 1, the behaviour of the system during the early times shows a very good agreement with the predicted exponential growth, and does not depend on the number of particles, which only changes the initial magnetization resulting from the perturbation. The slow start that one can notice in Figure 1 is due to the fact that the initial random configuration does not belong to a pure eigenmode, so that the contribution of the other eigenvalues takes some time to become negligible in front of the growing one. The same phenomenon can be viewed the instantaneous growth rates displayed in Figure 6.

Bi-clustered quiet start

The previous case involved a finite-N analog of a homogeneous Vlasov force-free equilibrium. We can construct another finite-N equilibrium with zero magnetization by uniformly distributing N/2 particles in a cluster of size ∆θ centered on a given position, and by settling the positions of the N/2 remaining ones by rotating the first cluster by π with

∀ k ∈ 1, ..., N 2 ,      θ * k = -∆θ + 4k∆θ N θ * N/2+k = θ * k + π p * k = 0 (18)
Figure 2 shows an example of such a bi-cluster configuration. The same stability analysis as in the homogeneous quiet start case can be performed. However, one expects now the growth rate to depend on ∆θ, the clusters' size. Equation ( 11) is still valid, but the matrix A is no longer circulant as it was in the simple quiet start case. A can be rewritten under the form

A = L -L -L L , ( 19
)
where L is a N/2 × N/2 matrix with coefficients

L ij = 1 N cos 4∆θ N (i -j) . (20)
Therefore, the characteristic polynomial of A reads det A -

λ 2 I N = (-2λ 2 ) N/2 det L - λ 2 2 I N/2 . (21)
Or, equivalently,

χ A (λ 2 ) = (-2λ 2 ) N/2 χ L λ 2 2 ( 22
)
This decomposition allows us to focus on the smaller matrix L. Unfortunately, L is not circulant either, but is a Toeplitz matrix. Indeed, one can write L ij = L |i-j| . A work performed by Treichler [19] showed that for a m × m Toeplitz matrix generated from the coefficients t k = cos(kω), the only two non-zero eigenvalues are

ν± (m, ω) = 1 2 m ± sin (mω) sin(ω) . ( 23
)
Equation ( 23) allows us to obtain the eigenvalues {ν k } of the matrix

L ν ± = 1 N ν± N 2 , 4∆θ N = 1 4 ± sin(2∆θ) 2N sin (4∆θ/N ) . (24)
Equations ( 11), ( 22) and ( 24) give the growth rate as

γ BCQS = 2ν + = 1 2 + sin(2∆θ) N sin(4∆θ/N ) , (25)
that is, in the large N limit,

γ BCQS = 1 √ 2 1 + sin(2∆θ) 2∆θ , (26)
up to O(N -2 ) terms. Figure 3 shows the comparison between the numerically computed growth rates and the theoretical prediction of Equation (25). When ∆θ = π/2, particles are uniformly distributed on the circle as in the previous equilibrium and the growth rate is exactly 1/ √ 2 with no finite-N correction. Otherwise, the growth rate depends on the number of particles but we checked that the difference with the asymptotic result ( 26) is already very small for N above 10 particles.

Extension to symmetric random initial configurations

In this section, we will only consider symmetric equilibria prepared to give M = 0 in the following way: we distribute N/2 particles at random on a partition of a given length and put the N/2 remaining particles by shifting the random ones by π. Using the appropriate indexation of particles (see (18)), the calculation of growth rates amounts to determine the largest eigenvalue of the N/2 × N/2 matrix L defined by

∀(i, j) ∈ 1, ..., N 2
2 , L ij = cos (θ i -θ j ) , (27)
where the particle positions {θ i } are distributed according to some f 0 (θ), and where the 1/N normalization factor has been voluntarily omitted so that the coefficients' distribution does not depend on the number of particles. We used a method based on Random Matrix Theory to calculate the growth rate's expectation. When the random coefficients verify L ij f 0 = µ > 0, an extension of Wigner's law [20] states that the largest eigenvalue is asymptotically approximated by

ν = 2 N N/2 i,j L ij + σ 2 µ + o 1 √ N , (28)
where [20] states that ν has a normal distribution of expectation 1 + (N/2 -1) µ + σ 2 /µ and bounded finite variance 2σ 2 . Using Equations ( 11) and ( 22), one finds the mean squared growth rate as

σ 2 = L 2 ij f 0 . Since L ii f 0 = 1, Theorem 2 in Reference
λ 2 = 2 N ν = 2 N 1 + N 2 -1 µ + σ 2 µ . ( 29
)
The value of λ 2 is hence distributed according to a normal law of variance 8σ 2 /N 2 because of the 2/N rescaling factor. The expectation of the growth rate reads then

γ = N -1 1 0 √ x exp - N 2 16σ 2 x -λ 2 2 dx, ( 30
)
where N is the normalization factor. This approach simplifies the calculation of the growth rate, since no more effort in the diagonalization of matrix A has to be done. However, the most restrictive applicability condition of Equation ( 30) is the particular symmetry of the initial state, which allows the use of Equation [START_REF] White | Gravitational Dynamics[END_REF]. The conditions on the probability distribution of L ij are not as limiting as the required symmetry of the initial state.

Random uniform bi-cluster

The equilibrium configuration generated by a waterbag distribution of parameter ∆θ yields a state topologically close to Eqs. (18). Therefore, we expect the growth rates for these random states to be close to the ones given in Equation (25). However, the latter Equation cannot be the asymptotic form of Equation ( 30) since the highest eigenvalue fluctuates around 1 + (N/2 -1) µ + σ 2 /µ, which corresponds to the eigenvalue of the deterministic matrix a ij = µ for i = j, a ii = 1, that is completely different from the deterministic matrix L defined by Equation (20).

One needs to verify the applicability of Equation ( 28) before calculating the expectation of the growth rate with Equation (30). µ and σ 2 are accessible through the following formulae

µ = cos (θ i -θ j ) f 0 = f 0 (θ i )f 0 (θ j ) cos (θ i -θ j ) dθ i dθ j , (31)
σ 2 + µ 2 = cos 2 (θ i -θ j ) f 0 = f 0 (θ i )f 0 (θ j ) cos 2 (θ i -θ j ) dθ i dθ j . ( 32
)
Writing χ(X) the characteristic function of the set X, the waterbag distribution of parameter ∆θ has the probability density

f 0 (θ) = 1 2∆θ χ([-∆θ, ∆θ]). ( 33
)
Hence, we have

µ = sinc 2 (∆θ) , (34)
σ 2 = 1 2 + sin 2 (2∆θ) 16∆θ 2 -µ 2 . ( 35
)
Clearly, µ > 0 for ∆θ < π/2 and σ 2 is finite, which means that Equation (30) holds. Figure 4 shows the numerical fits for randomized bi-clustered initial configurations and the theoretical mean value given by Equation (30). We also plotted the growth rate (25) Plot of the growth rates from the uniformly randomized bicluster initial configuration with respect to ∆θ. The circles correspond to averaged exponential fits from the numerically integrated magnetization, the solid curve corresponds to Equation (30). The dashed curve corresponds to Equation (25). Each numerical calculation has been performed with 1000 particles, and each circle results from the average of 8 runs. We clearly see a very good agreement with the theoretical prediction.

corresponding to the deterministic bicluster to show that, as expected, its behaviour is close to the expectation of the random one for a wide range of ∆θ.

Random Gaussian bi-cluster

In this subsection we show another example of use of Equation ( 30) for a more difficult case. The particles are no longer distributed according to a waterbag density, but with a Gaussian one. We define f 0 (θ) by

f 0 (θ) = σ θ √ 2π erf π 2σ θ √ 2 -1 exp - θ 2 2σ 2 θ
As shown on Figure 5, the agreement between the experimental average and the random matrix theory is very good. Plot of the growth rates as a function of σ θ for the random Gaussian bicluster equilibria. The circles correspond to averaged exponential fits from the numerically integrated magnetization, while the curve corresponds to Equation (30). Each numerical calculation has been performed with 1000 particles, and each circle results from the average of 8 runs. Here again, the agreement with the random matrix theory prediction is very good.

In the general case, relaxing the symmetry assumption in the preparation of the finite-N equilibria yields a vanishing average value of the matrix elements of A consistently with a vanishing magnetization. This prevents the application of the Theorem 2 of [20] and Equation (30) no longer holds. To overtake this difficulty, we now turn to the more usual continuous approach.

Linear theory within the continuous approach

Writing f (θ, p, t) the distribution function, one has the formal Vlasov equation ∂f ∂t

+ p ∂f ∂θ + E(θ, t) ∂f ∂p = 0. ( 39
)
One can define the density of particles n(θ, t) and the mean velocity field v(θ, t) through

n(θ, t) = +∞ -∞ f (θ, p, t)dp and v(θ, t) = 1 n(θ, t) +∞ -∞ p f (θ, p, t)dp.(40)
The HMF force field E(θ, t) is then given by

E(θ, t) = π -π sin (α -θ) n(α, t)dα.
(41)

The early evolution of the cold HMF can be reduced to a fluid description. The hierarchy of the moments of the Vlasov equation can be stopped at the first order since the temperature vanishes. By taking the moment of order zero, one immediately obtains

∂n ∂t + ∂(nv) ∂θ = 0. ( 42
)
Multiplying Equation ( 39) by p and integrating over p, one has

∂(nv) ∂t + +∞ -∞ p 2 ∂f ∂θ dp -E(θ, t)nv = 0. ( 43
)
Injecting Equation (42) in Equation ( 43) and dividing by n yields

∂v ∂t + v ∂v ∂θ -E(θ, t) = 0, ( 44
)
in the cold case, for which the mean square of the momentum equals the square of the mean velocity. Considering the stationary solution given by some n 0 (θ) yielding a zero magnetization, i.e. having a zero m = 1 Fourier component, and v = 0, one puts v = δv(θ) exp(iωt) and n = n 0 (θ) + δn(θ) exp(iωt) in Equations ( 41), ( 42

This is

ωδn k + k m n 0,m δv k-m = 0, ( 45
) iωδv ±1 = ± π ω δn ±1 , (46)
and δv m = 0 for m = ±1. This gives finally

ω 2 δn k + kπ (n 0,k-1 δn 1 -n 0,k+1 δn -1 ) = 0. ( 47
)
The dispersion relation is thus given by det M (ω) = 0, where M is generally an infinite matrix of elements

M kℓ = ω 2 δ k,ℓ + kπ (δ ℓ,1 n 0,k-1 -δ ℓ,-1 n 0,k+1 ) , ( 48
) with (k, ℓ) ∈ Z 2 .
Let us consider the finite-size 2N +1 square matrix M (N ) , which coefficients coincide with

M ij , ∀(i, j) ∈ [[-N, N ]] 2 . For N = 1, we have M (1) =    ω 2 + πn 0,0 0 -πn 0,-2 0 ω 2 0 -πn 0,2 0 ω 2 + πn 0,0    . ( 49
)
The condition det M = 0 is then fulfilled when

ω 2 + πn 0,0 2 = π 2 n 0,2 n 0,-2 , ( 50
)
but since n(θ, t) is real, n 0,-2 = n 0,2 where the bar denotes the complex conjugate. Therefore, one can take the square root in Equation ( 50), leading to

ω 2 ± = -πn 0,0 ± π |n 0,2 | . ( 51
)
It is then easy to show by recurrence that ∀N, det M (N +1) = ω 4 det M (N ) , so that Equation (51) gives the only non-vanishing roots to the general dispersion relation. This allows to obtain the growth rate for the cold unmagnetized HMF in the infinite N limit as

γ = 1 + 2π|n 0,2 | √ 2 , ( 52
)
where we used the fact that n 0 (θ) is normalized giving n 0,0 = 1/2π. It is easy to check that the N → ∞ growth rate ( 26) may be obtained from Equation (52) in the uniform bicluster configuration, for which 2πn 0,2m = sinc(2m∆θ) and n 0,2m+1 = 0, ∀m ∈ Z. The formula (52) was successfully tested for a variety of spatially inhomogeneous equilibria with M = 0 (See the plots of the instantaneous growth rates Ṁ /M with respect to time in Figure 6).

Growth rate

n0 ∝ 1 + cos(2θ)

√ 3/2 1/ √ 2 n0 ∝ 1 + cos(3θ) n0 ∝ 1 + cos(4θ) n0 ∝ 1 + cos(5θ) Figure 6.
Plot of the numerical growth rate Ṁ /M with respect to time for different initial equilibria. When the magnetization is almost exponential, Ṁ /M is practically constant and matches the theoretical value given by Equation (52), represented by the horizontal lines. We clearly see that when the second harmonic of n 0 is zero, the growth rates are equal to 1/ √ 2 as in the homogeneous case. All the runs have been performed using 10 4 particles.

We shall now briefly discuss the large-time nonlinear features of the HMF model in the light of Lynden-Bell's picture of violent relaxation.

6. Final discussion on the cold HMF case : an example of violent relaxation

The concept of violent relaxation

The regularity of the observed luminosity profiles of elliptical galaxies suggests that they have reached some equilibrium states. As the two-body, collisional, characteristic timescale is much larger than the estimated ages of galaxies, the concept of "violent relaxation" was introduced by Lynden-Bell in his famous 1967 paper [21,[START_REF] White | Gravitational Dynamics[END_REF] as a collisionless scenario to account for the rapid evolution of the galaxies toward quasistationary or quasi-equilibrium states. On the basis of statistical arguments, Lynden-Bell gave an expression for the coarse grained distribution function of galaxies in a meta-equilibrium state. Several studies gave evidence of discrepancies between the numerically obtained stationary distributions and Lynden-Bell predictions [START_REF] Severne | [END_REF]24] for the self-gravitating system, which may result from insufficient mixing properties for the application of this statistical theory. As far as the HMF model is concerned, a recent application of Lynden-Bell's theory [25] gave however satisfactory quantitative predictions. In any case, that is in spite of the controversies related to the strict application of Lynden-Bell statistics, there appears to be nowadays a common agreement on the scenario of the relaxation process for N -body long-range systems. It is supposed to be divided in two parts: the properly speaking violent relaxation part, namely a rapid, collisionless, evolution of the system towards a quasi-stationary metastable state, followed by a much slower thermalization phase towards the state predicted by equilibrium statistical mechanics (see e.g. Ref. [26] for recent related results on generic long-range systems). The HMF model provides a simple, yet nontrivial, long-range system to consider this issue.

Linear instability and violent relaxation in the HMF model

Extensive numerical studies of the HMF model are now available in the literature, mostly for two types of initial states: waterbag initial distribution functions in positions and momenta and Maxwellian initial distribution functions in momenta with possibly waterbag initial distribution functions for positions. Long-time discrepancies between time averaged observables and their ensemble predictions have been reported only with initial waterbag distributions in momenta, in connection with the emergence of so-called quasi-stationary states (see e.g. [11,12,13]).

Linear theory, that has been up to now mostly formulated within the Vlasov framework, and for spatially homogeneous states -with the noticeable exception of very recent extensions to inhomogeneous states [27,28] -, can be used as a guideline to discuss relaxation properties in the spirit of Lynden-Bell's picture. For instance, following the derivation of Vlasov linear theory given in Ref. [9], the spatially homogeneous cases with initial waterbag distributions in momenta are unstable with growth rates equal to 1/2 -3T where T is the initial temperature associated to the waterbag. Consequently, for T > 1/6 the system becomes linearly stable. This case corresponds to an energy density U equal to 7/12 [9,12]. This value happens to coincide with the energy threshold value above which pathological relaxation behaviours have been reported for spatially homogeneous, waterbag in momenta, initial distributions functions.

In the following, we shall discuss the conditions under which the modulus of the mean-field of the linearly unstable HMF model may saturate nonlinearly at a value that is close, but yet different, to its ensemble prediction.

Nonlinear saturation of the mean-field

Let us consider the fluid model introduced in Section 5. The validity of the zero temperature approximation is limited to the initial stage of the instability but this is not critical here as our point is just to establish the threshold condition at which nonlinear effects come into play to stop the growth of the initially vanishing magnetization.

In terms of Fourier components, Equation (44) reads

dδv m (t) dt + i ℓ ℓδv m-ℓ (t)δv ℓ (t) = E m (t), (53)
with E m (t) = 0 for m = ±1, and

E 1 (t) = E * -1 (t) = i (M x -iM y ) /2.
Nonlinear effects are non longer negligible when the nonlinear contribution balances the other terms, namely when the nonlinear mode couplings term balances the linear term. On m = 1, this yields

dδv 1 (t) dt ∼ ℓ ℓδv 1-ℓ (t)δv ℓ (t) ∼ δv -1 δv 2 (54) with, on m = 2, dδv 2 (t) dt = -i ℓ ℓδv 2-ℓ (t)δv ℓ (t) ∼ -iδv 2 1 (55)
Eq. ( 55) translates the fact that the m = 2 mode is nonlinearly triggered, at about twice the linear growth rate, and will be the first mode to emerge from the otherwise essentially m = 1 instability. The m = 0 velocity perturbation remains identically constant. We have

δv 2 (t) ≃ -i t 0 δv 1 (s) 2 ds ≃ -i δv 2 1 (0) 2γ exp (2γt)
so that, replacing this in Eq. ( 54), the nonlinear saturation takes place when

2γ 2 ∼ δv 2 1 ≃ δv 2 1 (0) exp (2γt) (56)
together with the linear balance coming from Equation (53),

γδv 1 ∼ E 1 = i 2 (M x -iM y ) . (57)
Eventually the last two orderings (56) and (57) give the order of the modulus of the magnetization at the nonlinear saturation, M sat , as

M sat ∼ 2 3/2 γ 2 . ( 58
)
This threshold is qualitative. In order to obtain a more quantitatively valid estimate, which is not an easy task at all, one should have to take into account, in particular, the fact that the growth rate does not remain equal to its linear value up to nonlinear saturation. Let us use the common plasma physics terminology for wave-particle interaction and introduce the trapping time of the particles in the mean-field potential well, namely the inverse of the bounce frequency ω b = √ M sat at the nonlinear saturation. This is just the characteristic timescale for particles oscillating in the HMF potential well, and does not depend on N . This is a so-called nonlinear timescale as the mean-field is supposed to be initially vanishingly small. The linear timescale is obviously given by the e-folding time, namely by the inverse of the linear growth rate γ. Then Equation (58) translates the fact that nonlinear saturation takes place when both timescales balance, namely for γ ∼ ω b . As a consequence of Equation (58), one expects nonlinear saturation to take place close to equilibrium predictions provided that the linear growth rate is of the order of the ensemble average of the nonlinear frequency, namely provided that

γ ∼ √ M µ . (59)
On the contrary, if the linear growth rate γ is small enough so that γ ≪ √ M µ , the value of the nonlinear magnetization threshold will be significantly below the value predicted by equilibrium statistical mechanics.

Application to the cold beam case

The cold case, on which we have just focused, illustrates well this scenario. Actually, large-time simulations show that the modulus of the magnetization, for instance, quickly converges towards a saturated state with a value that is close to its ensemble prediction. Let us check that the condition (59) is indeed satisfied for the cold unmagnetized HMF. This corresponds to an energy density U = 1/2. According to equilibrium statistical dynamics [29], this gives M µ = M c = 1/ √ β ≡ √ T where the inverse of the temperature is given implicitly by I 1 ( √ β)/I 0 ( √ β) = 1/ √ β. Numerically, the ensemble average of the magnetization for the cold HMF is then M c ≃ 0.62 and the equilibrium temperature is T ≃ 0.39. Since γ = 1/ √ 2, one effectively has γ ∼ ω b µ , where we have used ensemble equivalence.

It is interesting to note that similar observations have been reported for another long-range, mean-field, wave-particle Hamiltonian model starting with a cold beam of particles [30].

The relaxation process

Figure 7a shows that the cold unmagnetized HMF does actually experience an initial violent relaxation phase. The figure makes however apparent that this phase is followed by a much slower thermalization phase that is needed for phase space sweeping and complete convergence towards equilibrium statistical predictions. Actually, in the early nonlinear saturation, the magnetization oscillates about 0.60, which is below the ensemble equilibrium value. The drift towards the equilibrium statistical predictions takes place on a much longer timescale than the violent relaxation timescale.

More precisely, putting δM 0 the initial infinitesimal perturbation of the modulus of the magnetization, that may contain some N -dependence for the finite-N HMF model, the time ∆t needed to reach the nonlinear saturation threshold is of the order of γ -1 ln(γ 2 /δM 0 ). Assuming that δM 0 ∝ N -1/2 , one gets ∆t ∝ log N , which is confirmed by the plot on Figure 7b. The ensuing timescale needed for complete thermalization of the system has a much stronger N -dependence, as it is basically a binary collisional timescale that diverges with N . Therefore, in the Vlasov N → ∞ limit, the system would have been trapped in the QSS corresponding to nonlinear saturation. The specific thermalization process will be more closely examined in a forthcoming study.

Figure 1 .

 1 Figure 1. Numerical plot of the magnetization with respect to time. As predicted by Equation (17), the slope does not depend on the number of particles N , and is in very good agreement with the theoretical value of 1/ √ 2, plotted in thin dashed line.

Figure 2 .

 2 Figure2. Plot of some finite-N bi-clustered equilibrium where each particle has a zero momentum and faces its symmetric on the circle, providing a zero magnetization.

Figure 3 .

 3 Figure3. Plot of the growth rate as a function of ∆θ for slightly perturbed bicluster initial configurations. The circles correspond to exponential fits from the numerically integrated magnetization, and the solid curve corresponds to Equation(25). Each numerical calculation has been performed with 1000 particles.

Figure 4 .

 4 Figure 4. Plot of the growth rates from the uniformly randomized bicluster initial configuration with respect to ∆θ. The circles correspond to averaged exponential fits from the numerically integrated magnetization, the solid curve corresponds to Equation(30). The dashed curve corresponds to Equation(25). Each numerical calculation has been performed with 1000 particles, and each circle results from the average of 8 runs. We clearly see a very good agreement with the theoretical prediction.

Figure 5 .

 5 Figure 5. Plot of the growth rates as a function of σ θ for the random Gaussian bicluster equilibria. The circles correspond to averaged exponential fits from the numerically integrated magnetization, while the curve corresponds to Equation(30). Each numerical calculation has been performed with 1000 particles, and each circle results from the average of 8 runs. Here again, the agreement with the random matrix theory prediction is very good.

  ) and (44). Expanding δv(θ) and δn(θ) in Fourier series, one obtains the linear system iω m δn m exp(imθ) + i m ℓ (m + ℓ) n 0,m δv ℓ exp [i (m + ℓ) θ] = 0, iω m δv m exp(imθ) + iπ m δn m δ -1,m e -iθδ 1,m e iθ = 0.

Figure 7 .

 7 Figure 7. (a) Time evolution of the magnetization starting from a gently perturbed homogeneous quiet start of the cold HMF. (b) N dependency of ∆t, defined as the time needed for the magnetization to first reach the value 0.7 that roughly corresponds to nonlinear saturation.