The quasilinear theory in the approach of long-range systems to quasi-stationary states

# The quasilinear theory in the approach of long-range systems to quasi-stationary states

Alessandro Campa and Pierre-Henri Chavanis Complex Systems and Theoretical Physics Unit, Health and Technology Department, Istituto Superiore di Sanità, and INFN Roma1, Gruppo Collegato Sanità, 00161 Roma, Italy
Laboratoire de Physique Théorique, Université Paul Sabatier, 118 route de Narbonne 31062 Toulouse, France
###### Abstract

We develop a quasilinear theory of the Vlasov equation in order to describe the approach of systems with long-range interactions to quasi-stationary states. We derive a diffusion equation governing the evolution of the velocity distribution of the system towards a steady state. This steady state is expected to correspond to the angle-averaged quasi-stationary distribution function reached by the Vlasov equation as a result of a violent relaxation. We compare the prediction of the quasilinear theory to direct numerical simulations of the Hamiltonian Mean Field model, starting from an unstable spatially homogeneous distribution, either Gaussian or semi-elliptical. In the Gaussian case, we find that the quasilinear theory works reasonably well for weakly unstable initial conditions (i.e. close to the critical energy ) and that it is able to predict the energy marking the out-of-equilibrium phase transition between unmagnetized and magnetized quasi-stationary states. Similarly, the quasilinear theory works well for energies close to the instability threshold of the semi-elliptical case , and it predicts the out-of-equilibrium transition at . In both situations, at energies lower than the out-of-equilibrium transition the quasilinear theory works less well, the disagreement with the numerical simulations increasing by decreasing the energy. In that case, we observe, in agreement with our previous numerical study [A. Campa and P.-H. Chavanis, Eur. Phys. J. B 86, 170 (2013)], that the quasi-stationary states are remarkably well fitted by polytropic distributions (Tsallis distributions) with index (Gaussian case) or (semi-elliptical case). In particular, these polytropic distributions are able to account for the region of negative specific heats in the out-of-equilibrium caloric curve, unlike the Boltzmann and Lynden-Bell distributions.

## I Introduction

Systems with long-range interactions exhibit peculiar properties both at equilibrium and out-of-equilibrium. Among the equilibrium features that are not present in short-range systems, probably the most striking are the ensemble inequivalence and the negative specific heat in the microcanonical ensemble. Even richer is the phenomenology one observes out-of-equilibrium: long-lived quasi-stationary states (QSSs), non-mixing behavior, phase transitions between non-Boltzmannian distributions, etc houches (); assisebook (); campabook (); physrep (). Many results have been established, or illustrated, in the context of a toy model called the Hamiltonian Mean Field (HMF) model hmf (). Although progress has constantly been made in the last couple of decades, there are still points very poorly understood. One of these issues concerns the approach to QSSs.

In the statistical dynamics study of long-range systems, one of the main tools is the Vlasov equation that represents the interaction through a mean-field. For this reason the Vlasov equation describes a “collisionless” dynamics. It can be shown that this is a very good approximation that improves more and more when the number of components increases (the Vlasov equation becomes essentially exact when the number of particles ). More precisely, the larger the system the larger the time range for which the true dynamics can be studied with the Vlasov equation. Eventually, “collisional” effects, i.e. finite effects, will take over, driving the system towards Boltzmann-Gibbs equilibrium. However, in some cases, especially in the astrophysical context, the time scale of these “collisional” effects can be orders of magnitudes larger that the age of the Universe. It is therefore important to characterize the various stages of the Vlasov dynamics. This equation admits both an infinite number of conserved quantitites (the Casimirs) and an infinite number of stationary states. The QSS reached by the system after an initial transient is just one of these stationary states, obviously one which is stable with respect to perturbations. The infinite number of stationary states of the Vlasov equation makes the prediction of which one is selected by the system very difficult ybbdr (); incomplete (). Actually, it is fair to say that, theoretically, this is still a completely unsolved problem.

