Initial Conditions

September 2002

1  Mass Models

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⁻⁴; 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 > r_t >> a:

ρ(r) =        ρa  a4 rγ (a + r)(4−γ)  , r ≤ rt ρ* exp(−r/r*) (r/r*)−2  , r > rt

(1)

where ρ_a is the density at radius r = a. The parameters ρ_* and r_* can be chosen so that both ρ(r) and dρ/dr are continuous at r = r_t. In practice, we assumed that at r_t >> a the logarithmic slope is equal to its asymptotic value of −4, and set r_* = r_t / 2 and ρ_* = 4 e²ρ_a a⁴ r_t^−γ (a+r_t)^(γ−4).

Figure 1: Profiles of tapered gamma models with γ = 0.5, 1, 1.5, 2. All have M = 1 and r1/2 = 0.25. Solid lines are density profiles ρ(r) (left-hand scale). Dashed lines are cumulative mass profiles M(r) (right-hand scale).

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

2  Realization

We use Eddington's formula (Binney & Tremaine 1987, p. 237), in the form

f(E) = 8−1/2 π−2    d dE ⌠ ⌡ 0 E  dΦ (Φ− E)−1/2    dρ dΦ  ,

(2)

to calculate the distribution function f = f(E) for a given density profile. Here the integrand depends on the derivative of the density expressed as a function of the potential: ρ = ρ(Φ). There is no requirement that the density and potential actually be related by Poisson's equation. This provides several freedoms; first, it permits the construction of composite models in which several components co-exist in their mutual gravitational potential, and second, it allows the potential to be modified to allow for smoothing (a.k.a. `softening') in the N-body simulation. The first freedom is not used in the present study (it would be useful if one wanted to include dark halos). The second freedom is crucial; it allows N-body realizations of power-law density profiles to be set up in nearly perfect equilibrium with the gravitational field actually used in the N-body experiments.

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

Mε(r) = [1 + (3/2)(κ/γ) (ε/r)κ](−γ/κ) M(r)  .

(3)

This fitting formula was developed by considering the effects of Plummer smoothing on some power-law density profiles; here κ is a parameter, typicaly between 1.0 and 2.0, which adjusts the shape of the fit. We take ε = 0.005 in the present experiments; this is small enough that the slopes of the central cusps are resolved tolerably well. The resulting M_ε(r) is then used to compute the gravitational potential:

dΦ dr =  Mε(r) r2  ,

(4)

where we assume that Φ→ 0 as r → ∞.

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 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 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² f(0.5 v² + Φ(r_i)) using von Neumann rejection.

Last modified: September 16, 2002
http://www.ifa.hawaii.edu/~barnes/research/cusp_mergers/init_cond [.ps]