On the accuracy of symplectic integrators for secularly evolving planetary systems

TL;DR

This study evaluates the accuracy of symplectic integrators in long-term planetary system simulations, revealing that certain metrics and higher-order methods better capture secular frequencies by addressing specific error terms.

Contribution

It introduces a framework to identify error sources in symplectic integrators and demonstrates the superiority of higher-order methods for secular frequency accuracy.

Findings

Symplectic correctors do not improve secular frequency accuracy despite smaller energy errors.

A specific shadow Hamiltonian term causes artificial precession affecting frequency calculations.

Higher-order symplectic methods outperform standard ones in secularly evolving systems.

Abstract

Symplectic integrators have made it possible to study the long-term evolution of planetary systems with direct N-body simulations. In this paper we reassess the accuracy of such simulations by running a convergence test on 20Myr integrations of the Solar System using various symplectic integrators. We find that the specific choice of metric for determining a simulation's accuracy is important. Only looking at metrics related to integrals of motions such as the energy error can overestimate the accuracy of a method. As one specific example, we show that symplectic correctors do not improve the accuracy of secular frequencies compared to the standard Wisdom-Holman method without symplectic correctors, despite the fact that the energy error is three orders of magnitudes smaller. We present a framework to trace the origin of this apparent paradox to one term in the shadow Hamiltonian. Specifically, we find a term that leads to negligible contributions to the energy error but introduces non-oscillatory errors that result in artificial periastron precession. This term is the dominant error when determining secular frequencies of the system. We show that higher order symplectic methods such as the Wisdom-Holman method with a modified kernel or the SABAC family of integrators perform significantly better in secularly evolving systems because they remove this specific term.

Tables

Table 1. Table 1: Secular frequencies obtained by our integrations and other authors. All frequencies are given in units of arc-seconds per year. The frequencies obtained in this work have been rounded to the same number of digits as in Laskar et al. (2011) for easier comparison.

	$g_1$	$g_2$	$g_3$	$g_4$	$g_5$	$g_6$	$g_7$	$g_8$
This work	5.59	7.422	17.360	17.914	4.257519	28.2457	3.088043	0.673301
Laskar et al. (2011), La2010a	5.59	7.453	17.368	17.916	4.257482	28.2449	3.087946	0.673019
Spalding et al. (2018)	5.71	7.44	17.19	17.76	4.26	28.25	3.09	0.67

Equations

$$\{g,h\}\equiv \sum _i\left(\frac{\partial g}{\partial q_i}\frac{\partial h}{\partial p_i}-\frac{\partial g}{\partial p_i}\frac{\partial h}{\partial q_i}\right),$$

$${\dot{q}}_i=\frac{\partial H}{\partial p_i}=\{q_i,H\}\text{and}{\dot{p}}_i=-\frac{\partial H}{\partial q_i}=\{p_i,H\}.$$

$$\dot{g}=\{g,H\}.$$

$${\mathcal{L}}_{H}g\equiv \{g,H\}.$$

$$\dot{y}={\mathcal{L}}_{H}y$$

$$\dot{y}={\mathcal{L}}_{A}y+{\mathcal{L}}_{B}y.$$

$${\Psi }_t(y_0)={\phi }_{t/2}^{[A]}\circ {\phi }_t^{[B]}\circ {\phi }_{t/2}^{[A]}(y_0).$$

$${\phi }_t^{[H]}(y_0)=\exp \left(t{\mathcal{L}}_H\right)\text{Id}(y_0)=\exp \left(t\left({\mathcal{L}}_A+L_B\right)\right)\text{Id}(y_0).$$

$$\displaystyle\varPsi_{t}(y_{0})=\exp\left(\frac{t}{2}{\mathcal{L}}_{A}\right)\exp\Big{(}t{\mathcal{L}}_{B}\Big{)}\exp\left(\frac{t}{2}{\mathcal{L}}_{A}\right)\text{Id}(y_{0}).$$

$$\displaystyle\varPsi_{t}(y_{0})=\exp\Big{(}t{\mathcal{L}}_{A}+t{\mathcal{L}}_{B}-\frac{t^{3}}{24}\left[{\mathcal{L}}_{A},\left[{\mathcal{L}}_{A},{\mathcal{L}}_{B}\right]\right]$$

$$\displaystyle-\frac{t^{3}}{12}\left[{\mathcal{L}}_{B},\left[{\mathcal{L}}_{A},{\mathcal{L}}_{B}\right]\right]+\mathcal{O}(t^{5})\Big{)}\text{Id}(y_{0}),$$

$$\left[{\mathcal{L}}_A,{\mathcal{L}}_B\right]g(y)={\mathcal{L}}_{\{A,B\}}g(y).$$

$$z_i(t)=e_i(t)\exp \left(i{\varpi }_i(t)\right)i=1,2,\dots ,8,$$

$$\mathbf{u}(t)=S^{-1}\mathbf{z}(t).$$

$$⟨e^{i{\omega }_{i,\max }t},u_i(t)⟩=\int e^{i{\omega }_{i,\max }t}{u}_i(t)h(t)dt,$$

$$E_{\mathrm{floor}}\approx 2^{-53}\sqrt{N_{\mathrm{step}}},$$

$${\mathcal{E}}_t(y_0)\equiv {\phi }_t(y_0)-{\Psi }_t(y_0).$$

$${\mathcal{E}}_t$$

$${\mathcal{C}}_t^{AAB}$$

$${\mathcal{C}}_t^{AB}$$

$${\mathcal{C}}_t^{BA}$$

$${\mathcal{C}}_t^{AAB}$$

$${\mathcal{C}}_t^{AAB}(y_0)$$

$$\displaystyle=\Big{(}\text{Id}+t^{3}\left[{\mathcal{L}}_{A},\left[{\mathcal{L}}_{A},{\mathcal{L}}_{B}\right]\right]\text{Id}+\mathcal{O}(t^{4})\Big{)}(y_{0}).$$

$$\displaystyle{\mathcal{E}}_{t}(y_{0})=\Big{(}{\mathcal{C}}^{AAB}_{\alpha t}+{\mathcal{C}}^{BAB}_{\beta t}-2\cdot\text{Id}+\mathcal{O}(t^{4})\Big{)}(y_{0}),$$

$$t^3{\mathcal{L}}_{\left\{A,\left\{A,B\right\}\right\}}\text{Id}(y_0)=\left({\mathcal{C}}_t^{AAB}-\text{Id}\right)(y_0)+\mathcal{O}(t^4)$$

$${\phi }_t^{\left[\left\{A,\left\{A,B\right\}\right\}\right]}(y_0)=\frac{1}{t^2}{\mathcal{C}}_t^{AAB}(y_0)+\mathcal{O}(t^2).$$

$${\Delta }_{t}E(y_0)=H(y_0)-H\left({\mathcal{C}}_{\alpha t}^{AAB}(y_0)\right)+\mathcal{O}(t^2),$$

$${\Delta }_t\varpi (y_0)=\varpi (y_0)-\varpi \left({\mathcal{C}}_{\alpha t}^{AAB}(y_0)\right)+\mathcal{O}(t^2).$$

