Three-Body Problem

The setup

Take three point masses m₁, m₂, m₃ at positions r₁, r₂, r₃. Each pulls on the other two through Newton's law of universal gravitation, an inverse-square force directed along the line joining them. The acceleration of body i is the vector sum of the pulls from the other two:

d²rᵢ/dt²  =  Σ   G·mⱼ·(rⱼ − rᵢ) / |rⱼ − rᵢ|³
           j≠i

That is three coupled second-order vector ODEs — nine scalar second-order equations, or eighteen first-order equations once you split each into position and velocity. It looks innocent. The two-body version of exactly this equation is the cleanest exact solution in all of physics. Add one body and it becomes a problem that has resisted full analytic solution for over three centuries.

Why two bodies is easy and three is not

The two-body problem has a trick: switch to the relative coordinate r = r₂ − r₁ and the center of mass. The center of mass drifts at constant velocity (it carries no dynamics), and the relative motion reduces to a single body of reduced mass μ = m₁m₂/(m₁+m₂) orbiting a fixed center under an inverse-square force. That central-force problem conserves energy and angular momentum, which is exactly enough to integrate it down to a conic section — an ellipse, parabola, or hyperbola. Kepler's laws fall straight out.

With three bodies the trick fails. You can still remove the center-of-mass motion, but the leftover relative problem has no extra symmetry to exploit. Count the bookkeeping: the full system has 18 first-order variables. The classical conserved quantities are total energy (1), total linear momentum (3), total angular momentum (3), and the uniform motion of the center of mass (3), plus time-translation lets you drop one more — ten integrals in all. To solve a system by quadrature you need as many independent integrals as equations. Ten is not enough for eighteen, and there is no eleventh.

Feature	Two-body problem	Three-body problem
First-order equations	12 (6 after COM)	18
Independent integrals	Enough to integrate	Only 10 — too few
General solution	Closed form (conic section)	None known / proven impossible in elementary terms
Long-term behaviour	Periodic, fully predictable	Generically chaotic
Sensitivity to ICs	None (regular orbits)	Exponential (positive Lyapunov exponent)

The proof there is no closed form

This is not a failure of cleverness — it is a theorem. In 1887 Heinrich Bruns proved that the only algebraic integrals of the three-body problem (functions built from the coordinates and momenta) are the ten classical ones. Three years later Henri Poincaré, working on a prize problem set for King Oscar II of Sweden, extended the result: there are no new single-valued analytic integrals either. The system is, in the technical sense, non-integrable.

In 1912 Karl Sundman found a genuine convergent power series for the three-body problem (excluding the measure-zero set of triple collisions). It is a real solution — but the series converges so slowly that getting modest accuracy would need on the order of 10⁸⁰⁰⁰⁰⁰⁰ terms. It is a triumph of analysis and completely useless for computation. The practical conclusion stands: to know where three gravitating bodies go, you integrate numerically.

Sensitive dependence: the birth of chaos

Poincaré's deeper discovery was geometric. Studying the trajectories near the periodic solutions, he found that stable and unstable manifolds intersect in a "homoclinic tangle" — an infinitely folded web that makes nearby orbits diverge wildly. This was the first clear description of deterministic chaos, decades before Lorenz rediscovered it in weather models.

Quantitatively, two trajectories starting a distance δ₀ apart separate as

δ(t) ≈ δ₀ · e^(λ t)

where λ is the largest Lyapunov exponent. The Lyapunov time τ = 1/λ sets the prediction horizon. After about 23 Lyapunov times (the natural log of 10¹⁰), an initial error in the tenth decimal place has grown to order one and the prediction is worthless. Because no measurement is infinitely precise, long-term forecasting of a chaotic three-body system is fundamentally impossible — not hard, impossible.

The exact solutions that do exist

"No general solution" does not mean "no solutions." Special, highly symmetric families are known exactly:

• Euler collinear (1767). Three bodies stay on a common straight line that rotates rigidly; the ratios of their separations are fixed by the masses. Three such configurations exist; they are the L1, L2, L3 Lagrange points of the restricted problem.
• Lagrange equilateral (1772). Three bodies sit at the vertices of an equilateral triangle that rotates and breathes while staying equilateral. The stable versions are the L4 and L5 points — where the Sun–Jupiter system parks the Trojan asteroids.
• Figure-eight choreography (Moore 1993; Chenciner & Montgomery 2000). Three equal masses chase one another around a single figure-eight curve, each lagging the next by one-third of a period. Remarkably it appears stable. It launched a whole industry of finding "choreographies."

Solution	Discovered	Geometry	Where it shows up
Euler collinear	1767	Rotating straight line	L1 / L2 / L3
Lagrange equilateral	1772	Rotating equilateral triangle	L4 / L5, Trojan asteroids
Figure-eight	1993 / 2000	Single figure-eight curve	Equal-mass choreography
Broucke–Hénon family	1970s	Periodic relative orbits	Numerical catalogs

The restricted three-body problem

The most useful approximation lets one body be massless. A spacecraft or asteroid feels the gravity of two heavy bodies (say the Sun and Earth) but does not perturb them, so the heavies orbit each other on fixed circular Keplerian paths. Work in the frame that rotates with them and the problem gains a conserved quantity — the Jacobi integral C — and five equilibrium points appear. Those are the Lagrange points L1–L5.

This model is mission-design bread and butter. The James Webb Space Telescope orbits the Sun–Earth L2 point, about 1.5 million km from Earth; SOHO watches the Sun from L1; spacecraft ride the "interplanetary transport network" of low-energy manifolds threading through these points. The effective potential in the rotating frame, including the centrifugal term, is

