# Symmetry Breaking in d-Dimensional Self-gravitating Systems

###### Abstract

Systems with long-range interactions, such as self-gravitating clusters and magnetically confined plasmas, do not relax to the usual Boltzmann-Gibbs thermodynamic equilibrium, but become trapped in quasi-stationary states (QSS) the life time of which diverges with the number of particles. The QSS are characterized by the lack of ergodicity which can result in a symmetry broken QSS starting from a spherically symmetric particle distribution. We will present a theory which allows us to quantitatively predict the instability threshold for spontaneous symmetry breaking for a class of -dimensional self-gravitating systems.

###### pacs:

05.20.-y, 05.70.Ln, 05.45.-aLord Rayleigh was probably the first to make an observation that long-range forces can lead to symmetry breaking Re82 (). Rayleigh was studying stability of conducting spherical fluid droplets carrying charge . He discovered that when exceeds a certain critical threshold , droplets becomes unstable to symmetry breaking perturbations, elongating and eventually breaking up, emitting jets of fluid that carry away a significant fraction of the charge Ta64 (). Rayleigh instability is now the basis for various technological applications, such as electrospraying and electrospinning. It also helps to understand the conformational structure of charged polymers, such as polyampholytes KaKa95 (). For self-gravitating systems a similar instability has been observed in gravitational simulations AgMe90 (). It has been found that an initially spherically symmetric self-gravitating system can become unstable, leading to formation of structures of reduced symmetry AgMe90 (). This radial orbit instability is believed to be important for the formation of elliptical galaxies AtMa13 ().

There is, however, a fundamental difference between the Rayleigh instability of charged spherical droplets and the instability of spherically symmetric self-gravitating systems. Since the droplets are in (canonical) thermodynamic equilibrium, their shape must correspond to the minimum of the Helmholtz free energy — in fact, even for somewhat below a spherical shape is already metastable, with the global minimum corresponding to a strongly prolate ellipsoid AiGa62 (). The thermal fluctuations, however, are too small to overcome the barrier that separates the metastable minimum from the global one, so that the spherical shape persists up to the Rayleigh threshold. On the other hand, gravitational systems are intrinsically microcanonical — isolated from environment Gr01 (); IsCo01 (); ChIs02 (). In the thermodynamic limit, such long-range systems do not evolve to thermodynamic equilibrium, but become trapped in quasi-stationary states (QSS), the life time of which diverges with the number of particles Campa09 (). The QSS are characterized by the broken ergodicity, making equilibrium statistical mechanics inapplicable BeTe12 (). To explore spontaneous symmetry breaking of systems with long-range forces, therefore, requires a completely different approach TeBe12 (). In this Letter we will present a theory which allows us to quantitatively predict the thresholds of symmetry breaking instabilities for systems with long-range interactions. The results of the theory will be compared with extensive molecular dynamics simulations.

To present the theory, we will study a class of self-gravitating systems of particles of mass in an infinite -dimensional space. The interaction potential between the particles is , where is the gravitational constant. We will work in thermodynamic limit, and , while the total mass remains fixed. The initial particle distribution is assumed to be a uniform spherically symmetric water-bag in both configuration and velocity space,

(1) |

where is the Heavyside step function and is the surface area of a -dimensional unit sphere, and is the gamma function. Since the initial water-bag distribution is not a stationary solution of the collisionless Boltzmann (Vlasov) equation the systems will evolve with time. We are interested to discover under what conditions Eq.(1) becomes unstable to small non-axisymmetric perturbations.

It is convenient to define dimensionless variables by scaling the distance, time, velocities, the gravitational potential, and the energy with respect to: , , , and , respectively. This is equivalent to setting . The particle dynamics is governed by Newton’s equations of motion

(2) |

where the dot stands for the time derivative and , , is the particle position. In thermodynamic limit, the correlations between the particles can be ignored, so that the force acting on a particle located at is , where is the mean gravitational potential which satisfies the Poisson equation

(3) |

where is the particle number density.

