Measuring Exactly How Unpredictable Three Gravitating Bodies Are

Three Bodies, One Force, No Formula

Two gravitating bodies solve exactly. Newton did it, and every physics student redoes it: an ellipse, a period, a closed and perfectly predictable curve running forever. Add one more body, keep the same single force, the same clean inverse-square law, and the exact solution vanishes. Not approximately. Not for lack of a clever substitution. Henri Poincaré proved in 1890 that the three-body problem admits no general closed-form solution of the kind Newton found for two [1], and the reason he gave is now one of the founding facts of chaos theory: the equations are perfectly deterministic and the outcomes are not usably predictable, because trajectories that start almost on top of each other pull apart at a rate that compounds. Good modern surveys of everything that sentence is shorthand for exist [10][17], and reading one is not a substitute for measuring the rate yourself, which is what the rest of this article does.

That is the sentence every popular account of chaos reaches for. It is also, on its own, a claim with no number attached. How fast do two nearby three-body trajectories actually separate? Fast compared to what? For how long can a prediction be trusted before it is worth exactly nothing? Those are measurement questions, and the club has no orbiting bodies to point a telescope at, so we built the orbiting bodies instead: a gravity simulator, validated twice against solutions nobody disputes, then pointed at the question everyone quotes and almost nobody quantifies.

What follows builds the integrator, checks it twice against answers nobody argues with, measures three real Lyapunov exponents and one prediction horizon, runs 1,020 Monte Carlo trials to see how ejection depends on energy, reports the one place a first attempt at this failed outright, and closes with the sharpest objection we could bring against our own numbers.

What We Built, and What It Cannot See

Three point masses, mutual Newtonian gravity, $G=1$, confined to a plane. Body $i$ accelerates as

Two dimensions rather than three is a real restriction, not a shortcut taken for convenience alone: a planar three-body system already has four degrees of freedom of relative motion, which is enough for full chaos, and every figure below is drawable because of it. A genuinely three-dimensional triple would differ quantitatively. We expect it would not differ in kind, and we have not tested that expectation.

Two published, closed orbits validate the machine before any random number enters the picture (§3): the two-body Kepler ellipse and the figure-eight choreography discovered by Moore [2] and proven to exist by Chenciner and Montgomery [3]. Neither orbit ever passes near a close encounter, so a fixed time step advanced by a 4th-order symplectic integrator, built the standard way [7] by composing three velocity-Verlet leapfrog half-steps with weights $w_1 = 1/(2-2^{1/3})$ and $w_0 = 1-2w_1$, handles both beautifully. Every validation number in §3 comes from that fixed-step integrator with no softening at all: the bare $1/r^2$ law, the one Newton wrote down.

Generic three-body motion is not so polite. Two bodies can swing arbitrarily close together, the force blows up as $1/r^2$, and a fixed step simply cannot resolve how fast that force is changing during the swing. §5 tells the story of finding this out by having it fail loudly. From §4 onward we use an adaptive step tied to the local two-body free-fall time at whichever pair is currently closest,

capped at a maximum step so that a system with no close pairs is not integrated needlessly finely. The step size now depends on the state, which means the method is no longer exactly symplectic in the strict sense: a truly symplectic adaptive scheme needs a time transformation, of the kind Kustaanheimo and Stiefel built for the two-body problem in 1965 [9] and that Wisdom and Holman later adapted for planetary N-body work [8], and we did not implement one. What we did instead is print the conserved quantities on every run and let the numbers argue for themselves, which §3 and §4 do.

A body is called ejected once it sits beyond twenty length units from the center of mass of the other two, still receding, with positive specific orbital energy relative to that pair treated as a single point mass. That is a standard simplification in scattering experiments, not the true three-body energy partition, and §9 names where it could mislead at the margin.

Checking the Machine Before Trusting It

Two validations ran before any random triangle was drawn, because an integrator with a sign error or a missing factor of two produces confident, internally consistent, wrong numbers all the way through a study like this one, and nothing about the later randomness would ever reveal it on its own.

Validation 1: the two-body Kepler orbit. Two equal masses on an eccentric ellipse ($e=0.5$), a third, massless body parked far away so the exact same three-body code runs with no special-casing. Twenty thousand steps per period, twenty-five sample points compared against the closed-form solution of Kepler's equation $E - e\sin E = M$, solved by Newton's method to $10^{-13}$. Worst relative position error over the whole orbit: 5.06×10⁻¹². Worst velocity error: 3.87×10⁻¹². Energy after one full period agreed with energy at the start to 7.5×10⁻¹⁴, close to the floor of double-precision arithmetic itself.

