A prediction with an expiry date

Ask whether the solar system is stable and the honest answer has changed shape twice. Newton thought it was not, and expected divine adjustment; the problem he could not close is the one that has no closed form for three bodies. Laplace and Lagrange proved that it was, to first order in the planetary masses, and the proof stood for a century and a half. Then it was computed, and what came out was neither answer: the system is chaotic and bounded at the same time, which is a combination the question was not built to accept.

The two halves of that have to be taken separately, because they are about different quantities and the confusion between them is the whole difficulty.

Sensitivity is a rate, and it has a unit #

The definition of chaos that does any work is arithmetic: nearby trajectories separate exponentially. Two states differing by $\delta_0$ come to differ by $\delta_0 e^{\lambda t}$, and $\lambda$ — the Lyapunov exponent — has units of inverse time. Its reciprocal, the Lyapunov time, is the interval over which an error multiplies by $e$.

Everything follows from that one number by a logarithm. If the initial error is $\delta_0$ and the largest tolerable error is $\Delta$, the prediction stops being useful after

$$T = \frac{1}{\lambda}\,\ln\frac{\Delta}{\delta_0}.$$

The consequence people find surprising is buried in that logarithm. Improving the initial data buys time additively, not multiplicatively. Reduce $\delta_0$ by a factor of ten and $T$ grows by $\ln 10$ Lyapunov times — about 2.3 of them, regardless of where it started. For the inner planets, at five million years per e-folding, a hundredfold improvement in every measurement of every planetary position buys about twenty-three million years of extra reach. Not a hundredfold; twenty-three million years.

That is why the horizon cannot be pushed away by better instruments. To predict Mercury’s longitude at the age of the Sun, five billion years hence, would take $\ln(\Delta/\delta_0) = 1000$, and no measurement of anything has ever had $10^{434}$ significant figures.

The elements do not do what the phase does #

Now the other half, and it is what makes the solar system’s chaos survivable.

A planet’s state can be described by six numbers, and they divide sharply. Five of them — semi-major axis, eccentricity, inclination, and the two angles fixing the orbit’s orientation — describe the shape and placement of the ellipse. The sixth describes where the planet is on it, which is the phase.

Chaos in the solar system is overwhelmingly in the phase — the sixth number, the one Kepler’s equation exists to compute. The shape parameters wander too, but they wander within limits: the semi-major axes in particular are protected by a result that goes back to Laplace, since to first order in the masses they have no secular terms at all and can only oscillate. The elements that genuinely wander are the others: the eccentricities and inclinations are less well protected and do drift, over hundreds of millions of years, which is where the whole question actually lives.

So an integration of the solar system produces a statement of the form: this planet’s semi-major axis stays within a per cent of its present value for five billion years, and its longitude is unknown after about a hundred million. Both come out of one run and neither contradicts the other. The system is stable in the quantities that decide whether it is a solar system, and unpredictable in the quantity that decides where anything is.

Where the amplifier is #

Chaos in a gravitating system needs a mechanism, and the mechanism is resonance — specifically, the overlap of resonances.

An isolated resonance is not chaotic. Two bodies whose periods are in a small-integer ratio exchange energy in a way that repeats, and the resonant angle librates — swings back and forth about an equilibrium instead of circulating. That is a stable, bounded, entirely predictable motion, and it is what keeps the Galilean moons in their chain.

Bring a second resonance close enough that their zones of libration overlap, and a trajectory near the boundary of one is also near the boundary of the other. It can be captured by either, released, captured again — and which happens depends on the phase at the moment of encounter, which is precisely the quantity being amplified. Chirikov’s criterion makes the condition quantitative: chaos when the sum of the two resonance half-widths exceeds their separation.

In the inner solar system the relevant resonances are secular rather than mean-motion: not ratios of orbital periods, but ratios of the far slower precession rates of the orbits themselves. The one implicated in Mercury’s fate is $g_1 - g_5$, the difference between Mercury’s own apsidal precession and Jupiter’s, which happens to be very nearly zero. An angle whose rate is nearly zero can be pushed either way by very little, and Mercury’s eccentricity is what gets pushed.

What was actually measured #

Nothing about chaos is observed directly. What is observed is a set of positions, and everything else is inference through an integration — which makes it worth being exact about which numbers are data.