We define the ”envelope” of the position and velocity particle distributions to be and , respectively. The corresponds to the average over all the particles. Note that in the reduced units, and for all , but as the dynamics evolves, it is possible for the symmetry between the different directions to become broken. Our goal is to determine the equations of evolution for Sa71 (). Taking two time derivatives of and one of and using the equations of motion, Eq. (2), we obtain

(4) |

and

(5) |

To calculate the averages appearing in Eqs. (4) and (5), we need to know the mean-gravitational potential. We suppose that the originally spherically symmetric homogeneous distribution can become distorted into an ellipsoidal shape with the semi-axis and uniform density . Using the ellipsoidal coordinate system Ke53 () the gravitational field inside a -dimensional ellipsoid with the semi-axis can be calculated explicitly to be,

(6) |

where

(7) |

Furthermore, for a -dimensional ellipsoid with a uniform mass distribution it can be shows that . Substituting these results in Eqs. (4) and (5), we obtain a closed set of coupled equations

(8) |

and

(9) |

We define the ”emittance” in the ’th direction as =. Taking a time derivative of and using the equations (8) and (9), it is possible show that the are the constants of motion, . Using this observation the set of Eqs (8) and (9) reduces to

(10) |

For the initial water-bag distribution, Eq. (1), .

The virial theorem requires that a stationary gravitational system in dimensions must have , where and are the total kinetic and potential energies, respectively. For the initial water-bag distribution and the potential energy is , so that the virial condition reduces to . Although the initial water-bag distribution is not a stationary solution of the collisionless Boltzmann (Vlasov) equation, we expect that if the virial condition is satisfied, the system will not exhibit strong envelope oscillations. This is indeed what has been observed for gravitational systems in and TeLe11 (); JoWo11 (); TeLe10 (); LePa08a (). On the other hand if the initial distribution does not satisfy the virial condition, the particle distribution will undergo violent oscillations which will lead to QSS with a core-halo structure TeLe11 (); TeLe10 (). To measure how strongly the initial distribution deviates from the virial condition, we define a viral number . With this definition the emittance becomes .

Let us first consider a uniform spherically symmetric mass distribution of radius , i.e., for . In this case the integral in Eq. (7) can be evaluated analytically to yield , and the equation of evolution for the radius of the sphere becomes

(11) |

with and . We see that in agreement with the earlier discussion, if the initial distribution satisfies the virial condition, , the sphere’s radius remains constant for all time, for any . For , this equilibrium is stable because a small deviation from will result in small periodic oscillations of . On the other hand, for the equilibrium is unstable, and any will lead to either collapse or an unbounded expansion of the particle distribution. These conclusions are in agreement with the old observation of Paul Ehrenfest, who first noted that there are no stable orbits for Newtonian gravity in Eh20 ().

To investigate the possible symmetry breaking of an initially spherically symmetric mass distribution we need, therefore, to only consider and . For , the integral in Eq. (7) can be performed analytically yielding . Eq. (10) then simplifies to

(12) |

The symmetry breaking occurs if an initially vanishingly small fluctuation grows as a function of time. To study this instability, it is convenient to introduce new variables

(13) |

where is the average of ’s and is the asymmetry along the th direction. Clearly ’s are related by . Hence, for there is only one independent asymmetry variable . To locate the region of instability, we perform a linear stability analysis of Eqs. (12). Noting that , to leading order in , Eq. (12) simplifies to

(14) |

while the dynamics of to this order is

(15) |

The dynamics of is driven by the oscillations of . In particular, if the virial condition is satisfied and , the is a stable fixed point of Eq. (14). Therefore if , for small initial asymmetry, will not grow in time. However, if the initial distribution does not satisfy the virial condition, will oscillate and may drive a parametric resonance which can make unstable. This is precisely what is observed in numerical integration of Eqs. (14) and (15). We find that for sufficiently small (or large) , the amplitude of oscillations grows without a bound. Note that in Eq. (14) the instability occurs as a consequence of a fluctuation either in the velocity ( and ), the position (), or as a combination of both. For sufficiently small (or large) , we find that any small fluctuation in the initial particle distribution is amplified by the dynamics. Of course, in practice the growth of will be saturated by the Landau damping TeLe11 (); TeLe10 () and will result in a QSS with a broken rotational symmetry.

