On the statistical convergence of N-body simulations of the Solar System

TL;DR

Long-term N-body Solar System simulations with fixed timesteps risk yielding unphysical results if the timestep is too large. The authors perform a statistical convergence study by running 1200 five-gigayear integrations across timesteps from $3$ to $60$ days using a high-order symplectic integrator with general relativistic corrections, and they analyze instability rates, Mercury's secular frequency $g_1$, and energy conservation. They find that timesteps up to $dt=32$ days produce statistically converged ensembles for key observables, with an instability rate near $1\%$ and Mercury's secular frequencies matching high-accuracy references; larger timesteps begin to introduce non-converged behavior and numerical diffusion. These results validate much of the literature's statistical conclusions and provide practical guidance for efficient, trusted long-term Solar System simulations by ensuring dt-independence of the observed statistics.

Abstract

Most direct N-body integrations of planetary systems use a symplectic integrator with a fixed timestep. A large timestep is desirable in order to speed up the numerical simulations. However, simulations yield unphysical results if the timestep is too large. Surprisingly, no systematic convergence study has been performed on long (Gyr) timescales. In this paper we present numerical experiments to determine the minimum timestep one has to use in long-term integrations of the Solar System in order to recover the system's fundamental secular frequencies and instability rate. We find that timesteps of up to 32 days, i.e. a third of Mercury's orbital period, yield physical results in an ensemble of 5 Gyr integrations. We argue that the chaotic diffusion that drives the Solar System's long-term evolution dominates over numerical diffusion and timestep resonances. Our results bolster confidence that the statistical results of most simulations in the literature are indeed physical and provide guidance on how to run time and energy efficient simulations while making sure results can be trusted.

Figures (5)

• Figure 1: Top panel: the fraction of simulations that went unstable within 5 Gyr. There are 50 simulations in each timestep bin. Bottom panel: the time of instability (black dots) for individual simulations. The colours indicate the three regimes we identify: statistically converged (green), partially statistically converged (blue), not statistically converged (orange).
• Figure 2: The standard deviation of the secular frequency $g_1$ as a function of time for different timestep bins. The colours indicate the three regimes and are the same as in Fig. \ref{['fig:instability-times']}.
• Figure 3: Top: Cumulative distributions of the secular frequencies $g_1$ through $g_5$ in simulations of the Solar System after 5 Gyr with different timesteps. The black line corresponds to integrations using of the high order WHCKL integrator. Bottom: p-value as a function of timestep for distributions of $g_1$ through $g_5$. We compare the distributions from each timestep bin to the distribution obtained with WHCKL and a 6 day timestep. The dotted line corresponds to $p=0.05$.
• Figure 4: Cumulative distribution of relative energy errors at the end of the 5 Gyr integration for each timestep bin.
• Figure 5: The relative energy error as a function time for a randomly chosen simulation from each timestep bin. Also plotted is a line proportional to the square root of time which is the expected time dependence for a simulation that follow Brouwer's law.