14. Axisymmetric & Triaxial Models

Astronomy 626: Spring 1997

Axisymmetric two-integral models of many elliptical galaxies don't fit the observed kinematics, implying that these systems have a third, non-classical integral. Schwarzschild's method permits the construction of both axisymmetric and triaxial models which implicitly depend on a third integral. But minor orbit families and chaotic orbits both limit the range of axial ratios permitted for `cuspy' elliptical galaxies; systems with steep power-law cusps probably evolve toward more axisymmetric shapes over tens of orbital times.

14.1 Axisymmetric Models

Numerical calculations show that most orbits in `plausible' axisymmetric potentials have a third integral, since their shapes are not fully specified by the classical integrals, energy and the z-component of angular momentum. No general expressions for this non-classical integral are known, although for nearly spherical systems it is approximated by |J|, the magnitude of the total angular momentum, while in very flattened systems the energy invested in vertical motion may be used.

Despite the existence of third integrals in most axisymmetric potentials, it is reasonable to ask if models based on just two integrals can possibly describe real galaxies. In such models the distribution function has the form

```(1)         f = f(E,J_z) .
```
One immediate result is that the distribution function depends on the R and z components of the velocity only through the combination v_R^2+v_z^2; thus in all two-integral models the velocity dispersions in the R and z directions must be equal:
```            _______   _______
(2)         v_R v_R = v_z v_z ,
```
This equality does not hold in the Milky Way; for disk stars, the radial dispersion is about twice that in the vertical direction (MB81). Thus our galaxy can't be described by a two-integral model. For other galaxies, however, the situation is not so clear, and a two-integral model may suffice.

Distribution functions

From f(E,J_z) to rho(R,z): Much as in the spherically symmetric case described before, one may adopt a plausible guess for f(E,J_z), derive the corresponding density rho(R,Phi), and solve Poisson's equation for the gravitational potential. Perhaps the most interesting example of this approach is a series of scale-free models with r^-2 density profiles (Toomre 1982); however, these models are somewhat implausible in that the density vanishes along the minor axis.

From rho(R,z) to f(E,J_z): Conversely, one may try to find a distribution function which generates a given rho(R,z). This problem is severely underconstrained because a star contributes equally to the total density regardless of its sense of motion about the z-axis; formally, if f(E,J_z) yields the desired rho(R,z), then so does f(E,J_z)+f_o(E,J_z), where f_o(E,J_z) = -f_o(E,-J_z) is any odd function of J_z. The odd part of the distribution function can be found from the kinematics since it determines the net streaming motion in the phi direction (BT87, Chapter 4.5.2(a)).

An `unbelievably simple' and analytic distribution function exists for the mass distribution which generates the axisymmetric logarithmic potential (Evans 1993). This potential, introduced to describe the halos of galaxies (Binney 1981, BT87, Chapter 2.2.2), has the form

```                  1  2      2     2   2    2
(3)         Phi = - v_0 ln(R_c + R + z  / q ) ,
2
```
where v_0 is the velocity scale, R_c is the core scale radius, and q is the flattening of the potential (the mass distribution is even flatter). The corresponding distribution function has the form
```                          2       4 E            4 E            2 E
(4)         f(E,J_z) = A J_z exp(-----) + B exp(-----) + C exp(-----) ,
v_0^2          v_0^2          v_0^2
```
where A, B, and C are constants. Evans also divides this distribution function up into `luminous' and `dark' components to obtain models of luminous galaxies embedded in massive dark halos; his results illustrate a number of important points, including the non-gaussian line profiles which result when the luminous distribution function is anisotropic.

But even if kinematic data is available, this approach is not very practical for modeling observed galaxies. The reason is that the transformation from density (and streaming velocity) to distribution function is unstable; small errors in the input data can produce huge variations in the results (e.g. Dejonghe 1986, BT87). A few two-integral distribution functions are known for analytic density distributions, and recent developments have removed some mathematical obstacles to the construction of more models (Hunter & Qian 1993).

Jeans-equation models

Since we can't construct distribution functions for real galaxies, consider the simpler problem of modeling observed systems using the Jeans equations. If we assume that the underlying distribution function depends only on E and J_z we can simplify the Jeans equations, since the radial and vertical dispersions must be everywhere equal; thus

```            d     _______   nu  _______   ___________         dPhi
(5)         -- nu v_R v_R + -- (v_R v_R - v_phi v_phi) = - nu ---- ,
dR              R                                  dR

d     _______                                     dPhi
(6)         -- nu v_R v_R                              = - nu ---- .
dz                                                 dz
```
At each R one can calculate the mean squared velocity in the R direction by integrating Eq. 6 inward from z = infinity; the mean squared velocity in the phi direction then follows from Eq. 5.

The Jeans equations don't tell how to divide up the azimuthal velocities into streaming and random components. One popular choice is

