Direct N-body ensembles establish that Mercury's orbit has a percent-level probability of destabilizing within the next 5 Gyr through the g1−g5g_1-g_5 secular resonance, but the statistics are unconverged at early times: in the roughly 9,601 direct simulations available, too few instability events occur before about 1.3 Gyr to estimate probabilities, and the reported probabilities shift substantially with modeling choices such as integrator splitting, timestep, relativistic precession, and the solar quadrupole moment. Rare-event sampling with a hand-built reaction coordinate (Quantile Diffusion Monte Carlo using the Mercury-Venus Minimum Orbit Intersection Distance range) and a machine-learned convolutional neural network classifier can both predict future instability gigayears in advance, but the physical origin of the predictive signal is unidentified in both cases, the learned classifier has not been deployed in an operational splitting-and-killing pipeline, and ensemble-convergence criteria are only partially established. The open problem is to deliver converged, uncertainty-quantified, physically interpretable long-term instability statistics for the Solar System: probabilities of Mercury destabilization (and of the rarer subsequent destabilization of all terrestrial planets) as a function of future time over 0-5 Gyr, obtained by validated rare-event sampling and/or reduced secular or stochastic models, with an identified dynamical mechanism for the predictive signal and quantified robustness to model and integrator choices.

Background

Whether the Solar System is dynamically stable over the remaining main-sequence lifetime of the Sun (about 5 Gyr) is answered only probabilistically, because planetary motion is chaotic: the Lyapunov time, the e-folding time over which initially close orbits diverge, is only about 5 Myr for the inner planets, so individual trajectories lose predictability while statistical ensembles remain meaningful. In secular dynamics, the equations of motion are averaged over the orbital periods, and each orbit is described by slow frequencies gig_i (perihelion precession) and sis_i (node precession). Mercury's perihelion precession rate g1g_1 chaotically diffuses toward Jupiter's nearly constant g5g_5; entry into this g1−g5g_1-g_5 secular resonance transfers angular momentum between planets, can drive Mercury's eccentricity toward unity, and leads to close encounters with Venus followed, typically, by collision with Venus or the Sun (Laskar 1994; Laskar 2008). Laskar (2008) estimated from secular ensembles that the probability of such a Mercury destabilization within 5 Gyr is of order 10−210^{-2}. Direct N-body integrations, which solve Newton's equations for the Sun and planets without averaging, confirmed this at the percent level: Laskar and Gastineau (2009) integrated 2,501 solutions for 5 Gyr and found that about 1% develop a large Mercury eccentricity increase (probability 8.0±1.8×10−38.0\pm1.8\times10^{-3}), and in one solution the subsequent angular-momentum exchange destabilizes all four terrestrial planets about 3.34 Gyr from now, with possible collisions of Mercury, Mars, or Venus with the Earth. Zeebe (2015) added 1,600 direct simulations but concluded that the ensemble results were statistically inconclusive.

The central obstacle is rarity. In the combined 9,601 direct N-body simulations now available (Laskar & Gastineau 2009; Zeebe 2015; Abbot et al. 2023), only a handful of Mercury instabilities occur before 3 Gyr (2 out of 2,501 and 1 out of 1,600), and the logarithm of the instability frequency is linear in time only between 1.3 and 5 Gyr; below about 1.3 Gyr, where event probabilities fall toward 10−510^{-5}, the direct-ensemble statistics are unconverged, and simply running vastly more simulations is computationally prohibitive. Rare-event sampling addresses exactly this regime. In splitting-and-killing schemes such as Quantile Diffusion Monte Carlo (QDMC), an ensemble of simulations is periodically ranked by a reaction coordinate, a function of the system state approximating the probability of reaching the rare event (ideally the committor); simulations that progress well are cloned (split), ones that regress are terminated (killed), and bookkeeping of this genealogy keeps the resulting probability estimates unbiased. Abbot et al. (2021) applied QDMC to the REBOUND N-body code using the range of the Mercury-Venus Minimum Orbit Intersection Distance (MOID, the minimum possible separation of the two orbits) over a 400 Myr window as reaction coordinate, obtaining unbiased estimates (∼10−4\sim10^{-4} for an instability within 2 Gyr) at up to 100 times lower cost, while stating they have no complete explanation for why this coordinate predicts instability. A convolutional neural network classifier trained on 5-Gyr ensembles subsequently ranked the unstable member of a 50-member pool correctly up to 3.5 Gyr before the event (Abbot et al. 2024, Research Notes of the AAS), suggesting machine-learned reaction coordinates; however, standard explainability tools (SHAP, LIME) could not identify the physical signal being learned, and the classifier has not been run inside a splitting pipeline.

