View source for N-body simulation

{{short description|Simulation of a dynamical system of particles}}
{{About|a topic in physics|the automobile platform|GM N platform}}
{{Hatnote|For engineering problems and simulations involving many components, see [[Multibody system]] and [[Multibody simulation]].}}
{{DISPLAYTITLE:''N''-body simulation}}
{{Use dmy dates|date=May 2020}}

[[File:Galaxy cluster sim.png|thumb|An ''N''-body simulation of the cosmological formation of a cluster of galaxies in an expanding universe|upright=1.36]]
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 problem]] for other applications). ''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]].<ref name=Trenti>{{cite journal| first1 = Michele | last1 = Trenti | first2 = Piet | last2 = Hut|title=N-body simulations (gravitational)| year = 2008 |journal=Scholarpedia| volume = 3 | issue = 5 | pages = 3930 | doi = 10.4249/scholarpedia.3930 | bibcode = 2008SchpJ...3.3930T | doi-access = free }}</ref> In [[physical cosmology]], ''N''-body simulations are used to study processes of non-linear [[structure formation]] such as [[galaxy filament]]s and [[galaxy halo]]s from the influence of [[dark matter]]. Direct ''N''-body simulations are used to study the dynamical evolution of [[star cluster]]s.

==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 [[interstellar cloud|gas cloud]] cannot afford to have a particle for each atom or molecule of gas as this would require on the order of {{val|e=23}} 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.

==Dark matter simulation==

[[Dark matter]] plays an important role in the formation of galaxies. The time evolution of the density f (in phase space) of dark matter particles, can be described by the collisionless [[Boltzmann equation]]

<math display="block">\frac{d f}{d t} = \frac{\partial f}{\partial t}+ \mathbf{v} \cdot\nabla f -  \frac{\partial f}{\partial \mathbf{v}} \cdot\nabla\Phi </math>

