View source for N-body simulations (gravitational)

==Introduction==

The underlying dynamics relevant in the astrophysical context for of a
system of N particles interacting gravitationally is typically
Newton's law plus, in case, an external potential field (see however
below for a discussion of N-body simulations in general
relativity). The force <math>\vec{F}_i</math> acting on particle
<math>i</math> of mass <math>m_i</math> is:

:<math>\label{eq:newton}
\vec{F}_i = - \sum_{j \ne i}  G \frac{m_i m_j (\vec{r}_i-\vec{r}_j)}{|\vec{r_i}-\vec{r_j}|^3 } - \vec{\nabla} \cdot \phi_{ext}(\vec{r}_i),  
</math>


where <math>G=6.67300 \cdot 10^{-11}</math> <math>m^{3}</math>
<math>kg^{-1}</math> <math>s^{-2}</math> is the gravitational
constant, and <math>\phi_{ext}</math> is the external potential. The problem
is thus a set of non-linear second order ordinary differential
equations relating the acceleration <math>\partial^2 \vec{r_i} /
\partial t^2 = \vec{F}_i /m_i</math> with the position of all the
particles in the system.

Once a set of initial condition is specified (for example the initial
positions <math>\vec{r}_i</math> and velocities <math>\vec{v}_i \equiv
\partial \vec{r}_i / \partial t </math> of all particles) it exists a
unique solution, analytical only for up to two bodies, while larger N
require numerical integration (e.g. see Press et al. 2007). However
special care must be employed to ensure both accuracy and
efficiency. In fact, the gravitational force
(eq. \eqref{eq:newton} presents a singularity when the distance of
two particles approaches 0, which can lead to arbitrarily large
relative velocities. In addition, given the non-linear nature of the
equations, the singularities are movable, that is they depend on the
specific choice of initial conditions. In contrast, all singularities
in linear ordinary differential equations are independent of initial
conditions and thus easier to treat. Therefore constant timestep
methods are unable to guarantee a given accuracy in the
case of gravitational dynamics and lead to unphysical accelerations
during close encounters, which in turn may create unbound stars. A
shared, adaptive timestep scheme can correctly follow a close
encounter, but the price is paid in terms of efficiency as all the
other particles of the system are evolved on the timescale of the
encounter, which may be several orders of magnitude smaller than the
global timescale, resulting essentially in a freezing of the
system.

The singularity may be avoided by introducing a softening length in
Eq. \eqref{eq:newton} (e.g. see Aarseth 1963), that is by modifying the gravitational
interaction at small scales. For example:

:<math>  
  \vec{F}_i = - \sum_{j \neq i} \frac{G m_i m_j
    (\vec{r}_i - \vec{r}_j)}{(|\vec{r}_i - \vec{r}_j|^2 + \epsilon^2)^{3/2}} ,
</math>

where <math>\epsilon > 0</math> is the softening length, that is a
typical distance below which the gravitational interaction is
suppressed. To minimize the force errors and the global impact of the
softening for distances larger than <math>\epsilon\ ,</math> finite size
kernels that ensure continuous derivatives of the force may be
employed (e.g., see Dehnen 2001). This strategy effectively
suppresses binary formation and strong gravitational interactions, but
at the price of altering the dynamics of the system.

The computational complexity of the numerical solution of a N-body
system for a fixed number of timesteps scales as <math>N^2\ ,</math> as the
evaluation of the force on each particle requires to take into account
contributions from all other members of the system. For example,
considering a single state of the art cpu core (speed <math>\approx 5</math>
GFlops), a single force evaluation through a direct method would
require about 1 second for a system with <math>N=10^4</math> particles (assuming
10 floating point operations per pair of particles) and more than a
week for <math>N=10^7\ .</math>

The arbitrarily large dynamic range in the unsoftened dynamics and the
expensive evaluation of the force have led to the development of a
wide number of numerical techniques aimed at obtaining a reliable
numerical solution with the minimum amount of computational resources,
depending on the astrophysical problem of interest. Here we start by
discussing the different astrophysical contexts where N-body
simulations are routinely employed and we then present the state of
the art techniques for these problems.