Introduction

Imagine trying to simulate a galaxy. Every star exerts a gravitational pull on every other star. With nn stars, that is n×(n−1)n \times (n-1) pairs — roughly n2n^2 force calculations per time step. At one million stars you need a trillion calculations every frame. At a billion stars the number becomes astronomical in both senses of the word.

This is the N-body problem: given nn point masses, compute the net gravitational force on each one so you can advance their positions by a tiny time step and repeat. It underlies everything from planetary mechanics to cosmological simulations.

The naive O(n2)O(n^2) algorithm is exact but hopeless at scale. In 1986 Josh Barnes and Piet Hut published a clever approximation: instead of summing every pairwise interaction, group distant particles together and treat the group as a single mass at its center of mass. Their algorithm runs in O(nlog⁡n)O(n \log n) — fast enough to simulate millions of bodies on ordinary hardware.

The key data structure is a quadtree (in 2-D) or octree (in 3-D): a recursive spatial subdivision that lets you decide, for each particle, which distant clusters are far enough away to be summarized safely.

Try It

The simulation below spawns particles in a rotating disk and evolves them under gravity. Toggle Show quadtree to see how the algorithm partitions space — each cell summarizes all its particles as a single mass when a particle is far enough away.

<div class="controls">
  <label>θ (theta): <input id="theta" type="range" min="0.1" max="1.5" step="0.05" value="0.5"> <span id="theta-val">0.50</span></label>
  <label><input id="show-tree" type="checkbox"> Show quadtree</label>
  <button id="btn-pause" type="button">Pause</button>
  <button id="btn-reset" type="button" class="ghost">Reset</button>
</div>
<canvas id="canvas" width="480" height="320"></canvas>
<div class="stats" id="stats"></div>

