N-body simulations (gravitational)

Post-publication activity

Curator: Piet Hut

Contributors:

0.78 -

Michele Trenti

0.09 -

Eugene M. Izhikevich

0.04 -

Juan Pablo Carbajal

0.04 -

Alessandra Celletti

Tobias Denninger

Søren Bertil F. Dorch

Introduction

The underlying dynamics relevant in the astrophysical context for of asystem of N particles interacting gravitationally is typicallyNewton's law plus, in case, an external potential field (see howeverbelow for a discussion of N-body simulations in generalrelativity). The force $\vec{F}_i$ acting on particle$i$ of mass $m_i$ is:

\[\tag{1}\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), \]

where $G=6.67300 \cdot 10^{-11}$ $m^{3}$$kg^{-1}$ $s^{-2}$ is the gravitationalconstant, and $\phi_{ext}$ is the external potential. The problemis thus a set of non-linear second order ordinary differentialequations relating the acceleration $\partial^2 \vec{r_i} /\partial t^2 = \vec{F}_i /m_i$ with the position of all theparticles in the system.

Once a set of initial condition is specified (for example the initialpositions $\vec{r}_i$ and velocities $\vec{v}_i \equiv\partial \vec{r}_i / \partial t $ of all particles) it exists aunique solution, analytical only for up to two bodies, while larger Nrequire numerical integration (e.g. see Press et al. 2007). Howeverspecial care must be employed to ensure both accuracy andefficiency. In fact, the gravitational force(eq. (1) presents a singularity when the distance oftwo particles approaches 0, which can lead to arbitrarily largerelative velocities. In addition, given the non-linear nature of theequations, the singularities are movable, that is they depend on thespecific choice of initial conditions. In contrast, all singularitiesin linear ordinary differential equations are independent of initialconditions and thus easier to treat. Therefore constant timestepmethods are unable to guarantee a given accuracy in thecase of gravitational dynamics and lead to unphysical accelerationsduring close encounters, which in turn may create unbound stars. Ashared, adaptive timestep scheme can correctly follow a closeencounter, but the price is paid in terms of efficiency as all theother particles of the system are evolved on the timescale of theencounter, which may be several orders of magnitude smaller than theglobal timescale, resulting essentially in a freezing of thesystem.

The singularity may be avoided by introducing a softening length inEq. (1) (e.g. see Aarseth 1963), that is by modifying the gravitationalinteraction at small scales. For example:

\[ \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}} ,\]

where $\epsilon > 0$ is the softening length, that is atypical distance below which the gravitational interaction issuppressed. To minimize the force errors and the global impact of thesoftening for distances larger than $\epsilon\ ,$ finite sizekernels that ensure continuous derivatives of the force may beemployed (e.g., see Dehnen 2001). This strategy effectivelysuppresses binary formation and strong gravitational interactions, butat the price of altering the dynamics of the system.

The computational complexity of the numerical solution of a N-bodysystem for a fixed number of timesteps scales as $N^2\ ,$ as theevaluation of the force on each particle requires to take into accountcontributions from all other members of the system. For example,considering a single state of the art cpu core (speed $\approx 5$GFlops), a single force evaluation through a direct method wouldrequire about 1 second for a system with $N=10^4$ particles (assuming10 floating point operations per pair of particles) and more than aweek for $N=10^7\ .$

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

Astrophysical domains and timescales

N-body simulations are applied to a wide range of differentastrophysical problems so that the most appropriate technique to usedepends on the specific context, and in particular on the timescaleand collisionality of the problem.

Timescales, Equilibrium and Collisionality

A system of N particles interacting gravitationally with total mass Mand a reference dimension R (for example the radius containing half ofthe total mass) reaches a dynamic equilibrium state on a timescalecomparable to a few times the typical time $T_{cr}$ needed for a particle tocross the system $(T_{cr} \approx 1/\sqrt{GM/R^3})\ .$ This is theresponse time needed to settle down to virial equilibrium, that is$2K/|W|=1\ ,$ where $K$ is the kinetic energy of the system\[K=1/2 \sum_{i=1,N} m_i |\vec{v}_i|^2 \ ,\] and W is its potentialenergy\[W = - 1/2 \sum_{i \ne j} G m_i m_j/|\vec{r}_i-\vec{r}_j|\] (assuming no external field). If the system is initially out of equilibrium,this is reached through mixing in phase space due to fluctuations ofthe gravitational potential, a process called violentrelaxation (Lynden-Bell 1967).