To precisely locate the instability threshold it is simplest to consider a small fluctuation with and . Since the is driven by the periodic oscillations of , to study this instability we must work in the Poincaré section PoHi93 (); SiRe06 ().

Consider a displacement vector from the fixed point, . From Eq. (14), we see that its dynamics is governed by , where

(16) |

and the dynamics of is given by Eq. (15) with . If we now define a mapping that relates to its initial condition by , and substitute this into the evolution equation for , we obtain

(17) |

where is the identity matrix. In order to determine the stability of fixed point, we simultaneously integrate Eqs. (11) and (17) over one period of the oscillation of (i.e., between two consecutive points in the Poincaré map), and determine the eigenvalues of the mapping matrix . If the absolute value of any eigenvalue is larger than 1, then fixed point will be unstable. We find that the asymmetric instability occurs for and for . A more detailed analysis shows that it is produced by a pitchfork bifurcation and is of second order. In Fig. 1 we compare the predictions of the theory with the results of extensive molecular dynamics simulations performed using the state-of-the-art gravitational oriented massively parallel GADGET2 code Gad (), which has been appropriately modified to integrate gravity in two dimensions. At the particles are distributed in accordance with Eq. (1). To force the symmetry breaking to occur along the -axis, a small perturbation in this direction is introduced. We then monitor the moments and as the dynamics evolves. Fig. 1 shows the evolution of the moments for two different virial numbers. We find that for the symmetry is broken while for the spherical symmetry is unaffected by the initial perturbation. This is in close agreement with the predictions of the present theory. Similar symmetry breaking transition is also found for large virial numbers, see Fig. 2. Since the transitions are continuous, it is difficult to precisely locate the thresholds of instability using molecular dynamics simulations.

In Fig. 2 we show the snapshots of two QSS to which the system relaxes after a few oscillations. In agreement with the theory, depending on the virial number, one of the QSS is spherically symmetric while the other one is not.

For the integral in Eq. (7) cannot be performed in terms of simple analytical functions, and must be evaluated numerically. To locate the instability, we once again make use of the variables defined in Eq. (13) and expand Eq. (10) to linear order in . For , there are two independent variables and . Numerical integration of these equation shows, once again, existence of an instability for small and large virial numbers. To precisely locate the instability we fix . To linear order the dynamics of equations for and then decouples and becomes identical. This means that we can study the stability using a single variable. The matrix that determines the evolution of the displacement vector from fixed point now takes the form

(18) |

where is given by Eq. (11) with . Substituting this matrix in Eq. (17) and adopting the procedure analogous to the one used before, we find that the fixed point , becomes unstable for and . Fig. 3 shows two snapshots of the evolution of a 3D gravitational systems. As predicted by the theory, both for small and large virial numbers the spherical symmetry of the initial distribution is broken by the parametric resonances.

For 3d systems finite angular momentum can also lead to breaking of the spherical symmetry. This, however, is not the case in 2d. Furthermore, in our simulations the initial particle distribution has very small angular momentum — in the thermodynamic limit it will be exactly zero. The rotation of the system is, therefore, very slow, while the instability happens very quickly, showing that the residual angular momentum does not play any role for the symmetry breaking studied in this paper.

