Astronomy & Celestial Mechanics
the sky, computed
Gravitation

A probability for a system with one future

The Solar System will do one thing, and an integration of it is wrong about which thing within a hundred million years. What long integrations report instead is a probability — one per cent that Mercury's orbit goes unstable before the Sun dies — and a probability of a single event is a strange object. It is a count over an ensemble of futures that differ by less than anything measurable, which makes sense only because chaos makes those futures independent; and how many of them are needed is a binomial question with an unforgiving answer for rare events.

Assumes Numerical integration, Chaos and Patched conics.

16 min read 6 figures Measuring the unreachableSmall perturbations, large consequences

The best integrations of the Solar System over the five billion years the Sun has left report that Mercury’s orbit has about a one per cent chance of becoming eccentric enough to cross Venus’s, collide with it or with the Sun, or destabilise the inner planets. The number comes from an ensemble: 2,501 integrations of the full Solar System, each started from initial conditions differing from the others by a fraction of a millimetre in Mercury’s position, of which about twenty went unstable.

The idea that a deterministic system has a probability needs some care. There is one Solar System and it will do one thing. Its equations have no random term in them. And the integrations that produce the probability are, individually, wrong about where everything is within a hundred million years, because the system is chaotic and every integration’s errors, however carefully controlled, grow exponentially from then on. This essay asks what the ensemble is measuring, why that measurement survives the errors that destroy each member, and how many members it takes — using the simplest chaotic system that diffuses, so that every number can be checked.

Sixty futures that began a trillionth apart. The momentum of the standard map at K = 3, a kicked rotor whose phase space is almost entirely chaotic, for 60 members started 10⁻¹² apart in angle, against iteration. For 41 iterations they are one trajectory; then the separation, growing by a factor of about 2.0 each step, reaches the size of the system and they become independent random walks. Three of the 60 cross the threshold at |p| = 90 within 400 iterations (dark), the stand-in for an orbit that goes unstable. Which members cross is decided by differences far below anything that could be measured; how many cross is a property of the dynamics, and it is the only thing the ensemble can be asked for.
Fig. 1 The momentum of the standard map at K = 3 for sixty members started 10⁻¹² apart, against iteration. For 41 iterations they are one trajectory; then they become independent random walks. Three of the sixty cross the threshold at |p| = 90 within 400 iterations (dark), the stand-in for an orbit gone unstable.

A model with one dial

The standard map is the simplest system in which the questions can be asked cleanly. It describes a rotor kicked once per period by a force that depends on its angle: the momentum pp changes by Ksin⁡xK\sin x at each kick, and the angle xx advances by the new momentum. It is a Poincaré section of a real Hamiltonian system, and for kicks stronger than about K=1K = 1 its phase space is dominated by a single connected chaotic sea, in which the angle is scrambled within a few kicks and the momentum wanders like a random walk.

The wandering momentum plays the role of Mercury’s eccentricity. In the Solar System the eccentricity diffuses slowly through a web of secular resonances whose frequencies the linear theory computes badly; occasionally it reaches values where the inner planets’ orbits begin to cross, and then the system goes unstable quickly. Here the momentum diffuses, and an “instability” is its crossing a fixed threshold. The question for both is the same: what is the probability that the crossing happens within a given time?

The hero figure is one answer and the whole difficulty. Sixty members started a trillionth of a radian apart follow one trajectory for forty iterations, indistinguishable at the resolution of the plot; then they come apart and wander independently. Three of the sixty cross the threshold within four hundred iterations. Which three is decided by differences of a trillionth in the starting angle — differences no measurement could ever make. That three of them did is a property of the dynamics.

A random walk made by a deterministic rule