Once the system is in dynamic equilibrium a long term evolution ispossible, driven by two-body relaxation. Energy is slowly exchangedbetween particles and the system tends to evolve toward thermodynamicequilibrium and energy equipartition. The timescale $T_{rel}$ for this processdepends on the number of particles and on the geometry of the system\[T_{rel} \propto N/log(0.11 N) T_{cr}\] (e.g. see Spitzer1987). N-body systems such as galaxies and dark matter halos have arelaxation time much longer than the life of the Universe and are thusconsidered collisionless systems. Smaller systems, such asglobular and open clusters, are instead collisional, as the relaxationtime is shorter than their age. Two body relaxation is also suppressedwhen one particle in the system dominates the gravitational potential,such as in the case of solar system dynamics, where planets areessentially quasi-test particles.

Close encounters between three or more particles not only contributeto energy exchange, but can also lead to the formation of boundsubsystems (mainly binaries). The formation and evolution of a binarypopulation is best followed through direct, unsoftened, N-bodytechniques.

A self-gravitating N-body system made of single particles has anegative specific heat, that is it increases its kinetic energy as aresult of energy losses (Lynden-Bell & Wood 1968). This is a consequence ofthe virial theorem and qualitatively it is analogous to theacceleration of a Earth artificial satellite in presence ofatmospheric drag. A negative specific heat system is thermodynamically unstable and overthe two body relaxation timescale it evolves toward a gravothermalcollapse, creating a core-halo structure, where the core progressivelyincreases its concentration, fueling an overall halo expansion. Thecollapse is eventually halted once three body interactions lead to theformation of binaries. The so called "core collapsed globularclusters" are considered to be formed as a result of this mechanism.

Mean field approach: the Boltzmann equation

A system of N particles interacting gravitationally defines a 6N+1dimensional phase space given by the N position and velocity vectorsassociated to each particle at each time t. The solution of the N-bodyproblem defines a trajectory in this phase space. If the number ofparticles is large enough, that is if the two body relaxation time islong compared to the time-frame one is interested in, then astatistical description of the problem is possible. This allows us to pass from a 6N+1 to a 6+1 dimensionphase space. The idea is to construct a mean fielddescription of the dynamical system in terms of a single particledistribution function $f(\vec{x},\vec{v},t) \ ,$ where$f(\vec{x},\vec{v},t) d^3x d^3v $ is proportional to theprobability of finding a particle in a 6D element of volume $d^3x d^3v $ centered around position $\vec{r}$ andvelocity $\vec{v}$ at time t. Within this simplified framework theknowledge of the distribution function uniquely defines all theproperties of the system. The dynamic is described by thecollisionless Boltzmann equation, which derives essentially from theLiouville theorem:

\[\frac{D f}{D t} = \frac{\partial f}{\partial t} + \vec{v} \cdot \frac{\partial f}{\partial \vec{x}} - \frac{\partial \phi_T}{\partial \vec{x}} \cdot \frac{\partial f}{\partial \vec{v}}= 0,\]

where the total potential field $\phi_T = \phi_{ext}(\vec{x},t)+\phi(\vec{x},t)$ is the sum of an external potential plus theself-consistent field $\phi(\vec{x},t)$ defined from thedistribution function itself through the solution of the Poissonequation:

Mean Field Approach: analogies and differences with fluid dynamics

The velocity moments of the Boltzmann Equation define a set ofequations known as the Jeans Equations (e.g. Binney & Tremaine2007). The first three equations of the set are formally identical tothe Navier-Stokes equations for a self-gravitating gas and, like inthe fluid-dynamics analogy, express the conservation of mass, momentumand energy. Therefore the numerical algorithms developed to follow thedynamics of N-body systems find a wide application also in the contextof fluid-dynamics, with one important example being the SmoothedParticle Hydrodynamics (SPH) method (Gingold & Monaghan 1977). Thefundamental difference between the two cases is that the Jeansequations are derived in the limit of a collisionless system, while theNavier-Stokes equations assume a highly collisional system, with themean free path of a particle approaching zero. For fluids, this leads tothe definition of an equation of state, which closes the Navier-Stokesequations. The Jeans Equations are instead an infinite open set, wherethe n-th velocity moment depends on the n-th+1 moment.