Full text

1 Introduction

The range of timescales present in planetary systems is truly astronomical, from a few days to billions of years. Specialized integrators are therefore needed to accurately calculate trajectories over the entire lifetime of planetary systems. In recent years, the development of such accurate numerical methods and the increase in computing power have made very long-term integrations possible (Laskar & Gastineau, 2009).

The mixed variable symplectic Wisdom & Holman (1991) integrator (WH) and higher order variants thereof are one of the preferred tools for this purpose (Wisdom et al., 1996; Laskar & Robutel, 2001). Recently, Rein et al. (2019b) implemented many of these higher order symplectic integrators in the freely available REBOUND package. This present paper uses several of these methods to perform a convergence study for integrations of the Solar System. Typically authors report the performance of a method by monitoring how well it can conserve energy, as energy should be conserved exactly in a conservative system (Laskar & Robutel, 2001; Wisdom, 2018; Rein & Tamayo, 2015). Although the energy error is a metric involving all planets, it is dominated by only a few. We show that some of these methods, despite reducing the energy error by three orders of magnitude, do not improve the accuracy of the apsidal precession frequencies that drive the solar system’s long-term evolution. In this paper, we present a framework which lets us resolve this apparent paradox.

We introduce the notation and describe the numerical setup and algorithms in Sect. 2. This includes a short overview of the symplectic integrators that we use, as well as an outline of the modified Fourier transform that we use to determine the secular frequencies of the system. The results of our convergence study of long-term integrations are presented in Sect. 3. We then develop a framework which lets us understand the origin of errors in symplectic methods in Sec. 4. We conclude with a summary and discussion of the implications in Sec. 5.

2 Methods

We use the freely available $N$-body package REBOUND (Rein & Liu, 2012) to perform simulations of the Solar System with all 8 planets. We model the planets as point masses and treat the Earth-Moon system as one particle in our simulation. The initial conditions are taken from the NASA Horizons database, provided by the Solar System Dynamics Group of the Jet Propulsion Laboratory111https://ssd.jpl.nasa.gov. The effects of general relativity are approximated by including a $1/r^{3}$ term in the potential (Nobili & Roxburgh, 1986) using the implementation provided by REBOUNDx (Tamayo et al., in prep222Code and documentation are already available at https://reboundx.readthedocs.io).

We integrate the simulations for 20 million years. Although the direction in which we integrate does not matter, we here integrate backwards so that we can more easily compare the secular frequencies to those obtained by Laskar et al. (2011). We use the Simulation Archive (Rein & Tamayo, 2017) to record 22000 snapshots throughout each integration. Since our goal is to compare simulations with different timesteps to very high precision, we have to ensure that we can record the snapshots at exactly the same times for all simulations without the need for any interpolation or change of the timestep during the integration. We achieve this by choosing the timesteps such that all snapshots occur at integer multiples of the timestep (the same procedure was used by Laskar et al., 2011). Furthermore, we count the number of timesteps taken rather than keeping track of the integration time itself. This ensures that the round-off errors originating in adding small numbers (the timestep) to a large number (the current simulation time) do not affect the analysis.

We use the IAS15 integrator, a high order Gauß-Radau integrator (Rein & Spiegel, 2015), to obtain a true solution to which we can then compare all other simulations. It is important to point out that we only consider it to be the true solution of our model, but not the true solution of the real Solar System. All we need for our convergence study is a well defined model that contains the most important dynamical effects determining the evolution of the Solar System. Other simulations, in particular those by Laskar et al. (2011), are better representations of the real evolution of the Solar System because they use more accurate initial conditions and include important physical effects such as tidal evolution, stellar mass loss, and perturbations from minor planets. We come back to why we have confidence that the IAS15 integrator indeed provides the true solution to our model later, and also comment on it in Appendix A.

We run simulations with the exact same initial conditions using the IAS15 integrator, the standard WH integrator, the WH integrator with first symplectic correctors of order 17, WHC (Wisdom et al., 1996), the WH integrator with the lazy implementer’s kernel method, WHCKL (Wisdom et al., 1996), and the SABA integrator with four function evaluations and correctors using the lazy implementer’s method, SABACL4 (Laskar & Robutel, 2001; Rein et al., 2019b).

Next, we introduce the notation that we use for operator splitting methods and then summarize the main features of the various integrators mentioned. The differences between the various symplectic integrators and specifically their implementation in REBOUND are also described in more detail in Rein et al. (2019b).

2.1 Operator Splitting Schemes

Let us begin by defining the Poisson bracket $\{g,h\}$ of two functions $g$ and $h$ of the canonical coordinates $(q_{i},p_{i})$ as

[TABLE]

which allows us to write down Hamilton’s equations as

[TABLE]

The time derivative of any function $g$ that depends only on $p$ and $q$ can then be written succinctly as

[TABLE]

Let us further introduce ${\mathcal{L}}_{H}$, an operator which can act on any phase space function $g$. We call this operator the Lie derivative with respect to the Hamiltonian $H$ and define it as

[TABLE]

Assume that the dynamical system we are interested in has $3N$ coordinates and $3N$ momenta. The Lie derivative allows us to write Hamilton’s equation in a compact way as

[TABLE]

where $y(t)\equiv(q_{1}(t),\ldots,q_{3N}(t),p_{1}(t),\ldots,p_{3N}(t))$ are the canonical coordinates and momenta.

Lastly, let us define the formal solution operator $\varphi^{[H]}_{t}(y_{0})$. This operator returns the solution to the differential equation $\dot{y}={\mathcal{L}}_{H}y$ at time $t$ with initial conditions $y_{0}$ given at $t=0$.

Because we do in general not know how to write down this abstract solution operator, we use an operator splitting method to approximate it. The idea behind an operator splitting method is to rewrite the differential equation in Eq. 5, to which we do not know the solution, as

[TABLE]

By virtue of the linearity of the Lie derivative, it follows that this equation is identical to Eq. 5 if the have $H=A+B$. This idea of expressing $H$ as a sum of $A$ and $B$ is often referred to as splitting the Hamiltonian. Let us refer to the evolution operators corresponding to the differential equations $\dot{y}={\mathcal{L}}_{A}y$ and $\dot{y}={\mathcal{L}}_{B}y$ as ${\varphi}^{[A]}_{t}$ and ${\varphi}^{[B]}_{t}$ respectively. Then, we can write down a second order (leap-frog) operator splitting scheme,

[TABLE]

As we will see below, this is a second order method which approximates the evolution operator ${\varphi}^{[H]}_{t}$. Note that in the above notation, the $\varphi$s are operators and we apply the right most operator first.

Next, we use a formal trick to express any solution operator $\varphi^{[X]}$ in terms of its corresponding Lie derivative operator ${\mathcal{L}}_{X}$ and an exponential (defined as an infinite series in the usual way), for example,

[TABLE]

