N-Body Gravity Simulation: How Planets Dance

Newton's Law of Universal Gravitation

Every mass in the universe attracts every other mass with a force given by Newton's law:

F = G × m₁ × m₂ / r²

Where G is the gravitational constant (6.674 × 10⁻¹¹ N m² kg⁻²), m₁ and m₂ are the masses of the two bodies, and r is the distance between their centres. The force is attractive and acts along the line connecting the two bodies.

For two bodies, this gives a perfectly solvable system. The bodies orbit their common centre of mass in ellipses — Kepler's laws fall out directly. Add a third body and everything changes.

Why N-Body Is Hard: The Three-Body Problem

In 1887, Henri Poincaré proved that the three-body gravitational problem has no general closed-form solution. You cannot write down an equation that gives the position of three gravitating masses at time t given their initial conditions. The system is generally chaotic — small differences in initial conditions lead to completely different long-term trajectories.

This is not a failure of mathematics. It is a fundamental property of the equations. Three bodies gravitationally interacting can exhibit every mode of behaviour: stable orbits, quasi-periodic motion, resonances, chaotic wandering, and ejection of one body to infinity.

The only practical way to predict N-body evolution is numerical integration: take small time steps, compute forces, update velocities and positions, repeat. This is what the simulation does — and what every astrophysics code from solar system models to galaxy simulations does.

The Naive Algorithm: O(n²) Complexity

The straightforward approach to computing gravitational forces is to check every pair of bodies:

for each body i:
    force[i] = (0, 0, 0)
    for each body j ≠ i:
        r = position[j] - position[i]
        dist = length(r)
        force[i] += G * mass[i] * mass[j] / dist² * normalize(r)

This requires n × (n−1) force evaluations per time step — O(n²) complexity. For 100 bodies: 9,900 evaluations. For 1,000 bodies: 999,000 evaluations. For 10,000 bodies: nearly 100 million evaluations. For a real galaxy with 100 billion stars, the naive algorithm is completely infeasible.

Softening to Prevent Singularities

There is a numerical problem lurking in the force equation: when two bodies get very close, r approaches zero and the force diverges to infinity. In reality, stars are not point masses — close encounters involve complex physics. In simulation, we add a softening length ε:

F = G × m₁ × m₂ / (r² + ε²)

This caps the maximum force and prevents numerical explosions during close encounters. The softening length is a simulation parameter — too large and the physics becomes unrealistic; too small and close encounters destabilise the integrator.

Barnes-Hut: Reducing to O(n log n)

In 1986, Josh Barnes and Piet Hut published an algorithm that reduces N-body complexity from O(n²) to O(n log n) — a transformative improvement. The key insight: a distant cluster of many bodies can be approximated as a single body at the cluster's centre of mass. The gravitational error from this approximation is small when the cluster is far away.

The Barnes-Hut algorithm builds an octree (a tree that recursively divides 3D space into eight octants at each node). Each internal node stores the total mass and centre of mass of all bodies in its subtree.

When computing the force on body i, the tree is traversed. At each node, a decision is made based on the ratio s/d — where s is the node's spatial size and d is the distance from body i to the node's centre of mass. If s/d < θ (the opening angle threshold, typically 0.5 to 1.0), the node is treated as a single point mass. Otherwise, the node's children are opened and evaluated individually.

function compute_force(body i, node):
    if node is leaf:
        add direct force from node's body
    else:
        s = node size
        d = distance from i to node's centre of mass
        if s/d < θ:
            add approximate force from node's mass
        else:
            for each child of node:
                compute_force(body i, child)

With θ ≈ 0.5, the error in force estimates is below 1% while computation drops from hours to seconds for the same particle count. Modern cosmological simulations use variants of Barnes-Hut to simulate billions of dark matter particles.

Leapfrog Integration: Stability Over Accuracy

Knowing the forces is only half the problem. You also need to integrate the equations of motion — update positions and velocities — in a way that remains stable over many thousands of time steps.

Simple Euler integration (x += v*dt; v += a*dt) seems natural but has a fatal flaw for orbital mechanics: it does not conserve energy. The total energy of the system drifts upward over time, and orbits slowly spiral outward. Over long simulations, planets escape their stars.

The leapfrog integrator solves this by staggering position and velocity updates by half a time step:

// kick-drift-kick (KDK) form:
v_half = v + a * (dt/2)      // half-step velocity
x = x + v_half * dt           // full-step position
a = compute_forces(x)         // forces at new position
v = v_half + a * (dt/2)       // second half-step velocity

The leapfrog is symplectic — it exactly conserves a modified Hamiltonian close to the true one. In practice this means energy oscillates around the correct value rather than drifting, allowing simulations to run for billions of time steps without accumulating error.

Higher-order integrators (Runge-Kutta 4, Yoshida 6th order) are more accurate per step but require more force evaluations. For gravitational N-body, leapfrog with a small time step is usually the best trade-off.

Real Astrophysics Applications

N-body simulation is a core tool of modern astrophysics:

• Planetary system formation — simulating the solar system's early history, when hundreds of planetesimals collided and merged over millions of years to form the planets.
• Galaxy collisions — the Milky Way will collide with the Andromeda galaxy in about 4.5 billion years. N-body simulations from the 1970s onward have predicted the merger's timeline and structure.
• Globular clusters — dense spherical collections of 100,000 to 1 million stars, simulated to understand stellar dynamics, black hole formation, and escape rates.
• Large-scale structure — cosmological simulations like Millennium and IllustrisTNG track billions of dark matter particles from the Big Bang to the present, reproducing the cosmic web of filaments, voids, and galaxy clusters.
• Asteroid dynamics — near-Earth asteroid trajectories are computed via N-body integration including perturbations from all eight planets, the Moon, and Pluto.

The Stability of Our Solar System

Is our own solar system stable? This turns out to be a surprisingly difficult question. Numerical integrations run for billions of simulated years show that the inner planets (Mercury through Mars) have a small but nonzero probability of entering chaotic orbits. Mercury's orbit could become unstable on timescales of a few billion years, potentially colliding with Venus or being ejected.

The outer planets are more stable but not perfectly so. The four-body resonance between Jupiter, Saturn, Uranus, and Neptune means small perturbations can accumulate. The solar system as we know it may be one of many possible configurations that happened to persist long enough for complex life to evolve.