The momentum’s wandering is the part of the model that carries over most directly, so it is worth seeing why it is a random walk at all. At each kick the momentum changes by Ksin⁡xK\sin x, and the angle at which the kick lands is the angle after the previous step, advanced by the momentum. When the momentum is large and the kicks are strong, a small change in the momentum moves the next angle by a large amount modulo 2π2\pi, so successive angles are effectively uncorrelated and the kicks Ksin⁡xK\sin x add like independent random steps of mean zero and mean square K2/2K^2/2. The momentum’s spread then grows as the square root of the number of kicks, with a diffusion coefficient near K2/4K^2/4 per kick — corrected, for moderate kicks, by correlations between successive angles that the simple estimate ignores.

Mercury’s eccentricity behaves analogously over much longer times. It is driven by the slow oscillations of the other planets’ orbits, which act on it through combinations of secular frequencies, and when several of those combinations are close to resonance at once the drive becomes erratic and the eccentricity diffuses. The integrations show it wandering by a few hundredths over a billion years for most members, and occasionally reaching the region near 0.6 where its orbit begins to cross Venus’s. The probability being estimated is the probability of a first passage — of a random walk reaching a barrier within a time — and a first-passage probability is small and rises slowly at first, because the walk needs time to spread far enough to reach the barrier at all.

That shape is visible in the curves of the ensembles below: no member crosses for the first hundred iterations, a trickle begins, and the cumulative fraction rises roughly steadily after that. A probability quoted for one horizon says little about another: in the model, the first quarter of the horizon contributes almost nothing, and a probability for four hundred iterations cannot be divided by four to give one for a hundred. The same is true of Mercury’s one per cent over five billion years, which cannot be read as a fifth of a per cent per billion.

How long the members stay together

A separation that grows until it is the size of the system. The distance between 6 pairs of trajectories of the standard map, each pair started 10⁻¹² apart at a different point, against iteration, on a logarithmic scale. Each grows exponentially, with fluctuations, at a mean rate of 0.69 per iteration — a factor of 2.00 each step, faster than the ln(K/2) = 0.41 that the formula for very strong kicks gives — until it saturates at the size of the phase space. Twelve decades take 40 iterations: every ten-fold improvement in the initial conditions buys 3.3 more iterations of prediction. For the inner Solar System the corresponding rate is one e-fold in about five million years, and twelve decades last about a hundred and forty million.
Fig. 2 The separation of six pairs of trajectories started 10⁻¹² apart, against iteration, on a logarithmic scale. Each grows exponentially at a mean rate of 0.69 per iteration — a factor of 2 per step — until it saturates at the size of the phase space. Twelve decades take 40 iterations; each ten-fold improvement in the start buys 3.3 more.

Two members that start close together separate exponentially, at a rate set by the system’s largest Lyapunov exponent. For the standard map at this kick strength the separation doubles at every step, with fluctuations — some pairs pass through regions where the stretching is weaker and fall behind for a while — and twelve decades of separation take about forty steps. After that the pair’s separation is the size of the phase space itself, and the two members are as unrelated as two trajectories started anywhere in the chaotic sea.

The logarithm is what makes prediction hopeless and ensembles possible. Each factor of ten in the precision of the starting point buys a fixed number of extra steps — here three and a third — so a calculation started with a billion times better knowledge lasts only thirty more steps. For the inner Solar System the corresponding rate is one e-fold in about five million years, and a perturbation of a millimetre in Mercury’s position, which is a part in 101410^{14} of its orbit, grows to the size of the orbit in about a hundred and sixty million years. That is the horizon beyond which no integration can say where a planet is, and it is short compared with the five billion years being asked about.

The same logarithm makes the members independent. Members started a trillionth apart are clones for forty iterations and independent samples after that, and if the question concerns what happens over four hundred iterations, the first forty hardly matter. In the Solar System ensembles, members started fractions of a millimetre apart are independent after a hundred and fifty million years, and the five-billion-year question is asked over thirty such times.

A count over futures

The probability is then simply a count: the fraction of independent members for which the event happens. Its error is binomial, and that is where rare events become expensive.