To see that this is indeed a solution to the differential equation in Eq. 5, note that the series expansion of the exponential provides a series of differential operators acting on $\text{Id}(y_{0})$. The series is simply the Taylor expansion of the solution ${\varphi}^{[H]}$ around $t=0$. We can use this formalism to also express our splitting method from Eq. 7 in terms of exponentials of Lie derivative operators, i.e.

[TABLE]

This allows us to use the symmetric Baker-Campbell-Hausdorff (BCH) identity to combine the three exponentials of non-commuting operators into one exponential:

[TABLE]

where the commutator $\left[{\mathcal{L}}_{A},{\mathcal{L}}_{B}\right]$ is short hand for ${\mathcal{L}}_{A}{\mathcal{L}}_{B}-{\mathcal{L}}_{B}{\mathcal{L}}_{A}$.

If we compare Eq. 8 with Eq. 10, we can see that the differential equations which our splitting scheme is solving are not those originating from the Hamiltonian $H$. Although the first two terms are equal, Eq. 10 also contains an infinite series of commutators of Lie derivatives.

Note that by virtue of the Jacobi identity, we have

[TABLE]

This allows us to identify the commutators of Lie derivatives in Eq. 10 as Poisson brackets of the Hamiltonians $A$ and $B$. We can then construct a shadow Hamiltonian $\tilde{H}$ which is the original Hamiltonian $H$ plus an infinite series of Poisson brackets. For example, the two terms at order $\mathcal{O}\left(t^{3}\right)$ in Eq. 10 correspond to the Poisson brackets $t^{2}\{A,\{A,B\}\}$ and $t^{2}\{B,\{A,B\}\}$ in the shadow Hamiltonian.

Note that all of these terms contain a factor of $t^{n}$ with $n\geq 2$. In the limit of $t\rightarrow 0$ all error terms are therefore converging333There is an issue of whether the series actually converges (see Wisdom, 2018). This is not important for the discussion in this paper. to 0. In practice we will choose a timestep smaller than the shortest dynamical timescale in the problem. We indicate that we take a small but finite timestep by replacing $t$ with $dt$ from now on.

2.2 Wisdom-Holman Mapping

The Wisdom-Holman (WH) integrator (Wisdom & Holman, 1991) is using the scheme given in Eq. 7 together with a specific splitting of the Hamiltonian $H$. In the WH integrator, part $A$ describes the evolution of planets on Keplerian orbits around the Sun and part $B$ describes the much smaller interplanetary forces. The exact solutions $\varphi^{[A]}$ and $\varphi^{[B]}$ of the equations of motions associated with $A$ and $B$, $\dot{y}={\mathcal{L}}_{A}y$ and $\dot{y}={\mathcal{L}}_{B}y$, can be calculated efficiently. $\varphi^{[A]}$ is referred to as the Kepler step and involves solving Kepler’s equation. This can be done using a series expansion which can be truncated once machine precision is reached. $\varphi^{[B]}$ is referred to as the kick or interaction step. It updates the velocities of planets by evaluating the planet-planet interaction potential.

Let us estimate the size of the additional terms in the shadow Hamiltonian of the WH method. Because the WH integrator splits the Hamiltonian in a dominant Keplerian part $A$ and a perturbation part $B$, we can introduce a dimensionless parameter to keep track of this separation of scales by formally replacing $B$ with $\epsilon B$, where $\epsilon\ll 1$. In the Solar System, the gravitational forces from other planets are roughly a factor of $10^{-3}$ smaller than the force from the Sun, thus in this particular case $\epsilon\approx 10^{-3}$. Using this substitution, the two terms in the shadow Hamiltonian involving $\{A,\{A,B\}\}$ and $\{B,\{A,B\}\}$ are of order $\mathcal{O}\left(\epsilon dt^{2}\right)$ and $\mathcal{O}\left(\epsilon^{2}dt^{2}\right)$, respectively. Because a factor of $\epsilon$ is present in all error terms, the WH integrator performs well in weakly perturbed systems where $\epsilon$ is small.

2.3 Wisdom-Holman Mapping with Symplectic Correctors

The standard WH method is a second order method in $dt$. Symplectic methods can also be constructed to higher order in the timestep (Yoshida, 1990). To improve the accuracy for a fixed $dt$, one can thus simply choose a higher order method to remove the dominant error terms in Eq. 10. However, this requires more force evaluations per timestep and is therefore slower. Symplectic correctors were introduced by Wisdom et al. (1996) as an alternative to improve the accuracy with effectively no overhead. They can be understood as a clever transformation between the original problem and integrator variables which removes terms of order $\mathcal{O}(\epsilon)$ from the shadow Hamiltonian. Because the transformations between real and mapping variables only need to be performed when outputs are needed (not every timestep), they incur negligible computational cost. For example a third-order corrector removes the $\{A,\{A,B\}\}$ term which is of order $\mathcal{O}(\epsilon dt^{2})$. A fifth-order corrector removes all the terms up to and including those of order $\mathcal{O}(\epsilon dt^{4})$, and so forth. In this paper, we use the symplectic correctors of order 17 and refer to the method as WHC. The now dominant term in the shadow Hamiltonian for small timesteps is $\{B,\{A,B\}\}$, which is of order $\mathcal{O}\left(\epsilon^{2}dt^{2}\right)$.

2.4 Wisdom-Holman Mapping with a modified Kernel

If we wish to further improve the WHC method, we have to remove the now leading order term, $\{B,\{A,B\}\}\sim\mathcal{O}(\epsilon^{2}dt^{2})$. It goes beyond the scope of this paper to describe in detail how Wisdom et al. (1996) achieve this using a modified kernel. What is important for our analysis is that this new method, WHCKL, has now a leading order term of $\mathcal{O}(\epsilon^{2}dt^{4})$ by getting rid of the $\{B,\{A,B\}\}\}$ term. As summarized in Rein et al. (2019b), there are different ways the kernel can be implemented. In this paper, we use the lazy implementer’s method. In comparison to the modified kick method WHCKM, it allows for the inclusion of non-Newtonian forces. And in comparison to the composition method WHCKC, it is significantly faster (Rein et al., 2019b).

2.5 SABA integrator family

We also test the symplectic integrators of the SABA family (Laskar & Robutel, 2001). Specifically, we use the SABACL4 method which uses four force evaluations and one corrector per timestep (Rein et al., 2019b). The correctors in SABACL4 are responsible for removing the $\{B,\{A,B\}\}$ term in the shadow Hamiltonian. Note that although both Wisdom et al. (1996) and Laskar & Robutel (2001) use the word corrector to describe parts of their integrators, the derivation, interpretation, and implementation of the correctors differ significantly.

The SABACL4 method has a leading order error term of $\mathcal{O}(\epsilon dt^{6}+\epsilon^{2}dt^{4})$. For small timesteps the $\mathcal{O}(\epsilon^{2}dt^{4})$ term dominates. Most importantly for our later discussion, this method got rid of the $\{B,\{A,B\}\}$ error term in the shadow Hamiltonian.

2.6 Frequency Analysis