In parallel, reduced models pursue physical interpretability. A compact forced secular model of the inner planets reproduces Mercury's eccentricity statistics, the g1−g5g_1-g_5 destabilization, and a Mercury dynamical half-life of order 40 Gyr (Mogavero & Laskar 2021); computer-algebra analyses identify the resonance web responsible for the chaos (Mogavero & Laskar 2022), and truncation experiments show that degree-4 terms in eccentricity and inclination leave the model practically stable over 5 Gyr while degree-6 terms drive destabilization (Hoang, Mogavero & Laskar 2022); the full Lyapunov spectrum reveals quasi-integrals of motion that constrain the chaotic diffusion (Mogavero, Hoang & Laskar 2023). Phenomenological stochastic models treat the wander of g1g_1 as anomalous diffusion: plain diffusion underpredicts sub-5-Gyr instability probabilities by factors of 3-10,000, and a subdiffusive model (diffusion slower than Brownian) partially repairs this (Abbot et al. 2024, ApJ). Instanton theory, which predicts the most-likely path to a rare event, demonstrates path concentration in a simplified deterministic model of Mercury's resonance exit (Woillez & Bouchet 2020), with the extension to the real Solar System unresolved. Finally, the estimated probabilities are model-relative rather than absolute: numerical stepsize errors alter Mercury's possible orbits and instability probability (Hernandez, Zeebe & Hadden 2021); the choice of symplectic integrator splitting matters, since democratic heliocentric coordinates at multi-day timesteps introduce artificial precession that suppresses Mercury instabilities (Rein, Dey & Tamayo 2026; Hernandez 2026); general-relativistic precession reduces the 5-Gyr destabilization probability by a factor of about 60 (Brown & Rein 2023); the solar quadrupole moment and asteroid treatment trade off against each other in long-term integrations (Zeebe & Kocken 2025); and stellar flybys add external perturbations (Brown & Rein 2022; Kaib & Raymond 2025).

Problem Statement

