Long-term direct N-body simulations of planetary systems almost universally use symplectic integrators with a fixed timestep, yet the timestep below which ensemble statistics become physically trustworthy had been established only for the Solar System. Rein, Brown & Kanda (2025) showed that for the secular-dominated Solar System, fixed timesteps up to 6.5 days are statistically converged over 5 Gyr ensembles, and even 32-day timesteps are only partially converged (the ensemble instability fraction rising from about 1% to about 2%), because the system's physical chaotic diffusion dominates numerical diffusion and timestep resonances. They explicitly cautioned that these thresholds rely on the Solar System's separation of secular and orbital timescales and should not be applied to planetary architectures with strong mean-motion resonances or different mass distributions. The open problem is to determine the timestep and integrator conditions under which long-term ensemble simulations of such resonant and otherwise non-Solar-System architectures are statistically converged, and to establish the mechanism — in particular whether timestep resonances and numerical diffusion, no longer shielded by a timescale separation, can dominate the physical chaotic diffusion on which the Solar System's surprising robustness rests.

Background

Direct N-body integration — numerically solving Newton's equations of gravitation for all bodies simultaneously — is the standard tool for studying the long-term evolution of planetary systems. Almost all such simulations use symplectic integrators, which preserve the Hamiltonian phase-space structure and therefore exhibit no long-term drift in energy. The most common family is the Wisdom–Holman (WH) splitting, which divides the Hamiltonian into analytically solvable Keplerian two-body parts and a perturbing interaction part, and advances each with a fixed timestep dtdt chosen by the user. A larger timestep makes simulations cheaper, but if it is too large the results are unphysical: in particular, so-called timestep resonances — spurious resonances that occur when the number of timesteps per orbit is close to a small rational ratio — can artificially pump eccentricities (Wisdom 2015). Short-term studies in clean test problems suggest one may need on the order of 20 timesteps per pericenter passage for reliability (Wisdom 2015; Hernandez, Zeeze & Hadden 2022). Whether such criteria are the right ones for billion-year integrations of full planetary systems was, surprisingly, untested before Rein, Brown & Kanda (2025). The relevant dynamics differ between system types. Mean-motion resonances (MMRs) are commensurabilities between orbital periods; secular resonances are commensurabilities between the precession frequencies of the orbital ellipses, and drive the slow exchange of angular momentum that governs the Solar System's long-term chaos. The Solar System is chaotic with a Lyapunov time (the timescale for infinitesimal perturbations to grow by a factor ee) of only about 5 Myr, so only ensemble statistics — not individual trajectories — are predictable on gigayear timescales. After the exponential divergence saturates, the system undergoes chaotic diffusion, a slow random walk of orbital elements; ensemble studies (Laskar & Gastineau 2009) find a roughly 1% probability of a catastrophic encounter or ejection within 5 Gyr. A simulation ensemble is called statistically converged when its statistics (instability fraction, instability-time distribution, fundamental secular frequencies) are independent of the timestep. Few-body studies had shown that ensemble statistics can converge even when individual trajectories do not, provided the timestep and energy error are small enough (Boekholt & Portegies Zwart 2015; Trani et al. 2024, for the three-body problem), and Rein, Brown & Tamayo (2019) showed that the traditional energy-error metric is unreliable as an accuracy criterion for secularly evolving systems, identifying an artificial precession term in the shadow Hamiltonian (the slightly modified Hamiltonian that a symplectic integrator exactly preserves) as the dominant error. Rein, Brown & Kanda (2025) then performed the first systematic statistical convergence study for the full Solar System: 1200 integrations of 5 Gyr with timesteps log-spaced between 3 and 60 days, using the WHFast implementation of the Wisdom–Holman method and the higher-order corrected WHCKL variant, both from the open-source REBOUND integration package. They found three regimes — statistically converged for dt≤6.5dt\le 6.5 days, partially statistically converged up to 32 days (instability fraction raised from about 1% to about 2%, instabilities still occurring late), and clearly unconverged beyond. Their explanation is that the Solar System's secular timescales are at least two orders of magnitude longer than the timescales on which timestep resonances act (Tamayo & Hadden 2025), and the physical chaotic diffusion washes out numerical artifacts; the system also lacks strong mean-motion resonances. It is exactly this protective separation of timescales that fails for other planetary architectures, which is why the authors caution that their thresholds 'should not be blindly applied' to resonant systems — and why the corresponding question for those systems is open.