* { box-sizing: border-box; }
body { margin: 0; background: #0b0d14; font-family: system-ui, sans-serif; color: #c8d4e0; }
.controls { display: flex; align-items: center; gap: .8rem; flex-wrap: wrap; padding: .5rem .6rem; background: #131622; border-bottom: 1px solid #1e2538; font-size: .85rem; }
label { display: flex; align-items: center; gap: .3rem; cursor: pointer; user-select: none; }
input[type=range] { width: 90px; accent-color: #5b9bd5; }
button { font: 600 13px system-ui; padding: .3rem .8rem; border-radius: 6px; cursor: pointer; border: 1px solid #5b9bd5; background: #5b9bd5; color: #fff; }
button.ghost { background: transparent; color: #5b9bd5; }
canvas { display: block; width: 100%; max-width: 480px; }
.stats { font-size: .78rem; padding: .35rem .6rem; background: #131622; color: #7a91a8; min-height: 1.6em; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Barnes-Hut N-body simulation. Build a quadtree each step, then walk it
const N = 300;
const G = 0.4;
const DT = 0.4;
const SOFT = 1.5;

const canvas = document.getElementById('canvas');
const ctx = canvas.getContext('2d');
const W = canvas.width, H = canvas.height;

// with the opening-angle criterion (s/d < theta) to approximate distant
let theta = 0.5;
let showTree = false;
let paused = false;

const thetaSlider = document.getElementById('theta');
const thetaVal = document.getElementById('theta-val');
const showTreeCb = document.getElementById('show-tree');
const btnPause = document.getElementById('btn-pause');
const btnReset = document.getElementById('btn-reset');

thetaSlider.addEventListener('input', () => {
  theta = parseFloat(thetaSlider.value);
  thetaVal.textContent = theta.toFixed(2);
});
showTreeCb.addEventListener('change', () => { showTree = showTreeCb.checked; });
btnPause.addEventListener('click', () => {
  paused = !paused;
  btnPause.textContent = paused ? 'Resume' : 'Pause';
});
btnReset.addEventListener('click', init);

// clusters as single masses. O(n log n) vs O(n^2) naive.
let particles = [];

function init() {
  particles = [];
  const cx = W / 2, cy = H / 2, r = Math.min(W, H) * 0.38;
  for (let i = 0; i < N; i++) {
    const angle = Math.random() * 2 * Math.PI;
    const dist = r * Math.sqrt(Math.random()) * 0.9 + r * 0.05;
    const x = cx + dist * Math.cos(angle);
    const y = cy + dist * Math.sin(angle);
    // tangential velocity for a circular orbit (Kepler approximation)
    const v = Math.sqrt(G * N * 0.5 * dist / r) * 0.6;
    const vx = -Math.sin(angle) * v + (Math.random() - 0.5) * 0.3;
    const vy = Math.cos(angle) * v + (Math.random() - 0.5) * 0.3;
    particles.push({ x, y, vx, vy, ax: 0, ay: 0 });
  }
  if (paused) { paused = false; btnPause.textContent = 'Pause'; }
}

// Quadtree node: bounding box + aggregated mass + center of mass
function makeNode(x, y, w, h) {
  return { x, y, w, h, mass: 0, cx: 0, cy: 0, children: null, particle: null };
}

// Insert a particle into the quadtree, splitting leaves as needed
function insert(node, p) {
  if (node.particle === null && node.children === null) {
    node.particle = p;
    node.mass = 1; node.cx = p.x; node.cy = p.y;
    return;
  }
  if (node.children === null) {
    const hw = node.w / 2, hh = node.h / 2;
    node.children = [
      makeNode(node.x, node.y, hw, hh),
      makeNode(node.x + hw, node.y, hw, hh),
      makeNode(node.x, node.y + hh, hw, hh),
      makeNode(node.x + hw, node.y + hh, hw, hh)
    ];
    const old = node.particle; node.particle = null;
    insertIntoChild(node, old);
  }
  insertIntoChild(node, p);
  node.mass++;
  node.cx = (node.cx * (node.mass - 1) + p.x) / node.mass;
  node.cy = (node.cy * (node.mass - 1) + p.y) / node.mass;
}

function insertIntoChild(node, p) {
  const mx = node.x + node.w / 2, my = node.y + node.h / 2;
  const qi = (p.x >= mx ? 1 : 0) + (p.y >= my ? 2 : 0);
  insert(node.children[qi], p);
}

// Recursively accumulate force on particle p from quadtree node
function addForce(node, p, out) {
  if (node.mass === 0) return;
  if (node.particle === p) return;
  const dx = node.cx - p.x, dy = node.cy - p.y;
  const dist2 = dx * dx + dy * dy + SOFT * SOFT;
  // opening-angle criterion: treat this node as a point mass
  if (node.children === null || node.w / Math.sqrt(dist2) < theta) {
    const inv = G * node.mass / (dist2 * Math.sqrt(dist2));
    out.ax += dx * inv;
    out.ay += dy * inv;
    return;
  }
  for (const child of node.children) addForce(child, p, out);
}

// Draw quadtree cell boundaries for visualization
function drawTree(node) {
  if (node.mass === 0) return;
  ctx.strokeRect(node.x, node.y, node.w, node.h);
  if (node.children) node.children.forEach(drawTree);
}

let interactions = 0;

function step() {
  // Rebuild the quadtree from scratch each time step
  const root = makeNode(0, 0, W, H);
  for (const p of particles) insert(root, p);

  // Compute accelerations via Barnes-Hut tree walk
  interactions = 0;
  for (const p of particles) {
    const out = { ax: 0, ay: 0 };
    addForce(root, p, out);
    p.ax = out.ax; p.ay = out.ay;
    interactions++;
  }

  // Leapfrog-style velocity-Verlet integration + wrap-around boundary
  for (const p of particles) {
    p.vx += p.ax * DT; p.vy += p.ay * DT;
    p.x += p.vx * DT; p.y += p.vy * DT;
    if (p.x < 0) p.x += W; if (p.x > W) p.x -= W;
    if (p.y < 0) p.y += H; if (p.y > H) p.y -= H;
  }
}

function draw() {
  ctx.fillStyle = 'rgba(11,13,20,0.35)';
  ctx.fillRect(0, 0, W, H);

  if (showTree) {
    const root2 = makeNode(0, 0, W, H);
    for (const p of particles) insert(root2, p);
    ctx.strokeStyle = 'rgba(91,155,213,0.25)';
    ctx.lineWidth = 0.5;
    drawTree(root2);
  }

  ctx.fillStyle = '#c8d4e0';
  for (const p of particles) {
    ctx.beginPath();
    ctx.arc(p.x, p.y, 1.2, 0, 2 * Math.PI);
    ctx.fill();
  }
}

function loop() {
  if (!paused) { step(); }
  draw();
  document.getElementById('stats').textContent =
    i18nFmt('%n% particles · θ = %theta% · Barnes-Hut O(n log n)', { n: N, theta: theta.toFixed(2) });
  requestAnimationFrame(loop);
}

init();
loop();
Result
Console

Notice how the quadtree cells near a particle are small and numerous (high resolution), while distant regions are large single boxes (one interaction instead of hundreds). The theta parameter θ\theta controls the trade-off: smaller values are more accurate but slower; larger values are faster but rougher. The default θ=0.5\theta = 0.5 is the value Barnes and Hut recommended in their original paper.

The Real Complexity

Naive N-body is O(n2)O(n^2) per time step: every particle loops over every other particle. Doubling nn quadruples the work.

Barnes-Hut breaks that with two observations:

  1. Build a quadtree. Recursively subdivide the simulation region into four quadrants until each leaf holds at most one particle. This costs O(nlog⁡n)O(n \log n) in the average case.
  2. Apply the opening-angle criterion. For each particle pp, walk the tree. At each internal node (a cell of width ss at distance dd from pp), check whether s/d<θs / d < \theta. If yes, the whole cell's mass can stand in for all its particles — one interaction instead of potentially thousands. If no, recurse into the children.

The fraction of nodes that survive the opening-angle test is O(log⁡n)O(\log n) on average, so the total cost per particle is O(log⁡n)O(\log n), giving the full step O(nlog⁡n)O(n \log n).

The approximation error in the force on particle pp is bounded by a term proportional to θ2⋅(s/d)2\theta^2 \cdot (s/d)^2, so halving θ\theta roughly quarters the error while increasing runtime. In practice θ∈[0.5,1.0]\theta \in [0.5, 1.0] gives results close enough for astrophysical simulations.

Compare this to the fast Fourier transform: both algorithms replace a seemingly irreducible O(n2)O(n^2) task with O(nlog⁡n)O(n \log n) by exploiting structure in the problem — spatial proximity here, frequency decomposition there.

Where It Matters

Anywhere you sum a long-range interaction over many particles, Barnes-Hut or one of its cousins shows up:

  • Astrophysical simulations: galaxy formation, dark-matter halos, stellar clusters — Barnes-Hut (and its successor the fast multipole method) made large-scale cosmological simulations feasible on real machines.
  • Molecular dynamics: electrostatic forces between atoms follow the same 1/r21/r^2 law as gravity. The particle-mesh Ewald and fast multipole methods are direct descendants of the same idea.
  • Force-directed graph layouts: tools like D3.js use Barnes-Hut to repel graph nodes, turning a layout that would freeze on large graphs into one that runs in the browser.
  • Smoothed-particle hydrodynamics (SPH): fluid simulation where each "particle" carries a small blob of fluid; neighbor queries use the same spatial tree.
  • Machine learning: tt-SNE uses Barnes-Hut to approximate gradient repulsions among nn data points, enabling visualization of datasets with tens of thousands of items.

The common thread is long-range interaction + spatial locality: nearby things need full precision, distant things can be grouped. Whenever those two conditions hold, a spatial tree can replace O(n2)O(n^2) with O(nlog⁡n)O(n \log n).

To see a related O(nlog⁡n)O(n \log n) idea in a completely different domain, compare with the fast Fourier transform.

Conclusion

Barnes and Hut's insight was not to compute less — it was to compute smartly. A distant cluster of a thousand stars can be replaced by one number (its total mass) and one point (its center of mass) without meaningfully affecting the trajectory of a star far away. The quadtree automates the decision of what counts as "far away."

The result is an algorithm that scales to millions of particles on a laptop, to billions on a supercomputer — all because of a single ratio check: s/d<θs / d < \theta.

That ratio check is an example of a broader pattern in algorithm design: controlled approximation. Rather than demanding exact answers at every step, accept a bounded error and win orders of magnitude in speed. The same philosophy drives fast Fourier transform algorithms, compressed sensing, and the approximate nearest-neighbor searches that power modern recommendation systems.

Gravity is, in the end, a tractable problem — not because the physics simplifies, but because the geometry does.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/n-body-barnes-hut/Content licensed under CC BY-NC 4.0.