In the equation,<math>\mathbf{v}</math> is the velocity, and Φ is the gravitational potential given by [[Poisson's Equation]]. These two coupled equations are solved in an expanding background Universe, which is governed by the [[Friedmann equations]], after determining the initial conditions of dark matter particles. The conventional method employed for initializing positions and velocities of dark matter particles involves moving particles within a uniform Cartesian lattice or a glass-like particle configuration.<ref name="Oxford University Press (OUP)1995">{{cite journal | last1 = C.M.Baugh | last2 =  E.Gaztañaga | last3 =  G. Efstathiou | title = A comparison of the evolution of density fields in perturbation theory and numerical simulations - II. Counts-in-cells analysis | journal = Monthly Notices of the Royal Astronomical Society | date = 1995 | eissn = 1365-2966 | doi = 10.1093/mnras/274.4.1049 | doi-access = free | pmid = | url = | arxiv = astro-ph/9408057 }}</ref> This is done by using a linear theory approximation or a low-order [[perturbation theory]].<ref name="Jenkins2010">{{cite journal | last1 = Jenkins | first1 = Adrian | title = Second-order Lagrangian perturbation theory initial conditions for resimulations | journal = Monthly Notices of the Royal Astronomical Society | date = 21 April 2010 | volume = 403 | issue = 4 | pages = 1859–1872 | issn = 0035-8711 | eissn = 1365-2966 | doi = 10.1111/j.1365-2966.2010.16259.x | doi-access = free | pmid = | url = | arxiv = 0910.0258 | bibcode = 2010MNRAS.403.1859J }}</ref>

==Direct gravitational ''N''-body simulations==

[[File:N-body Simulation.ogv|thumb|N-body simulation of 400 objects with parameters close to those of [[Solar System]] planets]]
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 gravitational ''N''-body simulations were carried out by [[Erik Holmberg (astronomer)|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 {{tmath|O(N)}} effort.<ref name="Holmberg1941">{{Cite journal | last=Holmberg | first=Erik | title=On the Clustering Tendencies among the Nebulae. II. a Study of Encounters Between Laboratory Models of Stellar Systems by a New Integration Procedure | year = 1941 | journal=The Astrophysical Journal|volume=94|issue=3|pages=385–395|doi=10.1086/144344| bibcode=1941ApJ....94..385H | doi-access=free}}</ref> The first purely calculational simulations were then done by [[Sebastian von Hoerner]] at the [[Astronomical Calculation Institute (University of Heidelberg)|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.{{citation needed|date=November 2017}}

==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 (Universe)|scale factor]]) in a [[comoving coordinates|comoving coordinate]] system, which causes the particles to slow in comoving coordinates (as well as due to the [[redshift]]ing 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 6''N'' [[ordinary differential equation]]s 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 ''N''<sup>2</sup>, 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 theory|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 [[S. Rao Kosaraju|Kosaraju]] yield optimal O(''n''&nbsp;log&nbsp;''n'') time per iteration with fixed dimension.

===Particle mesh method===

Another possibility is the [[Particle Mesh|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 surrounding 2x2 vertices of the mesh. The potential energy Φ can be found with the [[Poisson equation]]

<math display="block">\nabla^2\Phi=4\pi G{\rho},\,</math>

where ''G'' is [[Newton's constant]] and <math>\rho</math> is the density (number of particles at the mesh points). The [[fast Fourier transform]] can solve this efficiently by going to the [[frequency domain]] where the Poisson equation has the simple form

<math display="block">\hat{\Phi}= -4\pi G\frac{\hat{\rho}}{k^2},</math>

where <math>\vec{k}</math> is the comoving wavenumber and the hats denote Fourier transforms. Since <math>\vec{g} = -\vec{\nabla}\Phi</math>, the gravitational field can now be found by multiplying by <math>-i\vec{k}</math> 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 [[figure of the Earth|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 10<sup>9</sup> [[solar mass]]es. This can introduce problems with short-range interactions between the particles such as the formation of two-particle [[binary (astronomy)|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 [[Softening|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.

===Softening===
'''Softening''' is a numerical trick used in N-body techniques to prevent numerical [[divergence]]s when a particle comes too close to another (and the [[force]] goes to infinity). This is obtained by modifying the regularized [[gravitational potential]] of each particle as

<math display="block">\Phi = - \frac{1}{\sqrt{r^2 + \epsilon^2}},</math>

(rather than 1/r) where <math>\epsilon</math> is the softening parameter. The value of the softening parameter should be set small enough to keep [[simulation]]s realistic.

==Results from ''N''-body simulations ==

''N''-body simulations give findings on the large-scale dark matter distribution and the structure of dark matter halos. According to simulations of cold dark matter, the overall distribution of dark matter on a large scale is not entirely uniform. Instead, it displays a structure resembling a network, consisting of voids, walls, filaments, and halos. Also, simulations show that the relationship between the concentration of halos and factors such as mass, initial fluctuation spectrum, and cosmological parameters is linked to the actual formation time of the halos.<ref name="NavarroFrenkWhite1997">{{cite journal | last1 = Navarro | first1 = Julio F. | last2 = Frenk | first2 = Carlos S. | last3 = White | first3 = Simon D. M. | title = A Universal Density Profile from Hierarchical Clustering | journal = The Astrophysical Journal | date = December 1997 | volume = 490 | issue = 2 | pages = 493–508 | issn = 0004-637X | eissn = 1538-4357 | doi = 10.1086/304888 | pmid = | url = | arxiv = astro-ph/9611107 | bibcode = 1997ApJ...490..493N }}</ref> In particular, halos with lower mass tend to form earlier, and as a result, have higher concentrations due to the higher density of the Universe at the time of their formation. Shapes of halos are found to deviate from being perfectly spherical. Typically, halos are found to be elongated and become increasingly prolate towards their centers. However, interactions between dark matter and [[baryons]] would affect the internal structure of dark matter halos. Simulations that model both dark matters and baryons are needed to study small-scale structures.

==Incorporating baryons, leptons and photons into simulations==

Many simulations simulate only [[cold dark matter]], and thus include only the gravitational force. Incorporating [[baryon]]s, [[lepton]]s and [[photon]]s 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 and Tate<ref>{{cite conference |author=John H. Reif |author2=Stephen R. Tate| citeseerx = 10.1.1.38.6242 |series=Lecture Notes in Computer Science | title = The Complexity of N-body Simulation |book-title=Automata, Languages and Programming | year = 1993 | pages = 162–176 }}</ref> 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 (physics)|ray tracing]].

== Example simulations ==

=== Common boilerplate code ===
The simplest implementation of N-body simulations where <math display="inline">n\geq 3</math> is a naive propagation of orbiting bodies; naive implying that the only forces acting on the orbiting bodies is the gravitational force which they exert on each other. In [[object-oriented programming]] languages, such as [[C++]], some [[boilerplate code]] is useful for establishing the fundamental mathematical structures as well as data containers required for propagation; namely [[State vector (navigation)|state vectors]], and thus [[Vector (mathematics and physics)|vectors]], and some fundamental object containing this data, as well as the mass of an orbiting body. This method is applicable to other types of N-body simulations as well; a simulation of point masses with charges would use a similar method, however the force would be due to attraction or repulsion by interaction of electric fields. Regardless, acceleration of particle is a result of summed force vectors, divided by the mass of the particle:

<math display="block">\vec{a} = \frac{1}{m} \sum\vec{F}</math>

An example of a programmatically stable and [[Scalability|scalable]] method for containing kinematic data for a particle is the use of fixed length arrays, which in optimised code allows for easy memory allocation and prediction of consumed resources; as seen in the following C++ code:
<syntaxhighlight lang="c++" line="1">
struct Vector3
{
    double e[3] = { 0 };

    Vector3() {}
    ~Vector3() {}

    inline Vector3(double e0, double e1, double e2)
    {
        this->e[0] = e0;
        this->e[1] = e1;
        this->e[2] = e2;
    }
};

struct OrbitalEntity
{
    double e[7] = { 0 };

    OrbitalEntity() {}
    ~OrbitalEntity() {}

    inline OrbitalEntity(double e0, double e1, double e2, double e3, double e4, double e5, double e6)
    {
        this->e[0] = e0;
        this->e[1] = e1;
        this->e[2] = e2;
        this->e[3] = e3;
        this->e[4] = e4;
        this->e[5] = e5;
        this->e[6] = e6;
    }
};
</syntaxhighlight>Note that {{Code|OrbitalEntity|C++}} contains enough room for a state vector, where:

*<math display="inline">e_{0} = x</math>,  the projection of the objects position vector in Cartesian space along <math>\left[1\;0\;0\right]</math>
*<math display="inline">e_{1} = y</math>,  the projection of the objects position vector in Cartesian space along <math>\left[0\;1\;0\right]</math>
*<math display="inline">e_{2} = z</math>,  the projection of the objects position vector in Cartesian space along <math>\left[0\;0\;1\right]</math>
*<math display="inline">e_{3} = \dot{x}</math>,  the projection of the objects velocity vector in Cartesian space along <math>\left[1\;0\;0\right]</math>
*<math display="inline">e_{4} = \dot{y}</math>,  the projection of the objects velocity vector in Cartesian space along <math>\left[0\;1\;0\right]</math>
*<math display="inline">e_{5} = \dot{z}</math>,  the projection of the objects velocity vector in Cartesian space along <math>\left[0\;0\;1\right]</math>

Additionally, {{Code|OrbitalEntity|C++}} contains enough room for a mass value.

=== Initialisation of simulation parameters ===
Commonly, N-body simulations will be systems based on some type of [[equations of motion]]; of these, most will be dependent on some initial configuration to "seed" the simulation. In systems such as those dependent on some gravitational or electric potential, the force on a simulation entity is independent on its velocity. Hence, to seed the ''forces'' of the simulation, merely initial positions are needed, but this will not allow propagation- initial velocities are required. Consider a planet orbiting a star- it has no motion, but is subject to gravitational attraction to its host star. As a time progresses, and time ''steps'' are added, it will gather velocity according to its acceleration. For a given instant in time, <math>t_{n}</math>, the resultant acceleration of a body due to its neighbouring masses is independent of its velocity, however, for the time step <math>t_{n+1}</math>, the resulting change in position is significantly different due the propagation's inherent dependency on velocity. In basic propagation mechanisms, such as the symplectic euler method to be used below, the position of an object at <math>t_{n+1}</math> is only dependent on its velocity at <math>t_{n}</math>, as the shift in position is calculated via

<math display="block">\vec{r}_{t_{n+1}}=\vec{r}_{t_{n}}+\vec{v}_{t_{n}} \cdot \Delta t</math>

Without acceleration, <math display="inline">\vec{v}_{t_{n}}</math> is static, however, from the perspective of an observer seeing only position, it will take two time steps to see a change in velocity.

A solar-system-like simulation can be accomplished by taking average distances of planet equivalent point masses from a central star. To keep code simple, a non-rigorous approach based on [[Orbital elements|semi-major axes]] and mean velocities will be used. [[Random-access memory|Memory space]] for these bodies must be reserved before the bodies are configured; to allow for scalability, a [[Memory allocation|malloc]] command may be used:
<syntaxhighlight lang="c++" line="1">
OrbitalEntity* orbital_entities = malloc(sizeof(OrbitalEntity) * (9 + N_ASTEROIDS));

orbital_entities[0] = { 0.0,0.0,0.0,        0.0,0.0,0.0,      1.989e30 };   // a star similar to the sun
orbital_entities[1] = { 57.909e9,0.0,0.0,   0.0,47.36e3,0.0,  0.33011e24 }; // a planet similar to mercury
orbital_entities[2] = { 108.209e9,0.0,0.0,  0.0,35.02e3,0.0,  4.8675e24 };  // a planet similar to venus
orbital_entities[3] = { 149.596e9,0.0,0.0,  0.0,29.78e3,0.0,  5.9724e24 };  // a planet similar to earth
orbital_entities[4] = { 227.923e9,0.0,0.0,  0.0,24.07e3,0.0,  0.64171e24 }; // a planet similar to mars
orbital_entities[5] = { 778.570e9,0.0,0.0,  0.0,13e3,0.0,     1898.19e24 }; // a planet similar to jupiter
orbital_entities[6] = { 1433.529e9,0.0,0.0, 0.0,9.68e3,0.0,   568.34e24 };  // a planet similar to saturn
orbital_entities[7] = { 2872.463e9,0.0,0.0, 0.0,6.80e3,0.0,   86.813e24 };  // a planet similar to uranus
orbital_entities[8] = { 4495.060e9,0.0,0.0, 0.0,5.43e3,0.0,   102.413e24 }; // a planet similar to neptune
</syntaxhighlight>
where {{Code|N_ASTEROIDS|c++}} is a variable which will remain at 0 temporarily, but allows for future inclusion of significant numbers of asteroids, at the users discretion. A critical step for the configuration of simulations is to establish the time ranges of the simulation, <math>t_{0}</math> to <math>t_\text{end}</math>, as well as the incremental time step <math>dt</math> which will progress the simulation forward:
<syntaxhighlight lang="c++" line="1">
double t_0 = 0;
double t = t_0;
double dt = 86400;
double t_end = 86400 * 365 * 10; // approximately a decade in seconds
double BIG_G = 6.67e-11; // gravitational constant
</syntaxhighlight>
The positions and velocities established above are interpreted to be correct for <math>t = t_{0}</math>.

The extent of a simulation would logically be for the period where <math>t_{0}\leq t < t_\text{end}</math>.

=== Propagation ===
An entire simulation can consist of hundreds, thousands, millions, billions, or sometimes trillions of time steps. At the elementary level, each time step (for simulations with particles moving due to forces exerted on them) involves

* calculating the forces on each body
* calculating the accelerations of each body (<math>\vec{a}</math>)
* calculating the velocities of each body (<math>\vec{v}_{n} = \vec{v}_{n-1} + \vec{a}_{n} \cdot \Delta t</math>
* calculating the new position of each body (<math>\vec{r}_{n+1} = \vec{r}_{n} + \vec{v}_{n} \cdot \Delta t</math>
The above can be implemented quite simply with a [[while loop]] which continues while <math>t</math> exists in the aforementioned range:
<syntaxhighlight lang="c++">
while (t < t_end)
{
    for (size_t m1_idx = 0; m1_idx < 9 + N_ASTEROIDS; m1_idx++)
    {       
        Vector3 a_g = { 0,0,0 };

        for (size_t m2_idx = 0; m2_idx < 9 + N_ASTEROIDS; m2_idx++)
        {
            if (m2_idx != m1_idx)
            {
                Vector3 r_vector;

                r_vector.e[0] = orbital_entities[m1_idx].e[0] - orbital_entities[m2_idx].e[0];
                r_vector.e[1] = orbital_entities[m1_idx].e[1] - orbital_entities[m2_idx].e[1];
                r_vector.e[2] = orbital_entities[m1_idx].e[2] - orbital_entities[m2_idx].e[2];

                double r_mag = sqrt(
                        r_vector.e[0] * r_vector.e[0]
                      + r_vector.e[1] * r_vector.e[1]
                      + r_vector.e[2] * r_vector.e[2]);

                double acceleration = -1.0 * BIG_G * (orbital_entities[m2_idx].e[6]) / pow(r_mag, 2.0);

                Vector3 r_unit_vector = { r_vector.e[0] / r_mag, r_vector.e[1] / r_mag, r_vector.e[2] / r_mag };

                a_g.e[0] += acceleration * r_unit_vector.e[0];
                a_g.e[1] += acceleration * r_unit_vector.e[1];
                a_g.e[2] += acceleration * r_unit_vector.e[2];
            }
        }

        orbital_entities[m1_idx].e[3] += a_g.e[0] * dt;
        orbital_entities[m1_idx].e[4] += a_g.e[1] * dt;
        orbital_entities[m1_idx].e[5] += a_g.e[2] * dt;
    }

    for (size_t entity_idx = 0; entity_idx < 9 + N_ASTEROIDS; entity_idx++)
    {
        orbital_entities[entity_idx].e[0] += orbital_entities[entity_idx].e[3] * dt;
        orbital_entities[entity_idx].e[1] += orbital_entities[entity_idx].e[4] * dt;
        orbital_entities[entity_idx].e[2] += orbital_entities[entity_idx].e[5] * dt;
    }
    
    t += dt;
}
</syntaxhighlight>
Focusing on the inner four [[rocky planet]]s in the simulation, the trajectories resulting from the above propagation is shown below:
[[File:Simulated trajectories of four rocky planets with dt 86400.gif|center|512x512px]]

==See also==
{{Portal|Physics}}
*{{Annotated link|Millennium Run}}
*{{Annotated link|Large-scale structure of the cosmos}}
*{{Annotated link|GADGET}}
*{{Annotated link|Galaxy formation and evolution}}
*{{Annotated link|Natural units}}
*{{Annotated link|Virgo Consortium}}
*{{Annotated link|Barnes–Hut simulation}}
*{{Annotated link|Bolshoi cosmological simulation}}

==References==
{{reflist}}

===Further reading===
* {{cite journal | first = Sebastian | last = von Hoerner | author-link = Sebastian von Hoerner | language = de | title = Die numerische Integration des n-Körper-Problemes für Sternhaufen. I | journal = Zeitschrift für Astrophysik | volume = 50 | pages = 184 | bibcode = 1960ZA.....50..184V | date = 1960 }}
* {{cite journal | first = Sebastian | last = von Hoerner | author-link = Sebastian von Hoerner | language = de | title = Die numerische Integration des ''n''-Körper-Problemes für Sternhaufen. II | journal = Zeitschrift für Astrophysik | volume = 57 | pages = 47 | bibcode = 1963ZA.....57...47V | date = 1963 }}
* {{cite book|title=Gravitational ''N''-body Simulations: Tools and Algorithms| first = Sverre J. | last = Aarseth | author-link = Sverre Aarseth | date = 2003 | publisher = [[Cambridge University Press]] | isbn = 978-0-521-12153-8 }}
* {{cite journal | first = Edmund | last = Bertschinger | title = Simulations of structure formation in the universe | journal = Annual Review of Astronomy and Astrophysics | volume = 36 | issue = 1 | pages = 599–654 | date = 1998 | doi = 10.1146/annurev.astro.36.1.599 | bibcode = 1998ARA&A..36..599B }}
* {{cite book|title=Galactic Dynamics| first1 = James | last1 = Binney | first2 = Scott | last2 = Tremaine | date = 1987 | publisher = [[Princeton University Press]] | isbn = 978-0-691-08445-9 }}
*{{Cite book
 | last1 = Callahan | first1 = Paul B.
 | last2 = Kosaraju | first2 = Sambasiva Rao 
 | contribution = A decomposition of multidimensional point sets with applications to ''k''-nearest-neighbors and n-body potential fields (preliminary version)
 | publisher = ACM
 | title = STOC '92: Proc. ACM Symp. Theory of Computing
 | date = 1992}}.

{{Numerical integrators}}

[[Category:Physical cosmology]]
[[Category:Gravity]]
[[Category:Simulation]]
[[Category:Cosmological simulation]]
[[Category:Articles containing video clips]]
[[Category:Computational physics]]
[[Category:Particles]]