The initial state is in general a non-stationary, or unstable stationary, state of the Vlasov equation. In that case, the system evolves very rapidly, with the dynamics governed by the Vlasov equation, towards a QSS. Simulations show that, as one could argue on the basis of the existence of an infinite number of stationary states, the selected QSS strongly depends on the initial state of the system. A predictive theory was built by Lynden-Bell lb () in an astrophysical context. According to this theory, the QSS reached by the system is the one that maximizes an entropy functional while preserving the energy and all the Casimirs. Although perfectly defined, this problem is amenable to a computable solution only when the distribution of the initial state has a small number of phase levels, i.e., it is piecewise constant. However, even when the computation is possible, the comparison with the simulations shows that the prediction of the Lynden-Bell theory works only in some cases. The explanation of its failure111Although we shall focus in this paper on situations where the Lynden-Bell theory does not provide a good description of the QSSs, we would like to emphasize that there exist situations where the Lynden-Bell theory works remarkably well and is able to predict the existence of out-of-equilibrium phase transitions and re-entrant phases that would not have been detected without the help of this theory; see in particular Refs. precommun (); stan () and Fig. 36 of Ref. epjb2013 () for successful predictions of the Lynden-Bell theory. has to be sought in the absence of a complete mixing of the dynamics (under given constraints), a property that is implicitly assumed in the theory incomplete (). Even worse, in self-gravitating systems, e.g. elliptical galaxies, the Lynden-Bell distribution has infinite mass, thus it cannot represent the actual QSS of the system.

In the attempt to find alternatives, we have proposed epjb2013 () to compare the QSS found in numerical simulations of the HMF model with polytropic distributions (sometimes called Tsallis distributions tsallisbook ()). These distributions are critical points of functionals of the one-particle distribution function that generalize the usual Boltzmann entropy. This by itself does not give to polytropes a privileged status with respect to other possible stationary states of the Vlasov equation. The reasons to employ them can be summarized as follows. They have been introduced long ago in astrophysics, where they are called stellar polytropes, in the attempt to describe the non-Boltzmannian distributions observed btnew (). Being the critical points (at fixed particle number and energy) of functionals of the form , where is a convex function, they determine distributions of the form , where is the individual energy, with thlb (). The stability of these Vlasov stationary states has to be determined in each case (see, e.g., the comments in this respect in Ref. jsm2010 (); epjb2010 ()). The function corresponding to the Boltzmann entropy is . Under suitable assumptions the function that gives rise to polytropic distributions can be obtained from a modification of the Lynden-Bell theory that takes into account the absence of complete mixing assise (); epjb2010 (). In a numerical study epjb2013 (), we have shown that in some cases in which the prediction of the Lynden-Bell theory fails to obtain the distribution of the QSSs reached by the system after the violent relaxation, polytropic distributions are a good approximation of the QSSs. However, we have found that the quality of the approximation worsens when the initial distribution is only weakly unstable.

In the case of weak instability of initially homogeneous distributions, another type of approximation could be more suitable, i.e. the quasilinear (QL) approximation.222The QL theory is well-known in plasma physics for the Coulombian interaction nicholson (). We apply it here to a new situation, namely the HMF model, in which the interaction between particles is attractive. This approach is based on the assumption that, although the initial distribution is not Vlasov stable, nevertheless its evolution towards a Vlasov stable stationary state is such that it is all the time only slightly inhomogeneous. The assumption will be made more precise below. In this paper we have tested the validity of the QL theory in the HMF model, motivated by the fact that the simulations performed in epjb2013 () to study the polytropic approximation have shown that, when the initial distribution is only weakly unstable and the polytropic approximation is not good, the QSS is only weakly inhomogeneous. Therefore, we make the hypothesis that the distribution is only slightly inhomogeneous all the time and compare the prediction of the QL theory to the results of numerical simulations.