Validation 2: the figure-eight choreography. Three equal masses chase each other around a single figure-eight curve, a genuine periodic solution of the unrestricted three-body problem rather than a disguised two-body case [2][3]. We used the eight-decimal initial condition in wide circulation and commonly quoted alongside a period near 6.3259 [4]. Integrating for exactly that literature period left a residual of 1.5×10⁻³, too large to call a clean return, and small enough to suspect the published digits rather than our arithmetic, since energy conservation over the same stretch was already sitting at 4.6×10⁻¹⁴. We measured our own return time directly: integrated a little past the nominal period, recorded the distance back to the start at every step nearby, and fit a parabola through the three closest points. That gave T = 6.32656, a difference of +0.0101% from the literature figure, consistent with the widely circulated initial condition being accurate to only eight or nine significant figures rather than being exact. Re-integrating cleanly at this self-measured period gave a residual of 1.4×10⁻⁴, which we call a pass against a tolerance of 5×10⁻⁴, a tolerance set by the initial condition's own precision, not by the integrator: over twenty full periods (126.5 time units) the relative energy drift never exceeded 8.5×10⁻¹⁴, and it does not grow with time, which is exactly the signature a symplectic method is supposed to leave.

Plain Arithmetic

No simulation in this section, only the arithmetic a Lyapunov exponent is actually made of, because the phrase "grows exponentially" gets thrown around a great deal more often than it gets worked out on paper.

Take the random-triple median exponent from §6, λ = 1.326 time units⁻¹. A separation growing as $e^{\lambda t}$ doubles whenever $\lambda t = \ln 2$, so the doubling time is

Half a time unit to double. Ten doublings, $2^{10}=1024$, take $10 \times 0.523 = 5.23$ time units and turn one part in a thousand of initial uncertainty into essentially the whole system. Twenty doublings, $2^{20}\approx 1.05\times10^6$, take 10.46 time units, which is, up to rounding, exactly the ten-and-a-bit-time-unit prediction horizon §7 reports for a one-in-a-million starting error. That agreement is not a coincidence dressed up as one: the horizon formula in §7, $t_{\text{horizon}}=(1/\lambda)\ln(L/\delta)$, is precisely $(1/\lambda)$ multiplied by $\ln(10^6)$, and $\ln(10^6)=6\ln10=13.82$, so $t_{\text{horizon}} = 13.82/1.326 = 10.42$. Two ways of doing the same multiplication landed on the same answer, which is the entire point of plain arithmetic sections: if they hadn't matched, something upstream would be wrong.

Compare the figure-eight's own rate, λ = 0.0374. The same twenty-doubling stretch, $20\ln2/0.0374$, takes 370 time units, fifty-eight full periods of the orbit itself. A random triple forgets its starting condition roughly thirty-five times faster than the exact figure-eight solution does. Same force law, same code, same two bodies of arithmetic, wildly different answers, because the underlying dynamics are not the same at all.

Working Notes From the Club Table

Three Configurations, Three Very Different Answers

A Lyapunov exponent is measured here the standard way, due to Benettin and colleagues [11]: run a reference trajectory and a copy nudged by $\delta_0=10^{-8}$ in a random direction of the full twelve-dimensional phase space, advance both by the same adaptive steps, and every half a time unit measure how far apart they have grown, accumulate $\ln(D/\delta_0)$, then rescale the perturbed copy back down to $\delta_0$ along the same direction before continuing. The running sum divided by elapsed time converges to $\lambda$.

The figure-eight itself, §3's periodic solution. Ten independent perturbation directions, each carried for twenty-five periods (158 time units). Mean λ = 0.0374 ± 0.0014 (SEM, n=10). A Lyapunov time of nearly 27 time units, over four periods of the orbit's own motion. This is not zero. The figure-eight is known to be linearly stable in the strict Floquet sense [3][4], and a genuinely non-chaotic periodic orbit would return a rate indistinguishable from zero within our numerical noise. Ours does not: a small, real, positive rate survives a perturbation off the exact periodic solution, which is a much weaker and more honest statement than "the figure-eight is chaotic" and a much stronger one than "it is stable."

