In the first chapter, we briefly introduce some essential results about the thermodynamics and dynamics of long-range system and in particular, self-gravitating systems. We present notably a derivation of the Vlasov equation through the BBGKY hierachy and discuss its limitations.

Introduction

Interactions are canonically characterized as short-range or long-range on the basis of the fundamental distinction which arises in equilibrium statistical mechanics between interactions for which the energy is additive and those for which it is nonadditive. For a system of particles interacting via two body interactions with a pair potential V (r), the system is then long-range if and only if V (r) decays at large distances slower than one over the separation r to the power of the spatial dimension d. In the last decade there has been considerable study of this class of interactions (for reviews, see e.g. [1,2]). One of the very interesting results about systems in this class which has emerged -essentially through numerical study of different models -is that, like for the much studied case of gravity in astrophysics, their dynamics leads, from generic initial conditions, to so-called quasi-stationary states (QSS): macroscopic non-equilibrium states which evolve only on time scales which diverge with particle number. Their physical realizations arise in numerous and very diverse systems, ranging from galaxies and "dark matter halo" in astrophysics and cosmology (see e.g. [START_REF] Binney | Galactic Dynamics[END_REF]), to the red spot on Jupiter [START_REF] Bouchet | [END_REF], to laboratory systems such as cold atoms [5], and even to biological systems [6]. Theoretically these states are commonly interpreted in terms of a description of the dynamics of the system by the Vlasov equation of which they represent stationary solutions. The non linear Vlasov equation describes the evolution of a smooth one particle distribution, in the limit of a large system N ! +1. In this limit, the Vlasov equation is exact. For finite systems, the Vlasov equation only describe its evolution up to a finite time diverging with the system size. These results are valid for conservative system in the micro-canonical framework.

In nature, however, real systems like those mention above, are in general finite, not isolated, and the entities (or particles) that compose the system are macroscopic, have internal degrees of freedom with short range interactions which might not be conservative. Thus the question naturally arises, is whether the presence of such QSS is valid and more generally, how the dynamics of long-range system is modified, beyond the idealised limit of Vlasov. The questions at the center of the thesis are the following: Under what conditions on a pair interaction does one expect the QSS to exist for isolated (Hamiltonian) systems? Do we expect QSS to exist also when such Hamiltonian system is subjected to additional dissipative forces? Do we expect QSS to exist in the presence of small stochastic perturbations to the system?

INTRODUCTION

In a second chapter, we explore the conditions on a pair interaction for the validity of the Vlasov equation to describe the dynamics of an interacting N particle system in the large N limit. Using a coarse-graining in phase space of the exact Klimontovich equation for the N particle system, we evaluate, neglecting correlations of density fluctuations, the scalings with N of the terms describing the corrections to the Vlasov equation for the coarse-grained one particle phase space density. Considering a generic interaction with radial pair force F (r), with F (r) ⇠ 1/r at large scales, and regulated to a bounded behaviour below a "softening" scale ", we find that there is an essential qualitative difference between the cases < d and > d, i.e., depending on the the integrability at large distances of the pair force. In the former case, the corrections to the Vlasov dynamics for a given coarse-grained scale are essentially insensitive to the softening parameter ", while for > d the amplitude of these terms is directly regulated by ", and thus by the small scale properties of the interaction. This corresponds to a simple physical criterion for a basic distinction between long-range (  d) and short range ( > d) interactions, different to the canonical one (  d + 1 or > d + 1 ) based on thermodynamic analysis. This alternative classification, based on purely dynamical considerations, is relevant notably to understanding the conditions for the existence of so-called quasi-stationary states in long-range interacting systems. This chapter follows very closely the content of an article submitted to Physical Review E [7].

In a third chapter, which is an extended report on the study of which the highlight results were published in Physical Review Letters [8], we consider long-range interacting systems subjected to additional dissipative forces. Using an appropriate mean-field kinetic description, we show that models with dissipation due to a viscous damping or due to inelastic collisions admit "scaling quasi-stationary states", i.e., states which are quasi-stationary in rescaled variables. A numerical study of one dimensional self-gravitating systems confirms both the relevance of these solutions, and gives indications of their regime of validity in line with theoretical predictions. We underline that the velocity distributions never show any tendency to evolve towards a Maxwell-Boltzmann form.

In the last chapter, which is a report on work in progress (and not yet published), we study two different toy models to explore how different kinds of stochastic perturbation affect the dynamics of long-range interaction system, and, in particular, the QSSs which are characteristic of them. In the models used, both extensions of the model with inelastic collision of the previous chapter, the perturbation is associated with the collision events, which can lead to gain or a loss of kinetic energy. For a first model, the collision are built such that on average over the realizations, the total energy in conserved, and in a second model, they are such that the system relax to a state where the loss and the gain of energy balance. For both models the perturbation drives the system to evolve, through a family of QSS, into a nonequilibrium stationary state. This final stationary distribution, which is itself also a stable QSS, does not depend on the initial conditions but does depend strongly on the microscopic detail of the perturbation.

Thus our conclusion is that QSS appear to be a very robust feature of longrange systems, and their existence or occurrence is not conditioned in general on the idealization that the system is isolated.

Chapter 1

Introduction to long-range interacting systems 1 Thermodynamics and dynamics

Definition

A long-range interacting system is usually defined as collection of entities (particles, spins, vortices,..) interacting with a pair potential which is not integrable at r ! +1. As we explain below, these systems have the property of non additivity and this property has strong consequences on the equilibrium statistical mechanics of the system.

To be a little more precise, and explain the origin of this definition let us estimate, in d dimensions, the potential energy e p of a particle placed at the center of sphere of radius R enclosing a volume V = ⌦ d R d where other particles are homogeneously distributed (where ⌦ d is the appropriate constant for the dimension d). We consider that particles interact with a power low central potential defined as

(r) = g |r| ↵ , (1.1)
where g is the coupling constant of the model and ↵ the exponent characterizing the interaction. We use a cut-off ", to regulate the possible short-scale divergence of (for ↵ > 0) i.e. we do not consider the contributions to e p from particles located within a sphere of small radius " around the particle considered. We will return below to this point to discuss further the role and necessity for this cut-off. In the continuous limit, the individual potential energy is expressed as:

e p = Z R " d d rn 0 g r ↵ = 0 J⌦ d Z R " drr 1 = n 0 g⌦ d d ↵ (R d ↵ " d ↵ ) , if ↵ 6 = d
where n 0 = N V is the particle density. The following distinction then follows.

• For ↵ > d, the potential is integrable: the energy per particle remains finite if we take R ! +1 at fixed n 0 (i.e. the usual thermodynamic limit). In this limit, the contribution to e p due to far away particles is negligible compared to that due to the ones close to it. Further, the total energy, E = R V d d rn 0 e p , 3 is extensive: it scales linearly in the thermodynamic limit with the volume V .

For such system, the energy per particle is an intensive quantity and it is well defined in thermodynamic limit. This class of system, which we will refer to, following the canonical usage, as short-range interacting system, includes most of the familiar systems studied in physics, and in particular those usually considered in the context of the study of phase transitions and critical phenomena. In this scope one usually defines the exponent of the potential with the parameter = ↵ d, and the studies are usually restricted to the case > 0. To avoid ambiguity, we mention that in this latter community the term "long-range" is then often used to qualify a universality class of systems for which the critical exponents characterizing the divergence of physical quantities, affected by the nature of the interaction, are different to the one of the mean field limit and induce longrange correlations [2,9]. This typically happens in the range d/2 < < (d) where (d) is a critical exponent for which the precise d-dependence (in d  4) is still discussed [10,11].

• For ↵ 6 d, the potential energy per particle (or density of potential energy) e p diverges for R ! +1 in the usual thermodynamic limit. The potential is non-integrable. This is by definition the case of long-range interacting systems 1 . One cannot neglect the contribution of e p due to "far-away" particles. The energy per particle scaling as e p / R d ↵ / V 1 ↵/d , is no more intensive and the total potential energy of the system U / V 2 ↵/d is super-extensive: it increases more rapidly than linearly with the volume.

The long-range interactions concerns a large variety of systems [START_REF] Dauxois | Dynamics and Thermodynamics of Systems with Long-Range Interactions[END_REF] such as plasma [13], self-gravitating systems [START_REF] Binney | Galactic Dynamics[END_REF]14], 2d turbulent flow [15], cold atoms [5,16], biological systems [6]. Further, many theoretical toy models such as the Hamitonian mean field (HMF) model(e.g. [1] and reference therein), the Blume-Emery-Griffiths (BEG) model [17], the 1d self-gravitating [18][19][20] or self-gravitating ring model [21] have been proposed to study the physics of long-range systems.

The distinctive thermodynamics of long-range systems

Various crucial differences with respect to short-range systems emerge immediately when one applies equilibrium statistical mechanics to long-range systems.

The first and fundamental one is that related to the definition given above. When we study the equilibrium thermodynamics of any system we must have extensivity of the total energy E as N ! +1 and V ! +1 (in order to have an N independent energy density). In order to recover this property, we must define the thermodynamic limit differently. There are various mathematically equivalent ways of doing this. One, the so-called Kac prescription, is to scale the coupling constant g with the size of the system as V (1 ↵/d) with the density n 0 = N V is fixed, and the coupling constant decreases in increasing the system size.

Alternatively, in the context of 3d self-gravitating system (for which g = Gm 2 ) one can use the "dilute" limit with g / 1 N 2 and n 0 / V (1+ d ) . Recently, another scaling limit have been proposed in [22], the author propose a scaling limit for a self-gravitating system of hard spheres keeping constant the packing fraction. Such a procedure restores the extensivity of the system; but nevertheless it remains non-additive i.e. the interaction energy of any part of the system with the whole is not negligible with respect to the internal energy of the given part. We detail further in the next section this point and discuss its particular consequence on the thermodynamics of long-range interacting systems.

A second issue which present itself immediately in analysing the equilibrium thermodynamics of long-range systems is one which occurs for any attractive power law potential 1/r ↵ with ↵ > 0. An interesting example, of this class of system includes, 3d self-gravitating system whose thermodynamics is discussed further below. Just as for short-range systems whose exponents are also positive, such systems require a regularisation at short distances to avoid collapse. Indeed if only two particles get closer and closer one another, these two can make the potential energy of the entire system diverges to 1. The total energy being conserved, the kinetic energy of the system increases and the number of accessible states diverges in the velocity space. As a consequence, the micro-canonical partition function:

⌦(E) / Z d dN xd dN x (E H({x N , v N })), (1.2)
which enumerates the number of microscopic states at a fixed energy E, diverges if and only if N 3 [23]. Moreover, the Boltzmann entropy

S(E) = k B ln(⌦(E)) (1.3)
increases without bound as well; and without entropy maxima an equilibrium state cannot exist. Therefore the introduction of a short-scale regulator is needed to bound below the potential energy. For spin systems defined on a lattice like the BEG models, this divergence is naturally regulated. For particle system, one may also consider fermions and made use of the Pauli exclusion principle as it a been done for a selfgravitating system [24], or again hard core particle as in [22]. In general, an ad-hoc cut-off is added to bound the potential. We will discuss in chapter 2, how the large scale dynamics of a long-range system is affected by the choice of the small scale regularisation.

Finally, just as for systems, the number of accessible states (Eq. (1.2)) will also diverge if the system is not confined, unless the pair potential is attractive enough to "self-confine" the particles (as is the case, as we will see, for example for 1d selfgravitational system). Thus, one has often to impose periodic boundary condition, or confine the system in a box.

Non-additivity

The central difficulty, as we have underlined, characteristic of any long-range system, arises from the non-integrable nature of the long-range potential at large distance. If the Kac prescription restore the extensivity of the long-range interacting system Figure 1.1: Schematic picture of a system partitioned into two parts, in red is the boundary between the two subsystems. and makes possible the existence of a thermodynamic limit, such system, however, remains non-additive. Let us consider first a short-range system, and partition the system in two subsystems 1 and 2 (as in Fig. 1.1) of respective energies E

1 and E 2 . The total energy of the system is given by

E = E 1 + E 2 + E
12 where E 12 is the energy of interaction between the two subsystems. For a short-range system, a particle, situated close to the boundary between the two subsystems, interacts only with particles close to the other side of the boundary, while the interactions with the particle in the bulk of the other subsystem are negligible. Therefore, the energy of interaction E

12 is proportional to the surface of contact between the two volumes. Then, in thermodynamic limit E 12 becomes negligible with respect to E 1 and E 2 and in this limit E = E 1 + E 2 : the system is additive. On the contrary for a long-range system, any particle in one of the two subsystem interacts with every single particles of the other subsystem. The energy of interaction E 12 is no longer negligible at the thermodynamic limit and

E 6 = E 1 + E 2 .

Ensemble inequivalence