In Section II we give the derivation of the QL approximation for a generic mean field system composed of particles moving on a circle. In Section III, we apply the theory to the HMF model, where the particular form of the interaction allows a simplification of the expressions. In Sections IV, V and VI we show the results for a Gaussian and semi-elliptical initial condition and determine the domain of validity of the QL theory and the polytropic fit.

## Ii The QL approximation

We describe here the derivation of the QL approximation for the evolution of the one-particle distribution function as determined from the Vlasov equation. In view of the application to the HMF model, we consider a generic one dimensional system on the circle, i.e., a system with Hamiltonian

 H=N∑i=1pi2+1N∑∑1≤i

where , and the potential , which has periodicity, is supposed to be continuous (the mass of the particles has been put equal to without loss of generality). The Vlasov equation for the one-particle distribution function is

 ∂f(θ,p,t)∂t+p∂f(θ,p,t)∂θ−∂Φ(θ,t;f)∂θ∂f(θ,p,t)∂p=0, (2)

where the mean field potential is given by

 Φ(θ,t;f)=∫dp′∫2π0dθ′V(θ−θ′)f(θ′,p′,t). (3)

For large long-range systems with the Vlasov equation is a very good approximation for times of order houches (); assisebook (); campabook (); physrep (). Functions that depend only on are particular stationary solution of the Vlasov equation.333It is easy to see that the mean field potential is constant in for a uniform distribution, i.e., a function that does not depend on . The stability of these stationary solutions is studied by putting and linearizing the Vlasov equation around , i.e., keeping only the terms at most linear in . We then obtain the following linear equation for :

 ∂f1(θ,p,t)∂t+p∂f1(θ,p,t)∂θ−∂Φ(θ,t;f1)∂θ∂f0(p)∂p=0, (4)

which is valid as long as . This is the linearized Vlasov equation. The proper frequencies of the dynamics determined by this equation can be studied by inserting solutions where the time dependence is of the form . The function is stable if all the solution for have a non positive imaginary part. In that case, the distribution function can change only because of the effects of “collisions” (finite effects) that are not considered in the Vlasov equation physrep (); epjp (). We recall that in the case of one dimensional systems the time range of validity of the Vlasov equation grows like (for inhomogeneous systems) or with a higher power of (for homogeneous systems) physrep (); epjp (). If, on the contrary, some of the frequencies have a positive imaginary part, is Vlasov unstable, and the distribution will depart from rapidly.

There is an approximation that can work if the positive imaginary parts of the proper frequencies are small enough (the scale with which to measure this smallness will be precised later): the so-called QL approximation. We first derive the corresponding equations, then we describe their solutions, and finally we explain in which case one can expect that the approximation is good. We begin by defining the angle average of the distribution :

 f0(p,t)=12π∫2π0dθf(θ,p,t). (5)

The angle averaged distribution is then used to define by

 f(θ,p,t)=f0(p,t)+f1(θ,p,t). (6)

Clearly the angle average of is zero. Substituting Eq. (6) in the Vlasov equation (2) we have

 ∂f0(p,t)∂t+∂f1(θ,p,t)∂t+p∂f1(θ,p,t)∂θ−∂Φ(θ,t;f1)∂θ(∂f0(p,t)∂p+∂f1(θ,p,t)∂p)=0. (7)

We can derive an equation for the time variation of by taking the angle average of the last equation. We obtain

 ∂f0(p,t)∂t−12π∂∂p∫2π0dθ∂Φ(θ,t;f1)∂θf1(θ,p,t)=0. (8)

Using this equation to substitute the time derivative of in Eq. (7), we have

 ∂f1(θ,p,t)∂t + p∂f1(θ,p,t)∂θ−∂Φ(θ,t;f1)∂θ(∂f0(p,t)∂p+∂f1(θ,p,t)∂p) (9) + 12π∂∂p∫2π0dθ∂Φ(θ,t;f1)∂θf1(θ,p,t)=0.

Let us now suppose that the following disequality is satisfied during the time evolution

 f1(θ,p,t)≪f0(p,t). (10)