```            _____ 2      ___________   _______
(7)         v_phi   = k (v_phi v_phi - v_R v_R) ,
```
where k is a free parameter (Satoh 1980). The dispersion in the phi direction is then
```                 2      ___________   _____ 2
(8)         sigma_phi = v_phi v_phi - v_phi   .
```
Note that if k = 1 the velocity dispersion is isotropic and the excess azimuthal motion is entirely due to rotation, while for k < 1 the azimuthal dispersion exceeds the radial dispersion.

Application to elliptical galaxies

Jeans equations models have been constructed of a number of elliptical galaxies (Binney, Davies, & Illingworth 1990, van der Marel, Binney, & Davies 1990, van der Marel 1991). The procedure is to:

1. Observe the surface brightness Sigma(x',y');
2. Deproject to get the stellar density nu(R,z), assuming an inclination angle;
3. Compute the potential Phi(R,z), assuming a constant mass-to-light ratio;
4. Solve the Jeans equations for the mean squared velocities;
5. Divide the azimuthal motion into streaming and random parts;
6. Project the velocities back on to the plane of the sky to get the line-of-sight velocity and dispersion v_los(x',y') and sigma_los(x',y');
7. Compare the predicted and observed kinematics.
The inclination angle, mass-to-light ratio, and rotation parameter k are unknowns to be determined by trial and error.

Some conclusions following from this exercise are that:

• Isotropic oblate rotators (k = 1) generally don't fit, even though some galaxies lie close to the expected relation between v_0/sigma_0 and epsilon;
• Some galaxies (e.g. NGC 1052) are well-fit by two-integral Jeans equation models;
• The models predict major-axis velocity dispersions in excess of those observed in most galaxies;
• Consequently, most of the galaxies must have a third integral, or may even be triaxial.

14.2 Triaxial Models

The general problem of modeling triaxial galaxies is illustrated with a couple of special cases. In separable models the allowed orbits are fairly simple, and the job of populating orbits so as to produce the density distribution is well-understood. In scale-free models the allowed orbits are more complex, and it's not clear if such models can be in true equilibrium.

Schwarzschild's Method

Schwarzschild (1979) invented a powerful method for constructing equilibrium models of galaxies without explicit knowledge of the integrals of motion. To use this method,

1. Specify the mass model rho(x) and find the corresponding potential;
2. Construct a grid of K cells in position space;
3. Chose initial conditions for a set of N orbits, and for each one,
• integrate the equations of motion for many orbital periods, and
• keep track of how much time the orbit spends in each cell, which is a measure of how much mass the orbit contributes to that cell;
4. Determine non-negative weights for each orbit such that the summed mass in each cell is equal to the mass implied by the original rho(x).

The last step is the most difficult. Formally, let P_i(c) be the mass contributed to cell c by orbit i; the task is then to find N non-negative quantities Q_i such that

```                   -- N
\
(9)         M(c) =  ) Q_i P_i(c)
/
-- i = 1
```
simultaneously for all cells, where M(c) is the integral of rho(x) over cell c. This may be accomplished by taking N > K, so as to obtain a reasonably rich set of `basis functions', and using any one of a number of numerical techniques, including linear programming (Schwarzschild 1979), non-negative least squares (Pfenniger 1984), Lucy's method (Newton & Binney 1984), or maximum entropy (Richstone & Tremaine 1988).

In general Eq. 9 has many solutions, reflecting the fact that many different distribution functions are consistent with a given mass model. Some methods allow one to specify additional constraints so as to select solutions with special properties (maximum rotation, radial anisotropy, etc.).

Separable Potentials

In three dimensions, a separable potential permits four distinct orbit families:

• box orbits,
• short-axis tube orbits,
• inner long-axis tube orbits, and
• outer long-axis tube orbits.

The time-averaged angular momentum of a star on a box orbit vanishes; such orbits therefore do not contribute to the net rotation of the system. Short-axis and long-axis tube orbits, on the other hand, preserve a definite sense of rotation about their respective axes; consequently, their time-averaged angular momenta do not vanish. The total angular momentum vector of a non-rotating triaxial galaxy may thus lie anywhere in the plane containing the short and long axes (Levison 1987).

Using Schwarzschild's method, it is possible to numerically determine orbit populations corresponding to separable triaxial models (Statler 1987). A somewhat more restricted set of models can be constructed exactly; these models make use of all available box orbits, but only those tube orbits with zero radial thickness (Hunter & de Zeeuw 1992). Apart from the choice of streaming motion, thin tube models are unique. One use of such models is to illustrate the effects of streaming motion by giving all tube orbits the same sense of rotation; the predicted velocity fields display a wide range of possibilities (Arnold, de Zeeuw, & Hunter 1994). Non-zero streaming on the projected minor axis is a generic feature of such models; a number of real galaxies exhibit such motions and thus must be triaxial (Franx, Illingworth, & Heckman 1989b).

Scale-Free Potentials

In scale-free models, box orbits tend to be replaced by minor orbital families known as boxlets (Gerhard & Binney 1985, Miralda-Escude & Schwarzschild 1989). Each boxlet family is associated with a closed and stable orbit arising from a resonance between the motions in the x and y directions.

The appearance of boxlets instead of boxes poses a problem for model building because boxlets are `centrophobic' (meaning that they avoid the center) and do not provide the elongated density distributions of the box orbits they replace. As a result, the very existence of scale-free triaxial systems is open to doubt (de Zeeuw 1995).

