### 2. METHODS

Simple estimates imply that relaxation times of galaxies are many orders of magnitude longer than a Hubble time (e.g. Binney & Tremaine1987). Consequently, the dynamics of stars in galaxies is well-described by the collisionless Boltzmann equation

ðf / ðt + v · ðf / ðv - · ðf / ðv = 0, 1.

where f = f(x, v, t) is the distribution function normalized so that f dx dv is the number of stars in the phase-space volume dx dv centered on the point (x, v). The potential, , includes the self-consistent field generated by the stars as well as the gravitational influence of dark matter and gas. While the nature of the dark matter remains enigmatic, observational constraints suggest that it too obeys Equation 1.

N-Body Methods

In the purely stellar-dynamical limit, the interaction of galaxies is a well-posed problem, in the sense that the evolution is determined entirely by the Boltzmann equation. Unfortunately, we don't know how to obtain relevant analytic solutions to this equation, except in the limit of weak interactions (e.g. Weinberg 1986, 1989) or for rather special systems (e.g. Sridhar & Nityananda 1990), so numerical methods must be used. Rather than solve the Boltzmann equation with a finite-difference technique in six-dimensions, a Monte-Carlo approach is employed. The initial f(x,v) is used as a probability distribution function to pick phase-space coordinates xi and vi for particles i = 1, . . . , N; these particles are then integrated along the characteristic curves of Equation 1:

dxi / dt = vi , 2.

dvi / dt = -i . 3.

Despite the enormous simplification entailed by Equations 2 and 3, the evolution of interacting galaxies is quite costly to compute, owing mainly to the difficulty of calculating the self-consistent potential. In practice, the mass in real galaxies is much more finely divided than in any computer model, so the model's self-consistent potential is much noisier than that of actual galaxies. This gives rise to an unwanted diffusion in phase-space which can alter the structure of simulated galaxies on uncomfortably short time-scales.

Broadly speaking, the various known algorithms for computing the potential of a system of N particles can be categorized as either ``action-at-a-distance'' or ``field'' methods. The former explicitly treat interactions between individual particles while the latter do so only indirectly through the contributions of particles to the global gravitational field. The simplest action-at-a-distance technique is direct summation. A common expression for the gravitational potential at the location of particle i is

(ri) = - G j i (mj / [ | ri - rj |2 + 2]1/2) , 4.

where is the ``softening parameter.'' Direct summation is flexible but has an asymptotic cpu cost per step scaling as ~ O (N2), limiting its practical use to small-N systems. Nevertheless, some early investigations of the dynamics of mergers were performed using this algorithm (e.g. White 1978, 1979, 1980; Gerhard 1981; Quinn 1982) and much larger simulations with N > 10,000 may ultimately be feasible on special-purpose machines (e.g. Sugimoto et al. 1990).

A variety of methods have been proposed to reduce the cost of computing the self-consistent potential. For some problems, it is expedient to neglect self-gravity altogether and evolve the system in a potential which is known a priori, such as in test-particle methods (e.g. Schwarzschild 1979; Quinn 1984), the restricted three-body method (e.g. Pfleiderer & Siedentopf 1961; Toomre & Toomre 1972), or semi-restricted N-body codes (e.g. Lin & Tremaine 1983; Quinn & Goodman 1986; Hernquist & Weinberg 1989). Generally speaking, however, the dynamics of interacting galaxies will not be represented faithfully without including self-gravity (Barnes 1988).

The majority of the known self-consistent methods have historically not been useful for studying systems with small filling factors in three-dimensions such as interacting and merging galaxies. (We refer the reader to the excellent monographs by Sellwood [1987] and Hockney & Eastwood [1988] for a discussion of the limitations of grid methods.) An exception is the technique known as the ``hierarchical tree'' method (Appel 1985; Jernigan 1985; Barnes & Hut 1986) which recently has had a significant impact on the study of colliding galaxies. All variants of this technique are effectively gridless and retain the essential ingredients of the action-at-a-distance formulation. Particles are processed similarly to direct summation, but the potential from distant groups of particles is approximated using low-order multipole expansions. If the required manipulations are performed using tree-structured data the cpu cost per step scales as ~ O (N log N). In fact, it is possible to reduce the cost even further to ~ O (N) by using tree-structured data in the context of a field representation of the potential (Greengard & Rokhlin 1987; Greengard 1988). A distinct approach, used in ``expansion codes,'' relies on the use of basis function expansions to compute the potential from the known density field sampled by the particles (Aarseth 1967; Clutton-Brock 1972a, b; van Albada & van Gorkom 1977; Villumsen 1982, 1983; White 1983a; McGlynn 1984; Bontekoe & van Albada 1987).

As emphasized by Sellwood (1987), the choice of method is largely dictated by considerations of efficiency. Contrary to claims which persist in the literature, action-at-a-distance and field methods suffer from comparable degrees of relaxation for the same spatial resolution and particle number (Hernquist & Barnes 1990; Hernquist & Ostriker 1992). Therefore, there is no objection in principle to using direct methods to evolve collisionless systems, provided that a sufficiently large N can be used to suppress relaxation effects to the extent required.

While present-day galaxies typically contain less than 10% of their luminous mass in this form, gas is critical to our understanding of galactic phenomena such as starbursts and active galactic nuclei, because it can dissipate and form stars. Galactic gas obeys the ordinary conservation laws for a compressible fluid (e.g. Landau & Lifshitz 1959). However, given our present level of understanding of the interstellar medium of even our own galaxy, the dynamics of galactic gas is ill-posed. Because of cooling and feedback from supernovae, interstellar gas is highly inhomogeneous and comprises several distinct phases (McKee & Ostriker 1977). We lack the computational hardware and theoretical knowledge to model the different phases simultaneously, and some compromises must necessarily be made.

One strategy, used to advantage by some workers, is to focus on global issues where the detailed small-scale structure of the gas is probably not crucial. In the limit where the gas can be represented as a continuous medium, its equations of motion can be solved either by finite-difference algorithms or by particle-based techniques. A promising example of the latter is the method known as smoothed-particle hydrodynamics (SPH; e.g. Lucy 1977; Gingold & Monaghan 1977; Monaghan 1992). Since SPH is Lagrangian, it does not suffer from grid-based limitations on spatial resolution or global geometry and is especially well-suited for studying interacting systems. In addition, SPH can be generalized so that it is adaptive in both space and time (Hernquist & Katz 1989); these refinements offer significant advantages when dealing with large density contrasts (Hernquist 1989a; Barnes & Hernquist 1991).

A phenomenological treatment is sometimes useful in the opposite extreme when the gas is modeled as a set of discrete clouds (e.g. Negroponte & White 1983). In the ``sticky-particle'' method, particles representing the gas evolve similarly to collisionless particles but undergo inelastic collisions amongst themselves, thereby dissipating kinetic energy. Computationally, this approach is simpler than those which solve the equations of motion for continuum gas, but requires a number of ad hoc parameters whose relationship to actual physical processes is problematic. Given our ignorance of the true state of interstellar gas it is difficult to argue for or against one of the two extremes embodied by, e.g. the SPH or sticky-particle formalisms.