Then Eq. (9) can be approximated by keeping only the first order terms in , i.e.

 ∂f1(θ,p,t)∂t+p∂f1(θ,p,t)∂θ−∂Φ(θ,t;f1)∂θ∂f0(p,t)∂p=0, (11)

which is a linear equation in similar to Eq. (4), but in which is not constant in time. The procedure now is to write the solution of this linear equation in , which depends on , and to insert it in Eq. (8), to obtain a closed equation for .444We note the formal similarity between Eqs. (8) and (11) and those arising in the kinetic theory of systems with long-range interactions leading to the Lenard-Balescu equation (see Eqs. (12) and (13) in epjp ()). Both are based on a “quasilinear approximation”. However, there are crucial physical differences. We are considering here the collisionless evolution () of a Vlasov unstable distribution function while the Lenard-Balescu kinetic theory is concerned with the collisional evolution (induced by finite effects) of a Vlasov stable distribution. Thus, our kinetic theory is based on the Vlasov equation (2) instead of the Klimontovich equation (see Eq. (3) in epjp ()). Our kinetic theory is also different from the kinetic theory of the Vlasov equation (based on a maximum entropy production principle (MEPP) csr () or on another quasilinear approximation cg ()) leading to the equilibrium Lynden-Bell distribution. The kinetic theory in csr (); cg () implicitly assumes an efficient mixing while our theory, as we shall see, assumes a weak mixing. We remark the conceptual difference between Eq. (4) on the one hand, and the couple of equations given by Eqs. (8) and (11) on the other hand. In the former case is a stationary solution of the Vlasov equation (2) and, as such, it is by definition constant in time even if from the proper frequencies of Eq. (4) we could find that it is unstable. In the latter case, even if solving Eq. (11) we find that does not have exponential growth, in general defined by Eq. (8) depends on time.

To analyze the solution of Eq. (11) we first consider Eq. (4). Then, we generalize our results to Eq. (11). The solution of Eq. (4) is best obtained by decomposing in Fourier components. To study its -th Fourier component, we substitute in Eq. (4) the plane wave expressions

 f1(θ,p,t) = ak(p)ei(kθ−ωt), (12) Φ(θ,t) = bkei(kθ−ωt). (13)

As usual, the physically meaningful solutions are, separately, the real and the imaginary parts of these expressions. The wavenumber takes all integer values from to . Since and are real, in their Fourier expansion (as in Eq. (30) below) there will be both components with and . Substituting Eq. (12) and Eq. (13) in Eq. (4) we have

 −iωak(p)+ipkak(p)−ikbk∂f0∂p=0. (14)

Obviously and are related. Using the Fourier decomposition of the potential

 V(θ)=12π∞∑k=−∞Vkeikθ, (15)

where

 Vk=∫2π0dθV(θ)e−ikθ, (16)

we get

 bk=Vk∫dpak(p). (17)

With the natural assumption we have that is real and that . Substituting Eq. (17) in Eq. (14) this equation becomes

 (pk−ω)ak(p)−kVk∂f0∂p∫dp′ak(p′)=0. (18)

We have to look for which values of this equation admits solutions which are not identically zero. Since this equation is linear, we can freely impose the normalization condition

 ∫dp′ak(p′)=1, (19)

so that the equation becomes

 (pk−ω)ak(p)−kVk∂f0∂p=0. (20)

As a matter of fact this equation is satified (for ) for any real value of , that we denote by , by putting

 ak(p;ωR)=kVkP∂f0∂ppk−ωR+|k|c(k,ωR)δ(pk−ωR), (21)

where denotes that the principal value is taken when integrating over . The normalization condition for provided by Eq. (19) implies that

 c(k,ωR)=1−kVkP∫dp∂f0∂ppk−ωR. (22)

Therefore the real values of are not solutions of a dispersion relation, and thus they are not proper frequencies of the linear equation (20) in the usual sense. These solutions are present also when the Fourier component is zero; in that case is a simple delta function and . Let us note that . For the only solution of Eq. (20) is .

