Initial Conditions
September 2002
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
falloff 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 > r_{t} >> a:
 (1) 

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 halfmass radius r_{1/2} = 0.25; the taper radius r_{t} = 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
form
 (2) 
The effect of smoothing is to take the pointmass of each body and
spread it out in a little cloud; the resulting density field generates
the Nbody potential via an accurate numerical solution of Poisson's
equation. Our Nbody 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
 (3) 
 (4) 
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 r_{t} = ∞ and ε = 0, it reproduces Hernquist's (1990) analytic solution for the distribution function of a γ = 1 model. And as will be shown, it permits Nbody models to be set up quite close to equilibrium.
Once f(E) has been calculated, actual generation of Nbody initial conditions is fairly straightforward. The body with index i is placed at a random point on a shell of radius r_{i}, where r_{i} is chosen by generating a random number X between 0 and total mass M and setting r_{i} = r_{M}(X), the inverse of the mass profile M(r). Each body is given a velocity with random direction and magnitude v_{i}, where v_{i} is chosen from the speed distribution v^{2} f(0.5 v^{2} + Φ(r_{i})) using von Neumann rejection.
Last modified: September 16, 2002
http://www.ifa.hawaii.edu/~barnes/research/cusp_mergers/init_cond
[.ps]