Barnes-Hut: Simulating 100,000 Gravitating Bodies in Real Time

The n² wall

Gravity is a pairwise force. To advance an N-body system by one step you need the acceleration of every body, and the direct sum gives it exactly: for each body, add the pull of every other body. That is n(n−1)/2 pairs — O(n²) work per step. A thousand bodies means half a million interactions per step and runs comfortably at 60 fps in a browser. A hundred thousand bodies means five billion interactions per step, and you are now waiting minutes per frame.

The Barnes–Hut algorithm (Josh Barnes and Piet Hut, published in Nature in 1986) breaks the wall with one physical observation: a distant cluster of stars pulls almost exactly like a single star at its centre of mass. You do not need to know that the Andromeda galaxy has a trillion stars in it to compute its pull on the Sun — you need its total mass and its centre of mass. Barnes–Hut turns that observation into a data structure and gets O(n log n) per step.

The tree

Recursively subdivide the simulation volume. In 2D each square splits into four quadrants — a quadtree; in 3D each cube splits into eight octants — an octree. You keep subdividing a cell until it contains at most one body. Every internal node caches two aggregates computed on the way back up: the total mass of the bodies below it, and their centre of mass.

insert(node, body):
  if node is empty leaf:            node.body = body; return
  if node is leaf with one body:    subdivide(node); reinsert the old body
  node.mass += body.mass            // aggregates updated on the way down
  node.com  = weighted average of centres of mass
  insert(child_containing(body.pos), body)

Building the tree costs O(n log n) for a reasonably uniform distribution — each of the n insertions descends a tree of depth ~log n. The tree is rebuilt from scratch every step; that is simpler than updating it and, in practice, the build is a small fraction of the total cost compared with the force traversal.

The θ criterion

To get the force on one body, walk the tree from the root. At each node compare the node's side length s with the distance d from the body to the node's centre of mass. If the node is small enough and far enough away — that is, if

s / d < θ           // θ = opening angle, typically 0.5

— treat the whole node as a single point mass at its centre of mass, add that one contribution and do not descend further. Otherwise open the node and recurse into its children. A leaf is always evaluated directly, and a body never attracts itself.

force(body, node):
  if node is a leaf:
    if node.body !== body: return pointForce(body, node.body);
    return 0;
  const s = node.size, d = distance(body.pos, node.com);
  if (s / d < theta) return pointForce(body, node);   // one term, whole subtree
  let f = 0;
  for (const c of node.children) if (c) f += force(body, c);
  return f;

Because each body only descends a logarithmic number of levels before the criterion is satisfied, the traversal costs O(log n) per body and O(n log n) per step overall. That is the entire algorithm.

The trade-off you are choosing with θ

θ is a single knob between exactness and speed. θ = 0 never satisfies the criterion, so every node is opened, every pair is evaluated and you are back to the exact O(n²) direct sum. Increasing θ opens fewer nodes: fewer interaction terms, faster steps, larger approximation error in the force. Values around 0.5 are the conventional default in astrophysics, with 0.3 for accuracy-critical work and 1.0 for a fast visual demo where a slightly wrong tidal field is unnoticeable.

Two caveats worth knowing. First, the monopole (centre-of-mass) approximation ignores the quadrupole moment of the cell, so error grows with θ faster than a naive estimate suggests; some codes add the quadrupole term to buy accuracy back cheaply. Second, Barnes–Hut error is not random noise — it is correlated with the tree geometry, so a body sitting near a cell boundary can see its force jitter as the tree is rebuilt. Softening the potential (replacing 1/r² with 1/(r² + ε²)) is standard practice anyway: it stops the force from exploding when two bodies pass close, which would otherwise demand an impossibly small time step.

How the N-body simulation here uses it

The simulation on this site builds a fresh quadtree each frame, walks it once per body with θ = 0.5 and a softening length of a few pixels, and then advances the state with Leapfrog — a symplectic integrator, because a gravitational system integrated with RK4 or explicit Euler bleeds or gains energy and the galaxy visibly collapses or evaporates. The force approximation and the integrator are chosen together: there is no point spending four force evaluations per step on RK4 when each force is itself an approximation controlled by θ.

each frame:
  1. build quadtree from current positions   O(n log n)
  2. accumulate mass + centre of mass upward O(n)
  3. traverse per body with θ = 0.5          O(n log n)
  4. Leapfrog kick–drift–kick                O(n)

For completeness: Barnes–Hut is not the only escape from O(n²). The Fast Multipole Method expands the field in multipoles on both sides of the interaction and reaches O(n) for a prescribed accuracy, at the price of a considerably heavier implementation; particle-mesh methods solve for the potential on a grid with an FFT and are the usual choice for cosmological volumes. For an interactive canvas with tens of thousands of bodies, a quadtree and one well-chosen θ is the sweet spot.

Frequently asked questions

Try it live

Everything above runs in your browser — open N-Body Gravity and change the parameters while it is running. Nothing is installed, nothing is uploaded, the whole model lives in one tab.