The Pythagorean problem, §5's cautionary tale. Masses 3, 4 and 5 at rest at the vertices of a 3-4-5 right triangle, first integrated by Burrau in 1913 [5] and carried to its dramatic conclusion by Szebehely and Peters in 1967 [6]: repeated close encounters and an eventual violent ejection. Ten perturbation directions to t=65, adaptive step throughout: mean λ = 0.472 ± 0.007 (SEM, n=10), a Lyapunov time under three time units. The closest triple encounter across all ten runs came within 0.0004 length units, a hundred times tighter than a fixed step at dt=0.001 could survive without the energy error already reported in §5.

Published numbers to test our own against are scarce for this specific measurement, since most three-body chaos results in the literature are qualitative demonstrations rather than quoted exponents, but Boekholt and Portegies Zwart's survey of N-body integration reliability is the closest quantitative touchstone available: running the Pythagorean problem on many independent codes and precisions, they find the numerical trajectories diverge from each other at a rate set by the encounter's own dynamical timescale, order the local free-fall time, while still agreeing on the qualitative outcome [15]. Our own Pythagorean Lyapunov time of 2.12 time units, against a free-fall time at the closest encounter of a small fraction of a time unit, sits at exactly that order of magnitude, the sanity check a real bug would have failed. For scale in the other direction: the real Solar System has a measured Lyapunov time of roughly five million years [13], later confirmed with independent long-term integrations [14]. That is not a different phenomenon from what we measured here in a few time units. It is the same instability at a vastly gentler setting, because the planets sit enormously farther from a genuine close encounter, relative to their own orbital scales, than three freshly dropped random bodies ever do [16].

A Number That Refuses to Behave

Fourteen of the forty random trials exceeded twice the median. Two sat above thirteen. If those forty numbers came from a well-behaved bell curve, that should not happen; a proper Gaussian puts almost nothing two standard deviations out, let alone ten times the typical value. Instead the distribution has a heavy right tail, and the tail has a specific physical cause: a trial's exponent is a time average, and if the trajectory happens to pass unusually close during a triple encounter partway through, that one brief, violent stretch of exponential stretching dominates the whole average, the way one lottery winner can dominate the mean income of an entire street.

The mean of forty such numbers is not a stable estimate of anything. Add or remove one extreme trial and it moves by a large fraction of itself; our own mean, 2.545, sits nearly twice the median. The median does not have that problem, because it only asks which value sits in the middle of a sorted list, and a handful of enormous values cannot drag a rank-based statistic around the way they drag an average. That is why §6's table reports both and why every later section that needs one number uses 1.326, the median, not 2.545, the mean.

How Far Ahead Can You See

A Lyapunov exponent is not by itself a forecast. Turned around, it gives a horizon: the time before an initial measurement error of size $\delta$ grows to the size of the system itself, $L$,

Using the random-triple median, $\lambda=1.326$, and $L=1$ in these units:

Every row down that first column trades four more orders of magnitude of measurement precision for a fixed 3.47 more time units of foresight. That is the logarithm doing exactly what a logarithm does, and it is the whole reason chaos defeats brute-force precision as a strategy: going from one-in-a-hundred to one-in-a-trillion, a ten-billion-fold improvement in how well you know where you started, buys the same fixed increment of extra warning as going from one-in-a-hundred to one-in-ten-thousand did. Nobody builds a measuring instrument ten billion times better to buy a linear return like that.

Who Gets Thrown Out, and When

A second question, separate from how fast nearby trajectories diverge: given a random bound triple, how often does the system eventually throw one body out for good? We swept a single parameter, $\alpha$, the ratio of each body's initial speed to the local virial speed, from 0.3 (cold, deeply bound) to 1.3 (already unbound), drew 170 independent random triangles at each of six settings, and integrated each for up to 20 time units, checking at every step whether one body had crossed twenty length units from the other two's center of mass while still receding with positive relative energy.

Two things stand out. First, the bound systems ($\alpha \le 0.9$, all E<0) do not show a clean monotonic trend with energy: the fraction wanders between 15% and 26% with overlapping error bars, which given only 170 trials per bin is consistent with a genuinely flat or weakly-varying relationship in that regime, not with a strong one we failed to resolve. Second, the picture changes sharply as energy crosses zero: 62% at $\alpha=1.1$ (still formally bound, E=−0.74, but close to the threshold) and 81% once unbound outright. That is the expected physics: an already-unbound system barely needs the chaotic dynamics to separate a body, while a deeply bound one needs a rare, favourable sequence of close encounters to hand one body enough kinetic energy to escape at all, and 20 time units is not always enough for that rare sequence to occur.

