Astronomy & Celestial Mechanics
the sky, computed
Gravitation

A rare event reached by cloning the near misses

Counting how often an ensemble of integrations goes unstable works for an event of one in a hundred and becomes hopeless for one in a million, because a ten per cent answer needs a hundred events and so a hundred million futures. Splitting gets there by a different route: wait for the members that come closest, copy them, and run the copies on. Each stage measures a common event and the rare one is a product of common ones — provided the copies forget their parents and the stages know how much time is left.

Assumes Numerical integration and Chaos.

16 min read 5 figures Small perturbations, large consequencesScales that defeat intuition

The question whether the Solar System will stay as it is has no single answer, because its future is chaotic: two integrations started a fraction of a millimetre apart diverge completely within a hundred million years. What it does have is a probability, and the probability is measured by running many integrations from slightly different starts and counting how many end in disaster. The largest such ensemble, 2,501 integrations of the planets over five billion years, found that Mercury’s orbit becomes eccentric enough to collide with Venus or the Sun in about one per cent of them.

Counting works for that, barely. The relative error of a probability estimated by counting is (1−p)/(pN)\sqrt{(1-p)/(pN)}, so a one per cent event counted among 2,501 members is known to about twenty per cent of itself. But the questions that matter more are rarer. Could the outer planets go unstable? Could the Earth be ejected, or collide with Mars? Could anything happen in the next billion years rather than the next five? These are events at a part in a thousand, a million, or less, and counting them to ten per cent needs a hundred events — a hundred thousand integrations for one in a thousand, a hundred million for one in a million. At the cost of a five-billion-year planetary integration, that is not a computing problem; it is a different century.

The way out is to stop waiting for the rare event and start steering toward it. The hero figure shows how.

Splitting, in one picture

The system is a kicked rotor — the standard map, p′=p+Ksin⁡xp' = p + K\sin x, x′=x+p′x' = x + p' — at a kick strength K=10K = 10, strong enough that almost all of its phase space is chaotic and the momentum pp performs a random walk. It is the same model that stands in for the slow chaotic diffusion of Mercury’s eccentricity, with the event being ∣p∣|p| passing a threshold within a fixed horizon of four hundred iterations — an orbit going unstable before the age of the Sun runs out.

Multilevel splitting divides the distance to the event into stages. A handful of members are run from the start until each either reaches the first level or runs out of time. The ones that reach it are copied, back up to the original number, and each copy is restarted from its parent’s crossing point with a tiny nudge to its angle. Chaos does the rest: within a few iterations the nudged copies have independent futures, and they are run on to the next level. The fraction that passes each level is a probability of tens of per cent, which a few dozen members measure well. The rare event’s probability is the product of the fractions.

In the hero figure sixteen members start near zero momentum. Two reach the first level; their copies spread, and one lineage carries on through the next three, ending at the threshold of 450 just before the horizon. The fractions passing the four levels are 0.13, 0.06, 0.63 and 0.88, and their product, 4×10−34\times10^{-3}, is the estimate. That small run is a rough estimate, but it reached an event that counting with sixteen members would essentially never see.

Where counting stops and splitting does not

The method is unbiased when done correctly: averaged over many runs, the product of fractions equals the true probability. That claim is worth checking rather than trusting, and the check is to compare it with counting where counting can still be done.

Counting stops where steering keeps going. The probability that the kicked rotor's momentum passes a threshold within 400 iterations, against the threshold, estimated two ways. Counting among 100,000 independent members (4·10⁷ iterations in all) finds 0.052 at 250, 0.016 at 300, 0.0036 at 350, 6.8·10⁻⁴ at 400, 6·10⁻⁵ at 450, 2·10⁻⁵ at 500, nothing at 550, nothing at 600. Splitting with 300 members per level, at an average cost of 3·10⁵ iterations per threshold — 132 times less — agrees with the count within the scatter of a single run — a few tens of per cent — wherever the count has thirty or more events, and carries on past where it has none, to 1.3·10⁻⁷ at 600. Dashed rings: counts of zero, drawn at the floor of one event in 100,000.
Fig. 1 The probability that the rotor’s momentum passes a threshold within 400 iterations, estimated by counting among 100,000 independent members (4×1074\times10^7 iterations, rings) and by splitting with 300 members per level (3×1053\times10^5 iterations per threshold, dots). The two agree within the scatter of a single splitting run wherever the count has thirty or more events. The count finds nothing above 500; splitting continues to 1.3×10−71.3\times10^{-7} at 600. Dashed rings mark counts of zero at the floor of one event in 100,000.