We follow the procedure of Laskar (1988, 1990, 1993, 2003) to determine the secular frequencies of the Solar System using a modified Fourier transform. Specifically we will compute the $g$-modes of the Solar System, which are related to the periastron and eccentricity of the planets. A full description of the algorithm is given in the references above (see also Šidlichovský & Nesvorný, 1996). Here, we only present a short outline of the algorithm.

We start by calculating the orbital elements in Jacobi coordinates for all planets at each of our snapshots. For simulations that use a symplectic corrector, we need to apply the correctors beforehand. This is done by calling the synchronize routine in REBOUND. We then compute the complex eccentricities,

[TABLE]

where $e_{i}(t)$ and $\varpi_{i}(t)$ are the eccentricity and the longitude of periastron of the $i$-th planet at time $t$.

Next, we perform a linear transformation to a basis that approximately corresponds to the proper modes. Although not absolutely necessary, this transformation makes it easier to identify the dominant secular frequencies in each planet’s orbital parameters. We use the $S$ matrix published by Laskar (1990) to do the transformation to proper modes,

[TABLE]

Since the Solar System’s secular frequencies are changing very slowly over time, this transformation should be thought of as an approximation. For our purpose, the precision of this approximation is not relevant as long as we use the same transformation for the entire convergence study.

On each proper mode $u_{i}(t)$, we then perform a Fourier transform and find the frequency $\omega_{i,\rm fft}$ with the maximum power. This frequency can only be determined to within one Nyquist frequency using a standard Fourier transform. We would need an unreasonably large number of samples ($10^{10}$) to achieve the accuracy we need for our convergence study ($10^{-10}$). To get around this, we search in the neighbourhood of $\omega_{i,\rm fft}$ for the frequency $\omega_{i,\rm max}$ which maximizes the following scalar product

[TABLE]

where $h(t)$ is the Hanning window function. The integral is evaluated at the discrete sampling points using the trapezoidal rule. We refer to the frequency $\omega_{i,\rm max}$ as the secular frequency $g_{i}$. Whereas one can continue with an orthogonalization procedure to extract other frequencies in each proper mode, we here only determine one frequency per mode.

3 Long term integrations of the Solar System

In this section we present the results of our 20 Myr integrations of the Solar System. With these results as a motivation, we then develop a formalism to explain them in Sec. 4.

3.1 Comparing our secular frequencies to previous work

As a consistency check, we start by comparing the secular frequencies obtained with our IAS15 simulation to previous work of Laskar et al. (2011) and Spalding et al. (2018). The results are listed in Table 1. We can reproduce the frequencies obtained by Laskar et al. (2011) with high precision, with the largest discrepancy of about 0.5% being in $g_{2}$, the frequency associated with Venus444Given the low eccentricity of Venus, it makes sense that this mode is particularly sensitive to changes in the initial conditions.. The frequencies of the outer planets are all reproduced to four decimal digits or better.

We do not expect to reproduce the frequencies of Laskar et al. (2011) exactly as these authors use a different set of initial conditions, treat the Moon as a separate object, include a model of tidal dissipation as well as the effects of some minor planets. Whereas the treatment of these additional effects are important to accurately model the long-term evolution of the Solar System, the high level of agreement we find even without including these effects gives us confidence that our simulations are an accurate model which captures the most important dynamics in the Solar System.

A further difference between the studies is that the frequencies obtained in this work come from a 20 Myr integration; the work of Laskar et al. (2011) use a 20 Myr integration for $g_{1}$ through $g_{4}$ and a 50 Myr integration for $g_{5}$ through $g_{8}$; and the work of Spalding et al. (2018) use a 450 Myr integration for all their frequencies.

We note that although Spalding et al. (2018) also uses the REBOUND integrator package and initial conditions from NASA Horizons, we were not able to reproduce their frequencies as accurately as those of Laskar et al. (2011). Possible explanations for this discrepancy might be their treatment of general relativity or a difference in their data reduction. Given that a precise integration was not necessary for the conclusions of Spalding et al. (2018) and that we agree with Laskar et al. (2011) to such high precision, we have confidence in this being the right solution.

3.2 Convergence study

Next, we present the results of our convergence study of the WH, WHC, WHCKL, and SABACL4 integrators in our 20 Myr integrations of the Solar System. For that, we vary the timestep from $0.05$ days, to $100$ days and measure various different error metrics for each simulation. For small enough timesteps, all integrators should converge to the true solution.

Fig. 1 shows the maximum relative energy error in our simulations. Since we are integrating a conservative system the energy should be conserved exactly by an ideal integrator and we can therefore measure the energy error directly without the need to know the exact solution. This is why many studies only look at the energy error. We can see that the energy error in simulations using the WH integrator converges following a $dt^{2}$ power-law until it reaches machine precision at $dt\approx 1$ day. We overplot the expected floor originating from IEEE754 double floating point round-off errors,

[TABLE]

where $N_{\rm step}$ is the number of timesteps required for a 20 million year integration (Brouwer, 1937; Quinn & Tremaine, 1990). We cannot expect to achieve a better accuracy than that because floating point numbers on a computer are represented using a finite number of digits (about 16 in decimal, 53 in binary). The above estimate assumes that one rounding error of $\sim 10^{-16}$ occurs after every timestep, ignoring other operations during the timestep. Clearly, every integrator has many operations per timestep but this provides a lower limit.

The energy error in simulations using the WHC and WHCKL integrators converge much faster and reach machine precision at $dt\approx 8$ days (the two curves for WHC and WHCKL are on top of each other). Similarly, the energy error in the SABACL4 simulations reaches machine precision at $dt\approx 20$ days. Note that the SABACL4 integrator uses four sub-timesteps and therefore four force evaluations per timestep (plus one corrector step also requiring a force evaluation). If we were to take this into account and plot the error as a function of the effective timestep, corresponding to one force evaluation per timestep, then the SABACL4 curve would shift to the left.

From this plot, we might conclude that the WHC, WHCKL, and SABACL4 integrators are thus very similar in performance and several orders of magnitude better than standard WH for intermediate timesteps.

As a reference, we also show in Fig. 1 the relative energy error in the simulation using the IAS15 integrator. Since IAS15 is an adaptive scheme, we plot the average timestep on the horizontal axis. Note that IAS15 uses compensated summation for some crucial parts of the calculation and therefore achieves an accuracy better than the IEEE754 floor. The fact that the energy error in the IAS15 simulation is significantly smaller than for the other integrators gives us confidence that we can consider it the true solution of our model for all practical purposes (see also Appendix A).

Finally, note that the energy error is really not a particularly good metric for the convergence study in the first place. There is only about one decade in the timestep, from $dt=10$ days to $dt=1~{}$day, that we can effectively use. For timesteps larger than 10 days, the methods are not well converged because we barely resolve the shortest timescale in the problem, Mercury’s orbital period of 88 days. For timesteps smaller than 1 day, all methods show an energy error which is dominated by round-off error, suggesting that if we further decrease the timestep, we do not increase the accuracy anymore (this conclusion is wrong as we show below). Part of the reason why the energy error is not a good metric for Solar System simulations is that it is a global error metric, but the Solar System has a wide range of orbital periods and a wide range of planet masses. Jupiter’s mass is more than 5000 times larger than that of Mercury, but Mercury’s orbital period is 50 times shorter. Mercury has the fewest timesteps along its orbit, but because of its low mass, the errors in its trajectory contribute comparatively little to the total energy of the system. Nevertheless, it is Mercury’s orbit which is most likely to undergo macroscopic instabilities over billions of years (Laskar & Gastineau, 2009).