Astrophysical domains

Based on the previous considerations about collisionality andtimescales, four main astrophysical domains for N-body simulations canbe identified, each requiring a different numerical technique toguarantee maximum performance and accuracy:

Celestial mechanics (solar and extrasolar planetary systems). Herea single body dominates the gravitational field and all the otherobjects move almost like test particles, subject to reciprocalperturbations. In this framework very high accuracy is required tocorrectly evaluate the perturbative terms and to avoid being dominatedby numerical noise such as time discretization and round-offs errors.

Dense Stellar systems, such as open clusters and globularclusters. These collisional systems made of components of roughlyequal mass present a rich dynamics, with multiple close encounters ofstars. The evolution requires to be followed on a relaxation timescalewith a correct description of the short range interactions.

Sphere of influence of a massive BH at the center of a stellarsystem. The sphere of influence of a BH is the volume within which thegravity of the BH dominates over that of the other particles. Thesituation resembles that of solar system dynamics, but here given thevery high density of stars two body encounters are frequent, makingthe problem a difficult hybrid between the two previous cases. Inaddition, Post Newtonian physics may need to be included if highaccuracy is required in the proximity of the BH.

Galaxy dynamics and cosmology. Galaxies, and especially darkmatter halos, are constituted by a very large number of particles, sothat their dynamics can be well described in terms of a meanfield. Close encounters are not important and softening is usuallyemployed in these N-body simulations to avoid the unphysical formation of binaries. Within this class, Self-Interacting Dark Matter Particles need a special mention: if dark matter halos aremade of Weakly Interacting Massive Particles, then their dynamics canbe modified by non-gravitational self-interactions, especially effective at the centerof cuspy dark halos. The dynamics of such a system is described by theCollisional Boltzmann Equation, which can be approximately solved using Fokker-Plankmethods.

Figure 1: Visualization of a simulation of galaxy collision using stars (blue) and SPH gas (yellow) particles. Credit: Frank Summers - STScI, reproduced with permission. Original images and movies available at http://terpsichore.stsci.edu/~summers/viz

Newtonian gravity: methods

