Symmetry Breaking in
d
Dimensional SelfGravitating Systems
Renato Pakter,
1
Bruno Marcos,
1,2
and Yan Levin
1
1
Instituto de Fı´ sica, Universidade Federal do Rio Grande do Sul, Caixa Postal 15051, CEP 91501970,Porto Alegre, Rio Grande do Sul, Brazil
2
Universite´ de Nice SophiaAntipolis, CNRS, Laboratoire J.A. Dieudonne´ , UMR 7351, Parc Valrose,06108 Nice Cedex 02, France
(Received 14 June 2013; revised manuscript received 18 November 2013; published 6 December 2013)Systems with longrange interactions, such as selfgravitating clusters and magnetically conﬁnedplasmas, do not relax to the usual BoltzmannGibbs thermodynamic equilibrium, but become trappedin quasistationary states (QSS) the lifetime of which diverges with the number of particles. The QSS arecharacterized by the lack of ergodicity which can result in a symmetry broken QSS starting from aspherically symmetric particle distribution. We will present a theory which allows us to quantitativelypredict the instability threshold for spontaneous symmetry breaking for a class of
d
dimensionalselfgravitating systems.
DOI: 10.1103/PhysRevLett.111.230603 PACS numbers: 05.20.
y, 05.45.
a, 05.70.Ln
Lord Rayleigh was probably the ﬁrst to make an observation that longrange forces can lead to symmetry breaking [1]. Rayleigh was studying the stability of conductingspherical ﬂuid droplets carrying charge
Q
. He discoveredthat when
Q
exceeds a certain critical threshold
Q
c
, droplets become unstable to symmetry breaking perturbations,elongating and eventually breaking up, emitting jets of ﬂuid that carry away a signiﬁcant fraction of the charge[2]. Rayleigh instability is now the basis for various technological applications, such as electrospraying and electrospinning. It also helps to understand the conformationalstructure of charged polymers, such as polyampholytes [3].For selfgravitating systems a similar instability has beenobserved in gravitational simulations [4]. It has been foundthat an initially spherically symmetric selfgravitating system can become unstable, leading to formation of structures of reduced symmetry [4]. This radial orbit instabilityis believed to be important for the formation of ellipticalgalaxies [5].There is, however, a fundamental difference between theRayleigh instability of charged spherical droplets and theinstability of spherically symmetric selfgravitating systems. Since the droplets are in (canonical) thermodynamicequilibrium, their shape must correspond to the minimumof the Helmholtz free energy—in fact, even for
Q
somewhat below
Q
c
, a spherical shape is already metastable,with the global minimum corresponding to a stronglyprolate ellipsoid [6]. The thermal ﬂuctuations, however,are too small to overcome the barrier that separates themetastable minimum from the global one, so that thespherical shape persists up to the Rayleigh threshold. Onthe other hand, gravitational systems are intrinsicallymicrocanonical—isolated from environment [7–9]. In the
thermodynamic limit, such longrange systems do notevolve to thermodynamic equilibrium but become trappedin quasistationary states (QSS), the lifetime of whichdiverges with the number of particles [10]. The QSS arecharacterized by the broken ergodicity, making equilibrium statistical mechanics inapplicable [11]. To explorespontaneous symmetry breaking of systems with longrange forces, therefore, requires a completely differentapproach [12]. In this Letter wewill present a theory whichallows us to quantitatively predict the thresholds of symmetry breaking instabilities for systems with longrangeinteractions. The results of the theory will be comparedwith extensive molecular dynamics simulations.To present the theory, we will study a class of selfgravitating systems of
N
particles of mass
m
in an inﬁnite
d
dimensional space. The interaction potentialbetween theparticles is
V
ð
r
Þ¼
Gm
2
=
½ð
2
d
Þ
r
d
2
, where
G
is thegravitational constant. We will work in thermodynamiclimit,
N
!1
and
m
!
0
, while the total mass
M
Nm
remainsﬁxed.The initialparticle distributionisassumed tobe a uniform spherically symmetric waterbag in both conﬁguration and velocity space,
f
0
ð
r
;
v
Þ¼
d
2
C
2
d
r
dm
v
dm
ð
r
m
r
Þ
ð
v
m
v
Þ
;
(1)where
is the Heaviside step function,
C
d
¼
2
d=
2
=
ð
d=
2
Þ
is the surface area of a
d
dimensional unitsphere, and
ð
x
Þ
is the gamma function. Since the initialwaterbag distribution is not a stationary solution of thecollisionless Boltzmann (Vlasov) equation, the systemswill evolve with time. We are interested in discoveringunder what conditions Eq. (1) becomes unstable to smallnonaxisymmetric perturbations.It is convenient to deﬁne dimensionless variables byscaling the distance, time, velocities, gravitational potential, and energy with respect to
r
0
¼
r
m
,
t
0
¼
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
r
dm
=GM
p
,
v
0
¼
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
GM=r
d
2
m
p
,
c
0
¼
GM=r
d
2
m
, and
E
0
¼
GM
2
=r
d
2
m
,respectively. This is equivalent to setting
r
m
¼
G
¼
M
¼
1
.PRL
111,
230603 (2013) PHYSICAL REVIEW LETTERS
week ending6 DECEMBER 2013
00319007
=
13
=
111(23)
=
230603(5) 2306031
2013 American Physical Society
The particle dynamics is governed by Newton’s equationsof motion
€
r
¼r
c
ð
r
;t
Þ
;
(2)where the dot stands for the time derivative and
r
¼
P
i
x
i
^
e
i
,
i
¼
1
;
...
;d
, is the particle position. In the thermodynamic limit, the correlations between the particles canbe ignored, so that the force acting on a particle located at
r
is
F
¼r
c
ð
r
;t
Þ
, where
c
ð
r
;t
Þ
is the mean gravitationalpotential which satisﬁes the Poisson equation
r
2
c
¼
C
d
n
ð
r
;t
Þ
;
(3)where
n
ð
r
;t
Þ
is the particle number density.We deﬁne the ‘‘envelope’’ of the position and velocityparticle distributions to be
X
i
ð
t
Þ¼
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ð
d
þ
2
Þh
x
2
i
i
q
and
V
i
ð
t
Þ¼
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ð
d
þ
2
Þh
v
2
i
i
q
, respectively. The
hi
correspondsto the average over all the particles. Note that in thereduced units,
X
i
ð
0
Þ¼
1
and
V
i
ð
0
Þ¼
v
m
for all
i
, but asthe dynamics evolves, it is possible for the symmetrybetween the different directions to become broken. Ourgoal is to determine the equations of evolution for
X
i
ð
t
Þ
[13]. Taking two time derivatives of
X
2
i
ð
t
Þ
and one of
V
2
i
ð
t
Þ
and using the equations of motion, Eq. (2), we obtain
_
X
2
i
þ
X
i
€
X
i
¼
V
2
i
ð
d
þ
2
Þ
x
i
@
c
@x
i
(4)and
V
i
_
V
i
¼ð
d
þ
2
Þ
_
x
i
@
c
@x
i
:
(5)To calculate the averages appearing in Eqs. (4) and (5),
we need to know the meangravitational potential. Wesuppose that the srcinally spherically symmetric homogeneous distribution can become distorted into an ellipsoidalshape with the semiaxis
f
X
i
g
and uniform density
n
ð
r
;t
Þ¼
d=C
d
Q
i
X
i
ð
t
Þ
.Using the ellipsoidal coordinate system [14],the gravitational ﬁeld inside a
d
dimensional ellipsoid withthe semiaxis
f
X
i
g
can be calculated explicitly to be
@
c
@x
i
¼
d
2
x
i
g
i
ð
X
1
;
...
;X
d
Þ
;
(6)where
g
i
ð
X
1
;
...
;X
d
Þ¼
Z
1
0
ds
ð
X
2
i
þ
s
Þ
Q
d j
¼
1
ð
X
2
j
þ
s
Þ
1
=
2
:
(7)Furthermore, for a
d
dimensional ellipsoid with a uniformmass distribution, it can be shown that
h
x
2
i
i¼
X
2
i
=
ð
d
þ
2
Þ
.Substituting these results in Eqs. (4) and (5), we obtain a
closed set of coupled equations:
_
X
2
i
þ
X
i
€
X
i
¼
V
2
i
d
2
X
2
i
g
i
ð
X
1
;
...
;X
d
Þ
(8)and
V
i
_
V
i
¼
d
2
X
i
_
X
i
g
i
ð
X
1
;
...
;X
d
Þ
:
(9)We deﬁne the ‘‘emittance’’ [15] in the
i
th direction as
2
i
ð
t
Þð
d
þ
2
Þ
2
½h
x
2
i
ih
_
x
2
i
ih
x
i
_
x
i
i
2
¼
X
2
i
V
2
i
_
X
2
i
X
2
i
. Takinga time derivative of
2
i
ð
t
Þ
and using Eqs. (8) and (9), it is
possible to show that the
i
ð
t
Þ
are the constants of motion,
i
ð
t
Þ¼
i
ð
0
Þ
i
. Using this observation, the set of Eqs. (8) and (9) reduces to
€
X
i
¼
2
i
X
3
i
d
2
X
i
g
i
ð
X
1
;
...
;X
d
Þ
:
(10)For the initial waterbag distribution, Eq. (1),
2
i
ð
0
Þ¼
v
2
m
.The virial theorem requires that a
stationary
gravitational system in
d
dimensions must have
2
K
¼ð
2
d
Þ
U
,where
K
and
U
are the total kinetic and potential energies,respectively. For the initial waterbag distribution,
K
¼
v
2
m
d=
½
2
ð
d
þ
2
Þ
and the potential energy is
U
¼
d=
½ð
2
d
Þð
d
þ
2
Þ
, so that the virial condition reduces to
v
m
¼
1
. Although the initial waterbag distribution is not astationary solution of the collisionless Boltzmann (Vlasov)equation, we expect that if the virial condition is satisﬁed,the system will not exhibit strong envelope oscillations.This is indeed what has been observed for gravitationalsystems in
d
¼
1
, 2, and 3 [16–19]. On the other hand, if
the initial distribution does not satisfy the virial condition,the particle distribution will undergo violent oscillationswhich will lead to QSS with a corehalo structure [16,18].
To measure how strongly the initial distribution deviatesfrom the virial condition, we deﬁne a viral number
R
0
ð
2
K=
ð
ð
2
d
Þ
U
Þ
Þ¼
v
2
m
. With this deﬁnition the emittancebecomes
2
i
ð
t
Þ¼
R
0
.Let us ﬁrst consider a uniform spherically symmetricmass distribution of radius
R
ð
t
Þ
, i.e.,
X
i
ð
t
Þ¼
R
ð
t
Þ
for
i
¼
1
;
...
;d
. In this case the integral in Eq. (7) can beevaluated analytically to yield
g
i
¼
2
R
d
=d
, and the equation of evolution for the radius of the sphere becomes
€
R
¼
R
0
R
3
1
R
d
1
;
(11)with
R
ð
0
Þ¼
1
and
_
R
ð
0
Þ¼
0
. We see that, in agreementwith the earlier discussion, if the initial distribution satisﬁes the virial condition,
R
0
¼
1
, the sphere’s radiusremains constant for all time,
R
ð
t
Þ¼
1
for any
d
. For
d
3
, this equilibrium is stable because a small deviationfrom
R
0
¼
1
will result in small periodic oscillations of
R
.On the other hand, for
d
4
the equilibrium is unstable,and any
R
0
1
will lead to either collapse or anunbounded expansion of the particle distribution. Theseconclusions are in agreement with the old observation of Ehrenfest, who ﬁrst noted that there are no stable orbits forNewtonian gravity in
d
4
[20].To investigate the possible symmetry breaking of aninitially spherically symmetric mass distribution, weneed, therefore, to only consider
d
¼
2
and 3. For
d
¼
2
,the integral in Eq. (7) can be performed analyticallyPRL
111,
230603 (2013) PHYSICAL REVIEW LETTERS
week ending6 DECEMBER 2013
2306032
yielding
g
i
ð
X
1
;X
2
Þ¼
2
=X
i
ð
X
1
þ
X
2
Þ
. Equation (10) thensimpliﬁes to
€
X
i
¼
i
X
3
i
2
X
1
þ
X
2
; i
¼
1
;
2
:
(12)The symmetry breaking occurs if an initially vanishinglysmall ﬂuctuation grows as a function of time. To study thisinstability, it is convenient to introduce new variables,
X
i
ð
t
Þ¼
X
ð
t
Þþ
i
ð
t
Þ
;
(13)where
X
¼ð
P
i
X
i
Þ
=d
is the average of
X
i
’s and
i
is theasymmetry along the
i
th direction. Clearly
i
’s are relatedby
P
i
i
¼
0
. Hence, for
d
¼
2
there is only one independent asymmetry variable
¼
1
¼
2
. To locate theregion of instability, we perform a linear stability analysisof Eq. (12). Noting that
ð
21
22
Þ
O
ð
Þ
, to leading orderin
, Eq. (12) simpliﬁes to
€
þ
3
ð
21
þ
22
Þ
2
X
4
ð
t
Þ
¼ð
21
22
Þ
2
X
3
ð
t
Þ
;
(14)while the dynamics of
X
ð
t
Þ
to this order is
€
X
¼
21
þ
22
2
X
3
1
X :
(15)The dynamics of
is driven by the oscillations of
X
ð
t
Þ
. Inparticular, if the virial condition is satisﬁed and
21
¼
22
¼
R
0
¼
1
, the (
¼
0
,
_
¼
0
) is a stable ﬁxed point of Eq. (14). Therefore, if
R
0
1
, for small initial asymmetry,
ð
t
Þ
will not grow in time. However, if the initialdistribution does not satisfy the virial condition,
X
ð
t
Þ
willoscillate and may drive a parametric resonance which canmake
ð
t
Þ
unstable. This is precisely what is observed innumerical integration of Eqs. (14) and (15). We ﬁnd that
for sufﬁciently small (or large)
R
0
, the amplitude of
ð
t
Þ
oscillations grows without a bound. Note that in Eq. (14)the instability occurs as a consequence of a ﬂuctuationeither in the velocity [
ð
0
Þ¼
0
and
1
2
], the position[
ð
0
Þ
0
], or as a combination of both. For sufﬁcientlysmall (or large)
R
0
, we ﬁnd that any small ﬂuctuation inthe initial particle distribution is ampliﬁed by the dynamics. Of course, in practice the growth of
ð
t
Þ
will besaturated by the Landau damping [16,18] and will result
in a QSS with a broken rotational symmetry.Topreciselylocatetheinstabilitythreshold,itissimplestto consider a small ﬂuctuation with
ð
0
Þ
0
and
1
¼
2
.Since the
ð
t
Þ
is driven by the periodic oscillations of
X
ð
t
Þ
,to study this instability we must work in the Poincare´section [21,22].
Consider a displacement vector from the (
¼
0
,
_
¼
0
) ﬁxed point,
X
ð
t
Þ¼ð
;
_
Þ
. From Eq. (14), we seethat its dynamics is governed by
_
X
¼
M
X
, where
M
¼
0 1
3
R
0
X
4
0
!
;
(16)and the dynamics of
X
ð
t
Þ
is given by Eq. (15) with
21
¼
22
¼
R
0
. If we now deﬁne a mapping
M
ð
t
Þ
that relates
X
ð
t
Þ
to its initial condition by
X
ð
t
Þ¼
M
ð
t
Þ
X
ð
0
Þ
,and substitute this into the evolution equation for
X
, weobtain
_
M
¼
M
M
;
with
M
ð
0
Þ¼
I;
(17)where
I
is the identity matrix. In order to determine thestability of (
¼
0
,
_
¼
0
) ﬁxed point, we simultaneouslyintegrate Eqs. (15) and (17) over one period
R
of theoscillation of
X
ð
t
Þ
(i.e., between two consecutive pointsin the Poincare´ map), and determine the eigenvalues of themapping matrix
M
ð
R
Þ
. If the absolutevalue of any eigenvalue is larger than 1, then (
¼
0
,
_
¼
0
) ﬁxed point willbe unstable. We ﬁnd that the asymmetric instability occursfor
R
0
<
0
:
255893...
and for
R
0
>
2
:
55819...
. A moredetailed analysis shows that it is produced by a pitchfork bifurcationand is of second order. In Fig. 1 we compare thepredictions of the theory with the results of extensivemolecular dynamics simulations performed using thestateoftheart gravitational oriented massively parallel
GADGET2
code[23],whichhasbeenappropriatelymodiﬁedto integrate gravity in two dimensions. At
t
¼
0
the particles are distributed in accordancewith Eq. (1). To force thesymmetry breaking to occur along the
x
axis, a smallperturbation in this direction is introduced. We then monitor the moments
h
x
2
i
and
h
y
2
i
as the dynamics evolves.Figure 1 shows the evolution of the moments for twodifferent virial numbers. We ﬁnd that for
R
0
¼
0
:
16
thesymmetry is broken, while for
R
0
¼
0
:
36
the sphericalsymmetryisunaffectedbytheinitialperturbation.Thisisinclose agreement with the predictions of the present theory.A similar symmetry breaking transition is also found for
0.1 0.2 0.3 0.4 0.5 0 100 200 300 400
(a)
0 0.5 1 1.5 2 0 25 50 75 100
0.1 0.2 0.3 0.4 0.5 0 100 200 300 400
(b)
0 0.5 1 1.5 2 2.5 3 0 25 50 75 100
FIG. 1 (color online). The evolution of
x
and
y
moments of themass distribution,
h
x
2
i
(red solid curve) and
h
y
2
i
(blue dashedcurve) obtained using molecular dynamics simulations for a 2Dsystem with
N
¼
8000
. A small asymmetry in the
x
direction isintroduced in the initial particle distribution. For initial distribution with
R
0
¼
0
:
36
(a), the system relaxes to a QSS with aspherical symmetry (see also Fig. 2), while for
R
0
¼
0
:
16
(b),spherical symmetry is broken. Similar behavior is found for
R
0
above the upper critical threshold; see Fig. 2. The inset shows theevolution of the virial number. Both the symmetric and theasymmetric QSS are fully virialized,
R
¼
1
.
PRL
111,
230603 (2013) PHYSICAL REVIEW LETTERS
week ending6 DECEMBER 2013
2306033
large virial numbers; see Fig. 2. Since the transitions arecontinuous, it is difﬁcult to precisely locate the thresholdsof instability using molecular dynamics simulations.In Fig. 2 we show snapshots of two QSS to which thesystem relaxes after a few oscillations. In agreement withthe theory, depending on the virial number, one of the QSSis spherically symmetric while the other one is not. For
d
¼
3
the integral in Eq. (7) cannot be performed in termsof simple analytical functions, and must be evaluatednumerically. To locate the instability, we once againmake use of the variables deﬁned in Eq. (13) and expandEq. (10) to linear order in
i
. For
d
¼
3
, there are twoindependent variables
1
and
2
. Numerical integration of these equations shows, once again, the existence of aninstability for small and large virial numbers. To preciselylocate the instability, we ﬁx
21
¼
22
¼
23
¼
R
0
. To linearorder the dynamics of equations for
1
and
2
thendecouples and becomes identical. This means that we canstudy the stability using a single
ð
t
Þ
variable. The matrixthat determines the evolution of the displacement vectorfrom (
¼
0
,
_
¼
0
) ﬁxed point now takes the form
M
¼
0 1
R
15
R
0
5
R
4
0
!
;
(18)where
R
ð
t
Þ
is given by Eq. (11) with
d
¼
3
. Substitutingthis matrix in Eq. (17) and adopting the procedure analogous to the one used before, we ﬁnd that the ﬁxed point(
¼
0
,
_
¼
0
) becomes unstable for
R
0
<
0
:
388666...
and
R
0
>
1
:
61133...
. Figure 3 shows two snapshots of the evolution of a 3D gravitational systems. As predictedby the theory, for both small and large virial numbers thespherical symmetry of the initial distribution is broken bythe parametric resonances.For 3D systems ﬁnite angular momentum can also leadto breaking of the spherical symmetry. This, however, isnot the case in 2D. Furthermore, in our simulations theinitial particle distribution has very small angularmomentum—in the thermodynamic limit it will be exactlyzero. The rotation of the system is, therefore, very slow,while the instabilityhappens veryquickly, showing that theresidual angular momentum does not play any role for thesymmetry breaking studied in this Letter.It is interesting to compare and contrast the Rayleighinstability of charged conducting droplets and the instabilityof selfgravitating systems. While the Rayleigh instability isa true thermodynamictransition, the gravitationalsymmetrybreaking is not. When the charge on a droplet exceeds thecritical value
Q
c
, it will undergo a ﬁrstorder transition to aprolate ellipsoid. On the other hand, the instability of a selfgravitating system is a purely dynamical phenomenon, arising from a parametric resonance that drives an asymmetricmode of oscillation. The magnitude of the instability issaturatedbythenonlinearLandaudamping[24]whichleadsto the formation of a nonequilibrium corehalo QSS. If theinstability occurs, the broken ergodicity [11] prevents thesymmetry from being restored. In
d
¼
2
, a selfgravitatingsystemwithaﬁnitenumberofparticleswilleventuallyrelaxto thermodynamic equilibrium in which the distributionfunction will have the usual BoltzmannGibbs form [18]and the meangravitational potential will once again bespherically symmetric. The relaxation time to equilibrium,however, diverges with
N
, so that in practice a sufﬁcientlylarge system (such as an elliptical galaxy) will never evolveto equilibrium, but will stay in a nonequilibrium stationarystate forever [25]. For such systems, once the instabilityoccurs, the symmetry will remain irrevocably broken.This work was partially supported by the CNPq,FAPERGS, INCTFCx, and by the U.S.AFOSR underGrant No. FA95501210438. Numerical simulationshave been performed at the cluster of the SIGAMMhosted at ‘‘Observatoire de Coˆte d’Azur,’’ Universite´ deNice–Sophia Antipolis.
[1] L. Rayleigh, Philos. Mag.
14
, 184 (1882).[2] G.I. Taylor, Proc. R. Soc. A
280
, 383 (1964).[3] Y. Kantor and M. Kardar, Phys. Rev. E
51
, 1299 (1995).[4] L. Aguilar and D. Merritt, Astrophys. J.
354
, 33 (1990).
42 0 2 442 0 2 4
(a)
4020 0 20 404020 0 20 40
(b)
FIG. 2. Snapshots of the
x

y
particle distribution in a QSS at
t
¼
200
. In (a)
R
0
¼
2
and the symmetry remains unbroken,while in (b)
R
0
¼
6
:
25
and the QSS has a broken rotationalsymmetry. Note that the ﬁnal particle distribution in both caseshas a characteristic corehalo structure.
105 0 5 10105 0 5 10
(a)
10050 0 50 10010050 0 50 100
(b)
FIG. 3. Snapshots of the
x

y
projection of the 3D particledistribution (
N
¼
20000
) at
t
¼
25
. In (a)
R
0
¼
0
:
5
, the symmetry remains unbroken, while in (b)
R
0
¼
0
:
01
there is aspontaneous symmetry breaking. Note that because of the particle evaporation a 3D system does not relax to a QSS.
PRL
111,
230603 (2013) PHYSICAL REVIEW LETTERS
week ending6 DECEMBER 2013
2306034