Another type of solutions of Eq. (20) can be obtained for particular complex values of . They are given by

 ak(p;ω)=kVk∂f0∂ppk−ω, (23)

provided the normalization condition (19) is satified. Integrating both sides with respect to we see that the complex value of must satisfy

 D(k,ω)≡1−kVk∫dp∂f0∂ppk−ω=0. (24)

The function is called the response dielectric function. There are isolated complex values of such that (dispersion relation). These are the proper frequencies of Eq. (20). There are two important relations to note. Since is real, then if we have that also and . Therefore, if for a given wavenumber there is a complex proper frequency , there is also the complex conjugate proper frequency , while for the wavenumber there will be the proper frequencies and . Correspondingly, we have and .

From the Plemelj formula

 limη→0+1x−x0∓iη=P1x−x0±iπδ(x−x0) (25)

we can obtain the value of for real values of . If and are the real and imaginary parts of , we have

 limωI→0±D(k,ω)=1−kVkP∫dp∂f0∂ppk−ωR±iπk∂f0∂p∣∣∣p=ωR/k=c(k,ωR)±iπk∂f0∂p∣∣∣p=ωR/k. (26)

The different limit from above and below the real line shows that Eq. (24) defines two different analytic functions of , one in the upper plane and one in the lower plane. The two limits coincide only for the particular values of for which the last term in the right-hand side of Eq. (26) is zero, i.e., for equal to the values of where is minimum or maximum. If it happens that for these particular values we also have , then these real frequencies belong to the proper frequencies, since they satisfy . Therefore among the isolated solutions there can also be real such that and . In this case the normalized is given by Eq. (23) with and without the necessity of the principal value.

It is possible to show case () that, for any given , the functions in Eq. (21) and in Eq. (23) constitute a complete orthogonal system (with the usual scalar product), when in the former runs over the real axis and in the latter runs over the discrete solutions of the dispersion relation .555Equation (4) is not self-adjoint, and for each solution or there is a corresponding solution or of the adjoint equation. The orthogonality of the system is expressed by the fact that each given solution of Eq. (4) is orthogonal to all the solutions of the adjoint equation, except the one corresponding to it case (). Therefore, any function can be expanded as

 h(p)=∫dωRα(ωR;k)ak(p;ωR)+∑jβj(k)ak(p;ωj(k)) (27)

with proper coefficients and . As just mentioned, the sum runs on the discrete values (generally both complex and real) satifying . Therefore, if the initial value of of Eq. (4) is , the following evolution is given by

 f1(θ,p,t)=eikθ[∫dωRα(ωR;k)ak(p;ωR)e−iωRt+∑jβj(k)ak(p;ωj(k))e−iωj(k)t]. (28)

The corresponding evolution of is

 ∂Φ(θ,t)∂θ=ikVkeikθ[∫dωRα(ωR;k)e−iωRt+∑jβj(k)e−iωj(k)t]. (29)

A general initial value of can be expanded in Fourier series as

 f1(θ,p,0)=12π∞∑k=−∞fk(p)eikθ, (30)

where the term with can be assumed to be absent (it simply corresponds to a constant in time), and where . Expanding any in the form given by Eq. (27), the evolution of is given by

 f1(θ,p,t)=12π∞∑k=−∞eikθ[∫dωRα(ωR;k)ak(p;ωR)e−iωRt+∑jβj(k)ak(p;ωj(k))e−iωj(k)t]. (31)

The reality of this expression follows from and when the proper frequencies for and are numbered so that for the same we have for and for . Finally, the evolution of is

 ∂Φ(θ,t)∂θ=12π∞∑k=−∞ikVkeikθ[∫dωRα(ωR;k)e−iωRt+∑jβj(k)e−iωj(k)t], (32)

which is also real.