The history of N-body simulations starts with a pioneering attempt byHolmberg (1941), who followed the evolution of a 37 particle system,where the force was calculated using lightbulbs and galvanometers(taking advantage of the same $r^{-2}$ scaling ofelectromagnetic and gravitational interactions). Computer simulationsstarted in the early sixties using up to 100 particles (e.g. see vonho*rner 1960 and Aarseth 1963) and had their full bloom in theeighties with the development of fast and efficientalgorithms to deal with collisionless systems, such as particle-meshcodes (see Hockney & Eastwood 1988 and references therein) and thetree method (Barnes & Hut 1986). At the same time regularization techniqueswere developed to deal with close encounters and binary dynamics in the case ofdirect simulations of a collisional system (e.g. see Aarseth's NBODY-Xcode series based on KS and chain regularization - Aarseth 2003and references therein). These algorithm advancements were coupled withtremendous progresses in the hardware, with the cpu speed growingexponentially. In addition to parallelization of serial codes, thefield advanced also thanks to special purpose hardware, such as theGRAPE (Makino et al. 1997). Today's (2008) N-body simulations are performed with upto $N=10^5$ (e.g. see Baumgardt & Makino 2003) for direct integration codes over a two-body relaxation timescale and upto $N=10^{10}$ for collisionless dynamics/cosmology (e.g. see theMillennium Run - Springel et al. 2005). In the context of planetary dynamics, self-gravitating systems of disk/ring particles with $N\approx 10^6 $ can be followed over hundreds of dynamical times (e.g. Richardson et al. 2000). Major breakthroughs are also expected in the nearfuture thanks both to the next generation GRAPE-DR and to doubleprecision graphic processing units, which provide extremely costcompetitive high performance computing capabilities.

Figure 2: Visualization from the Millennium Run Simulation, the largest published cosmological simulation to date. The figure shows the projected density field for a 15 Mpc/h thick slice of the redshift z=0 output. The overlaid panels zoom in by factors of 4 in each case, enlarging the regions indicated by the white squares. Credit: Volker Springel and the Virgo Consortium team, reproduced with permission. Original image available at http://www.mpa-garching.mpg.de/galform/virgo/millennium/index.shtml

Direct methods

Direct methods do not introduce approximations in the solution of theequations of motions and thus deliver the highest accuracy at theprice of the longest computation time, of order$O(N^2)$ per timestep. Integration is performed using adaptive(individual) timesteps and commonly a fourth order Hermiteintegrator. Close encounters and bound subsystems are treated exactlyin terms of Kustaanheimo-Steifel transformations. These essentially consist in transformations of coordinates using a perturbativeapproach over the analytical two body solution. If more than twoparticles are strongly interacting between each other, a chain regularization strategy(Mikkola 1990) can be used, which consists in recasting the problem interms of a series of separate Kustaanheimo-Steifel interactions. Astate of the art, publicly available, serial direct N-body integratoris Aarseth's NBODY6. Even with this specialized software, the numberof particles that can be effectively followed for timescalescomparable to the Hubble time is limited. For example, if one isinterested in the dynamical evolution of globular clusters, currentlyabout $N=20000$ is the practical limit for a serial run, assuch a run takes about $1000$ cpu hours. A run with$10^6$ particles carried out for a similar number ofrelaxation times $T_{rel}$ would require about$10^8$ cpu hours. The algorithm can be parallelized, but inpractice load imbalances may saturate the gain in efficiency,so some of the most cpu demanding simulations have been carried out on special purposehardware, such as the GRAPE, where the chip architecture has beenoptimized to compute gravitational interactions, thus deliveringTeraflops performance.

Tree codes

The tree code method (Barnes & Hut 1986) provides a fast, generalintegrator for collisionless systems, when close encounters are notimportant and where the force contributions from very distantparticles does not need to be computed at very high accuracy. In fact,with a tree code, small scale, strong interactions are typically softened (but see McMillan & Aarseth 1993), whilethe potentials due to distant groups of particles are approximated bymultipole expansions about the group centers of mass. The resultingcomputation time scales as $O(N log(N))$ but the approximations introduces small force errors. The long-range force errors are controlled by a singleparameter (the opening angle) that determines how small and distant agroup of particles must be to use the approximation. This strategyworks well to keep the average force error low, but a worst casescenario analysis highlights that unbound errors can arise for rare,but astrophysically reasonable configurations, such as that of theclassic "exploding galaxy" (Salmon & Warren 1994). In addition, forceerrors from the tree code may lead to violation of momentumconservation. Typical implementations of the tree code expand thepotentials to quadrupole order and construct a tree hierarchy ofparticles using a recursive binary splitting algorithm. The tree doesnot need to be recomputed from scratch at every timestep, savingsignificant cpu time. Systems with several hundred thousands ofcollisionless particles can be easily simulated on a GFlopsworkstation for a Hubble time using this method.

Fast Multipole Methods

A standard tree code implementation does not take advantage of thefact that nearby particles will be subject to a similar accelerationdue to distant groups of particles. The Fast Multipole Method(Greengard & Rokhlin 1987) exploit this idea and uses a multipoleexpansion to compute the force from a distant source cell within asink cell. This additional approximation of the gravitational interactionwas claimed to reduce the complexity from $O(N log(N))$ to$O(N)\ ,$ but the exact scaling seems implementationdependent and has been debated in the literature (e.g. see Dehnen 2000and references therein). One advantage of the fast multipole method isthat the symmetry in the treatment of sink and source cells withrespect to the multipole expansion can guarantee an exact conservationof the momentum. Recent successful implementations of fast multipolecodes or hybrids with tree code scheme, include Dehnen's Cartesianexpansion scheme (the GyrfalcON code- Dehnen 2000) and PKDGRAV (Stadel 2001).

Particle-mesh codes

The particle mesh method represents another route to speed up directforce evaluation for collisionless systems. In this case thegravitational potential of the system is constructed over a gridstarting from the density field and by solving the associated Poissonequation. Particles do not interact directly between each other butonly through a mean field. The method essentially softens thegravitational interactions at small scales, that is below the celllength. The density field is constructed using a kernel to split themass of the particles to the grid cells around the particleposition. The simplest choice is to assign all the mass to a singlecell, but this leads to significant force fluctuations, which can bereduced using a cloud in cell (8 points) or a triangular shaped cloud(27 points) kernel. The Poisson equation is typically solved using a FastFourier Transform, but other grid methods such as successive overrelaxation can also be used - e.g. see Bodenheimer et al. (2007). The deriving force, defined on the grid, is thenassigned back to the particles using the same kernel employed for thedensity field construction, in order to avoid spurious selfforces. The method achieves a linear complexity in the number ofparticles and ($O(N_g log(N_g)$) in the number of gridcells (this latter scaling is that of the FFT method). The price topay is in terms of short range accuracy as the force is a poorapproximation of Newton's law up to several grid spacing ofdistance.

Adaptive Mesh Refinement method

The dynamic range of particle-mesh codes can be increased by using anadaptive rather than a static grid to solve the Poisson Equation. Inthe Adaptive Mesh Refinement (AMR) method the grid elements areconcentrated where a higher resolution is needed, for example aroundthe highest density regions. One possibility to obtain an adaptiveresolution is to first construct a low-resolution solution of thePoisson Equation and then to progressively refine regions where thelocal truncation error (estimated through the Richardsonextrapolation) is highest. A multigrid structure needs to take intoaccount issues such as matching the solution at the grid interfaces.AMR codes are well suited for cosmological simulations (e.g. see the ENZO code, Bryan & Norman1998).

Self consistent field methods

A variant over the Particle Mesh code is the expansion of the densityand potential of the system in terms of a basis of orthogonaleigenfunctions. Clutton-Brock (1972) was one of the first to applythis idea in stellar dynamics, while a modern implementation is thatof Hernquist & Ostriker (1992). This method guarantees at fixedcomputational resources a higher accuracy than the tree code and theparticle mesh algorithms, provided that the set of basis function isappropriately selected. This limits in practice a general applicationof the method, which remains however very competitive for the study ofthe dynamical stability of collisionless systems constructed fromdistributions functions models.

P3M and PM-Tree codes

In order to increase the force resolution of particle mesh codes ithas been proposed to couple a mean field description on large scaleswith a direct, softened, treatment of the gravitational interactionson distances of the order of or below a few grid spacing (e.g. see Springel 2005). Thismethod is called $P^3M\ :$ Particle-Particle-Particle-Meshand efficiently increases the dynamic range of the parent PMalgorithm. However in presence of strong clustering a large number ofparticles will interact directly between each other, slowing downsignificantly the computation to $O(N^2)\ .$ This problem canbe resolved by using adaptive meshes, so that the spatial resolutionis refined in regions of high density. Adaptive $P^3M$codes have a computational cost which scales as $O(Nlog(N))\ ,$ like in a tree code. Finally another possibility is toresort to a tree code for the short range force evaluation leading toa hybrid PM-Tree scheme. These methods are generally extremely wellsuited for cosmological simulations, for example see Gadget2 (Springel 2005).

Celestial mechanics codes

Computational Celestial Mechanics refers to a series of methodstargeted at studying the dynamics of small N systems ($ N\lesssim 20 $). The smallest non trivial N is N=3, that is thethree body problem, which has many applications ranging from spaceflight to planets satellite motions and to binary-single starsencounters. Celestial mechanics requires extremely high precisiongiven the chaotic nature of the N-body problem. Numerical methods arebased on the use of local system of coordinates, to fight round-offerrors in systems with a wide dynamic range, such as in the study ofstar-planet-satellite problems, as well as on the variationalequations formalism and on perturbation theory to take advantage ofthe analytical, unperturbed motion of planets in the gravitation fieldof their star (e.g. see Beutler 2005). In this context symplectic integrators are widely used (e.g. see Wisdom & Holman 1991; Leimkuhler & Reich 2005).

Mean Field Methods

As an alternative to particle based N-body methods, the dynamics of asystem of particles interacting gravitationally can be followed bysolving the time dependent Boltzmann Equation coupled with theself-consistent Poisson equation.

Grid based solvers for the Collisionless Boltzmann Equation

This approach can take advantage of standard computational methodsdeveloped to solve partial differential equations, such as successiveover-relaxation and conjugate gradient methods. However it requires tosolve a highly dimensional (6D+time) non-linear system of partialdifferential equations. In general, the bottleneck is thus the verylarge amount of memory needed (for example, Terabites just to have amoderate resolution grid with 100 elements in each dimension). Howeverthis method is competitive if the astrophysical problem of interestpresents symmetries that reduce the number of dimensions needed in themodel. For example, in the case of globular cluster dynamics a verygood approximation can be obtained via a 3 dimensional model byassuming spherical symmetry in the position space (1D) and radialanisotropy in the velocity space (2D).

Fokker-Planck and Monte Carlo methods

These methods solve the collisional Boltzmann equation starting from agiven distribution function and by following test particles in the sixdimensional position-velocity phase space. At each timestep thevelocity of the particles is perturbed by random fluctuationsaccordingly to the assumed form for the collision operator $C[f]\ ,$ which depends on computed cross sections for two, three andfour body encounters. The complexity of Monte Carlo codes is linearwith the number of particles and thus a realistic number of particlescan be used for simulations of collisional systems with $ N>10^5$ with a serial code. The method is ideal forexploring grids of initial conditions, after proper validationthrough comparison with direct integration (e.g. see Heggie et al. 2006).

Beyond Newton: strong gravitational fields

In presence of a strong gravitational field, such as that in theproximity of the event horizon of a black hole, N-body simulationscannot be based on Newtonian physics, but must take intoaccount a general relativity framework. As a numerical solution of theEinstein equation is extremely challenging, Post-Newtonianapproximations are used when the gravitational field does notdeviate too much from the Newtonian case. Post-Newtonian correctionsare typically good enough to treat most astrophysical problems of thedynamics of stars around a black hole. A full general relativityframework is only required to study the merging and gravitationalwaves emission of two black-holes (e.g. see Baker et al. 2006).

Hardware

An alternative approach to increase the efficiency of numericalsolution of the N-body problem is to optimize the hardware. For directsimulations this approach can be very effective, thanks to the factthat the bottle neck of computation is just the evaluation of thegravitational force, which has a very simple expression. Along thisroute the GRAPE (GRavityPipE) concept has been extremelyeffective. The basic idea is to optimize a hardware pipeline tocompute $ (\vec{r_i}-\vec{r_j})/| \vec{r_i}-\vec{r_j}|^3\ .$ This special purpose hardware can then be interfaced with ageneral purpose computer, which takes care of all the other numericaloperations required to solve the equations of motions. With theGRAPE-6, the largest simulation on a collisional timescale publishedto date has N=131028 (Baumgardt & Makino 2003).

Another recent promising hardware development is the possibility touse Graphic Cards (GPUs) to carry out the cpu intensive forceevaluation. The performance of current generation of GPUs appears tobe superior (in terms of Flops/$ ratio) to that of the GRAPE6 series(Portegies-Zwart et al. 2007).

Simulation environments

In addition to the availability of stand-alone codes, several software environments have been created that contain various tools to set up initial conditions, run simulations, and analyze and visualize their results. Some examples are NEMO, Starlab, ACS, MUSE.

References

Aarseth, S. 1963, MNRAS, 126, 223
Aarseth, S. 2003, "Gravitational N-Body Simulations: Tools and Algorithms", Cambridge University Press
Baker, J.G. et al. 2006, ApJ, 653, 93
Barnes, J.E. and Hut, P. 1986, Nature, 324, 466
Baumgardt, H. and Makino, J. 2003, MNRAS, 340, 227
Beutler G. 2005, "Methods of Celestial Mechanics", Springer
Binney J. & Tremaine S. 1987, "Galactic Dynamics", Princeton University Press
Bodenheimer, P., Laughlin, G.P., Rozyczka, M. and Yorke, H.W. 2007, "Numerical Methods in Astrophysics: An Introduction", Taylor & Francis
Bryan G.L. and Norman, M.L. 1998, ApJ, 495, 80
Clutton-Brock, M. 1972, Ap&SS, 16, 101
Dehnen, W. 2001, MNRAS, 324, 273
Dehnen, W. 2000, ApJL, 536, 39
Gingold, R.A. and Monaghan, J.J. 1977, MNRAS, 181, 375
Greengard, L. & Rokhlin, V. 1987, J. comput. Phys., 73, 325
Heggie, D.C., Trenti, M. and Hut, P. 2006, MNRAS, 368, 677
Hernquist, L and Barnes, J.E. 1990, ApJ, 349, 562
Hockney, R.W. and Eastwood, J.W. 1988, "Computer Simulation Using Particles", Taylor & Francis
Holmberg, E. 1941, ApJ, 94, 385
Leimkuhler, B. and Sebastian R. 2005, "Simulating Hamiltonian Dynamics", Cambridge University Press
Lynden-Bell, D. 1967, MNRAS, 136, 101
Lynden-Bell, D. and Wood, R. 1968, MNRAS, 138, 495
Makino, J., f*ckushige, T., Koga, M. and Namura, K. 2003, PASJ, 55, 1163
Portegies-Zwart, S.F., Belleman, R.G. and Geldof, P.M. 2007, New Astronomy, 12, 641
Press, W.H., Teukolsky, S.A., Vetterling, W.T. and Flannery B.P. 2007, "Numerical Recipes", Cambridge University Press
Richardson, D.C., Quinn, T., Stadel, J. and Lake, G. 2000, Icarus, 143, 45
Salmon J.K. and Warren M.S. 1994, J. Comp. Phys., 111, 136.
Spitzer, L. 1987, "Dynamical Evolution of Globular Clusters", Princeton University Press
Springel, V. 2005, MNRAS, 364, 1105
Springel, V. et al. 2005, Nature, 435, 629
Stadel J. 2001, PhD. Thesis, University of Washington
von ho*rner, S. 1960, Z. Astrophys. 50, 184
Wisdom, J. and Holman, M. 1991, AJ, 102, 1528

Internal references

Teviet Creighton and Richard H. Price (2008) Black holes. Scholarpedia, 3(1):4277.

James M. Stone (2007) Computational astrophysics. Scholarpedia, 2(10):2419.

James Meiss (2007) Dynamical systems. Scholarpedia, 2(2):1629.

Eugene M. Izhikevich (2007) Equilibrium. Scholarpedia, 2(10):2014.

Kendall E. Atkinson (2007) Numerical analysis. Scholarpedia, 2(8):3163.

Philip Holmes and Eric T. Shea-Brown (2006) Stability. Scholarpedia, 1(10):1838.

David H. Terman and Eugene M. Izhikevich (2008) State space. Scholarpedia, 3(3):1924.

Alain Chenciner (2007) Three body problem. Scholarpedia, 2(10):2111.

External links

Open source codes

Aarseth's direct integration codes: http://www.ast.cam.ac.uk/~sverre/web/pages/nbody.htm

ACS, a collection of tools and introductory texts: http://www.artcompsci.org/

ENZO, a cosmological AMR code: http://lca.ucsd.edu/portal/software/enzo

Gadget2, a cosmological PM-tree+SPH code (massively parallel): http://www.mpa-garching.mpg.de/gadget/

Mercury (a mixed variable symplectic integrator code for planetary dynamics): http://www.arm.ac.uk/~jec/home.html

MUSE, a software framework for simulations of dense stellar systems: http://muse.li/

NEMO collection (includes particle-grid and tree codes): http://bima.astro.umd.edu/nemo/

Starlab (including the direct integration Kira code): http://www.ids.ias.edu/~starlab/starlab.html

Pictures and movies from N-body simulations

Cosmological Simulations and Galaxy Interactions (includes movies and snapshots of the Millennium Run): http://www.mpa-garching.mpg.de/galform/data_vis/index.shtml

Cosmological Simulations, Galaxy Interactions, Globular Cluster evolution, Outreach movies: http://terpsichore.stsci.edu/~summers/viz/

N-body simulations (gravitational) - Scholarpedia (2024)

Contents