# 15. Instabilities in Spheroidal Systems

## Astronomy 626: Spring 1997

Galaxies should be stable equilibrium solutions of the CBE. Spheroidal, pressure-supported systems with highly anisotropic velocity distributions are often unstable, relaxing to less flattened and less anisotropic configurations in a few dynamical times. Such instabilities may explain the absence of elliptical galaxies flatter than E6 or E7.

No all equilibria are stable. One example is a pencil balanced on end; this is equilibrium position, but the slightest perturbation will cause it to fall over. The balanced pencil is unstable. A second example is a `public-address' system, consisting of a microphone, an amplifier, and a loud-speaker; a signal going around this loop can grow in amplitude with each cycle, producing ear-splitting feedback. The PA system is overstable.

Instabilities also occur in stellar systems. Examples of potentially unstable stellar systems include:

• Homogeneous systems (Jeans instability);
• Disk systems (ring & bar instabilities);
• Spherical systems (Henon, radial-orbit, & tangential-orbit instabilities);
• Axisymmetric systems (radial-orbit, ring, off-center, & bending instabilities).
The Jeans instability plays a key role in structure formation in cosmology. Instabilities in rotating disk systems play important roles in constraining and determining the structure of disk galaxies. This lecture will focus on instabilities in pressure-supported spherical and axisymmetric systems. These instabilities may provide important constraints on the structure of elliptical galaxies.

## 15.1 Methods of Analysis

### Linear Stability Analysis

Suppose that f_0(r,v) is a time-independent solution of the CBE, and that Phi_0(r) is the self-consistent potential generated by this solution. To determine if this solution is stable, consider perturbed solutions of the form

```(1)         f(r,v,t) = f_0(r,v) + epsilon f_1(r,v,t)
```
where epsilon is a small parameter and f_1(r,v,t) is an arbitrary function. Substitute this perturbed distribution function into the time-dependent CBE, and keep only terms of first order in epsilon; the result is the linearized CBE,
```            d f_1       d f_1                d f_1                d f_0
(2)         ----- + v . ----- - grad Phi_0 . ----- - grad Phi_1 . ----- = 0 .
d t         d r                  d v                  d v
```
Here the first three terms describe the evolution of the perturbed distribution function in the unperturbed field Phi_0, while the last term describes the effect of the perturbation in the potential, Phi_1(r,t), on the unperturbed distribution. Here Phi_1 is given by Poisson's equation,
```                                    /
(3)         div grad Phi_1 = 4 pi G | dv f_1(r,v,t) .
/
```
Eqs. 2 and 3 together describe the evolution of the perturbed system under the assumption that the perturbations are small.

To discover if an instability exists -- or to prove that the systems is stable -- one must consider a complete set of perturbations f_1(r,v,t). The choice of the set of basis functions depends on the geometry of the system. For homogeneous systems, a Fourier expansion is appropriate; this approach is useful in deriving the Jeans instability.

Having devised a complete set of perturbations, the trick is to find linear combinations which grow. A linear combination can be represented as a vector, and time evolution as a matrix acting on this vector. Then combinations which grow are eigenvectors of this matrix. For example, if the growing solutions depend on time like

```                   - i omega t
(4)         f_1 ~ e            ,
```
with omega = i gamma where gamma is a real number, then the perturbation grows exponentially and the system is unstable; if omega = omega_R + i gamma, where omega_R is also real, then the perturbation oscillates as it grows and the system is overstable.

### Numerical Simulations

Due to problems in chosing good basis sets, the linear analysis just described is a fairly daunting task. An alternate approach is to feed an N-body realization of the equilibrium system into an N-body code, and run the system to see if it remains in equilibrium. This approach may rely on particle noise to `seed' growing perturbations, although one also has the option of introducing perturbations `by hand'. It's worth emphasizing that such numerical simulations only reveal fairly gross and violent instabilities; slowly-growing perturbations may be lost in the fluctuations due to discreteness and consequently can't be detected.

Nonetheless, N-body simulations have played an important role in the investigation of dynamical instabilities in stellar systems. Among things, they have led to the

• discovery of new instabilities,
• treatment intractable cases, and
• study of non-linear evolution.

## 15.2 Spherical Systems

### Isotropic Systems: f = f(E)