That last clause matters. Every fraction in this table is a finite-time measurement bounded at 20 time units, not a claim about eventual fate. Classical results in celestial mechanics hold that almost every bound, aperiodic triple system eventually disrupts given unlimited time; a system that survives our window is reported here as "not yet ejected by t=20", never as "permanently stable," and the true asymptotic ejection fraction for the bound bins is almost certainly higher than what a 20-time-unit window can show.

The Strongest Objection We Could Make

Here is the sharpest challenge available against this whole study, and it goes after the method rather than any single number: from §2 onward, everything quantitative in this article ran on an integrator that has abandoned the one property, exact symplecticity, that its own docstring promised at the start. We wrote a 4th-order symplectic method, validated it beautifully, then switched to a state-dependent adaptive step the moment real chaos showed up, and a symplectic map with a step size that depends on the state is no longer, strictly, a symplectic map. If the whole point of using a Yoshida integrator was long-term conservation, and we threw that guarantee away exactly when the dynamics got interesting, what is left standing?

Two answers, and we think both are honest rather than convenient. First, the guarantee we gave up was never really available for this problem in the first place: a fixed-step symplectic method conserves its shadow Hamiltonian only while the force stays smooth on the scale of one step, and a close three-body encounter violates that constantly, which is exactly what §5's 350-fold energy blowup demonstrates. Trading an unavailable exact guarantee for a measured one is not a downgrade if the exact guarantee was never actually being delivered. Second, and more fundamentally, exact trajectory-following is not what a chaotic simulation needs to be useful at all. Sauer, Grebogi and Yorke showed that a numerically computed chaotic trajectory, even one that has already diverged completely from the true orbit with those exact initial conditions, typically still shadows some other true trajectory of the same system closely for long stretches, a result they titled, pointedly, "How long do numerical chaotic solutions remain valid?" [12]. What a chaos study can honestly claim is not "this is the one true orbit" but "this is a real orbit of the same system, and the statistics we measured from many such orbits, Lyapunov exponents, ejection rates, are properties of the system rather than of one lucky or unlucky numerical path." Every headline number in this article is exactly that kind of statistic: an average or a median over many independent trials, never a claim about one irreplaceable trajectory.

That defense has a real limit, and we want to name it rather than let the citation do all the work. Shadowing results guarantee a nearby true orbit exists; they do not, by themselves, guarantee that our specific choice of adaptive rule and safety length reproduces the correct statistical distribution of outcomes across an ensemble, only that each individual trajectory is not numerical garbage. §10's sensitivity check, where a coarser step moved the ejection fraction from 0.22 to 0.34, is the closer, harder test of that distributional question, and it did not come back clean.

Where a Different Choice Would Change the Answer

Four choices in this study could plausibly have gone another way, and three of them measurably move the answer.

The step-size fraction η. Our default, η=0.005, gave an ejection fraction of 0.218 at α=0.9. A four-times-coarser step (η=0.02) gave 0.344 ± 0.050; a four-times-finer step (η=0.00125) gave 0.278 ± 0.047. Neither extreme agrees with the default within its own error bar, and the finer step did not simply converge onto the coarser one from the other side, which tells us 90 trials per setting is not yet enough to fully pin down where this number settles as η→0. This is, by a wide margin, the largest source of uncertainty in the whole Monte Carlo section, larger than any of the binomial standard errors reported in §9's table, and we would treat any of this study's ejection fractions as accurate to a factor of order 1.5, not to the two significant figures the table displays.

The safety length. Changing the numerical safety floor by a factor of 100 in either direction (10⁻⁴ and 10⁻⁸ in place of the default 10⁻⁶) moved the same fraction only to 0.256 and 0.278, both inside the default's own standard error. That confirms the safety length is doing what its name says, guarding against literal division by zero, and is not secretly acting as a physical softening the way the failed fixed-step attempt in §5 used one.

Two dimensions instead of three. Every result above is for planar motion. A body's escape route in three dimensions has one more direction to take, and encounters that would be a direct hit in the plane can pass at a nonzero perpendicular offset in three dimensions, generally softening the closest approaches somewhat. We expect this would lower both the ejection fractions and the closest recorded encounter distances, and we have not measured it.

