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

###### 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.

###### pacs:

## 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 failure^{1}^{1}1Although 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.^{2}^{2}2The 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

(1) |

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

(2) |

where the mean field potential is given by

(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.^{3}^{3}3It 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 :

(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 :

(5) |

The angle averaged distribution is then used to define by

(6) |

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

(7) |

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

(8) |

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

(9) | |||||

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

(10) |

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

(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 .^{4}^{4}4We 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

(12) | |||||

(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

(14) |

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

(15) |

where

(16) |

we get

(17) |

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

(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

(19) |

so that the equation becomes

(20) |

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

(21) |

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

(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

(23) |

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

(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

(25) |

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

(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
.^{5}^{5}5Equation
(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

(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

(28) |

The corresponding evolution of is

(29) |

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

(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

(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

(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

(33) | |||||

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.^{6}^{6}6If
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).^{7}^{7}7We 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
frequency^{8}^{8}8According
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

(34) |

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

(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.

(36) |

(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):

(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):

(39) |

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

(40) |

with

(41) |

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

(42) |

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

(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

(44) |

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

(45) |

The diffusion coefficient of the QL approximation is then

(46) |

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

(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

(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
becomes^{9}^{9}9Another 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 .

(49) |

and the diffusion coefficient

(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

(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

(52) |

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

(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

(54) |

We have

(55) |

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

(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

(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

(58) |

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

(59) |

Since

(60) |

we find that

(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

(62) |

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