\documentclass[11pt]{article}
\usepackage[margin=1in]{geometry}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{booktabs}
\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{compat=1.18}
\usepackage[colorlinks=true, linkcolor=black, citecolor=black, urlcolor=blue]{hyperref}

\definecolor{tbdeep}{HTML}{163F2C}
\definecolor{tbmid}{HTML}{276E49}
\definecolor{tbbright}{HTML}{3C9A68}
\definecolor{tbolive}{HTML}{7D8A3A}
\definecolor{tbgrey}{HTML}{7E8C84}

\title{Measuring Exactly How Unpredictable Three Gravitating Bodies Are}
\author{The Science Journaling Club}
\date{Volume 3, Issue 1 $\cdot$ Fall 2026 $\cdot$ Original Research}

\begin{document}
\maketitle

\begin{abstract}
Three point masses under Newtonian gravity are the textbook example of chaos,
and the claim is usually left asserted rather than measured. We built a
4th-order symplectic integrator, a composition of velocity-Verlet leapfrog
sub-steps with Yoshida's coefficients \cite{yoshida1990}, validated it against
the analytic two-body Kepler orbit to $5.06\times10^{-12}$ relative error and
against the published figure-eight periodic solution \cite{chencinermontgomery2000,simo2002}
to a residual of $1.4\times10^{-4}$ after one full period, then measured the
largest Lyapunov exponent of three configurations by Benettin's two-trajectory
renormalisation method \cite{benettin1980}. The exact figure-eight
choreography is almost not chaotic at all: $\lambda = 0.0374 \pm 0.0014$, a
Lyapunov time of nearly 27 time units, over four full periods. The 1913
Pythagorean three-body problem \cite{burrau1913,szebehely1967}, famous for a
violent triple encounter, gives $\lambda = 0.472 \pm 0.007$, a Lyapunov time
under three time units, but only once we replaced a fixed integration step
with one tied to the local free-fall time: the fixed step we tried first
reported an energy that had grown by a factor of 350, a failure this article
reports rather than hides. Forty randomly drawn bound triangles gave a
heavy-tailed spread of exponents, median 1.326, mean pulled up to 2.545 by
rare close encounters, translating to a prediction horizon of about 10.4 time
units before a one-part-in-a-million initial uncertainty grows to the size of
the system. Across 1{,}020 Monte Carlo trials at six initial energies, the
fraction of triples that ejected a body within 20 time units rose from
roughly 15--26\% for bound systems to 62\% near the escape threshold and 81\%
for unbound ones. A step-size sensitivity check moved the headline ejection
fraction from 0.22 to as high as 0.34 under a four-times-coarser step, the
single largest source of uncertainty in this study, reported plainly rather
than tuned away. Seed 20260914.
\end{abstract}

\section{Three bodies, one force, no formula}

Two gravitating bodies solve exactly, an ellipse and a period, closed and
predictable forever. Add a third body under the same single force and the
exact solution disappears: Poincar\'e proved in 1890 that the three-body
problem admits no general closed-form solution of the kind Newton found for
two \cite{poincare1890}, the founding fact of what became chaos theory.
Nearby trajectories separate at a rate that compounds, and that sentence is
usually left as a slogan. This article measures the rate. Good modern
surveys of the whole subject exist \cite{musielak2014,chenciner2007}; reading
one is not a substitute for measuring the exponent, which is what follows.

\section{What we built, and what it cannot see}

Three point masses, mutual Newtonian gravity, $G=1$, confined to a plane.
Body $i$ accelerates as
\begin{equation}
\ddot{\mathbf{x}}_i = \sum_{j \ne i} \frac{G m_j (\mathbf{x}_j - \mathbf{x}_i)}{r_{ij}^3}.
\label{eq:eom}
\end{equation}
A planar restriction is a real simplification, not a shortcut: it already
carries four degrees of freedom of relative motion, enough for full chaos,
and every figure below is drawable because of it. A genuinely
three-dimensional triple would differ quantitatively; we have not tested
that.

Two published, closed orbits validate the machine before any random number
enters (Section~\ref{sec:validate}): the two-body Kepler ellipse and the
figure-eight choreography discovered by Moore \cite{moore1993} and proven to
exist by Chenciner and Montgomery \cite{chencinermontgomery2000}. Neither
orbit passes near a close encounter, so a fixed-step 4th-order symplectic
integrator, built the standard way \cite{yoshida1990} by composing three
velocity-Verlet half-steps with weights $w_1=1/(2-2^{1/3})$,
$w_0=1-2w_1$, handles both with no softening at all.