Let us now compute the integral in of the product of and , the expression that appears in Eq. (8). Multiplying Eqs. (31) and (32) and integrating, we obtain

 ∫2π0dθ∂Φ(θ,t)∂θf1(θ,p,t) (33) = −i12π∞∑k=−∞kVk[∫dω′Rα(ω′R;−k)e−iω′Rt+∑nβn(−k)e−iωn(−k)t] × [∫dωRα(ωR;k)ak(p;ωR)e−iωRt+∑jβj(k)ak(p;ωj(k))e−iωj(k)t].

When is stable, there are no complex solutions of Eq. (24). This follows from the fact that its complex solutions come in complex conjugate pairs, and therefore if there are complex solutions there are necessarily unstable modes.666If is a solution of the dispersion relation (24), then is also a solution. If the first one is stable the second one is unstable, and vice versa. Thus, for a stable , the term with the sum in Eqs. (31) and (32) is either absent or with only real proper frequencies. Then, , and the integral of their product, Eq. (33), will not have terms with exponential growth. Actually, the interferences between the real frequencies will give rise to an exponential decay (the Landau damping).777We remark that, here, we are treating the linearized Vlasov equation à la Van Kampen vankampen (), i.e., with the Fourier transform in time. This treatment is the most natural one to find the proper frequencies of a linear problem. However, it is also possible to treat the problem à la Landau landau (), i.e., with the Laplace transform in time. The latter procedure is more suitable to study an initial value problem and to obtain the collective damped modes that result from the interference effects (in particular the Landau damping). In the treatment with the Laplace transform, the dielectric function is defined by Eq. (24) only in the upper plane; in the rest of the complex plane, it is defined by the analytic continuation. In the treatment à la Landau one sees how the Landau damping stems from the real frequencies of the treatment à la Van Kampen (see, e.g., nicholson () for mathematical details.) On the other hand, when there are complex proper frequencies, and then with positive imaginary parts (since they come in complex conjugate pairs), there will be terms with exponential growths. In that case, one could approximate the right-hand side of Eq. (33) by keeping only those terms. With the further assumption that for each there is only one complex proper frequency888According to the Nyquist theorem nyquisthmf (); nyquistaa (), for an unstable single humped distribution function , there is only one unstable mode. In more general cases, if there are several complex proper frequencies, one selects the largest one, i.e., the one that corresponds to the maximum growth rate. with (and its complex conjugate) we obtain

 ∫2π0dθ∂Φ(θ,t)∂θf1(θ,p,t)≈−i12π∞∑k=−∞kVk|β(k)|2ak(p;ω(k))e−i(ω(k)−ω∗(k))t. (34)

Using the expression of from Eq. (23) the previous expression can be transformed in

 ∫2π0dθ∂Φ(θ,t)∂θf1(θ,p,t)≈12π∑k>0k2Vk|β(k)|22ωI(k)(pk−ωR(k))2+ω2I(k)∂f0(p)∂pe2ωI(k)t. (35)

We now consider Eq. (11). For a linear equation with time dependent coefficients there are not proper frequencies in the usual sense, and the solutions are not of the same form as in Eqs. (12) and (13). However, we can try, for a given Fourier component, similar expressions, i.e.

 f1(θ,p,t) = ak(p,t)eikθe−i∫t0dt′ω(t′), (36)
 Φ(θ,t) = bk(t)eikθe−i∫t0dt′ω(t′). (37)

We substitute Eqs. (36) and (37) in Eq. (11) and we neglect the time derivatives of and , so that we obtain an equation similar to (14):

 −iωak(p,t)+ipkak(p,t)−ikbk(t)∂f0(p,t)∂p=0. (38)

With this adiabatic approximation, in which the time derivatives of and are neglected, the last equation is solved as before, with the only difference that and thus depend on time, and the time dependence of is only through . Without repeating the whole procedure, we write directly the expression for the integral of the product of and in the approximation that for each there is only one proper frequency with positive imaginary part (see footnote 8):

 ∫2π0dθ∂Φ(θ,t)∂θf1(θ,p,t)≈12π∑k>0k2Vk|β(k)|22ωI(k,t)(pk−ωR(k,t))2+ω2I(k,t)∂f0(p,t)∂pe2∫t0dt′ωI(k,t′). (39)