The probability of instability, from ensembles of three sizes. The fraction of an ensemble whose momentum has crossed the threshold by each iteration, for ensembles of 20, 200, 4,000 members started 10⁻¹² apart. The largest says the event has a probability of 1.0% by iteration 400. The 20-member ensemble says 5.0% — off by a factor of 4.9 — and the 200-member one 0.5%. Each curve is a valid sample of the same probability; a small one is a sample with a large error, and with a rare event the error is larger than the answer.
Fig. 3 The fraction of an ensemble whose momentum has crossed the threshold by each iteration, for ensembles of 20, 200 and 4,000 members. The largest gives 1.0% by iteration 400. The 20-member ensemble gives 5.0%; the 200-member one 0.5%.

The three curves are three valid estimates of the same probability, and they disagree by an order of magnitude. Twenty members can report only multiples of five per cent; one of them happened to cross, and the ensemble reports five times the truth. Two hundred members report half of it. Neither is wrong as a sample — a probability of one per cent gives one crossing in twenty members about sixteen per cent of the time and one crossing in two hundred about a quarter of the time — but neither is a measurement.

Members an ensemble needs, against how rare the event is. The relative error on a probability estimated by counting, √((1 − p)/(pN)), against the number of independent members, for events of probability 1.0% and 10%. It falls as one over the square root of the number of members and rises as one over the square root of the probability, so a rare event needs proportionally more: a ten per cent error on a 1.0% event needs 9,900 members, on a 10% event 900. The marked ensemble of 2,501 — the size of the largest published set of integrations of the whole Solar System over five billion years — measures a one per cent probability to ±20 per cent of itself. The count is all binomial: the members are independent only because chaos made them so.
Fig. 4 The relative error of a counted probability, (1−p)/(pN)\sqrt{(1-p)/(pN)}, against the number of independent members, for events of 1% and 10%. A 10% error on a 1% event needs 9,900 members; on a 10% event, 900. An ensemble of 2,501 measures a 1% probability to ±20% of itself.

The formula is the textbook binomial one, and its consequence for rare events is the point. The relative error falls as one over the square root of the number of members, and rises as one over the square root of the probability itself. A ten per cent event needs nine hundred members for a ten per cent error; a one per cent event needs ten thousand; a one-in-a-thousand event needs a hundred thousand. The largest ensemble of full Solar System integrations over five billion years, 2,501 members, measures its one per cent to about plus or minus a fifth of itself — between 0.8 and 1.2 per cent. That is a real measurement and a modest one, and it took years of computing on a large cluster to make.

There is a quieter assumption in the formula: that the members are independent. A thousand members started so close together that they had not separated by the time the event could happen would be a single member counted a thousand times, and the binomial error would be a fiction. Chaos is what licenses the count. The Lyapunov time that destroys prediction is also what turns an ensemble of nearly identical starts into an ensemble of independent futures, and an ensemble of a regular system — one with no positive Lyapunov exponent — would give exactly the wrong answer with exactly the same confidence.

Clones, and the cost of starting too close

The independence the count relies on is not automatic, and an ensemble designed carelessly can lose it. Members started so close together that their separation has not reached the size of the system by the time the question concerns are not independent: they are the same future, counted many times. For the standard map that happens if the horizon is shorter than about forty iterations, or if the starting separation is so small that twelve decades become twenty — and a starting separation below the precision of the arithmetic is worse still, because members that start closer than the last bit of a double-precision number start in the same place.

The effective number of independent members in an ensemble is therefore the number whose separations have grown to the size of the system, and the designers of the Solar System ensembles chose their spacing with this in mind: 0.38 millimetres in Mercury’s initial position, far below any measurement and far above the arithmetic’s resolution, with a hundred and fifty million years to grow before the questions start. Over five billion years the members are independent for more than ninety-five per cent of the interval. For a question about the next hundred million years — whether anything goes wrong soon — the same ensemble would contain far fewer independent futures than members, and its binomial error would understate the true one.

