The N-body problem is the age-old challenge in physics of predicting the motion of N individual bodies interacting with each other under a force law. This has been motivated historically by understanding the orbits of astronomical bodies on the sky, such as the various planets of the solar system. It remains of central importance today with astrophysical applications such the dynamics of stellar clusters and structure formation in the Universe through dark matter simulations such as the Millennium Run (pictured below.)
Unfortunately the problem can only be solved through first integrals in a very small number of restricted cases. Motivated by trying to understand the orbit of the Earth around the Sun, the two-body problem was solved as early as the 17th century. The restricted three-body problem, where the third body is less massive and moves in the same plane as the two first bodies, can be solved completely as well. That is however the extent of simple solutions to the N-body problem in classical mechanics. A general analytical solution for arbitrary N exists through convergent power series; in practice however the convergence is too slow for the solution to be implemented efficiently.
Luckily, computational physics and applied mathematics come to the rescue. N-body problems can be solved by numerically integrating the equations of motion. In a collection of N self-gravitating particles, every particle interacts with every other particle. As a result we need to calculate (N-1) force terms for each particle, leading to a total number of N(N-1) calculations to perform the global force budget of the system. This N2 scaling of the number of computations for such direct summation, or Particle-Particle (PP) methods, doesn’t bode well for numerical computation. One alternative is Particle-Mesh (PM) methods, where forces are interpolated on a mesh. Hybrid PP-PM methods have been used to combine greater speed and accurate treatment of high-density regions in the mesh.
Tree codes were developed in the 1980s as an alternative to PM methods. Tree algorithms divide the total particle volume into a hierarchical tree. Going down the tree when computing the force on a given particle, nearby bodies can be treated individually, while distant objects are grouped into massive pseudo-particles and their effect is approximated by their multipole moments. Such tree algorithms allow for computational costs of order N log N and are intrinsically much more adaptive to high clustering than PM methods.
One of the most common tree algorithms, the Barnes-Hut tree, is the one we will investigate here.
Applications of the N-body problem
Solving the N-body problem is at the heart of cosmological simulations such as the Millennium Run (pictured above), where millions of collisionless dark matter particles interact through gravity and form the structures that we see today in the Universe: dark matter halos hosting spiral galaxies, galaxy clusters, dwarf satellites, etc. An example tree code for such simulations is GADGET.
Another application of the N-body problem is the modeling of star clusters. Globular clusters are spherical collections of hundreds of thousands of old stars tightly bound by gravity. The formation mechanism and history of these objects is an active area of research. N-body simulations are widely used to shed light on their dynamical evolution.
For our project we simulate a globular cluster’s evolution over a time range corresponding to at least 10 dynamical times (the characteristic time scale for collapse or other changes to propagate through the system.) To create a realistic set of initial conditions, we initialize the distribution to follow a spherically symmetric Plummer density profile (H. C. Plummer 1911, Monthly Notices of the Royal Astronomical Society, 460, 71.) This is done by sampling radii from the cumulative mass distribution, and velocity magnitudes from the velocity distribution function through rejection sampling. A great reference can be found at the Art of Computational Science by Piet Hut and Jan Makino.