Attractor non-equilibrium stationary states in perturbed long-range interacting systems

1.8 MB·

by  Michael Joyce

Attractor non-equilibrium stationary states in perturbed long-range interacting
systems
Michael Joyce1, Jules Morand1,2,3,4, and Pascal Viot4
1 Laboratoire de Physique Nucl´eaire et de Hautes E´nergies, UPMC IN2P3 CNRS UMR 7585,
Sorbonne Universit´es, 4, place Jussieu, 75252 Paris Cedex 05, France
2 National Institute for Theoretical Physics (NITheP),Stellenbosch 7600, South Africa
3 Institute of Theoretical Physics, Department of Physics,
Stellenbosch University, Stellenbosch 7600, South Africa and
4 Laboratoire de Physique Th´eorique de la Mati`ere Condens´ee, UPMC,
CNRS UMR 7600, Sorbonne Universit´es, 4, place Jussieu, 75252 Paris Cedex 05, France
6
(Dated: January 15, 2016)
1
0 Isolated long-range interacting particle systems appear generically to relax to non-equilibrium
2
states (“quasi-stationary states” or QSS) which are stationary in the thermodynamic limit. A

n fundamentalopenquestionconcernsthe“robustness”ofthesestateswhenthesystemisnotisolated.
a In this paper we explore, using both analytical and numerical approaches to a paradigmatic one
J dimensional model, the effect of a simple class of perturbations. We call them “internal local

4 perturbations”inthattheparticleenergiesareperturbedatcollisionsinawaywhichdependsonly
1 on the local properties. Our central finding is that the effect of the perturbations is to drive all
the very different QSS we consider towards a unique QSS. The latter is thus independent of the

] initialconditionsofthesystem,butdeterminedinsteadbyboththelong-rangeforcesandthedetails
h
of the perturbations applied. Thus in the presence of such a perturbation the long-range system
c
evolves to a unique non-equilibrium stationary state, completely different to its state in absence of
e
the perturbation, and it remains in this state when the perturbation is removed. We argue that
m
this result may be generic for long-range interacting systems subject to perturbations which are
- dependent on the local properties (e.g. spatial density or velocity distribution) of the system itself.
t
a
t PACSnumbers: 05.20.-y,04.40.-b,05.90.+m
s
.
t
a
m I. INTRODUCTION libria of such systems, just as Maxwell Boltzmann (MB)
equilibriaareforshort-rangesystems. Theoreticallythey
-
d Systems of large numbers of interacting particles are are understood to be stationary solutions of the Vlasov
n subject in their physical analysis to a fundamental dis- equation which in principle describes these systems in
o tinction based on whether they are short-range or long- the relevant thermodynamic limit. Unlike MB equilib-
c
[ range, depending on the rapidity of the decay with sep- ria, they are infinitely numerous at given values of the
aration of the two body interaction potential. The dis- global conserved quantities, and the actual equilibrium
1 tinction in its canonical form arises from the presence or reached depends strongly on the initial condition of the
v absence of the property of additivity of the macroscopic system.
4
energy, which plays a fundamental role in equilibrium AbasicquestionwhicharisesaboutQSSinlong-range
4
6 statistical mechanics. While most of familiar laboratory interacting systems concerns the “robustness” of such
3 systemsstudiedinphysicsareshort-range—notablyany states. They are strictly defined only for isolated Hamil-
0 systemconstitutedofneutralatomsormolecules—there tonian systems and the question is whether they con-
1. arenumerousexamplesalsooflong-rangesystems, rang- tinue to exist when the system is not exactly isolated, or
0 ingfromself-gravitatingsystemsinastrophysicsandcos- exactly Hamiltonian, or both. For attempts to observe
6 mology,tovorticesinturbulentfluids,lasercooledatoms, these intriguing equilibria in laboratory systems, which
1 and even biological systems (for a review, see e.g.[1, 2]). are necessarily perturbed by and coupled to the external
:
v Study of various isolated long-range systems has shown world in some way, it is an essential to know whether
i that they evolve from generic initial conditions to mi- these states can be expected to survive. Studies of toy
X
croscopical non-Boltzmann equilibria, known as “quasi- models coupled to a thermal bath (see e.g. [10]) show,
r stationarystates”(QSS)becausetheyevolvetowardsthe unsurprisingly, that such a coupling sends the system to
a
system’struestatisticalequilibriumontime-scaleswhich itsthermalequilibrium,onatimescalewhichdependson
diverge with the number of particles (see e.g. [3–8]) The thecoupling. However, muchmoregenerally, ithasbeen
evolution to such states appears not to be characteristic suggested on the basis of study of the one dimensional
ofalllong-rangeinteractions,butonlyofthesub-classof HMF model (see e.g.[11, 12]) that QSS will disappear
these interactions for which the pair force (rather than in the presence of any generic stochastic perturbation to
pair potential) is non-integrable at large distances (see thedynamics. Astudyoftheeffectofexternalstochastic
[7, 9]). As these systems remain in the QSS indefinitely fields with spatial correlation applied to the same model
in the thermodynamic limit, these states can be consid- (see e.g. [13, 14].) shows however that interesting non-
ered to be the fundamental relevant macroscopic equi- equilibrium steady states can be obtained in this case.
2
Inthisarticleweexplorethequestionoftherobustness erature, and comment on the possible generality of the
of QSS in long-range systems to weak perturbations us- behavioursweobserveinabroadclassofperturbedlong-
ing a paradigmatic toy model of long-range interactions range systems.
—aonedimensionalself-gravitatingsystem—subjected
to a particular class of perturbations, which we refer to
as “local internal perturbations”: the perturbations to
II. MODELS
the purely self-gravitating dynamics occur when parti-
cles collide, and are described by simple collision rules
for the colliding particles, which may be stochastic or A. The sheet model
deterministic. The perturbations can thus be considered
to model physical effects which come into play at very We consider a system of identical particles of mass m
small scales e.g. due to very short scale forces and/or moving in a one dimensional space and interacting by a
internal degrees of freedom. Differently to the study of force independent of their separation, i.e. the force on a
[13, 14] there are therefore no external forces, and in- particle i due to a particle j is
deed we will build our collision rules so that they con-
serve total momentum. The choice of the two specific F =gm2sgn(x −x ), (1)
ij i j
two models we study is then guided by simplicity: once
momentum conservation is imposed, a non-trivial colli-
where g is the interaction strength. The model is known
sionruleinonedimensioncannotconserveenergy(asan
as the “sheet model” because these particles in one di-
elastic collision gives rise simply to exchange of particle
mensionareequivalenttoinfinite,infinitelythin,parallel
velocities). As we wish to focus here on the effect on
sheets moving in three dimensions interacting by New-
the effect of the perturbations on QSS, which have fixed
tonian gravity, in which case gm = 2πΣG, where G is
energy,wechoosetodefineperturbationswhichcancon-
Newton’s constant and Σ is the mass per unit surface
serve average energy and lead (in principle) to a steady
area of the sheets. This model dates back at least to
state. The two models we study are then simple choices
the early study of Camm[15] and has been studied quite
withthisproperty,correspondingtocollisionrulesdrawn
extensively by numerous authors since (see e.g. [16–18]
from the literature on granular gases arising from simple
and references therein).
considerationsofenergybalance. Intheonedimensional
ForafinitesystemofN particles,thetotalforceacting
self-gravitating model these perturbations given by non-
on the ith particle at any time is simply given as
trivial collision rules are both very simple to implement
numerically, and, aswewillsee, alsoadmitastraightfor-
F =gm2[N+−N−], (2)
ward theoretical description for the kinetic theory in an i i i
appropriate mean-field limit.
where N+ and N− denote the numbers of particles on
Our main result is that in both models we study we i i
the right and on the left of ith particle, respectively.
observe an evolution from the initial QSS, which arises
Thedynamicsofthissystemhasbeenextensivelystud-
fromagiveninitialcondition,throughafamilyofstates,
ied(seee.g. [16–18])andshowsthatthetimeevolutionof
untilatrulystationarystateisreached. Theevolutionis
the system displays long-lived QSS before reaching equi-
through a continuum of QSS and the final state is also a
librium. More precisely the relaxation time associated
QSS: if the perturbation is removed the system remains
with QSS is diverging with the system size τ ∼ AN,
in this state. Further the final stationary state in each
where the coefficient A depends strongly on the initial
case appears to be an attractor for the perturbed long-
states [18].
range system, i.e., starting from different initial condi-
We consider now two variants of this model, in which
tions which evolve to very different QSS the perturba-
the particle collision still conserves the momentum but
tion drives them all finally to the same non-equilibrium
not the total energy. On the other hand, as we wish the
state. This “universal” QSS is thus determined essen-
systemtobeabletoattainastationarystateatconstant
tially by the detailed nature of the small perturbation
energy, we constrain the exact collision rules to allow
and the long-range force itself.
this. Bothcollisionlawsaretakenfromsimplemodelsof
The paper is organized as follows. We first define the granular systems which have been studied in the litera-
two models we study of self-gravitating particles per- ture (see references below).
turbed in a specific manner when particles collide. In
the following section we describe analytical approaches
to these models which are valid in an appropriate mean-
field limit. These give rise to kinetic equations which B. Model A: collisions with random coefficients of
restitution
allow us to determine a well defined large N limit. They
also provide predictions for the early time evolution of
the system. In the next section we present results of nu- Let us denote v = v − v the relative velocity of
ij i j
merical studies of the two models. In our conclusions we particles i and j which undergo a collision with precol-
compare our results with other relevant works in the lit- lisional velocities v and v . We adopt the rule that the
i j
3
postcollisional velocities, v∗ and v∗, are given by and therefore the energy change is
i j
1−c
vi∗ =vj + 2 vij (3) δK =m(cid:20)c2B −1v2 +c ∆|v |+∆2(cid:21) (9)
1−c 4 ij B ij
v∗ =v − v ,
j i 2 ij
where the coefficient of restitution c is a non-negative The term in ∆ in the collision rule thus leads to an
randomvariable. Equivalentlyitcorrespondstomomen- energy injection, which can be smaller or larger than the
tum conservation combined with the rule energy loss due to the inelastic term: more precisely, the
collision leads to an energy gain if |v | < v , and an
ij 0
v∗ =−cv . (4) energy loss if |v |>v , where
ij ij ij 0
Granular models of this kind, incorporating a random
2∆
coefficientofrestitution,havebeenintroducedby [19]in v = (10)
ordertoincludetheeffectofenergyinjectioninavibrated 0 1−cB
two-dimensional granular gas.
In each collision the change of the kinetic energy is isthevalueoftherelativevelocityforwhichthecollision
given by is elastic.
Thiscollisionruleistheonedimensionalversionofthat
m c2−1
δK = [(v∗)2−(v )2]= (v )2 (5) introducedintwodimensionsby[20]inaphenomenolog-
4 ij ij 4 ij
ical model of quasi-two dimensional experiments of ag-
i.e. the collision is inelastic if 0 ≤ c < 1 and super- itated granular particles: the particles are confined be-
elastic if c > 1. In order that the applied perturbation tween two horizontal plates, and the vibrating bottom
mayadmitstationarystates,wechoosecfromabimodal plate transfers the kinetic energy to the particles by col-
probability distribution in which c takes two values, c lisions [21–23]. In this quasi-two dimensional geometry,
A
and c˜ , with equal weight, with 0<c <1 and theperiodoftheverticalvibrationsismuchshorterthan
A A
the typical time scale of the horizontal dynamics. The
(cid:113)
c˜ = 2−c2 (6) collision rule, Eq. (7), then represents a time coarse-
A A
grained description of the energy transfer of particle-
The latter relation imposes that, for a collision at the bottomplatecollisionstohorizontalparticle-particlecol-
same initial relative velocity, the energy lost with c=c lisions.
A
is the same as the energy gained when c = c˜A. Thus, In this paper we have chosen this collision rule simply
in the ensemble of realizations of the stochastic pertur- because it provides a simple way, quite different to that
bations, the average energy is constant, with the energy in the first model, to obtain a non-trivial two body col-
loss of the inelastic collisions (0 < cA < 1) balanced by lision rule which can be expected to lead to a stationary
the energy gain in superelastic collisions (c˜A > 1). We state. More specifically if the particle velocities at colli-
expect that in this case the system may be able to reach sions are assumed to be uncorrelated, the kinetic energy
astationarystatewithconstantenergy(modulofiniteN of the system gives a direct measure of the typical rela-
fluctuations). tive velocity of colliding particles: (cid:104)(v −v )2(cid:105) = 2(cid:104)v2(cid:105).
i j i
The kinetic energy would then be expected to be driven
towards a value of order Nv2, as above this energy scale
0
C. Model B: inelastic collisions with energy energy will be dissipated while below it energy will be
injection injected. Indeed in the case in which gravity is turned
off,andtheparticlesareenclosedinaboxwithreflecting
Inthismodeltheself-gravitatingparticlesundergocol- walls,ifallparticleshavevelocity±v /2allcollisionsare
0
lisions specified by the following rule: elastic and the velocity distribution does not evolve at
all.
1−c
vi∗ =vj + 2 Bvij −(cid:15)ij∆ (7) Model B is in fact microscopically deterministic, while
Model A is explicitly stochastic. One other notable dif-
1−c
vj∗ =vi− 2 Bvij +(cid:15)ij∆ ferenceisthatthephasespacevolumeoccupiedbyparti-
cles involved in a collision strictly contracts in Model B,
where the constant c has a fixed value in the range 0< while it can contract or increase in Model A depending
B
c < 1 (i.e. as for an inelastic collision), ∆ is a positive on whether the collision is inelastic or elastic. Indeed for
B
constant(withdimensionsofvelocity)and(cid:15) =sgn(v ). acollisionwithcoefficientofrestitutioncineithermodel
ij ij
The collision manifestly still conserves total momentum, we have dv∗dv∗ =cdv dv which is always a contraction
i j i j
and corresponds to inModelB.Thispropertyleadstodistinctivefeaturesof
the long time behaviour of the models which we observe
v∗ =−sgn(v )[c |v |+2∆]. (8)
ij ij B ij below.
4
III. KINETIC THEORY for binary collisions, which is the sum of two contribu-
tions:
A. Mean field limit without collisions J [f,f](x,v,t)=G (x,v,t)−L(x,v,t). (15)
q q
where G is the gain term corresponding to collisions
For the purely self-gravitating model the dynamics in q
where a particle has a post-collisional velocity equal to
the appropriate large N mean field limit is described by
v,
theVlasovequation[9,24,25]. Thislimitisobtainedby
taking N → ∞ at fixed values of the total system mass N (cid:90) (cid:90)
G (x,v,t)= dv(cid:48)dv(cid:48)(cid:48)|v(cid:48)−v(cid:48)(cid:48)|f(x,v(cid:48),t)f(x,v(cid:48)(cid:48),t)
M and energy E. Denoting the mass density in phase q M
space, f(x,v,t), the Vlasov equation reads δ(v−qv(cid:48)−(1−q)v(cid:48)(cid:48)), (16)
∂ f(x,v,t)+J [f]=0 (11) and L(x,v,t) is the loss term corresponding to collisions
t V
where a particle with a velocity v undergoes a collision
where J [f] the Vlasov operator is at time t,
V
N (cid:90)
J [f]=v∂ f(x,v,t)+a¯(x,t)∂ f(x,v,t) (12) L(x,v,t)= f(x,v,t) dv(cid:48)|v(cid:48)−v|f(x,v(cid:48),t). (17)
V x v
M
where a¯(x) is the mean-field acceleration given by Notethatthelosstermdoesnotdependexplicitlyonthe
a¯(x,t)=g(cid:82) sgn(x−x(cid:48))f(x(cid:48),v(cid:48),t)dx(cid:48)dv(cid:48). The mass den- coefficient of restitution, and indeed we can write
sity obeys the normalization condition
J [f,f](x,v,t)=G (x,v,t)−G (x,v,t) (18)
q q 0
(cid:90) (cid:90)
Let us introduce a series expansion of the δ function
dxdvf(x,v,t)=M. (13)
in terms of the parameter q
where M is the total mass of the system. δ(v−qv(cid:48)−(1−q)v(cid:48)(cid:48))= (cid:88) (q(v(cid:48)(cid:48)−v(cid:48)))nδ(n)(v−v(cid:48)(cid:48))
QSSareinterpretedasstationarysolutionsofEq.(11). n!
n≥0
There are an infinite number of such solutions, includ-
(19)
ing as a particular case the statistical equilibrium of this where δ(n) denotes the nth derivative of the δ function.
model (see below).
The Boltzmann operator is then expressed as
N (cid:88)(cid:90)(cid:90) qn((v(cid:48)(cid:48)−v(cid:48)))n
J [f,f](x,v,t)= dv(cid:48)dv(cid:48)(cid:48)|v(cid:48)−v(cid:48)(cid:48)|
B. Model A q M n!
n≥1
δ(n)(v−v(cid:48)(cid:48))f(x,v(cid:48),t)f(x,v(cid:48)(cid:48),t) (20)
1. Collision operator in Boltzmann approximation
To determine whether the parameters characterizing
The Vlasov equation can be derived starting from the the collisions can be rescaled with N so that the colli-
BBGKY hierarchy and making the approximation that sion term remains well defined (and non-trivial) in the
the two point correlations can be neglected. The colli- mean-field (Vlasov) limit, we consider the limit q → 0
sionsinourmodelcanbetreatedinthesameapproxima- of the model i.e. the quasi-elastic limit. Physically this
tion, and are then described by a canonical Boltzmann is clearly the relevant limit: to obtain an N indepen-
operator. We thus expect our model in the mean field dentevolutioninpresenceofthecollisionsontime-scales
limit to be described by a kinetic equation characterizing the mean-field dynamics (e.g. the time
a particle typically takes to cross the system) one must
∂ f(x,v,t)+v.∂ f(x,v,t) clearly “compensate” the effect of the divergent growth
t x
(cid:88) of the number of collisions with N by making the effect
+∂ a¯[f](x,t).∂ f(x,v,t)= P(q)J [f,f](x,v,t).
x v q of each collision arbitrarily weak.
q
(14)
2. Expansion of kinetic equation
where, for convenience, we introduce the parameter q =
1−c to characterize binary collisions with a coefficient
2 Inserting Eq. (20) in the right hand side of Eq. (14)
of restitution equal to c, and J [f,f](x,v,t) is a collision
q and integrating by parts term by term we obtain (fol-
operatoraccountingforsuchcollisionswhichareassumed
lowing [26–28]) the collision operator given as a series of
tooccurindependentlywithaprobabilityP(q). Initially
differential operators:
wewillleaveP(q)undeterminedandthenreplaceitwith
thespecificbimodalformforModelAattheappropriate N (cid:88)
J [f]= [(cid:104)qn(cid:105)∂n(f(x,v,t)a [f](x,v,t))
point below. A M v n
n≥1
Assumingtheparticlestobepointlike,thecollisionop-
(21)
erator is a homogeneous Boltzmann operator accounting
5
with obtained for any variant of this model in which P(q) is
such that energy is conserved on average.
(cid:90) (v(cid:48)−v)n
a [f](x,v,t)= |v(cid:48)−v| f(x,v(cid:48),t)dv(cid:48) (22) We note further that in this derivation we have as-
n n! sumed implicitly that f(v) has the convergence prop-
erties required for the validity of the Taylor expansion,
wherethebrackets(cid:104)..(cid:105)indicateanaverageovertheprob-
which requires clearly sufficiently rapid decay of f(v) at
ability distribution P(q). Note that a [f](x,v,t) has a
1 largev toensurethefinitenessofthecoefficients. Indeed
simplephysicalmeaningasanaverageeffectiveforcedue
we see that while the finiteness of Eq.(16) requires only
to the collisions [29, 30].
that f(v) be integrable at large |v| (i.e f(v) ∼ 1/|v|α
Let us consider now the specific P(q) of model A:
with α > 2), the definiteness of the expression Eq.(27)
requires α > 4. As we will discuss below the latter as-
1 1
P(q)= 2δ(q−qA)+ 2δ(q−q˜A) (23) sumption turns out to break down at longer times in the
model.
where qA = 1−2cA and, from the average energy conserv-
ing condition Eq.(6),
3. Evolution of moments of velocity distribution
(cid:112)
1− 1+4(1−q )q
A A
q˜ = (24)
A 2 We now discuss some properties of the collision oper-
ator by considering the evolution of the moments of the
Expanding in qA (as qA → 0 in the quasi-elastic limit) velocity distribution. Multiplying both sides of the ki-
we have netic equation by vn, and integrating over v, we obtain
q˜ =−q +2q2 −4q3 +O(q4) (25) d (cid:90) (cid:90)
A A A A A [ρ(x,t)vn(x,t)]+ dvvnJ [f]= dvvnJ [f] (28)
dt V A
and thus, to leading order in powers of q we obtain
A
(cid:82)
where ρ(x,t)= dvf(x,v,t) is the spatial mass density,
(cid:104)qn(cid:105)=nqn+1 nodd and vn(x,t) is the n-th moment of the velocity distribu-
A
(cid:104)qn(cid:105)=qn neven. tion at x, i.e.,
A
(cid:90) f(x,v,t)
Defining now vn(x,t)= dvvnp (v,t) where p (v,t)=
x x ρ(x,t)
√
√ (1−c ) N (29)
γ ≡q N = A (26)
A A 2 It is straightforward to show, either directly from the
exact collision operator, or for each of the two terms in
we have, at leading order in 1/N, Eq. (28), that
(cid:90)
γ2
(cid:104)q(cid:105)=(cid:104)q2(cid:105)= A dvJA[f]=0
N
(cid:90)
3γ4
(cid:104)q3(cid:105)=3(cid:104)q4(cid:105)= A , dvvJA[f]=0
N2
while, for n > 4, (cid:104)qn(cid:105) decreases with N more rapidly which express, respectively, the conservation of particle
than 1/N2. Thus taking the mean-field limit N → ∞ number and conservation of momentum in the collisions.
Indeed it is straightforward to show that the left-hand
at constant γ , the full expansion of the collision term
A
side of Eq. (28) corresponds for n = 0 to the continuity
reducestothesumofthetwofirstderivativesoff(x,v,t):
equation, and n=1 to the Euler equation.
γ2 Forthecasen=2, integrationbypartsusingEq.(27)
JA[f]= MA [∂v(f(x,v,t)a1[f](x,v,t)) gives
+∂v2(f(x,v,t)a2[f](x,v,t))(cid:3) (27) (cid:90) γ2 (cid:90)
dvv2J [f]=2 A dv(va [f]+a [f])f(x,v,t) (30)
A M 1 2
Note that the non-linear structure of the integral colli-
sion operator Eq. (21) is conserved because the velocity-
from which it follows using Eqs. (22) that
dependent functions a (v) are functionals of f(x,v,t).
n
We note that the crucial relation leading to the result (cid:90)
(27) for this model is (cid:104)q(cid:105) = (cid:104)q2(cid:105), which is simply the dvv2JA[f]=0, (31)
condition of average energy conservation. Indeed from
Eq. (5)itfollowsthattheenergychangeinacollisionat whichexpressestheconservationofthekineticenergyby
anygivenrelativevelocityisproportionaltoq−q2. Thus the collisions. Thus the local pressure ρv2 can change
the same mean-field limit for the kinetic theory will be only due to the mean field gravitational force (through
6
the Vlasov flow term J [f]). Note that while the van- kinetic energy) can be neglected, we have
V
ishing of the zero and first moments hold for any P(q),
it can be verified from Eq. (21) that the second moment dβ2(t) (cid:39) γA2 (cid:82) dx(cid:82) dv(cid:82) dv(cid:48)|v−v(cid:48)|5f(x,v,t)f(x,v(cid:48),t)
vanishes only if (cid:104)(q−q2)(cid:105) = 0 which is, as noted above, dt (cid:0)(cid:82) dx(cid:82) dvv2f(x,v,t)(cid:1)2
just the condition of average energy conservation.
(36)
For any n≥2 it is straightforward to show that
We will test this prediction below for the case where
(cid:90) nγ2 (cid:90) (cid:90) the initial state is the statistical equilibrium of the sys-
dvvnJ [f]= A dv dv(cid:48)|v−v(cid:48)|(v−v(cid:48)) tem.
A 4M
×(cid:2)(n−3)(vn−1−v(cid:48)n−1)−(n−1)vv(cid:48)(vn−3−v(cid:48)n−3)(cid:3)
×f(x,v,t)f(x,v(cid:48),t) (32) C. Model B
from which we recover the previous result for n=2, and Following exactly the same approach we write the ki-
further find that the first non-zero moment is for n=4, neticequation,intheBoltzmannandmean-fieldapprox-
and it has the simple expression imations, for this model as
(cid:90) γ2 (cid:90) (cid:90)
dvv4J [f]= A dv dv(cid:48)|v−v(cid:48)|5f(x,v,t)f(x,v(cid:48),t). ∂tf(x,v,t)+JV[f]=JB[f,f](x,v,t). (37)
A M
(33) where J [f,f](x,v,t) is the collision operator with the
B
The fact that right hand of this expression is strictly same structure as Eq. (15), but with the gain operator
positive has an important consequence: if this kinetic now given by
equation is valid, the system cannot reach a stationary
state. Or, conversely, if the system reaches a stationary G (x,v,t)= N (cid:90) (cid:90) dv(cid:48)dv(cid:48)(cid:48)|v(cid:48)−v(cid:48)(cid:48)|f(x,v(cid:48),t)f(x,v(cid:48)(cid:48),t)
state,itmustbesuchthattheassumptionsnecessaryfor q M
the derivation of the kinetic equation break down. As δ(v−qv(cid:48)−(1−q)v(cid:48)(cid:48)−sgn(v(cid:48)−v(cid:48)(cid:48))∆) (38)
noted above we will see that our numerical study shows
that the system generically evolves to such a regime. In where q = 1−cB is a fixed positive parameter less than
2
fact, we will see that when the system reaches a sta- unity.
tionary state it is characterized by a non-Gaussian ve- Followingthesameargumentsasabove,itisclearthat
locity distribution with tails decaying as a power law. toobtainacollisionoperatorJ [f,f](x,v,t)whichisin-
B
Indeed the estimated exponent of the velocity distribu- dependentofN inthemean-fieldlimit,wemustconsider
tion is such that the 4th moment is not defined and thus a quasi-elastic limit, with q →0 and ∆→0 as N →∞.
the above equation is not applicable in this state. Further if the evolution to a stationary state is to be de-
Ontheotherhandifthesystemispreparedinaninitial scribed in such a limit, the energy of this state, which
state which is a QSS, and which does satisfy the condi- wehaveinferredmustbe∼Mv2,mustbeextensive(like
0
tions necessary for the validity of the derivation leading theenergyinthemeanfieldlimit)andthereforev must
0
to Eq. (27), we can use Eq. (33) to infer non-trivial betakenindependentofN. Now,since∆=qv ,holding
0
information about the temporal evolution at sufficiently v fixed and taking q → 0 indeed defines a quasi-elastic
0
short times. Indeed in this case limit.
d(cid:16) (cid:17) γ2(cid:90) (cid:90) Proceeding as in the previous case, we perform again
ρ(x,t)v4(x,t) = A dv dv(cid:48)|v−v(cid:48)|5f(x,v,t)f(x,v(cid:48),t) an expansion of the Boltzmann operator in powers of q
dt M
about q =0. This gives
(34)
N (cid:88)
In practice we measure the integrated quantity, i.e. the J [f,f]= qn∂n(a (x,v,t)f(x,v,t)) (39)
B M v n
rescaledfourthmomentofvelocity(thekurtosis),defined n≥1
by
where
M(cid:82) dx(cid:82) dvv4f(x,v,t)
β (t)= (35) 1(cid:90)
2 (cid:0)(cid:82) dx(cid:82) dvv2f(x,v,t)(cid:1)2 an[f](x,v,t)=n! |v(cid:48)−v|(v(cid:48)−v−sgn(v(cid:48)−v)v0)nf(x,v(cid:48),t)dv(cid:48)
(40)
i.e. the fourth rescaled moment of the full velocity dis-
tribution P(v,t) = 1 (cid:82) dxf(x,v). For a Gaussian dis- Defining now
M
tterimbupteiroant,utrheeofvathlueesoyfstβe2misacnodnestqaunatl,tiond3e.pTenhdeatnimteofetvhoe- γB ≡qN = (1−2cB)N (41)
lution β (t) characterizes deviation of the distribution
2
from a Gaussian distribution. and taking N → ∞ at constant γ , we obtain a finite
B
Integrating Eq. (34) over x and, assuming that the limit for the collision operator which corresponds to the
time evolution of the second moment (i.e. of the total meanfieldlimit. Inthiscaseonlytheleadinglinearterm
7
of the expansion contributes, and the effect of collisions [−V /2,V /2]. For the case of gravity only, which
0 0
corresponds to the presence of an effective velocity de- has no characteristic length scale, this is a one pa-
pendent force per unit mass a (x,v): rameter family of initial conditions which may be
1
convenientlycharacterizedfullybytheinitialvirial
γ
J [f]= B [∂ (f(x,v,t)a [f](x,v,t))] . (42) ratio R (and the particle number N), where the
B M v 1 0
virial ratio R (at any time) is defined as
We note that the diffusive term which was non-zero in
the mean-field limit of Model A thus vanishes for Model 2K
R= (44)
B. U
As for the model A, we can calculate the velocity mo-
where K is the kinetic energy and U the potential
ments of the collision operator J [f]. The first two mo-
B
energy. A virialized system this has R=1.
ments again vanish as a consequence of conservation of
particle number and momentum, while
• Thermal equilibrium: the statistical equilibrium
(cid:90) γ (cid:90) of the purely self-gravitating system in the micro-
dvv2J [f]=2 B dvva [f]f(x,v,t) (43)
B M 1 canonicaland canonicalensemble hasbeenderived
for any finite N by [34] and its mean-field limit
DifferentlytomodelAthisisnotzero,ingeneral: indeed
(derived earlier by [15]) is
the model does not necessarily conserve energy on aver-
atoger.eaOchnathsetaottihonerarhyanstdatweeinexwpheiccthtthheescyosltliesmiontoopbeeraatboler f(x,v)= √M e−σv22 sech2(cid:16)x(cid:17) (45)
2 πσΛ Λ
is zero, and in this case Eq. (43) will vanish.
withσ2 = 4E,Λ= 4E andE isthetotalenergy.
3M 3gM2
IV. NUMERICAL RESULTS
3. Units
A. Simulation method and units
Forourstudytheonlydimensionalparametersofrele-
1. Code
vancearethetime,and,inmodelB,thevelocity(because
of the parameter ∆). A natural choice of units for both
The molecular dynamics of a one dimensional self-
are those characteristic of the mean field dynamics. For
gravitating model is conveniently simulated using an
the time unit we choose
event-drivenalgorithmasbetweenparticlecollisionstra-
jectories can be calculated explicitly. Such an algorithm 1
is exact up to the machine rounding error in computing τdyn = √gρ , (46)
0
the solutions of quadratic equations giving the collision
times(see[16,18,31]andreferencestherein). Furtherthe where ρ is the initial mass density of the system, and
0
algorithm may be sped up using a “heap structure” [32] for the velocity
andbyupdatingpositionsateachsteponlytheparticles
involvedineachcollision. Itisstraightforwardtomodify (cid:114)
2E
this algorithm to implement the simple collision rules of vdyn = 3M0 (47)
ourtwomodelsinsteadofelasticcollisions(equivalentto
particlecrossings). Weuseamodifiedversionofthecode whereE istheinitialenergy. Withthisdefinitionv2 is
described in [28] (and greater detail in [31]) [33]. Model 0 dyn
thevelocitydispersionofavirializedsystemwithenergy
Aischaracterizedbythechoiceoftheparameterc <1,
A E .
0
and each collision is then chosen with probability 0.5 to
Previous studies (see e.g. [18]) of evolution from this
be inelastic (with c = c ) or superelastic (with c = c˜ ).
A A firstclassofinitialconditionsfortheself-gravitatingsys-
For model B is characterized fully by the values of c
A temshowthatthesystemevolves,inatimeoforder10−
and ∆.
100τ ,toQSSofwhichthepropertiesdependstrongly
dyn
on R . At longer times, of order (102−103)Nτ , the
0 dyn
different QSS all relax to thermal equilibrium [4].
2. Initial conditions
Forbothmodelswestudyevolutionstartingfromthree
4. Additional macroscopic observables
kinds of initial conditions:
• “Rectangular waterbag”: particles are randomly Tomonitorinasimplewaytheevolutionoftheglobal
distributed with uniform probability in a rect- properties of the system, we measure in addition to the
angular region of phase space, [−L /2,L /2]× energy and the virial ratio, that of a global parameter
0 0
8
which is a simple measure of the “phase space entangle-
1.4
ment” of the state of the system: N=512
1.3 N=256
N=128
(cid:104)|xv|(cid:105)
φ11 = (cid:104)|x|(cid:105)(cid:104)|v|(cid:105) −1 (48) 1.2
E0
/1.1
fi
As shown in [18] the only stationary solution of the E
›
Vlasov equation which is a separable function of posi- 1.0
tion and velocity is that corresponding to thermal equi-
0.9
librium. Thus if φ is constant and non-zero this indi-
11
cates that the system is in a QSS distinct from thermal 0.80 2000 4000 6000 8000 10000 12000 14000 16000 18000
equilibrium, and its amplitude can be taken roughly as
a measure of “proximity” to the latter. We also monitor
1.4 N=512
the evolution of the kurtosis β as defined in Eq. (35). N=256
2 1.2 N=128
2E0
/1.0
B. Model A 2E0.fi8
›
−0.6
fi
AsseeninsectionIIIB,therelevantparametercharac- 2E0.4
›
terizing the perturbations due to collisions in the mean-
0.2
fieldlimitisγ definedinEq. (26). Ifthislimitdescribes
A
0.0
accurately the dynamics of the system, the term arising 0 2000 4000 6000 8000 10000 12000 14000 16000 18000
from the perturbations is proportional to γ2, and thus
A t/τdyn
the time scale on which they are expected to modify the
evolution of the system is ∼ τ /γ2. In order to study
dyn A
thedesiredrangeofweakperturbation,andassesstheva-
FIG. 1. Model A: (Top panel) Dimensionless energy (cid:104)E(cid:105)/E
0
lidity of the mean-field limit, we will thus consider small (E is the initial energy) and (bottom panel) dimensionless
0
values of γA and vary N keeping γA fixed. We report variance ((cid:104)E2(cid:105)−(cid:104)E(cid:105)2)/E02 versus dimensionless time t/τdyn
here results for γ = 0.03 and γ = 0.1, and for N in averagedover100realizationswithγ =0.03, ofrectangular
A A A
the range N =128 to N =1024. We average our results waterbag initial conditions for N =128,256,512 particles.
over a large number of realizations in each case.
we have, however, significant finite N fluctuations devel-
1. Evolution of energy and virial ratio oping in all cases at times a few times τdyn/γA2, so we
might anticipate that such effects may begin to play a
significant role on these time scales.
We have constructed this model so that the fluctua-
Fig. 2 shows the evolution of the mean (upper panel)
tions in energy of the system should become arbitrarily
and standard deviation (lower panel) of the virial ra-
small for sufficiently large N. Indeed in the mean field
tio in the same ensemble of simulations as in the pre-
limit we have derived above the energy is exactly con-
vious figure. Because of the equilibrium initial condi-
served, and this limit evidently thus does not describe
tions, the system remains always, as we would expect,
effects associated with the energy fluctuations at finite
very close to virialized, with only finite small N fluctua-
N. In our numerical study, at finite N, we therefore
tions which (lower panel) clearly decrease monotonically
need to check whether, on the time scale simulated, the
as N increases. We will see below that the non-trivial
energy fluctuations are indeed small. Fig. 1 shows the
time dependence of these fluctuations are a reflection of
evolution as a function of time of the mean energy, for
themacroscopicevolutionofthesystemonthesametime
a model with γ = 0.03, in an ensemble of realizations
A
scales.
starting from waterbag initial conditions with the differ-
ent indicated N, over a time scale roughly an order of
magnitude greater than τ /γ2. We see that, on these
dyn A
time scales, that in all cases, the ensemble averaged en- 2. Macroscopic evolution due to perturbation
ergy is indeed very close to constant, but (lower panel),
the variance of the normalized energy (i) grows almost Fig. 3 shows the time evolution of the entanglement
linearlyintime,asindicatedbythedashedstraightlines parameterφ andtherescaledkurtosisβ (bothdefined
11 2
and(ii)monotonicallydecreasesasN increases. Thus,as above), averaged over 100 realizations of equilibrium ini-
wewouldexpect,takingN sufficientlylargeatanygiven tial conditions with N =1024 particles. As expected, at
time, we can in principle converge to arbitrarily precise t = 0, φ = 0 and β = 3: the equilibrium state phase
11 2
conservation of the energy in any single realization. In spacedensityisaseparablefunctionofspaceandvelocity,
our numerical simulations at (relatively small) finite N, and the velocity distribution has a Maxwell-Boltzmann
9
1.10 4.0
3.5
1.05 3.0
2) 2.5
1
5
N=1.00 φ11 21..05
(
R 1.0
0.95
0.5
0.0
0.90 0.5
0 1000 2000 3000 4000 5000 0 200 400 600 800 1000 1200
0.35 70
N=256
0.30 N=512 60
0.25 50
2 fi
R0.20 40
›
− β2
0.fi15 30
2R
0.›10 20
0.05 10
0.00 0
0 1000 2000 3000 4000 5000 0 200 400 600 800 1000 1200
t/τ t/τ
dyn dyn
FIG. 3. Model A: Evolution of the entanglement parameter
FIG. 2. Model A: (top panel) virial ratio R as a function of φ (toppanel)andofthedimensionlesskurtosisβ averaged
11 2
time t/τdyn averaged over 100 realizations with equilibrium over 100 realizations with γA = 0.1 and equilibrium initial
initial conditions and for γA = 0.03 and N = 512. (Bottom conditions for N = 1024. The black curves correspond to a
panel) variance of the virial ratio for N = 512 and N = 256 modifiedsimulationinwhichtheperturbationisswitchedoff
with the same initial conditions. for t≥300τ .
dyn
shape. The red curve corresponds in each case to the
evolution with γA = 0.1, while the black curve is the perturbatively.
evolutionwhenthisperturbationis“switchedoff”atthe
time t=300τdyn, i.e., starting from this time the evolu- • the parameter φ11, which in absence of the per-
tion is that of the purely self-gravitating system. Note turbation would remain stable at its initial value
that in line with what would be expected from mean- φ = 0, evolves to a final value φ ≈ 2−3 de-
11 11
field theory, the characteristic time for the macroscopic pending on N, i.e., towards a state in which the
evolution is about ten times shorter than in the data in correlation between position and velocity is ever
the previous figures, for the case γA =0.03. We see that stronger. The perturbation clearly drives the sys-
the evolution induced by the perturbation is through a tem away from thermal equilibrium.
continuum of QSS, i.e., at all times the system remains
very close to a stationary and stable state of the Vlasov
equation. Given that the perturbation is weak — in the
3. Validity of mean-field kinetic theory
sense that it perturbs the system macroscopically on a
time scale long compared to τ — this is indeed what
dyn
Let us now consider the degree to which the evolu-
one would expect. However, what is not evident, and to
tion in our simulations of the system are described well
be underlined, is that
by the mean-field limit of the kinetic equations derived
• the effect of the small perturbation is not, at suf- above. The most basic prediction of this theory is that,
ficiently long times, perturbative: the system is when we adopt the associated scalings of the parameters
clearly progressively driven very far from its ini- with N, we observe an evolution which is independent
tial state. This contrasts strongly to the behaviour of N. The upper two panels of Fig. 4 show respectively
onewouldexpectforfamiliarshortrangesystemsin the evolution of φ and β , in each case averaged over
11 2
whichasmallperturbationwouldbeexpected(e.g. 300 realizations of the model with γ = 0.03 starting
A
applyinglinearresponsetheory)toevolvetoanout from equilibrium initial conditions, and for the different
of equilibrium state close to the initial one, in the particle numbers indicated: N = 128,256,512. In both
sense that all macroscopic quantites are changed cases we observe that at sufficiently early times there is
10
a non-trivial evolution of the system which is very well 3.5
128
superimposed for the different N, indicating the validity 3.0 256
of the mean-field theory. Further the evolution of β at 2.5 512
2
these times agrees well with that predictedbythe mean- φ112.0
1.5
field kinetic theory at early time, shown as a straight
1.0
line obtained using Eq. (36) with f(x,v) taken equal to 0.5
the initial thermal equilibrium phase space density Eq. 0.0
0 2000 4000 6000 8000 10000
(45). At longer times however we see that the evolution
changes: foreachN,theevolutionbreaksaway,atatime
40
scalewhichappearsroughlytoincreasewithN,fromthe
30
common behaviour, and shows on a similar time scale a
β2
tendency to reach a plateau, indicating in principle the 20
attainment of a stationary state. By measuring the ve-
10
locity and spatial distributions below we will verify that
0
this is indeed the case, and in so doing also find the ex- 0 2000 4000 6000 8000 10000
planation for the N-dependence and the noisiness of the
0.65
evolution of φ and β at longer times: the stationary
11 2
sthtaetseetpoarwthiciuchlarthmeascyrsotsecmopeivcovlvaersiaibsleinsbfaecctomoneeiflolrdewfihniechd /vdyn00..6505
fi
in the mean field limit. This is the case because these v||0›.50
states are characterised by velocity and spatial distribu- 0.45
tions which have slowly decaying power-law tails at long 0.40
0 2000 4000 6000 8000 10000
distances, for which both (cid:104)|x|(cid:105) and (cid:104)v4(cid:105) diverge. Their
values in a finite simulation are then regulated by the t/τdyn
cut-off due to the finite particle number, and are thus
highly fluctuating. Shown in the lower panel of Fig. 4 is
FIG. 4. Model A: Entanglement parameter φ , dimension-
theevolutionof(cid:104)|v|(cid:105),which,incontrast,isawelldefined 11
less kurtosis of the velocity distribution β and (cid:104)|v|(cid:105)/v ,
quantity in the final state. In this case we see that the 2 dyn
averagedover300realizations,asfunctionsoftimet/τ for
dyn
evolution for different N in the mean-field scaling agrees
γ =0.03 and N =128,256,512 (from top to bottom).
A
well, and fluctuates little, right up to the time at which
the stationary state is attained. The small deviation at
the latest times for N =128 is a result of the large fluc- 5
tuations of the total energy in this case (see Fig. 1 ). Eq.
4 R0=1
3 R0=0.01
φ11
2
4. Dependence on initial conditions
1
0
Fig. 5 shows the evolution of φ and (cid:104)|v|(cid:105)/v for 0 2000 4000 6000 8000 10000
11 dyn
γ = 0.03 and N = 512 and with three different ini- 0.70
A
tial conditions: thermal equilibrium, and two rectangu- 0.65
lar waterbags, R0 = 1 and R0 = 0.01. The number of vdyn0.60
realizations are 300,100,200 respectively. At long times, / 0fi.55
these quantities evolve towards the same mean value, in- v||0›.50
dependently of the initial conditions. Note that for φ 0.45
11
fluctuations increase with time and are associated with 0.400 2000 4000 6000 8000 10000
the long tails of position and velocity distributions.
t/τ
Inthepresenceoftheperturbation,thesystemgoesat dyn
longtimes to anonequilibriumstationarystate whichis
an attractor of the dynamics, i.e., the perturbed dynam-
FIG. 5. Model A: φ (top) and of (cid:104)|v|(cid:105)/v (bottom) as a
icsofthislong-rangesystemhasa“universal”stationary 11 dyn
function of time t/τ , with γ =0.03 and N =512 for dif-
state. dyn A
ferent initial conditions: thermal equilibrium (averaged over
300 realizations), rectangular waterbags with R = 1 (100
0
realizations), and R =0.01 (200 realizations).
0
5. “Universal” final state
Let us consider now the properties of this apparently phase space distribution is indeed stationary and the
stationary state, and check in particular whether the same in the different cases.