Moreover, some scale-free potentials have irregular orbits; these have no integrals of motion apart from the energy E. In principle, such an orbit can wander everywhere on the phase-space hypersurface of constant E, but in actuality such orbits show a complicated and often fractal-like structure.

The scale-free elliptic disk is a relatively simple two-dimensional analog of a scale-free triaxial system. Because the model is scale-free, the radial dimension can be folded out when applying Schwarzschild's method; thus the calculations are fast (Kuijken 1993). The result is that self-consistent models can be built using the available boxlets, loops, and irregular orbits, but the minimum possible axial ratio b/a increases as the numerical resolution of the calculation is improved.

Similar results hold for scale-free models in three dimensions. Models have been constructed for triaxial logarithmic potentials with a range of axial ratios b/a and c/a (Schwarzschild 1993). Tubes and regular boxlets provide sufficient variety to produce models with c/a > 0.5, but not flatter. However, over intervals of 50 dynamical times, irregular orbits behave like `fuzzy regular orbits', and by including them it becomes possible to build near-equilibrium models as flat as c/a = 0.3. But these models are not true equilibria; over long times they will become rounder and less triaxial.

14.3 Rotation, Chaos, & Secular Evolution

Figure rotation adds a new level of complexity to the orbit structure of triaxial systems. A few models have been constructed using Schwarzschild's method, but little is known about the existence and stability of such systems. N-body experiments indicate that at least some such systems are viable models of elliptical galaxies. Rotation tends to steer orbits away from the center and so may lessen the effects of central density cusps.

Realistic potentials are likely to have some irregular or chaotic orbits, and there is no reason to think that such orbits are systematically avoided by processes of galaxy formation. Over 10^2 or more dynamical times, such orbits tend to produce nearly round density distributions (Merritt & Valluri 1996).

Consequently, it is likely that secular evolution over timescales of 10^2 dynamical times may be changing the structures of elliptical galaxies (Binney 1982, Gerhard 1986). The outer regions are not likely to be affected since dynamical times are long at large radii, but significant changes may occur in the central cusps where dynamical times are only 10^6 years (de Zeeuw 1995).

References

• Arnold, R.A., de Zeeuw, P.T., & Hunter, C. 1994, M.N.R.A.S., submitted.
• Binney, J.J. 1981, M.N.R.A.S. 196, 455.
• Binney, J.J. 1982, M.N.R.A.S. 201, 15.
• Binney. J.J., Davies R.L., & Illingworth, G.D. 1990, Ap.J. 361, 78.
• Dejonghe, H. 1986, Phys. Rep. 133, 218.
• de Zeeuw, P.T. 1995, in The Formation and Evolution of Galaxies, ed. C. Munoz-Tunon & F. Sanchez, p. 231.
• Evans, N.W. 1993, M.N.R.A.S. 260, 191.
• Franx, M., Illingworth, G.D. & Heckman, T.M. 1989b, Ap. J. 344, 613.
• Gerhard, O.E. 1986, M.N.R.A.S. 219, 373.
• Gerhard, O.E. & Binney, J.J. 1985, M.N.R.A.S. 216, 467.
• Hunter, C. & de Zeeuw, P.T. 1992, Ap.J 398, 79.
• Hunter, C. & Qian, E. 1993, M.N.R.A.S. 262, 401.
• Kuijken, K. 1993, Ap.J. 409, 68.
• Levison, H.F. 1987, Ap.J. 320, L93.
• Merritt, D. & Valluri, M. 1996, Ap.J. 471, 82.
• Miralda-Escude, J. & Schwarzschild, M. 1989, Ap.J. 339, 752.
• Newton, A. & Binney, J.J. 1984, M.N.R.A.S. 210, 711.
• Pfenniger, D. 1984, Astr.Ap. 141, 171.
• Richstone, D.O. & Tremaine, S.D. 1988, Ap.J. 327, 82.
• Satoh, C. 1980, P.A.S.J. 32, 41.
• Schwarzschild, M. 1979, Ap.J. 232, 236.
• Schwarzschild, M. 1993, Ap.J. 409, 563.
• Statler, T.S. 1987, Ap.J. 321, 113.
• Toomre, A. 1982, Ap.J. 259, 535.
• van der Marel, R. 1991, M.N.R.A.S. 253, 710.
• van der Marel, R., Binney, J.J., & Davies, R.L. 1990, M.N.R.A.S. 245, 582.

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