Numerical integrations of the Solar System have been carried out for decades. Their results have been used, for example, to determine whether the Solar System is chaotic, whether Mercury’s orbit is stable, or to help discern Earth’s climate history. We argue that all of the past studies we consider in this work are affected by numerical chaos to different degrees, affecting the possible orbits and instability probability of Mercury, sometimes significantly. We show how to eliminate the effects of numerical chaos by resolving Mercury’s pericentre passage. We also show that several higher order symplectic maps do not exhibit significant differences in resolving pericentre passage of Mercury (at fixed time step), making their advantages suspect for calculating long-term orbits. Resolving pericentre passage affects a wide array of orbital numerical studies, like exoplanet studies, studies of the galactic centre, and other N-body problems.

The long-term stability of the solar system is an issue of significant scientific and philosophical interest. The mechanism leading to instability is Mercury’s eccentricity being pumped up so high that Mercury either collides with Venus or is scattered into the Sun. Previously, only three five-billion-year N-body ensembles of the solar system with thousands of simulations have been run to assess long-term stability. We generate two additional ensembles, each with 2750 members, and make them publicly available at https://archive.org/details/@dorianabbot. We find that accurate Mercury instability statistics can be obtained by (1) including only the Sun and the eight planets, (2) using a simple Wisdom–Holman scheme without correctors, (3) using a basic representation of general relativity, and (4) using a time step of 3.16 days. By combining our solar system ensembles with previous ensembles, we form a 9601-member ensemble of ensembles. In this ensemble of ensembles, the logarithm of the frequency of a Mercury instability event increases linearly with time between 1.3 and 5 Gyr, suggesting that a single mechanism is responsible for Mercury instabilities in this time range and that this mechanism becomes more active as time progresses. Our work provides a robust estimate of Mercury instability statistics over the next five billion years, outlines methodologies that may be useful for exoplanet system investigations, and provides two large ensembles of publicly available solar system integrations that can serve as test beds for theoretical ideas as well as training sets for artificial intelligence schemes.

N-body integrations are used to model a wide range of astrophysical dynamics, but they suffer from errors which make their orbits diverge exponentially in time from the correct orbits. Over long time-scales, their reliability needs to be established. We address this reliability by running a three-body planetary system over about 200 e-folding times. Using nearby initial conditions, we can construct statistics of the long-term phase-space structure and compare to rough estimates of resonant widths of the system. Our statistics are approximately consistent for a wide range of numerical methods, including a Runge-Kutta method, Wisdom-Holman method, symplectic corrector methods, and a method by Laskar & Robutel. "Improving" an integrator did not affect the phase space accuracy, but simply increasing the number of initial conditions did.