df df df
(1) -- + v . -- - grad Phi . -- = 0 .
dt dr dv
The gravitational field Phi(r,t) is given by Poisson's equation.
/
|
(2) div grad Phi = 4 pi G | dv f(r,v,t) .
|
/
Eqs. 1 & 2 may be viewed as a pair of coupled PDEs which together
completely describe the evolution of a galaxy.
The Collisionless Boltzmann and Poisson equations can be effectively solved with a Monte-Carlo technique which represents the distribution function f(r,v,t) as a set of N bodies.
Impractically large grids are required to solve the time-dependent 3-D Collisionless Boltzmann Equation by finite-difference methods. N-body simulation is basically a Monte-Carlo method of solving this equation, with the number of bodies, N, governing the accuracy of the method.
The basic idea behind Monte-Carlo methods is shown by the following procedure for approximating pi. Draw a square of area A, and inscribe within it a circle of area A_c; by simple geometry, A_c = pi A / 4. Now scatter n points independently and randomly within the square, and count the number n_c which fall within the circle. Since the expected number of points within any area is proportional to that area, the quantity 4 n_c / n approximates pi, with a fractional uncertainty of order 1/sqrt(n). As a method of calculating pi, this procedure is inefficient; for example, a trial with n = 262144 points yielded the estimate pi = 3.138 +/- 0.007. But the error in a Monte-Carlo calculation does not depend on the number of dimensions, but only on the number of points. In evaluating multidimensional integrals, Monte-Carlo methods can outperform other numerical techniques.
To represent the mass distribution function f(r,v,t_0) at some instant t_0 in a form suitable for Monte-Carlo calculations, one uses a set of N bodies, each possessing a mass m_i, position r_i, and velocity v_i, where i = 1...N. In effect, the continuous distribution function is replaced with a set of delta-functions:
-- N
\
(3) f(r,v) -> ) m_i delta(r - r_i) delta(v - v_i) .
/
-- i = 1
For this substitution to work, the expected mass of the bodies within
any phase-space volume V must be equal to the integral
of the distribution function over that volume; thus,
/ / -- \
| 3 3 / \ \
(4) | d r d v f(r,v) = ( ) m_i ) ,
| \ / /
/ V \ -- (r_i,v_i) in V /
where the angle brackets indicate an average over statistically
equivalent realizations, and the sum includes all bodies with
phase-space coordinates within the volume V.
The simplest way to initialize bodies in accord with Eq. 4 is to pick phase-space coordinates by treating f(r,v) as a probability distribution; that is, select (r_i,v_i) with probability proportional to f(r_i,v_i), and assign all bodies the same mass m_i = M/N, where M is the total mass. Since bodies are selected independently, the actual number within any given volume V will have a Poissonian distribution about the mean. This scatter is characteristis of Monte-Carlo method and limits the accuracy of the calculation.
This probabilistic point representation of the distribution function may be used to calculate integrals of f(r,v) over phase-space. Suppose that we wish to estimate the value of some observable q, defined as the integral
/
| 3 3
(5) q = | d r d v f(r,v) Q(r,v) .
|
/
Using Eq. 3, this becomes
-- N
\
(6) q = ) m_i Q(r_i,v_i) .
/
-- i = 1
The fractional uncertainty in any individual estimate of q is
of order 1/sqrt(N) if bodies are selected independently, just
as in the the estimate of pi above.
N-body representations are useful for other things besides Monte-Carlo integrations; in particular, they are easily projected into the future. This is accomplished by moving bodies along the phase-flow defined by
. .
(7) (r,v) = (v, - grad Phi)
This projection preserves the relation Eq. 4; thus starting with a
valid realization of f(r,v,t_0), the result is a realization
of f(r,v,t) at time t > t_0.
Now what to use for the potential, Phi(r,t)? In a self-consistent calculation the potential is an unknown, to be estimated from the N-body representation of the mass distribution by using the bodies as the source term for Poisson's equation; thus,
-- N
\
(8) div grad Phi = 4 pi G ) m_i delta(r - r_i) .
/
-- i = 1
This yields the standard Newtonian equations of motion for N
point masses. Formally there is no conceptual problem using these
equations in a self-consistent N-body simulation. But in practice the
singular potential wells associated with point masses create awkward
numerical problems. These problems can be very naturally avoided by
smoothing the N-body representation of the density field, for example
via the substitution
3 epsilon^2
(9) delta(r - r_i) -> ---- ----------------------------- ,
4 pi (|r - r_i|^2 + epsilon^2)^5/2
where epsilon is a parameter with dimensions of length. The
resulting equations of motion are
-- N
d(r_i) d(v_i) \ m_j (r_j - r_i)
(10) ------ = v_i , ------ = G ) ------------------------------- .
dt dt / (|r_j - r_i|^2 + epsilon^2)^3/2
-- j#i
Apart from the smoothing parameter epsilon these equations
are very similar to the N-body equations of motion in Lecture 4. However, it's important to remember
that bodies are tracers of the phase-space distribution, and not to be
identified with individual stars.
The smoothing procedure incorporated into these equations is commonly called `Plummer softening', since it effectively replaces each point mass with a little Plummer (1911) model; or equivalently, it replaces the potential of each point with the potential of a Plummer model. The latter interpretation leads to the phrase `softened potential', and to worries that one has invalidated the simulations. Such worries are laid to rest if the smoothing inherent in these equations is recognized as distinct from the gravitational force calculation.
Besides removing singularities in the equations of motion, smoothing suppresses small-scale fluctuations due to the discrete nature of N-body representations. Discreteness fluctuations are the bane of collisionless N-body simulation, so some smoothing is a good thing. But smoothing always comes at a price in spatial resolution, and no useful amount of smoothing will completely eliminate the effects of discreteness (Hernquist & Barnes 1990)!
The Monte-Carlo interpretation of N-body simulation places a premium on large values of N to reduce uncertainties. Thus N-body experimenters seek to run more bodies with the same fervor that observational astronomers seek to gather more photons. Force calculation is the most computation-intensive part of N-body simulation. Starting with Holmberg's (1941) optical computations, much ingenuity has gone into the rapid evaluation of forces in N-body systems.
No single force calculation method is optimal for all applications. Hierarchical methods will be discussed here in the greatest detail since they are useful for galactic encounter simulations, which typically involve irregular mass distributions, require high spatial resolution, and demand many bodies. But no method will ever make N-body calculation `cheap'; faster computers and better algorithms simply shift the focus of attention to larger problems.
Straightforward evaluation of the sum in Eq. 8 is robust, accurate, and completely general. As everyone knows, the cost of evaluating the force on all N bodies is O(N^2). To some extent, direct summation methods can beat the high cost of force calculation by (1) efficiently using individual time-steps (eg. Ahmad & Cohen 1973, Aarseth 1985), and (2) implementing the force calculation in hardware (eg. Sugimoto et al. 1990). Such improvements have kept direct summation remarkably competitive even as N has increased by several factors of 10.
A number of techniques are available for integrating sets of ordinary differential equations such as Eq. 8. One popular method for N-body simulations of collisionless systems is the leap-frog method, discussed in Lecture 2. The main drawback to the leap-frog method is that all bodies are advanced with the same time-step, which must be small enough to follow motions throughout the system. Thus when only a small percentage of bodies require very small time-steps, considerable computation is wasted. The speed-up of an ideal individual time-step algorithm, as compared to the leap-frog, is S = max[1/Delta t_i] / < 1/Delta t_i>, where the the maximum and average are over all bodies. In a spherical system composed of equal-mass bodies, each with a time-step proportional to the circular orbit period at its present radius, the speed-up depends on the density profile. For a Plummer (1911) profile, rho(r) proportional to (r^2+a^2)^-5/2, the speed-up is only S = 2.0867 since most of the mass lies within the `core'. The speed-up is formally infinite for a mass model with a density profile diverging as r -> 0, but simulated density profiles are never singular after smoothing by Eq. 9.
There are two kinds of uncertainties in N-body simulations of collisionless systems. First, the dynamical equations (Eq. 10) are integrated numerically with less-than-perfect accuracy. Second, even exact solutions of these equations do not correspond to exact solutions of Collisionless Boltzmann and Poisson equations. Numerical errors are fairly easy to measure and control, but smoothing and discreteness effects have subtle implications for the interpretation of N-body simulations.
Tests of N-body codes are hampered by a lack of exact solutions with which to compare the output of numerical integrations. To verify the correctness of the simulations, one should ideally show that all uncertainties can be made as small as desired, and that the results converge to a unique limit as the calculation is refined.
Errors in numerical solutions include approximations introduced by a hierarchical force calculation algorithm, truncation caused by using a finite time-step, and roundoff due to finite computer word-length. All can be treated as small perturbations introduced at every time-step; their cumulative effects can be gauged by monitoring the conservation of energy and momentum, or studied in more detail by running the same set of initial conditions with different time-steps and opening angles. Convergence testing shows that it is generally possible to constrain the uncertainty associated with numerical errors at a `reasonable' computational cost.
Due date: 2/16/04
Using 1000000 random numbers between 0 and 1, calculate pi using the Monte Carlo procedure described in Lecture 7. Estimate the error of the calculation.
Following the approach in Eqs. 5,6 of Lecture 7, write down the functions Q(r,v)and q for (a) the center-of-mass position, (b) the total momentum, (c) the total energy, and (d) the total angular momentum with respect to the origin.
Verify the conservation laws for energy, momentum, and angular momentum in the time evolution of the gravitational sphere using the Aarseth code in your N-body simulation. You should choose epsilon^2=0.25, eta=0.02 and G=1. As an initial condition, generate a sphere with 0.5 radius and with randomly distributed point masses of m=0.5 for each star. The center of mass of the system at t=0 should be located at x=0.5, y=0.5, z=0.5 with a non-vanishing velocity vx=vy=vz=1.5 in the direction of the origin. You should run N=2000 particles for tcrit=0.1 which allows substantial changes in the macroscopic evolution of your system due to the gravitational bouncing. Plot the three conserved quantities and the center-of-mass as a function of simulation time. Compare the theoretical predictions with the simulation results in your plots.