Ω(x,y) = ½ ω²(x² + y²) + G·m₁/r₁ + G·m₂/r₂

and the Lagrange points are its stationary points; L4 and L5 are local maxima that are nonetheless dynamically stable when the mass ratio exceeds ~24.96.

Integrating it numerically — code

In practice you integrate the equations directly. The single subtlety is the 1/r² force, which blows up during close encounters and demands either tiny time steps or a symplectic scheme that conserves energy over long runs. Here is a compact velocity-Verlet (leapfrog) integrator — symplectic, second order, and the standard starting point for gravitational N-body work:

// Three-body integrator (velocity Verlet, symplectic)
const G = 1;                       // work in units where G = 1
const bodies = [
  { m: 1, p: [-1,  0, 0], v: [ 0.347, 0.533, 0] },  // figure-eight ICs
  { m: 1, p: [ 1,  0, 0], v: [ 0.347, 0.533, 0] },
  { m: 1, p: [ 0,  0, 0], v: [-0.694,-1.066, 0] },
];

function accelerations(bs) {
  const a = bs.map(() => [0, 0, 0]);
  for (let i = 0; i < bs.length; i++) {
    for (let j = 0; j < bs.length; j++) {
      if (i === j) continue;
      const dx = bs[j].p[0] - bs[i].p[0];
      const dy = bs[j].p[1] - bs[i].p[1];
      const dz = bs[j].p[2] - bs[i].p[2];
      const r2 = dx*dx + dy*dy + dz*dz + 1e-9;   // softening avoids blow-up
      const inv = G * bs[j].m / (r2 * Math.sqrt(r2));
      a[i][0] += dx * inv;
      a[i][1] += dy * inv;
      a[i][2] += dz * inv;
    }
  }
  return a;
}

function step(bs, dt) {
  let a = accelerations(bs);
  for (let i = 0; i < bs.length; i++)        // half-kick + drift
    for (let k = 0; k < 3; k++) {
      bs[i].v[k] += 0.5 * a[i][k] * dt;
      bs[i].p[k] += bs[i].v[k] * dt;
    }
  a = accelerations(bs);
  for (let i = 0; i < bs.length; i++)        // second half-kick
    for (let k = 0; k < 3; k++)
      bs[i].v[k] += 0.5 * a[i][k] * dt;
}

// Total energy — should stay nearly constant for a good integrator
function energy(bs) {
  let ke = 0, pe = 0;
  for (let i = 0; i < bs.length; i++) {
    ke += 0.5 * bs[i].m * (bs[i].v[0]**2 + bs[i].v[1]**2 + bs[i].v[2]**2);
    for (let j = i + 1; j < bs.length; j++) {
      const dx = bs[j].p[0]-bs[i].p[0], dy = bs[j].p[1]-bs[i].p[1], dz = bs[j].p[2]-bs[i].p[2];
      pe -= G * bs[i].m * bs[j].m / Math.sqrt(dx*dx + dy*dy + dz*dz);
    }
  }
  return ke + pe;
}

const E0 = energy(bodies);
for (let n = 0; n < 100000; n++) step(bodies, 0.001);
console.log('energy drift:', (energy(bodies) - E0).toExponential(2));  // tiny

// Demonstrate chaos: perturb body 1 by one part in 10^9 and re-run
const twin = JSON.parse(JSON.stringify(bodies));
twin[0].p[0] += 1e-9;
for (let n = 0; n < 100000; n++) { step(bodies, 0.001); step(twin, 0.001); }
const sep = Math.hypot(
  bodies[0].p[0]-twin[0].p[0],
  bodies[0].p[1]-twin[0].p[1]);
console.log('separation grew from 1e-9 to', sep.toFixed(3)); // exponential blow-up

Run that perturbation experiment and the two "identical" systems, started a billionth of a unit apart, end up on visibly different trajectories. That is sensitive dependence made concrete.

Where the three-body problem shows up

• Mission design. Lagrange-point halo orbits (JWST at L2, SOHO at L1), gravity assists, and low-energy "ballistic capture" transfers all live inside the restricted three-body model.
• Solar System stability. Laskar's integrations give the inner planets a Lyapunov time near 5 million years and a ~1% chance of Mercury destabilizing over the Sun's remaining lifetime.
• Star clusters and galaxies. The full N-body problem is the three-body problem writ large; close triple encounters eject stars and harden binaries.
• Hierarchical triples. A close binary orbited by a distant third star drives Kozai–Lidov oscillations, swapping eccentricity and inclination — relevant to hot Jupiters and merging black holes.
• Chaos theory itself. Poincaré's three-body work is the historical root of dynamical-systems theory, strange attractors, and the modern study of chaos.
• Asteroid and Trojan dynamics. Bodies trapped near L4/L5, resonances, and the Kirkwood gaps in the asteroid belt.

Common misconceptions

• "No solution exists." Wrong — Sundman's series is a solution, and special families are exact. What's missing is a general closed form in elementary functions.
• "Chaos means random." The system is fully deterministic. Chaos is sensitive dependence, not randomness — the same initial conditions always give the same trajectory.
• "More bodies, more chaos." Chaos already appears with three; the obstacle is having too few conserved integrals, and that obstruction is already present at N = 3.
• "Lagrange points are gravity-free." They are points where gravity plus the centrifugal force balance in the rotating frame, not places with no gravity.
• "The Solar System is provably stable." It is only marginally, probabilistically stable; the planets are weakly chaotic over million-year timescales.
• "Energy conservation pins down the orbit." It constrains the motion to an energy surface but leaves a high-dimensional chaotic sea on that surface.