This equation can be written in a different form by defining a diffusion coefficient

 D(p,t)=∑k>02χk(t)ωI(k,t)(pk−ωR(k,t))2+ω2I(k,t) (40)

with

 ∂χk(t)∂t=2ωI(k,t)χk(t), (41)

where some constants have been incorporated in the definition of . Therefore, Eq. (39) becomes

 ∫2π0dθ∂Φ(θ,t)∂θf1(θ,p,t)≈D(p,t)∂f0(p,t)∂p. (42)

Substituting the right-hand side of this expression in Eq. (8) we obtain

 ∂f0(p,t)∂t=∂∂p(D(p,t)∂f0(p,t)∂p), (43)

which is of the form of a diffusion equation. This is the QL approximation for the evolution of the angle averaged distribution function .

Let us finally discuss when we expect that this approximation is good, i.e., when Eq. (10) is satisfied throughout the time evolution. If we assume that is of order , we can argue that the approximation is good as long as for each . This defines the smallness of also in relation to the initial values . We note that, in general, we expect that decreases in time since the dynamics drives the function towards a Vlasov (marginally) stable distribution.

## Iii Application to the HMF model

The expressions simplify for the HMF model where the potential is . Then, we have

 Vk=2πδk,0−π(δk,1+δk,−1). (44)

As a consequence, the dielectric function has a structure only for , and there is a collective dynamics only for these wavenumbers:

 D(k,ω)=1+πk(δk,1+δk,−1)∫dp∂f0∂ppk−ω. (45)

The diffusion coefficient of the QL approximation is then

 D(p,t)=2χ(t)ωI(t)(p−ωR(t))2+ω2I(t), (46)

where the evolution of is given by Eq. (41). It is sufficient to study the dispersion relation only for . It is given by

 1+π∫dp∂f0(p,t)∂pp−ω=0. (47)

Things simplify further if we consider a function that is initially even in and with only a single maximum at . These properties are conserved by the diffusion equation for (we will prove this below). It is then easy to show that if there is a proper frequency with a positive imaginary part (unstable case), its real part must be zero. Indeed, we can write the dispersion relation as

 1+π∫dpf′0(p,t)(p−ωR(t))(p−ωR(t))2+ω2I(t)+iπωI∫dpf′0(p,t)(p−ωR(t))2+ω2I(t)=0. (48)

The real and imaginary parts must be separately equal to zero. However, the last integral is negative (positive) definite if , and it is zero when . This can be seen as follows: (i) for it vanishes since is odd; (ii) for , exploiting that is odd and that it is negative (positive) for (), the sum of the values of the integrand for and is negative; (iii) analogously for this sum is positive. To be more explicit, consider for example the case . Then, for any , we have , i.e., . Thus, since for one has , the sum of the contributions of and to the last integral in Eq. (48) is negative. Therefore, since we must have , the dispersion relation becomes999Another proof of this result is given in Sec. 2.8. of nyquisthmf (). Considering an unstable single humped symmetric distribution , and using the Nyquist theorem, one can show that the dispersion relation (48) has a unique solution with ; it is such that .

 1+π∫dppf′0(p,t)p2+ω2I(t)=0, (49)

and the diffusion coefficient

 D(p,t)=2χ(t)ωI(t)p2+ω2I(t). (50)

Summarizing, for the HMF model, the integration of the diffusion equation (43) is performed with given by Eq. (50), where is the unique positive solution of Eq. (49), with single humped and even, and where the evolution equation for is

 dχ(t)dt=2ωI(t)χ(t). (51)