Antonov (1960, 1962) derived a necessary and sufficient condition for the stability of spherical, isotropic systems. This criterion, expressed in terms of a complex variational principle, can distinguish both stable and unstable systems. Some simple consequences of this criterion are that

• any system with df/dE < 0 is stable to non-spherical perturbations, and
• any system with df/dE < 0 and d^3 rho/d Phi^3 <= 0 is stable to all perturbations.
These secondary criteria are sufficient but not necessary conditions for stability; systems which violate them may still be stable. It's unknown if any spherical, isotropic systems are actually unstable. The best candidate is the polytrope with index n = 1/2, in which all stars have the same binding energy; this system exhibits radial oscillations which seem larger than expected given the number of bodies used in the simulations (Henon 1973).

### Anisotropic Systems: f = f(E,J)

Instability in systems composed of stars on exactly radial orbits was predicted by Antonov (1973), and confirmed numerically by Polyachenko (1981), who showed that such systems rapidly evolve into elongated, bar-shaped configurations.

Henon (1973) systematically investigated the dynamical stability of a class of models known as generalized polytropes, which have the distribution function

```                                  n-3/2  2m
{ K (E_1 - E)      J   , if E <= E_1
(5)         f(E,J) = {
{ 0 ,                    otherwise .
```
These models have finite radius since f -> 0 at some energy E = E_1 < 0; the parameters K and E_1 together fix the total mass and radius of the system. The parameters n >= 1/2 and m >= -1 govern the structure of the system. The energy distribution is controlled by n, while m determines the velocity anisotropy; in particular, the ratio of radial to tangential velocity dispersion is
```            sigma_r^2     1
(6)         --------- = ----- ,
sigma_t^2   1 + m
```
so m = -1, 0, infinity correspond to radial, isotropic, and tangential systems, respectively.

Using a spherically-symmetric N-body code, Henon (1973) found a dynamical instability in generalized polytropes `when n is low and the velocity distribution is radially elongated'. This instability takes the form of radial oscillations which grow in amplitude. As the system pulsates, bodies are scattered in binding energy; eventually the distribution function changes enough to shut off the instability.

At P. Hut's suggestion, I repeated Henon's experiments with an N-body code which did not enforce spherical symmetry and thus allowed non-spherical perturbations to grow as well (Barnes 1985). Being unaware of Antonov and Polyachenko's discovery of a non-spherical instability in systems dominated by radial orbits, I was initially suspicious of numerical bugs when my simulations of generalized polytropes with preferentially radial orbits (m < 0) evolved from spheres into triaxial, bar-like configurations. Only after much testing with different N-body codes could I defend the claim that this instability was real.

I also found that generalized polytropes with preferentially tangential orbits (m > 0), which are unusual in having hollow centers, are unstable in a different way. These systems exhibit non-spherical oscillations which grow in amplitude; bodies scattered into orbits of low angular momentum eventually fill in the hollow center of the model and shut down the instability. Only systems with d rho/dr > 0 are subject to this instability, so it's probably not relevant for realistic galaxies. The limiting case in which m -> infinity and the system consists of a thin spherical shell of bodies on circular orbits was subject to a linear stability analysis by J. Goodman (Barnes, Goodman, & Hut 1986). In this limit, all stars have the same orbital period, so any initial non-spherical perturbation will be recreated at the antipodal point half an orbit period later; gravity attracts more bodies to over-dense perturbations, which consequently grow.

To map the range of these instabilities, we ran a 10 by 14 grid of models in the (m,n) plane and measured changes in mass profile and ellipticity (Barnes, Goodman, & Hut 1986). These showed that Henon's spherically-symmetric instability was confined largely to systems with m + n < 1/2. On the other hand, the radial-orbit and tangential-orbit instabilities were completely insensitive to the value of n and even the sign of df/dE.

The astrophysical relevance of the radial-orbit instability was emphasized when it was found in several more-or-less realistic models of spherical galaxies. Merritt & Aguilar (1985) showed that anisotropic Jaffe (1983) models with distribution functions of the form f = f(E + J^2 / 2r_a^2) are unstable if the anisotropy radius r_a < 0.3 a, where a is the scale radius. And Merritt (1987) showed that an anisotropic model of M87, which reproduced the observed velocity profiles without invoking a central black hole (Newton & Binney 1984), would evolve into an elongated configuration.

A crude picture of the radial-orbit instability is to imagine exactly radial orbits as rigid rods which are constrained to pass through the center but can freely pivot about that point; a spherical system of such rods can move to a state of lower potential energy if the rods clump together about some axis determined by a slight initial overdensity. This picture relates the radial-orbit instability to the Jeans instability for mass points constrained to stay on the surface of a sphere.

A better picture of this instability was put forward by Palmer & Papaloizou (1987). In a general spherical potential, a highly elongated orbit will precess at a slow and constant rate. If a weak bar-like potential is added, the orbit will gain angular momentum as it comes into alignment with the bar, and lose angular momentum as it continues to precess past the bar. If the orbit's initial angular momentum is low enough, it will be trapped by the bar and confined to a box orbit, adding its mass to the mass already comprising the bar; hence the bar will grow.

Palmer & Papaloizou (1987) also showed that all models with f(E,J) diverging like a power law in J as J -> 0 are unstable. This implies that no rigorous stability criterion can be formulated in terms of the ratio of radial to total kinetic energy; systems with f(E,J) diverging as J -> 0 exist with arbitrarily small amounts of radial anisotropy, and all are formally unstable. Palmer & Papaloizou's criterion is sufficient but not necessary for instability; as they and others showed, radially anisotropic systems with finite f(E,0) can also be unstable (e.g. Dejonghe & Merritt 1988).

## 15.4 Axisymmetric Systems

N-body methods are probably the only viable way to assess the stability of three-dimensional systems. The enormous range of possible axisymmetric models makes any definitive survey of instabilities in non-spherical systems a daunting task. Here I will simply mention some results obtained for axisymmetric systems.

The construction of axisymmetric equilibrium systems for stability testing is tricky since most orbits in axisymmetric potentials possess a third, non-classical integral of unknown form. One possible set of axisymmetric models which can be derived analytically are the `shell-orbit' models described by Bishop (1987). These are built out of tube orbits of zero radial thickness; they are thus somewhat unusual and may be more prone to instabilities than models using orbits of finite radial thickness. Another way to construct axisymmetric models is Schwarzschild's (1979) linear programming technique. This approach permits model builders to use the full range of possible orbits. Levison & Richstone (1985) have used this method to construct E6 models with flat rotation curves and a wide range of kinematic properties.