Since the energy error makes it hard to isolate potential issues, we plot different quantities in Fig. 2. Specifically, we plot the semi-major axis, periastron, and secular frequency errors associated with each planet. Although all secular modes couple to each planet’s orbital elements to some degree via Eq. 13, we associated the secular frequency $g_{1}$ with Mercury, $g_{2}$ with Venus, and so forth, following traditional conventions. Unlike for the energy error, we need to compare the integrations to the true solution to calculate the error in the secular frequencies. We here consider the IAS15 integration as the true solution. The semi-major axis error is scaled relative to each planet’s initial semi-major axis.

Looking at the top row, we can see that the semi-major axis error for all planets follows qualitatively the results for the energy error. The WHC, WHCKL, and SABACL4 integrators perform significantly better than the standard WH integrator. This is expected as to leading order the energy error is simply the sum of the semi-major axes error weighted by planet mass.

However, looking at the middle and bottom row of Fig. 2, we can see that the periastron errors and the secular frequency errors behave very differently. Using those metrics, the accuracy of the WHC and WH integrators are now suddenly comparable for the inner Solar System, with WHC even having a slightly larger secular frequency error in $g_{1}$ (bottom left panel) than WH. The other two integrators in our sample, WHCKL and SABACL4, remain significantly better than the WH method.

This, the fact that the WHC method is performing so well when looking at the energy error, but so poorly555Relatively speaking, of course. when looking at quantities important for the secular evolution, is the main paradox which we will explain in Sec. 4.

Note that the secular frequency error seems to be a particularly well suited metric for this convergence study. We can now clearly see a second order power law for the WHC method when changing the timestep $dt$ over several orders of magnitude. Compare this to the energy error in Fig. 1, where we were not able to see a power law for WHC. Furthermore, note that the simulations which appeared to be dominated by round-off error for small timesteps in Fig. 1, do in fact continue to improve in accuracy if the timestep is reduced further.

3.3 High accuracy of secular frequencies

Before we go on to derive a formalism to explain these results, let us reiterate that we use several integrators which are based on conceptually different approaches to solving the evolution of planetary systems. IAS15 can be thought of as a general purpose high-order integrator, using compensated summation, variable timesteps, and not making use of the fact that we are integrating a Hamiltonian system. Then there is the SABACL4 integrator which can also be thought of as a very general purpose scheme for nearly integrable systems, not necessarily being derived for only Hamiltonian systems (Laskar & Robutel, 2001). And there is the family of WH integrators whose derivation was heavily influenced by the dynamical systems approach (Wisdom, 2018).

Above, we have shown evidence that all these different integrators agree on the secular frequencies of the Solar System to a precision of $10^{-9}$ arc-seconds per year for small enough timesteps. This might not be surprising. The integrators should converge to the same solution, after all. But given that the integrations differ significantly and consist of up to $10^{11}$ timesteps, this is nevertheless a reassuring result.

Note that determining the secular frequencies to extremely high precision does not immediately lead to new physical insight about the system. In fact, the frequencies have little obvious physical meaning beyond approximately 1 part in $10^{2}$ for the inner planets and 1 part in $10^{4}$ for the outer planets in a secular framework because they can be no longer considered constants at higher precision. That is why they are only quoted to a finite precision in Tab. 1. For example, if we were to measure the frequencies in an even slightly shorter or longer integration, they would change by a factor much larger than $10^{-9}$, the precision to which we have determine them. However, the precise physical interpretation of the secular frequencies does not matter for our convergence study. The way to think about the secular frequencies is not as fundamental constants of the Solar System. Rather, they are complicated combinations of action and angle variables which are important for the long-term evolution of the system. And because we can measure them so precisely in a numerical simulation, we can use them to verify the behaviour of our numerical methods in an idealized but well defined model of the Solar System. As we have shown, the observed behaviour of some methods in numerical experiments is not always what their formal properties might suggest.

4 Analysis

In this section we will explain why, compared to the standard WH method, the WHC method improves the energy error but not the secular frequency error. To do this, we need to introduce a bit of notation first. We use an operator based approach, with a notation similar to that in Rein et al. (2019b) and Rein et al. (2019a). For different perspectives, we refer the reader to books and papers by Hairer et al. (2006); Laskar & Robutel (2001), and Wisdom (2018) as well as references therein.

4.1 Error operator

We are interested in the error of our splitting methods. Let us therefore define an error operator as

[TABLE]

At every point in phase space, this operator returns the difference between the true solution ${\varphi}$ and the approximate solution of our splitting method $\varPsi$ after some time $t$. We can express this operator using exponentials from Eqs. 8 and 10, expanding them, and then matching terms of the same power in $t$. This quickly becomes tedious but can be done by hand for the leading order term that we’re interested in. To leading order, we have

[TABLE]

Note that all terms up to second order in $t$ cancel because we have a second order splitting $\varPsi$. To save space we have dropped the state on which the operator is acting on from the notation above.

4.2 Approximating the error operator using a composition

Next, consider the composition of operators defined by

[TABLE]

where

[TABLE]

and thus

[TABLE]

Note these operators look like commutators of group elements, in contrast to the commutators of Lie derivatives (which are part of an algebra). After repeatedly applying the BCH formula, and ignoring terms of order $\mathcal{O}(t^{4})$ and higher, we can rewrite this composition as

[TABLE]

Comparing this to Eq. 17 and using the analogous operator ${\mathcal{C}}^{BAB}$, we can explicitly calculate the error operator to leading order in $t$ via

[TABLE]

where $\alpha=1/\sqrt[3]{24}$ and $\beta=1/\sqrt[3]{12}$. Note that this allows us to calculate ${\mathcal{E}}$ without implementing any Lie brackets. We only need the evolution operators $\varphi^{[A]}$ and $\varphi^{[B]}$. Furthermore, if we are only interested in one of the leading order terms in Eq. 17 we can calculate each term individually. The calculation above can be repeated for different splitting methods. For example, the error operator of a first order method can be expressed to leading order by the composition ${\mathcal{C}}^{AB}_{t}$.

As previously discussed in Sect. 2, due to the Jacobi identity, we have $\left[{\mathcal{L}}_{A},{\mathcal{L}}_{B}\right]g(y)={\mathcal{L}}_{\{A,B\}}\,g(y)$, which allows us to identify the commutator of Lie derivatives with a corresponding Poisson bracket of Hamiltonians $A$ and $B$. Thus, we could also compute the operator ${\mathcal{C}}^{AAB}$ by first explicitly calculating the Poisson bracket $\{A,\{A,B\}\}$, then calculating, and finally solving Hamilton’s equations. In practice though, this can get tedious very quickly if calculated by hand for non-trivial Hamiltonians666An automatic differentiation algorithm could do this, but most likely not very efficiently as Lie derivatives not only contribute derivatives, but also sums over all coordinates.. On the other hand, the composition methods we introduced above are straightforward to implement using only two already existing operators. Note that these compositions are similar to the compositions used by Wisdom et al. (1996) to express their symplectic correctors777Similar to how these authors construct higher order symplectic correctors, we could improve our approximation of ${\mathcal{E}}$ by going to higher order. However, since we are interested in the leading order to identify relevant terms anyway, going to higher order is not of much use for us here..