The escape criterion's energy partition. Calling a body "ejected" once it clears twenty length units, receding, with positive energy relative to the other pair treated as one point mass, is standard practice in scattering experiments and an approximation nonetheless: right at the escape threshold, a body can cross that radius, still be weakly coupled to the pair it left, and later return. A stricter, more expensive criterion, tracking the full three-body energy partition until it stabilises, would likely reclassify a small number of borderline α=1.1 trials, the bin closest to the threshold, and would not meaningfully change the deeply bound or clearly unbound bins at either end of Figure 4.

Reproducing this

Everything above comes from one file and one command. Python 3.12 and numpy are required; no other package is imported.

python three-body-chaos.py > three-body-chaos-output.txt

Expect 5 to 8 minutes on a laptop; our archived run took 366.2 seconds. Master seed 20260914 is hard-coded at the top of the file, and every random stream, initial triangles, velocity directions, perturbation directions, is spawned from it through numpy's SeedSequence, so a rerun reproduces every number in this article exactly. The interactive model lets you drag the figure-eight's own starting positions and velocities, watch a faint second copy peel away from the first, and slide the initial uncertainty on a prediction-horizon calculator built from this study's own measured exponent.

References

python three-body-chaos.py > three-body-chaos-output.txt

1. 1.

   Poincaré, H. (1890). Sur le problème des trois corps et les équations de la dynamique. Acta Mathematica 13, 1–270.

2. 2.

   Moore, C. (1993). Braids in classical dynamics. Physical Review Letters 70(24), 3675–3679. doi:10.1103/PhysRevLett.70.3675

3. 3.

   Chenciner, A. & Montgomery, R. (2000). A remarkable periodic solution of the three-body problem in the case of equal masses. Annals of Mathematics 152(3), 881–902. doi:10.2307/2661357

4. 4.

   Simó, C. (2002). Dynamical properties of the figure eight solution of the three-body problem. In Contemporary Mathematics: Celestial Mechanics 292, 209–228. American Mathematical Society. doi:10.1090/conm/292/04926

5. 5.

   Burrau, C. (1913). Numerische Berechnung eines Spezialfalles des Dreikörperproblems. Astronomische Nachrichten 195(6), 113–118. doi:10.1002/asna.19131950602

6. 6.

   Szebehely, V. & Peters, C. F. (1967). Complete solution of a general problem of three bodies. The Astronomical Journal 72, 876–883. doi:10.1086/110355

7. 7.

   Yoshida, H. (1990). Construction of higher order symplectic integrators. Physics Letters A 150(5–7), 262–268. doi:10.1016/0375-9601(90)90092-3

8. 8.

   Wisdom, J. & Holman, M. (1991). Symplectic maps for the n-body problem. The Astronomical Journal 102, 1528–1538. doi:10.1086/115978

9. 9.

   Kustaanheimo, P. & Stiefel, E. (1965). Perturbation theory of Kepler motion based on spinor regularization. Journal für die reine und angewandte Mathematik 218, 204–219. doi:10.1515/crll.1965.218.204

10. 10.

    Musielak, Z. E. & Quarles, B. (2014). The three-body problem. Reports on Progress in Physics 77, 065901. doi:10.1088/0034-4885/77/6/065901

11. 11.

    Benettin, G., Galgani, L., Giorgilli, A. & Strelcyn, J. M. (1980). Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Parts 1 and 2. Meccanica 15, 9–20 and 21–30. doi:10.1007/BF02128236, doi:10.1007/BF02128237

12. 12.

    Sauer, T., Grebogi, C. & Yorke, J. A. (1997). How long do numerical chaotic solutions remain valid? Physical Review Letters 79, 59–62. doi:10.1103/PhysRevLett.79.59

13. 13.

    Laskar, J. (1989). A numerical experiment on the chaotic behaviour of the Solar System. Nature 338(6212), 237–238. doi:10.1038/338237a0

14. 14.

    Hayes, W. B. (2007). Is the outer Solar System chaotic? Nature Physics 3, 689–691. doi:10.1038/nphys728

15. 15.

    Boekholt, T. & Portegies Zwart, S. (2015). On the reliability of N-body simulations. Computational Astrophysics and Cosmology 2, 2. doi:10.1186/s40668-014-0005-3

16. 16.

    Di Cintio, P. & Trani, A. A. (2025). Partial suppression of chaos in relativistic three-body problems. Astronomy & Astrophysics 693, A53. doi:10.1051/0004-6361/202452678

17. 17.

    Chenciner, A. (2007). Three body problem. Scholarpedia 2(10), 2111. doi:10.4249/scholarpedia.2111