Problem Statement

Determine the conditions — in particular the minimum fixed timestep dtdt — under which long-term direct N-body ensemble simulations of planetary architectures other than the secular-dominated Solar System are statistically converged. The architectures of primary interest are systems with strong mean-motion resonances (for example compact resonant chains of super-Earths, i.e. systems of several planets of a few Earth masses on pairwise near-commensurable orbital periods) and systems with substantially different mass distributions, for which the Solar System's separation between secular timescales and the timestep does not exist. Concretely, the problem requires: (i) ensemble experiments integrating representative such architectures over physically long timescales (gigayears, or the equivalent in dynamical times for compact systems) with a range of fixed timesteps, identifying the boundary between statistically converged, partially statistically converged, and unconverged regimes as judged by physically meaningful observables — instability fraction and instability-time distribution, resonance-angle and libration statistics, and fundamental or resonant frequencies; (ii) establishing the mechanism, i.e. whether and when timestep resonances and timestep-dependent numerical diffusion dominate over the physical chaotic diffusion that rescues the Solar System's large-timestep results, and how the answer depends on resonance strength and mass distribution; and (iii) stating the practical consequences — which timestep, integrator, and coordinate choices can be trusted, and at what computational cost, for long-term stability studies of such systems. As the source emphasizes, the Solar System thresholds (dt≤6.5dt\le 6.5 days for full convergence, 32 days as the partial-convergence ceiling) should not be blindly applied to these architectures; the answer must be established for them directly.

References

  1. Rein, Hanno; Brown, Garett; Kanda, Mei (2025). 'On the statistical convergence of N-body simulations of the Solar System.' The Open Journal of Astrophysics 8. DOI: 10.33232/001c.154745. arXiv:2507.04987. https://arxiv.org/abs/2507.04987
  2. Rein, Hanno; Dey, Kavi; Tamayo, Daniel (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://arxiv.org/abs/2601.08019
  3. Hernandez, David 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://arxiv.org/abs/2603.24456
  4. Laskar, J.; Gastineau, M. (2009). 'Existence of collisional trajectories of Mercury, Mars and Venus with the Earth.' Nature 459, 817. DOI: 10.1038/nature08096. https://doi.org/10.1038/nature08096
  5. Wisdom, J. (2015). 'Resolving the pericenter.' The Astronomical Journal 150, 127. DOI: 10.1088/0004-6256/150/4/127. https://doi.org/10.1088/0004-6256/150/4/127
  6. Hernandez, D. M.; Zeeze, R. E.; Hadden, S. (2022). 'Stepsize errors in the N-body problem: discerning Mercury's true possible long-term orbits.' Monthly Notices of the Royal Astronomical Society 510, 4302. DOI: 10.1093/mnras/stab3664. https://doi.org/10.1093/mnras/stab3664
  7. Boekholt, T.; Portegies Zwart, S. (2015). 'On the reliability of N-body simulations.' Computational Astrophysics and Cosmology 2, 2. DOI: 10.1186/s40668-014-0005-3. https://doi.org/10.1186/s40668-014-0005-3
  8. Trani, Alessandro A.; Leigh, Nathan W. C.; Boekholt, Tjarda C. N.; Portegies Zwart, Simon (2024). 'Isles of regularity in a sea of chaos amid the gravitational three-body problem.' Astronomy & Astrophysics 689, A24. DOI: 10.1051/0004-6361/202449862. arXiv:2403.03247. https://arxiv.org/abs/2403.03247
  9. Rein, H.; Brown, G.; Tamayo, D. (2019). 'On the accuracy of symplectic integrators for secularly evolving planetary systems.' Monthly Notices of the Royal Astronomical Society 490, 5122. DOI: 10.1093/mnras/stz2942. arXiv:1908.03468. https://arxiv.org/abs/1908.03468
  10. Tamayo, Daniel; Hadden, Samuel (2025). 'A Unified, Physical Framework for Mean Motion Resonances.' The Astrophysical Journal 986, 11. DOI: 10.3847/1538-4357/adc1c4. https://doi.org/10.3847/1538-4357/adc1c4