To summarize, the composition ${\mathcal{C}}^{AAB}$ that we have introduced above is a way to approximate the evolution operator $\varphi_{t}^{\left[\left\{A,\left\{A,B\right\}\right\}\right]}$ corresponding to evolution under the Hamiltonian $\left\{A,\left\{A,B\right\}\right\}$ to leading order in $t$ at every point in phase space $y_{0}$ without the need to calculate any derivatives, no matter how complicated the Hamiltonian is. Explicitly, from Eq. 19 we have

[TABLE]

and therefore

[TABLE]

Note that the factor of $t^{-2}$ acts like a denominator in a finite difference. If the result of this operator is dominated by round-off errors, we can easily rescale $t$ on the right hand side of Eq. 22 to avoid these errors.

4.3 Using the error operator

We can use the composition operators to approximate the evolution of any function of the phase space variables under any Poisson bracket. Here, we want to estimate the evolution under the Poisson brackets which appear at leading order in the shadow Hamiltonian of a given splitting scheme.

For example, we can approximate the change in energy that occurs after one step with the WH method $\varPsi_{t}$ due to the term $[{\mathcal{L}}_{A},[{\mathcal{L}}_{A},{\mathcal{L}}_{B}]]$ in Eq. 17 as follows

[TABLE]

where $H$ is the Hamiltonian, now used as a function operating on the phase space variables and returning the energy as a scalar. And $\alpha$ is a constant taking into account the constant coefficients which appear in Eq. 17. In fact we can calculate how any quantity changes due to the $[{\mathcal{L}}_{A},[{\mathcal{L}}_{A},{\mathcal{L}}_{B}]]$ term. For example the error in the periastron $\varpi$ of a planet after a step of size $t$ is given by

[TABLE]

Compare this to simply calculating the Poisson bracket $\left\{A,\left\{A,B\right\}\right\}$ explicitly. This would give us the change in energy. But to calculate the change in any other quantity, we would need the extra steps of calculating the corresponding equations of motions and then solving them.

4.4 One planet and general relativistic precession

Let us start exploring how these tools might help us understand the longterm behaviour of symplectic integrators by looking at a very simple test problem consisting of just one planet orbiting a central mass together with general relativistic corrections in the form of an additional $1/r^{3}$ term in the potential (Nobili & Roxburgh, 1986). We treat the general relativistic corrections as a perturbation in part $B$ of the Hamiltonian. In this case they are the sole perturbation as there are no other planets. The dominant Keplerian motion of the planet is described by $A$. In the particular case we are looking at, Mercury and the Sun, the general relativistic corrections are small, leading to a precession rate of 0.5 arc-seconds per year (Park et al., 2017).

None of the symplectic methods described above will be able to solve this problem exactly because $\varphi^{[A]}$ and $\varphi^{[B]}$ don’t commute and, among others, the Poisson brackets $\{A,\{A,B\}\}$ and $\{B,\{A,B\}\}$ are non-zero. Let us focus on the standard WH method for now. In Fig. 3, we plot the one-step as well as the cumulative errors of an integration with the WH integrator and a timestep of about 10% of the planet’s orbit. Note that in contrast to the one-planet Kepler problem without GR, here neither the phase nor the periastron are exact action or angle variables, but combinations thereof. As one can see in the first panel, both the one-step and the cumulative energy errors oscillate but show no sign of any secular increase, as one might expect for a symplectic method. The same conclusion could be reached by looking at the periastron errors shown in the third panel, which do seem to average out, and no long-term drift is observed. However, one can see in the middle panel, the one step phase errors do not average out, leading to the absolute value of the cumulative phase error to increase linearly with time over long timescales. In a system such as the Solar System which evolves primarily due to secular interactions on long timescales, it is important to get the periastron resolved accurately. The phases on the other hand matter in systems which are dominated by mean motion resonances or encounters.

Typically, one would estimate the error of the method by looking at Eq. 10 and conclude that the leading order term is $\mathcal{O}(\epsilon dt^{2})$ coming from $\{A,\{A,B\}\}$. And indeed, this does provide a good estimate of the energy error. However, this is not the only metric we might want to use. For example, we might reasonably ask how the errors in the planet’s other orbital elements behave. As mentioned above, one way to estimate these errors would be to analytically calculate the evolution under $\{A,\{A,B\}\}$ and $\{B,\{A,B\}\}$. For this simple one planet test case this is straightforward to do, at least to leading order. But let us use the composition method introduced above which allows us to quickly estimate the errors due to the terms in Eq. 10 using any metric we want by using the already implemented operators $\varphi^{[A]}$ and $\varphi^{[B]}$.

We plot the estimate for the periastron error due to the two leading order term in Eq. 10 in Fig. 4. We calculate these estimates using the ${\mathcal{C}}^{AAB}$ and ${\mathcal{C}}^{BAB}$ operators. The top panel shows that the these estimates match the actual errors (Fig. 3) almost perfectly (the curves are hardly visible because they are on top of each other). The contributions from ${\mathcal{C}}^{BAB}$ are significantly smaller than those from ${\mathcal{C}}^{AAB}$. Note the different scale on the vertical axes, approximately one factor of $\epsilon$. Nevertheless, we can observe a qualitatively different behaviour for the errors coming from ${\mathcal{C}}^{BAB}$: these errors do not cancel out over one orbital period and lead to a secular growth in the cumulative error.

We can use these results to predict the accuracy of the WH integrator over very long timescales without any further integrations. Not only that, we can also correctly predict the long-term behaviour of the WH integrator with symplectic correctors and the higher order methods WHCKL and SABACL4 using only the results shown in Fig. 4.

Over short timescales, the secular error growth from $\{B,\{A,B\}\}$ is not important because it is so much smaller. However, if we integrate for long enough, this will eventually become the dominant contribution to the periastron error. Just by looking at our calculation over five orbital periods, we can estimate the growth rate and thus estimate when this happens. For this specific problem, the error due to $\{B,\{A,B\}\}$ will become dominant after $\approx 10^{6}$ orbital periods.