The mirror-image mistake is spreading the members too widely, across initial conditions the observations exclude. In a single connected chaotic region this costs nothing, as the next figure shows, but a real phase space has structure, and members started outside the observed uncertainty may begin in a region with different behaviour.

A probability of what, exactly

Every member is started inside a tiny region around the measured initial conditions. The probability the ensemble reports is the fraction of that region whose futures cross the threshold. It would be natural to expect the answer to depend on how big the region is — on how well the initial conditions are known.

A probability that does not care how well the start is known. The fraction of 3,000 members crossing the threshold within 400 iterations, against the width of the range of initial angles they are spread over, from 10⁻¹⁰ radians to a radian; bars are one standard error. Over ten decades of width — from a cloud far smaller than any measurement to one a radian across — the probability stays near 1.3%, scattering by about its counting error. It does not change even for the widest cloud, because the chaotic sea is one connected region and any cloud inside it is mixed across the same region within a few dozen iterations; it would change only if the cloud reached parts of the phase space with different long-term behaviour — an island of regular motion, or a different momentum. The ensemble measures a property of the region the system is in, not of the uncertainty in where it is, which is why a probability can be computed for a system that has only one future.
Fig. 5 The fraction of 3,000 members crossing the threshold within 400 iterations, against the width of the range of starting angles, from 10⁻¹⁰ radians to a radian; bars are one standard error. Across ten decades the probability stays near 1.3%, scattering by about its counting error.

It does not. A cloud of starting angles a ten-billionth of a radian across and one a radian across give the same probability within the counting error. Inside a connected chaotic region, every cloud is stretched and folded across the whole region within a few Lyapunov times, and what it remembers of its own size after that is nothing. The probability is a property of the region the system is in — the chaotic sea, with its particular rate of diffusion — and not of the uncertainty in where within it the system started.

That is the answer to the question of what the probability means. It is not a statement about ignorance of the initial conditions, which would change as the measurements improved. It is a statement about the dynamics of the region the Solar System occupies, which no improvement in measurement would change: any initial condition consistent with the observations, measured to a millimetre or to a nanometre, leads into the same sea, and the same fraction of the sea’s trajectories reach the threshold within five billion years. The Martian obliquity, whose past can only be given as odds, is the same kind of quantity read backwards.

The independence of the width breaks down at two edges. A cloud large enough to reach a different region — an island of regular motion inside the sea, or a part of the sea with a different diffusion rate — would change the answer, because then the initial conditions would genuinely matter. And for a chaotic orbit that clings to the edge of an island, the time needed to forget the start can be much longer than the Lyapunov time suggests, and an ensemble whose members are all still stuck in the same sticky layer counts one history many times.

Arithmetic as one more perturbation

The members of an ensemble are not exact solutions of the equations. Each is computed with a finite step and finite precision, and those errors grow exponentially just like the differences in initial conditions. It would be reasonable to worry that a probability computed from members that are each wrong is itself wrong.

The same ensemble in two kinds of arithmetic. The fraction of 4,000 members crossing the threshold by each iteration, members started a hundred-thousandth of a radian apart, computed in double precision (blue) and with every step rounded to single precision (orange) — an error of about one part in ten million added every iteration, a billion times larger than double precision's own. Member by member the two calculations are different trajectories within a few tens of iterations: of the members that cross in either, only 2 cross in both, and 97.0% of all members have the same fate in both only because most never cross. The fractions agree: 1.5% and 1.5%, within the counting error of ±0.2%. Rounding acts on a chaotic system as one more perturbation, and an ensemble's statistics are indifferent to it — which is why long integrations can be trusted for probabilities and never for individual futures.
Fig. 6 The fraction of 4,000 members crossing the threshold, computed in double precision (blue) and with every step rounded to single precision (orange). Member by member the two calculations diverge within a few tens of iterations — of the members that cross in either, only two cross in both — but the fractions agree: 1.5% and 1.5%.

