In physics and astronomy, an N-body simulation is a simulation of a dynamical system of particles, usually under the influence of physical forces, such as gravity (see n-body problem). N-body simulations are widely used tools in astrophysics, from investigating the dynamics of few-body systems like the Earth-Moon-Sun system to understanding the evolution of the large-scale structure of the universe. In physical cosmology, N-body simulations are used to study processes of non-linear structure formation such as galaxy filaments and galaxy halos from the influence of dark matter. Direct N-body simulations are used to study the dynamical evolution of star clusters.
Nature of the particles
The
'particles' treated by the simulation may or may not correspond to
physical objects which are particulate in nature. For example, an N-body
simulation of a star cluster might have a particle per star, so each
particle has some physical significance. On the other hand, a simulation
of a gas cloud cannot afford to have a particle for each atom or molecule of gas as this would require on the order of 1023 particles for each mole of material (see Avogadro constant), so a single 'particle' would represent some much larger quantity of gas (often implemented using Smoothed Particle Hydrodynamics).
This quantity need not have any physical significance, but must be
chosen as a compromise between accuracy and manageable computer
requirements.
Direct gravitational N-body simulations
In direct gravitational N-body simulations, the equations of motion of a system of N
particles under the influence of their mutual gravitational forces are
integrated numerically without any simplifying approximations. These
calculations are used in situations where interactions between
individual objects, such as stars or planets, are important to the
evolution of the system.
The first direct N-body simulations were carried out by Erik Holmberg at the Lund Observatory
in 1941, determining the forces between stars in encountering galaxies
via the mathematical equivalence between light propagation and
gravitational interaction: putting light bulbs at the positions of the
stars and measuring the directional light fluxes at the positions of the
stars by a photo cell, the equations of motion can be integrated with effort. The first purely calculational simulations were then done by Sebastian von Hoerner at the Astronomisches Rechen-Institut in Heidelberg, Germany. Sverre Aarseth at the University of Cambridge (UK) has dedicated his entire scientific life to the development of a series of highly efficient N-body
codes for astrophysical applications which use adaptive (hierarchical)
time steps, an Ahmad-Cohen neighbour scheme and regularization of close
encounters. Regularization is a mathematical trick to remove the
singularity in the Newtonian law of gravitation for two particles which
approach each other arbitrarily close. Sverre Aarseth's codes are used
to study the dynamics of star clusters, planetary systems and galactic
nuclei.
General relativity simulations
Many simulations are large enough that the effects of general relativity in establishing a Friedmann-Lemaitre-Robertson-Walker cosmology are significant. This is incorporated in the simulation as an evolving measure of distance (or scale factor) in a comoving coordinate system, which causes the particles to slow in comoving coordinates (as well as due to the redshifting of their physical energy). However, the contributions of general relativity and the finite speed of gravity
can otherwise be ignored, as typical dynamical timescales are long
compared to the light crossing time for the simulation, and the
space-time curvature induced by the particles and the particle
velocities are small. The boundary conditions of these cosmological
simulations are usually periodic (or toroidal), so that one edge of the
simulation volume matches up with the opposite edge.
Calculation optimizations
N-body simulations are simple in principle, because they involve merely integrating the 6N ordinary differential equations defining the particle motions in Newtonian gravity. In practice, the number N of particles involved is usually very large (typical simulations include many millions, the Millennium simulation included ten billion) and the number of particle-particle interactions needing to be computed increases on the order of N2,
and so direct integration of the differential equations can be
prohibitively computationally expensive. Therefore, a number of
refinements are commonly used.
Numerical integration is usually performed over small timesteps using a method such as leapfrog integration.
However all numerical integration leads to errors. Smaller steps give
lower errors but run more slowly. Leapfrog integration is roughly 2nd
order on the timestep, other integrators such as Runge–Kutta methods can have 4th order accuracy or much higher.
One of the simplest refinements is that each particle carries
with it its own timestep variable, so that particles with widely
different dynamical times don't all have to be evolved forward at the
rate of that with the shortest time.
There are two basic approximation schemes to decrease the computational time for such simulations. These can reduce the computational complexity to O(N log N) or better, at the loss of accuracy.
Tree methods
In tree methods, such as a Barnes–Hut simulation, an octree
is usually used to divide the volume into cubic cells and only
interactions between particles from nearby cells need to be treated
individually; particles in distant cells can be treated collectively as a
single large particle centered at the distant cell's center of mass (or
as a low-order multipole
expansion). This can dramatically reduce the number of particle pair
interactions that must be computed. To prevent the simulation from
becoming swamped by computing particle-particle interactions, the cells
must be refined to smaller cells in denser parts of the simulation which
contain many particles per cell. For simulations where particles are
not evenly distributed, the well-separated pair decomposition methods of Callahan and Kosaraju yield optimal O(n log n) time per iteration with fixed dimension.
Particle mesh method
Another possibility is the particle mesh method in which space is discretised on a mesh and, for the purposes of computing the gravitational potential,
particles are assumed to be divided between the nearby vertices of the
mesh. Finding the potential energy Φ is easy, because the Poisson equation
where G is Newton's constant and is the density (number of particles at the mesh points), is trivial to solve by using the fast Fourier transform to go to the frequency domain where the Poisson equation has the simple form
where is the comoving wavenumber and the hats denote Fourier transforms. The gravitational field can now be found by multiplying by
and computing the inverse Fourier transform (or computing the inverse
transform and then using some other method). Since this method is
limited by the mesh size, in practice a smaller mesh or some other
technique (such as combining with a tree or simple particle-particle
algorithm) is used to compute the small-scale forces. Sometimes an
adaptive mesh is used, in which the mesh cells are much smaller in the
denser regions of the simulation.
Special-case optimizations
Several different gravitational perturbation algorithms are used to get fairly accurate estimates of the path of objects in the solar system.
People often decide to put a satellite in a frozen orbit.
The path of a satellite closely orbiting the Earth can be accurately
modeled starting from the 2-body elliptical orbit around the center of
the Earth, and adding small corrections due to the oblateness of the Earth, gravitational attraction of the Sun and Moon, atmospheric drag, etc.
It is possible to find a frozen orbit without calculating the actual path of the satellite.
The path of a small planet, comet, or long-range spacecraft can
often be accurately modeled starting from the 2-body elliptical orbit
around the sun, and adding small corrections from the gravitational
attraction of the larger planets in their known orbits.
Some characteristics of the long-term paths of a system of
particles can be calculated directly. The actual path of any particular
particle does not need to be calculated as an intermediate step. Such
characteristics include Lyapunov stability, Lyapunov time, various measurements from ergodic theory, etc.
Two-particle systems
Although
there are millions or billions of particles in typical simulations,
they typically correspond to a real particle with a very large mass,
typically 109 solar masses. This can introduce problems with short-range interactions between the particles such as the formation of two-particle binary
systems. As the particles are meant to represent large numbers of dark
matter particles or groups of stars, these binaries are unphysical. To
prevent this, a softened Newtonian force law is used, which does not
diverge as the inverse-square radius at short distances. Most
simulations implement this quite naturally by running the simulations on
cells of finite size. It is important to implement the discretization
procedure in such a way that particles always exert a vanishing force on
themselves.
Incorporating baryons, leptons and photons into simulations
Many simulations simulate only cold dark matter, and thus include only the gravitational force. Incorporating baryons, leptons and photons
into the simulations dramatically increases their complexity and often
radical simplifications of the underlying physics must be made. However,
this is an extremely important area and many modern simulations are now
trying to understand processes that occur during galaxy formation which could account for galaxy bias.
Computational complexity
Reif et al. prove that if the n-body reachability problem is defined as follows – given n
bodies satisfying a fixed electrostatic potential law, determining if a
body reaches a destination ball in a given time bound where we require a
poly(n) bits of accuracy and the target time is poly(n) is in PSPACE.
On the other hand, if the question is whether the body eventually reaches the destination ball, the problem is PSPACE-hard. These bounds are based on similar complexity bounds obtained for ray tracing.