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).

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

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 vHere the first three terms describe the evolution of the perturbed distribution function in the unperturbed field

/ (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

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

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.

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.

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

sigma_r^2 1 (6) --------- = ----- , sigma_t^2 1 + mso

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).

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 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 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).

- 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.

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) + aDerive the equations of motion for a test particle moving in this potential, and implement them in the

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 awhere

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)

Last modified: March 14, 1997