The Lyapunov time of the inner solar system, about five million years, is a computed quantity: two integrations from initial conditions differing by a chosen amount, and the slope of the logarithm of their separation. It has been computed independently by several groups with different integrators and different force models, and the agreement between them — a factor better than two — is the closest thing to a measurement available. The number is robust in the sense that it does not depend on which scheme produced it. It is not robust in the sense that it varies by a factor of a few depending on which planet’s which coordinate is used to define the separation.

The initial conditions are the data. Radar ranging to Venus and Mars, and Doppler tracking of orbiters at Mars, Jupiter and Saturn, give the inner planets’ positions to metres and the outer planets’ to kilometres. Metres, on orbits of $10^{11}$ metres, is a relative precision of $10^{-11}$ — and at a five-million-year Lyapunov time that is exhausted in about 125 million years. The best-determined initial conditions in all of dynamics buy an eighth of a per cent of the age of the system.

The one piece of genuinely observational evidence about the long-term behaviour is geological, and it is remarkable that it exists at all. Sedimentary sequences record the Milankovitch cycles — the beat frequencies of the Earth’s orbital and axial variations — and a long enough sequence dates them. The 405,000-year cycle from the $g_2 - g_5$ secular resonance between Venus and Jupiter is stable enough to be used as a clock back through the whole Mesozoic — a beat between two precession rates, read out of rock, and a relative of the eccentricity term in the Earth’s own insolation. Others are not: the 2.4-million-year cycle between Mars and the Earth is found in sediments at a period that has demonstrably changed, which is a direct geological observation of chaotic drift in a planetary orbit.

The arithmetic that made the answer possible #

Every statement in this essay is the output of an integration over five billion years, and a five-billion-year integration of the solar system was not possible with a general-purpose method. What made it possible is a choice of scheme, and the choice is worth setting out because it is the reason the numerical error can be separated from the physical divergence at all.

An ordinary integrator — a Runge–Kutta method, say — approximates the solution and accumulates energy error that grows with the number of steps, typically linearly. Over $10^{11}$ steps that error swamps everything, and no amount of care with the step size fixes it, because the growth is in the method rather than in the arithmetic.

A symplectic integrator does something different. It solves exactly a Hamiltonian system slightly different from the intended one, chosen so that its solutions stay close to the true ones. The energy error therefore does not grow: it oscillates about a small offset, indefinitely, however many steps are taken. The trajectory is wrong in phase and right in structure, which is exactly the division this essay is about.

The further refinement that made planetary integrations practical splits the Hamiltonian into a large Keplerian part and a small interaction part. Each planet’s motion about the Sun is advanced analytically, using the exact two-body solution, and only the planetary perturbations are stepped numerically. Since the perturbations are a thousandth of the Keplerian term, a step size of a fraction of the shortest orbital period gives accuracy that a direct method would need a step orders of magnitude smaller to match.

So the long integrations exist because the equations were rewritten rather than because computers got faster, and the rewriting exploits precisely the fact that the solar system is nearly integrable — the same near-integrability that makes the semi-major axes bounded. The tool and the result rest on the same property.

Where the model stops #

The Lyapunov exponent is a local quantity presented as a global one. It is defined as a limit along a single trajectory, and in a system with several distinct chaotic zones the number depends on where the trajectory spends its time. Quoted as “the Lyapunov time of the solar system” it is an average over an interval and a region, and the interval matters: measured over ten million years it comes out differently from the same calculation over a hundred.

Exponential growth is a statement about infinitesimal separations. The straight stretch in the first figure is only straight while the separation is small compared with the system; after that it saturates, because two bounded systems of the same energy can only be so far apart. Everything the exponent says is about the early part of the curve, and the interesting question — what happens after the horizon — is a question about the saturated regime, where the exponent says nothing.

Chaos does not imply instability, and the reverse is what everyone assumes. A trajectory can be exponentially sensitive and confined to a bounded region for ever. The solar system is thought to be exactly that in the outer planets, where the Lyapunov time is short but nothing ever leaves. The inner planets are the case where the boundedness is not certain — and the uncertainty is not about the chaos, which is established, but about whether the chaotic zone touches a configuration in which something crosses something else.

What the picture cannot show #