Rounding every step to single precision adds an error of a part in ten million at every iteration, a billion times larger than double precision’s own. Within a few tens of iterations every member of the rounded ensemble is following a different trajectory from its double-precision twin, and of the members that cross the threshold in either calculation, only two cross in both. The probabilities agree to within their counting error. In a chaotic region, a numerical error is one more small perturbation, indistinguishable in its effect from a different starting point; the rounded member is an exact trajectory of a slightly different start, and the ensemble of such members samples the same sea.

This is the deep reason long integrations of the Solar System can be trusted for statistics and never for particular futures, and it carries a condition. The errors must act like perturbations of the starting point, not like a change of the equations. A step adapted to the orbit can break the conservation laws the dynamics depends on and introduce a systematic drift — an energy that creeps, a precession that is slightly too fast — and a systematic error moves the whole sea rather than shuffling trajectories within it. That is why the integrators used for these ensembles are symplectic, with fixed steps, and why the error that grows like a random walk rather than linearly is the one they are designed to have.

The past has the same horizon

The horizon runs backwards as well as forwards. The positions of the planets fifty million years ago are as uncertain as their positions fifty million years from now, because the integration backwards is subject to the same exponential growth of errors. That matters in a field far from celestial mechanics. The layered sediments of the ocean floor record the slow cycles of the Earth’s orbit and tilt — the same cycles whose obliquity drives the ice sheets — and geologists date the layers by matching them to an astronomical solution computed backwards. Beyond about fifty million years the solutions computed from slightly different initial conditions disagree about the phase of the longest cycles, and the matching can no longer be done with one solution.

So the dating of old sediments is done the other way round: an ensemble of astronomical solutions, and the sediment record used to choose among them. The geological record then becomes a measurement of the Solar System’s past state, constraining which members of the ensemble are consistent with what actually happened. A probability over futures, read backwards, is a set of possible pasts, and the rocks select from it — which is the one place where a chaotic system’s single realisation can be checked against its ensemble at all.

What the model leaves out

The standard map is a stand-in. Its diffusion is fast and nearly uniform, the threshold is a clean line, and the time to cross it is a few hundred iterations. Mercury’s eccentricity diffuses through a web of resonances whose overlap varies from one part of the phase space to another, its instability is reached by a particular resonance with Jupiter’s perihelion precession rather than by a uniform drift, and the crossing time is billions of Lyapunov times — long enough that the eccentricity’s slow diffusion is better described by a stochastic equation fitted to the integrations than by any analytic rate.

The ensemble’s members also share their equations. Every member of a Solar System ensemble uses the same model: the same masses, the same treatment of general relativity, the same neglect of the asteroids beyond a few hundred, the same Sun losing mass at the same rate. Those are not perturbations within the sea; they decide where the sea’s edges are. General relativity, in particular, is not a small correction here: integrations without it find instabilities of Mercury about sixty times more often, because the relativistic precession moves Mercury’s perihelion frequency away from the resonance that destabilises it. A probability is conditional on the model, and the model’s errors are not counted by the binomial formula.

Still open: how many is enough for a rarer event

One per cent is a rare event, but not a very rare one. The questions that matter beyond it — whether the outer planets could ever go unstable, whether the Earth could be ejected, the probability of collisions in the next billion years rather than the next five — concern events at a part in a thousand or less, and the binomial count makes them a hundred times more expensive at fixed precision.

The way through is to stop counting and start steering. Methods developed for rare events in statistical physics run many short integrations, keep the ones moving towards the event and branch them, and reweight the results so that the probability of the rare outcome is estimated from far fewer full integrations. Applied to the Solar System, they have begun to estimate probabilities of events far rarer than any brute-force ensemble could see. Whether their reweighting is trustworthy for a system whose diffusion is not uniform, and whose rare paths may run through narrow channels the method never samples, is the question that will decide how small a probability of the Solar System’s future can be stated with confidence — and every answer will still be a count over futures, made meaningful only by the chaos that makes them independent.

Where this goes next