Generic three-body motion is not so polite. Two bodies can swing arbitrarily
close together, the force diverges as $1/r^2$, and a fixed step cannot
resolve how fast that force is changing during the swing.
Section~\ref{sec:notes} tells the story of finding this out by having it fail
loudly. From Section~\ref{sec:lyapunov} onward we use an adaptive step tied
to the local two-body free-fall time at whichever pair is currently closest,
\begin{equation}
\Delta t = \eta \cdot 2\pi\sqrt{\frac{r_{\min}^3}{G(m_i+m_j)}}, \qquad \eta = 0.005,
\label{eq:adaptive}
\end{equation}
capped at a maximum step. The step size now depends on the state, so the
method is no longer exactly symplectic: a genuinely symplectic adaptive
scheme needs a time transformation, of the kind Kustaanheimo and Stiefel
built for the two-body problem \cite{kustaanheimo1965} and that Wisdom and
Holman later adapted for planetary work \cite{wisdomholman1991}, which we did
not implement. We printed the conserved quantities on every run instead and
let the numbers argue for themselves.

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, a
standard simplification in scattering experiments and an approximation to
the true three-body energy partition (Section~\ref{sec:choices}).

\section{Checking the machine before trusting it}
\label{sec:validate}

\textbf{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
same three-body code runs unmodified. 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 orbit: $5.06\times10^{-12}$. Worst
velocity error: $3.87\times10^{-12}$. Energy after one full period agreed
with energy at the start to $7.5\times10^{-14}$.

\textbf{Validation 2: the figure-eight choreography.} We used the
eight-decimal initial condition in wide circulation, commonly quoted
alongside a period near 6.3259 \cite{simo2002}. Integrating for exactly that
literature period left a residual of $1.5\times10^{-3}$, too large to call a
clean return, while energy conservation over the same stretch already sat at
$4.6\times10^{-14}$, which pointed at the published digits rather than our
arithmetic. We measured our own return time directly: integrated 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. Re-integrating cleanly at this
self-measured period gave a residual of $1.4\times10^{-4}$, a pass against a
tolerance of $5\times10^{-4}$ 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\times10^{-14}$ and did not grow
with time.

\begin{figure}[htbp]
\centering
\begin{tikzpicture}
\begin{axis}[
  width=12.5cm, height=8cm,
  xlabel={$x$}, ylabel={$y$},
  xmin=-1.2, xmax=1.2, ymin=-0.42, ymax=0.42,
  grid=both, grid style={gray!18},
  legend pos=outer north east, legend cell align=left,
  legend style={font=\footnotesize, draw=gray!40},
  tick label style={font=\footnotesize},
  label style={font=\small},
  axis equal image
]
\addplot[tbdeep, thick, no markers] table[row sep=\\] {
x y \\
0.97004 -0.24309 \\
1.01987 -0.18654 \\
1.05343 -0.12802 \\
1.07324 -0.06895 \\
1.08091 -0.00983 \\
1.07709 0.04929 \\
1.06146 0.10840 \\
1.03272 0.16720 \\
0.98864 0.22465 \\
0.92623 0.27829 \\
0.84243 0.32341 \\
0.73574 0.35239 \\
0.60892 0.35574 \\
0.47062 0.32616 \\
0.33222 0.26366 \\
0.20139 0.17568 \\
0.07894 0.07247 \\
-0.03938 -0.03649 \\
-0.15982 -0.14260 \\
-0.28763 -0.23678 \\
-0.42402 -0.30884 \\
-0.56355 -0.34982 \\
-0.69540 -0.35693 \\
-0.80938 -0.33544 \\
-0.90080 -0.29473 \\
-0.97003 -0.24328 \\
-1.01987 -0.18674 \\
-1.05344 -0.12823 \\
-1.07326 -0.06917 \\
-1.08094 -0.01004 \\
-1.07713 0.04908 \\
-1.06151 0.10819 \\
-1.03278 0.16699 \\
-0.98872 0.22445 \\
-0.92633 0.27811 \\
-0.84255 0.32325 \\
-0.73587 0.35227 \\
-0.60905 0.35566 \\
-0.47073 0.32612 \\
-0.33231 0.26365 \\
-0.20145 0.17569 \\
-0.07898 0.07248 \\
0.03937 -0.03647 \\
0.15983 -0.14257 \\
0.28766 -0.23675 \\
0.42407 -0.30880 \\
0.56362 -0.34976 \\
0.69548 -0.35683 \\
0.80946 -0.33531 \\
0.90087 -0.29458 \\
0.97010 -0.24312 \\
};
\addlegendentry{body 1 (bodies 2, 3 trace the same curve)}
\addplot[only marks, mark=*, mark size=2pt, tbolive] coordinates {(0.97004,-0.24309)};
\addlegendentry{$t=0$}
\end{axis}
\end{tikzpicture}
\caption{The figure-eight three-body orbit, one self-measured period
($T=6.32656$), plotted from this study's own validated integrator. All three
equal masses trace the identical curve, offset from each other in time by a
third of a period.}
\label{fig:fig8}
\end{figure}

\section{Plain arithmetic}

Take the random-triple median exponent from Section~\ref{sec:lyapunov},
$\lambda=1.326$ time units$^{-1}$. A separation growing as $e^{\lambda t}$
doubles when $\lambda t=\ln 2$, so the doubling time is
\begin{equation}
t_{2\times} = \frac{\ln 2}{\lambda} = \frac{0.6931}{1.326} = 0.523 \text{ time units.}
\end{equation}
Ten doublings, $2^{10}=1024$, take $5.23$ time units and turn one part in a
thousand of initial uncertainty into essentially the whole system. Twenty
doublings, $2^{20}\approx1.05\times10^6$, take $10.46$ time units, matching
(up to rounding) the prediction horizon Section~\ref{sec:horizon} reports for
a one-in-a-million starting error: $t_{\text{horizon}}=(1/\lambda)\ln(10^6)
= 13.82/1.326 = 10.42$. Two routes to the same multiplication landed on the
same answer, which is the point of doing it twice.

The figure-eight's own rate, $\lambda=0.0374$, needs $20\ln2/0.0374=370$ time
units for the same twenty doublings, fifty-eight full periods of the orbit
itself. A random triple forgets its starting condition roughly 35 times
faster than the exact figure-eight solution does, from the same force law and
the same code.

\section{Working notes from the club table}
\label{sec:notes}

\noindent\textbf{Session 1.} Validated the two-body orbit in an afternoon.
Figure-eight residual came back at $1.5\times10^{-3}$ against a hoped-for
$10^{-6}$; energy conservation over the same run was already at $10^{-14}$,
which does not happen if the physics is wrong, so we suspected the published
digits rather than the code.

\noindent\textbf{Session 2.} Wrote a parabolic-fit period search instead of
guessing further. It found a period 0.01\% off the literature value: eight
significant figures of initial condition cannot close a loop to one part in a
million. Kept the whole detour in the article.

\noindent\textbf{Session 3.} First Monte Carlo run on the Pythagorean problem
with a fixed step and a small softening length: energy had grown by a factor
of 350 by $t=30$. Shrinking the step by ten made it worse, 887$\times$, the
tell that a fixed step near a close encounter is the wrong tool rather than a
slowly converging one. Switched to the adaptive step of
Equation~\ref{eq:adaptive}; same problem, same close encounter (0.0004 units),
energy drift $8.8\times10^{-4}$, using fewer total steps.

\noindent\textbf{Session 4.} Forty random Lyapunov trials: mean 2.54, median
1.33, two trials above 13. Nearly called the run broken a third time before
looking at the actual distribution instead of the summary statistic.

\noindent\textbf{Session 5.} Full run, seed 20260914, 366 seconds end to end.
The sensitivity check in Section~\ref{sec:choices} moved the headline
ejection fraction from 0.22 to 0.34 under a coarser step, a gap larger than
any statistical error bar in the study.

\section{Three configurations, three very different answers}
\label{sec:lyapunov}

A Lyapunov exponent is measured the standard way, due to Benettin and
colleagues \cite{benettin1980}: advance a reference trajectory and a copy
nudged by $\delta_0=10^{-8}$ in a random direction of the full
twelve-dimensional phase space, and every half a time unit measure the
separation $D$, accumulate $\ln(D/\delta_0)$, then rescale the perturbed copy
back to $\delta_0$ along the same direction before continuing. The running
sum divided by elapsed time converges to $\lambda$.

\begin{table}[htbp]
\centering
\caption{Largest Lyapunov exponent for three configurations.}
\label{tab:lyapunov}
\begin{tabular}{lrrr}
\toprule
configuration & $\lambda$ & Lyapunov time & $n$ \\
\midrule
figure-eight (exact periodic orbit) & $0.0374\pm0.0014$ & 26.8 & 10 \\
Pythagorean problem \cite{burrau1913,szebehely1967} & $0.472\pm0.007$ & 2.12 & 10 \\
random bound triples, median & 1.326 & 0.754 & 40 \\
random bound triples, mean & $2.545\pm0.478$ & 0.393 & 40 \\
\bottomrule
\end{tabular}
\end{table}

\textbf{The figure-eight itself.} Ten independent perturbation directions,
each carried for 25 periods (158 time units): mean $\lambda=0.0374\pm0.0014$
(SEM). The figure-eight is known to be linearly stable in the strict Floquet
sense \cite{chencinermontgomery2000,simo2002}, 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 solution, weaker than ``the figure-eight is
chaotic'' and stronger than ``it is stable.''

\textbf{The Pythagorean problem.} Masses 3, 4 and 5 at rest on a 3-4-5 right
triangle, first integrated by Burrau \cite{burrau1913} and carried to its
dramatic conclusion by Szebehely and Peters \cite{szebehely1967}. Ten
perturbation directions to $t=65$: mean $\lambda=0.472\pm0.007$ (SEM), a
Lyapunov time under three time units, closest triple encounter $0.0004$
length units.

Published numbers to compare against are scarce for this exact measurement.
Boekholt and Portegies Zwart's survey of N-body integration reliability finds
that in Pythagorean-type close encounters, independent codes and precisions
diverge from each other at a rate set by the encounter's own dynamical
timescale while still agreeing on the qualitative outcome
\cite{boekholt2015}. Our Lyapunov time of 2.12 against a closest-encounter
free-fall time far below one time unit sits at exactly that order of
magnitude. For scale, the real Solar System has a measured Lyapunov time near
five million years \cite{laskar1989}, later confirmed with independent
long-term integrations \cite{hayes2007}, the same instability operating on
planets that sit enormously farther from a genuine close encounter, relative
to their own orbital scales, than three freshly dropped random bodies ever do
\cite{dicintio2025}.

\textbf{Forty random bound triangles.} Mean energy $-2.86$ in these units,
all forty finite over 35 time units: mean $2.545\pm0.478$ (SEM), median
$1.326$, range $0.198$ to $14.468$. Fourteen of forty exceeded twice the
median; two exceeded thirteen. The distribution has a heavy right tail
because a trial's exponent is a time average, and one unusually close triple
encounter partway through can dominate the whole average. The mean is not a
stable estimate of anything here; every later section that needs one number
uses the median, $1.326$, not the mean.

\begin{figure}[htbp]
\centering
\begin{tikzpicture}
\begin{axis}[
  width=12.5cm, height=7.5cm,
  ylabel={Lyapunov exponent $\lambda$ (log scale)},
  ymode=log,
  symbolic x coords={figure-eight,Pythagorean,random (median)},
  xtick=data,
  ymin=0.01, ymax=20,
  grid=both, grid style={gray!18},
  tick label style={font=\footnotesize},
  label style={font=\small}
]
\addplot[only marks, mark=*, mark size=2.5pt, tbdeep,
  error bars/.cd, y dir=both, y explicit]
  coordinates {
  (figure-eight,0.0374) +- (0,0.0014)
  (Pythagorean,0.472) +- (0,0.007)
  };
\addplot[only marks, mark=*, mark size=2.5pt, tbbright,
  error bars/.cd, y dir=both, y explicit]
  coordinates {
  (random (median),1.326) +- (0,2.252)
  };
\end{axis}
\end{tikzpicture}
\caption{Largest Lyapunov exponent for three configurations, log scale.
Error bars are $\pm1$ SEM for the figure-eight and Pythagorean cases
($n=10$ each); the random-triple point shows the interquartile range
(0.756 to 3.578) rather than a standard error, since that spread is real
dynamical variation across trials.}
\label{fig:lyapunov}
\end{figure}

\section{How far ahead can you see}
\label{sec:horizon}

A Lyapunov exponent gives a prediction horizon: the time before an initial
error $\delta$ grows to the size of the system $L$,
\begin{equation}
t_{\text{horizon}} = \frac{1}{\lambda}\ln\!\left(\frac{L}{\delta}\right).
\label{eq:horizon}
\end{equation}
Using the random-triple median $\lambda=1.326$ and $L=1$:

\begin{table}[htbp]
\centering
\caption{Prediction horizon at the random-triple median exponent.}
\label{tab:horizon}
\begin{tabular}{rrr}
\toprule
$\delta$ & horizon (time units) & horizon / Lyapunov time \\
\midrule
$10^{-2}$ & 3.47 & 4.61 \\
$10^{-4}$ & 6.95 & 9.21 \\
$10^{-6}$ & 10.42 & 13.82 \\
$10^{-8}$ & 13.89 & 18.42 \\
$10^{-10}$ & 17.37 & 23.03 \\
$10^{-12}$ & 20.84 & 27.63 \\
\bottomrule
\end{tabular}
\end{table}

Every four more orders of magnitude of measurement precision buys a fixed
3.47 more time units of foresight. A ten-billion-fold improvement in how well
you know the starting state buys the same fixed increment as a
hundred-fold improvement did, because the horizon grows only as the
logarithm of precision.

\section{Who gets thrown out, and when}

A separate question: given a random bound triple, how often does the system
eventually eject a body? We swept $\alpha$, the ratio of each body's initial
speed to the local virial speed, from 0.3 to 1.3, drew 170 independent random
triangles at each of six settings, and integrated each for up to 20 time
units, checking whether a body crossed twenty length units from the other
two's center of mass while still receding with positive relative energy.

\begin{table}[htbp]
\centering
\caption{Ejection fraction by initial energy, 170 trials per row.}
\label{tab:ejection}
\begin{tabular}{rrrrr}
\toprule
$\alpha$ & mean $E$ & ejected & fraction & SE \\
\midrule
0.3 & $-3.107$ & 45 & 0.265 & 0.034 \\
0.5 & $-2.836$ & 26 & 0.153 & 0.028 \\
0.7 & $-2.211$ & 25 & 0.147 & 0.027 \\
0.9 & $-1.611$ & 37 & 0.218 & 0.032 \\
1.1 & $-0.741$ & 106 & 0.624 & 0.037 \\
1.3 & $+0.468$ & 138 & 0.812 & 0.030 \\
\bottomrule
\end{tabular}
\end{table}

\begin{figure}[htbp]
\centering
\begin{tikzpicture}
\begin{axis}[
  width=12.5cm, height=7cm,
  xlabel={mean initial energy $E$}, ylabel={ejection fraction},
  xmin=-3.5, xmax=1.0, ymin=0, ymax=1,
  grid=both, grid style={gray!18},
  tick label style={font=\footnotesize},
  label style={font=\small}
]
\addplot[tbdeep, thick, mark=*, mark size=2pt,
  error bars/.cd, y dir=both, y explicit]
  coordinates {
  (-3.107,0.265) +- (0,0.034)
  (-2.836,0.153) +- (0,0.028)
  (-2.211,0.147) +- (0,0.027)
  (-1.611,0.218) +- (0,0.032)
  (-0.741,0.624) +- (0,0.037)
  (0.468,0.812) +- (0,0.030)
  };
\addplot[dashed, tbgrey, thick, no markers] coordinates {(0,0) (0,1)};
\addlegendentry{$E=0$}
\end{axis}
\end{tikzpicture}
\caption{Fraction of 170 random triples per bin that ejected a body within
20 time units, against mean initial energy. Error bars are binomial standard
errors. Bound systems ($E<0$) sit in a noisy 15--26\% band before the
fraction climbs sharply once energy approaches and crosses zero.}
\label{fig:ejection}
\end{figure}

The bound bins do not show a clean monotonic trend: the fraction wanders
between 15\% and 26\% with overlapping error bars, consistent with a flat or
weakly-varying relationship given only 170 trials per bin. The picture
changes sharply once energy crosses zero: 62\% at $\alpha=1.1$ and 81\% once
unbound outright, the expected physics, since an unbound system barely needs
the chaotic dynamics to separate a body while a deeply bound one needs a rare
favourable sequence of close encounters. Every fraction here is a
finite-time measurement bounded at $t=20$, not a claim about eventual fate;
classical results hold that almost every bound, aperiodic triple eventually
disrupts given unlimited time, so a system surviving our window is reported
as ``not yet ejected,'' never as permanently stable.

\section{The strongest objection we could make}

The sharpest challenge here targets the method: from
Section~\ref{sec:lyapunov} onward, every quantitative result ran on an
integrator that abandoned exact symplecticity, the one property its own
construction promised. Two answers. First, the guarantee we gave up was never
really available for this problem: 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, exactly what
the 350-fold energy blowup in Section~\ref{sec:notes} demonstrates. Second,
exact trajectory-following is not what a chaotic simulation needs to be
useful. Sauer, Grebogi and Yorke showed that a numerically computed chaotic
trajectory typically still shadows some other true trajectory of the same
system closely for long stretches, even after diverging completely from the
one with those exact initial conditions \cite{sauer1997}. What a chaos study
can honestly claim is that the statistics measured from many such
trajectories, exponents, ejection rates, are properties of the system rather
than of one lucky or unlucky numerical path, which is exactly the kind of
statistic every headline number in this article is. That defense has a real
limit: shadowing guarantees a nearby true orbit exists, not that our specific
adaptive rule reproduces the correct statistical distribution of outcomes
across an ensemble. Section~\ref{sec:choices}'s sensitivity check, where a
coarser step moved the ejection fraction from 0.22 to 0.34, is the harder
test of that question, and it did not come back clean.

\section{Where a different choice would change the answer}
\label{sec:choices}

\textbf{The step-size fraction $\eta$.} Default $\eta=0.005$ gave an ejection
fraction of 0.218 at $\alpha=0.9$. A four-times-coarser step ($\eta=0.02$)
gave $0.344\pm0.050$; a four-times-finer step ($\eta=0.00125$) gave
$0.278\pm0.047$. Neither extreme agrees with the default within its own error
bar, and the finer step did not converge onto the coarser one from the other
side, so 90 trials per setting is not enough to fully pin down where this
number settles as $\eta\to0$. This is the largest source of uncertainty in
the whole Monte Carlo section, and we would treat any ejection fraction here
as accurate to a factor of order 1.5, not to the two significant figures the
tables display.

\textbf{The safety length.} Changing the numerical safety floor by
$100\times$ in either direction moved the same fraction only to 0.256 and
0.278, both inside the default's own standard error, confirming it guards
against division by zero rather than acting as physical softening.

\textbf{Two dimensions instead of three.} A body's escape route in three
dimensions has one more direction to take; encounters that are direct hits in
the plane can pass at a nonzero offset in three dimensions. We expect this
would lower both the ejection fractions and the closest recorded encounter
distances, and have not measured it.

\textbf{The escape criterion.} Calling a body ejected at twenty length units
with positive two-body relative energy is standard in scattering experiments
and an approximation: a body can cross that radius, remain weakly coupled,
and later return. A stricter criterion tracking the full three-body energy
partition would likely reclassify a small number of borderline
$\alpha=1.1$ trials without changing the bins at either end of
Figure~\ref{fig:ejection}.

\section{Reproducing this}

Everything above comes from one file and one command. Python 3.12 and numpy
are required; no other package is imported.
\begin{verbatim}
python three-body-chaos.py > three-body-chaos-output.txt
\end{verbatim}
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 is spawned from it through numpy's \texttt{SeedSequence}, so a rerun
reproduces every number in this article exactly.

\begin{thebibliography}{99}

\bibitem{poincare1890}
Poincar\'e, H. (1890). Sur le probl\`eme des trois corps et les \'equations
de la dynamique. \emph{Acta Mathematica} \textbf{13}, 1--270.

\bibitem{moore1993}
Moore, C. (1993). Braids in classical dynamics. \emph{Physical Review
Letters} \textbf{70}(24), 3675--3679.
\href{https://doi.org/10.1103/PhysRevLett.70.3675}{doi:10.1103/PhysRevLett.70.3675}

\bibitem{chencinermontgomery2000}
Chenciner, A. and Montgomery, R. (2000). A remarkable periodic solution of
the three-body problem in the case of equal masses. \emph{Annals of
Mathematics} \textbf{152}(3), 881--902.
\href{https://doi.org/10.2307/2661357}{doi:10.2307/2661357}

\bibitem{simo2002}
Sim\'o, C. (2002). Dynamical properties of the figure eight solution of the
three-body problem. In \emph{Contemporary Mathematics: Celestial Mechanics}
\textbf{292}, 209--228. American Mathematical Society.
\href{https://doi.org/10.1090/conm/292/04926}{doi:10.1090/conm/292/04926}

\bibitem{burrau1913}
Burrau, C. (1913). Numerische Berechnung eines Spezialfalles des
Dreik\"orperproblems. \emph{Astronomische Nachrichten} \textbf{195}(6),
113--118.
\href{https://doi.org/10.1002/asna.19131950602}{doi:10.1002/asna.19131950602}

\bibitem{szebehely1967}
Szebehely, V. and Peters, C. F. (1967). Complete solution of a general
problem of three bodies. \emph{The Astronomical Journal} \textbf{72},
876--883.
\href{https://doi.org/10.1086/110355}{doi:10.1086/110355}

\bibitem{yoshida1990}
Yoshida, H. (1990). Construction of higher order symplectic integrators.
\emph{Physics Letters A} \textbf{150}(5--7), 262--268.
\href{https://doi.org/10.1016/0375-9601(90)90092-3}{doi:10.1016/0375-9601(90)90092-3}

\bibitem{wisdomholman1991}
Wisdom, J. and Holman, M. (1991). Symplectic maps for the n-body problem.
\emph{The Astronomical Journal} \textbf{102}, 1528--1538.
\href{https://doi.org/10.1086/115978}{doi:10.1086/115978}

\bibitem{kustaanheimo1965}
Kustaanheimo, P. and Stiefel, E. (1965). Perturbation theory of Kepler
motion based on spinor regularization. \emph{Journal f\"ur die reine und
angewandte Mathematik} \textbf{218}, 204--219.
\href{https://doi.org/10.1515/crll.1965.218.204}{doi:10.1515/crll.1965.218.204}

\bibitem{musielak2014}
Musielak, Z. E. and Quarles, B. (2014). The three-body problem.
\emph{Reports on Progress in Physics} \textbf{77}, 065901.
\href{https://doi.org/10.1088/0034-4885/77/6/065901}{doi:10.1088/0034-4885/77/6/065901}

\bibitem{benettin1980}
Benettin, G., Galgani, L., Giorgilli, A. and 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. \emph{Meccanica}
\textbf{15}, 9--20 and 21--30.
\href{https://doi.org/10.1007/BF02128236}{doi:10.1007/BF02128236},
\href{https://doi.org/10.1007/BF02128237}{doi:10.1007/BF02128237}

\bibitem{sauer1997}
Sauer, T., Grebogi, C. and Yorke, J. A. (1997). How long do numerical
chaotic solutions remain valid? \emph{Physical Review Letters} \textbf{79},
59--62.
\href{https://doi.org/10.1103/PhysRevLett.79.59}{doi:10.1103/PhysRevLett.79.59}

\bibitem{laskar1989}
Laskar, J. (1989). A numerical experiment on the chaotic behaviour of the
Solar System. \emph{Nature} \textbf{338}(6212), 237--238.
\href{https://doi.org/10.1038/338237a0}{doi:10.1038/338237a0}

\bibitem{hayes2007}
Hayes, W. B. (2007). Is the outer Solar System chaotic? \emph{Nature
Physics} \textbf{3}, 689--691.
\href{https://doi.org/10.1038/nphys728}{doi:10.1038/nphys728}

\bibitem{boekholt2015}
Boekholt, T. and Portegies Zwart, S. (2015). On the reliability of N-body
simulations. \emph{Computational Astrophysics and Cosmology} \textbf{2}, 2.
\href{https://doi.org/10.1186/s40668-014-0005-3}{doi:10.1186/s40668-014-0005-3}

\bibitem{dicintio2025}
Di Cintio, P. and Trani, A. A. (2025). Partial suppression of chaos in
relativistic three-body problems. \emph{Astronomy \& Astrophysics}
\textbf{693}, A53.
\href{https://doi.org/10.1051/0004-6361/202452678}{doi:10.1051/0004-6361/202452678}

\bibitem{chenciner2007}
Chenciner, A. (2007). Three body problem. \emph{Scholarpedia} \textbf{2}(10),
2111.
\href{https://doi.org/10.4249/scholarpedia.2111}{doi:10.4249/scholarpedia.2111}

\end{thebibliography}

\end{document}