Can the Solar System's long-term instability statistics be converged with validated, physically interpretable methods? The single overarching objective is a converged, uncertainty-quantified account of the probability that the Solar System's inner planets destabilize as a function of future time over the next 5 Gyr. The currently open, complementary components of this objective, posed separately by the primary literature, are:

  1. Converged statistics at early times. Provide estimates of the Mercury instability probability (Mercury's orbit destabilizing through the g1−g5g_1-g_5 secular resonance into a close encounter with Venus, collision with Venus or the Sun, or comparable disruption) as a function of future time over 0-5 Gyr, with quantified confidence intervals extending below the roughly 1.3 Gyr horizon where the existing combined direct N-body ensembles (about 9,601 members) contain too few events, together with explicit criteria distinguishing statistically converged, partially converged, and unconverged ensemble regimes for the statistics reported. The rarer subsequent full terrestrial-planet destabilization (possible collisions of Mercury, Mars, or Venus with the Earth) is the natural second statistic; because it is attested only once in the combined direct ensembles, a converged estimate for it is a desirable extension rather than a required deliverable.
  2. Validated accelerated estimation. Demonstrate accelerated estimators - splitting-and-killing rare-event schemes with hand-built or machine-learned reaction coordinates, and/or reduced secular or stochastic models - that are unbiased relative to direct simulation wherever direct statistics exist. Measured computational speed-up in an operational end-to-end pipeline (as achieved by Abbot et al. 2021 with the hand-built MOID-range coordinate) is a desirable but not required demonstration; validation against direct or rare-event statistics, including retrospective validation of a learned coordinate or validation of a reduced model, remains an admissible contribution.
  3. Physical interpretability. Identify the dynamical origin of the predictive signal in any learned or empirical reaction coordinate - which physical quantities, resonance structure, or diffusion properties it tracks - rather than leaving an unexplained correlation; connect this to the reduced-model description (secular resonances, quasi-integrals, anomalous diffusion of g1g_1, or instanton paths).
  4. Robustness to modeling choices. Quantify how the resulting probabilities depend on the dynamical model and integrator - timestep, Hamiltonian splitting, relativistic precession, solar quadrupole moment, asteroid treatment, and external perturbations - which are known to shift the 5-Gyr destabilization probability at the factor-of-60 level; the instability statistics are correspondingly tied to the model class and event definition that produced them.

These four components are individually open, jointly supporting research directions rather than a conjunctive acceptance rule: a contribution addressing any one of them advances the overarching objective.

The instability probabilities are model-relative: each refers to a particular model class, event definition, and time window, and is comparable with the published direct-ensemble event counts only within a matching model class and event definition.

References

  1. Laskar, J., & Gastineau, M. (2009). Existence of collisional trajectories of Mercury, Mars and Venus with the Earth. Nature 459, 817-819. DOI: 10.1038/nature08096. https://doi.org/10.1038/nature08096
  2. Laskar, J. (2008). Chaotic diffusion in the Solar System. Icarus 196, 1-15. DOI: 10.1016/j.icarus.2008.02.017. https://doi.org/10.1016/j.icarus.2008.02.017
  3. Laskar, J. (1994). Large-scale chaos in the solar system. Astronomy & Astrophysics 287, L9. NASA ADS bibcode 1994A&A...287L...9L. https://articles.adsabs.harvard.edu/pdf/1994A%26A...287L...9L
  4. Zeebe, R. E. (2015). Dynamic stability of the Solar System: statistically inconclusive results from ensemble integrations. The Astrophysical Journal 798, 8. DOI: 10.1088/0004-637X/798/1/8. https://doi.org/10.1088/0004-637X/798/1/8
  5. Webber, R. J., Plotkin, D. A., O'Neill, M. E., Abbot, D. S., & Weare, J. (2019). Practical rare event sampling for extreme mesoscale weather. Chaos 29, 053109. DOI: 10.1063/1.5081461. https://doi.org/10.1063/1.5081461
  6. Abbot, D. S., Webber, R. J., Hadden, S., Seligman, D., & Weare, J. (2021). Rare Event Sampling Improves Mercury Instability Statistics. The Astrophysical Journal 923, 236. DOI: 10.3847/1538-4357/ac2fa8. arXiv:2106.09091. https://doi.org/10.3847/1538-4357/ac2fa8
  7. Abbot, D. S., Hernandez, D. M., Hadden, S., Webber, R. J., Afentakis, G. P., & Weare, J. (2023). Simple Physics and Integrators Accurately Reproduce Mercury Instability Statistics. The Astrophysical Journal 944, 190. DOI: 10.3847/1538-4357/acb6ff. arXiv:2212.14844. https://doi.org/10.3847/1538-4357/acb6ff
  8. Abbot, D. S., Webber, R. J., Hernandez, D. M., Hadden, S., & Weare, J. (2024). Mercury's Chaotic Secular Evolution as a Subdiffusive Process. The Astrophysical Journal 967, 121. DOI: 10.3847/1538-4357/ad3e5f. arXiv:2306.11870. https://doi.org/10.3847/1538-4357/ad3e5f
  9. Abbot, D. S., Laurence-Chasen, J. D., Webber, R. J., Hernandez, D. M., & Weare, J. (2024). AI can identify Solar System instability billions of years in advance. Research Notes of the AAS 8, 3. DOI: 10.3847/2515-5172/ad18a6. arXiv:2401.12180. https://doi.org/10.3847/2515-5172/ad18a6
  10. Woillez, E., & Bouchet, F. (2020). Instantons for the destabilization of the inner Solar System. Physical Review Letters 125, 021101. DOI: 10.1103/PhysRevLett.125.021101. arXiv:1910.04005. https://doi.org/10.1103/PhysRevLett.125.021101
  11. Mogavero, F., & Laskar, J. (2021). Long-term dynamics of the inner planets in the Solar System. Astronomy & Astrophysics 655, A1. DOI: 10.1051/0004-6361/202141007. arXiv:2105.14976. https://doi.org/10.1051/0004-6361/202141007
  12. Mogavero, F., & Laskar, J. (2022). The origin of chaos in the Solar System through computer algebra. Astronomy & Astrophysics 662, L3. DOI: 10.1051/0004-6361/202243327. arXiv:2205.03298. https://doi.org/10.1051/0004-6361/202243327
  13. Hoang, N. H., Mogavero, F., & Laskar, J. (2022). Long-term instability of the inner Solar System: numerical experiments. Monthly Notices of the Royal Astronomical Society 514, 1342-1350. DOI: 10.1093/mnras/stac1299. arXiv:2205.04170. https://doi.org/10.1093/mnras/stac1299
  14. Mogavero, F., Hoang, N. H., & Laskar, J. (2023). Timescales of Chaos in the Inner Solar System: Lyapunov Spectrum and Quasi-integrals of Motion. Physical Review X 13, 021018. DOI: 10.1103/PhysRevX.13.021018. arXiv:2305.01683. https://doi.org/10.1103/PhysRevX.13.021018
  15. Hernandez, D. M., Zeebe, R. E., & Hadden, S. (2021). Stepsize errors in the N-body problem: discerning Mercury's true possible long-term orbits. Monthly Notices of the Royal Astronomical Society 510, 4302-4307. DOI: 10.1093/mnras/stab3664. arXiv:2111.08835. https://doi.org/10.1093/mnras/stab3664
  16. Brown, G., & Rein, H. (2022). On the long-term stability of the Solar system in the presence of weak perturbations from stellar flybys. Monthly Notices of the Royal Astronomical Society 515, 5942-5950. DOI: 10.1093/mnras/stac1763. https://doi.org/10.1093/mnras/stac1763
  17. Brown, G., & Rein, H. (2023). General relativistic precession and the long-term stability of the Solar system. Monthly Notices of the Royal Astronomical Society 521, 4349-4355. DOI: 10.1093/mnras/stad719. arXiv:2303.05567. https://doi.org/10.1093/mnras/stad719
  18. Zeebe, R. E., & Kocken, I. J. (2025). Reduced Solar Quadrupole Moment Compensates for Lack of Asteroids in Long-term Solar System Integrations. The Astronomical Journal 170, 20. DOI: 10.3847/1538-3881/add5dc. arXiv:2505.00187. https://doi.org/10.3847/1538-3881/add5dc
  19. Rein, H., Dey, K., & Tamayo, D. (2026). Democratic Heliocentric Coordinates Underestimate the Rate of Instabilities in Long-term Integrations of the Solar System. The Astrophysical Journal 1004, 197. DOI: 10.3847/1538-4357/ae7363. arXiv:2601.08019. https://doi.org/10.3847/1538-4357/ae7363
  20. Hernandez, D. M. (2026). Artificial Precession and Instability in Solar System and Planetary Simulations: Analytic and Numerical Results. The Astrophysical Journal 1005, 25. DOI: 10.3847/1538-4357/ae6dad. arXiv:2603.24456. https://doi.org/10.3847/1538-4357/ae6dad
  21. Kaib, N. A., & Raymond, S. N. (2025). The influence of passing field stars on the solar system's dynamical future. Icarus 439, 116632. DOI: 10.1016/j.icarus.2025.116632. arXiv:2505.04737. https://doi.org/10.1016/j.icarus.2025.116632