A step size below which a smaller step is worse. The error left in a long integration against the step size, for methods of three different orders, with both contributions drawn. The falling lines are truncation error, whose slope on these axes is exactly the order of the method. The rising line is round-off, identical for all three because it is a property of the arithmetic and not of the algorithm: every operation loses a few bits, the losses are independent, and they accumulate as the square root of the number of steps — which is why its slope is −1/2 and why it rises as the step shrinks. Each method's total has a minimum, at a step of 1.0e-6, 5.0e-6, 9.7e-4 for orders 1, 2, 4. Below that minimum every halving of the step costs time and makes the answer worse. That is the practical reason a solar-system integration is not run at an arbitrarily fine step, and it is a reason with nothing to do with computer time. An error that grows like a random walk A long integration accumulates two errors with opposite habits. One falls when the step is made smaller and grows in proportion to the time; the other grows when the step is made smaller and accumulates as a square root. Which of the two dominates decides whether a billion-year integration means anything. A separation of 10⁻¹⁰ becomes 0.0285. The distance between two copies of Burrau's problem, started 10⁻¹⁰ apart in one coordinate of one body and then stepped in lockstep — one loop, one step size, both states advanced by it — for 60 time units. The vertical axis is logarithmic, so the straight stretch is exponential growth, and its slope is the Lyapunov exponent: 0.3413 per time unit, fitted over the 604 samples lying between ten times the starting offset and a tenth of the system's own size. That is an e-folding every 2.93 time units, so the separation multiplies by ten every 6.75 and by 7.8·10⁸ across the whole run. It has not saturated within the drawn interval: the growth rate over the last tenth of the run is still 0.6002 per time unit, 176% of the fitted exponent, and the separation has reached only 0.90% of the system's own size. The exponent is what sets a prediction horizon, and this run is drawn short of it on purpose: the straight line is the measurement, and the flat part that follows is arithmetic about how far apart two bounded systems can get. A prediction with an expiry date The inner solar system's Lyapunov time is about five million years, so a centimetre of error becomes an orbit in a hundred million. The ephemeris dies while the system survives — because the elements stay bounded when the phase does not, and only one of those is what stability means. The divisor is 0.129″/yr and the theory's own error is 0.24. Six frequencies of the secular solar system on one logarithmic axis, in arcseconds per year. The top two are the pair whose near-equality is the whole story: the perihelia of Mercury and Jupiter separate at 1.333″ a year, the nodes of Mercury and Venus at 1.462, and the difference of those two differences is 0.129 — a resonant argument that turns once every 10.0 million years. A term with that argument in the disturbing function acts in one direction for five million years at a stretch, which is what pumps Mercury's eccentricity, and it is the reason the inner solar system's Lyapunov time is what it is. The bottom three bars are why this figure exists. The divisor is smaller than the corrections the theory that computes it leaves out. Relativity contributes 0.4298″ a year to g₁ alone — the same 43 arcseconds a century that broke Newtonian gravity — which is 3.3 times the divisor; the fourth-order terms in the eccentricity that Laplace–Lagrange truncates come to about 0.24″, which is 1.8 times it; and the second-order solution computed on this page gets 0.35″, missing the published value by more than the value itself. A theory cannot bound what it cannot resolve. Laplace's proof that the eccentricities stay bounded is a proof about a system whose frequencies are constants, and the frequency that decides the question is not one. The bound that holds only in the linear theory Laplace proved the planetary eccentricities bounded, and the proof is a proof about a linearised system with constant frequencies. One combination of those frequencies is nearly zero — and it is smaller than the terms the linearisation threw away, which is why the stability of the solar system is a probability rather than a theorem.

About the same objects

Not linked from either essay — found by the objects both name.

What links here

Essays that link to this one from their own argument.

The objects this essay names

Each one links to every other essay that touches it.

Binomial statisticsChaosEnsembleLyapunov exponentNumerical integrationRoundoff errorSolar system stabilityStandard map