The dimension. A three-body plot is two spatial dimensions of a six-dimensional phase space; the real problem is fifty-four dimensional at minimum. Every drawing here is a shadow, and the structures that matter — resonance zones, separatrices, the web that connects them — are geometrical objects in that full space with no two-dimensional representative.

Time. The divergence figure runs for sixty time units, which is chosen so that the exponential stretch is visible. The solar-system statements are about $10^{9}$ years, which is $10^{10}$ orbits of Mercury, and no figure on any page can show that many of anything.

The probability. The modern result is a distribution over an ensemble of runs, and a plot of one run — however carefully integrated — is one sample from it. A picture of a trajectory that survives says nothing about whether the trajectory that does not is one run in a hundred or one in a million, and the whole content of the modern answer is that number.

What a probability over an ensemble means #

The result quoted below is a percentage, and a percentage attached to a single system deserves a moment’s care, because there is only one solar system and it will do whatever it does.

The probability is not a statement about the solar system’s propensity to destabilise. It is a statement about a set of solar systems, all of which are consistent with every measurement anybody has made, differing among themselves by amounts far below the observational uncertainties. One per cent of that set destabilises. Which member of the set this one is, nobody knows and nobody can find out, because the differences between them are smaller than anything measurable and grow to dominate only over the interval in question.

That is an unusual kind of ignorance and it is worth distinguishing from the ordinary kind. Normally a probability quoted for a unique event stands in for a lack of information that better measurements could remove. Here better measurements cannot remove it: the sensitivity is exponential, so reducing the initial uncertainty by any factor a laboratory could achieve moves the horizon by a few Lyapunov times and leaves the five-billion-year question exactly where it was.

The practical response has been to treat the ensemble as the object of study rather than the trajectory. What is computed and reported is the distribution of outcomes — how often Mercury’s eccentricity exceeds a threshold, how often an orbit crossing occurs, how often a collision follows — and those distributions are stable: run a different ensemble with a different integrator and the percentages agree, even though no individual trajectory does.

The stable quantity is the statistic and the unstable one is the trajectory, which is the same division the whole essay has been making, applied one level up.

The answer stopped being a yes or a no #

Laskar’s 2009 result, from 2,501 integrations over five billion years with initial conditions scattered inside the observational uncertainties, is the standing answer. In roughly one per cent of the runs Mercury’s eccentricity grows past 0.7 and it crosses Venus’s orbit; in a handful of those, a collision or an ejection follows, and in one, the Earth is involved.

That is not a prediction, and its form is the point. The system is not being said to be stable or unstable; a probability is being attached to an outcome, over an ensemble of solar systems all consistent with the one that is observed. The observations do not determine which of those solar systems this is, and no observations could, because the differences between them are far below anything measurable and grow to dominate over the interval in question.

The surprising part is how much survives that. The semi-major axes stay put in essentially all the runs. The outer planets are untouched in all of them. Ninety-nine per cent of the ensemble looks, at the end, very much like the beginning. A system can be exponentially unpredictable in every one of its phases and still be recognisably itself five billion years later — for as long, at least, as the star at the centre keeps the mass it has —, and the reason is the division the first three figures make: the phase is what diverges, and the phase is not what the system is.

One more integration extends the divergence past the point where any prediction survives.

Where the ladder goes next #

Later rungs on this anchor: the Kolmogorov–Arnold–Moser theorem, which says that most of the invariant tori survive a small perturbation and is the modern replacement for Laplace’s proof. Chirikov’s overlap criterion, quantitatively, and the resonance web it implies. The transition to global chaos as the perturbation grows. Arnold diffusion, the slow leak along resonance channels in more than two degrees of freedom, which is the mechanism nobody has yet demonstrated in a real planetary system. Chaotic obliquity, and why the Earth’s is stable only because of the Moon. And the astronomical time scale in geology, which is the one place where a chaotic dynamical prediction can be checked against a physical record millions of years long.

Poincaré found the mechanism in 1890 while trying to win a prize for proving the opposite. He had submitted a proof of stability, the error was caught in press, and what replaced it — homoclinic tangles, and a figure he said he would not attempt to draw — is the beginning of this subject. The result that founded chaos theory was a corrected mistake about the solar system, and the correction took a century of arithmetic to turn into a number.