It is interesting to compare and contrast the Rayleigh instability of charged conducting droplets and the instability of self-gravitating systems. While the Rayleigh instability is a true thermodynamic transition, the gravitational symmetry breaking is not. When the charge on a droplet exceeds the critical value , it will undergo a first order transition to a prolate ellipsoid. On the other hand, the instability of a self-gravitating system is a purely dynamical phenomenon, arising from a parametric resonance that drives an asymmetric mode of oscillation. The magnitude of the instability is saturated by the non-linear Landau damping LePa08 () which leads to the formation of a non-equilibrium core-halo QSS. If the instability occurs, the broken ergodicity BeTe12 () prevents the symmetry from being restored. In , a self-gravitating system with finite number of particles will eventually relax to thermodynamic equilibrium in which the distribution function will have the usual Boltzmann-Gibbs form TeLe10 () and the mean-gravitational potential will once again be spherically symmetric. The relaxation time to equilibrium, however, diverges with , so that in practice a sufficiently large system (such as an elliptical galaxy) will never evolve to equilibrium, but will stay in a non-equilibrium stationary state forever Ma13 (). For such systems once the instability occurs, the symmetry will remain irrevocably broken. This work was partially supported by the CNPq, FAPERGS, INCT-FCx, and by the US-AFOSR under the grant FA9550-12-1-0438. Numerical simulations have been performed at the cluster of the SIGAMM hosted at “Observatoire de Côte d’Azur”, Université de Nice – Sophia Antipolis.

## References

- (1) Lord Rayleigh, Phil. Mag. 14, 184 (1882).
- (2) G. I. Taylor, Proc. R. Soc. London Ser. A 280, 383 (1964).
- (3) Y. Kantor and M. Kardar, Phys. Rev. E 51, 1299 (1995).
- (4) L. Aguilar and D. Merritt, Astrophysical Journal, 354, 33-51 (1990).
- (5) E. Athanassoula, R. E. G. Machado and S. A. Rodionov, Mon. Not. R. Astron. Soc. 429, 1949-1969 (2013).
- (6) G. Ailam and I. Gallily, Phys. Fluids 5, 575 (1962).
- (7) D. H. E. Gross Microcanonical Thermodynamics: Phase transitions in Small Systems (Lecture Notes in Phys., vol. 66 World Scientific, Singapore (2001)).
- (8) I. Ispolatov and E. G. D. Cohen, Phys. Rev. E 64, 056103 (2001).
- (9) P. H. Chavanis and I. Ispolatov Phys. Rev. E 66, 036109 (2002).
- (10) A. Campa, T. Dauxois, and S. Ruffo, Phys. Rep. 480, 57 (2009); Y. Levin, R. Pakter, F.B. Rizzato, T. N. Teles, and F. P. da C. Benetti, Phys. Rep. (2013) http://dx.doi.org/10.1016/j.physrep.2013.10.001.
- (11) F. P. da C. Benetti, T. N. Teles, R. Pakter, and Y. Levin, Phys. Rev. Lett. 108, 140601 (2012).
- (12) T. N. Teles, F. P. da C. Benetti, R. Pakter, and Y. Levin, Phys. Rev. Lett. 109, 230601 (2012).
- (13) F. J. Sacherer, IEEE Transactions on Nuclear Science 18 1105 (1971).
- (14) M. Reiser, Theory and Design of Charged Particle Beams (Wiley-InterScience, New York, 1994).
- (15) T. N. Teles, and Y. Levin, and R. Pakter, Mon. Not. R. Astron. Soc. 417, L21, (2011).
- (16) M. Joyce and T. Worrakitpoonpon, Phys. Rev. E 84, 011139 (2011).
- (17) T.N. Teles, Y.Levin, R. Pakter, and F.B. Rizzato, J. Stat. Mech. P05007 (2010).
- (18) Y. Levin, R. Pakter, and F.B. Rizzato, Phys. Rev. E 78, 021130 (2008).
- (19) O. D. Kellogg, Foundations of Potential Theory, (Dover, New York, 1953).
- (20) P. Ehrensfest, Annalen der Physik 61, 440 (1920).
- (21) C. Polymilis and K. Hizanidis , Phys. Rev. E 47, 4381 (1993).
- (22) W. Simeoni Jr., F. B. Rizzato, and R. Pakter, Phys. Plasmas 13, 063104 (2006).
- (23) V. Springel, Mon. Not. R. Astron. Soc. 364, 1105 (2005).
- (24) Y. Levin, R. Pakter and T. N. Telles, Phys. Rev. Lett. 100, 040604 (2008).
- (25) B. Marcos, Phys. Rev. E, 88, 032112 (2013).