Initial Data

September 2002

The galaxy models used in these experiments have a central bulge, an exponential disk, and a spherical dark halo. Disks are actually represented in two different ways; as flat, rotating structures corresponding to the usual notion of a disk, or as spherical, isotropic distributions with the same cumulative mass profiles as an exponential disk. When the latter disk representation is used, distribution functions for all components may be calculated exactly using Abel integrals.

1  Mass models

For ease of comparison with earlier work, we initially adopted galaxy models similar to those used by Velazquez & White (1996), with spherical bulges following Hernquist (1990) profiles, exponential & isothermal disks, and non-singular isothermal halos (Hernquist 1993). But a problem came to light in trying to realize the halo using an isotropic distribution function; the r⁻¹ profile of the inner bulge creates a potential with a nonzero gradient as r → 0, and in such a potential well a halo with a constant-density core can only be realized with an anisotropic distribution function. We therefore replaced the isothermal halo profile with an NFW profile; the latter has the same profile as the bulge as r → 0, so it can be realized with an isotropic distribution function.

For numerical convenience, the bulge profile was tapered off smoothly at large radii:

ρb(r) =          ρba  ab4 r (ab+r)3  , r ≤ bb ρbb exp[2 − 2r/bb] (r/bb)−2  , r > bb

(1)

where ρ_ba is the bulge density at the scale radius r = a_b, and ρ_bb = ρ_ba a_b⁴ b_b⁻¹(a_b+b_b)⁻³ is the density at the taper radius r = b_b.

The NFW halo has a logarithmically diverging mass, so it must be tapered at large radii; we used the form adopted by Springel & White (1998):

ρh(r) =          ρha  ah3 r (ah+r)2  , r ≤ bh ρhb exp[Q − r/ah] (r/bh)γ  , r > bh

(2)

where ρ_ha is the halo density at the scale radius r = a_h, ρ_hb = ρ_ha a_h³ b_h⁻¹ (a_h+b_h)⁻² is the density at the taper radius r = b_h, the constant Q = b_h/a_h, and γ = b_h / (b_h + a_h) − 0.5.

The disk's density profile depends on cylindrical radius R and height z:

ρd(R,z) =  Md 4 πz0 Rd2 e− R/Rd sech2(z/z0)  .

(3)

Parameters for the galaxy model are listed below. For each component, M is the mass, a is the scale radius, and b is the taper radius.

These quantities are reported in arbitrary units with G = 1. Any such system of units is internally self-consistent, and the results may be rescaled as needed. To scale this model to the Milky Way, for example, we may equate the disk's scale length R_d and mass M_d to the scale length and mass of the MW's disk, R_d = 3.5  kpc and M_d = 6 ×10¹⁰ M_\odot.

Figure 1: Left: density profiles for each component of the galaxy model. Heavy curve is NFW halo; dotted curve is truncated isothermal halo. Bulge dominates density at small radii. Right: mass profiles for each component of the galaxy model.

Fig. 1 shows density and mass profiles for each component of the galaxy model. Here the disk is represented by a spherical configuration which has the same mass profile M(r). At all radii the spherically-averaged disk density and mass profiles fall below the corresponding halo profiles. This is not the case in the V&W model - there the constant density of the halo core allows the disk to dominate (slightly) at intermediate radii.

2  Spherical realizations

Initial data are generated from isotropic distribution functions f_b(E), f_s(E), and f_h(E), which represent the bulge, `spherical disk', and halo, respectively. We calculated these distribution functions using an Abel equation which relates component c's density profile ρ_c(r) and distribution function f_c(E) to the potential generated by all components, Φ(r). To allow for the finite resolution of the N-body force calculation, we smoothed the total density profile ρ(r) with a Plummer kernel before evaluating the potential function Φ(r). This smoothing procedure employs a semi-empirical fitting function with a free parameter κ; the tests below were run to establish the optimal value of this parameter.

The N-body realizations used in these tests had N_b = 20480 bulge particles, N_s = 65536 `spherical disk' particles, and N_h = 229376 halo particles. Plummer smoothing with a length scale of ε = 0.05 length units was adopted. Forces were calculated with a modified tree-code. Bodies were advanced with a time-centered leap-frog; for these initial tests a time-step ∆t = 1/16 was used.

κ = 1.15	κ = 1.25	κ = 1.30	κ = 1.35

Figure 2: Evolution of binding energies in realizations generated using different values of κ. Each image is linked to a short animation showing times t = 0 to t = 8. Horizontal axis shows average binding energy (Ei(t) + Ei(0)) / 2 of each body i, ranging from −5.625 (left) to 0.0 (right). Vertical axis shows relative change in binding energy 2 (Ei(t) − Ei(0)) / (Ei(t)+ Ei(0)), ranging from −10% (bottom) to +10% (top). Yellow, blue, and red represent bulge, disk, and halo components.

Fig. 2 and the associated animations shows how bodies evolve in binding energy E_i during a fairly short run of 8 time units. As the animations make clear, bodies are scattered in binding energy by stochastic events and by large-scale potential fluctuations; the former give rise to a diffusive process, while the latter create coherent patterns of motion. These coherent motions are very nearly absent in the κ = 1.35 model, indicating that this model is closest to equilibrium.

3  Disk realizations

Initial data for the disk is generated using a routine which sets up a self-gravitating disk in approximate equilibrium with the external potential due to the bulge and halo. Initial data for the bulge and halo are generated using the same Abel integral technique described above; to calculate the required distribution functions, we assume that the total potential is spherical. Of course, the actual potential is aspherical because of the disk, but this does not seem a significant source of error in the present experiments.

Figure 3: Evolution of surface density in disk realizations. Each image is linked to an animation showing times t = 0 to t = 64; time t is shown at upper right. Grey-scale shows departure from underlying exponential profile; black is half density, white is twice density. Left: disk models with z0 = 0.2. Right: disk models with z0 = 0.1.

Figure 4: Evolution of binding energies in disk realizations. Each image is linked to an animation showing times t = 0 to t = 64; time t is shown at upper right. Left: disk models with z0 = 0.2. Right: disk models with z0 = 0.1.

Last modified: September 24, 2002
http://www.ifa.hawaii.edu/~barnes/research/halo_noise/inital_data.html [.ps]