The figure carries two results. Where both methods work — thresholds up to about 400, probabilities above a part in a thousand — splitting agrees with counting to within its own run-to-run scatter, at a hundred and thirty times less computation per threshold. And where counting fails, splitting goes on: at a threshold of 600 the probability is about one in eight million, an event that a hundred thousand members would see with probability one in eighty, and that splitting estimates from a few hundred thousand iterations of the map.

The agreement in the overlap is what licenses the extrapolation. A splitting estimate cannot be checked against counting at one in ten million, by definition; it can be checked at one in a thousand, using the same code and the same levels, and its unbiasedness there is evidence it is unbiased further out. That is the same logic every rare-event method in physics relies on — validate where the brute-force answer exists, then trust the method where it does not.

A cost that grows with the logarithm

The cost comparison is the reason the method exists.

Iterations a ten per cent answer costs, two ways. The number of iterations needed to estimate the probability of crossing a threshold to ten per cent of itself, against that probability, both logarithmic. Counting needs about 100/p independent members of 400 iterations each, so its cost rises as 1/p (rings). Splitting's cost is measured from 12 independent repetitions at each threshold and scaled to a ten per cent spread (dots): 6.2·10⁵ at p = 0.016 and 8·10⁶ at p = 2.2·10⁻⁷, where counting would need 2.5·10⁶ and 1.8·10¹¹. The advantage is 4-fold for the commonest event drawn and 22432-fold for the rarest; it keeps growing because each extra decade of rarity adds a level or two to splitting and a factor of ten to counting.
Fig. 2 Iterations needed for a ten per cent estimate, against the probability of the event, both logarithmic. Counting needs about 100/p100/p members of 400 iterations each (rings). Splitting’s cost, measured from twelve independent repetitions at each threshold and scaled to a ten per cent spread (dots), is 6.2×1056.2\times10^5 at p=0.016p = 0.016 and 8×1068\times10^6 at p=2.2×10−7p = 2.2\times10^{-7}, where counting would need 2.5×1062.5\times10^6 and 1.8×10111.8\times10^{11}. The advantage is four-fold for the commonest event and over twenty-thousand-fold for the rarest.

Counting’s cost is proportional to 1/p1/p: every factor of ten in rarity is a factor of ten in members. Splitting’s grows far more slowly. Each extra decade of rarity adds a level or two, and each level costs about the same, so in the ideal case the cost grows with the logarithm of 1/p1/p — in practice with its square, because the errors of successive stages compound. Over the five decades drawn, counting’s cost rises by a factor of seventy thousand and splitting’s by a factor of thirteen.

For the commonest event in the figure, the methods are within a factor of four of each other, and counting is simpler and has no tuning to get wrong. The advantage appears only where it is needed, which is the sense in which splitting is not a better way of doing what counting does but a way of asking questions counting cannot.

The same budget, spent two ways

A comparison at equal cost shows what the advantage means for a single study.

