Phase transitions in self-gravitating systems and bacterial populations with a screened attractive potential

1.1 MB·English

by  P. H. Chavanis

Phase transitions in self-gravitating systems and bacterial populations
with a screened attractive potential
P.H. Chavanis and L. Delfini1
1
Laboratoire de Physique Th´eorique (IRSAMC), CNRS and UPS,
Universit´e de Toulouse, F-31062 Toulouse, France
WeconsiderasystemofparticlesinteractingviaascreenedNewtonianpotentialandstudyphase
transitions between homogeneous and inhomogeneous states in the microcanonical and canonical
0
ensembles. Like for other systems with long-range interactions, we obtain a great diversity of
1
microcanonical and canonical phase transitions depending on the dimension of space and on the
0
importanceofthescreeninglength. WealsoconsiderasystemofparticlesinNewtonianinteraction
2
in the presence of a “neutralizing background”. By a proper interpretation of the parameters, our
n studydescribes(i)self-gravitatingsystemsinacosmologicalsetting,and(ii)chemotaxisofbacterial
a
populations in theoriginal Keller-Segel model.
J

2 PACSnumbers: 47.10.A-,47.15.ki
1

] I. INTRODUCTION is modified so as to allow for the existence of spatially
h
homogeneous distributions at equilibrium. This will be
c
e calledthe modified Newtonian model. In that model, the
m Many biological species like bacteria, amoebae, en- source of the potential is the deviation between the ac-
- dothelial cells, or even ants interact through the phe- tual density ρ(r,t) and the average density ρ. This is
at nomenon of chemotaxis [1]. These organisms secrete a similar to the effect of a “neutralizing background” in
chemical substance (like a pheromone) that has an at-
t plasmaphysics [10]. This model canbe derivedfromthe
s tractive (or sometimes repulsive) action on the organ-
. Keller-Segel model in the limit of vanishing degradation
t ismsthemselves. Thisphenomenonisresponsibleforthe
a of the chemical [11]. It also appears in cosmology, due
self-organization and morphogenesis of many biological
m to the expansion of the universe, when we work in the
species. Ithasalsobeenproposedasaleadingmechanism
- comoving frame [12]. It is therefore interesting to con-
d for the formation of blood vessels during embriogenesis siderthis formofinteractionatagenerallevelandstudy
n [2]. Onatheoreticalpointofview,chemotaxiscanbede-
the corresponding phase transitions. We shall also com-
o scribedbytheKeller-Segel[3]modeloritsgeneralizations
pare them with the ones obtained within the ordinary
c [4]. The Keller-Segel model consists in a drift-diffusion
[ equationfortheevolutionofthedensityofbacteriaρ(r,t) Newtonian model [13–30] (see a review in [31]).

1 coupled to a reaction-diffusion equation for the evolu- The paper is organized as follows. In Sec. II, we
v tion of the secreted chemical c(r,t). In certain approxi- discuss several kinetic models taken from astrophysics,
2 mations, the reaction-diffusion equation is replaced by a plasma physics and biology for which our study applies.
4
Poissonequation. In that case, the Keller-Segel (KS) [3] We consider either isolatedsystemsdescribedby the mi-
9
1 model becomes isomorphic to the Smoluchowski-Poisson crocanonicalensemble(fixedenergyE)ordissipativesys-
. (SP)system[5]describingself-gravitatingBrownianpar- tems describedby the canonicalensemble (fixedtemper-
1 ticles(see,e.g.,[6]foradescriptionofthisanalogy). The atureT). We characterizetheir equilibriumstates in the
0
KS model and SP system have been studied thoroughly mean field approximation: in the microcanonical ensem-
0
1 in applied mathematics (see Refs. in [7]) and in theoret- ble(MCE),theymaximizetheentropyatfixedmassand
: ical physics (see Refs. in [5]). energy and in the canonical ensemble (CE) they min-
v
However, the original KS model [3] also allows for the imize the free energy at fixed mass. In Sec. III, we
i
X possibilitythatthechemicalsuffersadegradationprocess specifically consider the case of a Newtonian interaction
r which has the effect of reducing the range of the interac- with a neutralizing background. We study phase tran-
a
tion. In that case, the Poisson equation is replaced by a sitions between homogeneous and inhomogeneous states
screened Poisson equation [8]. In the gravitational anal- depending on the dimension of space. In d = 1, the
ogy,thisamountstoreplacingthegravitationalpotential systempresentscanonicalandmicrocanonicalsecondor-
by a screened gravitational potential, i.e. an attractive der phase transitions. In d = 2, the system presents
Yukawa potential. In that case, there exists interest- an isothermal collapse in CE (zeroth order phase transi-
ingphasetransitionsbetweenspatiallyhomogeneousand tion)andafirstorderphasetransitioninMCE.Ind=3,
spatially inhomogeneous equilibrium distributions. This the system presents an isothermal collapse in CE and a
is a physicalmotivationto considerthe thermodynamics gravothermal catastrophe in MCE (zeroth order phase
of N-body systems interacting via an attractive Yukawa transitions). In Sec. IV, we perform a similar study
potential [9]. This will be called the screened Newto- fortheattractiveYukawapotentialwithscreeninglength
nianmodel. Weshallalsoconsiderarelatedmodelwhere k−1. In d = 1, there exists a canonical tricritical point
0
the interaction is not screened but the Poisson equation (k ) R = √2π 4.44 and a microcanonical tricritical
0 c
≃
2
point (k ) R 11.8, where R is the system size. If with
0 m
≃
k < (k ) , the system presents canonical and micro-
0 0 c
f f
canonical second order phase transitions. In that case, S = k ln drdv, (4)
B
the ensembles are equivalent. If (k ) < k < (k ) , − m m
0 c 0 0 m Z
the system presents a canonical first order phase transi-
tionandamicrocanonicalsecondorderphasetransition.
In that case, there exists a region of negative specific M = ρdr, (5)
heats in MCE and the ensembles are inequivalent. If Z
k > (k ) , the system presents canonical and micro-
0 0 m
canonical first order phase transitions. In d = 2 and
v2 1
d=3,the phasetransitionsaresimilartothosereported E = f drdv+ ρ(r,t)u(r,r′)ρ(r′,t)drdr′
2 2
for the modified Newtonian model. In Sec. V, we study Z Z
the dynamical stability of the homogeneous phase and + ρV dr, (6)
analytically determine the critical point (E∗,T∗) that
c c Z
markstheonsetofinstabilityofthehomogeneousbranch
and the starting point of the bifurcated inhomogeneous where ρ(r,t) = f(r,v,t)dv is the spatial density. In-
branch. Direct numerical simulations associated with troducing the mean field potential
R
these phase transitions will be reported in a forthcom-
ing paper. Φ(r)= u(r,r′)ρ(r′)dr′+V(r), (7)
Finally, it may be noted that the phase transitions
Z
reported in this paper share analogies (but also differ-
the energy can also be written
ences) with phase transitions observed in the Hamil-
tonian mean field (HMF) model [32–36], the spherical
1 1
mass shell (SMS) model [37], the Blume-Emery-Griffiths E = fv2drdv+ ρ(Φ+V)dr. (8)
2 2
(BEG)model[38],theinfinite-rangeattactiveinteraction Z Z
(IRAI) model [39], the self-gravitating Fermi gas (SGF) Weshallbeinterestedinglobalandlocalentropymax-
model [28], the self-gravitating ring (SGR) model [40] ima. Let us first determine the critical points of entropy
and the one-dimensional static cosmology (OSC) model atfixedmassandenergywhichcancelthefirstordervari-
[41]. ations. Introducing Lagrange multipliers, they satisfy
1
δS δE αδM =0. (9)
II. KINETIC MODELS AND STATISTICAL
− T −
EQUILIBRIUM STATES
The variations are straightforward to evaluate and we
1. Isolated systems obtain the mean field Maxwell-Boltzmann distribution
We consideranisolatedsystemofN particlesininter- f =Ae−βm v22+Φ , (10)
action described by the Hamiltonian equations (cid:0) (cid:1)
whereβ =1/k T andΦ(r)isgivenbyEq. (7). Integrat-
dri ∂H dvi ∂H ingovertheveBlocity,thefindthatthedensityisgivenby
m = , m = , (1)
dt ∂vi dt −∂ri the mean field Boltzmann distribution
where ρ=A′e−kmBΦT. (11)
1
H = mv2+m2 u(r ,r )+m V(r ). (2)
2 i i j i This criticalpoint is a (local)entropymaximum atfixed
i i<j i
X X X mass and energy iff
We assume that the particles interact through a binary
potential u(r,r′) that is symmetric with respect to the (δf)2 1
interchange of r and r′, and that they also evolve in a δ2J =− 2mf drdv− 2β δρδΦdr≤0, (12)
fixed external potential V(r). Since the system is iso- Z Z
lated, with strict conservation of energy and mass, the for all perturbations δf that conserve mass and energy
proper statistical ensemble is the microcanonical ensem- at first order. In Appendix A, we provide an equivalent
ble [9]. Inthis paper, we shalluse a mean fieldapproach but simpler condition of stability in the microcanonical
[63]. Inthemicrocanonicalensemble,thestatisticalequi- ensemble [see inequality (A12)].
librium state is obtained by maximizing the Boltzmann The time evolution of the distribution function
entropy at fixed mass and energy. We thus have to solve f(r,v,t) is governedby a kinetic equation of the form
the maximization problem
∂f ∂f ∂f ∂f
max S[f] E[f]=E, M[f]=M , (3) +v Φ = , (13)
f { | } ∂t · ∂v −∇ · ∂r (cid:18)∂t(cid:19)coll
3
where with
v2 1
Φ(r,t)= u(r,r′)ρ(r′,t)dr′+V(r), (14) F = f drdv+ ρ(r,t)u(r,r′)ρ(r′,t)drdr′
2 2
Z Z Z
f f
+ ρV dr+k T ln drdv. (18)
is the time-dependent mean field potential. The l.h.s. is B
m m
anadvectiveoperator(Vlasov)inphasespace. Ther.h.s. Z Z
is a “collision” operator like the Boltzmann operator in We shallbe interestedbyglobalandlocalminimaoffree
thekinetictheoryofgasesorliketheLandau(orLenard- energy. Let us first determine the critical points of free
Balescu) operator in plasma physics or stellar dynamics. energy at fixed mass which cancel the first order varia-
The “collision” operator in Eq. (13) takes into account tions. Introducing a Lagrange multiplier, they satisfy
the development of correlations between particles. It
δF +αTδM =0. (19)
can have a more or less complicated form but it satis-
fies general properties associated with the first and sec-
The variations are straightforward to evaluate and we
ond principles of thermodynamics: (i) it conserves mass
obtain the mean field Maxwell-Boltzmann distribution
and energy; (ii) it satisfies an H-theorem for the Boltz-
(10) and the mean field Boltzmann distribution (11) as
mann entropy (4), i.e. S˙ 0 with an equality iff f is in the microcanonical ensemble. This critical point is a
≥
the Maxwell-Boltzmann distribution (10). Furthermore, (local) minimum of free energy iff
the Maxwell-Boltzmann distribution is dynamically sta-
ble iff it is a (local) entropy maximum at fixed mass and 1 k T (δf)2
δ2F = δρδΦdr+ B drdv 0, (20)
energy. These generalproperties canbe checkeddirectly 2 m 2f ≥
fortheBoltzmannequation,fortheLandauequation,for Z Z
the Lenard-Balescuequation and for the BGK operator. for all perturbations δf that conserve mass.
Therefore, the kinetic equation (13) is consistent with In the mean field approximation, the evolution of the
the maximization problem (3) describing the statistical distribution function f(r,v,t) is governed by a kinetic
equilibriumstateofthesysteminMCE.Ifweneglectthe equation of the form
collisionsfor sufficiently shorttimes, Eq. (13) reduces to
∂f ∂f ∂f ∂ ∂f
the Vlasov equation which can experience a complicated +v Φ = D +ξfv , (21)
∂t · ∂v −∇ · ∂r ∂v ∂v
processof collisionlessviolentrelaxationtowardsa quasi (cid:18) (cid:19)
stationary state (QSS) [43].
coupled to the mean field potential (14). This is called
the mean field Kramers equation. The mean field
Kramers equation conserves mass and satisfies an H-
2. Dissipative systems in phase space theorem for the Boltzmann free energy (18), i.e. F˙ 0
≤
withanequalityifff istheMaxwell-Boltzmanndistribu-
We consider a dissipative system of N Brownian par- tion (10). Furthermore, the Maxwell-Boltzmann distri-
ticles in interactiondescribedby the Langevinequations butionisdynamicallystableiffitisa(local)minimumof
freeenergyatfixedmass. Therefore,thekineticequation
dr ∂H (21)isconsistentwiththeminimizationproblem(17)de-
i
m = , (15)
dt ∂v scribing the statisticalequilibriumstate of the systemin
i
CE.
Remark: the critical points in MCE and CE are the
dv 1 ∂H same because the variational problems (3) and (17) are
i = ξv +√2DR (t), (16)
dt −m∂r − i i equivalentatthelevelofthefirstordervariations(9)and
i
(19). However,theyarenotequivalentatthe levelofthe
where H is the Hamiltonian defined by Eq. (2), ξv second order variations (12) and (20) because of the dif-
i
is a friction force and R (t) is a white noise satis−fying ferent class of perturbations to consider. Therefore, we
i
R (t) =0and Rµ(t)Rν(t) =δ δ δ(t t′). Thediffu- canhaveensembles inequivalence[22,31,44,45]. Infact,
h i i h i j i ij µν −
sioncoefficientDandthefrictioncoefficientξ arerelated the condition of canonical stability (17) provides a suf-
toeachotheraccordingtotheEinsteinrelationξ =Dβm ficient condition of microcanonical stability (3). Indeed,
where β = 1/(k T) is the inverse temperature. Since if inequality (20) is satisfied for all perturbations that
B
this systemisdissipative,the properstatisticalensemble conservemass, then it is a fortiori satisfiedfor perturba-
is the canonicalensemble [9]. In the canonicalensemble, tions that conserve mass and energy, so that inequality
the statistical equilibrium state is obtained by minimiz- (12) is satisfied. Therefore, canonical stability implies
ing the Boltzmann free energy F[f] = E[f] TS[f] at microcanonical stability:
−
fixedmass. Wethushavetosolvetheminimizationprob-
(17) (3). (22)
lem
⇒
However, the converse is wrong in case of ensembles in-
min F[f] M[f]=M (17)
f { | } equivalence.
4
3. Dissipative systems in physical space and the time-dependent density ρ(r,t) is solution of the
Smoluchowskiequation(28). Using Eq. (29), we can ex-
press the free energy (18) as a functional of the density
In the strong friction limit ξ + , we can formally
neglecttheinertialtermdv /dti→nEq.∞(16)andweobtain and we obtain the free energy (25) up to some unimpor-
i
tant constants.
the overdamped Langevin equations
Remark 2: it is shown in Appendix A that the max-
dr 1 ∂H imization problems (17) and (24) are equivalent in the
ξ dti =−m∂r +√2DRi(t). (23) sense thatf(r,v) is solutionof (17)iff ρ(r) is solutionof
i
(24). Thus, we have
Thestatisticalequilibriumstateofthissystem(described
by the canonicalensemble [9]) is obtained by solving the (17) (24). (30)
⇔
minimization problem
As a consequence, the Maxwell-Boltzmann distribution
min F[ρ] M[ρ]=M , (24) f(r,v) is dynamically stable with respect to the mean
ρ { | } field Kramers equation (21) iff the corresponding Boltz-
mann distribution ρ(r) is dynamically stable with re-
with
spect to the mean field Smoluchowskiequation (28). On
1 the other hand, according to the implication (22), the
F = ρ(r,t)u(r,r′)ρ(r′,t)drdr′ Maxwell-Boltzmann distribution f(r,v) is dynamically
2
Z stable with respect to the kinetic equation (13) if it is
ρ ρ
+ ρV dr+k T ln dr. (25) stable with respect to the mean field Kramers equation
B
m m
Z Z (21), but the reciprocal is wrong in case of ensembles
inequivalence.
Writing the variational principle as
δF +αTδM =0, (26)
4. The Keller-Segel model of chemotaxis
we obtain the mean field Boltzmann distribution (11).
This criticalpoint is a (local)minimum offree energy at The Keller-Segel model [3] describing the chemotaxis
fixed mass iff of biological populations can be written as
1 k T (δρ)2 ∂ρ
δ2F = δρδΦdr+ B dr 0, (27) = (D ρ χρ c), (31)
2 m 2ρ ≥ ∂t ∇· ∇ − ∇
Z Z
for all perturbations δρ that conserves mass.
1 ∂c
In the mean field approximation, the evolution of the =∆c k2c+λρ, (32)
densityprofileρ(r,t)isgovernedbyakineticequationof D′ ∂t −
the form
whereρistheconcentrationofthebiologicalspecies(e.g.
bacteria)andcistheconcentrationofthesecretedchem-
∂ρ 1 k T
B
= ρ+ρ Φ , (28) ical. The bacteria diffuse with a diffusion coefficient D
∂t ∇· ξ m ∇ ∇
(cid:20) (cid:18) (cid:19)(cid:21) and undergo a chemotactic drift with strength χ along
the gradient of chemical. The chemical is produced by
coupled to the mean field equation (14). This is called
the bacteria at a rate D′λ, is degraded at a rate D′k2
the mean field Smoluchowski equation. The mean field
and diffuses with a diffusion coefficient D′. We adopt
Smoluchowskiequation(28)conservesmassandsatisfies
Neumann boundary conditions [3]:
an H-theorem for the Boltzmann free energy (25), i.e.
F˙ 0 with an equality iff ρ is the Boltzmann distribu- c n=0, ρ n=0, (33)
≤
tion (11). Furthermore, the Boltzmann distribution is ∇ · ∇ ·
dynamically stable iff it is a (local) minimum of free en- where n is a unit vector normal to the boundary of the
ergyatfixedmass. Therefore,thekineticequation(28)is domain. The drift-diffusion equation (31) is similar to
consistentwiththeminimizationproblem(24)describing the mean field Smoluchowski equation (28) where the
the statistical equilibrium state of the system in CE. concentration of chemical c(r,t) plays the role of the
Remark1: theSmoluchowskiequation(28)canalsobe potential Φ(r,t). Therefore−, there exists many analogies
deduced from the Kramers equation (21) in the strong between chemotaxis and Brownian particles in interac-
friction limit [46]. For ξ,D + and β = ξ/Dm fi- tion [6]. In particular, the effective statistical ensemble
nite, the time-dependent dist→ribut∞ion function f(r,v,t) associated with the Keller-Segel model is the canonical
is Maxwellian ensemble. The steady states of the Keller-Segel model
are of the form
f(r,v,t)= β2mπ d/2ρ(r,t)e−βmv22 +O(1/ξ), (29) ρ=AeDχc, (34)
(cid:18) (cid:19)
5
which is similar to the Boltzmann distribution (11) with Thisisvalidintheabsenceofdegradationofthechemical
an effective temperature T = D/χ. The Lyapunov andforsufficientlylargedensitiesρ ρ. Thismodelcan
eff
≫
functional associated with the KS model is [4]: be usedin particular to study chemotactic collapse. The
corresponding free energy is
1
F = ( c)2+k2c2 dr ρcdr
2λ ∇ −
Z (cid:2) +T(cid:3) ρZlnρdr. (35) F =−21 ρcdr+Teff ρlnρdr. (42)
eff Z Z
Z
It is similar to a free energy F = E T S in thermo- In that model, the boundary conditions (33) must be
eff
−
dynamics,whereE istheenergyandS istheBoltzmann modified [64] and we must impose that c 0 at infinity
→
entropy. The KS model conserves mass and satisfies an like for the gravitational potential in astrophysics. Fur-
H-theorem for the free energy (35), i.e. F˙ 0 with an thermore,wemustimposethatthenormalcomponentof
≤
equalityiffρistheBoltzmanndistribution(34). Further- the current vanishes on the boundary: (D ρ χρ c)
more, the Boltzmann distribution is dynamically stable n = 0 so as to conserve mass. In that ∇case−, the∇KS·
iff it is a (local) minimum of free energy at fixed mass. model is isomorphic to the Smoluchowski-Poisson (SP)
In that context, the minimization problem system describing self-gravitating Brownian particles in
the overdamped limit [29].
min F[ρ,c] M[ρ]=M , (36)
ρ,c { | }
determines a steady state of the KS model that is dy-
namically stable. This is similar to a condition of ther-
5. Physical justification of the canonical ensemble for
modynamical stability in the canonical ensemble.
systems with long-range interactions
Let us consider some simplified forms of the Keller-
Segel model that have been introduced in the literature:
(i)InthelimitoflargediffusivityofthechemicalD′ In statistical mechanics, the canonical distribution is
+ at fixed k2 and λ, the reaction-diffusion equatio→n usually derived by considering a subpart of a large sys-
(3∞2) takes the form of a screened Poissonequation [8]: tem and assuming that the rest of the system plays the
role of a thermostat [48]. However, this justification im-
∆c k2c= λρ, (37) plicitly assumes that energy is additive. Since energy is
− −
non-additive for systems with long-range interactions, it
and the free energy becomes is sometimes concluded that the canonical ensemble has
no foundation to describe systems with long-rangeinter-
1
F = ρcdr+T ρlnρdr. (38) actions[49]. Infact,thisisnotquitetrue[9]. Wecangive
eff
−2 two justifications of the canonical ensemble for systems
Z Z
with long-range interactions:
In that case, the KS model is isomorphic to the Smolu-
chowski equation (28) with an attractive Yukawa poten- (i) The canonical ensemble is relevant to describe a
tial (64). system of particles in contact with a thermal bath of a
(ii) In the limit of large diffusivity of the chemical different nature [9]. This is the case if we consider a
D′ + and a vanishing degradation rate k2 = 0, system of Brownian particles in interaction described by
the→reacti∞on-diffusion equation (32) takes the form of a thestochasticequations(15)-(16). Theparticlesinteract
modified Poisson equation [11]: through a potential u(r,r′) that can be long-range, but
they also undergo a friction force and a stochastic force
∆c= λ(ρ ρ), (39) that are due to other types of interaction (they model
− −
in general short-range interactions). As we have seen,
where ρ = M/V is the average density, and the free en- this system is described by the canonical ensemble. It
ergy becomes does not correspond to a subsystem of a larger system,
but simply to a system as a whole with long-range and
1
F = (ρ ρ)cdr+T ρlnρdr. (40) short-range interactions [65].
eff
−2 −
Z Z (ii) Since canonical stability implies microcanonical
In that case, the KS model is isomorphic to the Smolu- stability[44],theconditionofcanonicalstabilityprovides
chowski equation (28) with a modified Poisson equation a sufficient condition of microcanonical stability. In this
(43). sense,thecanonicalstabilitycriterion(seeSecs. II 2and
(iii)Someauthorshavealsoconsideredasimplemodel II 3) can be useful even for an isolated Hamiltonian sys-
of chemotaxis where the reaction-diffusion equation (32) tem (see Sec. II 1) because if we can prove that this
is replaced by the Poisson equation [47]: system is canonically stable, then it is granted to be mi-
crocanonically stable. This remark also applies to other
∆c= λρ. (41) ensembles (grand canonical, grand microcanonical,...).
−
6
III. THE MODIFIED NEWTONIAN MODEL is the CE since the KS model has a canonical structure.
Thismodelhasbeenstudiedbyappliedmathematicians,
In this section, we discuss phase transitions that ap- starting with J¨ager & Luckhaus [11], but they have not
pear in the modified Newtonian model. performed the type of study that we are developing in
this paper.
In view of these different applications, we shall study
A. Physical motivation of the model thismodelinthemicrocanonicalandcanonicalensembles
in any dimension of space.
We consider a system of particles interacting via a
mean field potential Φ(r,t) that is solution of the modi-
B. The modified Emden equation
fied Poissonequation
∆Φ=S G(ρ ρ), (43) In the modified Newtonianmodel, the statistical equi-
d
−
libriumstateisgivenbytheBoltzmanndistribution(44)
where ρ=M/V is the averagedensity (conserved quan- coupled to the modified Poisson equation (43). We look
tity). At statistical equilibrium, the density is given by for spherically symmetric solutions because, for non ro-
the Boltzmann distribution tating systems, entropy maxima (or minima of free en-
ergy) are spherically symmetric. Introducing the cen-
ρ=Ae−βmΦ. (44)
tral density ρ = ρ(0), the central potential Φ = Φ(0),
0 0
the new field ψ = βm(Φ Φ ) and the scaled distance
0
We have used the notations of astrophysics (where G is −
ξ = (S Gβmρ )1/2r, the Boltzmann distribution (44)
d 0
theconstantofgravityandS thesurfaceofaunitsphere
d can be rewritten
in d dimensions) in order to make the connection with
ordinary self-gravitating systems where Eq. (43) is re- ρ=ρ e−ψ(ξ). (46)
0
placed by the Poisson equation ∆Φ = S Gρ. However,
d
this model can have application in other contexts as ex- Substituting this relation in the modified Poisson equa-
plained below. We assume that the system is confined tion (43), we obtain the modified Emden equation
in a finite domain (box) and we impose the Neumann
boundary conditions 1 d dψ
ξd−1 =e−ψ λ, (47)
ξd−1dξ dξ −
Φ n=0, ρ n=0, (45) (cid:18) (cid:19)
∇ · ∇ ·
where λ=ρ/ρ plays the role of the inverse centralden-
where n is a unit vector normal to the boundary of the 0
sity. Since Φ′(0)=0 for a spherically symmetric system,
box (the explicit expression of the potential in d = 1 is
the boundary conditions at the origin are
given in Appendix B). This model admits spatially ho-
mogeneous solutions (ρ=ρ and Φ=0) at any tempera- ψ(0)=ψ′(0)=0. (48)
ture. Italsoadmitsspatiallyinhomogeneoussolutionsat
sufficiently low temperatures. We shall study this model TheordinaryEmdenequation[53]isrecoveredforλ=0,
in arbitrary dimensions of space d with explicit compu- i.e. for very large central densities with respect to the
tations for d = 1,2,3. This model has different physical average density. The function e−ψ(ξ) is plotted in Figs.
applications: 1 and 2 for different values of λ and different dimensions
(i) It describes self-gravitating systems in a cosmolog- of space d. It presents an infinity of oscillations. For
ical setting [12]. Due to the expansion of the universe, d = 1, the oscillations are undamped and their period
when we work in the comoving frame, the Poissonequa- is given by Eq. (E11). For d 2, the oscillations are
tion takes the form of Eq. (43) where the potential is damped and the function ψ(ξ)≥tends to the asymptotic
produced by the deviation between the actual density value lnλ for ξ + .
ρ(r,t) and the mean density ρ. In cosmology, we must We−assume tha→t the∞system is enclosed in a spheri-
also account for the scale factor a(t) but if we consider cal box of radius R. The normalized box radius α =
timescalesthatareshortwithrespecttotheHubbletime (S Gβmρ )1/2R is determined by the boundary condi-
d 0
H−1 = a/a˙, we can ignore this time dependence. This tion Φ′(R)=0 that becomes
modelhasbeenstudiedbyValageas[41]ind=1withpe-
riodicboundaryconditions. Inthatcontext,therelevant ψ′(α)=0. (49)
ensemble is the MCE since the system is isolated.
(ii) By a proper reinterpretation of the parameters, For a given value of λ, we need to integrate the modi-
thefieldequation(43)describestherelationbetweenthe fiedEmdenequation(47)-(48)untilthepointξ =αsuch
concentrationofthe chemicalandthedensityofbacteria that ψ′(α) = 0. Since the function ψ(ξ) presents an in-
in the Keller-Segel model (39). In that case, the most finite number of oscillations, this determines an infinity
physicaldimensionisd=2andtheboundaryconditions of solutions α (λ), α (λ),... that will correspond to dif-
1 2
areofthe form(45). Furthermore,therelevantensemble ferent branches in the following diagrams. Once α (λ)
n
7
4 400
λ > 1 d=1 d=1
3 300
n=3
−ψ 2 η 200
e
n=2
1 100
λ < 1
n=1
0 0
0 10 20 30 40 50 -15 -10 -5 0 5
ξ ln(1/λ)
FIG.1: Thefunctione−ψ ford=1andλ=0.5<1(bottom) FIG. 3: Inverse temperature η as a function of the central
or λ=2>1 (top). In d=1, theoscillations are undamped. density 1/λ for the first three branchesin d=1.
3
concentrations, the Newtonian solution provides a good
d=2
2.5 λ > 1 approximationofthemodifiedNewtoniansolutioninthe
core (see Appendix E).
2
ψ-e1.5 C. The temperature
1
λ < 1 Wemustnowrelatethenormalizedcentraldensity1/λ
to the temperature T. Recalling that ρ = M/V with
0.5
V = 1S Rd, we obtain
d d
0
0 25 50 75 100 ρ dM 1 GMmβ 1
ξ λ= = =d . (50)
ρ S Rdρ Rd−2 α2
0 d 0
FIG.2: Thefunctione−ψ ford=2andλ=0.5<1(bottom) Introducing the normalized temperature
or λ = 2 > 1 (top). In d 2, the oscillations are damped.
≥
The case d=3 (not represented) is similar. βGMm
η , (51)
≡ Rd−2
we find the relation
is determined, the density profile is given by Eq. (46).
1
The density profile is extremum at the center and at the η = λα2. (52)
d
boundary. Onthe n-th branch,the density profile shows
n “clusters” corresponding to the oscillations of e−ψ(ξ). Recalling that α=α (λ) for the n-th branch, this equa-
n
Close to the origin, the density increases for λ>1 while tiongivestherelationbetweentheinversetemperatureη
itdecreasesforλ<1. Thehomogeneousstateψ =0cor- andthe centraldensity 1/λfor the n-thbranch. In Figs.
responds to λ = 1. This solution is degenerate because 3, 4 and 5, we plot the inverse temperature η as a func-
the boundary condition (49) is satisfied for any α. tionofthecentraldensity1/λforthefirstthreebranches
Remark: Whenλ 0,correspondingtolargevaluesof n=1,2,3 in different dimensions of space d=1,2,3.
→
thecentraldensity,weexpecttoobtainresultssimilarto Let us discuss the asymptotic behaviors of the tem-
those obtained for the usual Newtonian model since the perature (we only describe the first branch n = 1) and
differentialequation(47)reduces to the ordinaryEmden compare with the Newtonian model (see, e.g., [29]):
equation. However, the results are different because the In d = 1: for the ordinary Newtonian model, the
•
boundary conditions are not the same. In the Newto- series of equilibria is parameterized by α, which is a
nian model, the force at the boundary is non zero (for measure of the central density. When α + , the
→ ∞
a spherically symmetric system, according to the Gauss distribution tends to a Dirac peak ρ = Mδ(x) and the
theorem,wehaveΦ′(R)=GM/Rd−1)whileinthe mod- inverse temperature η + . When α 0, the dis-
→ ∞ →
ified Newtonian model the force at the boundary is zero tribution is homogeneous and the inverse temperature
(Φ′(R)=0). Therefore,strictlyspeaking,theNewtonian η 0. For the modified Newtonian model, the se-
→
and the modified Newtonian models behave differently ries of equilibria is parameterized by the central den-
even when ρ + . Nevertheless, for large central sity 1/λ. When 1/λ + , the distribution tends to
0
→ ∞ → ∞
8
100 In d=3: for the ordinary Newtonian model, the se-
d=2 rie•sofequilibriais parameterizedby α. When α + ,
→ ∞
the distribution tends to the singular isothermal sphere
75 n=3 ρ (r) =1/(2πGβmr2) and the inverse temperature η
s
→
η = 2. The curve η(α) displays damped oscillations
s
around this value. When α 0, the distribution is ho-
η 50 →
mogeneous and the inverse temperature η 0. For the
→
modified Newtonian model, the series of equilibria is pa-
n=2
25 rameterized by 1/λ. When 1/λ + , the distribution
→ ∞
is concentrated at the center and we numerically find
ηc=4 n=1 that η 3.05... (the value is different from the Newto-
→
0 nianresultη =2 due to different boundaryconditions).
-15 -10 -5 0 5 10 s
ln(1/λ) The curve η(λ) displays damped oscillations around this
value. Whenλ=1,thedistributionis homogeneousand
FIG. 4: Inverse temperature η as a function of the central η = η∗ = 1x2 6.7302445 (see Appendix F). When
density 1/λ for thefirst three branchesin d=2. 1/λ c0, th3e d1is≃tribution is concentrated at the bound-
→
ary and we numerically find that η + .
→ ∞
50
d=3
D. The energy
40
n=3 Wemustalsorelatethenormalizedcentraldensity1/λ
30 to the energy E. The total energy is given by (see Ap-
η pendix C):
20
n=2 v2 1
E = f drdv+ (ρ ρ)Φdr. (53)
2 2 −
10 Z Z
n=1
Using the Maxwell-Boltzmann distribution (10), the ki-
0-10 -5 0 5 10 15 20 netic energy is simply
ln(1/λ)
d
K = Nk T. (54)
FIG. 5: Inverse temperature η as a function of the central 2 B
density 1/λ for thefirst three branchesin d=3.
UsingthemodifiedPoissonequation(43)andanintegra-
tion by parts, the potential energy can be written
a Dirac peak ρ = Mδ(x) and η → +∞ with the same W = 1 ( Φ)2dr. (55)
asymptoticbehaviorasintheNewtonianmodel(seeAp- −2S G ∇
d Z
pendix E). When λ = 1, the distribution is homoge-
neous and η = η∗ = π2 9.8696044 (see Appendix F). The total energy E =K+W is therefore given by
c ≃
When 1/λ 0, the distribution tends to a Dirac peak
ρ= M(δ(x→R)+δ(x+R))concentratedatthe boxand E = dNk T 1 ( Φ)2dr. (56)
η +2 . − 2 B − 2SdGZ ∇
→ ∞
In d = 2: for the ordinary Newtonian model, the
• Introducing the dimensionless variables defined previ-
series of equilibria is parameterized by α. When α
+ , the distribution tends to a Dirac peak ρ = Mδ(→r) ously, recalling that r =ξR/α, and introducing the nor-
∞ malized energy
and the inverse temperature tends to η = 4. When
c
α 0, the distribution is homogeneous and the inverse ERd−2
→
temperature η 0. For the modified Newtonian model, Λ , (57)
→ ≡− GM2
the series of equilibria is parameterized by 1/λ. When
1/λ + , the distribution tends to a Dirac peak ρ = we obtain
Mδ(→r) an∞d η η = 4 (since the density is very much
c
concentrated,→theboundaryconditionsdonotmatterand d 1 1 α dψ 2
Λ= + ξd−1dξ. (58)
we recover the same results as in the Newtonian case). −2η 2η2αd−2 dξ
Z0 (cid:18) (cid:19)
When λ = 1, the distribution is homogeneous and η =
η∗ = 1j2 7.3410008 (see Appendix F). When 1/λ Recalling that α = α (λ) and η = η (λ) for the n-th
c 2 11 ≃ → n n
0, the distribution is concentrated at the boundary and branch, this equation gives the relation between the en-
η + . ergy Λ and the central density 1/λ for the n-th branch.
→ ∞
9
0.2
d=1 n=1 0 n=3
0.15
-0.1
n=2
0.1
-0.2
Λ Λ
0.05 n=2
n=3 -0.3
d=3 n=1
0
-0.4
-0.05
-15 -10 -5 ln(1/λ) 0 5 10 -15 -10 -5 0ln(1/λ)5 10 15 20
FIG.6: EnergyΛasafunctionof thecentraldensity1/λfor FIG.8: EnergyΛasafunctionofthecentraldensity1/λ for
thefirst threebranches in d=1. thefirst threebranchesin d=3.
0.4 d=2 n=1 riesofequilibriais parameterizedby α. When α + ,
the distribution tends to a Dirac peak ρ = Mδ→(r) a∞nd
the energy Λ + . When α 0, the distribution is
→ ∞ →
homogeneousandtheenergyΛ . Forthemodified
0.2 →−∞
Newtonian model, the series of equilibria is parameter-
Λ
n=2 ized by 1/λ. When 1/λ + , the distribution tends
to a Dirac peak ρ = Mδ(→r) an∞d Λ + . When λ = 1
0 n=3 the distribution is homogeneous an→d Λ=∞Λ∗ = 1/η∗
c − c ≃
0.13622121. When 1/λ 0, the distribution is con-
− →
centrated at the boundary and we numerically find that
-0.2 Λ 0.1.
-15 -10 -5 0 5 10
ln(1/λ) →
In d = 3: for the ordinary Newtonian model, the
•
series of equilibria is parameterized by α. When α
FIG.7: EnergyΛasafunctionof thecentraldensity1/λfor →
+ , the distribution tends to the singular isothermal
thefirst threebranches in d=2. sp∞hereρ (r)=1/(2πGβmr2)withenergyΛ =1/4. The
s s
curve Λ(α) undergoes damped oscillations around this
value. When α 0, the distribution is homogeneous
→
In Figs. 6, 7 and 8, we plot the normalized energy and the energy Λ . For the modified Newtonian
→ −∞
Λ as a function of the central density 1/λ for the first model, the series of equilibria is parameterized by 1/λ.
threebranchesn=1,2,3indifferentdimensionsofspace When1/λ + ,thedistributionisconcentratedatthe
→ ∞
d=1,2,3. center and we numerically find that Λ 0.38... (the
→ −
Letusconsidertheasymptoticbehaviorsoftheenergy value is different from the Newtonian result Λs = 1/4
(we only describe the first branch n = 1) and compare due to different boundary conditions). The curve Λ(λ)
with the Newtonian model (see, e.g., [29]): undergoes damped oscillations around this value. When
In d=1: for the ordinary Newtonian model, the se- λ = 1 the distribution is homogeneous and Λ = Λ∗c =
rie•s of equilibria is parameterized by α, which is a mea- −3/(2ηc∗)≃−0.22287452. When λ→0, the distribution
sure of the central density. When α + , the distri- is concentratedatthe boundaryandwe numericallyfind
bution tends to a Dirac peak ρ=Mδ→(x) a∞nd the energy that Λ 0.05.
→
Λ 0. When α 0, the distribution is homogeneous
→ →
and the energy Λ . For the modified Newtonian
→ −∞
model, the series of equilibria is parameterized by the
central density 1/λ. When 1/λ + , the distribution E. The entropy and the free energy
→ ∞
(1)
tends to a Dirac peak ρ = Mδ(x) and Λ Λ = 1/6
max
→
(see Appendix D). When λ = 1 the distribution is ho- Finally,werelatethecentraldensity1/λtotheentropy
mogeneous and Λ = Λ∗c = −1/(2ηc∗) ≃ −0.0506606. SandtothefreeenergyF. UsingEqs. (4),(10)and(11),
When 1/λ 0, the distribution tends to a Dirac peak the entropy is given by
ρ= M(δ(x→R)+δ(x+R))concentratedatthe boxand
2 −
Λ Λ(1) . d ρ ρ
→ max S = Nk lnT k ln dr. (59)
In d=2: for the ordinary Newtonian model, the se- 2 B − B m m
• Z
10
0 -1
-0.5 λ→ ←0→λ → ∞ d=1 ηc=4 d=2 → λ→ ∞
-2 n=1
-1 n=1
n=2 n=2
λ=1
j-1.5 j-3
λ=1
-2
-4
-2.5 →
λ→ 0
-3 -5
0 20 40 60 80 100 120 0 10 20 30 40 50 60
η η
FIG.9: Freeenergy J asafunctionoftheinversetemper- FIG. 10: Free energy J as a function of the inverse tem-
NkB NkB
ature η in d = 1. Note that the branches λ < 1 and λ > 1 peratureη in d=2.
coincide.
Substituting Eq. (46) in Eq. (59), and introducing the -2 -3.15
dimensionless variables defined previously, we get -3.2
S d -3 j-3.25
= lnβ lnρ0 λ→ 0 -3.3
NkB −2 − j-4 → λ=1
ρ R d α -3.35
+ 0 S ψe−ψξd−1dξ, (60) 2.7 2.8 2.9η3 3.1 3.2
d -5
Nm (cid:18)α(cid:19) Z0 n=1
→ λ→ ∞
up to some unimportant constants. Using α = -6 d=3
n=2
(S Gβmρ )1/2Rtoexpressρ intermsofαandintroduc-
d 0 0
ing the normalized temperature (51), we finally obtain -70 5 10 15 20 25 30 35
η
S d 2
= − lnη 2lnα FIG. 11: Free energy J as a function of the inverse tem-
NkB − 2 − NkB
1 1 α peratureη in d=3.
+ ψe−ψξd−1dξ, (61)
ηαd−2
Z0
up to some unimportant constants. Using the previous
results,thisexpressionrelatesthe entropyS/Nk tothe curves will be helpful in the next section to analyze the
B
central density 1/λ. The free energy is F =E TS. In phasetransitionsinthecanonicalandmicrocanonicalen-
−
thefollowing,itwillbemoreconvenienttoworkinterms sembles respectively.
of the Massieu function J = S k βE (by an abuse of
B Remark: Since δS = k βδE, the extrema of entropy
− B
language, we shall often refer to J as the free energy).
S(λ) and energy E(λ) coincide. Since the series of equi-
We have
libria E(λ) exhibits damped oscillations for 1/λ +
→ ∞
J S in d = 3 (see Fig. 8), this implies that the curve S(λ)
= +ηΛ. (62) will also exhibit damped oscillations at the same loca-
Nk Nk
B B
tions. Correspondingly, S(E) will present some “spikes”
Usingthepreviousresults,thisexpressionrelatesthefree for 1/λ + in d=3 (see inset of Fig. 15). Similarly,
→ ∞
energy J/Nk to the central density 1/λ. since δJ = Ek δβ, the extrema of free energy J(λ)
B B
−
In Figs. 9, 10 and 11, we have plotted the free en- and temperature β(λ) coincide. Since the series of equi-
ergy J/Nk as a function of the inverse temperature η libriaβ(λ)undergoesdampedoscillationsfor1/λ +
B
→ ∞
(parameterizedby the centraldensity 1/λ) in d=1,2,3. in d = 3 (see Fig. 5), this implies that the curve J(λ)
In Figs. 12, 13, 14 and 15, we have plotted the entropy will also exhibit damped oscillations at the same loca-
S/Nk as a function of the energy Λ (parameterized by tion, and that the curve J(β) will present some “spikes”
B
the central density 1/λ) in d = 1,2,3. In these figures, for 1/λ + in d = 3 (see inset of Fig. 11). In ad-
→ ∞
the solid lines without label refer to the homogeneous dition, the curve J(β) presents a minimum for η 24.7
≃
phase. The solid lines with label n = 1 refer to the corresponding to E = 0. Similar behaviors were previ-
firstinhomogeneousbranch. The dashed lines with label ously observed in the model of self-gravitating fermions
n= 2 refer to the second inhomogeneous branch. These [28, 31].