An analog of the radial-orbit instability is seen in axisymmetric models constructed using orbits with low angular momentum J_z (Palmer, Papaloizou, & Allen 1990, Levison, Duncan, & Smith 1990). Oblate and prolate systems composed of such orbits evolve into triaxial structures; this suggests that the natural endpoint of the radial-orbit instability is a triaxial configuration.

### Oblate Systems

Oblate Bishop models as flat as E6 or flatter are subject to an axisymmetric instability, with particles on adjacent orbits clumping together into thin, curving cylinders (Merritt & Stiavelli 1990). A similar axisymmetric instability has long been known in disk systems composed of stars on nearly circular orbits; such a disk will break up into a set of nested rings (Toomre 1964). It seems likely that this instability occurs because of the small radial velocity dispersion of these `shell-orbit' models.

Bishop models even as round as E1 are subject to a `m = 1' instability which displaces the density maximum from the system's center of mass (Merritt & Stiavelli 1990). This instability is not peculiar to Bishop models; it's also seen in Levison & Richstone models with low radial velocity dispersions (Levison, Duncan, & Smith 1990).

### Prolate Systems

Prolate Bishop models more elongated than E6 are subject to bending instabilities, temporarily becoming either banana or S-shaped before relaxing to less elongated configurations (Merritt & Hernquist 1991). This instability is probably related to the ``firehose'' instability in a thin, infinite sheet of stars (Toomre 1966). It seems that this instability is not peculiar to `shell-orbit' models but is likely to occur in any prolate system more elongated than about E7 (Merritt & Sellwood 1994).

## References

• Antonov, V.A. 1960, Astr. Zh. 37, 918 (Soviet Astron. 4, 859).
• Antonov, V.A. 1962, Vestnik Leningrad Univ. 19, 96.
• Antonov, V.A. 1973, in Dinamika Galaktik i Zvezdnykh Skoplenii (Alma-Ata: Nauka).
• Barnes, J. 1995, in Dynamics of Star Clusters, eds. J. Goodman & P. Hut, p. 297.
• Barnes, J., Goodman, J., & Hut, P. 1986, Ap. J. 300, 112.
• Bishop, J. 1987, Ap. J. 322, 618.
• Dejonghe, H. & Merritt, D. 1988, Ap. J. 328, 93.
• Henon, M. 1973, Astr. Ap. 24, 229.
• Jaffe, W. 1983, MNRAS 202, 995.
• Levison, H., Duncan, M.J., & Smith, B.F. 1990, Ap. J. 363, 66.
• Levison, H. & Richstone, D. 1985, Ap. J. 295, 349.
• Merritt, D. 1987, Ap. J. 319, 55.
• Merritt, D. & Aguilar, L.A. 1985, MNRAS 217, 787.
• Merritt, D. & Hernquist, L. 1991, Ap. J. 376, 439.
• Merritt, D. & Sellwood, J.A. 1994, Ap. J. 425, 551.
• Merritt, D. & Stiavelli, M. 1990, Ap. J. 358, 399.
• Newton, A.J. & Binney, J. 1984, MNRAS 210, 711.
• Palmer, P.L. & Papaloizou, J. 1987, MNRAS 224, 1043.
• Palmer, P.L., Papaloizou, J., & Allen, A.J. 1990, MNRAS 243, 282.
• Polyachenko, V.L. 1981, Soviet Astr. Lett. 7, 79.
• Schwarzschild, M. 1979, Ap. J. 232, 236.
• Toomre, A. 1964, Ap. J. 139, 1217.
• Toomre, A. 1966, in Notes on the 1966 Summer Study Program in Geophysical Fluid Dynamics at the Woods Hole Oceanographic Institution, p. 111.

## Homework

Due date: 3/20/97

These problems explore the picture of radial-orbit instability due to Palmer & Papaloizou (1987). The first two ask you to study an orbit in a Hernquist potential with a weak bar; the third seeks to show that a collection of such orbits will yield a bar-shaped mass distribution.

15. In Cartesian coordinates (X,Y), the potential of a Hernquist (1990) model with mass M and scale radius a is

```                                G M
(7)         Phi(X,Y) = - ------------------- .
sqrt(X^2 + Y^2) + a
```
Derive the equations of motion for a test particle moving in this potential, and implement them in the accel routine of the program leapint that you got in Lecture 6. Setting G = 1, M = 1 and a = 1/2, plot the orbit of a test particle starting at position (X,Y) = (1,0) with velocity (v_X,v_Y) = (0,v_0), where v_0 is a small value compared to the circular velocity at radius 1. Verify that the resulting orbit is a rosette, and adjust v_0 so that the test particle reaches apocenter about 20 to 50 times while precessing once around the origin. Show the resulting orbit, along with your value for v_0. (Hint: set the `number of points', now more properly called the `number of dimensions', n = 2; use two elements of the vector x for the coordinates (X,Y), and two elements of the vector v for the velocities (v_X,v_Y); set mstep = 32768, nout = 32, and dt = 1.0/256.0 to run for a sufficient time while outputting a reasonable amount of data).

16. Starting with the system in problem 15, add a weak bar-like potential aligned with the X axis,

```                                G M            G m  2
(8)         Phi(X,Y) = - ------------------- + --- Y  ,
sqrt(X^2 + Y^2) + a    a
```
where m << M parametrizes the strength of the bar. Re-run the initial conditions for the test particle orbit from problem 15 in this new potential, and adjust the parameter m so that the orbit can no longer precess completely around the origin but rather remains trapped within some angle of the X axis. Plot the resulting orbit.

17. In problem 16 you probably noticed that the orbit precesses fastest when it's aligned with the bar; consequently the orbit contributes less to the density along the bar than it does elsewhere. But the bar as a whole is built from an ensemble of orbits which oscillate through different angles with respect to the X axis, and the sum of all these orbits gives a density which is highest along the bar. Illustrate this by considering a simpler system: an ensemble of simple harmonic oscillators, with oscillation amplitudes x_max uniformly distributed between 0 and 1. Assuming the oscillators have random phases, what is the density distribution along the x axis?

Joshua E. Barnes (barnes@galileo.ifa.hawaii.edu)