The above calculation was done for the standard WH integrator. However we can also use it to predict the long-term periastron error of other integrators. If we use WHC with symplectic correctors, we remove the $\{A,\{A,B\}\}$ term. This leads to a significant improvement of the integrations’s accuracy over short timescales because the only contribution to the periastron error is now the $\{B,\{A,B\}\}$ term which is several orders of magnitude smaller (compare the top and bottom panels of Fig. 4). But on long timescales, the $\{B,\{A,B\}\}$ term, which is still present in the shadow Hamiltonian even when using the WHC integrator, will grow in exactly the same way as it did for the WH integrator without correctors. Thus, over long timescales, we do not expect any benefit from using symplectic correctors in reducing the periastron error. We test this by running the same simulation for $10^{7}$ orbits and comparing the periastron error for the WH and WHC integrators to an IAS15 simulation. The results are shown in Fig. 5. We can see that the correctors indeed improve the periastron error significantly, but only on short timescales. On long timescales, $\gtrapprox 10^{6}$ orbital periods, the errors are identical for the WH and WHC integrators. Given that Mercury’s orbital period is only 88 days, $10^{6}$ orbits correspond to about 250 kyrs, i.e. a short timescale compared to both the lifetime of the system and the timescale on which the dynamics of the Solar System change significantly (a few million years).

We also show the periastron error for the WHCKL integrator in Fig. 5. This integrator does not have the $\{B,\{A,B\}\}$ term in the shadow Hamiltonian and we therefore expect it to perform significantly better. A linear growth of the periastron error is still visible at late times, but this is now suppressed by two orders of magnitude. Given that the now leading error term is two orders of magnitude smaller and the timestep is 10% of the orbital period, we can deduce that the now leading order term must be of order $\mathcal{O}\left(\epsilon^{2}dt^{4}\right)$.

4.5 Two planets

We now repeat the line of argument from above, but this time with two planets and without general relativistic corrections. The two planets correspond to Mercury and Saturn. Saturn contributes about $7.3^{\prime\prime}$ per century to the precession rate of Mercury (Park et al., 2017). Now part $A$ of the Hamiltonian corresponds to the Keplerian motion of the planets and part $B$ corresponds to planet-planet interactions.

The same one step and cumulative errors of the periastron as before are shown in Fig. 6. And as before, we can see that the contributions to the periastron error due to $\{A,\{A,B\}\}$ average out to leading order. The contributions due to $\{B,\{A,B\}\}$ on the other hand do not average out and we can see a secular growth, clearly visible after only 60 Mercury periods. Thus, even though the errors due to $\{B,\{A,B\}\}$ are significantly smaller than those due to $\{A,\{A,B\}\}$, we expect that they become dominant after $\approx 10^{5}$ Mercury periods, i.e. $25$ kyrs.

Note that once again, we were able to predict when which error dominates without actually running a long-term integration. To verify this, we run the two-planet case for $10^{6}$ Mercury period ($250$ kyr) and plot the periastron error in Fig. 7 for simulations using the WH, WHC, and WHCKL integrators. The results confirm our estimates. On short timescales, the symplectic correctors significantly improve the periastron error. However, after about a hundred thousand years, the simulations using correctors have the same error in the periastron as those without. We conclude that symplectic correctors do not improve the accuracy of secularly evolving planetary systems over long timescales.

We also show the periastron error for the WHCKL integrator in Fig. 7. As before, since this integrator does not have the $\{B,\{A,B\}\}$ term in the shadow Hamiltonian, it performs significantly better. Periastron errors are still present due to higher order terms, but they are significantly suppressed.

4.6 Implications for long-term integrations of the Solar System

We can now come back to explaining the results of the long-term integration of the entire Solar System presented in Sect. 3. As compared to the simple test problems shown above, the full Solar System is a much more complex dynamical system. Nevertheless, we can now qualitatively understand why the WHC integrator performs so poorly in reproducing the secular frequencies compared to higher order methods despite its good energy conservation properties.

This behaviour is simply a result of the $\{B,\{A,B\}\}$ term in the shadow Hamiltonian. Although this term leads to only oscillatory behaviour in the semi-major axis error of a planet as well as the global energy error, it leads to a secular drift in the periastron error. If we use a timestep of $8$ days, then after only $\approx 25$ kyrs, the error in the periastron of Mercury is dominated by the cumulative contributions from the $\{B,\{A,B\}\}$ term for both the WH and WHC integrator. And since the periastron is used to determine the secular frequencies, we find no improvement in the secular frequencies despite the fact that other quantities such as the semi-major axes are more precise.

The integrators WHCKL and SABACL4 on the other hand get rid of the $\{B,\{A,B\}\}$ term and can therefore achieve a significantly better accuracy of the secular frequencies.

As one can see from Fig.2, the secular frequencies, and therefore the dynamics, in the outer Solar System can be improved by using symplectic correctors (bottom row). Specifically, the secular frequencies associated with Jupiter and Saturn, $g_{5}$ and $g_{6}$, improve by about one order of magnitude. Note that this is still significantly less than the improvement in the semi-major axes (top row) which is of order $\epsilon\sim 10^{-3}$. The frequencies $g_{7}$ and $g_{8}$ improve even more. We attribute this difference in the convergence of secular frequencies in the inner and outer Solar System to two effects. First, the orbital timescales are longer in the outer Solar System. Thus, it will take longer for the error contributions from the $\{B,\{A,B\}\}$ term to become dominant for the outer planets. After all, over short timescales, the correctors lead to an improvement (see Fig. 7) and we here only integrate for 20 Myrs. Second, the dynamics in the outer Solar System are also physically different from the dynamics in the inner Solar System. For example, the semi-major axes of the outer planets change significantly ($\Delta a/a\sim\epsilon^{2}\sim 10^{-6}$) over the timespan of our integrations (Milani et al., 1987). Since symplectic correctors improve the semi-major axes of all planets, we can expect them to also somewhat improve the secular frequencies in the outer Solar System if we determine them to such extreme precision as we do. Note that this shows that secular frequencies are indeed a very useful metric to use in a convergence study given that they can capture these effects.

There is a simple physical explanation as to why removing the $\{B,\{A,B\}\}$ term is important for the accurate determination of secular frequencies, but removing the $\{A,\{A,B\}\}$ term has little effect according to our results. First note that the secular precession frequencies are to first order functions of the planets’ semi-major axes and masses only (see e.g. Exercise 7.1 in Murray & Dermott, 2000). However, the semi-major axes of the planets are constant to within $\mathcal{O}\left(\epsilon\right)$ in a secularly evolving system as already shown by Laplace (1776). Because the $\{A,\{A,B\}\}$ term is $\mathcal{O}\left(\epsilon\right)$, its effects on the semi-major axes average out (at least to leading order in $\epsilon$). Thus, this term does also not affect the secular frequencies over long timescales (again at least to leading order in $\epsilon$).

5 Conclusions

In this paper, we have shown that the energy error alone is not a sufficient metric when measuring the accuracy of long-term planetary systems which evolve due to secular interactions. Specifically, we have presented evidence showing that the Wisdom-Holman integrator with symplectic correctors does not improve the accuracy of secular frequencies in the inner Solar System compared to the WH integrator without symplectic correctors. This is despite the fact that the symplectic correctors improve the energy error by three orders of magnitude.

We have presented a framework that helps us understand this behaviour. We are now able to trace the origin of the poor performance to a specific term in the shadow Hamiltonian. If we use a higher-order symplectic integrator such as the Kernel method of Wisdom et al. (1996) or the SABACL4 method of Laskar & Robutel (2001) which get rid of this term, then the integrators’ performance is improved significantly.