The same budget, spent two ways. 40 independent splitting estimates of the probability that |p| passes 500 within 400 iterations, each with 300 members per level and about 3.5·10⁵ iterations, plotted as the ratio of each estimate to their mean on a logarithmic scale (bars), against what counting would report for the same budget: 885 members, a Poisson count with a mean of 0.02 events (stems; the column at the left is a count of zero). The splitting estimates scatter about their mean of 1.9·10⁻⁵ by 42 per cent. The count reports zero 98 per cent of the time and, when it reports anything, at least 59 times the truth. The budget is the same; the information is not.
Fig. 3 Forty independent splitting estimates of the probability that ∣p∣|p| passes 500 within 400 iterations, each costing about 3.5×1053.5\times10^5 iterations, as the ratio of each estimate to their mean on a logarithmic scale (bars), against what counting would report with the same budget: 885 members, a Poisson count with a mean of 0.02 events (stems, with the column at the left for a count of zero). The splitting estimates scatter by 42 per cent about their mean of 1.9×10−51.9\times10^{-5}. The count reports zero 98 per cent of the time and, when it reports anything, at least 59 times the truth.

With the same computing, counting almost always reports that the event never happens, and on the rare occasions it sees one, it reports a probability nearly sixty times too high, because one event among 885 is the smallest non-zero answer it can give. Neither outcome is wrong in the statistical sense — averaged over many studies the count is unbiased too — but no single study learns anything. The splitting runs scatter by forty per cent around the right answer, every one of them. A forty per cent estimate of a one-in-fifty-thousand probability is a useful number; a zero is not.

The clock is part of the distance

There is a subtlety in how the levels are placed, and it turned out to be the whole of the method’s skill. The obvious choice is to put the levels at fixed values of ∣p∣|p|: 60, 120, 180, and so on up to the threshold. A member that reaches ∣p∣=300|p| = 300 is closer to 500 than one at 200, so it is cloned.

But the event has a horizon. A member that reaches 300 at iteration 390 has ten iterations left and almost no chance of going further, while one that reaches 250 at iteration 50 has plenty. Fixed levels clone both the same way. The late arrivals produce copies that fail at the next level, the early ones produce copies that succeed, and the estimate comes to depend on which few lineages happened to arrive early — a high-variance gamble, unbiased on average and erratic in any one run.

A rare crossing reached in stages. Multilevel splitting on a kicked rotor (the standard map at K = 10), aimed at the event |p| > 450 within 400 iterations. 16 members are run until each reaches the first level or runs out of time; those that reach it are copied back up to 16, each copy restarted from its parent's crossing point with its angle nudged by a thousandth of a radian, and run on toward the next. The levels (dashed) are fixed values of |p|, the naive choice: a member that reaches one late, with little time left, is cloned exactly like one that reached it early. 5 levels; the fractions passing them are 0.88, 0.38, 0.25, 0.19, 0.00, and their product, zero in this small run, is the estimate. Each stage is a common event; the rare one is never waited for.
Fig. 4 The same sixteen members aimed at the same threshold, with levels at fixed momenta of 90, 180, 270 and 360 (dashed). The first level is easy and the members that reach it are cloned regardless of when they got there; the late arrivals’ copies run out of time at the next level, the fractions passing fall to 0.88, 0.38, 0.25, 0.19 and then zero, and the run ends with no estimate at all.

The run in the figure fails for exactly this reason. Its early levels are crossed by many members, most of them late in the horizon, and each stage’s survivors are progressively older; by the fourth level the few that remain have too little time to reach the fifth. With more members such a run would reach the threshold eventually, and on average it would give the right answer, but the answer would rest on a few lineages that happened to be early.

The fix is to measure progress toward the event in units that include the time left. The hero figure’s levels are not straight lines but curves that rise as the horizon approaches: a member at momentum ∣p∣|p| with τ\tau iterations remaining is scored by (∣p∣−P)/2Dτ(|p|-P)/\sqrt{2D\tau}, roughly how many standard deviations of random walk it still needs, with DD the walk’s diffusion rate. Reaching a level now means being likely to finish, not merely being far along.

Levels that know about the clock. The scatter of single splitting estimates — the standard deviation of 12 independent runs of 300 members, divided by their mean — against the threshold, with levels set on |p| alone (rings) and with levels that also count the time left (dots). Both are unbiased; the naive levels clone members that reach a level with too little time left to go further, so most of their copies fail at the next stage and the estimate rests on a few lucky lines of descent. At 600 the naive scatter is 118 per cent of the answer and the clock-aware one 44. The choice of what counts as progress toward the event is the whole of the method's skill.
Fig. 5 The scatter of a single splitting estimate, as a share of the answer, against the threshold, with levels set on ∣p∣|p| alone (rings) and with levels that count the time left (dots). Both are unbiased. At thresholds up to 400 they behave alike; at 600 the naive levels scatter by 118 per cent of the answer and the clock-aware ones by 44.

