N-Body Gravity: Simulating the Cosmos in Your Browser

The N-body problem is one of the oldest in physics: given N point masses interacting through gravity, predict their future positions. For N = 2, Newton solved it exactly in the seventeenth century — the two-body problem has a closed-form solution in terms of ellipses, parabolas, and hyperbolas. For N = 3, Poincaré proved in 1889 that no general closed-form solution exists. For N = 10,000 — the scale of our galaxy simulation — you need an efficient algorithm and careful numerical tricks.

The Brute Force Approach: O(N²)

The straightforward way to simulate N gravitating bodies is to compute, for every body, the gravitational force from every other body. Newton's law of gravitation gives the force between bodies i and j as:

F_ij = G m_i m_j / |r_i − r_j|² in the direction from i toward j.

For each of N bodies, you sum N − 1 pairwise forces. That's O(N²) force calculations per timestep. For N = 1,000, that's one million calculations per frame. For N = 10,000, it's one hundred million. At 60 frames per second, a brute-force 10,000-body simulation requires six billion force calculations per second — too slow for real-time in a browser.

Brute force is nonetheless useful for small N — fewer than about 500 bodies — where it is maximally accurate and simple to implement. Our simulation uses brute force for small particle counts and switches to Barnes-Hut above a threshold.

Barnes-Hut: O(N log N) with a Quadtree

The Barnes-Hut algorithm, published by Josh Barnes and Piet Hut in 1986, reduces the complexity to O(N log N) using a key approximation: distant groups of bodies can be treated as a single body at their centre of mass.

The algorithm builds a quadtree (in 2D) or octree (in 3D) that recursively subdivides space into cells. Each cell stores the total mass and centre of mass of all bodies within it. To compute the force on body i, you traverse the tree. For each node, if the node is far enough away relative to its size — specifically, if the node's width s divided by its distance d satisfies s/d < θ (a tuning parameter, typically 0.5–1.0) — you treat the entire node as a single body at its centre of mass. Otherwise, you recurse into the node's children.

Algorithm	Complexity	N = 1,000	N = 10,000
Brute Force	O(N²)	1M ops	100M ops
Barnes-Hut	O(N log N)	~10K ops	~133K ops
Fast Multipole	O(N)	~1K ops	~10K ops

The θ parameter controls the accuracy-speed trade-off. θ = 0 gives exact brute-force results. θ = 1.0 gives significant approximation errors. θ = 0.5 provides a good balance — force errors of less than 1% for most bodies while keeping the simulation interactive.

Gravitational Softening: Preventing Singularities

Newton's gravity law has a singularity: as two bodies get very close, the force between them grows without bound (1/r² → infinity as r → 0). In a real galaxy, actual close encounters between stars are extraordinarily rare because stars are tiny relative to the distances between them. In a simulation, particles can get arbitrarily close, causing numerical explosions — one timestep sees enormous force, produces enormous velocity, and the particle flies off to infinity.

The standard fix is gravitational softening: replace the denominator r² with r² + ε², where ε is the softening length. This caps the maximum force at a finite value when two particles are within distance ε of each other. The force now looks like:

F = G m_i m_j / (r² + ε²)

Physically, this is equivalent to smearing each point mass into a sphere of radius ε. For values of r much larger than ε, the force is indistinguishable from point masses. For very small r, the force is capped at G m_i m_j / ε².

Choosing ε is an art. Too large and you suppress real dynamical effects — particles that should slingshot each other instead pass through smoothly. Too small and close encounters dominate the timestep, slowing the simulation drastically. A typical choice for a galaxy simulation is ε equal to about 1% of the mean inter-particle spacing.

Time Integration: The Leapfrog Method

Given the forces, how do we advance particle positions and velocities through time? Simple Euler integration — update velocity, then position, using current force — is fast but loses energy over time, causing orbits to slowly spiral outward. For an N-body simulation run for thousands of timesteps, this energy drift ruins the results.

The leapfrog integrator, also called the Störmer-Verlet method, solves this with a staggered update. Velocities are evaluated at half-integer timesteps while positions are evaluated at integer timesteps:

• v_n+1/2 = v_n-1/2 + a_n Δt
• x_n+1 = x_n + v_n+1/2 Δt

This tiny change — velocities and positions no longer updated at the same time — makes the integrator symplectic: it exactly conserves a slightly modified Hamiltonian (energy function) of the system. In practice, this means total energy oscillates around the true value rather than drifting — orbits stay stable for millions of timesteps without becoming unphysical.

For even better accuracy, higher-order symplectic integrators like Yoshida's fourth-order method or the classic fourth-order Runge-Kutta (non-symplectic but highly accurate per step) are used in high-precision astronomical simulations. Our browser simulation uses leapfrog for its simplicity, stability, and excellent energy conservation.

Galaxy Formation and Emergent Structure

Run enough gravitating particles with appropriate initial conditions and remarkable structure emerges spontaneously. Give a cloud of particles some initial angular momentum and let gravity collapse it: the cloud flattens into a disc, clumps form under Jeans instability, and spiral arms emerge from differential rotation. These are not programmed features — they arise from N simple particles following Newton's law.

Real cosmological simulations like the Illustris project and EAGLE use the same fundamental algorithms — Barnes-Hut or Fast Multipole for gravity, SPH or moving mesh for gas dynamics — scaled to billions of particles on supercomputer clusters. The qualitative behaviours they produce match observed galaxy morphologies remarkably well. Our browser simulation is a miniature version of the same computational physics that professional cosmologists use to understand the universe's large-scale structure.