The use of composition operators to study a specific term in the shadow Hamiltonian has multiple advantages over an order of magnitude estimate of Poisson brackets. First, it lets us study terms that grow with time but are initially small and obscured by much larger terms that might not grow with time. Second, we can study the effect of any error term on any arbitrary quantity in a simulation, i.e. not just the energy error. This allows us to study cases where the error in the energy is oscillatory, but the error in other quantities such as the periastron error is not, leading to an increasing error with time. The composition method is very simple to implement as it only uses already implemented functions which are part of the standard Wisdom-Holman integrator and is compatible with any position dependent potential.

We have further shown that long-term simulations of the Solar System can achieve an extremely high level of agreement in the secular frequencies across many different integrators. This makes secular frequencies a great tool to study the accuracy of symplectic integrators. If we use higher order methods which get rid of the dominant error terms for the periastron precession of planets, then we can achieve a remarkable level of precision. With a timestep of $10$ days, corresponding to only 8 timesteps per shortest dynamical timescale (Mercury’s orbital period), we can calculate the secular frequency associated with Mercury to a precision of $10^{-9}$ arc-seconds per year! However, it is important to reiterate that the purpose of this work was not an attempt to produce a list of the most accurate secular frequencies for the Solar System, but to study the performance of various symplectic integrators. For the former, we refer the reader to Laskar et al. (2011) who use a much more realistic model of the Solar System than we do. Further, note that simply determining the secular frequencies to extremely high precision does not lead to new physical insight about the system. In fact, the frequencies are not well defined in secular theory beyond approximately 1 part in $10^{2}$ for the inner planets and 1 part in $10^{4}$ for the outer planets.

Providing evidence that numerical methods converge and that physical results do not depend on nuisance parameters such as the timestep or the specific numerical integrator used are crucial for any numerical experiment. Most studies that check the convergence of $N$-body simulations monitor integrals of motion such as the energy. We argue that it is wise to go beyond this and check that the quantities which are physically relevant for the evolution of the system are converged as well. In the case of secularly evolving systems such as the Solar System, the secular frequencies are a natural choice. As we have shown, having an energy error which is dominated by round-off, does not necessarily imply that a simulation is converged. Reducing the timestep even further can continue to improve the accuracy, not of the energy error of course, but of other quantities.

Finally, let us comment on how chaos fits into this picture. It is well known that the Solar System is chaotic with the shortest Lyapunov timescale being about 5 Myr (Laskar, 1989). Let us consider two of our 20 Myr simulations with two different high order integrators. If the timestep is very small, say $dt=0.1$ days, then the simulations are converged to machine precision, $10^{-16}$. Even if the these simulations use the exact same model and the exact same initial conditions, they will differ after only one timestep because the two integrators perform slightly different operations and this leads to a round-off error due to finite floating point precision. Let us assume that the two simulations differ by $10^{-16}$ after one timestep. We can then think of one of these simulations as a shadow simulation of the other, using slightly perturbed initial conditions. If the system is chaotic, this difference will grow exponentially. Our integration span four $e$-folding times, thus the initial difference of $10^{-16}$ will increase to about $5\cdot 10^{-15}$ by the end of the simulation. But since we introduce round-off errors at every timestep (not just in the first one) the round-off errors grow with time as well. By the end of the simulation the round-off error is $\approx 10^{-10}$, see Eq. 15 and Fig. 1. We obtain the secular frequencies to a precision of $10^{-9}$. Thus, the frequencies are clearly not affected by chaos. The precision is either limited by round-off errors for small timesteps, or discretization errors coming from the integrator for large timesteps as shown in Fig. 2. But this is only true for this specific simulation. If we integrate for much longer, then the precision to which we can obtain the secular frequencies will eventually be limited by chaos. To summarize, there are four different contributions to the divergence of trajectories in our simulations with symplectic integrators:

Discretization errors in the form of oscillating error terms which average out over time and lead to a constant error, $\mathcal{O}\left(1\right)$. 2. 2.

Discretization errors where the errors do not average out, growing linearly with time, $\mathcal{O}\left(t\right)$. 3. 3.

Floating point round-off errors, growing like a random walk in predominantly action-like variables in an unbiased implementation, $\mathcal{O}\left(\sqrt{t}\right)$, and growing linearly in a biased implementation, $\mathcal{O}\left(t\right)$. For predominantly angle-like variables this can also lead to errors growing as $\mathcal{O}\left(t^{1.5}\right)$ or $\mathcal{O}\left(t^{2}\right)$, respectively. 4. 4.

Chaos, with the seed coming from round-off errors, leading to an exponential divergence, $\mathcal{O}\left(e^{t}\right)$.

Which contribution dominates depends on the system at hand, the integrator and timestep used, and last but not least the length of the integration because the various contributions depend very differently on time, i.e. $\mathcal{O}\left(1\right)$, $\mathcal{O}\left(\sqrt{t}\right)$, $\mathcal{O}\left(t\right)$, $\mathcal{O}\left(t^{1.5}\right)$, $\mathcal{O}(t^{2})$, $\mathcal{O}\left(e^{t}\right)$. The tools described in this paper help to determine which of these contributions is limiting the accuracy of trajectories in a specific $N$-body simulation, and in some cases, how to improve the accuracy.

Acknowledgments

We thank an anonymous referee for helpful comments. This research has been supported by the NSERC Discovery Grant RGPIN-2014-04553 and the Centre for Planetary Sciences at the University of Toronto Scarborough. Support for this work was provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51423.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. This research was made possible by the open-source projects Jupyter (Kluyver et al., 2016), iPython (Pérez & Granger, 2007), and matplotlib (Hunter, 2007; Droettboom et al., 2016).

Appendix A Taking IAS15 as the true solution

Let us comment on why Fig. 2 provides further evidence that the IAS15 simulation is indeed much more accurate than any of the other simulation in our sample. This is not to say that the IAS15 algorithm is superior. The difference is partly due to the fact that IAS15 uses compensated summation to reduce round-off errors whereas the other integrators do not have this capability in our implementation (see also Laskar & Robutel 2001 and Wisdom 2018 who use these algorithms with some form of extended precision).

Note that all error metrics for all integrators in Fig. 2 are eventually dominated by round-off errors at small timesteps. Specifically, note that the curves are noisy as expected from a random walk. Because we are comparing all simulations shown in the plot to the same IAS15 simulation, this noisy pattern would not show up if the IAS15 simulation were less accurate than the WHFast simulations, or if the round-off error came from the IAS15 simulation. If the IAS15 simulation had a larger error than the other simulations, then all the curves for the different integrators would be on top of each other on the left side of Fig. 1 and 2.

This gives us further confidence that our assumption of using the IAS15 simulation as the true solution is appropriate. Alternatively, we could have used one of the higher order symplectic integrations with the smallest timestep as the reference simulation. This would still allow us to determine when the error is dominated by round-off errors. But we would not be able to determine the power law in the round-off error dominated regime.