The mass models used in these collisions are based on the ``gamma
model'' (Dehnen 1993; Tremaine et al. 1994), which has a density
profile of the form ρ ∝ r−γ (a+r)(γ−4),
where a is a scale length and γ ∈ [0, 3) parametrizes the
logarithmic slope at r << a. These models include the Jaffe (1987)
model (γ = 2) and the Hernquist (1990) model (γ = 1).
At r >> a the density profile falls off as r−4; while this
fall-off is rapid enough to insure that the total mass is finite, it's
a bit inconvenient for numerical work. The profiles used for the
simulations are therefore tapered at radii r > rt >> a:
Fig. 1 shows density and mass profiles for the tapered gamma models considered here. All have been scaled to have total mass M = 1 and half-mass radius r1/2 = 0.25; the taper radius rt = 4. This scaling, together with the assumption that the gravitational constant G = 1, defines the units used for the numerical experiments.
We use Eddington's formula (Binney & Tremaine 1987, p. 237), in the
The effect of smoothing is to take the point-mass of each body and
spread it out in a little cloud; the resulting density field generates
the N-body potential via an accurate numerical solution of Poisson's
equation. Our N-body code
adopts Plummer smoothing (Aarseth 1963), in effect convolving the
density profile with a Plummer (1911) model which has length scale
ε. In terms of the cumulative mass profile M(r), Plummer
smoothing may be approximated by
Given a density profile (1), the corresponding cumulative mass profile M(r) can be evaluated analytically (indeed, the particular form of the tapering function was chosen to simplify this calculation). We tabulate ρ(r) and M(r) on a fine grid in logr, and use (3) to get the smoothed mass profile Mε(r) and (4) to tabulate the corresponding potential. Then ρ is expressed as a function of Φ and numerically differentiated. The integral in (2) is numerically evaluated for a grid of E values between Φ(r = 0) and 0, and the results are numerically differentiated to obtain f(E). Despite two numerical differentiations, this procedure is stable. With rt = ∞ and ε = 0, it reproduces Hernquist's (1990) analytic solution for the distribution function of a γ = 1 model. And as will be shown, it permits N-body models to be set up quite close to equilibrium.
Once f(E) has been calculated, actual generation of N-body initial conditions is fairly straightforward. The body with index i is placed at a random point on a shell of radius ri, where ri is chosen by generating a random number X between 0 and total mass M and setting ri = rM(X), the inverse of the mass profile M(r). Each body is given a velocity with random direction and magnitude vi, where vi is chosen from the speed distribution v2 f(0.5 v2 + Φ(ri)) using von Neumann rejection.
Last modified: September 16, 2002