We remark that Eq. (43) is a diffusion equation with a time (and velocity) dependent diffusion coefficient given, for the HMF model, by Eq. (50), which is different from as long as does not vanish. If it happens that vanishes at a certain time, then the solution of the diffusion equation reaches a stationary state. Actually, as already noted, we expect that the diffusion equation will drive the function towards a stationary state by monotonically decreasing up to . We do not have a general proof of this, but we can show that, if at time the imaginary pulsation is positive, it has to tend to for increasing time [see Eq. (62) below]. We cannot prove that decreases monotonically; however, we can give an approximate argument, described in Appendix A, to show that this is likely to be the case. The argument is based on the result presented below in Eq. (59), and on the fact, found in the numerical integration of the diffusion equation and shown in the next Section, that the function keeps approximately the same functional form while it varies as driven by the diffusion equation.

It is useful to make the following evaluation. If we suppose that, at a given time , is a Gaussian

 f0(p,t)=12π√β(t)2πe−12β(t)p2, (52)

which is normalized such that , then the dispersion relation (49) becomes

 1−√β38π∫dpp2e−12βp2p2+ω2I=0, (53)

where we have not indicated the time dependence of and . The integral in the last expression can be expressed in term of the error function

 erf(x)=2√π∫x0dye−y2. (54)

We have

 1−12β+ωI√πβ38e12βω2I⎡⎣1−erf⎛⎝√β2ωI⎞⎠⎤⎦=0. (55)

For a general distribution function , we define the kinetic temperature by

 T(t)=⟨p2⟩=2π∫+∞−∞dpp2f0(p,t). (56)

This is the variance of the distribution. For the Gaussian distribution (52), we have . Therefore, Eq. (55) relates the growth rate to the temperature for the Gaussian distribution (52). Let us suppose that at time the function is a Gaussian. Its evolution, as determined by Eq. (43), together with Eqs. (50) and (51), will not maintain this structure, i.e., we cannot hope that there will simply be a change with time of the parameter . However, in our simulations we will see that the true evolution does not depart very much from a Gaussian, at least for the cases considered.

We now prove that an initial which is even in and with only a single maximum at will keep these properties during the dynamics determined by Eq. (43). We have seen that, as long as is of this form, then , and the diffusion coefficient is given by Eq. (50). This is an even function of , with a single maximum at . Therefore, the right-hand side of Eq. (43) is even in . Furthermore, taking the derivative of Eq. (43) with respect to we have

 ∂∂t∂f0∂p=∂2D∂p2∂f0∂p+2∂D∂p∂2f0∂p2+D∂3f0∂p3. (57)

If develops a point with , at that moment, by continuity, we will have and . Then, we will have . This proves that the function will remain with a single maximum at .

We also prove that the average kinetic energy increases with time. In fact, we have

 ∂∂t∫dpp22f0(p,t) = ∫dpp22∂f0(p,t)∂t =∫dpp22∂∂p(D(p,t)∂f0(p,t)∂p) = −∫dppD(p,t)∂f0(p,t)∂p. (58)

Substituting the expression of the diffusion coefficient (50) in Eq. (III), and using the dispersion relation (49), we obtain

 ∂∂t∫dpp22f0(p,t)=−2χ(t)ωI(t)∫dpp∂f0(p,t)∂pp2+ω2I(t)=2χ(t)ωI(t)π>0. (59)

Since

 ϵkin(t)=2π∫dpp22f0(p,t)=12T(t), (60)

we find that

 dTdt=8χ(t)ωI(t)>0, (61)

and we conclude that the kinetic temperature increases with time. Therefore, the velocity distribution has the tendency to spread. The increase of the kinetic temperature also qualitatively confirms that the system tends to be less unstable as time goes on.

As anticipated, we can prove that must necessarily tend to . In fact, the relation

 ∂∂t∫dpf0(p,t)2 = 2∫dpf0(p,t)∂f0(p,t)∂t =2∫dpf0(p,t)∂∂p(D(p,t)∂f0(p,t)∂p) = −2∫dpD(p,t)(∂f0(p,t)∂p)2≤0, (62)

is a sort of -theorem for the functional similar to the enstrophy in two-dimensional turbulence.101010This inequality is actually true for any -function thlb () of the form where is convex () since . Since and since is bounded from below by