The difference is small for moderately rare events and decisive for very rare ones, where most of the stages are near the end of the horizon and the naive levels fill up with members that have no time left. Since the cost to reach ten per cent goes as the square of the scatter, the clock-aware levels are about seven times cheaper at the rarest threshold.

The function that decides what counts as progress is called the importance function, and choosing it is where the physics enters a method that is otherwise pure statistics. A perfect importance function would be the probability of success itself, and with it splitting would have zero variance; a poor one makes splitting no better than counting. The art is in finding a quantity that tracks the probability of success well enough, using what is known about the dynamics.

A trick older than the problem

The idea is older than computers that could use it. It was devised at Los Alamos in the late 1940s for the transport of neutrons through shielding, where the question was how many neutrons penetrate a thick wall — a rare event, because each neutron is scattered or absorbed many times on its way. Herman Kahn and Theodore Harris described the scheme in 1951 as “splitting and Russian roulette”: a neutron that reaches a deeper layer of the shield is split into several, each carrying a fraction of its weight, and one that wanders back toward the source is subjected to a game of chance that removes most such neutrons while promoting the survivors’ weight. The estimate is unbiased because the weights keep track of what was duplicated and what was removed.

The planetary version differs in one respect that matters, and the neutron version has no counterpart to it. A neutron’s future is random; split it, and the copies diverge because the scatterings they meet are drawn independently. A planetary system’s future is deterministic, and a copy of an orbit is the same orbit. The copies diverge only because the system is chaotic and they are nudged — so the chaos that makes the long-term future unpredictable is also what makes the method work. In a regular system splitting would produce identical copies and gain nothing, and in a regular system it would not be needed.

There is a related point about the integrations themselves. The members are integrated with symplectic methods because those keep the energy error bounded over billions of years, and with a fixed step because an adaptive step breaks the symplectic property. Neither choice is affected by splitting: a clone is integrated exactly as its parent was, from a state the integrator produced, so whatever errors the integrator makes are made consistently across the whole tree.

Two conditions the copies have to meet

The method’s unbiasedness rests on two conditions, and both are physical rather than statistical.

The copies must forget their parent. A clone restarted with a tiny nudge shares its parent’s future until the nudge has grown to the size of the system — the time set by the Lyapunov exponent. If the next level is reached before that, the copies all reach it together or all fail together, and the stage’s fraction measures one future many times. In the rotor the nudge is a thousandth of a radian, and at K=10K = 10 the separation reaches the size of the system in a few iterations, quickly compared with the time between levels. With a nudge of 10−1210^{-12} the same method, on the same map, gives estimates whose scatter is about twice as large at every threshold, because the copies stay together for seventeen iterations — long enough that sibling copies often cross the next level, or fail to, as a group — which is the same failure that an ensemble started too close together suffers, appearing inside the method.

The state must carry the whole future. A clone is restarted from its parent’s crossing point, which is assumed to contain everything that decides what happens next. For the rotor that is true: the position, the momentum and the time are the whole state. For a planetary system integrated with a symplectic method it is also true — positions, velocities, the time — but for a model with hidden variables, such as a planet whose spin or internal state is not part of the integration, cloning from the orbit alone would be wrong.

What it means for the planets

The Solar System’s ensembles have not yet been run with splitting at scale, and the reasons are instructive. A planetary integration over five billion years is a hundred billion steps; the levels have to be placed on a quantity that predicts instability, such as Mercury’s maximum eccentricity over the last few million years; the clones have to be nudged by amounts small enough to stay within the observational uncertainty and large enough to decorrelate in a fraction of the time between levels; and the cost of a wrong importance function, which the rotor only hints at, is measured in years of computing.