An important result of the statistical mechanics built for short-range systems is the equivalence of statistical ensemble which means that mean quantities converge to the same value in the thermodynamic limit. Indeed this result allows one, for example, to study the phase diagram of a system in the ensemble of its choice. Short-range interacting systems are usually study in the canonical ensemble, because calculations are easier than in the micro-canonical ensemble. The latter is essential, however, at a fundamental level to construct the canonical ensemble [1,[START_REF] Dauxois | Dynamics and Thermodynamics of Systems with Long-Range Interactions[END_REF][START_REF] Brunet | Cours de Physique Statistique[END_REF]. We now briefly recall the latter to highlight how the distinction on the basis of additivity is important between short-range and long-range interacting systems.

One considers an isolated system with energy E and divides it into two parts: a "small" part of energy E 1 that will be the system we want to study the statistics and a "large" part of energy E 2 E

1 playing the role of the thermal bath. The probability that the system 1 has an energy E 1 , for an additive system is given by

p(E 1 ) = 1 ⌦(E) Z (E 1 + E 2 E)⌦ 2 (E 2 )dE 2 (1.4) = ⌦ 2 (E E 1 ) ⌦(E) (1.5)
where ⌦ 2 is the micro-canonical partition function of the subsystem defined similarly than ⌦ in Eq. (1.2). Then, using the expression of the micro-canonical entropy (Eq. (1.3)), and Taylor expanding with E >> E 1 , we have:

p(E 1 ) = 1 ⌦(E) e 1 k B S(E E 1 ) (1.6) = 1 ⌦(E) e 1 k B S(E)+ 1 k B E 1 @S @E | E +o(E 2 1 )
(1.7)

' 1 Z e E 1 , (1.8)
with Z the partition function and where we have define the inverse temperature of the bath:

= 1 k B @S @E E .
(1.9)

We thus recover the canonical distribution for the system 1. We clearly see (Eq. (1.4)) that the property of additivity is crucial to build the canonical ensemble. For a short-range system, of course, this construction become strictly valid at the thermodynamic limit. For a long-range system, this construction is not mathematically valid even at the thermodynamic limit. Hence non additive systems, and in particular long-range one, are expect to have a very particular behaviour when they are in contact with a thermal reservoir. The micro-canonical ensemble is clearly well defined and the only issue is how to take the limit N ! +1. There is no formal barrier to calculate, for long-range system, in the canonical ensemble, but its physical meaning must be scrutinized carefully. For many long-range model, the calculation have been done in both ensembles, defining formally the free energy 2and the canonical partition function.

For some long-range systems, it turns out that the equilibrium results in both ensembles are actually equivalent. This is, for example, the case of the 1d selfgravitating system [20,[START_REF] Worrakitpoonpon | Relaxation of isolated self-gravitating systems in one and three dimensions[END_REF][START_REF] Sire | [END_REF], or of the HMF model [1]. For these two models, the thermal equilibrium can be work out in both ensembles and the equilibrium statistical mechanics is well defined. Moreover the dynamical studies have shown that the two long-range models effectively relax to the predicted equilibria which has, in particular, a Gaussian velocity distribution [1,28].

However for other long-range systems, like in 3d self-gravity [START_REF] Sire | [END_REF]29] or for BEG model [17], the calculation shows that the ensembles are not equivalent: the equilibrium solution found in maximising the entropy in the micro-canonical ensemble, or in minimizing the free energy in the canonical ensemble, are different in a substantial region of the space of the parameters.

Non concave entropy and negative specific heat

It has been understood and shown rigorously [17,30] that the fundamental explanation of the non-equivalence of ensemble is the presence of a non-concave region in the entropy-energy curve in the micro-canonical ensemble. To be more precise, the thermodynamic potential in the micro-canonical ensemble (S(E)) and the one in the canonical ensemble (the free energy F ( )) are related by a Legendre-Fenchel transform [START_REF] Touchette | Legendre-Fenchel transforms in a nutshell[END_REF]. We denote this operation with the subscript ⇤. The entropy S and the free-energy F are thus related by: S ⇤ ( )

. = F ( ) = inf E { E S(E)} (1.10)
Similarly one defines the Legendre transform of the free-energy which is the canonical entropy:

S can . = S ⇤⇤ (E) . = ( F ) ⇤ (E) = inf { E F (E)} (1.11)
In points where the curve are differentiable, this transformation implies the relation = @S @E (1.12) which is the micro-canonical definition of the temperature. In Fig. 1.2,3 one can see that when the entropy is a concave function of the energy E, the Legendre-Fenchel transform is self-inverse. This means that if one considers the free-energy F ( ) which is the thermodynamic potential in the canonical ensemble and that one wants to go back in the micro-canonical ensemble, using a Legendre transform, then the canonical entropy S can = F ⇤ equals the micro-canonical one S. In this case then, the ensembles are equivalent.

In Fig. 1.3, we see that this is not the case for a non-concave entropy. This gives rise to a non-differentiable point is the transformed function S ⇤ = F (middle panel). Moreover when one wants to transform the free energy into an entropy, the function obtained, is the concave envelope of the original entropy and in this case the Legendre-Fenchel transform is not self-inverse: S ⇤⇤ 6 = S (right panel).

Hence, for systems with a convex intruder in the entropy-energy curve, the ensembles are not equivalent. For such systems, the micro-canonical ensemble may show richer behaviour than the canonical one. For example, considering the specific heat: related to the entropy by the relation

C v = @E @T , (1.13)
@ 2 S @E 2 = 1 T 2 1 C v . (1.14)
In the region where the curve S(E) is convex, the second derivative is positive, and then C v < 0. A counter intuitive consequence is that when ones increase the energy of the system, it may cool down. Moreover non-concave entropy is related to a rich variety of phenomena such as the appearance of first-order phase transitions as well as metastable states in the canonical ensemble [START_REF] Bouchet | [END_REF] and also the possibility of ergodicity breaking [33]. In [START_REF] Bouchet | [END_REF] the authors propose a complete classification of micro-canonical phase transitions in long-range interacting, their link to canonical ones, and of the possible situations of ensemble non-equivalence.

The results from equilibrium statistical mechanics aim to predict the final macroscopic state of the system. However, to determine if the system does effectively relax to such states, the dynamics of the system has to be studied.

Dynamics of long-range systems

Qualitative description

In a long-range interacting, the contribution to the force acting on a given particle cannot be approximated by that due to particles in its neighbourhood. The resultant strong coupling of numerous degrees of freedom allows important collective effects in the many body dynamics. The detailed motion of every single particle is, in practice, inaccessible except by numerical simulation, and what we attempt to describe analytically are the behaviours of macroscopic quantities such as the one particle phase space density.

Studies of the dynamics of long-range interacting system (e.g. [34,35] or for review [1]) -both numerical and theoretical -have lead to the conclusion that it is characterized broadly by two distinct regimes: starting from a generic initial condition, a long-range system of many bodies evolves first, on a time scale (⌧ mf or ⌧ dyn ) characterizing the evolution under the mean field, towards a macroscopically stationary state which does not coincide with the equilibria derived from the equilibrium thermodynamics. This early time process is often referred to as "violent relaxation" following the nomenclature originally given by Lynden-Bell for self-gravitating systems [36]. These so-called quasi-stationary states are believed to describe notably globular clusters, galaxies [START_REF] Binney | Galactic Dynamics[END_REF], the red spots of Jupiter [START_REF] Bouchet | [END_REF] and are reproduced in simulations (e.g.: for the HMF model [34], 1d self-gravitating system [18], 2d euler equations [37]) for various long-range systems. They are interpreted theoretically as stationary solutions of the Vlasov equation which we derive in the next section. This equation is a non linear equation for the single particle distribution function, and takes the mathematical form of a Boltzmann equation without collisional term. In the astrophysics community (see e.g. [START_REF] Binney | Galactic Dynamics[END_REF]) it is usually called the "collisionless" Boltzmann equation. It has been rigorously shown [38], for a certain class of potential (with an exponent ↵  d 1 as we will see in chapter 2) that it gives a very good approximation, at a macroscopic level, of the N body dynamics at least on time scale ⌧ V ⇠ ln(N )⌧ mf On a longer time scale diverging also with the particle number but (in general much longer than ⌧ V ), numerous long-range interacting systems have been observed to relax to thermal equilibrium [1,28]. This second phase evolution is no longer described within the framework of the Vlasov equation and additional terms accounting for "collisions" (or correlations, as we will discuss below) have to be included to describe it.

Derivation of the Vlasov equation

In this section we present a derivation of the Vlasov equation (see e.g. [39,[START_REF] Balescu | Equilibrium and Non-Equilibrium Statistical MechanicsWiley-Interscience publication[END_REF]) through the BBGKY hierarchy and a short discussion of what this equation represents.

-The BBGKY hierachy -We consider a system of N indistinguishable point particles of identical masses (or charges) m; we denote x i and v i respectively the position and the velocity of the i th particle in a d dimensional space. At time t, the information about an ensemble of realisations of the system is given by the N particle probability distribution function:

f N (x N , v N , t) . = f N (x 1 , ...x N , v 1 , ..., v N , t). (1.15)
It gives the probability of finding N particles in the infinitesimal volume element

S N i=1 [x i , x i + dx i ] [ [v i , v i + dv i ] in the 2dN dimensional space of configurations.
For s = 1, ..., N, we define the reduced s-particle phase space probability:

f s (x s , v s , t) . = f s (x 1 , ...x s , v 1 , ..., v s , t) = Z P N (x N , v N , t)dx (N s) dv (N s) , (1.16)
with the notation dz (N s) . = dz s+1 ...dz N . The conservation of the probability give the sum rule

Z f s (x s , v s , t)dx s dv s = 1, 8s = 1, ..., N, (1.17)
where dx s dv s . = Q s i=1 dx i dv i . As there is no sink or source in the space of configurations, a consequence of the divergence theorem (in 2dN -dimension) is that f N obeys the continuity equation:

@f N (x N , v N , t) @t + N X i=1  @x i @t • @f N (x N , v N , t) @x i + @v i @t • @f N (x N , v N , t) @v i = 0. (1.18)
where • denotes the scalar product and we use the notation @ @z ⌘ r z for the gradient. Now if the N particles interact via a pair potential (|x y|), the Hamiltonian of this system is given by:

H = 1 2 m N X i=1 v 2 i + N X i=0 X j>i (|x i x j |). (1.19)
and we implicitly assume a regularization at r ! 0, if needed (we will come back to this crucial point). For a Hamiltonian system, the equations of motion are:

@x i @t = @H @mv i = v i (1.20) @v i @t = @H @x i = N X i6 =j @ (|x i x j |) @x i .
Inserting then in Eq. (1.18), we obtain the Liouville equation:

@f N (x N , v N , t) @t + N X i=1  v i • @f N (x N , v N , t) @x i (1.21) X j6 =i @ (|x i x j |) @x i • @f N (x N , v N , t) @v i # = 0
Integrating Eq. (1.21) over velocities and position of N s particle we obtain the so called Bogoliubov-Born-Green-Kiriwood-Yvon (BBGKY) hierarchy and assuming that f N ! 0 when |v i |, |x i | ! +1, we have 8s = 1...N :

@f s @t + s X i=1 v i • @f s @x i s X i6 =j @ (|x i x j |) @x i • @f s @v i (1.22) = (N s) s X i=1 Z dx s+1 dv s+1 @ (x i x s+1 ) @x i • @f s+1 @v i ,
where the factor (N s) comes from integrals over (assumed) indistinguishable particles. One sees that each equation for f s is coupled to the next one for f s+1 . This set of equations contain exactly the same information as the Liouville equation Eq. (1.21).

-The Vlasov equation -It takes the from of a closed equation for the equation for the single particle phase space density f 1 (x 1 , v 1 ). In order to obtain it, one needs to close the BBGKY hierarchy at the level of f 1 with assumptions that we discuss here.

CHAPTER 1. INTRODUCTION TO LONG-RANGE INTERACTING SYSTEMS

The first equation of the BBGKY hierachy (Eq. (1.22)) is:

@f 1 (x 1 , v 1 , t) @t + v 1 • @f 1 (x 1 , v 1 , t) @x 1 (1.23) (N 1) Z @ (|x 1 x 2 |) @x 1 • @f 2 (x 1 , x 2 , v 1 , v 2 , t) @v 1 dx 2 dv 2 = 0.
Now let us decompose the two-body distribution function:

f 2 (x 1 , x 2 , v 1 , v 2 , t) = f 1 (x 1 , v 1 , t)f 1 (x 2 , v 2 , t) + g(x 1 , x 2 , v 1 , v 2 , t), (1.24)
where g is the connected or irreducible, two-point correlation function. To close the hierarchy i.e. to obtain an equation for the single particle distribution function f 1 only, we simply assume that g(x

1 , x 2 , v 1 , v 2 , t) = 0.
In other words, we assume that the pair correlation and higher order correlations are negligible, or suppose what Boltzmann called molecular chaos in the context of a dilute gas. If we suppose the initial condition, indeed, satisfies,

f N (x 1 , x 2 , v 1 , v 2 , 0) = Q N i=1 f 1 (x i , v i , 0
), one may ask until which time t this assumption is a realistic description of the system. This property for a N body system to remains uncorrelated during its dynamical evolution is called propagation of chaos. With this assumption, and in the limit of large N , the last term can be rewritten:

(N 1) @ @v 1 Z @ @x 1 f 2 (x 1 , x 2 , v 1 , v 2 , t)dx 2 dv 2 ' @f 1 @v 1 • @ ¯ [f 1 ](x 1 , t) @x 1 , (1.25)
where the mean field potential ¯ is:

[f 1 ](x 1 , t) = (N 1) Z (|x 1 x 2 |)f 1 (x 2 , v 2 , t)dx 2 dv 2 .
(1.26)

Let us now take the limit N ! +1 at fixed f 1 (x, v) (i.e. increasing the particle density f = Nf 1 in phase space). We can obtain consistently an N independent limit if we take (N 1) ⇠ C st , in other words, rescale the coupling of the interaction by 1 N . This correspond to the Kac prescription, discussed earlier, in which the coupling is scaled in order to obtain an extensive energy 4 . Indeed E ⇠ N 2 ⇠ N when ⇠ 1 N . We obtain then the Vlasov equation in its dimensionless form:

@f 1 (x, v, t) @t + v • @f 1 (x, v, t) @x @ ¯ [f 1 ](x, t) @x • @f 1 (x, v, t) @v = 0, (1.27)
with the mean field potential:

[f 1 ](x, t) = Z (|x x 0 |)f 1 (x 0 , v 0 , t)dx 0 dv 0 . (1.28)
For this equation, we recall that f 1 represents the single particle probability distribution normalised to unity.

One may alternatively rewrite the equation in terms of the mass density in the phase space f (x, v, t) = Mf 1 (x, v, t). In the limit of large N with a potential scaling as ⇠ 1

N and the mass m ⇠ 1 N , such that the total mass M = mN remains fixed, we obtain the (dimensional) Vlasov equation:

@f @t + v • @f @x @ [f ] @x • @f @v = 0, (1.29)
with the mean field potential:

[f ](x) = 1 M Z (|x x 0 |)f (x 0 , v 0 , t)dx 0 dv 0 . (1.30)
In this last expression Eq. (1.29), we see that the Vlasov equation is a continuous description of the particle dynamics in the limit of an infinite number of infinitely light particles. The mean field limit can thus be described as a fluid limit. It can be shown that the Vlasov dynamics has an infinite number of conserved quantities called Casimir invariants:

C[f ] = Z H(f (x, v, t))dxdv, (1.31)
where H is any function. The Vlasov equation admits an infinite number of stationary solutions [1,[START_REF] Chavanis | [END_REF]. The continuous dynamics described by the Vlasov equation is a major subject of research in itself and also addressed in various different fields (astrophysics, plasma physics,..). Indeed, it is a complex non-linear equation and its resolution, just as the N body dynamics, is very costly numerically (and even more so because it is an equation for the 2d dimensional phase space). This standard derivation, just given, of the Vlasov equation is not completely satisfactory: indeed, it remains quite unclear in what circumstances the essential approximation made of neglecting correlation will apply to a system in the class of Hamiltonian systems. It is known that for hard-core or short-range potential (e.g. van der Waals), the system is not well described by the Vlasov equation but instead by the Boltzmann equation with a collision term accounting for short-range correlations. Looking more carefully at the derivation, we have implicitly supposed, for example in Eq. (1.23), that integrals are not divergent. Also, for a system for which a short-scale regulation of the pair force is needed, the present derivation does not allows one to understand the effects of this regularization. We have simply stated that if we make certain approximations we obtain this equation, but have not justified the approximation with any rigour.

One of the aims of the chapter 2 of the thesis is precisely to improve our understanding of the range of application of the Vlasov equation.

Violent relaxation: the establishment of QSS

As we have stated, studies of long-range systems show that they evolve generically through a phase of violent relaxation and organise in a QSS. Both this dynamical evolution and the QSSs are believed to be described by the Vlasov equation and the QSSs are stationary states of the latter5 . As noted above, there are an infinite

CHAPTER 1. INTRODUCTION TO LONG-RANGE INTERACTING SYSTEMS
number of stationary solutions of the Vlasov equation. The question then arise naturally of how the systems selects which stationary state it relaxes to and to what extent this does or does not depend on the detail of the initial conditions for the system. Indeed, it is known that the macroscopic distribution of the QSS reached depends on the initial conditions given to the system. As noted, the early time Figure 1.4: Illustation of the "violent relaxation". The plot shows a particles position in phases space at different time (N body simulations of 1d gravity [42]) relaxation can be studied within the framework of the Vlasov equation.

Lynden-Bell proposed a statistical approach to calculate the phase space distribution of QSSs [36] in the framework of the Vlasov equation that we briefly describe here. A detailed derivation can be found for example in [1,[START_REF] Worrakitpoonpon | Relaxation of isolated self-gravitating systems in one and three dimensions[END_REF]36,42]. The Lynden-Bell theory is similar to the usual Boltzmann statistics, but instead of working with the particles, Lynden-Bell studied the distribution of the phase space density level ⌘ at a microscopic level in phase space. This approach assumes, just as for Boltzmann statistics, the existence of ergodicity with respect to the phase space configuration accessible to the Vlasov dynamics i.e. consistentlty with its conserved quantities. The phase space is divided into P "macrocells" and each of them subdivided into ⌫ "microcells" of volume h d . The volume fraction occupied by the level ⌘ inside the macrocell i is

⇢(x, v) = n i ⌫ , (1.32)
where n i is the number of microcells inside a macrocell i occupied by the level ⌘. The volume fraction is related to the distribution function by ⇢(x, v) = f (x, v)/⌘. Then with simple combinatorics, one can work out the number of possible "microstates" W (n i ) and defined the entropy S lb = k B ln(W (n i )). In the limit in which the variaton of ⇢(x, v) between macrocells are infinitesimal, one can write the entropy as:

S lb = k B Z d d xd d v h d [⇢(x, v)⇢(x, v) + (1 ⇢(x, v))ln(1 ⇢(x, v))].
(1.33)

Then maximising the latter entropy with the constraint that the energy and the number of particle is conserved, one finds the following solution, for the case of a equation induces a filamentation of the phase space. Looking at a macroscopic level, with time going on, the filaments, so thin, becomes indistinguishable. The time dependence is thus progressively transferred to small and smaller scales.

single of a single level of density:

f lb (x, v) = ⌘ 1 + e (e(x,v) µ) , (1.34)
where e(x, v) = 1 2 mv 2 + m ¯ (x) is the single particle energy, and and µ are the Lagrange multipliers related to the two constraints (interpreted as the "inverse temperature" and "the chemical potential").

It has been found through numerical studies that the LB theory can provide a very accurate description of the relaxed QSS, in certain parts of the initial condition space of specific models: in the HMF [43], and in model with 2d vortices [44]. However more generally it is clear that the Lynden-Bell theory is inadequate, in particular for system starting from state very far from a QSS. This have been shown for the 1d self-gravitating system [18,42] and for the HMF [45]. In 3d gravity notably, an ubiquitous property of QSS is to have non-isotropic velocity distribution and even non isotropic spatial distribution arising from spherically symmetric initial condition and such state are not predicted by the theory of Lynden-Bell.

The main reason why the theory does not work for any initial condition is that it supposes a mixing and ergodicity in the phase space, and these two requirements are in general not verified by long-range systems. Other approaches to predict the property of QSS have been explored notably by [42]. The model constructed by the authors, have given, in various cases (HMF, 1d self-gravity and others), framework to understand well these states based on an analysis of the dynamics of the phase of violent relaxation. In summary, during violent relaxation, the mean field potential is characterised by quasi-periodic oscillations. These breathing modes have been notably shown within the framework of kinetic equation (beyond Vlasov), in the context of cold-atoms [5]. It is possible, therefore, for some particles to enter in resonance with the oscillations and gain large amounts of energy at the expense of the collective motion. This process is known as Landau damping. The Landau damping diminishes the amplitude of the oscillations and leads to the formation of a halo of highly energetic particles which surround the high density core [19]. The phenomenon of Landau damping, first described heuristically by Landau [START_REF] Landau | Physical kinetics[END_REF] in the context of plasma physics has been recently formulated rigorously in the famous work of Villani and Mouhot [47,48] in particular limits.

Evolution of QSS: beyond Vlasov equation

The state attained by a long-range system is a quasi-stationary state because, for a large but finite N particle system, a second relaxation occurs, on a much longer time scale ⌧ R eventually bringing the system toward the thermal equilibrium (if well defined). The time scale, in general also diverges with N and its precise dependence varies with the model.

For 3d self-gravitating systems, Chandrasekahar [49] estimated the time scale for thermal relaxation to be ⌧ R ⇠ N/ln(N )⌧ dyn . His calculation, based on an estimation of the number of "collisions" a test particle has in the time it crosses the system (⇠ ⌧ dyn ). This estimate of the relaxation rate due to two-body collisions (for 3d gravity) has been generalised to a broad class of pair interactions in dimensions d 2 in [50]. The latter calculation also shows that the main contribution to this rate for a long-range system is due to small angle deflections (or "soft" collisions) rather than due to "hard" collisions giving large deflections.

While the mean field dominates the trajectories of particles, many small deflections lead eventually to significant deviations from the collective dynamics. The physics of such"collisions" is not included into the Vlasov equation describing only the collective effects, and additional term can be derived to account for them.

One common approach is to add Fokker-Plank to the Vlasov equation [START_REF] Binney | Galactic Dynamics[END_REF]14]. To compute the diffusion coefficient, the microscopic detail of the mechanism of collision is needed with all the associated hypotheses underlying the calculation of Chandrasekar. This have been done notably in the context of 2d gravity in [51].

The Lennard-Balescu equation, coming from plasma physics, is another attempt to build such terms from a consistent derivation. To obtain it, one has to close the BBGKY hierarchy at the level of the two body correlation function f 2 . In order to do so, one usually expands the pair density function as

f 2 (z 1 , z 2 ) = f (z 1 )f (z 2 )+g(z 1 , z 2 )
as discussed above (we use the simpler notation z i = (x i , v i ) and f = f 1 ), and f 3 as

f 3 (z 1 , z 2 , z 3 ) = f (z 1 )f (z 2 )f (z 3 ) + f (z 1 )g(z 2 , z 3 ) (1.35) + f (z 2 )g(z 1 , z 3 ) + f (z 3 )g(z 1 , z 3 ) + h(z 1 , z 2 , z 3 ).
with h the connected three-point correlation function.

Then assuming that h ⇠ 0 for large N , one obtains two coupled equation, closed in f and g. Then solving formally the equation for g, it is possible to obtain, for an homogeneous state (which is relevant in plasma physics but not for self-gravitating system), an expression for the additional term of the Vlasov equation in f which take into account the finite N effects. The additional term take finally the form of Fokker-Planck terms with a diffusion coefficient which is a functional of f . In practice to solve this equation is difficult and it remains unclear whether it correctly describe the evolution of relaxation of the system.

Long-range systems with stochastic perturbations

The Vlasov equation applies to strictly conservative systems in a micro-canonical framework, and the question inevitably arises of the robustness of QSSs, which are stationary solution of this equation, beyond this idealized limit.

In [52] the authors consider a long-range interacting system driven by external stochastic forces acting collectively on all the particles constituting the system. Given a long-range Hamiltonian H (Eq. (1.19)), they consider the following equation of motion: qi = @H @p i and ṗi = @H @q i ↵p i + p ↵⇠(q i , t), (1.36) where ⇠(q, t) is a statistically homogeneous Gaussian stochastic force field ( or "noise") with a zero mean and without temporal correlation i.e h⇠(q, t)i = 0 and h⇠(q i , t)⇠(q j , t)i = C(|q i q j |) (t t 0 ), (1.37) where ↵ controls the strength of the force and C is the (isotropic) correlation function.

CHAPTER 1. INTRODUCTION TO LONG-RANGE INTERACTING SYSTEMS

Starting from a generalised Liouville equation for the N -particle distribution function f N (Eq. (1.15)), they derive a kinetic theory starting from a generalized BBGKY hierarchy which takes the stochastic perturbation into account. Because the stochastic force is correlated in space it induces correlations on the system and, to treat this at a non-trivial level, the kinetic theory has to take into account the connected pair correlation g between particle in the derivation the kinetic equation. To obtain Lennard-Balescu equation, their generalised BBGKY hierarchy is closed at the level of g by neglecting the connected three-point and higher order correlation functions. They obtain thus a set of 2 equations, with both finite N effects and the stochastic effect. Considering the limit N 1 and N↵ 1, in which the characteristic time ⌧ stoch = 1 ↵ ⌧ dyn of the external perturbation is much shorter that the one of the finite N effects ⌧ fN ⇠ N⌧ dyn , the author then treat the limit ↵ ⌧ 1. The authors able to show, in particular, that detail balance is respected (and the system evolved to a Gaussian velocity distribution) if and only if the noise is white i.e. uncorrelated in space. In presence of correlated noise, the author shows that the system relaxes toward a non-equilibrium stationary state (NESS). The numerical simulations of the HMF model reported in these article give not only very good agreement with the kinetic theory developed for homogeneous state but also, for some particular value of the parameter of the perturbation, the system exhibits a very interesting bistable behaviour: the system switches abruptly and intermittently between a homogeneous and an inhomogeneous state.

In [53,54] the authors consider a perturbation of the long-range dynamics with a diffusion process à la Langevin in order to attempt an operative description of a canonical ensemble for long-range interacting. There studies correspond to the one in [52] but for an uncorrelated noise. The author report also, in particular, that the system relaxes to non-equilibrium stationary state that are QSS.

In an other study [55] the HMF model is perturbed with a energy conserving three-particle collision dynamics. At a given frequency, the velocities of three particles, chosen randomly within the whole system, are exchanged by a periodic permutation. The authors showed in this case the QSSs are "destroyed" by the stochastic dynamics, and that the perturbation makes the system relax faster the thermal equilibrium.

In the spirit of these studies, in the chapter 3 and chapter 4 we consider longrange systems subjected to various kinds of perturbations, different to those described above.

2 Self-gravitating systems

3d self-gravity

The gravitational N body problem is evidently the most notable example of a longrange interacting system. The problem can be relevant at different scales in astrophysics and in cosmology, e.g. planetary system, stars a globular cluster, stars in galaxies, galaxies in galaxy cluster. Considering such a system in an infinite space (or with periodic boundary conditions) is also a pertinent model in the context of the problem of non-linear structure formations in cosmology.

Results from equilibrium statistical mechanics

The study of the thermodynamics of three dimensional self-gravitating systems is a specific example in which the small scale regulation of the divergence mentioned above is relevant. One also has to enclose the system of N particles of masses m in a spherical box of radius R. Then, in the mean field limit N ! +1 with ER/GM 2 = C st where M is the mass of the whole system and E the energy or N ! +1 with GM m/R = C st with = (k B T ) 1 the inverse temperature), it is possible to determine the phase diagram in the canonical and micro-canonical ensemble. A result it is important to underline is that the phase diagram obtained depends on the small scales cut-off introduced to regularise the potential at r = 0 [56]. More precisely it depends not only on the nature of this cut-off (one can use a softened potential [57] or consider fermionic particles [24]) but also on its value. Without regularisation, nevertheless the problem can be solved formally using the mean field limit [14] as there are, for example, local maxima of entropy which are well defined independently of the cut-off. A famous result of the calculation is that we obtain a caloric curve spiralling around a point. It presents then a non concave region and the ensembles are not equivalent in general. With a short-scale regularisation, this is also true and lead to a very rich thermodynamics [56] [24].

When the energy in the micro-canonical ensemble or the temperature in the canonical ensemble drops below a certain critical value E c or T c , respectively, the corresponding thermodynamic potentials undergo a discontinuous jump [14,56]. If no short-range cut-off is introduced, the discontinuous jump is infinite and the entropy and free energy diverge, and then no extrema can be found. This makes all normal (non-singular) states of the self-attractive system metastable with respect to such a collapse; the collapse energy E c is in fact an energy below which the metastable states cease to exist [57]. Above this energy, isothermal sphere type solutions [14] exist. If, on the other hand, a short-range cutoff is introduced, the entropy and free energy jumps are finite. In this case, as a result of the collapse, the system goes into a non-singular state with a dense core and a sparse halo. The precise nature of the state depends on the details of the short-range behaviour of the potential [24]. There is an energy E pt for which both core-halo and isothermal systems have the same entropy; above this energy the collapsed state becomes metastable and at some higher energy it ceases to exist. It is possible to regard the energy E pt as that where a true (first order) phase transition occurs. The value of E pt , and in general the form of the entropy-energy curve, is highly sensitive to the details of the short-range cutoff [58].

These results thus shows a rich variety of equilibrium states. However both astrophysical observation and numerical simulations show that, as we have discussed above, a self-gravitating system is rapidly trapped into a QSS (in astrophysics literature they more usually referred to as "collisionless equilibria"). While numerical simulations show that the these states evolve on time scales in very good agreement with the predictions of Chandrasekhar described above [59][60][61], whether the equilibria derived from the mean field limit of equilibrium thermodynamics are actually ever attained by such system remains, to our knowledge, an unanswered question. The Virial theorem and virialization.

As explained above a system starting from a phase space distribution which is not a time independent solution of the Vlasov equation undergoes a rapid evolution toward a Vlasov stable state (QSS).

One of the most noteworthy properties of such state is that they are so-called virial equilibria. This results is based on the "virial theorem" which we briefly discuss now for the self-gravitating case (but which can be generalised for any potential). We consider a system of self-gravitating particles in 3d interacting only through an exact 1/r potential. We introduce the inertia tensor [START_REF] Binney | Galactic Dynamics[END_REF]:

I µ⌫ = N X i=1 mx i,µ x i,⌫ , (1.38)
where x i,µ is the µ th component of the i th particle. The second time derivative of this quantity is:

Ï = m N X i=1 (ẍ i,µ x i,⌫ + x i,µ ẍi,⌫ + 2 ẋi,µ ẋi,⌫ ) . (1.39)
Inserting the equations of motion of the particles in the last expression, one easily obtains

Ï = 2m N X i=1 v i,µ v i,µ Gm 2 N X i=1 (x i,µ x i,µ )(x j,⌫ x i,⌫ ) |x i x j | 3 . (1.40)
Taking the trace of the last expression one recognises that the first term on the right hand side is four times the kinetic energy and the second twice the potential energy. One thus obtains the Lagrange identity:

1 2 Ï = 2K + U (1.41)
If now we assume that the system is, at a macroscopic level, in a stationary state, one obtains:

2K + U = 0 (1.42)
Thus, up to fluctuations due to finite particle number, we expect any such system to obey the relation (1.42). It can alternatively be directly derived from the Vlasov equation assuming only stationarity of f . We will illustrate this below for the case of a 1d self-gravitating system. It is straightforward also to generalise this result to any potential 1/r ↵ and to a system confined in a box, associated with a pressure:

2K + ↵U = 3P V 8↵ 6 = 0 (1. 43
)
where P is the pressure and V the volume of the system. For a system in thermal equilibrium at temperature T , K = 3 2 Nk B T , or more generally one can defined a kinetic temperature in this way (by 3 2 Nk B T = hv2 i) Then, if the pressure term can be neglected (e.g. in any QSS of an unbound selfgravitating system), we have, for the gravitational case:

2K + U = 0, thus E = K + U = K = 3 2 Nk B T , implying dE dT = k B
above. Thus when a self-gravitating system loses energy it heats up. By losing heat, the system grows hotter and continues to radiate energy. This is the mechanism that prevails in the internal region of stars [24]. For instance, at a stage where no more nuclear fuel is available, the core contracts and becomes hotter, giving its energy to the outer part which expands and becomes colder.

1d self-gravity

We introduce here briefly the 1d self-gravitating system, which is a paradigmatic model for long-range interacting systems. Besides the fact that it provides a natural case to study as a toy model for 3d gravity, this 1d model has the very nice feature that the particles trajectories can be analytically integrated between crossing which means that they can be simulated "exactly" as describe in detail in chapter 3.

Definition of the model

We consider a system of N particles of identical masses m = 1 moving in 1d and interacting through the gravitational pair interaction which obey the 1d Poisson equation:

@ 2 x (x) = 2g (x), (1.44)
with g the coupling constant. The model can be also seen as infinite self-gravitating sheets in three dimensions of surface mass density ⌃, which lead to the identification has g = 2⇡G⌃ with G the (3d) gravitational constant. Then, in a system of N such particles, the force exerted on a particle i by all the other is:

F i = g N X j6 =i sgn(x j x i ) , (1.45)
which can be written as

F i = g [N > (x i ) N < (x i )] (1.46)
where N > ( N < ) the number of particle on the left (right) of the particle i. Thus the force is constant other than when particles cross. We note further that compared to 3d gravity, 1d gravity differs notably in that (1) the force (and potential) is regular at r = 0 and thus no short-scale regularization is required and (2) the potential diverges at large separation and is therefore a confining potential. No box is therefore necessary to treat the equilibrium thermodynamics.

Virial theorem

In this section we derive the virial theorem for 1d gravity starting from the Vlasov equation. Integrating the 1d Vlasov equation over v, and using the assumption f (x, v, t) ! 0 for v, x ! ±1, one obtains the continuity equation:

@⇢(x, t) @t + @ (⇢(x, t)v(x, t)) @x = 0. (1.47)
Now, if we multiply the Vlasov-Poisson equation by x and v and integrate with respect the two variables over the whole space we obtain

Z xv @f @t dxdv + Z xv 2 @f @x dxdv Z xv @ @x @f @v dxdv = 0 . (1.48)
Considering first the integration over v, and integrating by parts the last term, this last equation is rewritten:

Z x @ (⇢(x, t)v(x, t)) @t dx + Z x @(⇢(x, t)v 2 (x, t)) @x dx + Z x⇢(x, t) @ @x dx = 0. (1.49)
Integrating by parts and using the continuity equation (1.47), the first term gives:

1 2 d dt Z x 2 @⇢(x, t)v(x, t) @x = 1 2 d 2 dt 2 Z x 2 ⇢(x, t) = 1 2
d 2 I(t) dt 2 .
(1.50)

with I the moment of inertia. The second term is

Z x @(⇢(x, t)v 2 (x, t)) @x dx = Z ⇢(x, t)v 2 (x, t) = 2K(t). (1.51)
Finally, the third term equals may be written as

Z x⇢(x, t) @ ¯ (x) @x dx = 1 M Z Z ⇢(x, t)gm x(x x 0 ) |x x 0 | ⇢(x 0 , t)dxdx 0 . (1.52)
Noting that the symmetric part, by exchange x $ x 0 , of

= x(x x 0 ) |x x 0 | is Sym = 1 2 |x
x 0 | and as ⇢(x, t)⇢(x 0 , t) is symmetric under this exchange, only the symmetric part survives in the integral and the third term equals the mean field potential energy:

1 2M Z Z ⇢(x 0 , t)gm|x x 0 |⇢(x, t)dx 0 dx = U (t). (1.53)
Collecting all terms one obtains the Lagrange indentity 1 2

d 2 I dt 2 = 2K(t) U (t). (1.54)
Hence for stationary states of the Vlasov equation, i.e. for QSS, we obtain the virial theorem for 1d gravity: 2K = U . In this case a virialised state verify E = 3K = 3 2 U = C ste . Further we note the specific heat is always positive, and indeed, it turns out that as we saw, the thermodynamic ensembles are equivalent.

Thermodynamics

As we have noted, this model does not present the small scale divergence problem; moreover, as we noted, it is not necessary to enclose the system in a box. The equilibrium thermodynamics has been completely solve by Rybicki in 1971 [20] and reproduced in greater detail in [START_REF] Worrakitpoonpon | Relaxation of isolated self-gravitating systems in one and three dimensions[END_REF]. Remarkably, the calculation give an exact expression, even for a finite N , of the equilibrium solution in both the canonical and micro-canonical ensemble. The Hamiltonian of the system is

H(x, v) = K(v) + U (x) = N X i 1 2
mv 2 i + gm 2 2 N X j=1 N X k6 =j |x j x k |; (1.55)
where x = x 1 , ...x N and v = v 1 , ..., v N . In this calculation, we start from the definition of the equilibrium distribution in the canonical ensemble

f c (x, v) = (zN !) 1 Z Z dx N dv N (x) (v)e H(v,x) 1 N N X i=1 (x x i ) (v v i ) , (1.56)
where (x) (v) is the constraint on the center of mass to remains fixed at the origin, and z the corresponding partition function. The factor N ! arises because we assumed indistinguishable particles. With this assumption the distribution function

1 N P N i=1 (x x i ) (v v i
) can be replaced by (p p N ) (x x N ) under the integral. Hence, an important consequence, is that the equilibrium solution is separable:

f c (x, v) = ⇢ c (x)⇥ c (v), (1.57)
with

⇢ c (x) = (QN !) 1 Z d d x (x)e U (x) (x x N ), (1.58)
and

⇥ c (v) = (R) 1 Z d d v (v)e K(v) (v v N ), (1.59)
where Q and R are the corresponding partition functions which normalise the two distributions with the relation z = QR. And then by direct calculation, using Fourier transform and complex analysis one can show that

⇥ c (v) = s N 2⇡m(N 1) e B mv 2 2(N 1) , (1.60)
and

⇢ c (x) = N gm 2 N X j=1 A N j e N gm 2 j|x| , (1.61)
with

A N j = j( 1) j+1 [(N 1)!] 2 (N 1 j)!(N 1 + j)! .
(1.62)

The calculation of the equilibrium solution in the micro-canonical ensemble is then obtained performing an inverse Laplace transform of the canonical solutions [20]. Then taking the limit N ! +1 at fixed mass M and energy E using the characteristic momentum and length scale: 2 = 4m 2 E 3M and ⇤ = 4E 3gM 2 , one obtains:

f eq (x, v) = 1 2 p ⇡ 1 ⇤ sech 2 ⇣ x ⇤ ⌘ e ( v ) 2 (1.63)

Dynamics

The dynamics of 1d self-gravitating system is qualitatively typical of that of any long-range interacting system as we have outlined it above. In particular, as we will detail and illustrate with results from our simulation in chapter 3, the system reaches first a given QSS on a time scale of order ⌧ dyn after a period of violent relaxation. The properties of the QSS depends on the initial condition of the system. Then on much longer times scales, proportional to N in units of ⌧ dyn , these QSS relax to the thermal equilibrium (given by Eq. (1.63) above). One notable feature is that this relaxation is very slow, in the sense that ⌧ relax ' zN ⌧ dyn where z is a numerical factor of order 10 2 10 3 which depends on the QSS [28]. The calculation of Chandrasekhar for 3d is not generalizable to this case as there are no real collisions in 1d because particles just cross on another and indeed the mechanism for this relaxation to thermal equilibrium remains an open problem for this system.

Chapter 2

Finite N corrections to Vlasov dynamics We saw in the previous chapter that one of the most interesting feature of the dynamics of a long-range interacting system it that they are trapped into QSS. A basic question is whether the appearance of these out of equilibrium stationary states -and more generally the validity of the Vlasov equation to describe the system's dynamics -applies to the same class of long-range interactions as defined by equilibrium statistical mechanics, or only to a sub-class of them, or indeed to a larger class of interactions. In short, in what class of systems can we expect to see these quasi-stationary states? Are they typical of long-range interactions as defined canonically? Or are they characteristic of a different class?

To answer these questions requires establishing the conditions of validity of the Vlasov equation, and specifically how such conditions depend on the two-body interaction. In the literature there are, on the one hand, some rigourous mathematical results establishing sufficient conditions for the existence of the Vlasov limit. It has been proven notably [38,[START_REF] Spohn | Large scale dynamics of interacting particles[END_REF][START_REF] Hauray | [END_REF] that the Vlasov equation is valid on times scales of order ⇠ logN times the dynamical time, for strictly bound pair potentials decaying at large separations r slower than r (d 2) . On the other hand, both results of numerical study and various theoretical approaches, based on different approximations or assumptions, suggest that much weaker conditions are sufficient, and the timescales for the validity of the Vlasov equation can be much longer. In the much studied case of gravity, notably, a treatment originally introduced by Chandrasekhar, [49] and subsequently refined by other authors (see e.g. [35,59,64,65] in which non-Vlasov effects are assumed to be dominated by incoherent two body interactions gives a time scale ⇠ N/logN times the dynamical time for the validity of the Vlasov equation, at least close to stationary solutions representing quasi-stationary states, and this in absence of a regularisation of the singularity in the two-body potential. Theoretical approaches in the physics literature derive the Vlasov equation and kinetic equations describing corrections to it (for a review see [1,66]) either within the framework of the BBGKY hierarchy [START_REF] Balescu | Equilibrium and Non-Equilibrium Statistical MechanicsWiley-Interscience publication[END_REF], as we saw in the previous chapter, or starting from the exact Klimontovich equation for the N body system [START_REF] Klimontovich | The Statistical Theory of Non-equilibrium Processes in Plasma˜International series of monographs in natural philosophy[END_REF][START_REF] Nicholson | [END_REF]. These approaches are both widely argued (e.g. see e.g. [1,2,66,69,70], to lead generally to lifetimes of quasi-stationary states of order ⇠ N times the dynamical time for any softened pair potential , except in the special case of spatially homogeneous quasi-stationary states one dimension.

CHAPTER 2. FINITE N CORRECTIONS TO VLASOV DYNAMICS

In this chapter, we address the question of the validity of Vlasov dynamics using an approach starting from the exact Klimontovich equation. Instead of considering, as is often done (see e.g. [1,2,69] , an average over an ensemble of initial conditions to define a smooth one particle phase density, we follow an approach (described e.g. in [71]) in which such a smoothed density is obtained by performing a coarse-graining in phase space. This approach gives the Vlasov equation for the coarse-grained phase space density when certain terms are discarded. We study how the latter "non-Vlasov" terms depend on the particle number N , on the scales introduced by the coarse-graining. In particular we develop this study analyzing the dependence on the large and small distance behavior of the two body potential. Our analysis leading to the scaling behaviours of these terms is based only one very simple -but physically reasonable -hypothesis that we can neglect all correlations in the (microscopic) N body configurations other than those coming from the mean (coarse-grained) phase space density.

The main physical result we highlight is that, under this simple hypothesis, the coarse-grained dynamics of an interacting N -particle system shows a very different dependence on the pair interaction at small scale depending on how fast the interaction decays at large distances: for interactions of which the pair force is integrable at large scales the coarse-grained N body dynamics is highly sensitive to how the potential is softened at much smaller scales, while for pair forces which are non-integrable the opposite is true. Correspondingly, while the Vlasov limit may be obtained for any pair interaction which is softened suitably at small scales, the conditions on the short-scale behaviour of the interaction are very different depending on whether its large scale behaviour is in one of of these two classes. This result provides a more rigorous basis for a "dynamical classification" of interactions as long-range or short-range, which has been introduced on the basis of simple considerations of the probability distribution of the force on a random particle in a uniform particle distribution in [72], and found also in [50] to coincide with a classification based on the dependence on softening of collisional time scales using a generalisation of the analysis of Chandrasekhar for the case of gravity.

The chapter is organized as follows. In the next section we derive the equation for the coarse-grained phase space density and write in a simple form the non-Vlasov terms our subsequent analysis focusses on. In section 2 we first explain our central hypothesis concerning the N -body dynamics, and then apply it to evaluate the statistical properties of the non-Vlasov terms. In the following section we then determine the scaling behaviours of these expressions, i.e., how they depend parametrically on the relevant parameters introduced, and in particular on the two parameters characterising the two body interaction -its large scale decay and the scale at which it is softened. In section 4 we use these expressions to identify the dominant contributions to the non-Vlasov terms, which turn out to differ depending on how rapidly the interaction decays at large scales. In the following section we present more complete exact results for the one dimensional case and the comparison with a simple numerical simulation. We then summarize our results and conclusions, discussing in particular the central assumptions and the dependence of our findings on them.

A Vlasov-like equation for the coarse-grained phase space density

We summarize a standard approach used to justify the validity of the Vlasov equation for long-range interacting systems alternative to the one described in the chapter 1 (BBGKY). The approach involves using a coarse-graining, in phase space, of the full N body dynamics and leads to an evolution equation for the coarse grained phase space density which consists of the Vlasov terms, plus additional terms. This equation, and the specific form of the non-Vlasov terms we derive, is the starting point for our analysis in the subsequent sections. We follow closely at the beginning the presentation and notation of [71].

We consider a d-dimensional system of N particles of identical mass m = 1 interacting only through the a generic two body force, denoting g(x) the force on a particle at x exerted by another one at the origin.

At any time t, the N particles have phase space positions which we denote {(x i , v i )} i=1..N , and the microscopic (or fine-grained, or Klimontovich) one particle phase-space density is simply the distribution

f k (x, v, t) = N X i=1 (x x i (t)) (v v i (t)). (2.1)
Likewise the microscopic one particle density distribution in coordinate space is

n k (x) = Z f k d d v = N X i=1 (x x i ), (2.2)
The full evolution of the N body system can be written in the form of the so-called Klimontovich equation for the microscopic phase space density:

@f k @t + v @f k @x + F[n k ](x) @f k @v = 0, (2.3)
where

F[n k ](x) = Z ⌦ g(x x 0 )n k (x 0 )d d x 0 = N X i=1 g(x x i ), (2.4)
is the exact force at point x (due to all particles). The detail derivation of this equation is given in [1], this equation is exact and contain the same information than the Hamilton's equation or the Liouville equation (cf. chapter 1). The only assumption made in deriving this equation from the equations of motion of the individual particles is that the force g(x) is bounded as

x ! 0. Introducing a top-hat window function W (z = z 1 , . . . , z d ), W (z) = ⇢ 1, if |z| < 1, 0, otherwise, (2.5)
we define the coarse-grained phase space density:

f 0 (x, v, t) = Z d d x 0 d x d d v 0 d v W ✓ x x 0 x ◆ W ✓ v v 0 v ◆ f k (x 0 , v 0 , t), (2.6)
where x and v are the characteristic sizes of the coarse-graining cell in position and, respectively, velocity space. We will denote by C(x, v) the coarse graining cell centred at (x, v), which thus has a phase space volume d

x d v , and

N c (x, v, t) = ⇡ d ( (d/2 + 1)) 2 d x d v f 0 (x, v, t), (2.7)
is the number of particles in C(x, v). We suppose that the coarse graining cell is always much smaller (in both real and velocity space) than the characteristic size of the system, but sufficiently large to contain a large number of particles, i.e., x ⌧ L x and v ⌧ L v where L x and L v the characteristic size of the system in coordinate and velocity space respectively, and

1 ⌧ N c (x, v, t) ⌧ N . (2.8)
Likewise we define the coarse-grained spatial density as

n 0 (x, t) = Z d d vf 0 (x, v, t) = Z d d x 0 x W ✓ x x 0 d x ◆ n k (x 0 ). (2.9)
To obtain a Vlasov-like equation we integrate Eq. ( 2.3) over the coarse-graining cell over x and v. We obtain, for the first term;

Z Z d d x 0 d d v 0 @f k @t W ✓ x 0 x x ◆ W ✓ v 0 v v ◆ = @f 0 (x, v, t) @t . (2.10)
For the second term using the Einstein summation convention for repeated index, for j = 1 . . . d, we have:

Z Z d d x 0 d d v 0 d x d v v 0 j . @f k @x 0 j W ✓ x 0 x x ◆ W ✓ v 0 v v ◆ (2.11) = Z Z d d x 0 d d v 0 d x d v v j . @f k @x 0 j W ✓ x 0 x x ◆ W ✓ v 0 v v ◆ + Z Z d d x 0 d d v 0 d x d v (v 0 j v j ). @f k @x 0 j W ✓ x 0 x x ◆ W ✓ v 0 v v ◆ = Z Z d d x 0 d d v 0 d x d v f k v j . @ @x 0 j W ✓ x 0 x x ◆ W ✓ v 0 v v ◆ Z Z d d x 0 d d v 0 d x d v f k (v 0 j v j ). @ @x 0 j W ✓ x 0 x x ◆ W ✓ v 0 v v ◆ .
To obtain the last equality we have used an integration by parts, in which the term

h f k W ⇣ x 0 x x ⌘ W ⇣ v 0 v v ⌘i +1 1 vanishes. Then noting that @ @x 0 j W ✓ x 0 x x ◆ = @ @x j W ✓ x 0 x x ◆ , (2.12)
CHAPTER 2. FINITE N CORRECTIONS TO VLASOV DYNAMICS the derivatives @ @x i can be taken out of the integrals and we obtain for the second term

v j . @f 0 @x j + @ @x j Z Z C(x,v) d d x 0 d d v 0 d x d v (v 0 j v j )f k = v @f @x + @ @x 2 4 1 N c (x, v, t) X i2C(x,v) v i (t) v 3 5 .
(2.13) Then following the same procedure, adding F 0 (x, t) F 0 (x, t) under the integral and again integrating by parts the third term gives:

F 0 (x) • @f 0 @v + @ @v 2 4 1 N c (x, v, t) X i2C(x,v) F i (t) F 0 (x, t) 3 5 .
(2.14)

Collecting all expressions derived above, we obtain the following equation:

@f 0 @t + v • @f 0 @x + F 0 (x) • @f 0 @v = @ @x [f 0 ⇠ v ] @ @v [f 0 ⇠ F ] , (2.15)
where

F 0 (x, t) = Z ⌦ g(x x 0 )n 0 (x 0 , t)d d x 0 , (2.16)
is the force at the point x due to the coarse-grained distribution (which we will identify as the mean-field force). Furthermore

⇠ v (x, v, t) = 2 4 1 N c (x, v, t) X i2C(x,v) v i (t) 3 5 v , (2.17)
where v i is the velocity of particle i, and

⇠ F (x, v, t) = 2 4 1 N c (x, v, t) X i2C(x,v) F i (t) 3 5 F 0 (x, t) , (2.18)
where F i is the exact force acting on the particle i, i.e., F i = F(x i , t) where we simplify the above notation:

F(x, t) . = F[n k ](x) = Z ⌦ g(x x 0 )n k (x 0 , t)d d x 0 . (2.19)
Thus ⇠ v is the "velocity fluctuation" in the cell (x, v) around the coarse-grained velocity v, i.e., difference between the arithmetic mean of particles velocities in the cell (i.e. the velocity of the center of mass) and the velocity at the centre of the coarse-graining cell, and ⇠ F (x, v, t) is the "force fluctuation" around the coarsegrained force, i.e. the difference between the arithmetic mean of the exact forces acting on each particle in the cell (equal to the force on the centre of mass of the particles in the cell) and the force at the centre of the cell due to the coarse-grained particle distribution.

If the right hand side of Eq. (2.15) is set equal to zero, we recover, given Eq. (2.16), the Vlasov equation for the coarse-grained phase space density f 0 (x, v, t). 29

Establishing the validity of the Vlasov equation in an appropriate limit thus requires showing that the terms on the right-hand side may indeed be taken to zero in this limit. For a real system, for which N is finite and the typical number of particles in a coarse-graining cell is finite, Eq. (2.15) is not closed for the coarse-grained phase space density, but rather coupled to the fine-grained density through the terms on the right-hand side. If it is possible to define the Vlasov limit for the system, these terms will represent at any finite (but large) N , small perturbations to the pure Vlasov evolution of the coarse-grained distribution associated with the "graininess" of the system which, under suitable hypothesis, are responsible for the relaxation of the system to the thermodynamic equilibrium. In the rest of this article we focus on these terms, and establish their scaling with N (or, equivalently, as a function of the characteristic scales of the coarse-graining cell), given certain simplifying hypotheses.

Statistical evaluation of the finite N fluctuating terms

We now focus on the two terms on the right-hand side of Eq. (2.15). Their direct evaluation is clearly impossible as in principle it requires knowledge of the full fine-grained phase space density. We can, however, determine how they scale as a function of relevant parameters by using a statistical approach: given a coarsegrained distribution f 0 (x, v) we can consider in each cell of the single particle phase space, the fine-grained distribution to be characterized by an ensemble of realizations of particle distributions having f 0 (x, v) as constant mean density. If we know the statistical properties of this ensemble, we can then, in principle, calculate those of the fluctuating terms on the right-hand side of Eq. (2.15). In particular, we can then consider how the amplitudes of these statistical quantities depend on the relevant parameters.

We make here the most simple possible hypothesis about this ensemble for the fine-grained phase space density: we suppose that it corresponds to the ensemble of realisations of an inhomogeneous Poissonian point process with mean density given by the coarse-grained phase space density f 0 (x, v). In other words we assume that the particles are randomly distributed, without any correlation, inside each coarsegrained cell, with a mean density which varies from cell to cell following f 0 (x, v). The density-density correlations are thus fully described by the one point distribution f 0 (x, v), and all other correlations, of the fluctuations around this mean density, are neglected. Physically this means we retain the "pure" finite N (discreteness) effects arising from the fluctuations of the mean density, but neglect the correlation of these fluctuations. Alternatively one can consider that we proceed by assuming complete ignorance of the distribution of particles below the coarse-graining scale, other than the information furnished about it by the coarse-grained density itself. This allows us to close, at least in a statistical formulation, Eq. (2.15) for f 0 (x, v). Indeed this hypothesis is arguably the most natural one to make in seeking to obtain a criterion for the validity of the Vlasov equation which is based only on the coarse graining density f 0 (x, v) and the properties of the pair interaction itself. Formally we can state our assumption to be that the relevant terms can be evaluated by considering an ensemble realisations of a point process with the N CHAPTER 2. FINITE N CORRECTIONS TO VLASOV DYNAMICS particle probability distribution in phase space given by (see, e.g., [73])

P N (x 1 , v 1 ; ..; x N , v N ) = N Y i=1 f 0 (x i , v i ) N , (2.20)
assuming the coarse-grained phase space density f 0 (x, v) to be a smooth function. In practice it is convenient to perform the calculation with finite coarse-grained cells in which f 0 (x, v) is fixed, and the particles are then distributed randomly in each coarse-grained cell.

Mean and variance of ⇠ v

We first evaluate the mean and variance of ⇠ v , as defined by Eq. (2.17) above, assuming now the defined properties of the ensemble of realisations. Given that the latter assigns randomly the velocities of the particles in the coarse-grained cell, which is centred at (x, v), we evidently have *

1 N c (x, v, t) X i2C(x,v) v i (t) + = v , (2.21)
where h...i denotes the ensemble average. Therefore we have, as would be expected,

h⇠ v i = 0 . (2.22)
Calculating the variance of ⇠ v gives straightforwardly that

⌦ ⇠ 2 v ↵ = 1 N c (x, v) d 12 2 v , (2.23)
where we can take

N c = f 0 (x, v) d x d
v to be the average number of particles in the cell, given that N c = hN c i + N c where by hypothesis

N c ⇠ p N c ⌧ N c 1 . Given that N c (x, v
) is large, and v i are considered independent and identically distributed variables, ⇠ v is thus simply a Gaussian random variable with mean zero and variance given by Eq. (2.23).

Mean and variance of ⇠ F

Let us evaluate now the first two moments of ⇠ F . To evaluate these quantities in the ensemble defined above, it is useful first to note that

F(x, v, t) . = 1 N c (x, v, t) X i2C(x,v) F i (t) (2.24) = 1 N c Nc X i=1 N Nc X I=1 g(x i y I )
where the i runs over the N c particles inside the coarse-graining cell C(x, v), and the index I over the other N N c particles outside the cell. This is the case because we assume g( x) = g(x) (making the sum over all the mutual forces of the pairs of particles in the cell C(x, v) vanishes). As the individual pair forces depend only on the spatial positions of particles, we need only specify, to calculate the ensemble average of powers of the force, the probability distribution for the spatial positions of these points. Given the writing of the force in (2.24), in which the sum is performed separately over the particles inside and outside the coarse-grained cell C(x, v) considered, it is convenient to write the ensemble average in a similar form. As noted above we can take N c to be fixed and equal to its average value, n 0 (x) d x , up to negligible corrections of order 1/ p N c . We can then write the N particle probability distribution as

P 0 N (x 1 , ..x Nc ; y 1 , ..y N Nc ) = Nc Y i=1 p(x i ) N Nc Y I=1 p(y I ), (2.25)
where p(x) is the one-point probability distribution function of the spatial position of a particle given the condition that it is contained in the cell, and p(y) is the onepoint probability distribution of the spatial position of a particle given the condition that it is outside the cell. Given that the N c particles are randomly distributed inside the cell, we have simply

p(x i ) = ⇢ d x , x i 2 C(x, v) 0, otherwise. (2.26)
We define further S(x), the support of W (x/ x ), i.e., the "stripe" in phase space with the same spatial coordinates as the phase space cell C(x, v) (see Fig. 2.1). Assuming the number of particles in the cell N c to be small compared to the total number of particles N , p(y) can be approximated everywhere simply by is the unconditional one point PDF for the spatial distribution obtained by integrating (2.20) over all but one space coordinate, i.e.,

p(y) = n 0 (y) N . (2.27)
Mean of ⇠ F Using (2.25) to calculate the ensemble average of the exact force exerted by all particles on those in a coarse-grained cell, we have

⌦ F(x, v, t) ↵ = (N N c ) Z ⌦ d d x 0 Z ⌦ d d y 0 p(x 0 )p(y 0 )g(x 0 y 0 ) = Z S(x) d d x 0 d x Z ⌦ d d y 0 n 0 (y 0 )g(x 0 y 0 ) = Z S(x) d d x 0 d x F 0 (x 0 ), (2.28)
where we have approximated N Nc N by 1. Thus

h⇠ F (x, v, t)i = Z S(x) d d x 0 d x [F 0 (x 0 ) F 0 (x)] . (2.29)
Assuming that we can neglect the variation of the mean-field F 0 in the coarsegraining cell, we obtain h⇠ F (x, v, t)i = 0.

(2.30)

Figure 2.1: Scheme of the coarse-graining of phase space (in 1 dimension). Illustration of the different regions corresponding to the domains of integration used our analysis of the integrals: C is the coarse-grained cell with centre at (x, v), S (in white) is the "stripe" enclosing all points in [x x 2 , x + x 2 ]; @ " S is the region containing particles within a distance " of the boundary of the stripe.

Variance of ⇠

F We now calculate ⌦ ⇠ 2 F (x, v, t) ↵ = ⌦ F2 (x, t) ↵ F 2 0 (x, t). (2.31) We first break h⇠ 2 F (x, v, t)i into two parts ⌦ F2 (x, t) ↵ = * 1 N 2 c X i F 2 (x i ) + + * 1 N 2 c X i,j;j6 =i F(x i ) • F(x j ) + , (2.32)
where all the sums are over the N c particles in the cell C(x, v). We will refer to the first term on the right hand side as the "diagonal" contribution, and the second term as the "off-diagonal" contribution to the variance: the first is the contribution to the variance due to the variance of the force on each particle of the cell, the second the contribution to the variance arising from the correlation of the forces on different particles in the cell. Further we can use Eq. (2.24) to split each term into two terms

⌦ F2 (x, t) ↵ = * 1 N 2 c Nc X i N Nc X I g 2 (x i x I ) + + * 1 N 2 c Nc X i N Nc X I,J;I6 =J g(x i x I ) • g(x i x J ) + + * 1 N 2 c Nc X i,j;i6 =j N N c X I g(x i x I ) • g(x j x I ) + + * 1 N 2 c Nc X i,j;i6 =j N N c X I,J;I6 =J g(x i x I ) • g(x j x J ) + ,
where i, j again denote sums over the N c particles inside the coarse-grained cell and I, J over the remaining N N c ⇠ N particles . Performing the ensemble average by integrating over the PDF (2.25), the result can be conveniently divided in two parts. The first and third terms give

1 N c Z ⌦ d d x 0 d x n 0 (x 0 ) Z S(x) d d y g 2 (y 0 x 0 ), (2.33)
and ✓

1 1 N c (x, v) ◆ Z ⌦ d d x 0 d x n 0 (x 0 ) Z S(x) d d y 0 g(y 0 x 0 ) 2 , (2.34)
respectively. Both the second and fourth terms can be expressed purely in terms of the mean-field, as

1 N c ✓ 1 1 N ◆ Z S(x) d d x 0 d x F 2 0 (x 0 ), (2.35)
and ✓

1 1 N c ◆✓ 1 1 N ◆Z S(x) d d x 1 d x Z S(x) d d x 2 d x F 0 (x 1 ) • F 0 (x 2 ) (2.36)
respectively.

Assuming that we can neglect the variation of the mean-field F 0 in the coarsegraining cell, we can perform the integrals in the last two expressions, and then obtain

⌦ ⇠ 2 F (x, v) ↵ = 1 N c (x, v) Z ⌦ d d x 0 n 0 (x 0 ) Z S(x) d d y 0 d x g 2 (y 0 x 0 ) + (1 1 N c (x, v) ) Z ⌦ d d x 0 n 0 (x 0 ) Z S(x) d d y 0 d x g(y 0 x 0 ) 2 1 N F 2 0 (x).
(2.37)

Our analysis below will focus essentially on the first two terms in this expression as they are those which describe the contribution to the fluctuating terms which are potentially sensitive to the small scale properties of the pair force g(x). We note that the first term comes from the "diagonal" part of the variance, and more specifically it represents the contribution to the variance arising from the force on a single particle in the cell due to a particle outside the cell: we will thus refer to it as the two body contribution. The second integral in (2.37), on the other hand, arises from the "off-diagonal" part of the variance, and more specifically it represents the contribution to the variance of the force on the cell due to the correlation of the force exerted on two particles inside the cell exerted by a particle outside the the cell: we will refer to it therefore as the three body contribution.

Parametric dependence of the fluctuations

We now analyse the expressions we have obtained, focussing on how their value depends parametrically on the relevant parameters we have introduced, notably the number of particles in the system (N ), its size (L x , L v ), the number of particles in the coarse-graining cells (N c ) and the coarse-graining scales ( x , v ). Further we will take the pair force to be given by g(x) = g|x| x for < 0, and by

g(x) = g ⇢ |x| x, |x| " " x, |x| < " , (2.38)
for > 0, where g is the coupling constant (and g < 0 for the case of an attractive interaction). We will focus then on how the amplitude of the fluctuations of the force depend also on the exponent and the characteristic length " at which the force is regularized at small scales. Indeed it is evident that we must introduce such a regularization of the pair force, as without it the integrals in (2.37) can be ill-defined.

Mean field Vlasov limit

The mean-field Vlasov limit is formulated by taking N ! 1 at fixed system size, and scaling g / N 1 so that the mean field F 0 remains fixed. Applying this procedure to the expressions we have obtained for h⇠ 2 F i and h⇠ 2 v i, both indeed converge to zero, for any non-zero ". In order to obtain these expressions we have, as noted, assumed also that the mean-field does not vary on the scale of the coarse-graining cell. Assuming that the characteristic scale of variation of the coarse-grained phase space density, and mean field, is the system size L x , for a finite coarse-graining cell we expect corrections to our expressions due to the variation of these quantities on the scale x which are suppressed at least by x /L x relative to those we have calculated. The coarse-grained phase space density f 0 (x, v, t) thus indeed obeys the Vlasov equation in the usual formulation of this limit, when the size of the coarse-graining cell is taken to be negligible with respect to the system size.

We now study more closely the approach to the Vlasov limit as characterized by the scaling behaviour of the corrections to it. We assume that these are given by those of the statistical quantities we have calculated, i.e. we take

|⇠ v | ⇠ p h⇠ 2 v i , |⇠ F | ⇠ q h⇠ 2 F i. (2.39) CHAPTER 2. FINITE N CORRECTIONS TO VLASOV DYNAMICS

Velocity fluctuations

Given the result (2.23) and that

N c ⇠ N L d x L d v d x d v we infer |⇠ v | ⇠ 1 p N ✓ L x x ◆ d/2 ✓ L v v ◆ d/2 v . (2.40)
There is evidently no dependence on the pair force in this term. As already noted we can recover the Vlasov limit by taking N ! 1. Further we can see that this result remains valid for arbitrarily small (but non-zero) values of the ratio x Lx so that variation of coarse-grained quantities on the coarse-grained scale can indeed be neglected.

Force fluctuations

The scaling of the last term on the right-hand side (2.37) is already explicited, representing simply a fluctuation of the force about the mean field of order 1/ p N times the mean field itself. In order to determine the dependences of the first two terms we need to analyse carefully that of the integrals. To do so we divide the domains of integration in the double integral into appropriate subdomains, isolating the region which may depends on the lower cut-off " in the pair force.

body contribution

We consider first the diagonal two body contribution, writing the double integral as

Z ⌦/S/@"S d d x 0 + Z @"S d d x 0 + Z S/@"S d d x 0 ⌦ ✓ n 0 (x 0 ) Z B(x 0 ,") d d y 0 d x g 2 (y 0 x 0 ) (2.41) +n 0 (x 0 ) Z S(x, )/B(x 0 ,") d d y d x g 2 (y 0 x 0 ) ◆ ,
where in this context ⌦ indicates the integration operation on the terms in parenthesis, and as illustrated for the one dimensional case in Fig. 2.1, the integral over

x 0 , over the whole of space (⌦) has been divided into three parts:

• @ " S: the set of points which are within a distance " of the boundary of the stripe S. The volume of this region is of order " d 1 x . • S/@ " S: the set of points belonging to S but not belonging to @ " S. For " ⌧ x , its volume is of order d x .

• ⌦/S/@ " S: the rest of space. and the integral over y 0 , over the stripe S has been divided into two parts:

• B(x 0 , "): the set of points in S which are a distance of less than " from the point x 0 . This region has a volume of order " d .

• S/B(x 0 , "): the rest of S; for " ⌧ x , its volume is of order d x .

We consider now one by one the terms in the integral written as in (2.41). The first term in the integration over x 0 excludes the region where the pair force is "dependent, and is over a volume of order L d x . Hence if we suppose n 0 ⇠ const, it gives a contribution which scales as

⇠ n 0 L d 2 x . (2.42)
For the second region of integration over x 0 , of volume of order " d 1 x , the region B(x 0 , ") of the integration over y 0 gives a contribution of order " 2 over a volume of order " d , and thus

⇠ n 0 (" d 1 x )" d d x " 2 .
(2.43)

In the region S(x)/B(x 0 , ") on the other hand, of volume of order d x , we have

⇠ n 0 (" d 1 x ) 2 x .
(2.44)

In the third region of integration over x 0 , of volume of order d x , we obtain again a contribution of order " 2 in the volume B(x 0 , "), and thus

⇠ n 0 " d 2 , (2.45)
while the second term region of the y 0 integration gives

⇠ n 0 d 2 x
.

(2.46)

body contribution

Proceeding in the same manner, we write the 3-body contribution to the variance of the force as Z ⌦/S/@"S

d d x 0 + Z @"S d d x 0 + Z S/@"S d d x 0 ⌦ ✓ n 0 (x 0 ) Z B(x 0 ,") d d y 0 d x g(y 0 x 0 ) (2.47) +n 0 (x 0 ) Z S(x, )/B(x 0 ,") d d y d x g(y 0 x 0 ) ◆ 2 .
The first term in the integration over x 0 gives then a contribution

⇠ n 0 L d 2 x .
Compared to the 2-body integral, the analysis of the remaining parts is essentially the same, except for one important difference: as the integration over y 0 is over the vector pair force, the integral over y 0 is zero when integrated in a sphere around x 0 ; in particular the integration over B(x 0 , ") vanishes when B(x 0 , ") is fully contained in S. This is the case for the integration region S/@ " S, which therefore does not depend on " and simply gives a contribution of order

⇠ n 0 d 2 x . (2.48)
For the second integration region in the integral over x 0 , the volume ⇠ " d 1 x , B(x 0 , ") is not fully contained in S and we have therefore a contribution from this CHAPTER 2. FINITE N CORRECTIONS TO VLASOV DYNAMICS part of the integration over y 0 of order " d d

x " , while the second part, which does not depend on ", gives a contribution of order x . Taking the square and multiplying by the volume of @ " S, we obtain three terms:

⇠ n 0 " d 1 x " 2d 2d x " 2 , (2.49)
from the square of the first term,

⇠ n 0 " d 1 x " d d x " x , (2.50)
from the cross term, and

⇠ n 0 " d 1 x 2 x , (2.51)
from the square of the second term.

4 Force fluctuations about the Vlasov limit: dependence on "

Gathering together the expressions derived above, we obtain, keeping only the leading divergence in " in each of the 2-body and 3-body contributions,

⌦ ⇠ 2 F (x, v) ↵ = g 2 N c C ⌦ n 0 L d 2 x + C S n 0 d 2 x + C @S n 0 " d 2 (2.52) + g 2 (1 1 N c ) C 0 ⌦ n 0 L d 2 x + C 0 S n 0 d 2 x + C 0 @S n 0 " 2d+1 2 (d+1) x 1 N |F[n 0 ](x)| 2 ,
where all C ⇤ and C 0 ⇤ are constants (Note that we have not included (2.50) because when it diverges, for > d + 1, the term retained is indeed more rapidly divergent.)

Depending on the values of and d different terms dominate. We consider each case.

Case < d 2

For this case, there are no divergences as " ! 0, and for " ⌧ x ⌧ L x , the dominant term from the two integrals is

g 2 n 0 L d 2 x ⇠ 1 n 0 L x g 2 n 2 0 L 2 x ⇠ 1 N |F 0 | 2 , (2.53)
and therefore we infer the scaling of the total force fluctuation is

|⇠ F | ⇠ 1 p N |F 0 |. (2.54)
As noted above we therefore obtain in this case the Vlasov limit taking N ! 1 with g ⇠ 1/N . Further we conclude that, at finite N , the fluctuations around the mean-field force are dominated by contributions coming from fluctuations of the density at the scale of the system size, which dominate those coming both from the scale x of the coarse-graining cell and those from the scale " at which the pair force is regularized.

Case

d 2 < < d + 1 2
In this range of , there is a divergence at " ! 0 in the contribution coming from the 2-body term, while the 3-body term remains finite. Keeping only the dominant contributions to the two integrals when " ⌧ x ⌧ L x , we obtain

h⇠ 2 F (x)i ⇠ C " g 2 n 0 (x) N c (x, v) " d 2 + C 0 g 2 n 0 (x) d 2 x , (2.55)
where C " and C 0 are constants. Given the divergence in " we see explicitly that in this case the Vlasov limit is obtained taking N ! 1 with g ⇠ N 1 at finite non-zero ", while the limit does not exist in the absence of this small scale regularisation of the pair force.

Using

N c (x, v) / f 0 (x, v) d x d
v , we can write the dominant fluctuations as

h⇠ 2 F (x)i ⇠ g 2 n 0 (x) d x ✓ C " " d 2 f 0 (x, v) d v + C 0 2d 2 x ◆ . (2.56)
This expression allows us to conclude, as anticipated, that there is a crucial difference between the following sub-cases: i) the range of in which the first term dominates, and ii) the range in which the second term dominates:

For d 2 < < d:
In this case the exponent of x in the second term inside the brackets in (2.56) is positive, and therefore when we increase x at fixed " (and fixed v ), this term dominates over the first one. More specifically when x "(f

0 (x, v) d v " d ) 1 2d 2
this term dominates, and

h⇠ 2 F i ⇠ g 2 n 0 (x) d 2 x (2.57)
Thus, even though the amplitude of the fluctuations depends on ", and diverges as " ! 0, for a sufficiently large coarse-graining cell, fluctuations become in practice effectively insensitive to the value of ", for a wide range of values which is such that the larger is x the smaller is the lower limit on ", with the latter vanishing as

x ! 1. Note that since

⇠ F ⇠ 1 p n 0 (x) d x ✓ x L x ◆ d |F 0 (x)| ⇠ 1 p N ✓ x L x ◆ d/2 |F 0 |, (2.58)
we can also neglect the final term in (2.37).

For d < < d + 1 2 : In this range, it is instead the exponent of x in the second term inside the brackets in (2.56) which is now negative, and as a consequence it is now the first term which dominates when we make the coarse-graining scale x large. We have, therefore

h⇠ 2 F i ⇠ g 2 n 0 f 0 d x d v " d 2 .
(2.59)

Further as this can be rewritten as

|⇠ F | ⇠ 1 p NN c ✓ " L x ◆ d/2 |F 0 |. (2.60)
It follows that for sufficiently small " this is the dominant contribution to the fluctuations. In this range therefore the leading contribution to the force fluctuations is directly dependent on ".

Case d + 1 2 <

In this case both the integrals giving the 2-body and 3-body contributions are divergent at small ", but the dominant divergence comes from the latter giving

h⇠ 2 F (x)i ⇠ g 2 n 0 (x) " 2d+1 2 d+1 x . (2.61)
This case is therefore like the previous case (d < < d + 1 2 ): the dominant contribution to the fluctuations is divergent as " ! 0.

Exact one dimensional calculation and numerical results

For a one dimensional system, d = 1, it is possible to perform explicitly the integrals in (2.37) to obtain exactly the expression of the force fluctuations. In order to illustrate our main result above, we can compare the expressions we obtain, and in particular their leading scaling behaviours, with what is obtained directly by measuring the force fluctuations in cells on realisations of a homogeneous Poisson particle distribution. As we are interested primarily in the " dependence of these fluctuations we consider solely the contribution to them from the part of the integral which is potentially sensitive to them. We can therefore write

h⇠ 2 F (x)i S = n 0 N c Z x 2 x 2 dx 0 Z x 2 x 2 dy 0 1 x g 2 (y 0 x 0 ) +n 0 (1 1 N c ) Z x 2 x 2 dx 0 " Z x 2 x 2 dy 0 1 x g(y 0 x 0 ) # 2 ,
where we use the subscript in h⇠ 2 F (x)i S to indicate that this is the contribution to the force variance sourced by particles in the phase space "stripe" S, and

g(x) = g ( x |x| +1 , |x| " x |x| " |x| < ". (2.62)
Integrating we obtain

h⇠ 2 F (x)i S = g 2 1 N c n 0 (x) x  1 1 " 2 2 4 1 2 x " 1 2 + 1 (1 2 )(1 ) 2 2 x
(2.63) 1)(3 2 )

+g 2 (1 1 N c ) n 0 (x) 2 x  2 (1 + 2 ) 2(
" 3 2 4 (1 ) 2 (2 )
" 1 ( 2 x ( x ") 2 ) 4 (1 ) (2 )(3 )
" ( 3 x ( x ") 3 ) + 4 (1 ) (2 )
" 1 ( x ") 2 4 (1 ) 2 ✓ x 2 ◆ 2 2 ✓ x 2 " ◆ 1 X n=0 ( 1) n 1 + 2n 1 n! ✓ 1 2 " x ◆ 2n + 2 (1 ) 2 (3 2 ) 3 2 x
, where (x) n is the Pochhammer symbol, (x) n = (x+n) (x) = x(x + 1)...(x + n 1). Expanding this expression in the limit "/ ! 0, and keeping only the leading "-dependent terms, we obtain:

h⇠ 2 F (x)i S = g 2 n 0 (x) N c  4 1 2 " 1 2 + 1 (1 2 )(1 ) 1 2 x
(2.64)

+g 2 n 0 (x)(1 1 N c )  2 (1 + 2 ) 2(1 )(3 2 ) 2 x " 3 2 + 4 (1 ) 2 1 x " 2 + C 1 2 x ,
where C is a constant depending on . Comparing this expression with (2.52) which we obtained for the d-dimensional case, we note that we indeed have agreement when we take d = 1. The term of order " 2 corresponds to the term of order " d+1 in (2.50) which we did not include in (2.52) because it is never the leading divergence. We now compare these analytical expressions with those obtained from a direct numerical estimation of the same quantity in a Poissonian realisation of a particle system. To do so we have distributed N = 72900 particles randomly in a "phase space box" (cf. illustration in Fig. 2.1) of side 27 v (in velocity) and 9 x (in position). Thus there are 243 cells, containing, on average, 300 particles. We then calculate the exact force on each particle in the system due to all particles in the stripe to which it belongs, i.e., F(x i ) = P j2S g(x i x j ). For each cell we then average this quantity to get F = 1

Nc

P

i2C F[n k ](x i ), and we finally estimate the variance of this quantity over the 243 cells.

Shown in Fig. 2.2 are the results of these numerical simulations compared with the theoretical results, Eq. (2.63), for a range of values of " and different values of . The agreement is excellent in all cases. The crosses indicate the points obtained from the numerical simulation described in the text and dotted lines are the corresponding analytical results, Eq. (2.63). We use units in which x = 1 and g = 1.

Further we have verified the validity of the expanded expression of the fluctuation Eq. (2.64). In Fig. 2.3 (upper panel) is shown, for example, the case = 1.25. In addition to the simulation results and the exact theoretical results Eq. (2.63), we have plotted for comparison the dominant divergence of the 2-body term which scales as " 1 2 . The perfect matching of the two curve at low " (full and dominant divergence) confirms the expression of the dominant power.

According to Eq. (2.64), the dominant power of the 3-body term is expected to be constant in this case and we have plotted also the full 3-body contribution to check it.

Except in the sector where " ⇠ x , it is well verified that the 3-body contribution is a constant, and present no divergence for = 1.25 < d/2 + 1.

In the lower panel of Fig. 2.3 we display the same quantities for the case = 1.75. The expected dominant power of the 2-body contribution matches with the full expression and the numerical results. The 3-body contribution is also plotted rescaled by the expected leading divergence " 3 2 . The curve obtain, except in the region where " ⇠ x , is a straight line and validate the expansion Eq. (2.64) for the 3-body contribution as well.

For the others cases, the results are in excellent agreement with the expression (2.64) for the small "/ x behaviour.

Thus we find results completely in line with our determination of the scalings in this limit for the d-dimensional case, leading to the different behaviours described in Section 4: for = 0.25, in the range < d/2, the force fluctuations are independent of "; for = 0.75 and = 1.25, in the range d/2 < < d/2 + 1, the divergent 2-body term ⇠ " 1 2 dominates at small ", but is then overtaken at larger " by the flat behaviour of the the dominant "-independent term in the 3-body term; for = 1.75 and = 2.5, both in the range d/2 + 1 < we have again at the smallest ✏ Let us now discuss further the physical significance of these results, and in particular how they justify the basic qualitative distinction between interactions as we have anticipated in the introduction. We have considered a generic N -particle system with Hamiltonian dynamics described by a two body interaction with a pair force ⇠ 1/r and regularized below a scale " to a constant value. Introducing a coarse-graining in phase space, one obtains an equation for the coarse-grained phase space density with terms corresponding to the Vlasov equation, in addition to "non-Vlasov" terms which are functionals of the microscopic phase density. We then made the hypothesis that the typical amplitude of these terms can be estimated by assuming this microscopic phase space density to be given by a realisation of an inhomogeneous Poissonian point process, in which the mean density is specified by the coarse-grained phase space density but the density fluctuations are uncorrelated. In other words we neglect all contributions due to correlation of density fluctuations with respect to the ones due to the simple products of mean densities. Doing so we have determined the scalings of these latter terms as a function of the different scales in the problem, of and of the spatial dimension d. If the assumed hypothesis is valid (for large N ), a limit in which the Vlasov equation applies is simply any one in which the derived non-Vlasov terms go to zero while the Vlasov terms tend to fixed finite values.

Our results show firstly that, for any interaction in the class we have considered, such a limit exists and can be defined in different ways: notably, by taking N to infinity at fixed values of the other scales, and in particular at fixed smoothing; or, alternatively, by taking the limit in which the coarse-graining scale, while always remaining small compared to the system size, is taken arbitrarily large compared to the other scales, and in particular the smoothing scale ". However, when this limit is taken, there is a important difference between the case in which > d or < d: for the former case the dominant correction term to the Vlasov limit is always strongly dependent on ", while in the latter case the dominant term, at fixed x , is independent of " in a wide range whose lower limit vanishes as x diverges. This different leading dependence on " corresponds to fluctuations which are sourced by quite different contributions in the two cases: for > d, the dominant fluctuation in the force on a coarse-graining cell comes from the contributions coming from forces on individual particles due to single particles which are "closeby", i.e., within a radius ": for < d, on the other hand, the dominant fluctuation comes from the coherent effect on two different particles in the cell coming from a single particle which is "far-away", i.e. at a distance of order the size of the coarse-graining cell, or, even for < d/2, of order the size of the system. Given that these amplitudes are those of the terms describing the corrections to the Vlasov evolution due to the large, but finite, number of particles in a real system, this means that the dynamics of the particle system at a coarse-grained scale in the two cases is dominated by completely different contributions: by the particle distribution at the smallest scales when > d, and by the particle distribution at the coarse-graining scale or larger when < d. Thus in the latter case we can decouple the dynamics at the coarse-graining scale from that at smaller scales (the interparticle distance, the scale of particle size), while for > d we cannot do so. Or, in the language of renormalisation theory, the former admit a kind of universality in which the coarse-grained dynamics is insensitive to the form of the interaction below this scale, while the latter do not. This is a basic qualitative distinction between the dynamics in these two cases, which corresponds naturally to what one call "short-range" or "long-range". In order to distinguish it from the canonical distinction based on thermodynamical considerations, following [72], we can refer to > d as dynamically short-range, and < d as dynamically long-range 2For what concerns quasi-stationary states the implications of this result and relevance of the classification are simple: for all cases one would expect that such states may exist (since the Vlasov limit exists) but the conditions for their existence, which requires that the time scales of their persistence be long compared to the system dynamical time, will be very different. For > d, their lifetime, which would be expected to be directly related to the amplitude of the non-Vlasov term, can be long only if the smoothing scale in the force is sufficiently large; in the case < d, their lifetime will be expected to be independent of ". We note that this is precisely in line with results of analytical calculations based on the Chandrasekhar approach to estimation of the relaxation rate, and the results of numerical simulations of systems of this kind reported in [50].

These results are all, as we have emphasized, built on our central hypothesis that correlation of density fluctuations, associated with the finite particle number, may be neglected. To formulate it we must define a smooth mean phase space density by introducing a coarse-graining scale, which is assumed to be a "mesoscopic scale": arbitrarily small compared to the system size, and yet large enough so that the phase space cell contains many particles. Our central hypothesis is not one of which we can prove the validity, but it is a consistent, simple and physically reasonable, analogous to that of "molecular chaos" in the derivation of the Boltzmann equation. It is also in line with the fact that when Vlasov approximation is valid, the system dynamics is defined in terms only of the mean density and fluctuations due to higher order density correlations are considered subdominant. We note above all that it leads to conclusions and behaviours which are very reasonable physically, even in the cases we have not focussed on but to which our analysis can be applied. For example, if we consider a hard repulsive core interaction without a smoothing, i.e. let " ! 0 for large, we infer that the force fluctuation on a coarse-grained cell diverges and that there is thus no Vlasov limit. Indeed in this case an appropriate two body collision operator is required to take into account the effects of interactions between particles.

Finally we note that the ensemble we have assumed to describe the fine-grained phase space density is a realisation of a Poisson process with mean density given by the coarse-grained space. This implies that we include (Poissonian) fluctuations of particle number at all scales: indeed these (finite particle number) fluctuations are the source of the fluctuation terms in the dynamical equations which we have analysed. One could consider that it might be more appropriate physically to take, in estimating the fluctuations induced at a given coarse-grained density, the ensemble in which the particle number is constrained to be fixed in all coarse-graining cells. In other words we could average, given a coarse-grained phase space density, only over the configurations in which the particles, of fixed number in each cell, are distributed randomly within the cells. It is straightforward to verify, with the appropriate modification of the average, that doing so can only change our results for what concerns the large scale contributions to the force fluctuations: the diverging behaviours at small scales we have focussed on arise from the fact that there is a finite probability for a particle to be arbitrarily close to another one, and the local value of the density will at most modify the amplitude but not the parametric dependence of this term. On the other hand, our determinations of the parametric dependence of the contributions to the force fluctuations from the bulk will be expected to depend on how the particle fluctuations at larger scales are constrained 3 .

In future work we plan to explore the possibility of developing the approach used in this article to understand and describe further the effect of the "non-Vlasov"collisional -terms on the evolution of a finite N system, i.e., to use it to develop a kinetic equation for the N particle system, alternative to those developed in the literature (see, e.g., [66,69] for a discussion and references). We would to clarify in particular the relation between this approach and that of Chandrasekhar which, as we have noted, when extended to a generic interaction has been shown [50] to give very consistent conclusions about the sensitivity of collisional relaxation of quasistationary states to the small scale regularisation. Likewise it would be interesting to try to test directly in numerical simulations for the validity of our central hypothesis about correlations, and also characterize using analytical methods the robustness of our central conclusions to the existence of different weak correlations of the density fluctuations.

7 Appendix 1: Alternative derivation of h⇠ 2 F i

We present in this appendix an alternative (but equivalent) calculation of the mean and variance of ⇠ F . We recall it's definition:

⇠ F (x, v, t) = F(x, v, t) F 0 (x, t), (2.65)
where we have defined:

F(x, v, t) = X i2C(x,v) Z ⌦ g(x i x 0 )n k (x 0 )d d x 0 . (2.66)
To evaluate the moments, it is useful first to introduce n C(x,v) the microscopic density of particle restricted to the cell C(x, v):

n C(x,v) (y 0 ) = Nc X i2C(x,v) (y 0 y i ), (2.67)
where y i are the position of the particle inside the cell. We introduce also the restricted Klimontovich density n k/C which does not include particle inside the cell such that:

n k (x 0 ) = n C(x,v) (x 0 ) + n k/C (x 0 ). (2.68)
Then we express F as:

F(x, v) = 1 N c Z ⌦ d d x 0 Z ⌦ d d y 0 n k/C (x 0 )n C(x,v) (y 0 )g(x 0 y 0 ), (2.69)
where the integral over n C (x 0 )n C (y 0 ) vanishes because of the symmetry of the two body force g( x) = g(x).

Inside the cell we distribute the particle with an homogeneous Poisson process. Out side the cell, by integration of the phase space probability (2.20) over the velocities v 1 , ...v N , this is given by

P(x 1 , x 2 .., x N ) = N Y i=1 p(x i ) (2.70)
where

p(x i ) = ⇢ d x , x i 2 C(x, v) n 0 (x i )
N Nc , otherwise.

(2.71)

In other words, we assume the N c particles in the cell to be distributed randomly in the cell, and the remaining N N c particles to be distributed according to the inhomogeneous Poisson process with mean density n 0 (x). Using the latter, the follows that

⌦ n k/C (x) ↵ = n 0 (x).
(2.72)

Further we compute with, Eq. (2.25) the mean of n C(x,v)

⌦ n C(x,v) (y) ↵ = N c (x, v) x .
(2.73)

Moreover to compute the variance of F, it is useful to compute the correlation functions. Outside the cell we have, by definition of an inhomogeneous Poisson process:

⌦ n k/C (x 0 )n k/C (x 00 ) ↵ = n 0 (x 0 ) (x 0 x 00 ) + (1 1 N N c )n 0 (x 0 )n 0 (x 00
).

(2.74)

Inside the cell the correlation are zero as well, thus for any y 0 , y 00 2 C(x, v) we have:

hn C (y 0 )n C (y 00 )i = N c d x (y 0 y 00 ) + N c (N c 1) 2d x , (2.75)
where the (x, v) dependence of N c is implicit. Finally, we have assume no correlation between particles inside and outside the cell, and then:

⌦ n k/C (x 0 )n C(x,v) (y 0 ) ↵ = ⌦ n k/C (x 0 ) ↵ ⌦ n C(x,v) (y 0 ) ↵ . (2.76)
Mean of ⇠ F Using the expression Eq. 2.69 and the mean of the densities computed above we obtain:

h F(x)i = Z S(x, x) d d y d x Z ⌦ d d x 0 n 0 (x 0 )g(y 0 x 0 ) = Z S(x, x) d d y d x F 0 (y 0 ), (2.77)
where S = S(x, x ) stands for an integration inside the ball centered in x of radius

x in real space, as represented in Fig. 2.1. Then the mean of the force fluctuation

⇠ F gives h⇠ F (x, v, t)i = Z S(x, x) d d y d x F 0 (y 0 ) F 0 (x, t). (2.78)
The first term represents the mean value of the mean field force F 0 over the coarsegaining cell. Though, the mean field is defined at a coarse-grained level, then under the scale x , we do not make any other assumption in neglecting its fluctuation. Thus We obtain :

h⇠ F (x, v, t)i = 0. (2.79)
Variance of ⇠

F

The mean field being constant over realisations we have:

h⇠ 2 F (x, v, t)i = h F2 (x, t)i F 2 0 (x, t). (2.80)
Using the expression (2.69), we express the variance of F as:

⌦ F2 (x, v) ↵ = 1 N 2 c (x, v) Z Z Z Z ⌦ n k/C (x 0 )n k/C (x 00 ) ↵ hn C (y 0 )n C (y 00 )i g(y 0 x 0 )g(y 00 x 00 )d d x 0 d d x 00 d d y 0 d d y 00 . (2.81)
Using this last equation as well as the correlation function of the density eq. (2.75) and eq. (2.74), we get four terms. Two of these terms can be written as a function of the mean-field force squared and we obtain the following expression for the variance of the force noise:

h⇠ 2 F (x, v)i = 1 N c (x, v) Z ⌦ d d x 0 n 0 (x 0 ) Z S(x, ) d d y d g 2 (y 0 x 0 ) + (1 1 N c (x, v) ) Z ⌦ d d x 0 n 0 (x 0 ) Z S(x, ) d d y 0 d g(y 0 x 0 ) 2 1 N N c F 2 0 (x).
(2.82)

Approximating N N c ⇠ N , we obtain the same result than Eq. (2.37).

Chapter 3

Long-range systems with weak dissipation 1 Introduction

We previously discussed, as a function of the range of the pair force, the validity of the Vlasov equation to describe the dynamics of a N body system, and the possible presence of QSS. These results apply a priori for a strictly Hamiltonian system and one may wonder what happens to such states when the system is subjected to non-Hamiltonian perturbations. In this chapter, we address the general issue of the robustness of QSS in presence of dissipation. We consider a long-range interacting system with two different classes of dissipative forces: on the one hand, a viscous damping force, which we will refer to as the viscous damping model (VDM); on the other hand, instantaneous inelastic, but momentum conserving, collisions, which we will refer to as the inelastic collisional model (ICM).

For the two classes of model we derive a kinetic equation in the weak dissipative limit taking these effects into account. The precise meaning of weak here and the derivation are detailed in section 2. In the next section, the kinetic approach for the two class of models shows that a long-range interacting system with dissipation may admit scaling solutions. These solutions describe the evolution of the system through a sequence of QSS, remaining virialised at all time. In spite of the shrinking of the system due to the dissipation, when the system follows such scaling quasi-stationary states the phase space distribution remains fixed up to a rescaling.

To check these theoretical predictions, numerical simulations of a 1d self-gravitating system have been performed. The latter is simulated first in the presence of a fluid damping force, and is thus an example of the VDM class of systems. Then, a second numerical investigation of the 1d self-gravitating system is considered with particles undergoing inelastic collisions in order to also test the predictions with an example of ICM. The simulation method used is detailed in section 4.

In section 5, the numerical results are presented. In order to illustrate and to complete the description of the dynamics of long-range interacting system systems given in the chapter 1, we first reproduce some known results for a 1d self-gravitating system without dissipation. Then the simulations with the two types of dissipation are presented and discussed in the context of the kinetic theory developed.

Mean field limit for long-range systems with dissipation

In this section, for the two class of models studied, we derive, via a BBGKY hierarchy generalised to non-Hamiltonian systems, a kinetic equation describing the evolution of the one particle phase space distribution, in the mean field limit. It takes the form of a Vlasov equation with an additional term describing the non-conservative effects of the given model. For the VDM, this first derivation is valid in any dimension.

For the ICM, for simplicity, the derivation of kinetic equation is presented in 1d. Subsequently, we will restrict ourselves, for both models, to the 1d case. We consider the case of N particles of identical mass m interacting through the long-range power law potential

(|r|) = g n|r| n , n < d, (3.1)
where g > 0 is the coupling constant. We consider attractive forces for simplicity.

For n < d, we will consider below the large N mean field limit and thus implicitly we assume as seen in the previous chapter that the singularity in the potential at r = 0 is regulated. However for n < d 1, even the leading correction to the Vlasov equation are expected to be insensitive to the value of such a cut-off (if it is "soft").

Dissipation through viscous damping

In addition to the long-range potential, we consider the general viscous damping force:

F vd = ⌘m|v| ↵ v , (3.2)
where ↵ and ⌘ are constants and v = v |v| is the unit vector parallel to the particle velocity. Without dissipative effects, the one particle distribution function in phase space f is governed, in the mean field limit, by the Vlasov equation. We saw, in chapter 1, that the Vlasov equation may be derived starting from the N particle probability distribution function f N (x N , v N , t) through the BBGKY hierarchy. We come back to this point in order to express the term corresponding to a general viscous damping force in the kinetic equation. To do so we generalise the BBGKY hierarchy derivation. Here we start also from the continuity equation in the N particle phase space:

@f N (x N , v N , t) @t + N X i=1  @x i @t • @f N (x N , v N , t) @x i + @v i @t • @f N (x N , v N , t) @v i = 0. (3.3)
For the VDM, the non-Hamiltonian equations of motion, are given by:

ma i = F i + F vd = m X i6 =j r i (x i x j ) ⌘mv ↵ i vi . (3.4)
Using Newton's second law in the continuity equation Eq. (3.3) we obtain the following generalised Liouville equation:

@f N @t + N X i=1 ( v i • @f N @x i X j6 =i r i (x i x j ) • @f N @v i ) = ⌘ N X i=1 @ @v i |v i | ↵ vi f N . (3.5)
Integrating over N 1 coordinates, and assuming f N (x N , v N , t) = P N i=1 f 1 (x i , v i , t) in the mean field limit, the left-hand side of the equation gives the three usual terms of Vlasov equation. Performing this integration on the right-hand side of Eq. (3.5), we have

⌘ N X i=1 Z dx 1 ..dx N 1 Z dv 1 ..dv N 1 @ @v i |v i | ↵ vi f N = ⌘ @ @v |v| ↵ vf 1 (x 1 , v 1 , t). (3.6)
For a long-range system with viscous damping, in the mean field limit, the evolution of the one particle distribution is therefore described by a modified Vlasov equation with an operator J V DM [f ] accounting for the viscous damping:

@f @t + v. @f @x + F[f ](x). @f @v = J V DM [f ] =⌘ @ @v .|v| ↵ vf (x, v, t), (3.7)
where

F[f ](x) = r (x) (3.8)
is the mean field force with

(x) = 1 M Z (|x x 0 |)f (x 0 , v 0 , t)dx 0 dv, (3.9)
the mean field potential.

Cosmological motivations

We mention first one motivation for this model is given in the context of non linear structure formation in cosmology. In this context, we can often make the approximation that particles, whether stars, galaxies or "dark matter" particles, interact only through Newtonian gravity. Indeed, we suppose that other fundamental forces are not relevant at the scale of a star cluster, a galaxy or galaxy cluster. The equations of motion for a finite system of N self-gravitating particles in an expanding universe, with identical masses m are:

d 2 x i dt 2 + 2H dx i dt = Gm a 3 N X i6 =j sgn(x i x j ) |x i x j | d 2 , (3.10)
where x i are the positions of the particles in "co-moving" coordinates, t the "cosmic time". a is the scale factor, G the gravitational constant and H = ȧ/a the Hubble constant.

Changing the time coordinate dt = a 3/2 d⌧ , equation (3.10) becomes

d 2 x i d⌧ 2 = (⌧ ) dx i d⌧ + F i with F i = Gm N X i6 =j sgn(x i x j ) |x i x j | (n+1) , (3.11)
where (⌧ ) depends on the scale factor a. Thus in these coordinates the expansion of the universe is expressed by a term equivalent to a fluid damping force in an euclidean space. Further for the reference Einstein De Sitter model with a / t 2/3 one has = C st .

Dissipation by inelastic collisions.

We consider now a 1d long-range system undergoing momentum conserving, dissipative collisions. We describe first the microscopic rules by which the velocities change for two particle collisions. To describe the effects of many dissipative collisions on the dynamics, we introduce the pseudo-Liouville operator, used in the literature on granular systems, which models the irreversible evolution of the system. Then we derive a kinetic equation for the evolution of the one particle density, in both the mean field (N ! +1) and the quasi-elastic limit (c R ! 1). We show that these limits can be taken consistently.

Collision rule

We consider an inelastic collision of a pair of particles 1, 2 of the same mass m. The particle are ordered with positions x 1 < x 2 . The coefficient of restitution c R , the ratio of relative velocity before and after the collision, is canonically defined by the relation:

(v ⇤ 2 v ⇤ 1 ) = c R (v 2 v 1 ), (3.12)
where v ⇤ denotes the post-collisional velocities and where c R is constant. Using the conservation of the total momentum, one deduces the following 1d collision rule:

v ⇤ 1 =v 2 + 1 c R 2 (v 1 v 2 ) (3.13) v ⇤ 2 =v 1 1 c R 2 (v 1 v 2 ).
For c R = 1, the Eq. (3.13) results in simple exchange of the particle's velocities: the collision is elastic; while, for an inelastic collisions 0  c R < 1. The last relation may be inverted and written in terms of pre-collisional velocities v ⇤⇤ :

v ⇤⇤ 1 =v 2 + 1 c 1 R 2 (v 1 v 2 ) (3.14) v ⇤⇤ 2 =v 1 1 c 1 R 2 (v 1 v 2
). An inelastic collision, thus results in a loss of energy per collision:

K 12 = 1 2 m(v ⇤2 1 + v ⇤2 2 v 2 1 v 2 2 ) = m 1 c 2 R 4 (v 1 v 2 ) 2 . (3.15)

Pseudo-Liouville equation

To derive a kinetic term for non-Hamiltonian hard-core collisions we use the collisional operator introduced by Ernst in [75] for elastic collisions and extended to treat inelastic collision in the literature on granular systems [76][77][78]. Because of the dissipative nature of the collisions, the dynamics is irreversible, and the operator is not self-adjoint. One may then construct two binary collision operators, one for forward time evolution T + and its adjoint for backward evolution: T + = T 6 = T + . We construct T + as follows. First, we define the two operators that change velocities according to the collision rules (3.13). Let O(v i , v j ) be an observable depending on the two velocities v i , v j involved in the collision. We define the operators b ⇤ and b ⇤⇤ through the following relations:

b ⇤ O(v i , v j ) = O(v ⇤ i , v ⇤ j ) , (3.16)
i.e. b ⇤ changes the pre-collisional velocities into post-collisional velocities, and

b ⇤⇤ O(v i , v j ) = O(v ⇤⇤ i , v ⇤⇤ j ) , (3.17)
which performs the reverse operation and accounts for the backward evolution. Second, we consider here point particles (for 1d system, the particle size is irrelevant and can be set to 0), with hard core potential. The latter acts only when particles are in contact and thus the spatial dependence of the operator is T + / (x i x j ). Moreover, we consider the change of the observable during an interval [t, t + t] including the collision event. Then T + / b ⇤ Id where Id is the operator identity. Finally, the binary collision operator is defined as:

T (ij) + = |v i v j | (x i x j )(b ⇤ Id), (3.18)
where the factor |v i v j | is a flux term accounting for collisions to occurred only for approaching particles [77]. The time evolution operator, which is the adjoint of T (i,j) +

:

T (ij) + = T (ij) = |v i v j | (x i x j )( 1 c 2 R b ⇤⇤ Id). (3.19)
The probability that three particles collide at the same instant is assumed to be zero. Then, one can write the pseudo-Liouville equation, with a right hand side accounting for all possible binary collisions:

( @ @t + L (N ) 0 )f (N ) (x N , v N , t) = X i<j T (i,j) + f (N ) (x N , v N , t), (3.20)
where

L (N ) 0 = P N i L i 0 = P N i v i . @ @r i + a i . @ @v i .

Kinetic equation

We now integrate Eq. (3.20) over (N 1) coordinates to obtain the first equation of the BBGKY hierachy. Then assuming molecular chaos (Stosszahl ansatz)

f 2 (x 1 , x 2 , v 1 , v 2 , t) = f (x 1 , v 1 , t)f (x 2 , v 2 , t)
, and integrating the right hand side of Eq. (3.20), we obtain the operator J ICM [f ] that accounts for the dissipative collisions. With the same assumption, the left hand side gives the three terms of the Vlasov equation. Then, in the mean field limit, the kinetic equation for ICM is given by: @f @t + v

@f @x + F [f ](x, t) @f @v = J ICM [f ] , (3.21)
where F [f ](x, t) is the mean field force, as defined previously, and,

J ICM [f ] = N M Z dv 1 dx 1 T (1,2) + f (x 1 , v 1 , t)f (x 2 , v 2 , t) , (3.22)
or in its explicit form:

J ICM [f ](x, v, t) = N M Z dv 1 |v v 1 |  f (x, v ⇤⇤ , t)f (x, v ⇤⇤ 1 , t) c 2 R f (x, v, t)f (x, v 1 , t) , (3.23)
where the factor 1/M comes from the normalisation of f here taken to be

Z f (x, v, t)dxdv = M . (3.24)
We note that, for both sides of the equation, the derivation requires the molecular chaos assumption: the correlation between particles is neglected. This approximation is valid at low density, where all collisions are supposed to be independent events. In other words, before each collision, the particles has no memory of the past collisions.

Weak dissipative limit.

So far, it has not been established that the collisional operator goes to a finite value when N ! +1. More specifically we need to determine a scaling of this term so that the energy dissipation per unit time is independent of N .

The mean field limit taken to obtain the Vlasov equation is defined by taking N ! +1 with m / 1/N and g / 1/N (so that the mean number density n 0 ⇠ N but the mass density ⇢ 0 = C st ). This limit fixed the mean field time scale:

⌧ mf = 1 p gn 0 .
(3.25)

It is the characteristic time scale for the mean field dynamic, and in particular for the establishment of a QSS. During this time scale, a particle typically "crosses" the system1 and the number of collision undergoes by the system is then / N 2 . Then using Eq. (3.15), during ⌧ mf the dissipation of energy E / 1 c 2 R 4 N 2 . In order to take the limit N ! +1 and c R ! 1 such that this amount remains fixed, we define

= (1 c R (N ))N 2 (3.26)
to be the order parameter for this kind of dissipation, and take both limit keeping this parameter constant. This limit corresponds, to the quasi-elastic ((1 c R ) ⇠ 1 N ! 0) limit usually taken in the literature of granular gas [79,80]. We will refer to it here as the weak dissipation limit.

In this limit where q = 1 c R 2 ⇠ 1 N ⌧ 1 we can perform a Kramers-Moyal expansion of the collisional term [80]. To do so we split the integral in Eq. (3.23) into a gain term G q depending on q and a loss term L.

J ICM = G q (x, v, t) L(x, v, t).

(3.27)

The loss term is the rate at which particles leave the phase space element dxdv around (x, v):

L(x, v, t) = N M Z dv 0 |v 0 v|f (x, v 0 , t)f (x, v, t) , (3.28)
and the gain term G q , the rate at which dvdx is populated by the collisions, may be expressed with a Dirac distribution imposing the constraint on the two pre-collisional velocities v 0 and v 00 to change according to the rules (3.14) such that one of the two post-collisional velocities is v:

G = N M Z Z dv 0 dv 00 |v 0 v 00 |f (x, v 0 , t)f (x, v 00 , t) (v qv 0 (1 q)v 00 ). (3.29)
Then, taking q ⌧ 1 we expand the Dirac distribution to second order in q:

(v qv 0 (1 q)v 00 ) = (v v 00 + q(v 00 v 0 )) (3.30) = (v v 00 ) + q(v 00 v 0 ) 0 (v 00 v 0 ) + 1 2 q 2 (v 00 v 0 ) 00 (v 00 v 0 ),
where 0 , 00 are the successive the distributional derivatives of the Dirac distribution, defined by their action on the test function '.

Z dx 0 (x)'(x) = Z dx (x)' 0 (x). (3.31)
Using the expanded form of and the last equation in the Eq. (3.29), at the same order of approximation q ⇠ 1 N than for the Vlasov equation, we obtain for G q :

G q = N M  f (x, v, t) Z |v 0 v|f (v 0 )dv 0 + 1 c R 2 @ @ v Z |v 0 v|(v 0 v)f (v)f (v 0 )dv 0 .
(3.32) Gathering up Eqs. (3.28) and (3.32) we obtain, in the weak dissipative limit, the approximate form of the collision operator which we will use in what follows:

J ICM [f ] = M @ @v Z a 1 (x, v, t)f (x, v, t) , (3.33)
where

a 1 = R dv 0 (v 0 v)|v 0 v|f (x, v 0 , t).

Scaling quasi-stationary states

For both classes of dissipation considered, we have now derived, from continuity equations in the N body phase space and shown that their evolution, is, in a well defined mean field limit, an additional term, denoted J[f ], on the right-hand side of the Vlasov equation:

@f (x, v, t) @t + v. @f (x, v, t) @x + F [f ](x). @f (x, v, t) @v = J[f ] , (3.34)
where

F [f ](x) = g M R dx 0 dv 0 sgn(x 0 x)|x 0 x| (n+1)
. For the case of the VDM, in 1d, this operator is:

J V DM [f ] = ⌘ @ @v (v ↵ f (x, v, t)) , (3.35)
while for the ICM, it is expressed, in the quasi-elastic limit, as:

J ICM [f ] = Nm @ @v  f (x, v, t) Z dv 0 (v 0 v)|v 0 v|f (x 0 , v 0 , t) . (3.36)
We now look for scaling solutions to this equation using the following ansatz:

f scal (x, v, t) = M x(t)v(t) f ✓ x x(t) , v v(t) ◆ , (3.37)
where f is a dimensionless function. In what follows, all dimensionless quantities will be written with a tilde over their symbols. Inserting the ansatz in Eq. (3.34), obtain:

M x(t)v(t) " 1 x(t) dx(t) dt f (y, z) + 1 v(t) dv(t) dt f (y, z) 1 x(t) dx(t) dt y @ f (y, z) @y (3.38) 1 v(t) dv(t) dt z @ f (y, z) @z # + 1 x2 (t) z @ f (y, z) @z + gM F [ f ](y) v2 (t)x (n+2) (t) @ f (y, z) @z = J[f scal ] ,
where y and z rescaled variables (y = x x(t) , z = v v(t) ), and

F [ f ](y) = Z sgn(y y 0 )|y y 0 | (n+1) f (y 0 , z)dy 0 dz , (3.39)
is the dimensionless mean-field force. Further, we have for the VDM:

J V DM [f scal ] = M ⌘v ↵ 2 (t) x(t) @(z ↵ f (y, z)) @z , (3.40)
and for the ICM:

J ICM [f scal ] = M x2 (t) @( ã1 (y, z) f (y, z)) @z , (3.41)
where

ã1 (y, z) = Z (z z 0 )|z z 0 | f (y, z 0 )dz 0 . (3.42)
The existence of a scaling solution imposes to choose x and v such that each term of Eq. (3.38) has the same time dependence, implying the following set of proportionality relations:

dx(t) dt 1 x(t) / dv(t) dt 1 v(t) / 1 x2 (t) / 1 v2 (t)x n+2 (t) / ( ⌘v ↵ 1 (t), for VDM v(t)
x(t) , for ICM . This means that the viral ratio We now assume further that the system is in a QSS in the limit that the dissipation is absent. With this assumption, and with Eq. (3.44), we have that

R = 2K nU / v2 (t)x n (t) , (3.45
z @ f (y, z) @y + gM v2 xn Ā(y) @ f (y, z) @z = 0 . (3.46)
i.e., the last two terms on the right hand side of Eq. (3.38) cancel. This also implies that the virial ratio R is unity. Physically this means we assume that all time dependence of the evolution arises solely from the dissipation. This corresponds to an adiabatic limit of weak dissipation in which the time scale on which the dissipation causes macroscopic evolution is arbitrarily long compared to the time scale associated with the mean-field dynamics. In other words, the system has always the time to virialize i.e. to redistribute a part of the loss of kinetic energy in a proportional loss of potential energy. The scaling solution thus excludes all non-trivial time dependence due to the mean-field dynamics beyond its effect in virializing the system, and in particular does not describe the phase of violent relaxation to virial equilibrium. We will evaluate further below the validity of this crucial approximation. Using Eq. (3.46) the constraints on x and v reduce now to the simple one:

dv(t) dt 1 v(t) / ( ⌘v ↵ 1 , for VDM v x ,
for ICM Now, using this two last relation of proportionality it follows that the only requirement on the scaling solution is

dv(t) dt 1 v(t) = A 0 ⌘v (t), (3.49)
with

= ↵ 1, (3.50)
for the VDM, and

dv(t) dt 1 v(t) = A 1 v(t)x(t) 1 = A 1 v , (3.51)
with

= (n + 2)/n, (3.52)
for the ICM, where A 0 and A 1 are dimensionless positive constants. Integrating these equations, we obtain

v(t) = v0 8 < : ⇣ 1 + sgn( ) t tc ⌘ 1 6 = 0 e t tc = 0 , (3.53)
where t c is a characteristic time scale and v0 is an arbitrary constant. The solutions for x(t) follow from Eq. (3.44).

x(t) = x0 Further, one can deduce the energy decay law of the system. Using the virial theorem 2K + nU = 0, the total energy E may be expressed as:

E = (1 2 n )K. (3.55)
Finally the energy evolves as

E(t) = E 0 8 < : ⇣ 1 + sgn( ) t tc ⌘ 2 6 = 0 e 2t tc = 0. (3.56)
For attractive pair potentials with n < 0, which is the class of long-range potentials we are considering here (in d = 1), the scaling solution therefore describes, when < 0, a system which undergoes a collapse in the finite time t c . This feature of the dynamics is similar to the one of a dissipative granular system with a zero pair potential between collisions. Let us recall that, for such a system, the energy decay law is given by the so-called Haff's law [81] which describe also a collapse in a finite time.

For = 0, the system undergoes a monotonic contraction characterized by the time t c , but never collapses. For the ICM, this is the case for a system interacting by a pair potential with exponent n = 2 (a harmonic potential). For the VDM, this corresponds to a linear fluid damping for which ↵ = 1.

We now return to the essential approximation we have made in deriving the scaling solution. This corresponds to assuming that ⌧ diss ⌧ mf , where ⌧ diss and ⌧ mf are the characteristic times for, respectively, the dissipation of the system energy E and the mean-field (Vlasov) dynamics (where ⌧ mf has been defined above in Eq. (3.25)). We now work out how ⌧ diss evolves with time in the scaling solution.

To do so, coming back to Eq. (3.34), we derive an expression for t c starting from the arbitrary time t = 0. With the assumption, on the time scale ⌧ diss , one can ignore the two last terms of the right hand side of modified Vlasov equation Eq. (3.34). This means that we neglect the rapid exchange between the potential and the kinetic energy of the long-range system, and so K / U / E / C st . Then multiplying the right-hand side of the Vlasov equation by v 2 and integrating over position and velocity coordinates, the only remaining term gives:

Z Z @f (x, v, t) @t v 2 dxdv = 2 dK(t) dt . (3.57)
with K the kinetic energy of the system. Further, multiplying by v 2 and integrating over position and velocity the left hand side for the two models, we obtain, for the VDM,

dK dt = ⌘ 2 Z |v| ↵+1 f (x, v, t)dvdx; (3.58)
and for the ICM, remarking the non symmetric part of the function

(v, v 0 ) = v(v v 0
) under the exchange v $ v 0 vanishes, we obtain:

dK dt = 4 Z Z |v v 0 | 3 f (x, v, t)f (x, v 0 , t)dvdv 0 dx. (3.59)
Substituting the scaling solution, in the two last equations, we obtain, for the VDM,

I K dv 2 (t) dt = ⌘|v(t)| +2 I V DM ; (3.60)
and for the ICM, with the use of relation Eq. (3.44), we get

I K dv 2 (t) dt = 1 2 |v(t)| 3
x(t)

I ICM = 1 2 |v(t)| +2 I ICM , (3.61)
where

I K = R z 2 fdydz, I V DM = R z ↵+1
fdydz, and I ICM = R R |z z 0 | f fdydzdz 0 are dimensionless finite integrals depending on f . Finally, remarking that

dv 2 dt t=0 / 1 t c v2 (0), (3.62)
is it easy to work out that, for the two models we respectively obtain the scaling of t c :

t c / 1 ⌘ v2 (0) |v(0)| +2 , (3.63)
where ⌘ is replaced by for the ICM. For both cases, t c diverges as the inverse of the strength of the dissipation as one would expect. t c represents the time scale for dissipation starting from the (arbitrary) time t = 0. In the scaling solution, the characteristic time for the dissipation starting from an arbitrary time t thus scales as

⌧ diss (t) / v (t) (3.64)
For a typical system size x(t), the mean field acceleration scales as x (n+1) , and ⌧ mf , the typical time for a particle to cross the system scales as;

⌧ mf / x n+2 2 . (3.65)
With the use of Eq. (3.44), it follows that

⌧ diss /⌧ mf / v[ + n+2 n ] (3.66)
Therefore if > c = n+2 n the ratio of these timescales increases as a function of time. In other words, if > c , the scaling solution drives the system to a regime in which the approximation underlying it becomes arbitrarily well satisfied: the separation of time scale ⌧ diss ⌧ mf becomes more and more accurate. In this case, then, we might expect that the scaling solution may be an attractor of the dynamics, while for < c , the opposite is the case and the scaling solution might describe a transient behavior. The case = c , which corresponds precisely to the ICM, is the marginal one. In this case, the ratio ⌧ diss /⌧ mf remains constant in the scaling solution, and one might thus expect this solution to be a transient which persists on a time scale dependent on this ratio.

Numerical study: 1d self-gravitating system

To illustrate and check the theoretical results, and in particular the relevance of the scaling QSS solutions we have derived from the different kinetic equations, we performed numerical simulations on a reference long-range interacting system toy model: the 1d self-gravitating system. In this model, as discussed in chapter 1, the force is constant between particle collisions and the N equations of motion can be integrated exactly. Then, to simulate the system an event driven algorithm can be implemented. No discretization of the equations of motion is necessary and the dynamics is exact, in the sense that the code is only limited by the round-off of the computer [82]. Moreover, this 1d model exhibits many dynamical features similar to 3d self-gravitating systems.

In this section, we present the numerical method used to implement the simulation of three 1d self-gravitating system:

• the model without dissipation, which we refer to as SGS. We use this model as a test of our numerical simulation by comparing with known results.

• the SGS with an additional linear fluid damping force which is an example of the VDM, corresponding to the case ↵ = 1 and > c .

• the SGS with inelastic collisions. This is an example of the ICM with = c = 1 (n = 1).

All three models are simulated numerically with an event driven method. In the algorithm, the time evolution consists in elementary steps of the dynamics : The algorithm can be divided into steps:

1. computation of the next collision times of each pair of particles;

2. determination of the pair with the shortest time to its next collision;

3. update of the colliding particles positions and velocities forward up to the next collision time; 4. change of particle's velocities following collision rules.

Steps 1 and 3 require the equations of motion to be integrated. For step 2, we use a basic sorting algorithm. The object "particle" has several attributes: the position, the velocity, the next collision time t col and the time up to which its position is updated. Indeed, in order to reduce the time of the simulation we do not move all the particles at each step, but only the two particles involved in the collision, and then also the two next neighbouring particles in order to compute the three new collision times [83]. Thus, at a given step, all particles are not at the same physical time. Then, in order to extract observables at a given physical time t, we update all the particles up to this time. At step 4, the appropriate collision rule is implemented (elastic for SGS and VDM, inelastic for ICM). To simulate all three models the code used is based on an existing code (adapted where appropriate) written by F. Sicard to simulate in particular cosmological toy models of this kind. This code uses periodic boundary conditions rather than open boundary condition appropriate for our model. However we will show below that in the limit that the ratio of the size of the system compared to the size of the periodic box goes to zero, the force in the periodic system converges to the expression of the force for a isolated system. If this ratio is sufficiently small, the force is an arbitrary good approximation to the one of the our model. Finally, in the last subsection we present the different distributions we used as initial conditions for our simulation, and explain how we generate them. The following section will be devoted to our numerical results.

SGS-1d self-gravity (without dissipation)

Equation of motion of an isolated SGS

We recall the equation of motion of the 1d SGS model (discussed also in chapter 1). We consider a finite system of N identical self-gravitating particles with masses set equal to 1 for convenience. The force acting on particle i is given by:

F N (x i ) = g N X j=1 sgn(x j x i ) = g[N > (x i ) N < (x i )], (3.67)
where N > (x i ) (and N < (x i ) ) is the number of particles on the right (the left) of x i . Between two collisions, the force on each particle is thus constant. The equations of motion of the system are therefore:

d 2 x i dt 2 = g[N > (x i ) N < (x i )]. (3.68)
Then between two collisions the motion can be integrated trivially and an event driven dynamics can be implemented. This simulation method has been used by other studies (e.g. [START_REF] Worrakitpoonpon | Relaxation of isolated self-gravitating systems in one and three dimensions[END_REF]82]). As mentioned and detailed further below, we used a slightly different simulation method with periodic boundary conditions.

Time unit of the simulation

In order to estimate the characteristic time of the dynamics and choose an appropriate time unit in the simulation, one considers the following situation. At t = 0, we imagine the N particles distributed on a regular lattice between x = 0 and x = L 0 with lattice spacing L 0 N = n 1 0 . The equations of motion are:

ẍj = g N X i6 =j sgn(x i x j ) = g(2j N ). (3.69)
Integrating these equations from t = 0 (x i (t = 0) = ( 1 2 + i)n 1 0 ) and solving for the time ⌧ cm at which all particles reach the center of mass L 0 /2, we obtain ⌧ cm = 1 p gn 0 . For convenience, following the choice of [18], we define the unit of time for our the simulations as twice this time:

⌧ dyn = 2 q g N L 0 . (3.70)
We shall see in section 5 that ⌧ dyn is the characteristic period of oscillation of the virial ratio: it is roughly the time scale in which the system converts potential energy into kinetic energy and inversely. It is also called the mean field time and the notation ⌧ mf may also be used.

Equation of motions for a periodic system

To simulate the dynamics, as noted we used a code designed to simulate a 1d SGS which is infinite and periodic. We used this code because it allows one to study in an analogous manner to cosmological simulation of the universe (in 3d), initial conditions in which matter is distributed almost uniformly throughout space (using e.g. small displacement from a regular lattice [84,[START_REF] Sicard | Out-of-equilibrium dynamics in infinite one-dimensional selfgravitating systems[END_REF]). In practice we do not report results for this case here. The definition of the force Eq. (3.67) is for a system containing a finite number of particles, but it is ill defined for an infinite periodic system. The code used is based on the following regularisation of the expression (3.67):

F 1 (x i ) = lim µ!0 X j6 =i e µ|x i x j | sgn(x i x j ). (3.71)
For an infinite periodic distribution with N points in each cell, it can be shown [84,[START_REF] Miller | [END_REF] that this expression can be written as:

F 1 (u i ) = 2gn L (u i hui) , (3.72)
where n L . = N L is the number density in the box, and u i are the positions of the particles ordered from i = 1, .., N in [0, L] relative to a regular lattice i.e. u i = x i s i where s i = (i 1 2 )n 1 L ; and hui = 1 N P i u i is the average particle displacement. In practice, and in what follows, we take hui = 0 or equivalently we set the center of mass of the system to be at the center of the box. Substituting u i = x i s i in the expression of the regularised force given by Eq. (3.72) and then noting that N < (x i ) = i 1 and N > (x i ) = N i, the force may be expressed as [87]:

F 1 (x i ) = g[N > (x i ) N < (x i )] 2gn L (x i hxi)), (3.73)
where hxi is the center of mass of the system. The first term of Eq. (3.73) is ⇠ gN while the second is ⇠ gn L L 0 = g↵N with L 0 = ↵L the size of the system and ↵ the ratio between the system size and the size of the prediodic box. Thus, if ↵ ⌧ 1, the second term ("background" term2 ) is negligible and converges to the expression (3.67).

Integration of the equations of motion.

From Eq. (3.72) we see that each particle is simply evolving is an inverted harmonic potential centred on its lattice site. The equations of motion, valid between two collisions, becomes:

d 2 u i dt 2 1 ⌧ 2 L u i = 0, (3.74)
with

⌧ L = 1 p 2gn L . (3.75)
In simulations, calculations are done in the dimensionless time t a = t ⌧ L . Between two collisions, integrating Eq. (3.74) between an initial time t = t 0 a and t a , we then update the position and the velocity of particles using:

u(t a ) = u(t a ) cosh(t a t 0 a ) + v(t 0 a ) sinh(t a t 0 a ) (3.76) x(t a ) = s i + u i (t a ) v(t a ) = v(t 0 a ) cosh(t a t 0 a ) + u(t s 0 ) sinh(t a t 0 a ). (3.77)

Next collision time

After each collisions between two particles, we need to compute the next collision time of each pairs of particles (i, i + 1). Starting from an arbitrary time t 0 a , the condition x i (t col ) = x i+1 (t col ) leads to the following polynomial in

X = e t col ( u + v)X 2 + 2n 1 L X + u v (3.78) where u = u i+1 (t 0 a ) u i (t 0 a ) and v = v i+1 (t 0 a ) v i (t 0 a ).
From the positive roots we deduce the next collision time. If their is no positive root, the next collision time is set to infinity.

VDM -1d self-gravity with viscous damping force

Equation of motions for periodic system

To simulate this model, we again use a code built to simulate an infinite periodic system. As for SGS, the code reproduces our model in the limit that the size of the system is negligible compare to the size of the periodic box. Like the code described above, this code was designed to simulate a toy model for cosmological structure formation, but now including an additional term to mimic the effect of universe in expansion (see section 2.1) which, in appropriate variables, is just a fluid damping term.

The equations of motion in this case are derived starting from [87]:

d 2 x i dt 2 + 2H dx i dt = g a 3 lim µ!0 X j6 =i e µ|x i x j | sgn(x i x j ) (3.79)
which is simply the 1d analogy of the equation of motion of 3d cosmological simulation where H = ȧ a is the "Hubble constant" function, and a(t) the scale factor. Again, using displacements from a regular lattice:

u i = x i ( 1 2 +i)n 1 L .
The equation of motion may be written as:

d 2 u i dt 2 + 2H du i dt = 1 a 3 2gn 0 (u i hui) (3.80)
with hui = 0 in practice in the simulation code.

Taking a(t) = (t/t 0 ) 2/3 with t 0 = 1 p 3gn

L an arbitrary time as in the reference "Einstein De Sitter" cosmology [START_REF] Sicard | Out-of-equilibrium dynamics in infinite one-dimensional selfgravitating systems[END_REF]87]), we substitute a in Eq. (3.80) where ⌧ L given in Eq. (3.75), we obtain:

d 2 u i dt 2 + 1 ⌧ L p 6 du i dt 1 ⌧ 2 L u i = 0 . (3.83)
We recognise the two terms of Eq. (3.74) describing an Hamiltonian SGS, with an additional term that account for the damping. We now suppose the system is confine in a region L 0 = ↵L in the center of the box and express Eq. (3.83) with the dimensionless time variable

t = t ⌧ dyn = 1 2 p 2↵ t ⌧ L , (3.84)
where ⌧ dyn is the time unit associated with the finite isolated system. and going back to x i = s i + u i coordinates; the last equation can be written:

d 2 x i d t2 + ↵ dx i d t = g[N > (x i ) N < (x i )] 2gN ↵x i , (3.85)
where

↵ = 2 p ↵/3 = ⌘⌧ dyn (3.86)
is the coefficient of dissipation in units of ⌧ 1 dyn . Again we see explicitly that in the limit ↵ ⌧ 1, the background term is negligible and the equation used in the code approximates arbitrary well the one of a finite system with fluid damping. We note also that as the dissipation leads to a system contraction, this approximation turns out to be more and more accurate with time. Moreover, in these units, the equations are independent of n 0 = N L 0 . Therefore with this method one can tune the parameter ↵ = n L n 0 and allows the experiment of different viscous coefficient ↵ .

Integration of the equations of motion.

Then, between two collisions, it is simple to integrate Eq. (3.81). Then, in VDM simulations, particles are moved forward up to the next collision time using the following expressions: and obtain the following quintic polynomial equation:

u i (t e ) = u i (t 0 e ) " 3
3 5 ( u + t 0 e v) Y 5 + n 1 L Y 3 + 1 5
(2 u 3 vt 0 e ) = 0 (3.88) where u = u i+1 (t 0 e ) u i (t 0 e ) and v = v i+1 (t 0 e ) v i (t 0 e ). The first derivative of the polynomial (3.88) is quadratic in Y 2 and the expression of its roots can be computed exactly as a function of its coefficients. With the use of the maxima found with the derivative, a standard algorithm is implemented to find systematically the roots of the polynomial. Adding the constraint t col e > 0 and t col e > t 0 e the roots are obtain using the dichotomy method which is always numerically stable [START_REF] Viot | Numerical Simulation in Statistical Physics Lecture[END_REF].

ICM -1d self-gravity with inelastic collisions

It is straightforward to modify the code for the 1d SGS so that it can also simulate the ICM. For the two previous model (SGS, VDM), a collision results in a simple exchange of the velocities of the particles. For this model, once the two particles i and i + 1 are moved at the same point using the equation of the 1d SGS without dissipation, we change their velocity according to the collision rules Eq. (3.13). The strength of the dissipation is tuned with = (1 c R )N/2, consistently with the previous discussion 2.2.

For both models with dissipation, the system contracts and the density becomes higher and higher, and the duration between two collisions decreases continuously. For the two dissipative models, the simulation is interrupted when the time between two collisions becomes smaller than the accuracy of the computer.

Initial conditions for numerical simulation

For the SGS, the system relaxes, in few times of ⌧ dyn , to a QSS. The nature of the QSS, i.e. its distribution function, depends on the nature probability distribution we use to generate initial conditions. In presence of dissipation, we expect the system to proceed in the same way through this phase of its evolution because of the separation of time scales. As our purpose is to examine the robustness such QSS, it is important to have a variety QSS, and thus variety of initial conditions.

We present, here, the three different initial conditions we have employed: the "rectangular water-bag" initial condition with different initial virial ratio, the microcanonical equilibrium, and another water-bag configuration called "Fermi degenerate".

-1 -Rectangular water-bag distribution A water-bag distribution is one which is single valued and non-zero only in a compact region of phase space. The rectangular water-bag initial condition in which this region is rectangular has been used to study the Vlasov dynamics [18,[START_REF] Worrakitpoonpon | Relaxation of isolated self-gravitating systems in one and three dimensions[END_REF]42]. It is simple, and macroscopic quantities can be easily computed.

To start from a rectangular water-bag distribution, we choose randomly and uniformly the (phase space) coordinates of N particles in the rectangular domain

D = [ L 0 /2, L 0 /2] ⇥ [ V 0 /2, V 0 /2].
The one particle probability distribution function is then:

P 1 (x i , v i ) = ⇢ 1 L 0 V 0 , if (x, v) 2 D, 0, otherwise. (3.89)
The normalization is such that the integral of P N (x 1 , v 1 , ..., x N , v N ) = ⇧ N i=1 P 1 (x i , v i ) over the 2N coordinates of phase space is unity.

Using this probability distribution function we can compute the mean initial kinetic energy over the realisations of the process,

hK 0 i = Z L 0 /2 L 0 /2 Z V 0 /2 V 0 /2 1 2 m N X i=1 v 2 i P N (v i )d N v i = 1 8 Nm V 2 0 3 , (3.90)
and the mean initial potential energy 3

hU 0 i = Z L 0 /2 L 0 /2 Z V 0 /2 V 0 /2 U (x i )P({x i })d N x i d N v i (3.93)
3 For the force given in Eq. (3.73) (with periodic boundary condition), the expression for the potential energy in coordinates where the centre of mass is zero is:

U (x i, t) = g X i,j2S X i6 =j |x i x j | g↵n 0 x 2 i (3.91)
and the mean initial potential energy is:

hU 0 i = gL 0 N 2 N 6 2g↵L 0 N 2 6 . (3.92)
Given that the size if the box is always much larger than the size of the system, the background term is negligible (↵ ⌧ 1) We have:

R 0 = 2hK 0 i hU 0 i = V 2 0 2g(N 1)L 0 (3.94)
R 0 is thus directly related to the aspect ratio of the rectangular water-bag. The initial energy E 0 is completely determined by R 0 . In the following results we use mainly two values of the initial virial ratio:

• R 0 = 1 for which 2K 0 ⇠ U 0 : "hot" initial state. • R 0 = 0.01 for which K 0 ⌧ U 0 : "cold" initial state. We do not consider the initial conditions where R 0 > 1 as there are known to relax into similar QSS as for R 0 < 1 [START_REF] Worrakitpoonpon | Relaxation of isolated self-gravitating systems in one and three dimensions[END_REF].

-2 -Equilibrium state

As we have discussed in chapter 1, in 1d, the thermal equilibrium of an open selfgravitating system is well defined. Further the micro-canonical and the canonical ensemble are equivalent. Maximising this in canonical or micro-canonical ensemble, one can compute the exact thermal equilibrium solution of the system in the meanfield limit [20,[START_REF] Worrakitpoonpon | Relaxation of isolated self-gravitating systems in one and three dimensions[END_REF](chapter 1):

f eq (x, v) = 1 2 p ⇡ 1 ⇤ sech 2 ⇣ x ⇤ ⌘ exp ✓ ⇣ v ⌘ 2 ◆ , (3.95)
with the two characteristic scales depending on the total energy of the system E:

2 = 4E 3mN , ⇤ = 4E 3g(mN ) 2 . (3.96)
For this solution of the Vlasov equation, the distribution of positions and velocities are separable and using the dimensionless variable ⇠ = x/ and ⌘ = v/ , we choose the N particle initial positions randomly with the probability distribution function

g(⇠) = 1 2 sech 2 (⇠) = ⇤⇢ eq (x); (3.97)
and the velocities randomly following the distribution 4 :

h(⌘) = 1 p ⇡ exp ⌘ 2 = ✓ eq (v). (3.101)
It has been demonstrated [28] that the thermal equilibrium is the only separable stationary solution of Vlasov equation. Moreover simulation have shown that finite N system relax toward this equilibrium on a time scale ⇠ N⌧ dyn ( [28] and reference therein.

-

-Fermi degenerate water-bag distribution

A consequence of Jean's Theorem [START_REF] Binney | Galactic Dynamics[END_REF], in 1d, is that any distribution which is only a function of e(x, v) = 1 2

mv 2 + m ¯ (x) , (3.102)
the individual energy of a particle, is a stationary solution of the Vlasov equation. In 1966, Hohl and Feix derived a non-equilibrium stationary solution of the 1d Vlasov equation [START_REF] Hohl | [END_REF]. The authors seek for a water-bag distribution of this kind. They sought a solution in which f is taken to be constant in the area A of phase space such that e(x, v) < e max , where e max is a maximum particle energy, and zero outside, i.e.:

f (x, v) ⌘ f (e) = ⇢ f 0 , if 0 < e  e max , 0, if e > e max . ( 3

.103)

It is then possible to derive an expression for the curve bounding A: defining v + (x) > 0 and v (x) = v + (x) < 0 to be the contour of the water-bag, symmetric 4 The random generator of the computer gives a random number rect [0,1] between 0 and 1 with uniform probability. In order to generate the spatial distribution we use the inverse transform method (e.g. [START_REF] Viot | Numerical Simulation in Statistical Physics Lecture[END_REF]). For this we use the anti-derivative of g of Eq. (3.97) which is G(⇠) = 1 2 tanh (⇠). Then compute it inverse:

G 1 (⇠) = 1 2 ln ✓ 1 + ⇠ 1 ⇠ ◆ . (3.98)
Then the distribution g(⇠) in position is generated by transforming rect [0,1] with G 1 :

⇠ = 1 2 ln ✓ 1+rect [0,1] 1 rect [0,1] ◆ (3.99)
For the Gaussian velocity distribution h(⌘), we make use of the Box Muller method [START_REF] Viot | Numerical Simulation in Statistical Physics Lecture[END_REF] and generate with the following expression: with respect to the x axes. Then the water-bag distribution solution of the stationary Vlasov equation may be written as:

⌘ = 1 p ⇡ q ln(rect [0,1] ) sin(2⇡rect [0,1] ) (3.100)
f (x, v) = f 0 ⇥(v v + (x))⇥(v v (x)), (3.104)
with ⇥(x) the step function (heavy-side). Then, inserting in the stationary Vlasov equation we obtain:

✓ v @ @x + F @ @v ◆ f 0 ⇥(v v + (x))⇥(v v (x)) = 0 (3.105)
with F the mean field force. Using the derivative of the distribution ⇥, the latter equation is equivalent to:

(v v + ) ✓ v @v + @x + F ◆ (v v ) ✓ v @v @x + F ◆ = 0 (3.106)
Then, integrating over v and using v (x) = v + (x) we obtain:

v ± @v ± @x F = 0 (3.107)
Further, F is given by the Poisson equation:

@F @x = 2g⇢(x) = 4gm Z v + (x) 0 f (x, v)dv = 4gmf 0 v + (3.108)
Inserting in Eq. (3.107), one gets the following differential equation for the contour:

@ 2 v + @x 2 + 1 v + ✓ @v + @x ◆ 2 4gf 0 = 0 (3.109)
This equation is non linear and must be solved numerically. Do to so, we employ an adaptive integration step method. Such an integration method is essential when the slope of the curve becomes singular. Then we fit the point obtained by this numerical integration with an ansatz function.

v a (x) = p a + bx + cx 2 + dx 4 + ex 5 (3.110)
where a, ..., e are constant. In Fig. 3.3 is shown the curve obtained numerically which describe the contour of the water-bag. Then we use the symmetry of the solutions (under x ! x and v ! v) to obtain the complete contour. To generate our initial condition we start by distributing uniformly point in a rectangular region enclosing the contour. Knowing the function v a it is simple to implement an acceptancerejection method to select N points inside the contour (Fig. 3.3).

Figure 3.3: The "Fermi degenerate" water-bag phase space distribution determined numerically. The blue crosses have been obtained in solving Eq. (3.109) and the red curve is a fit of the envelope with the anstaz Eq. (3.110).

In [START_REF] Hohl | [END_REF], the contour Eq. (3.109) is derived from an energy variational principle. They allow the contour to vary at fixed volume of phase space. Then looking for the water-bag configuration that minimises the total energy E, Eq. (3.107) is obtained. The Fermi degenerate distribution is thus the minimum energy configuration of the system for a water-bag distribution with a given phase space density.

Moreover, this state minimises the Lynden-Bell entropy discussed briefly in chapter 1 for a single level of phase space density. Indeed, this solution corresponds to the so-called degenerate limit ( ! +1) of the solution given by the theory of Lynden-Bell (chapter 1).

f (x, v) = 1 1 + e [e(x,v) µ]
(3.111) with µ = e max is the analogous of the Fermi energy in a fermionic gas. This is why this state is called "Fermi degenerate". Because of the property of stationarity ans stability, we found interesting to use this particular distribution as an initial condition of our simulations.

Simulation results for 1d self-gravity

We now present our numerical results. First, we reconsider the case of an isolated system and check that we reproduce results obtained in other studies. Indeed the dynamics of a finite 1d SGS has been the subject of numerous studies [19,28,[START_REF] Hohl | [END_REF]90]. In particular, it is important to verify that the numerical method used simulates accurately a finite system despite the presence of the background term. We will see that the simulations indeed reproduce the known results on violent relaxation toward QSS and on the thermal relaxation. Moreover, increasing the temporal resolution of observables we were able to exhibit the presence of "breathing modes" in the 1d SGS model. These modes have been analytically studied for SGS [91] and for other long-range interacting models in [5,[START_REF] Olivetti | Des systèmes à longue portée aux atomes froids[END_REF].

Second, we present the results of the VDM. We show that the decay law (Eq. (3.53)) fits perfectly the theoretical prediction from the kinetic equation. Moreover for the different initial conditions we show that the scaling QSS are a very good description of the dynamics.

Finally, numerical results of the ICM are displayed. For this model, the prediction of a system collapse in a finite time fits with the simulation data. However, for initial conditions with R 0 ⌧ 1, we observe significant deviations from the scaling solution at longer times. These deviations do not correspond to any tendency of the system to approach its thermal equilibrium, but on the contrary takes the system further from this latter state.

Macroscopic observables

We monitor the kinetic energy and the potential energy which is computed with Eq. (3.91). In addition, two macroscopic "order parameters" are convenient to study the dynamics.

The first one is the virial ratio already introduced above and defined by

R(t) = 2K(t) U (t) , (3.112)
where K and U are respectively the kinetic energy and the potential energy. A stationary solution of the Vlasov equation implies the virial ratio to be strictly 1.

Thus in a mean field (N = 1) description of the system, this ratio being unity is a necessary condition to be in a QSS. For our finite N system, we say the system is "virialised" when R is unity (up to small finite N fluctuations). The second observable is the separability ratio where h.i is an arithmetic mean over the N particles. This ratio is zero if and only if the system is in thermal equilibrium [START_REF] Worrakitpoonpon | Relaxation of isolated self-gravitating systems in one and three dimensions[END_REF]. Since 11 is defined as a ratio if the moments of the position and velocity distribution, in a QSS this ratio has a constant value.

Finally, we study also the spatial and velocity distributions. The histograms display below will represents the averaged statistic over many realisations of the dynamics to reduce finite-size effects.

SGS (without dissipation)

As described in previous studies and discussed in chapter 1, the dynamics of 1d SGS is described by a two stages of evolution during each of which the shape of the distribution function changes significantly. We describe first the results of the simulation for the shorter time ("collisionless", or "violent" ) relaxation. We then illustrate the long time relaxation toward the thermal equilibrium.

Violent relaxation to QSS Fig. 3.4 illustrates the complex dynamics of violent relaxation. Each panel shows the particle positions in phase space, at a given time. In each case the 100 realisations of the system are over-plotted. As time progresses the phase space structure "winds up" into a filamentary structure. As the dynamics goes on, the filaments stretch spiralling around the center and finally mix and become indistinguishable (last 2 panels). At a macroscopic level, the phase space distribution does not evolve significantly any more: the QSS is attained.

We now analyse the evolution of the macroscopic quantities introduced above during this early stage relaxation. In Fig. 3.5 the virial ratio is displayed as a function of time up to 150⌧ dyn . For both initial conditions, in ⇠ 10⌧ dyn , the parameter relaxes to the value 1 with fluctuations smaller than 10 2 . A blow-up of the oscillations is presented below for both cases. In the left panels we zoom in on the first few dynamical times; and in the right panels on later times once it has virialised. For both case, when the system is relaxing (left panels), the smallest period of oscillations is apparently exactly ⌧ dyn independently of the initial condition and of the number of particles.

Moreover for the hot initial condition the virial ratio oscillates manifesting beat like structure, suggesting the existence of two slightly different frequencies. One of the two frequencies is clearly the inverse of the time units: ⌧ 1 dyn = p gn 0 /2. Indeed we note that these features are not present for the cold case while the oscillations of period of ⌧ dyn are still present. After the system has virialised (right panels), the oscillations are still present with a period of order ⌧ dyn , increasing the number of particles does not appears to quench the fluctuations. Moreover, the results, averaged over 100 realisations of the initial distribution, explore largely the statistic of the initial rectangular portion of phase space. Thus it seems that these oscillations are representative of the (coarse grained) Vlasov dynamics.

Such "breathing modes" have been studied and seen in various kinds of simulations [5] and they have been analytically studied for self-gravitating system in [91]. These modes generate an oscillating potential with which some particles resonate. As discussed in chapter 1, in the core-halo theory for violent relaxation [42], this process is proposed to be the analogous to a Landau damping mechanism, responsible for the formation of the "halo" in the distribution displayed below. In Fig. 3.6 is presented a plot of 11 for the first 100 dynamical times starting from many different initial conditions. After ⇠ 10 20⌧ dyn , the parameter relaxes to a stationary value which depends on the initial condition considered. The values are consistent with the ones obtained in [18,[START_REF] Worrakitpoonpon | Relaxation of isolated self-gravitating systems in one and three dimensions[END_REF] for the rectangular water-bag case. For the water-bag initial states, which are not stationary solutions of the Vlasov equation, the relaxation of this parameter and of the virial ratio, clearly indicate the establishment of QSS. Also one remarks that as R 0 ! 0 the stationary value 75 11 increases. This indicates that colder initial conditions lead to a phase space distribution with a stronger correlation between coordinate and velocity space.

The Fermi degenerate and the thermal equilibrium initial conditions are realisations of a stationary solution of the Vlasov equation. We expect that they will not evolve macroscopically. This is indeed what is indicated by the observed behaviour

Relaxation to equilibrium

We briefly present results on much longer time evolution of this system. We can observe, as in other studies, a relaxation to thermal equilibrium. The physical process involving this relaxation acts on much longer time scale which recently have been shown to be of order of to ⇠ zN ⌧ dyn where z is a numerical factor of order 10 2 10 3 [28,[START_REF] Yawn | [END_REF]94]. The relaxation time has been shown to be strongly dependent on the initial virial ratio R 0 [28]. In Fig. 3.9 is shown an average over 50 realisations of a simulation of 128 particles starting from a cold intial state run for 5.10 5 ⌧ dyn . The relaxation takes place on the same time scale for order 10 5 ⌧ dyn , in good agreement with the results of [28] obtained using a different code. Figure 3.9: Thermal relaxation of a system of 128 particles (50 realisations). The energy and the virial ratio R and the separability ratio 11 , versus t/⌧ dyn for rectangular water-bag R 0 = 0.01.

VDM: SGS with viscous damping

The 1d SGS with viscous damping belongs to the large class of systems concerned by the theoretical predictions derived in section 3, corresponding to the VDM case with n = 1, and ↵ = 1 . For the latter case, the parameter of Eq. (3.53) is zero and we recall that, for this case, the scaling QSS admits the theoretical scaling:

v(t) = v 0 e t/tc (3.114)
for the velocities, where t c is the typical time scale of the dissipation given by Eq. (3.63). For this case therefore t c = 3 ⌘ and the energy E / v2 (t)

E(t) = E 0 e 2⌘t 3 = E 0 e 2 3 ↵ t ⌧ dyn (3.115)
where ↵ = ⌘⌧ dyn is the dimensionless parameter used to tune the dissipation in simulations (see section 4). We observe excellent agreement with the decay law: in Fig. 3.10 the total energy E, 3K and 2U are plotted versus t/⌧ mf , for the VDM with ↵ = 0.012. Initially, the kinetic and potential energy curves oscillate in phase opposition. These oscillation are damped, and in a short time the system virialises. In the lower panel, the semi-log plot the decay of the energies for the two cases benchmarks the decay law derived. The best fit is obtained with the parameter ↵ = ⌘⌧ mf to 10 4 . Figure 3.12: Time evolution of velocity (left panels) and spatial (right panels) distribution for VDM and ⌘⌧ dyn = 0.037. The initial conditions are a rectangular waterbag with R 0 = 1 (upper panels) and rectangular water-bag with R 0 = 0.01 (lower panels). Successive times t/⌧ dyn = 17, 35, 53, 70, 88, 106, 123, 141, 159, 176, 194, 212 are overplotted, the time arrow goes from blue to red. Simulations have been perform with 1024 particles and the results are averaged over 100 realisations of each initial condition.

The parameter

11 also relaxes in ⇠ 40⌧ dyn towards a stationary value for both cases R 0 = 1 and R 0 = 0.01. These two parameters are independent of the system size; Fig. 3.11 illustrate the presence of scaling stationary states. In the inset, we can see the frequency of oscillations of the virial ratio increases. Indeed, if L(t) is the typical size of the system at time t, the mean field time scale of the system ⌧ mf (t) = 2 p L(t)/gN decreases as the system gets smaller. Thus the frequency of the oscillations / 1/⌧ mf increases. We recall that the time unit used is ⌧ dyn = 2 p L 0 /gN with L 0 the initial length of the water-bag. In Fig. 3.12 is plotted on top of each other the normalised position and velocity distributions at different times respectively versus x/x(t) and v/v(t) where x and v are the measured standard deviations of the distributions. The superposition of the curves illustrates the accurate description of the evolution of the system by scaling QSS. The blue curves correspond to the early time when the system has not virialised yet, while the red curves represent latter times.

At early time, the behaviour of 11 and of the virial ratio during the period of violent relaxation does not differ a lot from the one of the SGS without dissipation and indicates again, up to a rescaling of the variables, the two models have the same evolution.

Figure 3.14: Time evolution after violent relaxation of velocity (left panels) and spatial (right panels) distribution for ICM and N = 0.01. We used rectangular water-bag initial condition with R 0 = 1 (top panels) and R 0 = 0.01 (N = 4096) (bottom panels) Time goes from blue to red at indicated times. Each simulations have been run for 100 realisations of initial conditions. the system virialises and reaches a (macroscopically) stable state. On longer time scales (t = 70, 140, 212, 280 ⌧ dyn ), the four bottom panels allows one to see clearly a contraction of the system due to the dissipation.

Figure 3.14 shows the evolution of the position and velocity distribution of ICM after that the system has virialized starting from a hot and a cold RWB initial condition. In order to compare the successive distributions, the spatial and velocity coordinates have been rescaled by the estimated variances x(t) and v(t) of the distributions. For RWB initial condition with R 0 = 1 (top panels), after the violent relaxation, the system is perfectly described by a scaling solution. For cold RWB initial condition (bottom panels), the scaling solutions gives an accurate description of the dynamics only at early times. At latter times one can see clear deviations from such a solution, with the distributions becoming increasingly peaked at their centres.

To be more precise, we now examine the simulation results for the order parameters introduced above. In Fig. 3.15 the energy, the virial ratio and the separability ratio 11 are displayed. Three simulations start from different RWB initial conditions. For all initial conditions, the system collapses in a finite time and the energy decay, fitted by (1 t tc ) with = 2.00 ± 0.01, matches the theoretical prediction (even in the case of the deviation from scaling solution).

In the bottom panel of Fig. 3.15, we observe that, for all initial conditions, the system virialises in ⇠ 40⌧ dyn . In the inset is a blow up of the virial ratio (for R 0 = 1) at early times. As for the first two models, the virial ratio manifests beats like features for the hot RWB initial condition. One also notes, as in the VDM, an increase of the oscillation frequency of the virial ratio, due to the contraction of the system as explained beyond.

In the middle panel of Fig. 3.15, for initial RWB R 0 = 1 (blue), the 11 stabilises in few ⌧ dyn and then remains constant. This confirms that in this case the dynamics of the system is well described by a scaling QSS. For the two colder initial conditions (R 0 = 0.1, 0.01), the 11 stabilises but then deviates from the constant behaviour. For these cases, as we have seen, scaling solutions are indeed an accurate description of the system only for a transient period and this is also reflected in the behaviour of 11 .

Figure 3.16: For the ICM, reduced total energy E(t)/E(0) versus t/⌧ dyn , for the indicated values of and R 0 = 0.01. The inset shows a zoom on the shorter time evolution.

In Fig. 3.16 we plotted the energy decay curves for different values of as a function of the rescaled time t/⌧ dyn . With this change of variables the curves superpose and, as expected by the theoretical treatment (cf. Eq. (3.63)), the collapse time scales as t c ⇠ 1 just as the time scale ⌧ diss associated with the dissipation. We also note that the collapse time t c depends slightly on the initial conditions (Fig. 3.15). The inset of Fig. 3.16 shows small deviations from the scaling behaviour at short times, associated with the virial oscillations during the initial violent relaxation.

The bottom panel of Fig. 3.17, shows a semi-log plot of 11 as a function of t, for = 0.01, 0.005, 0.001. We see that the early time relaxation is not affected by the variation of the strength of the perturbation up to ⇠ 10⌧ dyn . Indeed choosing sufficiently small, the two time scales of the system are well separated (⌧ dyn ⌧ ⌧ diss ). Thus the effect of the dissipation happens on longer time scale, as expected.

The upper panel of Fig. 3.17 shows the evolution of 11 for the ICM starting from R 0 = 0.01, R 0 = 0.1 and R 0 = 1, and each for the three different values of . As well as the energy, the evolution of 11 , for all initial conditions, collapses using the rescaled time t/⌧ dyn . This seems to indicate that the dependence of the dynamics on might be scaled out using this time coordinates. For initial conditions with R 0 = 0.1, 0.01, the deviation from the scaling solution leads to an increase of 11 , indicating that the system evolves into a state in which spatial and velocity coordinates are more and more correlated. However during this evolution, the virial ratio remains equal to 1 up to small fluctuation, indicating that the system remains virialised at each time. Thus we can infer that throughout the deviation from the scaling solution, the system follows, nevertheless, a succession of (quasi-)stationary states of the Hamiltonian dynamics (Vlasov dynamics). In order to verify this statement, we run a simulation in which the inelasticity of the collisions has been "turned-off" after time t of f = 110⌧ dyn . The results are displayed in Fig. 3.18. As we observe in the inset, the energy is conserved after t of f . The constant value of 11 after t of f indicates clearly that the evolution of the system throughout follows a succession of quasi-stationary state, exactly like an adiabatic process follows a succession of equilibrium states in thermodynamics. In other words, the inelastic collisions change the state of the system, but the long-range dynamics constrains it to remain in a virial equilibrium. 11 as a function of t/⌧ dyn for ICM with = 0.01. In this simulation the dissipation has been stopped at time t = 100⌧ dyn . The inset shows the total energy versus time In Fig. 3.21 is shown the evolution of 11 for all the different initial conditions we have used. We observe that for the thermal equilibrium initial condition, the scaling QSS gives an accurate description of the system evolution at early time. At later times, for this case, the deviation from the scaling dynamics measured by

In this chapter we have studied the effect of various kinds of dissipation on the dynamics of long-range interacting systems. For a Hamiltonian system, in the mean field limit, the dynamics is described by the Vlasov equation. This non linear equation for the phase space density admits an infinite number of stationary solutions. The stationary solutions of the Vlasov description are quasi-stationary states of the long-range particles dynamics.

For a 1d long-range interacting system, we showed that the QSS are robust to a large class of dissipative effects including fluid damping and, to lesser extent, inelastic collisions. This robustness has been characterized theoretically by the existence of scaling solutions of the appropriately modified Vlasov equation. This equation, derived under mean field and molecular chaos assumptions, has the same structure as the Vlasov equation with additional terms accounting for the dissipation. The existence of scaling solutions was derived under the assumption of separation of the two central characteristic time scales of the system: ⌧ diss , the time scale of the dissipation and ⌧ dyn , the time scale of the mean field dynamics. Under this assumption, the system dissipates energy sufficiently slowly that the system always remains very close to a QSS. After the initial phase of violent relaxation lasting about 40⌧ dyn , in which the QSS is established, the system undergoes a contraction on the longer time scale ⌧ diss ⌧ dyn . Under a well defined rescaling of the position and the velocity variables, the system keeps the same QSS distribution throughout the contraction.

The existence of scaling QSS has been confirmed by simulation of the N body long-range dynamics for two models of dissipation in the case of 1d gravity. For a fluid dissipative model, where the dissipation takes place homogeneously in the system, scaling QSS are found to be a very good description of the dynamics for any initial condition. However for a model with inelastic collisions, the dissipation depends on the frequency of collisions and thus on the density of the system. In this case, for inhomogeneous QSS, the dissipation is more important in the denser region of the system. In coordinates following the contraction of the system, this implies an effective energy transfer in the system, a change of phase space distribution and the system leave the scaling solutions on time scale of order of ⌧ diss . We have shown that this evolution drives the system towards more correlated states between position and velocity distributions, and thus further away from the thermal equilibrium.

Furthermore, we have shown the perturbed system remains at all times virialised. Thus, by switching the dissipation off once the system has left the scaling solution, we have shown that the system stays in the more correlated state indicating that this increase of the correlation induced by inelastic collision is indeed an adiabatic evolution through quasi stationary states of the isolated system. Hence, this collisional perturbation of the dynamics leads to a whole range of stationary solution of the dynamics of the unperturbed system.

The inelastic collisions added to the 1d system induces local contractions of the phase space volume occupied by the system and makes it collapse. Moreover, following this collapse in dimensionless variables, the inelastic collision induces local change of velocity distribution, while a 1d elastic collisions does not change it. For the quasi-elastic limit of the model, each collision induces an infinitesimal and local change in the velocity distribution. Even if this change is purely deterministic microscopically, on a coarse-grained level one would expect to be able to treat it as a stochastic process. However, the contraction of the system due to the energy dissipation makes it very difficult to explore further the evolution away from the scaling QSS, and in particular whether the perturbation drives it to a new stationary state at some time.

In the next chapter, we pursue the question raised by this interesting evolution toward more correlated QSS. To do so we introduce a similar class of perturbation arising from collisions in a long-range system, but such that the energy is conserved on average.

Chapter 4

Long-range systems with internal local perturbations

Introduction

In the previous chapter we have studied the effects of dissipation on a long-range system and, in particular, on QSSs which are one of their most striking feature. In the limit of the separation of the relevant time scales, a kinetic equation, derived for a large class of one dimensional long-range models, enabled the prediction of the existence of "scaling QSS". Evolving in such a solution the systems keeps the same phase space distribution up to a rescaling, and thus the dissipation has a somewhat trivial effect on the QSS which would appear in the isolated system. Simulations of a 1d self-gravitating system with different sorts of dissipation have largely confirmed that such states do indeed describe well the evolution of these systems starting from a range of different initial conditions.

However, in one of the two models (the ICM), we have clearly seen evolution deviating from scaling QSS, and interestingly, as we have noted, this evolution brings the system toward more correlated states rather than toward thermal equilibrium. Indeed, on the basis of work in the literature [95,96] (see chapter 1) we might have anticipated that the "internal" stochasticity of the model might lead to a relaxation to thermal equilibrium. To further our study of this observation -and, more generally, to explore how different kinds of perturbation affect the dynamics of long-range systems and, in particular, their QSS-we study in this chapter two different toy models.

The first model is a variant of the ICM in which particles undergo not just inelastic, but also super-elastic collisions (i.e. in which there is a net gain rather than a loss of energy). As we will detail, this is a generalisation of a granular model studied in literature [97]. In contrast with the ICM, the model has a dynamics which is explicitly stochastic: we chose randomly whether the collision is inelastic of super-elastic. Further we choose, essentially for simplicity, to conserve the energy on average (over the realisations). Indeed, in the ICM, as dissipation leads to a continuous contraction of the system, it becomes difficult numerically to follow the evolution of the system on a long time scale. This model can be considered as a simple toy model for real long-range interacting systems in which there are additional short-range interactions which come into play when particles are close to one another-and which lead to energy loss or gain to other degrees of freedom (e.g. through radiation or to internal degrees of freedom). Superelastic collisions can model thermally activated nanoclusters [98,99] and have also been reported for large-scale magnetized plasmoid in the heliosphere [100]. Granular models incorporating a random coefficient of restitution have been used in [97] to study the effect of energy injection of a vibrated 2d granular gas. In the latter study, more generally, different distributions of the values of the coefficient of restitution are considered. Further we note that such model could be considered to study the effects of errors introduced in the numerical integration of a long-range system. For example in 3d self-gravitating system the velocity change of particles when they pass very close to one another, in a"hard" two body collision, is known to be associated often with non-physical energy change due to poor integration of the trajectories [101,102].

The second model we consider is based on another granular model [103]. It is another variant of the ICM in which energy is injected by collisions with an additional constant velocity "kick". Again we consider a one dimensional self-gravitating system subjected to such collisions. In this case the system is expected to be able to attain a stationary state in which the energy is constant. This model in contrast with the first one is deterministic like the ICM.

We describe both models as a long-range systems subjected to "internal" and "local" perturbation: internal in the sense that the perturbation is associated with particular events in the system i.e. interaction between elements internal to the system (rather than with external agents); and local because the perturbation is localised in space and depends only on the state at this point of the system (and not on the state of the whole system). As we will discuss in our conclusion, we believe that this locality is central to the very different behaviour we observe compared to the model of [55] showing relaxation to thermal equilibrium.

The chapter is organised as follows. In the next section we consider the first model described above. We derive for it, under the assumption that the total energy is conserved, a modified Vlasov equation which accounts for the effect of the stochastic collisions in the quasi-elastic (perturbative) limit. We use it to show in particular that the thermal equilibrium of the 1d self-gravitating system is unstable when the perturbation is present. We then study the evolution of the energy of a single realisation, and show, using heuristic arguments, that, at any finite time, the energy variance, goes to zero as N does to infinity. This result means that the assumption of the energy conservation is indeed valid in the mean field -quasi-elastic limit. Numerical results for the model are then presented, first for the model in the absence of gravity and then with gravity

In the following section we investigate the second model through an essentially numerical study. In the case of the one dimensional granular gas (i.e. without gravity) we show that the system relaxes toward a bimodal velocity distribution. Further, the simulations reveal hydrodynamical modes that have been seen for this model in higher dimensions [103]. For the self-gravitating model we find that the system reach a steady state, independent of the initial conditions which evolves at longer time to state which is highly ordered microscopically.

The stochastic collision model (SCM)

After giving the microscopic definition of the collisions in the model, we develop a kinetic theory to describe its macroscopic evolution. We show that, in the quasielastic limit, the collisional operator added to the Vlasov equation can be written in a Fokker Planck form. This kinetic equation allows us to determine the large N scaling of the parameter characterising the perturbation in order to obtain a consistent mean field -quasi-elastic limit. An analysis of the moments of the collisional term is then presented. This permits us to write a dynamical equation for the moments of the phase space density function. From these equations, we deduce, in particular, the instability of the thermal equilibrium distribution. Finally, by an analysis of the random process undergone by the energy of a single realisation, we establish the validity of the assumption we made to build the kinetic equation. These results are general and valid for the model both with and without long-range interactions. In a second part of this section numerical results for the model with and without gravity are presented and discussed.

Microscopic description of the stochastic collisions

The stochastic collision model (SCM) is defined as follows. Denoting v ij the relative velocity of particles i and j with velocities v i and v j prior to the collision, the postcollisional velocities v ⇤ i and v ⇤ j are given by:

v ⇤ i =v j + 1 ↵ 2 v ij (4.1) v ⇤ j =v i 1 ↵ 2 v ij ,
where the coefficient of restitution ↵ is chosen randomly with equal probability from the two values c R and cR related by

cR = q 2 c 2 R . (4.2)
This relation ensures that the magnitude of the energy change, for a given relative velocity at collision, is the same in all collisions. The sign of the energy change is determined by whether we choose c R or cR :

• For an inelastic collision we use the usual restitution coefficient 0  c R < 1, such that the energy change after collision is

K inel = m 1 c 2 R 4 (v ij ) 2  0, (4.3)
• For a collision involving a gain of energy, so called super-elastic collision, we define, as a function of c R , the corresponding super-elastic coefficient of restitution cR = p 2 c 2 R such that the energy change after collision is

K superel = K inel 0 . (4.4)
With this definition, 1 < cR  p 2, and the quasi-elastic limit is such that both c R and cR go to unity.

We note that for this model the Jacobian of the transformation of the phase space volume element is the random variable ↵ and takes a value either less than or greater than unity, i.e. for a collision between two particles 1 and 2, the measure of the phase space dx

1 dx 2 ...dx N dv ⇤ 1 dv ⇤ 2 ..dv N = ↵dx 1 dx 2 ...dx N dv 1 dv 2 .
.dv N , is correspondingly either smaller or greater than before the collision1 .

Kinetic equation

In the previous chapter, the modified Vlasov equation for a model with inelastic collisions was derived from a generalised Liouville equation. We have truncated the BBGKY hierachy, neglecting the two point correlation function and obtained the modified Vlasov equation:

@ t f (x, v, t) + v.@ x f (x, v, t) @ x ¯ [f ](x, t).@ v f (x, v, t) = J [f, f ](x, v, t) . (4.5)
where the collision term J [f, f ](x, v, t) is the rate of change of the phase space density resulting from the collision rules. We saw in chapter 3 that it can be split into a gain and a loss term:

J [f, f ](x, v, t) = G (x, v, t) L(x, v, t) . (4.6)
The loss term measures the decrease of the population of particles inside a volume dxdv in the phase space around point (x, v) during the time interval dt:

L(x, v, t) = N M f (x, v, t) Z dv 0 |v 0 v|f (x, v 0 , t) . (4.7)
and the gain term i.e. the increase of the population in an infinitesimal volume during dt is:

G (x, v, t) = N M Z Z dv 0 dv 00 |v 0 v 00 |f (x, v 0 , t)f (x, v 00 , t) (v v 0 (1 )v 00 ) .
We define for convenience

2 = (1 c R ) 2 . ( 4.8)
Given that (1 ) = (1 c2 R )/4, the corresponding ˜ for a super-elastic collision, is a solution of the equation (4.4)

K superel (˜ ) = K inel ( ) () (1 ˜ )˜ = (1 ) . (4.9)
and therefore:

˜ = 1 2 r (1 ) + 1 4 . (4.10)
Subsequently, we define the collisional operator for super-elastic collisions as:

J ˜ [f, f ](x, v, t) = G ˜ (x, v, t) L[f, f ](x, v, t) . (4.11)
We have derived a collisional operator for the two kinds of collisions the particles are subjected to. In the microscopic dynamics, collisions happen at discrete time intervals and instantaneously. The equation (4.5), for the ICM, implicitly involves a coarse-graining in time on a time interval dt in which many collisions occur in a phase space element. This means in particular (as we will verify below) that we have implicitly assumed energy conservation. This assumption is valid exactly in the model only if we consider an average over realisations. We will see below that it also be taken to be valid for a single realisation in the quasi-elastic limit of interest to us. With these assumptions which will be discussed further in section 2.7, the kinetic equation for the SCM is written:

@ t f (x, v, t) + v.@ x f (x, v, t) @ x [f ](x, t).@ v f (x, v, t) = J SCM [f, f ](x, v, t) (4.12) . = 1 2 (J (x, v, t) + J ˜ (x, v, t)) .

Kramers Moyal expansion of the collision operator

We derive now the expression for the collisional operator in the quasi-elastic limit using a Kramers Moyal expansion of this operator. We take the limit of small and expand the two gain rates G ⇤ to second order respectively in and ˜ . Doing so, we obtain:

G = N M  f (x, v, t) Z |v 0 v|f (v 0 )dv 0 + @ v Z |v 0 v|(v 0 v)f (v)f (v 0 )dv 0 (4.13) + 1 2 2 @ 2 v Z |v 0 v|(v 0 v) 2 f (v)f (v 0 )dv 0 , (4.14)
and the same expression for G ˜ but with ˜ instead of . With the use of Eq. (4.9), one has the following relation:

+ ˜ = 2 + ˜ 2 = 2 + o( 2 ) . (4.15)
Gathering terms, and using the last relation, we find that in the quasi-elastic limit, the collisional operator takes a Fokker Planck form:

J SCM [f ] = 2 N M  @ v (f (x, v, t)a[f ](x, v, t)) + 1 2 @ 2 v (f (x, v, t)D[f ](x, v, t)) , (4.16)
with an acceleration and diffusion coefficient which are functions of the distribution function:

a[f ](x, v, t) = Z |v 0 v|(v 0 v)f (x, v 0 , t)dv 0 ; (4.17) and D[f ](x, v, t) = Z |v 0 v|(v 0 v) 2 f (x, v 0 , t)dv 0 . (4.18)
i.e. in the limit in which the time scale ⌧ SCM ⇠ 1 2 N ⌧ dyn is long compared to the mean field time scale ⌧ dyn , so that the system always remains very close to virialized.

Using results of the last section, we then obtain:

dhv 2 i dt = 0, (4.36)
i.e. then the kinetic (and also total) energy is conserved. We also have,

dhv 4 i dt = 2 (1 ) 2 N 2M Z |v 1 v 2 | 5 f (x, v 1 , t)f (x, v 2 , t)dxdv 1 dv 2 . (4.37)
The right hand side of Eq. (4.37) is positive definite and thus it predicts that hv 4 i increases monotonically in time. We infer that if the system reaches a stationary state hv 4 i must diverge in this case. Indeed we have implicitly assumed, in our derivation, that the integral on the right-hand side of Eq. (4.37) is finite and thus Eq. (4.37) can be used only if the velocity distribution decreases faster than v 6 .

Instability of thermal equilibrium

In this section we apply Eq. (4.37) to the particular case where the initial distribution is the thermal equilibrium. This shows the important result that the equilibrium of the isolated system is not one for our perturbed system. More specifically it gives us a prediction which we will be able to compare to the results from simulation.

Inserting the expression of the thermal equilibrium, given in chapter 1, for 1d selfgravitating system, on the right hand side of Eq. (4.37),

dhv 4 i dt = 2 (1 ) 2 N 2M Z |v 1 v 2 | 5 f eq (x, v 1 , t)f eq (x, v 2 , t)dv 1 dv 2 (4.38) = 2 (1 ) 2 N 2M Z ⇢ 2 eq (x)dx Z |v 1 v 2 | 5 e ( v 2 1 +v 2 2 2 ) 2 ⇡ dv 1 dv 2 = 2 (1 ) 2 N 2M 8 p 2 p ⇡ 5 Z ⇢ 2 eq (x)dx , where Z ⇢ 2 eq (x)dx = 1 4⇤ 2 Z sech 4 ( x ⇤ )dx = 1 3⇤ . (4.39)
On the other hand, for the granular case, if we consider a uniform spatial distribution with a Gaussian velocity distribution the calculation give the same result but for the spatial distribution:

Z ⇢ 2 eq (x)dx = Z +L 0 /2 L 0 /2 1 L 2 0 dx = 1 L 0 . (4.40)
We introduce the standard 4 th moment or kurtosis of the velocity distribution:

2 = hv 4 i hv 2 i 2 (4.41)
because we set v = hvi = 0 without loss of generality. This parameter, used in the simulations below, give a measure of the deviation from a Gaussian velocity distribution for which 2 = 3. For sufficiently short times we can deduce from Eq.(4.36) that hv 2 i(t) = C st . = 2 , then integrating Eq. (4.38) from 0 to t we predict the early evolution of 2 for this perturbation starting from equilibrium solution:

2 (t) = hv 4 i hv 2 i 2 = 3 + t ⇥ 2 (1 ) 2 N 2M 8 p 2 p ⇡ 5 Z ⇢ 2 eq (x)dx (4.42)
This last equation Eq. (4.42)should be a good approximation up to a finite time when the distribution deviates from the thermal equilibrium significantly.

Evolution of energy

In order to build the kinetic theory, we have assumed implicitly that, at each time interval dt, the kinetic energy gain and lost by the system cancels i.e. the energy is conserved exactly. Such a conservation of energy is only, evidently, valid in the ensemble average over all realisations of the stochastic dynamics while the energy in a given realisation does vary. However we show here that the variance, in the large N mean field -quasi-elastic limit, at any fixed time, goes to zero. Thus, in this limit we can expect the kinetic theory to be valid for a single realisation. We label, the successive collisions, k = 1...k t , with k t the number of collisions up to a time t. Defining = 1) , the k th collision (at time t k ) induces the microscopic energy change:

1 c 2 R 4 = (
E k = ⌘ k ( v k ) 2 (4.43)
where ⌘ k = ±1 are uncorrelated random variables and v k = v ij = v i v j where i and j the two particles involved in the collision k. In this model the energy thus performs a random walk with varying step-lengths. We now assume that (i) successive steps are uncorrelated i.e. h v k v 0 k i = k,k 0 for any two collisions k and k 0 ; and (ii) the velocities of colliding particles are uncorrelated, and more specifically that: h v 2 k i = 2hv 2 i(t k ). With these assumption we have a relation between the mean squared relative velocity and the kinetic energy (and total) energy of the system:

h v 2 k i = 4 K(t k ) N = a E(t k ) N ; (4.44)
where a is a numerical factor depending on the model considered:

• For a granular system, the total energy equals the kinetic energy and a = 4

• For a long-range system with a pair potential ⇠ 1/r n , we assume the system remains virialised, and a = 4n n+2

With these assumptions we have simply:

E k = ⌘ k a N E k , (4.45)
for the k th collision, where E k is the energy of the system at the time of the collision. It follows therefore that up to collision k t :

E kt = E 0 kt 1 Y k=0 (1 + ⌘ k a N ) . (4.46)
Taking the averaged over the ensemble of realisations (with ⌘ k = ±1 and uncorrelated) we obtain: h(1 + ⌘ k a /N ) >= 1 and

D 1 + ⌘ k a N 2 E = (1 + a 2 2 /N 2
) Thus, averaging Eq. (4.46) we get, hE kt i = E 0 , (4.47) i.e the averaged energy is indeed conserved and ⌧ 1, we have:

hE 2 kt i = ✓ 1 + a 2 2 N 2 ◆ kt E 2 0 . ( 4
hE 2 kt i = e ktlog ✓ 1+ a 2 SCM N 3 ◆ ' e kt a 2 2 SCM N 3 . (4.49)
Considering now a time scale on which energy perturbation (as measure by the variance) is small i.e.

k t a 2 SCM N 3 ⌧ 1. (4.50)
On such a time scale we can now estimate k t the number of collisions, in the approximation that the macroscopic evolution is unchanged by the perturbation: for the case without gravity we have k t / N⌧ granu (since ⌧ granu is the typical time between collisions), while for the case with gravity k t / N 2 ⌧ dyn (since ⌧ dyn is the typical time for a particle to cross the system). Thus for the first cases we need that t such that

a 2 2 SCM N 2 t ⌧ granu ⌧ 1 , (4.51)
for the energy to be conserved and for the case with gravity we need:

a 2 2 SCM N t ⌧ dyn ⌧ 1 . (4.52)
Thus in the large N limit, the variance of the energy can indeed always be taken to be arbitrary small and the kinetic theory developed should hold even for a single realisation (and not just for the ensemble average).

Correspondingly, for a large but finite number of particles N , or equivalently at small but fixed , the assumption of energy conservation in a realisation is a reasonable approximation up to a finite time t N , which scales as t N ⇠ N 2 ⌧ granu for a self-gravitating system and as t N ⇠ N⌧ dyn for a granular system.

Moreover, the variance is expected to grow linearly with time at short times after which it should grow exponentially. Therefore to run to longer time in the mean field limit (the limit of the kinetic theory) of the model, one needs to increase the system size.

In conclusion, the kinetic theory developed, and its results, are expected to hold, at least up to a finite time t growing with the particle number and this justify numerical simulation which we present the results in the following section. Finally, from Eq. (4.46), we can also deduce the full energy distribution. Considering the logarithm this equation is:

ln(E kt ) = ln(E 0 ) + kt 1 X j=1 ln(1 + a⌘ k ), (4.53)
i.e. the process is a biased random walk of step ln(1 ± a ). In the large N limit, we can neglect the small bias; indeed as ⇠ 1/N , we have ln(1 ± a ) ⇠ ±a . Then, in the large N limit, in the logarithm of the energy, the dynamics in energy space is a simple random walk. Therefore, the distribution P(E 0 |E, k t ) of the energy at collision k t , starting form energy E 0 , is log-normal:

P(E 0 |E, k t ) = E 0 E p 4⇡Dk t e ln 2 (E/E 0 ) 4Dk t , (4.54)
where the coefficient of diffusion is, for ⌧ 1, The superposition of the curves indicate, that the time evolution of the system is rescaled by 2 in agreement with Eq. (4.37). At early times, the linear behaviour of Eq. (4.42) is verified, while at longer times 2 apparently reaches a independent "plateau" plus residual large fluctuations. The left hand side of Eq. (4.37) being positive, one would expect, a monotonic growth of the kurtosis. 100 realisations of a 8192 particle simulation. After 20 ⌧ granu , the system apparently evolves into a state with a long-tailed velocity distribution. On the log-log plot on the left panel we find a good fit of tail of the distribution with ⇠ v  with  ' 4, 5 ± 0.5 which is in accordance with previous results [97]. In spatial coordinates (which we do not display) the system remain homogeneous as expected.

D = a 2 2 2 . ( 4.55)
For such a distribution, the 4 th moment is not defined and this is a very consistent explanation of the plateau behaviour of the parameter 2 (as a saturation phenomenon). Indeed for a finite number of particle (or realisations), the tails of the distribution is well sampled only up to finite value and the expected divergence of 2 is thus regulated by the finite particle number.

In summary, we have seen that for the SCM without gravity, the kinetic theory may be used to describe the averaged behaviour of a large (but finite) particle number system up to a finite time. An initial Gaussian velocity distribution is unstable in presence of these stochastic collisions, and in a short time, the system reaches an homogeneous stationary state with velocity distribution with a tail scaling as v 4 . The results are in good agreement with the previous studies [97]. Clearly the fluctuations of the virial ratio are below 1 (corresponding to a system which is virialised) and these fluctuations decrease in value and in amplitude as N increases. Thus as N ! +1, in the quasi-elastic limit, the system remains always exactly virialised and this is a strong indication that it evolves through a family of QSS. Fig. 4.6 shows the evolution of the separability ratio 11 . The upper panel shows it for two realisations and for the average over 300 realisations.

11 evolves on the same time scale as the variance of R, and reaches (as the variance of R) a well defined averaged "plateau" at longer times. The fluctuation from realisation to realisation of 11 are very large and their standard deviation is shown in the lower panel for N = 128 to 512. At earlier times it appears that the standard deviation decreases as N increases but at longer times it appears to be independent of N . This means that the system reaches a stationary virialised state but it is not clear that it converges to a well defined single value of 11 as N ! +1. This is further shown in Fig. 4.7 which shows for the same time axes, the evolution of both 11 (top) and 2 (bottom) for different N . In the lower panel is also plotted the theoretical prediction from the kinetic theory calculated in section 2.6.

We see therefore that the early time behaviour is independent of N as predicted by the kinetic theory. However at long times the evolution depends on N . In particular, the value of both 11 and 2 at "plateau" appears to increase with N , suggesting that they will diverge as N ! +1. This behaviour is an indication, just as we saw in the case without gravity, of a long tail distribution in the final stationary state.

In Fig. 4.8 where the evolution of both spatial and velocity distribution is displayed, we see that there is indeed such a behaviour. The tail of the velocity distribution scales as v  , with  = 3 ± 0.5, which appears slightly less step that from the case without gravity. The spatial distribution relaxes also toward a distribution with a power law tail. Though further numerical study would be required to establish this. Thus like in the case without gravity, the system is driven by the perturbation to change its velocity distribution dramatically, but in the presence of gravity this leads also to a complete reorganization of its spatial distribution. To check whether the virialized state is a QSS i.e. a stationary solution of the isolated system, we "turn-off" the stochastic perturbation at a given time. In Fig. 4.9 we present a simulation of 1024 particles with SCM = 0.1 up to t of f = 300⌧ dyn (i.e. for t > t of f , SCM = 0). In the top panel two realisations (in color) and the averaged energy over 100 realisations (in black) is shown. The gray lines after t of f are the (averaged) energy with the perturbation plotted as a comparison. We see that the energy is constant after t of f as expected. The two bottom panels show the separability ratio 11 (middle) and the kurtosis 2 (bottom) with the same color 115 convention. Up to t of f , R and 11 follow the evolution described above. After t of f , both parameter remain constant up to small fluctuation that are reduced by the average over realisation (black). We deduce that after this time, the system, now isolated, remains in the QSS to which the perturbation has driven it. As indicated by the constancy of the virial ratio, this confirms that evolution of the system is through a family QSSs.

Dependence on the initial conditions

Our results so far have been given for the thermal equilibrium initial condition, and indicate that at long times the system reaches a final QSS described by a long tailed velocity distribution. To see how this evolution depends on the initial conditions, we also have simulated the model starting from some of the different initial conditions described in chapter 3. Fig. 4.10 shows the evolution of 11 over 225 realisations, for a series of 128 particle simulations of SCM for SCM = 0.02, starting from thermal equilibrium, cold water-bag and hot water-bag initial conditions. Simulation have been performed over the time scale of the thermal relaxation of an isolated SGS described in chapter 3. The evolution of the 11 for an isolated SGS is also plotted as a comparison. We see that for all initial conditions, the 11 reaches the same "plateau" and never shows, even on long time scale, any tendency to relax to thermal equilibrium. The non equilibrium stationary state reached appear to be an attractor of the dynamics. Such a "stochastic attractor" have been reported in [52], in the context of a stochastically perturbed HMF.

Conclusion on SCM

Our numerical study of the SCM with and without gravity shows an early time behaviour in good agreement with the kinetic theory developed to describe it. In the presence of gravity, in the mean-field and quasi-elastic limit, the system, starting from any initial condition, evolves on a N independent time scale through a family of virialized state, each of which is a QSS of the unperturbed long-range system, until it reaches an apparently unique virialized state. This state is not the thermal equilibrium of the isolated model, and indeed is typically "further away" (in terms of correlations measured by 11 ) from the thermal equilibrium. Therefore we observed compelling evidence for the establishment of a "universal" non-equilibrium stationary state in this model. Despite the stochasticity of the dynamics, there is no tendency of the system to relax toward equilibrium. This clearly contrasts with the results obtained with a toy model in [95] constructed to test for the robustness of QSS to stochastic perturbations.

The BSR model

To go further in our study of the effect of "internal" perturbations on the dynamics of long-range systems, and in particular on their QSS, we have considered a second model introduced by Brito, Rosso & Soto [103] and also studied by Brey et al. [105,106] in the context of granular gas. In this case, particles undergo dissipative collisions and an energy influx, is necessary to maintain a steady state. In experiments of agitated granular particle this energy can be provided mechanically: the particles are distributed between two horizontal plates, and the bottom plate by a permanent shaking transfers the kinetic energy to the particles by collisions with the plate [107][108][109]. In this quasi-two dimensional geometry, the vertical vibrations are much faster that the time scale of the horizontal dynamic system and the energy provided vertically is transmitted horizontally to the system through collisions. The BRS model differs from the SCM in the way the injection of energy compensating the energy lost by inelastic collisions is modelled. Instead of super-elastic collisions, a fixed velocity "kick" is given to the colliding particles.

Compared to the SCM of the previous section, this model differs in several respects: Firstly, the energy is not conserved (even on average), but the system can attain a steady state of fixed macroscopic energy. Secondly, the long tails of the stationary state observed in the SCM are expected to be suppressed, because, as we will see in the collision rules, inelastic collisions dominate for high relative velocities. Thirdly, the model has, like the ICM, a strictly deterministic dynamics. It thus allows us to probe whether the explicit stochasticity of the SCM is really an essential feature for its qualitative behaviour, and notably the apparent existence of a universal non-equilibrium stationary state.

We first define the model and discuss some of its particular features. We then perform a numerical study, first considering the model without gravity (which has not been studied, in one dimension, in references [103,105]). We then study the full model with (1d) gravity. The results we report are a preliminary analysis of the model and we will discuss briefly, in our conclusions, how this analysis should be completed in future work.

BRS collisions

Let 1 and 2 be two particles of the system (such that x 1 < x 2 ) which collide with relative velocity v

12 = v 1 v 2 > 0.
The post-collisional velocities v ⇤ 1 and v ⇤ 2 are given by the following momentum conserving rules:

v ⇤ 1 =v 2 + v 12 (4.60) v ⇤ 2 =v 1 v 12 + with = 1 c R 2 , c R
the coefficient of restitution, and > 0 the constant velocity "kick". This is the one dimensional version of the rules in [103]. We note that in [105] the same collision rules are considered, except that the sign of is negative. In this case a particle gains velocity in the direction of the other colliding particle after the collision, which tends to increase the probability of re-collision phenomenon and leads rapidly to numerical difficulties. Thus we limit ourselves the original collision rules in [103]. We thus might expect that the system may relax, starting for a generic initial condition with no net velocity drift, to a uniform spatial distribution and a velocity distribution of the form:

g st (v) = 1 2 ( (v v st ) + (v + v st )). (4.66)
However it is clear that the velocity distribution will not evolve more generically if there are only particles with the two velocities ± 2 , irrespectively of the spatial distribution.

In such states and specifically for the velocity distribution given by Eq. (4.66) we have:

K st = 1 2 mh⌃ N i=1 v 2 i i = N 2 8 2 .
(4.67)

Thus the kinetic temperature of such a stationary state (or possible oscillating state in the configuration space with the same stationary velocity distribution) is expected to be proportional to 2 2 . It is interesting to compare this with the predicted kinetic energy derived for the 2d case. In [103] the stationary temperature is derived using hydrodynamic approached, and in [105] in any dimension by a kinetic approach with the use of modified collisional operator. For a Gaussian distribution, this temperature is:

T G s = mc 2 R ⇡ 4(1 c 2 R ) 2 " 1 + s 1 + 4(1 c 2 R ) ⇡c 2 R # 2 .
(4.68)

The later results does not a priori apply to the one dimensional case of the numerical simulation of the next section. Indeed, in [105] we remark that the coefficients of the Sonine expansion used to derived the results diverge for d = 1, indicating that the Gaussian approximation is very bad in one dimension. However for c R ! 1 the latter expression give the same scaling ( T G s ⇠

2 ) in the parameter which characterize the perturbation represented by the BRS collisions.

Numerical results for BRS model without gravity

Method and units

We present now the results of our numerical simulations of this model. The simulation method only requires a simple modification compared to the previous model, adapting the collision rule according to Eq. (4.60). The particles are placed again in a box of size L = N with elastic walls. We restrict the investigation to the quasi-elastic limit ⌧ 0 which is the relevant limit in the gravitational case, in which we wish to treat BRS collisions as a perturbation. We start from a spatially homogeneous distribution, and the velocities are distributed following a Gaussian. We work again in units with the particles masses m = 1, their number density n 0 = 1; the variance of the initial Gaussian velocity distribution v 0 = 1 and thus ⌧ granu = 1 n 1 0 v 0 = 1. 120 All simulations start form the same initial distribution which is homogeneous spatially and a Gaussian in velocity space, and all simulations have (by choice of units) initially the same kinetic energy K 0 . The yaxis is normalised to K st (Eq. (4.67)), the kinetic energy of the predicted stationary state. We observe that all the simulations converge very well at longer times to K = K st , indicating that the system indeed appears to go to this stationary state. The time axes is normalised by because in this unit, for all cases, the relaxation appear to happen on the same time scale. As expected starting from K K st i.e. in a state with h|v 12 |i 2 2 , we observe a rapid decay of the energy: indeed in this case the model should be well approximated at early time by the ICM; while for K ⌧ K st initially collisions predominantly inject energy since h|v 12 |i ⌧ 2 2 . Fig. 4.12 shows, for a single realisation with K 0 K st , the temporal evolution of the velocity distribution, with velocity now normalised for convenience in units of 2 = 1. The initial "hot" distribution cools down to a width of order unity and then rapidly reaches a clearly bimodal distribution. Peaks at v = ± 2 emerge and persist over the rest of the simulation time. There appears to be some evidence that these peaks narrow as time goes on. Longer simulation would be required to establish this more clearly. Indeed it would be interesting to determine the time scale for such relaxation if it does occur. Fig. 4.13 shows, for a single realisation of the same case, the temporal evolution of the velocity distribution, with a greater resolution in time, than in the previous figure. Although, to a first approximation, the distribution appears homogeneous, Figure 4.12: Distribution of the velocity v /2 for one realisation of 1d granular system with BRS collisions. The snapshots are taken at time t = 0, 5, 10, 15, ...100⌧ granu . The simulation is for N = 2048 particles with = 0.01 and = 0.01 it is evident that there are coherent time dependent fluctuations at large spatial scales. These appear to be the one dimensional equivalent of the hydrodynamic modes observed in the 2d version of this model in [103,105].

Results

In summary, in the 1d BRS model without gravity we observe that arbitrary initial states are driven efficiently to a state which has a very stable energy ⇠ relaxes to the same value, 0.15. This indicates that the BRS collisions drive the system toward a non equilibrium (non Gaussian) stationary state independent of the initial condition.

The next six figures show the spatial and velocity distributions and the phase space plots, for two chosen simulation from Fig. 4.14: the case = 0.01, = 0.0001 for "hot" and "cold" water-bag initial conditions (with E p hv 2 i(t), so that we can compare easily the shape of the distribution as the system shrinks. In Fig. 4.18 and Fig. 4.19 the same quantities are plotted for an initial water-bag initial condition with R 0 = 0.01. These figures illustrate and confirm the conclusion suggested by the plot of the energy and of 11 (Fig. 4.14). For both initial conditions, the two distributions follow a similar evolution: a first regime of rapid virialization of the system on a time of order 10 20 ⌧ dyn , and a second regime of evolution through "scaling" QSS (similar to ICM) up to 100 ⌧ dyn , which correspond more or less to the time scale of energy relaxation (Fig 4.14). Then both velocity and position distributions relax to an almost exactly constant distribution with very sharp edges. Like the initial water-bag distribution, the distributions subsequently attained are apparently homogeneous between the two [ x max , x max ] for position and [ v max , v max ] for velocity and equal zero outside these region (the value of x max and v max depends on / ).

The dynamical scenario with two regimes is confirmed: (1) an early time violent relaxation corresponding to the relaxation of the virial ratio (at t = 10 for the "cold" case one can see the noted "core-halo" state), (2) an evolution similar to ICM through "scaling QSS" . However these phase space plots reveal another feature of the evolution which could not be seen in the previous plots which appeared to indicate a stationary state is reached at t ⇠ 150⌧ dyn . In Fig. 4.20 the position in phase space of the 20 realisations of 128 particles are plotted on the top of each other for these two simulations at time t = 0, 10, ..240 ⌧ dyn . In the figure we see clearly that there is a further evolution (3) of the N body system in which the phase space distribution develops progressively an ordered aspect, with particles grouped on well separated curves.

This clearly reflects that the BRS collision which continue to produce a reduction of the phase space measure, making the system organise into a very ordered state. In this state the particles lie, in the phase space, on well define closed curves. For simulation of 128 particles, a central particle and five such orbits can be identified. (The five ellipses correspond in fact to 10 + 1 little bumps in the final velocity and spatial distribution in the previous figures). This "organisation" seems to emerge at the edge of the system and evolve progressively toward the center. This final state is to a very good approximation spatially homogeneous (Fig. 4.18 and 4.18). Then from the Poisson equation, one can deduce that the single particle potential energy is gn 0 x 2 . Thus, in this approximation, the one particle total energy of is e = 1/2 m v 2 + gn 0 x 2 , and at fixed particle energy the particles move on ellipses in the phase space. Then the five ellipses plus the central particle must correspond to six values of individual energy taken by the particles.

In Fig. 4.21 is shown the distribution of the individual energies e in units of E st the stationary total energy of the system. The individual energies of the particle are distributed in 6 clearly separated narrow regions of the energy space. This confirms the above analysis. In the left panel the 128 particles distribution at time 999⌧ dyn . In the right panel the distribution for a 256 particles simulation (t = 9980⌧ dyn ). We have plotted,on the top of each others, the position of the particles in the 1d phase space for 20 realisations (in black) and also a single realisation (red stars) in both case.

We have investigated briefly the dependence on particle number of the evolution towards this ordered state. Fig. 4.22 shows, at time 999⌧ dyn for 128 particles (left panel) and at t = 9980⌧ dyn for 256 particles (right panel), the phase space positions of the 20 particles of the 20 realisations as black points, while red stars represent a single realisation. We see that for one realisation (in a given orbit), particles are clearly distributed with a definite phase difference which was not evident on the previous plot. Second, the stationary phase space distribution for N = 256 particles (right panel) presents 8 ellipses (and no particle in the center) which are thus more densely packed than in the 128 particle case. Therefore, this attractor state (which does not depend on initial conditions) depends evidently on the number of particle. Thus unlike the first two phases of relaxation, we clearly cannot expect to describe using the mean field description of this ordered state and the relaxation to it. (We believe the latter should describe however very well the first two stages.) Further, the time scale of the relaxation to this state increases with the number of particles (the time of relaxation changes approximatively by an order of magnitude when we increase the number of particle by 2). Indeed, looking carefully the 256 particle case one sees small "defects" in the order of the one particle distribution, and also that the most central circles are not very well defined, indicating that the state is still "organising" at the time shown.

Finally, we have "turned-off" the perturbation once the system is in this state (at t = 1000⌧ dyn ) for a 128 particle simulation. We observed the system remains in this particular ordered state and found that it is stable at least up to 1000⌧ dyn . Thus, the state appears to be a stationary solution of the (Hamiltonian) N -body dynamics with elastic collisions.

Conclusion on BRS collisions.

The local and internal perturbation represented by BRS collisions appear to lead to a universal non-equilibrium stationary state. As in the SCM this states is very far from the thermal equilibrium of the 1d isolated self-gravitating system, but is very different to those observed in the SCM. For the granular system (without gravity), this state presents a bimodal velocity distribution that we were able to predict, and the spatial distribution showing hydrodynamic modes. The characterisation of the modes and a broader study of all the limits (in particular the limit c R ! 1) of this system would be an interesting development to complete the study. For the selfgravitating case, we find that the system reaches a steady state, again independent of the initial conditions. Once the system reaches the energy balance, particles organise precisely to a very stable, and ordered state. In this state particles lie on well defined, apparently elliptic, orbits in the phase space corresponding to different particle energy level. These results should be completed by further numerical and analytical study in order to establish that the two first stages of evolution can be described by a large N mean field limit. For the third phase, the relaxation to the final microscopically ordered phase occurs on a time scale dependent on N and also the properties of this "ordered" state could perhaps be understood analytically, maybe through some minimization procedures.

Conclusion

In this chapter, we have investigated the effects on the dynamics of long-range interacting systems of a certain class of "local internal" perturbations through the study of two one dimensional toy models. More specifically in both case we consider two perturbations of the dynamics inspired by granular studies, considering momentum conserving but non-Hamiltonian collisions.

We considered firstly the SCM, a non deterministic extension of the ICM, in which each collision is randomly inelastic or super-elastic. The collision rules are built such that on average, over the realisations of the stochastic process, the total energy remains fixed. In a second model, the BRS, the system is perturbed with inelastic collisions and driven by a constant energy input at each collision. The collisions, deterministic in this case, bring the system into a steady state where the gain and the loss of energy is balanced.

Both internal perturbation act microscopically and locally but completely modify the global organisation of the system: they push the long-range system far from its initial QSS. In that sense the QSSs are not robust to such perturbations. However, the evolution, in both cases, happens through a succession of virialised states, stationary solutions of the Vlasov equation. For both perturbations, the collisions, described by a non Hamiltonian term local in space and a functional of the velocity distribution added to the Vlasov equation, drive the system to a non-equilibrium stationary state (NESS). The stationary distribution does not depend on the initial conditions but depend strongly on the microscopic details of the perturbation. In both case we have checked that the NESS is a QSS i.e. it is a stationary solution of the unperturbed dynamics.

It is interesting to compare the latter results with the existing one in the literature. In [95] the stochastic force added to the long-range system which acts in a non-local manner on the system. Indeed the three particles of which the velocities are permuted, are chosen from the whole system. This process increases the mixing in velocity space, in a way which tends to destroy its space dependence and it drives the system to the thermal equilibrium. For the self-gravitating case, it has been verified that such a non-local perturbation indeed makes the system relax efficiently to the thermal equilibrium3 . In the models presented above (BRS,ICM,SCM), the perturbation, as we have underlined, has the property of being local, i.e the change in the velocity distribution depends only of the region of the system where the collision happens. It may be that such property is important for the stochastic perturbation to drive the system toward a NESS rather than thermal equilibrium.

Our study is complementary to that of reference [52]. As discussed in chapter 1, the authors add an external stochastic force (noise) acting also globally on the system. We can consider collisions as internal stochastic forces simply acting at the collision point, and what differs from [52], is that the force is coupled to the state of the system: this is appropriate to model the effects, for example, of others short-range interaction at play in the system. However, we found a similar feature of the dynamics as in the study [52] which is that the system indeed generically relaxes to a "stochastic attractor" independent of the initial conditions.

The latter results and the present study lead us to postulate that the relaxation toward thermal equilibrium is probably the exception for such perturb long-range systems and attained in very special circumstances. Conversely, QSSs with their non-trivial organisation in phase space appear to be a robust feature even when these systems are perturbed.

Figure 1 . 2 :

 12 Figure 1.2: Example of a concave entropy-energy curve (left panel) and its Legendre-Fenchel transform (right panel): the free-energy temperature curve.

Figure 1 . 3 :

 13 Figure 1.3: Sketch of a non concave entropy S verus the energy E (left panel), its Legendre transform S ⇤ ( ) = F ( )(middle panel), and the Legendre-Fenchel transform of the free-energy (full line, right panel) F ⇤ (E) = S ⇤⇤ (E) 6 = S(E).

CHAPTER 1 .

 1 INTRODUCTION TO LONG-RANGE INTERACTING SYSTEMS

CHAPTER 1 .

 1 INTRODUCTION TO LONG-RANGE INTERACTING SYSTEMS

CHAPTER 2 .

 2 FINITE N CORRECTIONS TO VLASOV DYNAMICS

Figure 2 . 2 :

 22 Figure 2.2: The variance of the force on particles in a phase space cell due to particles in the corresponding stripe, plotted as a function of ", for the different values of .The crosses indicate the points obtained from the numerical simulation described in the text and dotted lines are the corresponding analytical results, Eq. (2.63). We use units in which x = 1 and g = 1.

Figure 2 . 3 :

 23 Figure 2.3: The variance of the force on particles in a phase space cell due to particles in the corresponding stripe, versus ", for = 1.25 (upper panel) and = 1.75 (lower panel). The full expression is compare (see text) with the leading "-divergence of the 3 and 2 body contribution Eq. (2.64). The units are such that x = 1 and g = 1.

(3. 43 )

 43 Using the third relation, coming from the two last terms of the right hand side of Eq.(3.38) we have that:v2 xn = cste . (3.44)

  ) must be constant. The time derivative of Eq. (3.44) implies also the first relation of Eq.(3.43).

Figure 3 . 1 :

 31 Figure 3.1: Two examples of rectangular water-bag distributions, left panel for R 0 = 0.01 and right panel for R 0 = 1 both with 1024 particles. The units are such that the mean density n 0 = 1 and ⌧ dyn = 1.

Figure 3 . 2 :

 32 Figure 3.2: Phase space realisation for 1024 particles of the micro-canonical equilibrium.

Figure 3 . 4 :

 34 Figure 3.4: Phase space position 100 realisations of a 1024 particles system are plotted on top of each other for the times t/⌧ dyn = 0, 3.5, 10.5, 21.2, 31.8, 70, for the rectangular water-bag initial condition with R 0 = 0.01

Figure 3 . 5 :

 35 Figure 3.5: Evolution of the virial ratio of an isolated SGS with rectangular waterbag initial conditions (RWB) with initial virial ratio R 0 = 1(upper panels) and R 0 = 0.01 (lower panels) for N = 512, 1024, 2048, 4096 averaged over 100 realisations.

  of

Figure 3 . 7 :

 37 Figure 3.7: Spatial (right panels) and velocity (left panels) distribution function at time t = 0, 20, 40, 60, 80, 100 ⌧ dyn for a Fermi degenerate (upper panels) and thermal equilibrium (lower panel) initial conditions

Figure 3 . 8 :

 38 Figure 3.8: Spatial (right panels) and velocity (left panels) distribution function of the QSS reached after violent relaxation for an initial rectangular water-bag distribution with R 0 = 0.01 (upper panels) and R 0 = 1 (lower panels) at time t = 35, 70, 106 ⌧ dyn

Figure 3 . 10 :

 310 Figure 3.10: Normalised total energy per particle for VDM with ⌘⌧ dyn = 0.012, 0.037. In the upper panel are also plotted 3K and 2U as a function of time for ⌘⌧ dyn = 0.037. The lower panel displays the logarithm of the total energy for the two cases versus time. The different parameters are averaged over 50 realisations of the initial conditions.

Figure 3 . 15 :

 315 Figure 3.15: Normalised total energy, 11 and virial ratio versus t/⌧ dyn for the ICM with N = 0.01, for rectangular water-bag initial conditions with R 0 = 1, 0.1, 0.01 for 1024 particles (100 realisations). The inset shows a blow-up of the case R 0 = 1.

Figure 3 . 17

 317 Figure 3.17: Top: 11 versus t/⌧ mf for = 0.01, 0.005, 0.001 and R 0 = 0.01 (3 upper curves), R 0 = 0.1 (3 middle curves) and R 0 = 1 (3 lower curves). Bottom: 11 versus t/⌧ mf in a semi-log plot for R 0 = 0.01 and = 0.001 (red), 0.005 (green), 0.01 (blue).

  Figure 3.17: Top: 11 versus t/⌧ mf for = 0.01, 0.005, 0.001 and R 0 = 0.01 (3 upper curves), R 0 = 0.1 (3 middle curves) and R 0 = 1 (3 lower curves). Bottom: 11 versus t/⌧ mf in a semi-log plot for R 0 = 0.01 and = 0.001 (red), 0.005 (green), 0.01 (blue).

Figure 3

 3 Figure 3.18:11 as a function of t/⌧ dyn for ICM with = 0.01. In this simulation the dissipation has been stopped at time t = 100⌧ dyn . The inset shows the total energy versus time

Figure 4 . 1 :

 41 Figure 4.1: Evolution of the energy per particle average over 100 realisations, it estimated variance hK 2 i c , and the kurtosis 2 , SCM without gravity with = 0.01, for N = 1024, 8196 and with an initial homogeneuous state with Gaussian velocity distribution.

Fig. 4 .

 4 Fig. 4.2 shows the evolution of the kurtosis 2 , on a longer time scale than for the previous figure, for different values of versus the rescaled time t 2 /⌧ granu .The superposition of the curves indicate, that the time evolution of the system is rescaled by 2 in agreement with Eq. (4.37). At early times, the linear behaviour of Eq. (4.42) is verified, while at longer times 2 apparently reaches a independent "plateau" plus residual large fluctuations. The left hand side of Eq. (4.37) being positive, one would expect, a monotonic growth of the kurtosis. Fig.4.3 shows the evolution of the velocity distribution of the granular gas for

Fig. 4 .

 4 Fig. 4.2 shows the evolution of the kurtosis 2 , on a longer time scale than for the previous figure, for different values of versus the rescaled time t 2 /⌧ granu .The superposition of the curves indicate, that the time evolution of the system is rescaled by 2 in agreement with Eq. (4.37). At early times, the linear behaviour of Eq. (4.42) is verified, while at longer times 2 apparently reaches a independent "plateau" plus residual large fluctuations. The left hand side of Eq. (4.37) being positive, one would expect, a monotonic growth of the kurtosis. Fig.4.3 shows the evolution of the velocity distribution of the granular gas for

Figure 4 . 2 :

 42 Figure 4.2: Simulation of SCM without gravity with 2048 particles with varying: Evolution of the 2 averaged over 100 realisations versus 2 t/⌧ granu for = 0.01, 0.005, 0.0025.

Figure 4 . 3 :

 43 Figure 4.3: Plot (top) and log-log plot (bottom) of the velocity distribution for the SCM without gravity with = 0.01 at times (blue to red): t/⌧ granu = 0, 10, 20, .., 70

Figure 4 . 4 :

 44 Figure 4.4: Two top panels: Normalised total energy versus t/⌧ dyn for 5 different realisations (in color), the black line represents the mean over the 300 realisations of the process, for N = 128 (Top) and N = 512 (Middle). In the bottom panel the time evolution of the standard error for N = 128, 256, 512 particles ( SCM = 0.03).

Figure 4 . 5 :

 45 Figure 4.5: Top panels: Virial ratio as a function of time of a single realisation (in cyan), and for its mean over the 300 realisations (in black), for N = 512. Bottom panel: the evolution of the standard deviation over realisations for N = 128, 256, 512 particles ( SCM = 0.03).

Fig. 4 .

 4 Fig.4.5 shows, for the same model as in the previous figure, the evolution of the virial ratio as a function of time. The upper panel shows its value for an individual realisation (in color) and the averaged over 300 realisations for N = 512. The lower panel shows the standard deviation of R, for N = 128 to 512. Clearly the fluctuations of the virial ratio are below 1 (corresponding to a system which is virialised) and these fluctuations decrease in value and in amplitude as N increases. Thus as N ! +1, in the quasi-elastic limit, the system remains always exactly virialised and this is a strong indication that it evolves through a family of QSS. Fig.4.6 shows the evolution of the separability ratio 11 . The upper panel shows it for two realisations and for the average over 300 realisations.11 evolves on the same time scale as the variance of R, and reaches (as the variance of R) a well defined averaged "plateau" at longer times. The fluctuation from realisation to realisation of 11 are very large and their standard deviation is shown in the lower panel for N = 128 to 512. At earlier times it appears that the standard deviation decreases as N increases but at longer times it appears to be independent of N . This means that the system reaches a stationary virialised state but it is not clear that it converges to a well defined single value of 11 as N ! +1. This is further shown in Fig.4.7 which shows for the same time axes, the evolution of both 11 (top) and 2 (bottom) for different N . In the lower panel is also plotted the theoretical prediction from the

Figure 4 . 6 :

 46 Figure 4.6: Top panels: evolution of the separability ratio 11 for 2 realisations (in colors) and for the mean curve over the 300 realisations (in black), for N = 512. Bottom panel: Its standard deviation over realisations for N = 128, 256, 512 particle ( SCM = 0.03).

Figure 4 . 7 :

 47 Figure 4.7:11 (upper panel) and 2 (bottom panel) versus t/⌧ dyn of SCM with gravity for N = 128, 256, 512 particles with SCM = 0.03, starting from thermal equilibrium initial conditions. The black line on the bottom panel is the theoretical prediction of the initial slope starting from 2 = 3 at the initial time.

Figure 4 . 8 :Figure 4 . 9 :

 4849 Figure 4.8: Evolution of positions (right panels) and velocities distributions (left panels) for the SCM with gravity with SCM = 0.1 starting from an initial thermal distribution for 1024 particles and 100 realisations. The upper panels show log-log plots of the distributions.

Figure 4 .

 4 Figure 4.10: Parameter 11 versus t/⌧ dyn for simulations of 128 particles averaged over 225 realisations of the process. The three upper curves are simulations of SCM starting for different indicated initial conditions. The lower curve is the 11 of a simulation of a pure self-gravitating system relaxing to equilibrium on this time scale.

= v 2 . Thus we infer v 1 = ± 2 v 2 = ⌥ 2 .

 21222 velocities, v and v + , and further since the collision with the walls are elastic we must have v 1

Figure 4 . 11 :

 411 Figure 4.11: Log-log plot of the kinetic energy divided by K st (Eq. (4.67)) as a function of t/⌧ granu for a 2048 particles for different indicated value of and

2 2

 2 and with an almost stationary bimodal velocity distribution (probably converging to the exactly bimodal distribution of Eq. (4.66) ) and with persistent hydrodynamic modes about an average (in time) uniform spatial distribution.

Figure 4 . 13 :

 413 Figure 4.13: Spatial distribution for one realisation of 1d granular system with BRS collisions. The snapshots are taken at time t = 0, 1, 2, ...100⌧ granu . The simulation is for N = 2048 particles with = 0.01 and = 0.01 in a box L = N . The position are plotted here in units of L/2

Figure 4 . 14 :

 414 Figure 4.14: Simulation of a self-gravitating system with BRS collisions and with water-bag initial condition with R 0 = 1 and R 0 = 0.01. Averaged over 20 realisations, the total energy, R, 11 and 2 of a simulations with = 0.01 and = 0.0001 are plotted versus t/⌧ dyn for 128 particles

Figure 4 . 15 :

 415 Figure 4.15: E/( /4 ) versus t/⌧ dyn for one realisation of a simulation a 1-d selfgravitating with BRS collisions. The different curves are of several value of and and different initial conditions.

Figure 4 . 16 :

 416 Figure 4.16: Evolution of the spatial distribution versus x/ x (t) where x (t) is the estimated variance of the distribution at time t for 20 realisations of 128 particles with BRS collisions with = 0.01, = 0.0001 starting from a rectangular water-bag initial condition with R 0 = 1

Figure 4 . 17 :

 417 Figure 4.17: Evolution of the velocity distribution versusv/ v (t) ( v (t) = hv 2 i(t))for 20 realisations of 128 particles with BRS collisions with = 0.01, = 0.0001 starting from a rectangular water-bag initial condition with R 0 = 1

Figure 4 . 18 :

 418 Figure 4.18: Time evolution of the normalised spatial distribution for 20 realisations of 128 particles with BRS collisions with = 0.01, = 0.0001 starting from a rectangular water-bag initial condition with R 0 = 0.01

Figure 4 . 19 :Fig. 4 .

 4194 Figure 4.19: Time evolution of the normalised velocity distribution for 20 realisations of 128 particles with BRS collision with = 0.01, = 0.0001 starting from a rectangular water-bag initial condition with R 0 = 0.01

Figure 4 . 20 :

 420 Figure 4.20: Evolution of the particle position in phase space, for a simulation with 128 particles starting from rectangular water-bag with R 0 = 0.01 (the 25 top panels) and R 0 = 1 (bottom panels). On the top of each other the position of the particles of the 20 realisations of the self-gravitating -BRS dynamic are plotted at time t = 0, 10.., 240 ⌧ dyn . The parameters for BRS collisions are = 0.0001 and = 0.01.

Figure 4 . 21 :

 421 Figure 4.21: Distribution of the individual energy e in units of E st the stationary energy (value of the simulation) for a system of 128 self-gravitating particles with BRS collisions. The distribution in an average over 20 realisations and over a time window of with 10⌧ dyn (sampled at intervals of ⌧ dyn ).

Figure 4 .

 4 Figure 4.22: In the left panel the 128 particles distribution at time 999⌧ dyn . In the right panel the distribution for a 256 particles simulation (t = 9980⌧ dyn ). We have plotted,on the top of each others, the position of the particles in the 1d phase space for 20 realisations (in black) and also a single realisation (red stars) in both case.

  and express it using the time variable t e = t/t

	0 , we obtain:
	d 2 u i dt 2 e	+	4 3t e	du i dt e	=	3t 2 e 2	u i	(3.81)
	Making the change of time variable t ⌧ L =	r	3 2	ln(t e )	(3.82)

  To compute the collision time of two particles labelled by i + 1 and i initially at time t 0 e , we have to find the time t col e > t 0 e at which the solution of the equation satisfy u i+1 (t col e ) u i (t col e ) = n 1 L . We express this condition using the variable Y =

	5 x ⇣ ✓ t e t 0 e ◆ 2/3 + 2 5 ✓ t e t 0 e ◆ 1 # + t 0 e v i (t 0 e ) 3 5 " ✓ t e t 0 e ◆ 2/3 ✓ t e t 0 e ◆ 1 # (3.87) t 0 e e t col	⌘	1/3

i (t e ) = s i + u i (t e ) v i (t e ) = du i dt e with t 0 e , this initial time (or last collision time) in unit of t 0 66 CHAPTER 3. LONG-RANGE SYSTEMS WITH WEAK DISSIPATION Next collision time.

or "strong" long-range system in opposition to the "weak" long-range interaction to distinguish with the case ↵ > d ( > 0 )[2]

We note that the Kac prescription (or other equivalent) to make the energy extensive is important to define the free energy F = E T S. Indeed as the entropy S / N , one needs E / N so that both term defining F have the same scaling with N .

Figures taken form the slides of Hugo Touchette: http://www.maths.qmul.ac.uk/~ht/ archive/talks/safrica2011ensequivgen1_2.pdf

Here we assumed the system size is fixed, then this scaling is proportional to V as considered above.

This statement "seems",a priori, paradoxical because the macroscopic evolution described is irreversible, and one can easily check that the Vlasov equation is time reversible. The Vlasov

N . But dE dT is just, for a system in thermal equilibrium, its specific heat, and here we see it can be negative as mentioned

In a uniform Poisson process with mean density n 0 , the PDF P (N, V ) of the number of particlesN in a volume V is P (N, V ) = (n 0 V ) N e n0V/N ! (see, e.g.,[73])

Or, alternatively, if one adopts the terminology advocated by[2], in which the thermodynamic distinction is made between "strongly long-range" and "weakly long-range", our classification could be described as a distinction between "dynamically strongly long-range" and "dynamically weakly long-range" interactions.

A simple example is the case of one dimensional gravity, i.e., = 0 in d = 1. In this case Poissonian force fluctuations indeed diverge with system size[74] precisely as derived here. On the other hand, because the force is independent of distance, the force fluctuation in a cell coming from particles in other cells vanishes if the number of particles in these cells does not fluctuate.

For a system with elastic collisions, as noted in chapter 1, particles just exchange their velocities during collisions. One can then imagine a equivalent system of "ghost" particles, passing through one another. In this point of view, the force acting a ghost particle in discontinuous at each crossing and during ⌧ mf a particle undergoes a number of collision / N .

Originally, in the case of a infinite

3d system, the analogous force was regularised by subtracting the contribution due to the mean density (the background ) in a ad hoc manner (Jean's swindle). The derivation of the force from the expression(3.73) shows that it is a mathematically well defined regularization.

(t) = h|x(t)v(t)|i h|x(t)|ih|v(t)|i 1 (3.113)

is comparable to the one starting from cold RWB initial conditions. It is a striking feature of our results that the system, when its shows deviation from scaling QSS, does not show any tendency to evolve toward the thermal equilibrium (or a scaling version of this state), but on the contrary evolves away from this state towards more correlated states. To better understand this feature is one of the motivations of the next chapter.

In the view where the phase space of the non-Hamiltonian system is a manifold with a constant probability measure as proposed in[104], the manifold is locally dilated or contracted in a random manner at each collision.

N.B.: The definition of in this chapter is different the one in the previous chapter, where we used the notation q for this expression.

K st ).

Internship report of T. Epalle, 2013, under supervision of M. Joyce and D. Benhaiem

11 versus t/⌧ dyn for rectangular water-bag R 0 = 0, 0.01, 0.1, 1, thermal equilibrium and fermi initial state. The SGS simulations are performed with 1024 particle and the value of the two parameter averaged for over 50 realisations for each initial conditions of the separability ratio 11 in Fig. 3.6. For the thermal initial condition 11 remains equal to zero as expected, and for the Fermi degenerate state 11 is also constant and close to the one reached by the RWB R 0 = 1 initial condition. I In Fig. 3.7, the velocity and spatial distributions are displayed for times t = 0, 20, 40, 60, 80, 100 ⌧ dyn for Fermi degenerate and thermal equilibrium initial conditions. The superposition of the curves establishes that these two states are (quasi) stationary solution of the dynamics. In Fig. 3.8 shows at different times the distributions of a system starting from cold (R 0 = 0.01) and hot (R 0 = 1) water-bag distribution, once the system has relaxed to a QSS. The hot initial water-bag relaxes to distributions very close to the Fermi degenerate state, but differing only by additional small tails (halos) in the rectangular water-bag case. 11 and virial ratio as a function of time for VDM with ⌘⌧ dyn = 0.037, for RWB initial condition with R 0 = 1 (green) and R 0 = 0.01 (blue). The inset is a blow up of the virial oscillations for the case R 0 = 1 in the regime of violent relaxation. The different parameters are averaged over 50 realisations of the initial conditions.

ICM: SGS with dissipation through inelastic collisions

According to the discussion after Eq. (3.66), this is a limit case. Indeed this model is a particular case of the ICM class, with = c = 1, corresponding to the ratio between the two time scales ⌧ diss and ⌧ dyn which remains constant in the scaling QSS. Thus as we have discussed, for this case the scaling solution may be expected to represent well the behaviour of the system for a transient period. The four first panels of Fig. 3.13 show the early time evolution of the system (t = 0, 7, 14, 28 ⌧ dyn ). It appears to be just like the period of violent relaxation for a Hamiltonian long-range interacting system, and as we see below, during this period, Vlasov stable initial conditions. To explore further the domain of validity of the scaling-QSS as a function of the initial condition taken in this limit case of the theory, we run simulations with thermal equilibrium and Fermi degenerate initial condidtions (see sec. 4.4). Contrary to rectangular water-bag, these states are, by construction, stationary solutions of the Vlasov equation (Vlasov stable).

Rectangular water-bag initial conditions

In Fig. 3.19, the top panels show the evolution of the spatial and velocity distributions starting from a Fermi degenerate initial condition as a function of the rescaled variables x/x(t) and v/v(t). We observe that for this case also the scaling solution appears to be a accurate description of the dynamics. In the bottom panels, the same quantities are plotted for a simulation starting from the thermal equilibrium. In this case also, we observe clear deviation from the scaling dynamics at later times. For the two cases the evolution of the energy of the system is the same as for RWB initial condition and we do not display it here.

In Fig. 3.20 we show the evolution of the virial ratio and the separability ratio 11 of a simulation starting from the Fermi degenerate state. For comparison, we also display these quantities starting from an initial rectangular water-bag with R 0 = 1. In the inset we can verify that the system is indeed already virialised when starting from this initial condition. Further for this case 11 remains constant throughout the evolution, confirming that, for this initial condition also, the systems dynamics is well described by the scaling QSS. Large N scaling of Following the same reasoning as for the treatment of the ICM in the previous chapter, we consider how the parameter must be scaled in the N ! +1 mean field limit in which the Vlasov equation describes the evolution of the phase space density due to the long-range interaction. For the ICM, we saw that taking the limit 1 c R ⌧ 1, the Kramers Moyal expansion gives a term proportional to

= ICM . From Eq. (4.16) we see that, for the SCM, we need to take 2 N ⇠ C st as N ! +1 so that the collision term does not diverge in the large N limit. We thus introduce a rescaled parameter

In our numerical simulation with gravity we will keep this parameter fixed as we vary N .

Moments of the collision operator

We consider the following moments of the full collisional integral

We first consider the first term on the right hand side of the equation, and using the definition of the collisional operator we obtain:

) are the pre-collisional velocities. Splitting the last integral, changing the velocity coordinates of the first term with the post-collisional

) and regrouping terms we obtain:

Now we define

The rest of the integrand, |v

where S n is the symmetric part of n under the exchange v 1 $ v 2 defined by,

It is easy to compute the following expressions:

Inserting these results in (4.24), we can deduce the following expression for the moments of the collisional term:

which is an expression of the conservation of momentum in collisions, further, 1)

The equivalent expressions for the operator J ˜ are easily deduced. Then using the relation (1 ) = ˜ (1 ˜ ), we find:

which is simply the expression of energy conservation. For the fourth moment we have:

Temporal evolution of moments

We have now the expression of different moments of the collision operator. We now integrate the Vlasov equation over v n to get the time evolution of these moments.

The first term on left hand side of Eq. (4.12) gives:

In the case without gravity the mean field term in a Vlasov equation vanishes and the term proportional to @ x f vanishes if we assume the system to be spatially uniform. In the gravitating case we can, following the same reasoning as in the previous chapter, assume that these two terms also do not contribute in the limit of weak perturbation PERTURBATIONS

Numerical results for the granular gas

In this section we present out numerical results of SCM for the case where the particles do not interact other than at collisions. This is the 1d case of a model studied in [97]. We consider this model first to check that we can reproduce previous results, and also to probe the correctness of the kinetic theory in this case.

Method of simulation and units

We simulate the SCM without gravity with an event-driven dynamics in precisely the way described in the previous chapter for the ICM. The system is thus confined in a box of size L 0 with elastic boundaries conditions. We choose to work in units of m = 1, L 0 = N and ⌧ granu = 1

(equals thus unity) is the initial kinetic temperature.

In the code, we use the following definition of the perturbation parameter:

(i.e. we do not make use of SCM introduced above for the moment.)

We also measure the kurtosis of the velocity distribution 2 in Eq. (4.41). We shall also estimate the variance of the (kinetic) energy calculated as follows:

with K r is the energy of the realisation r, R the realisations number and hKi = 1 R P R r=1 K r the mean energy, in order to compare with of the theoretical variance, Eq. (4.48).

Numerical results

In the first panel of Fig. 4.1, we display the time evolution of the energy (averaged over realisations) of the SCM without gravity with = 0.01 for N = 1024, 8196. The mean energy is conserved to a good approximation up to a time which increases as N grows. The second panel of Fig. 4.1, shows the evolution of the estimated variance hK 2 i c .We clearly see that as N increases, the variance decreases and that the linear behaviour predicted in section 2.7 is well verified at early time. For longer times, larger fluctuations for both quantities appears. Fluctuation comes from the large deviation of a single realisation and one would need an average over more realisations to recover the average behaviour. In the third panel of Fig. 4.1, the kurtosis 2 of the velocity distribution is plotted as a function of t/⌧ granu . The straight line (in black) is a plot of the equation (4.42) with the parameter value of the simulation. The kurtosis starts from 3 (corresponding to the initial Gaussian distribution) as expected.

The fluctuations around the averaged curve becomes significant after t ' 1000⌧ granu . As expected the Eq. (4.42) describes well the deviation from the Gaussian distribution up to a finite time. These results indicate that the large N limit in which the kinetic theory was derived is an accurate description of the system at early time, for a large enough number of realisations.

Numerical results for 1d self-gravitating gas

Method of simulation and units

In order to simulate the SCM with gravity, we modify the code implemented for ICM by adding super-elastic collisions of Eq. (4.1). At each collision, with a random generator, the code chooses, with equal probability, if the collision is super-elastic or inelastic. We recall that this code used is valid to simulate a finite self-gravitating system only if particles remains confined in a small region compared to the size of the box. We choose units with m = 1 and the initial spatial density n 0 = 1. It is appropriate to adopt the same time unit as the SGS model of the last chapter: ⌧ dyn = 2 p gn 0 . As discussed in section 2.3 we characterize the collisions by SCM :

To probe the macroscopic evolution, we measure the parameter 11 and the virial ratio defined in chapter 3, and make also use of the kurtosis of the velocity distribution 2 in Eq. (4.41). For the each quantity Q (Q = E, 11 , R, ..) we define the estimated standard deviation:

with Q r the value of the quantity for the realisation r, R the number of numerical realisations and hQi the mean quantity. We use this standard error to check whether the mean over the finite number of realisations can be taken to be representative of the true mean.

Validity of the mean field description

The two upper panels of Fig. 4.4 show, for the SCM with SCM = 0.03, the temporal evolution of the energy, both for 5 individual realisations (in colour) and its average over 300 realisations (in black), for N = 128 (top panel) and 512 (lower panel). We see that for the former case the average is close to constant only at earlier times but not at longer times; while for the case N = 512, the conservation is within about 10% over the whole time range. The behaviour of the standard deviation, shown in the third panel, estimated from 300 realisations for each of N = 128, 256, 512, shows that the behaviour of fluctuation of the estimated mean is indeed just a reflection of the root of the variance which grows with time more or less monotonically. These behaviours are in line with our expectation following the discussion is section 2.7, and we note that the standard deviation of the energy appears (qualitatively) to grow in proportion to p t, which is consistent with the expression of the variance Eq. (4.48). We see clearly that, at fixed t, as N ! +1 we do indeed expect the energy to be conserved exactly as predicted.

Depending on the initial relative velocity of the two colliding particles, the absolute value of relative velocity between particles can increase or decrease. Indeed the kinetic energy change in a collision is:

The unique solution for colliding particles (with v 12 > 0) of the equation K BRS = 0 is:

Thus for 0 < v 12 < , K BRS > 0 the kinetic energy gain is positive, while for a relative velocity such that v 12 > the system loses energy. Note in particular that particles with velocity close to zero gain a velocity ± .

Stationarity

It is important to remark that the velocity kick does not depend on the precollisional velocity. Thus for a BRS collision the transformation of the phase space volume is just obtained from the inelastic part:

Each collision thus give a contraction of the occupied phase space region, and through the evolution the support of the probability distribution function is thus necessary reduced. In other words the relevant measure dv 1 ...dv N to describe the velocity space becomes arbitrary small with time, and selects more and more limited values of particle velocities. T his contracting measure is a property of this non-Hamiltonian system, just as it is of the ICM [104]. However while the ICM leads to a cooling which asymptotically brings the system to a zero energy state in which all particles are at rest, we expect the BRS model to lead to a stationary state at fixed non zero energy. (Indeed a stationary state with zero energy cannot exist because there is always then a net energy injection). Thus we conclude that if the BRS can reach a stationary state it must have a distribution with support of zero measure in the phase space of dimension 2dN .

From the collision rule (4.60), it is evident that the evolution may in fact lead to the "selection" of very particular values of the velocity and give a stationary distribution of zero measure: we have seen that if the relative velocity at collision v 12 has the value , the collision is elastic and thus the particles simply "exchange velocities" and the velocity distribution does not evolve i.e.

These considerations are valid for both the gravitational and granular case (without gravity). Now, in the case without gravity velocities change only at collision when they collide. Therefore if we have a distribution with v 12 = ± for all pairs of particles, the distribution is macroscopically stationary. This will be true if we have any two

Numerical results: BRS with gravity

Method and units

We now consider how the dynamics of a 1d self-gravitating system is modified by BRS collisions. To simulate this system we again use the code with a periodic box (see chapter 3), simply modifying appropriately the collision rules.

We work in units with particle masses m = 1, the initial density n 0 = L 0 N = 1 where L 0 is the initial region of the phase space (we also choose in this case L 0 ⌧ L, L the size of the periodic box, which is the requirement for the accuracy of the regularized gravity code) and time is measured again in units of

just as for SGS and ICM (chapter 3). As for the ICM, we expect ⌧ dyn to be the characteristic mean field time only at early time: as the system evolves away from its initial states its density can change greatly and thus the characteristic time of the mean field dynamics.

Preliminary results

In Fig. 4.14 is shown the average over 20 realisations of the energy, the virial ratio R, and the parameters 11 and 2 as a function of time for 128 particles with = 10 2 and = 10 4 . The blue curves are the results for "cold" water-bag initial distribution (R 0 = 0.01) and the green are for a "hot" water-bag initial condition with R 0 = 1 ( for details see chapter 3).

In the first panel, the temporal evolution of the total energy is shown in a semilog plot. For both simulations the energy relaxes toward the same stationary value for both initial conditions in approximately 150⌧ dyn . In the second panel we see that the system virialises in only ⇠ 10 20 ⌧ dyn and then remains very close to a virialised state.

Just as for the ICM, the inset is a blow up of the initial virial oscillations for the "hot" initial conditions. As might be anticipated, these oscillations have characteristic period decreasing in time. Indeed, as the energy decreases, the density n increases. and the period of the mean field oscillation 2 p gn(t) thus decreases.

The third panel and the fourth panel show the evolution, respectively, of the separability ratio 11 and the kurtosis of the velocity distribution 2 . Both simulations show, up to ⇠ 100⌧ dyn , a behaviour very close to that observed in the ICM for the same initial condition. After a first phase of (violent, in the case R 0 = 0.01) relaxation, lasting ⇠ 10 20 ⌧ dyn , both quantities settle to non-zero values, which are very close to a constant in the case R 0 = 1 but slowly increase in the case R 0 = 0.01. Indeed given that the starting energy E 0 is much greater than ⇠ 2 2 , the model is very well approximated initially by the ICM: then when the energy reaches a value close to 2 2 the behaviour of 11 and 2 rapidly changes as the system is rapidly (in ⇠ 10 20 ⌧ dyn ) apparently driven into a stationary state. After this second relaxation all parameters are stable and their oscillation are remarkably small (even for only one realisation), indicating that the system has reached a well defined stationary state, independent of the initial condition.