Splitting and its relatives are routine in chemistry, where they find the rare molecular rearrangements that take microseconds in simulations that step in femtoseconds, and in climate science, where they estimate the odds of heatwaves rarer than any in the record. Applying them to planetary dynamics is recent, and has begun on simplified models of the secular motion, where Mercury’s eccentricity is driven by a handful of slowly drifting frequencies and each integration is cheap. The step to full integrations is the step from the rotor here to the planets, and it is a step in cost and care rather than in principle. When it is taken, the questions that are now answered with “no instability was seen in 2,501 integrations” — the outer planets, the Earth, the next billion years — will be answered with numbers.

There is also a payoff beyond the number. A splitting run does not only estimate how likely an instability is; it produces a population of trajectories that actually reach it, sampled in proportion to how they occur. Counting delivers one or two such histories at best. Splitting delivers hundreds, and they show the route — which resonance Mercury’s eccentricity climbed through, how long it lingered near each level, whether the route is always the same or several routes contribute. For a question about a rare event, the typical way it happens is often worth as much as how often it happens, and it is visible only in a sample of events large enough to have a typical member.

What the rotor leaves out

The rotor is a stand-in. Its momentum diffuses like a random walk with a well-defined rate, which is what makes the importance function easy to write down. Mercury’s eccentricity diffuses irregularly, through resonances that open and close as the planets’ frequencies drift, and the best importance function for it is not known.

The horizon is sharp. The event must happen within four hundred iterations. For the planets the horizon is the Sun’s remaining main-sequence life, which is also sharp, but a question like “before the Sun’s luminosity has risen by ten per cent” has a softer edge.

The levels are fixed in advance. Adaptive versions of the method place each level where a fixed fraction of the members fail, which removes one choice but not the choice of what is being measured. The comparison of naive and clock-aware levels holds for both.

The scatter is measured, not derived. Each point in the cost figure is from twelve repetitions, so the scatters are themselves uncertain by about twenty per cent. The trends — logarithmic growth for splitting, inverse for counting, and the widening gap between naive and clock-aware levels — are much larger than that.

Still open: an importance function for a planetary system

The method converts a computing problem into a modelling problem: what quantity, measured now, best predicts that a planetary system will go unstable later? For the rotor the answer is the momentum and the time left. For the Solar System candidates include Mercury’s eccentricity, the angular-momentum deficit that the giant planets hold and the inner ones might borrow, and the distance of the secular frequencies from resonance. Each would give an unbiased estimate; each would give a different variance; and which gives the smallest is a question about how instability develops that is itself not settled. The probability of a rare catastrophe is reachable. What reaching it requires is knowing, in advance, what the approach to the catastrophe looks like.

Where this goes next

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. 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. 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 last curve across the cylinder, before and after it breaks. Surfaces of section for the standard map at K = 0.5, 0.9716, 1.4, each 22 trajectories iterated 41 times from a column of starting values. Nothing here is placed: every dot is an iterate. At K = 0.5 the islands are separate and the space between them is filled with curves that run all the way round in θ — a trajectory cannot get from one island to the next, and one started beside the hyperbolic fixed point wanders 2.91 in p — 0.46 of a cylinder — and stops. At K = 0.9716, which is Greene's threshold to four figures, the last of those curves is on the point of going and the same trajectory still only reaches 4.84, which is 0.77 of a cylinder. At K = 1.4 it covers 6.1 cylinders: the barrier is gone and there is nothing left to stop it. That transition is what a chaotic zone in the asteroid belt is, drawn without any asteroids. Where the chaos comes from A Lyapunov time says how long a prediction lasts. It does not say what destroyed it. The mechanism is two resonances whose libration widths overlap, and the transition can be watched happening on a surface of section as one number is turned up.

About the same objects

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

The objects this essay names

Each one links to every other essay that touches it.

Binomial statisticsChaosEnsembleLyapunov exponentNumerical integrationRare event samplingSolar system stabilityStandard map