High-order symplectic integrators for planetary dynamics and their implementation in rebound

ABSTRACT

Direct N-body simulations and symplectic integrators are effective tools to study the long-term evolution of planetary systems. The Wisdom–Holman (WH) integrator in particular has been used extensively in planetary dynamics as it allows for large time-steps at good accuracy. One can extend the WH method to achieve even higher accuracy using several different approaches. In this paper, we survey integrators developed by Wisdom et al., Laskar & Robutel, and Blanes et al. Since some of these methods are harder to implement and not as readily available to astronomers compared to the standard WH method, they are not used as often. This is somewhat unfortunate given that in typical simulations it is possible to improve the accuracy by up to six orders of magnitude (!) compared to the standard WH method without the need for any additional force evaluations. To change this, we implement a variety of high-order symplectic methods in the freely available N-body integrator rebound. In this paper, we catalogue these methods, discuss their differences, describe their error scalings, and benchmark their speed using our implementations.

gravitation, methods: numerical, planets and satellites: dynamical evolution and stability

1 INTRODUCTION

Astronomers have predicted the location of planets in the night sky since ancient times. However, it has only recently become possible to calculate the orbital evolution of planetary systems accurately over very long time-scales with the help of fast computers (Laskar & Gastineau 2009). Various different numerical algorithms, so-called integrators, have been used for that purpose. Because of the large separation of time-scales involved, from a day to billions of years, specialized integrators are needed to predict the orbital evolution of planetary systems over their entire lifetime that can correspond to up to 10¹² orbits.

There are many options when it comes to choosing an integrator. Some generic integrators achieve a good accuracy by being very high order. Examples include the Bulirsch–Stoer (Stoer & Bulirsch 2002) or the Gauß–Radau based IAS15 integrator (Rein & Spiegel 2015). A different approach is that of Wisdom & Holman (1991), which makes use of the fact that we are interested in the evolution of a dynamical system rather than simply the solution of a generic differential equation, and that the gravitational interactions between planets can be considered perturbations to otherwise Keplerian orbits. Their integrator, which we refer to the classical Wisdom–Holman integrator (WH) is therefore particularly well suited for planetary systems where orbits remain well separated. This is the case for many planetary systems including the Solar system. Different implementations of the WH integrator and various extensions of it (e.g. allowing close encounters) are freely available and have been used extensively by the astrophysics community (Chambers & Migliorini 1997; Duncan, Levison & Lee 1998; Rein & Tamayo 2015; Rein et al. 2019).

The classical WH integrator is a second-order method. If one requires greater accuracy one has two options: either reducing the time-step or choosing a different higher order method (which ideally also makes use of the fact that planet–planet interactions are small). Many simulations require an astronomically large number of time-steps to begin with¹ and it is not uncommon for simulations to run for a month or even a year in some cases (Laskar et al. 2011). Since reducing the time-step, and thus increasing the number of time-steps, leads to a greater computational cost there might be a significant advantage if one can use a higher order method that offers better accuracy at a fixed time-step

The higher order integrators that we are looking at in this paper all assume that planetary systems do not have close encounters and planetary orbits remain well separated. There are once again different approaches for obtaining such a high-order integrator. Each comes with many subtle choices that one can make along the way. In this paper, we focus on the advances driven by two groups, one led by Jack Wisdom (Wisdom, Holman & Touma 1996; Wisdom 2006, 2018) and the other by Jacques Laskar (Laskar & Robutel 2001; Blanes et al. 2013). Since their implementations are not easily available for the community, we have implemented their integrators within the rebound package, together with some new variants. We review the different approaches that these groups take in Section 2 and compare their integrators’ formal properties. We present our own implementations in Section 3. We then verify the accuracy of our implementations and measure their speed in numerical tests that we present in Section 4. We hope that making these integrators easily available within rebound will allow more people who are interested in planetary dynamics to use and build upon them.

2 SYMPLECTIC INTEGRATORS

2.1 Splitting methods

All integrators described in this paper are symplectic splitting schemes. In this section, we provide a short introduction to such methods. This is by no means a comprehensive discussion of this vast subject area. Our goal is not to rederive these methods but to introduce the notations that we will use later. We opt for a formal but nevertheless hopefully easy to understand approach for those with little background in this area. For different approaches and perspectives we refer the reader to papers and books by Laskar & Robutel (2001), Wisdom (2018), Hairer, Lubich & Wanner (2006) as well as references therein.

The differential equations of the N-body problem have special properties because they can be derived from a Hamiltonian H together with Hamilton’s equations. From a mathematical point of view, such differential equations have more structure than an arbitrary one. Many of these mathematical properties have direct physical consequences, like the conservation of energy and phase space volume. The integrators that we describe in this paper preserve some of these structures by construction, and one can argue that they are therefore preferable. We note, however, that even though we focus on conservative systems here, many of the splitting methods that we describe can also be applied to non-Hamiltonian systems (Tamayo et al. 2019).

Let us first define the Poisson bracket {g, h} of two functions g and h of the canonical coordinates (q_i, p_i) as

$$\begin{eqnarray*} \lbrace g,h\rbrace = \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). \end{eqnarray*}$$

(1)

The functions g and h are arbitrary, they can, for example, be a single coordinate, the energy, or the angular momentum. We can use the Poisson bracket to write down Hamilton’s equations as

$$\begin{eqnarray*} \dot{q}_i = \frac{\partial H}{\partial p_i} = \lbrace q_i,H\rbrace \quad \quad \text{and}\quad \quad \dot{p}_i = -\frac{\partial H}{\partial q_i} = \lbrace p_i,H\rbrace . \end{eqnarray*}$$

(2)

Thus, through the chain rule, the time derivative of any function g that depends only on p and q can be written succinctly as

$$\begin{eqnarray*} \dot{g} = \lbrace g,H\rbrace . \end{eqnarray*}$$

(3)

Assume that the planetary system we are interested in has 3N coordinates and 3N momenta. If we define a function y(t) = (q₁(t), …, q_3N(t), p₁(t), …, p_3N(t)), then solving the N-body problem corresponds to solving the 3N coupled differential equations |$\dot{y} = \lbrace y,H\rbrace$|⁠. Let us further define the Lie derivative with respect to the Hamiltonian H as

$$\begin{eqnarray*} {\mathcal {L}}_H y = \lbrace y,H\rbrace , \end{eqnarray*}$$

(4)

or, just writing down the operator without the operand,

$$\begin{eqnarray*} {\mathcal {L}}_H = \lbrace \cdot , H\rbrace . \end{eqnarray*}$$

(5)

Comparing with equation (3), we see that the Lie derivative with respect to the Hamiltonian H is an operator that yields the time derivative of any function of the phase space coordinates that it acts on. Finally, let us define the formal solution operator |$\varphi ^{[H]}_t(y_0)$|⁠. This operator returns the solution to the differential equation |$\dot{y} = \lbrace y,H\rbrace$| at time t with initial conditions y₀ given at t = 0.

In general we do of course not know how to write down this abstract solution operator. However, with the notation that we have introduced, we can now at least write down a formal expression. In terms of the Lie derivative defined above and an exponential defined in the usual way as an infinite series, a formal expression of the solution operator is

$$\begin{eqnarray*} {\varphi }^{[H]}_{t}(y_0) &=& \exp \left(t {\mathcal {L}}_{H} \right) \text{Id}(y_0) \nonumber \\ &=& \left(\text{Id} + t {\mathcal {L}}_H \text{Id}+\frac{1}{2} {\mathcal {L}}_H {\mathcal {L}}_H \text{Id} +\cdots \right) (y_0). \end{eqnarray*}$$

(6)

Note that the exponential in equation (6) is simply a formal way to write down an operator. It is acting on the function to its right, in the above case this is the identity Id. To see that this construct indeed provides a solution to the differential equation |$\dot{y} = \lbrace y,H\rbrace$| with initial conditions y₀, note that the series expansion of the exponential results in a series of differential operators acting on Id(y₀). This series is simply the Taylor expansion of the solution φ around t = 0.

Since it is not possible to write down an explicit solution for the N-body problem, we turn to operator splitting methods. One way to view an operator splitting method is to take a differential equation of the form

$$\begin{eqnarray*} \dot{y} = {\mathcal {L}}_H (y) \end{eqnarray*}$$

(7)

and rewrite it as

$$\begin{eqnarray*} \dot{y} = \left({\mathcal {L}}_A + {\mathcal {L}}_B\right)(y) = {\mathcal {L}}_A(y)+ {\mathcal {L}}_B(y) . \end{eqnarray*}$$

(8)

Note that due to the linearity of the Poisson bracket, we have |${\mathcal {L}}_{A+B} = {\mathcal {L}}_A + {\mathcal {L}}_B$|⁠. This implies that for Hamiltonian systems we can split the equations of motion by splitting the Hamiltonian H into two parts, i.e. H = A + B. Depending on the splitting, we get two new sets of differential equations that we can study independently:

$$\begin{eqnarray*} \dot{y} = {\mathcal {L}}_A (y) \quad \text{and} \quad \dot{y} = {\mathcal {L}}_B (y) . \end{eqnarray*}$$

(9)

In practice, we will choose a splitting that allows us to easily solve these new differential equations independently. In other words, we want to be able to explicitly write down the solution operators |${\varphi }^{[A]}_{t}(y_0)$| and |${\varphi }^{[B]}_{t}(y_0)$|⁠.

Using the notation from above, we have

$$\begin{eqnarray*} \exp \left(t {\mathcal {L}}_{H} \right) = \exp \left(t ({\mathcal {L}}_{A}+{\mathcal {L}}_B) \right), \end{eqnarray*}$$

(10)

but because these exponentials are operators, not scalars, we unfortunately cannot use standard rules for exponentials. Specifically, we typically have

$$\begin{eqnarray*} \exp \left(t ({\mathcal {L}}_{A}+{\mathcal {L}}_B) \right) \ne \exp \left(t {\mathcal {L}}_{A} \right) \exp \left(t {\mathcal {L}}_{B} \right). \end{eqnarray*}$$

(11)

For illustration purposes, imagine the above equation were true. We could then solve the differential equation in equation (7) by first applying the solution operator φ^[B] to our initial conditions for a time t (remember, operators act to the right), then take the result and apply the solution operator φ^[A] on it, again for time t. The result of the last operation would then correspond to the solution φ^[H] at time t. But because the equation above is not true, we need to rely on the Baker–Campbell–Hausdorff (BCH) identity,² to get a similar expression

$$\begin{eqnarray*} &&{\exp \left(t {\mathcal {L}}_{A} \right) \exp \left(t {\mathcal {L}}_{B} \right)} \nonumber \\ &&{= \exp \left(t {\mathcal {L}}_{A} + t {\mathcal {L}}_{B} + \frac{1}{2} t^2 \left[{\mathcal {L}}_A,{\mathcal {L}}_B \right] + \mathcal {O}\left(t^3\right)\right), } \end{eqnarray*}$$

(12)

where |${}[{\mathcal {L}}_A, {\mathcal {L}}_B] ={\mathcal {L}}_A {\mathcal {L}}_B - {\mathcal {L}}_B {\mathcal {L}}_A$| is a commutator. So if we apply the splitting procedure described above (corresponding to the left-hand side of equation 12), we do not solve the differential equation |$\dot{y} = {\mathcal {L}}_H (y)$|⁠, but some other differential equation (given by the argument in the exponential on the right-hand side of equation 12). Formally, we can write down the Hamiltonian corresponding to the differential equation that we are actually solving as

$$\begin{eqnarray*} H^{\prime } = A+B+\frac{1}{2} t \lbrace A,B\rbrace + \mathcal {O}\left(t^2\right), \end{eqnarray*}$$

(13)

where we have used

$$\begin{eqnarray*} \left[ {\mathcal {L}}_{A} , {\mathcal {L}}_{B} \right] = {\mathcal {L}}_{\lbrace A,B\rbrace }, \end{eqnarray*}$$

(14)

which follows from the Jacobi identity. Although we do not end up solving the Hamiltonian system we were trying to solve, the fact that we can write down a different Hamiltonian system that we seem to solve exactly is often quoted as why symplectic integrators have superior properties.³ We can study the properties of this ‘nearby’ Hamiltonian to see how far our approximation is from the true solution. As we can see from the equation above, for small t, the Hamiltonian is close to the original Hamiltonian and only differs by |$\frac{1}{2} t \lbrace A,B\rbrace + \mathcal {O}\left(t^2\right)$|⁠. The smaller a step we make, the closer we stay to the true Hamiltonian. We indicate that we take a small but finite step by replacing t with dt. To reach a time t we simply repeat the procedure n times such that |$t=\mathrm{ d}t\, n$|⁠.

2.2 Wisdom–Holman integrator

Wisdom & Holman (1991) split the N-body Hamiltonian H of a planetary system into two parts A and B such that A contains all terms describing the planets’ Keplerian motion and B contains the planet–planet interactions. Jacobi coordinates are particularly well suited for this. This coordinate choice has a further advantage, it keeps the momentum terms in the Hamiltonian in the form |$\sum _i p_i^2/(2m_i)$|⁠. For an analysis of other splittings, see Hernandez & Dehnen (2017) and Rein & Tamayo (2019).

Note that the Hamiltonian B only contains a position-dependent potential. We can easily calculate the corresponding equations of motion, |$\dot{y} = {\mathcal {L}}_B y$|⁠, and solve them exactly. The solution operator φ^[B] is simply a kick, i.e. a change of the particles’ momenta while keeping their positions fixed. Although it is not quite as straightforward, the solution φ^[A] to the equations of motion |$\dot{y}= {\mathcal {L}}_A y$| can also be computed relatively easily if one can efficiently solve Kepler’s equation (the solution is only required to within finite machine precision, see e.g. Rein & Tamayo 2015).

The motion of well-separated planets is described almost completely by only the Keplerian motion, A, and we can consider the interactions due to other planets, B, a perturbation. We can quantify this with a small number ϵ and keep track of it by formally replacing B in the Hamiltonian with ϵB such that the splitting is of the form H = A + ϵB.

The splitting scheme described in the previous section is a first-order scheme, the WH integrator uses a second-order one. To advance the solution by one time-step dt, it first applies the solution operator |${\varphi }_{\mathrm{ d}t/2}^{[A]}$| for half a time-step, then applies |${\varphi }_{\mathrm{ d}t}^{[B]}$| for a full time-step, followed by another operation of |${\varphi }_{\mathrm{ d}t/2}^{[A]}$|⁠. Let us formally write one time-step of this scheme as

$$\begin{eqnarray*} \text{WH} = e^{\frac{1}{2} A}e^{ B}e^{\frac{1}{2} A}. \end{eqnarray*}$$

(15)

As before, we can apply the BCH formula (twice in this case) to calculate the nearby Hamiltonian, i.e. the Hamiltonian of the system the integrator is actually solving (sometimes this is also referred to as the ‘shadow’ Hamiltonian), to be

$$\begin{eqnarray*} H^{\prime } &=& A + \epsilon B \nonumber - \epsilon \frac{\mathrm{ d}t^2}{24} \left\lbrace A,\left\lbrace A, B \right\rbrace \right\rbrace - \epsilon ^2 \frac{\mathrm{ d}t^2}{12} \left\lbrace B,\left\lbrace A, B \right\rbrace \right\rbrace \\ && + \mathcal {O}(\epsilon \mathrm{ d}t^4) + \mathcal {O}(\epsilon ^2 \mathrm{ d}t^4) + \cdots \end{eqnarray*}$$

(16)

Note that the first error terms appear at second order in dt. There are two terms, one is |$\mathcal {O}\left(\epsilon \mathrm{ d}t^2\right)$|⁠, the other is |$\mathcal {O}\left(\epsilon ^2 \mathrm{ d}t^2\right)$|⁠. For small ϵ, the first term dominates. Compared to other methods, for example the standard leap-frog integrator,⁴ the thing to note here is the factor of ϵ. For the Solar system, ϵ ∼ 10⁻³, making the WH integrator about 1000 times more accurate than the leap-frog integrator whose error terms have the same structure as those above, but do not include the factors of ϵ.

We refer to this integrator as the classical WH method.

2.3 Wisdom–Holman–Touma integrator family

2.3.1 Symplectic correctors

As we have seen above, the classical WH method has two error terms at second order in time, |$\mathcal {O}\left(\epsilon \mathrm{ d}t^2\right)$| and |$\mathcal {O}\left(\epsilon ^2 \mathrm{ d}t^2\right)$|⁠. To obtain a higher accuracy method, we would like to get rid of the leading terms. Instead of simply using a higher order splitting scheme (Yoshida 1993), Wisdom et al. (1996) describe a different approach where in their derivation they make use of the fact that the system we are integrating is a Hamiltonian system.

Wisdom et al. (1996) apply a corrector e^C step before and after each step to remove the leading term, |$\mathcal {O}\left(\epsilon \mathrm{ d}t^2\right)$|⁠. Specifically the method is

$$\begin{eqnarray*} \text{WHC} = e^{ C} e^{\frac{1}{2} A}e^{ B}e^{\frac{1}{2} A} e^{-C}. \end{eqnarray*}$$

(17)

Because e^−C is the inverse of e^C, these correctors only need to be applied at the beginning of the integration and at the end, or whenever an output is generated. But they do not need to be applied in-between time-steps since they cancel out. For long integrations this method is therefore effectively as fast as the classical WH method.

Deriving the e^C operator would go beyond the scope of this paper and we refer the reader to Wisdom et al. (1996) and Wisdom (2018). To summarize, the leading order error term (at order ϵ) in the classical WH method can be interpreted as arising from a mismatch in initial conditions between the real problem and the modified problem that the splitting scheme is solving. Applying the symplectic correctors corresponds to a canonical transformation to and from integrator coordinates that removes the |$\mathcal {O}(\epsilon)$| term.⁵ It is in general not possible to derive an exact expression for the corrector e^C, so it has to be approximated. We can do this approximation to any order we want and do not need to make it particularly efficient because the correctors are only applied at the beginning and end of the simulation, not during intermediate time-steps

Wisdom (2006) gives explicit expressions for e^C to various order dt^p in terms of the already implemented Keplerian motion and kick operators, e^A and e^B. After applying the correctors, the largest error terms are |$\mathcal {O}\left(\epsilon \mathrm{ d}t^{p+1}\right)$| and |$\mathcal {O}\left(\epsilon ^2 \mathrm{ d}t^2\right)$|⁠. We refer to this integrator as WHCp, where p is the order of the corrector. If we apply a high enough order corrector, for example a 17th order one, then |$\mathcal {O}\left(\epsilon \mathrm{ d}t^{18}\right)$| can be ignored and the dominant error term is |$\mathcal {O}\left(\epsilon ^2 \mathrm{ d}t^2\right)$| for all reasonable time-steps.

2.3.2 Kernel method

Using the 17th order correctors of Wisdom (2006), the leading order term of the WHC method is, for all practical purposes, |$\mathcal {O}\left(\epsilon ^2 \mathrm{ d}t^2\right)$|⁠. Naturally, we would like to get rid of this error term as well. The approach taken by Wisdom et al. (1996) is to change the kick step of the WH method, e^B. If we refer to |$e^{K}=e^{\frac{1}{2} A}e^{ B}e^{\frac{1}{2} A}$| as the kernel of the WH and WHCp methods and change it |$e^{\frac{1}{2} A}e^{ B^{\prime }}e^{\frac{1}{2} A}$|⁠, then with the right choice of B′ and one extra corrector step, one can, at least in theory, eliminate all |$\mathcal {O}\left(\epsilon ^2\right)$| terms. We refer to this ideal integrator as WHCCKI (WH + first corrector + second corrector + ideal kernel). If it were possible to implement this kernel (as well as the exact correctors), its leading error term would be |$\mathcal {O}\left(\epsilon ^3\mathrm{ d}t^4\right)$|⁠.

Unfortunately, the ideal B′ is another infinite series of Poisson brackets and we cannot (at least not for the N-body problem) hope to solve the evolution of the system under it exactly. Thus, we are once again required to rely on an approximation. However, because B′ itself is a power series in dt, we only need to match B′ to a finite power in dt to get rid of the leading order |$\mathcal {O}\left(\epsilon ^2 \mathrm{ d}t^2\right)$| term. Note that by doing so we still keep higher order terms involving ϵ², including the now leading term |$\mathcal {O}\left(\epsilon ^2 \mathrm{ d}t^4\right)$|⁠. Wisdom et al. (1996) describes three different methods for a practical implementation of such an approximation.

The first method involves a composition of e^A and e^B operators. One example of a kernel that has the right properties to eliminate the leading order term is

$$\begin{eqnarray*} e^K = e^{\frac{5}{8} A} e^{-\frac{1}{6} B} e^{-\frac{1}{4} A} e^{\frac{1}{6} B} e^{\frac{1}{8} A} e^{B} e^{-\frac{1}{8} A} e^{-\frac{1}{6} B} e^{\frac{1}{4} A} e^{\frac{1}{6} B} e^{\frac{3}{8} A}. \end{eqnarray*}$$

Note that the implementation of this kernel is straightforward if we already have the e^A and e^B operators implemented. However, in contrast to the symplectic correctors, the kernel operator needs to be applied at every time-step and its efficiency is therefore important. The force evaluation is typically the slowest part of an integrator because it scales as |$\mathcal {O}\left(N^2\right)$| for N particles, whereas all other parts scale as |$\mathcal {O}\left(N\right)$|⁠. The above kernel contains five e^B operators and thus five force evaluations per time-step. We can therefore expect this method to be approximately five times slower than the WH and WHCp methods. We will see below if the increase in accuracy outweighs the decrease of performance. We will refer to this method as WHCCKC (WH + first corrector + second corrector + composition kernel). Because the composition kernel is only accurate to order |$\mathcal {O}\left(\epsilon ^2 \mathrm{ d}t^4\right)$|⁠, we might ignore the second corrector that only acts at higher orders. We will refer to this method with only the first corrector applied as WHCKC (WH + first corrector + composition kernel).

The second method to approximate B′ is to calculate the leading order terms explicitly. For the N-body problem, the leading order term, B, corresponds to the standard kick, involving derivatives of the planet–planet potentials and some additional terms due to Jacobi coordinates (Rein & Tamayo 2015). The second term in B′ is a Poisson bracket of the form {B, {B, A}}. For a WH type splitting of the N-body Hamiltonian in Jacobi coordinates, B is only dependent on positions, and {B, {B, A}} evaluates to a combination of spatial derivatives of the accelerations (thus second derivatives of the potential) and leads to expressions similar in structure to those describing the jerk.⁶ Thus, it is possible to implement an efficient operator corresponding to the evolution under the Hamiltonian

$$\begin{eqnarray*} \widehat{B} = B + \mathrm{ d}t^2 \frac{1}{24} \lbrace B,\lbrace B,A\rbrace \rbrace . \end{eqnarray*}$$

(18)

In fact, it is possible to calculate |$e^{\widehat{B}}$| (the ‘modified kick’) without any additional square root evaluations beside those already needed to calculate e^B (Wisdom 2018). Thus, if we use the kernel

$$\begin{eqnarray*} e^K = e^{\frac{1}{2} A} e^{\widehat{B}} e^{\frac{1}{2} A}, \end{eqnarray*}$$

(19)

we can expect a performance almost identical to that of the standard WH method. We refer to method that uses this approximation of the kernel as WHCKM (WH + first corrector + modified kick kernel). If we decide to use the second correctors as well, the method becomes WHCCKM. As above, the use of second correctors is not expected to improve the accuracy because we already introduce error terms at order |$\mathcal {O}\left(\epsilon ^2 \mathrm{ d}t^4\right)$| by truncating the in principle infinite series of B′ to only two terms to get |$\widehat{B}$|⁠.

The third method to approximate B′ is referred to as the ‘lazy implementer’s method’ by Wisdom et al. (1996). Instead of calculating the Poisson bracket {B, {B, A}} explicitly, it uses a Taylor series approximation to estimate the modified kick by first calculating the acceleration a_j due to e^B (i.e. the unmodified kick), then shifting the positions with |$q_j\rightarrow q_j+\frac{\mathrm{ d}t^2}{12}a_j$|⁠, before calculating new accelerations |$a^{\prime }_j$| with the updated positions, and finally updating the momenta with |$p_j\rightarrow p_j+\mathrm{ d}t\, m_j\, a^{\prime }_j$|⁠. The intermediate position values are discarded. This trick requires two force evaluations for the kernel and the method is thus approximately two times slower than the standard WH method. Note that the Taylor series introduces a new error term at order |$\mathcal {O}\left(\epsilon ^3 \mathrm{ d}t^3\right)$|⁠. Depending on the relative size of ϵ compared to dt, this may or may not be the new dominant term. In either case, because the kick step is only updating the momenta, the scheme itself remains symplectic. Note that compared to the modified kick method, the lazy implementer’s method requires no extra work when it is used with any additional position dependent potential, i.e. due to general relativity⁷ or tides. We refer to this method as WHCKL (WH + first corrector + lazy implementer’s kernel). Although the use of the second corrector might not be of any advantage here either, for completeness, we refer to it as WHCCKL.

2.4 The SABA integrator family

A different approach to extend the WH integrator to higher order has been used by Laskar & Robutel (2001). To construct their SABA family of symplectic integrators, these authors use a Lie series approach to build composition methods similar to Yoshida (1993). An additional constraint that Laskar & Robutel (2001) enforce is that each sub-step in the composition method has only positive time-steps. They argue that this avoids large coefficients in the error terms, which might be beneficial for large time-steps

Specifically, their integrators SABA1, SABA2, SABA3, and SABA4, have error terms |$\mathcal {O}(\epsilon \mathrm{ d}t^2)$|⁠, |$\mathcal {O}(\epsilon \mathrm{ d}t^4+\epsilon ^2 \mathrm{ d}t^2)$|⁠, |$\mathcal {O}(\epsilon \mathrm{ d}t^6+\epsilon ^2 \mathrm{ d}t^2)$|⁠, and |$\mathcal {O}(\epsilon \mathrm{ d}t^6+\epsilon ^2 \mathrm{ d}t^2)$|⁠, respectively. They are compositions of the same e^A and e^B operators used in the classical WH method. For example, one SABA2 step can be written as

$$\begin{eqnarray*} \text{SABA2} &=& e^{c_1 A} e^{d_1 B} e^{d_2 A} e^{d_1 B} e^{c_1 A}, \end{eqnarray*}$$

(20)

where the coefficients c₁, d₁, and d₂ are positive constants. It is referred to as SABA2 because it involves two force evaluations per time-step. Similarly SABA3 has three force evaluations per time-step, and so forth. SABA1 is simply the classical WH method that uses one force evaluation.

As Laskar & Robutel (2001) point out and we will show later, the dominant term for small time-steps of the SABA methods is typically the |$\mathcal {O}(\epsilon ^2\mathrm{ d}t^2)$| term except for SABA1 (which is just WH). To remove this error term, Laskar & Robutel (2001) describe a corrector step for their integrators, leaving only terms of order |$\mathcal {O}(\epsilon \mathrm{ d}t^p + \epsilon ^2 \mathrm{ d}t^4)$| if used in conjunction with the SABAp integrators. We refer to these integrators as SABACp.

Is it important to stress that the SABACp correctors are very different from the ones described by Wisdom et al. (1996), both in terms of concept and implementation. For example, consider SABAC2 which we can write as

$$\begin{eqnarray*} \text{SABAC2} &=& e^{-\alpha C_S} e^{c_1 A} e^{d_1 B} e^{d_2 A} e^{d_1 B} e^{c_1 A} e^{-\alpha C_S}, \end{eqnarray*}$$

(21)

where α is a positive constant. The SABACp correctors |$e^{C_S}$| are simply given by the evolution of the system under the dominant term in the shadow Hamiltonian, {B, {B, A}} and the coefficient α is directly related to the coefficient of this term as it appears in the shadow Hamiltonian. Thus, the |$e^{\alpha C_s}$| operators explicitly remove one of the error terms accumulated over one time-step. Note that this is the same Poisson bracket we encountered in the kernel method WHCKM above. However, the interpretation is different. Further note that in contrast to the WHCp correctors that are only applied at the beginning and end of the integration, the SABACp correctors do not cancel out and have to be applied at every time-step.

To calculate these correctors, we need an operator corresponding to the evolution under the Hamiltonian {B, {B, A}}. In the case of the standard N-body problem with Newtonian gravity, the operator can be calculated analytically. In fact, it is the same calculation as for the modified kick in the WHCKM integrator except that we throw away the standard part of the kick and only keep the modifications. Thus, this corrector implementation takes about as much time as a normal force evaluation. This is the approach taken by Laskar & Robutel (2001).

However, since we have encountered different ways to calculate the kernel in the last section, it might come at no surprise that there is another option to calculate the |$e^{C_S}$| correctors. To our knowledge this method has not been described elsewhere. It uses the same idea as the lazy implementer’s kernel in WHCKL. However, rather than advancing the momenta by the modified kick, |$p_j\rightarrow p_j+\mathrm{ d}t\, m_j\, a^{\prime }_j$|⁠, we only advance them by the modification, i.e. |$p_j\rightarrow p_j+\beta \, \mathrm{ d}t\, m_j\, (a^{\prime }_j-a_j)$| with some constant β(α). As for WHCKL, this method has the disadvantage that it only approximates the evolution under {B, {B, A}} and will therefore introduce errors at higher order. But it also has the same advantage in that it works with any arbitrary position dependent force without the need to derive an analytic expression for {B, {B, A}}. We refer to this integrator as SABACL (SABA integrator with lazy correctors).

Laskar & Robutel (2001) also describe higher order SABA methods, the analogue SBAB integrators where the e^A and e^B operators are exchanged, as well as even higher order methods that can be constructed by compositions of lower order SABA and SBAB integrators. Since these integrators seem to be less efficient and useful for typical N-body simulations of planetary systems (the long and highly accurate integrations of Laskar et al. 2011, use the SABAC4 method), we do not consider them in this paper.

2.5 SABA integrators with negative time-steps

The dominant error term in the SABA4 integrator with correctors is of order |$\mathcal {O}\left(\epsilon ^2 \mathrm{ d}t^4\right)$| for small time-steps. Laskar & Robutel (2001) enforced that all sub-steps have positive time-steps. To achieve a higher order in dt, one needs to relax this condition. Doing so, Blanes et al. (2013) present methods with leading error terms of |$\mathcal {O}\left(\epsilon \mathrm{ d}t^{10} + \epsilon ^2 \mathrm{ d}t^4\right)$|⁠, |$\mathcal {O}\left(\epsilon \mathrm{ d}t^8 + \epsilon ^2 \mathrm{ d}t^6 + \epsilon ^3 \mathrm{ d}t^4\right)$|⁠, and |$\mathcal {O}\left(\epsilon \mathrm{ d}t^{10} + \epsilon ^2 \mathrm{ d}t^6 + \epsilon ^3 \mathrm{ d}t^4\right)$|⁠. The specific methods chosen by Blanes et al. (2013) have particularly small constants associated with their error terms. We refer to these methods as SABA(10,4), SABA(8,6,4), and SABA(10,6,4). They need 7, 7, and 8 force evaluations per time-step, respectively, but do not require a corrector step. Therefore these methods work with any position dependent forces. In contrast to the SABA methods in the last section, The negative time-steps make the method here harder to use with dissipative forces. Blanes et al. (2013) also discuss methods using heliocentric coordinates. Since these methods have almost identical efficiency compared to those using Jacobi coordinates, we do not further consider them here.

2.6 Comparison table for integrators

Whereas the WH method is well known by astronomers and the go-to method for long-term integrations of planetary systems, the higher order methods are used less often. The main motivation behind writing this paper is to clear up some of the mystification regarding the various higher order integrators and make them readily available for astronomers to use.

To do that we summarize all integrators that we have introduced above in Table 1. The first column lists the name of the methods and synonyms various authors have used for them. The second column lists the main reference(s) for each integrator. The third column lists the operators that need to be applied to start up and shut down an integrator. |$C^{(1)}_{p}$| indicates first symplectic correctors of order p and C⁽²⁾ indicates second symplectic correctors. The fourth column lists the operators that need to be applied at each time-step. Note that in some cases it is possible to combine the last operator of the previous step with the first operator of the current step. The symbol A corresponds to e^A, the evolution under Hamiltonian A that describes the Keplerian motion of planets. Similarly, B corresponds to e^B, the planet–planet interactions or the kick step. For integrators which use a variant of the modified kick operator, we denote |$\widehat{B}$|⁠. Similarly, for integrators which use a variant of the lazy implementer’s method, we denote |$\widehat{BB}$| (two B’s because it involves two force evaluations).

Table 1.

Comparison of higher order symplectic splitting integrator used in planetary dynamics. Integrators with particularly useful properties have been marked with ★. See the text for details.

Name and synonyms	Main references	Start up/shut down	One time-step	Cost	Only A, B	|$\mathcal {O}\left(\epsilon \, \mathrm{ d}t^?\right)$|	|$\mathcal {O}\left(\epsilon ^2\, \mathrm{ d}t^?\right)$|	|$\mathcal {O}\left(\epsilon ^3\, \mathrm{ d}t^?\right)$|	Implemented in rebound
WH ★ SABA1 (d) WHFAST (e) M2 (b)	(a), (e), (f)	–	|$A\, B\, A$|	1	|$\checkmark$|	2			|$\checkmark$|
WHCp CM2 (b)	(b), (c), (f)	|$C^{(1)}_{[p]}$|	|$A\, B\, A$|	1	|$\checkmark$|	p + 1	2		|$\checkmark$|(up to p = 17)
WHCCKI (ideal kernel)	(b)	|${C}^{(2)}_*\, {C}^{(1)}_*$|	|$A\, B_*\, A$| * = not possible			∞	∞	4	(not possible)
WHCKC (comp. kernel)	(b)	|$C^{(1)}_{[17]}$|	|$A\, (B\, A)^5$|	5	|$\checkmark$|	18	4		|$\checkmark$|
WHCKM (mod. kick kernel) CMM4 (b)	(b)	|$C^{(1)}_{[17]}$|	|$A\, \widehat{B}\, A$|	1		18	4		|$\checkmark$|
WHCKL ★ (lazy impl. kernel)	(b)	|$C^{(1)}_{[17]}$|	|$A\, \widehat{B\, B}\, A$|	2	|$\checkmark$|	18	4	3	|$\checkmark$|
WHCCKC (comp. kernel)	(b)	|$C^{(2)}\, C^{(1)}_{[17]}$|	|$A\, (B\, A)^5$|	5	|$\checkmark$|	18	4		|$\checkmark$|
WHCCKM (mod. kick kernel) WHCK (f)	(b) (f)	|$C^{(2)}\, C^{(1)}_{[17]}$|	|$A\, \widehat{B}\, A$|	1		18	4		|$\checkmark$|
WHCCKL (lazy impl. kernel)	(b)	|$C^{(2)}\, C^{(1)}_{[17]}$|	|$A\, \widehat{B\, B}\, A$|	2	|$\checkmark$|	18	4	3	|$\checkmark$|
SABA2	(d)	–	|$A\, B\, A\, B\, A$|	2	|$\checkmark$|	4	2		|$\checkmark$|
SABA3	(d)	–	|$A\, (B\, A)^3$|	3	|$\checkmark$|	6	2		|$\checkmark$|
SABA4	(d)	–	|$A\, (B\, A)^4$|	4	|$\checkmark$|	6	2		|$\checkmark$|
SABAC2	(d)	–	|$\widehat{B} A\, (B\, A)^2 \widehat{B}$|	3		4	4		|$\checkmark$|
SABAC3	(d)	–	|$\widehat{B} A\, (B\, A)^3 \widehat{B}$|	4		6	4		|$\checkmark$|
SABAC4	(d), (g)	–	|$\widehat{B} A\, (B\, A)^4 \widehat{B}$|	5		6	4		|$\checkmark$|
SABACL4 ★ (lazy impl. cor.)	This paper	–	|$\widehat{B\, B} A\, (B\, A)^4 \widehat{B\, B}$|	6		6	4	3	|$\checkmark$|
SABA(10,4)	(h)	–	|$A\, (B\, A)^7$|	7	|$\checkmark$|	10	4		|$\checkmark$|
SABA(8,6,4)	(h)	–	|$A\, (B\, A)^7$|	7	|$\checkmark$|	8	6	4	|$\checkmark$|
SABA(10,6,4) ★	(h)	–	|$A\, (B\, A)^8$|	8	|$\checkmark$|	10	6	4	|$\checkmark$|

Notes. (a) Wisdom & Holman (1991), (b) Wisdom et al. (1996), (c) Wisdom (2006), (d) Laskar & Robutel (2001), (e) Rein & Tamayo (2015), (f) Wisdom (2018), (g) Laskar et al. (2011), (h) Blanes et al. (2013).

We list the theoretical cost of each method in column five. It corresponds to the number of force evaluations per time-step and assumes all other operations take no time. If the last operator of the step involves a force evaluation and it can be combined with the first operator of the next step, then we only count it as one force evaluation. This column thus provides a runtime estimate relative to the classical WH method for a fixed time-step. We come back to this later, but note that two methods already stands out: WHCKM and WHCCKM. We should expect these integrators to be as fast as the classical WH method.

The sixth column lists if the method requires only the operators e^A and e^B. If this is the case, then the implementation is straightforward and just comes down to repeatedly applying e^A and e^B in the right order and for the right amount of time. Note that if other operators need to be implemented, then it is harder to use the method for simulations that include effects other than Newtonian gravity between all particles because these special operators need to be rederived analytically.⁸

The seventh, eighth, and ninth columns list the order of the leading error terms. A number k in the seventh column implies that there is an error term of order |$\mathcal {O}\left(\epsilon \mathrm{ d}t^k\right)$|⁠, a number k in the eighth column implies that there is an error term of order |$\mathcal {O}\left(\epsilon ^2 \mathrm{ d}t^k\right)$|⁠, and similarly for the ninth column. A column is left blank if the error is always dominated by the error given by the column to the left, regardless of the relative sizes of dt and ϵ.

The last column lists if the method is implemented in the rebound integrator package. We will describe the details of our implementations in the next section.

3 IMPLEMENTATION

3.1 WHFast extensions

We implement all methods described in Section 2.3 as extensions to the WHFast integrator in rebound . To do this we add two new parameters to the ri_whfast structure. The first is kernel that controls what kind of kernel is used. The default setting (DEFAULT) uses the standard WH method’s kernel. The other settings available are COMPOSITION, MODIFIEDKICK, and LAZY, corresponding to the integrators WHCKC, WHCKM, and WHCKL, respectively. The second new parameter is corrector2 that either turns the second correctors on or off. We also extend the first correctors already implemented in WHFast and controlled by the corrector parameter to 17th order.

The new kernel methods are currently implemented in a way that makes them somewhat more restricted in their usage than the basic WHFast algorithm. Specifically, the newly implemented kernels do not support variational equations (and therefore chaos indicators such as MEGNO), OpenMP parallelization, or force calculations using a BH tree code. For most cases where higher order symplectic methods might be used, these features are not essential. However, some of these features can be added at a later time if there is a need.

All kernel methods make use of the WHFast setting safe_mode which, when turned off, combines the last e^A operator in each time-step with the e^A operator at the beginning of the next time-step. If the safe mode is turned off, then Jacobi coordinates are only converted back to Cartesian coordinates at the end of the integration. In most cases it makes sense to turn-off the safe mode flag as it not only provides a speed-up, but also reduces round-off errors coming from continuous transformations to and from Jacobi coordinates.

If frequent outputs are required WHFast needs to apply the correctors and their inverses repeatedly. Round-off error can prevent the correctors and their inverses from cancelling out exactly. If this becomes a problem, one can turn on the keep_unsynchronized setting of WHFast which then only applies the correctors to generate an output in Cartesian coordinates, but continues the integration from a copy of the Jacobi coordinates it made before the correctors were applied.

rebound  allows the user to specify a routine to include additional forces that can be used to model effects due to general relativistic corrections, oblateness, or tides. As long as these forces are position dependent, all of the new kernel methods support additional forces, with the exception of MODIFIEDKICK.

Furthermore, all new kernels are compatible with the SimulationArchive (Rein & Tamayo 2017). Specifically, all kernel methods are bit-wise reproducible from any snapshot and are machine independent.

3.2 SABA

We implement the SABA integrator family described in Sections 2.4 and 2.5 as a new integrator saba in rebound . The parameter type in the ri_saba structure determines which specific SABA integrator is used. Currently, the three high-order integrators of Blanes et al. (2013), SABA(10,4), SABA(8,6,4), and SABA(10,6,4) are implemented. In addition, the integrators SABA1, SABA2, SABA3, and SABA4 with either no correctors, lazy correctors, or modified kick correctors are implemented. It is straightforward to extend the implementation to other variants should there be a need.

Both the keep_unsynchronized and the safe_mode parameters in the ri_saba structure works the same way as for WHFast . If the safe mode is turned off, it combines the e^A operators at the beginning and end of consecutive time-steps if no corrector is used. If a corrector is used, then the correctors at the beginning and end of a time-step are combined. As with WHFast , in most cases it makes sense to turn-off the safe mode flag to provide a speed-up and reduce round-off errors. Note that it is called safe mode because, when turned off, changing particle properties manually in-between time-steps requires extra scrutiny.

The modified kick corrector is not compatible with additional forces, however, SABACL integrators using the lazy correctors are, as are all the SABA integrator that do not use a corrector.

The SABA integrators as implemented in rebound rely on many of the internal WHFast methods and are therefore also compatible with the SimulationArchive, bit-wise reproducible, and machine independent.

4 TESTS

We run simulations of the outer Solar system for 10 kyr as a test case of our rebound implementations of the integrators introduced above. A python notebook to reproduce all of the test and figures can be found at https://github.com/hannorein/ReinTamayoBrown2019.

Fig. 1 shows the maximum relative energy error as a function of the time-step. We measure the energy error at 10 000 random times during each integration to avoid any aliasing. The left-hand panel shows the methods of Wisdom et al. (1996) whereas the right-hand panel shows those of Laskar & Robutel (2001) and Blanes et al. (2013). All methods are dominated by stepsize resonances (Rauch & Holman 1999) for time-steps larger than about 10 per cent of the shortest orbital time-scale in the problem. Note that for such large time-steps, the SABA type integrators perform better than others. The first reason for this is that the SABAp methods with p > 1 have more than one force evaluation during the time-step. Thus, the effective time-step in-between force evaluations is smaller than the time-step plotted on the horizontal axis. Of course, more force evaluations comes at the cost of being slower (see below). The second reason for the good behaviour of the SABA integrators at large time-steps is that the coefficients in the error terms are generally small even at higher order that is beneficial for large time-steps (Laskar & Robutel 2001).

Figure 1.

The maximum relative energy error as a function of the time-step in 10 kyr integrations of the outer Solar system using different symplectic integrators.

For more reasonable time-steps of a few per cent of the shortest orbital period, the methods show a power-law convergence as indicated by the terms in Table 1. In particular, note that the WH and WHC17 integrators follow a |$\mathcal {O}\left(\epsilon \mathrm{ d}t^2\right)$| and |$\mathcal {O}\left(\epsilon ^2 \mathrm{ d}t^2\right)$| power law, respectively. In other words, the error of WHC17 is reduced by one extra factor of ϵ compared to the standard WH integrator. From the plot, we can read off a value ϵ ∼ 10⁻³ for this test case, roughly equal to the Jupiter–Sun mass ratio. For small enough time-steps, once the |$\mathcal {O}\left(\epsilon \right)$| term is negligible, the SABA2, SABA3, and SABA4 integrators⁹ follow the same |$\mathcal {O}\left(\epsilon ^2 \mathrm{ d}t^2\right)$| power law as WHC17. In other words, going beyond SABA2 does not help to improve the accuracy significantly for small time-steps. Note however that SABA3 is slightly more accurate than SABA2, and SABA4 is slightly more accurate than SABA3. Once again, this can be explained because SABA3 evaluates the forces three times during a time-step, whereas SABA2 evaluates them only twice. Thus, the effective time-step of SABA3 is actually smaller than that of SABA2. The correctors of the SABAC integrators successfully remove the |$\mathcal {O}\left(\epsilon ^2 \mathrm{ d}t^2\right)$| terms, with higher order SABAC integrators having an advantage at large time-steps since the leading order errors fall off faster for higher order SABA integrators. The higher order integrators of Blanes et al. (2013) have even smaller energy errors.

Going back to the left-hand panel of Fig. 1, we can see that all kernel methods, WHCKC, WHCKM, and WHCKL, perform equally well. Thus, we conclude that any of the three approximations for the kernel work equally well for this problem, which is consistent with the results of Wisdom et al. (1996).

For very small time-steps, the integrators are dominated by numerical round-off error coming from the finite precision of floating point numbers. If all round-off errors are unbiased, then this error term behaves like a random walk and scales as |$\mathcal {O}(\mathrm{ d}t^{-1/2})$|⁠. This is known as Brouwer’s law (Brouwer 1937). This behaviour can be observed for all integrators in Fig. 1 that reach machine precision and is expected as all integrators use internally the unbiased WHFastKepler solver (Rein & Tamayo 2015). Worth noting is that the WHCKC method has a slightly larger round-off error than the WHCKM and WHCKL methods because it requires more e^A and e^B operators for one time-step.

To better compare the runtime performance of the integrators, we plot the maximum relative energy error as a function of the runtime in Fig. 2. The simulations were performed on an Intel Xeon CPU (E5-2620 v3, 2.40 GHz). The horizontal axis has been scaled so that it shows the number of hours required for a 1 billion year integration for this particular problem and CPU. Amongst all integrators, WHCKL is the fastest for moderate accuracy, ΔE/E ≳ 10⁻¹² (we plot it on both panels for comparison). It only requires 1 h to integrate the outer Solar system for 1 Gyr at a maximum relative energy error of 10⁻¹⁰. Since the dominant error term, |$\mathcal {O}\left(\epsilon ^2 \mathrm{ d}t^4\right)$|⁠, falls off as the fourth power of the time-step, an accuracy of 10⁻¹⁴ can be achieved in only 10 h. The WHCKM integrator with the modified kick step performs almost as well. For simulations that require extremely high accuracy, the higher order SABA integrators perform best. In particular, the SABA(10,6,4) integrator is more efficient than the WHCKL method when ΔE/E ≲ 10⁻¹² is required. It requires a runtime of 4 h to integrate the outer Solar system for 1 Gyr at a maximum relative energy error of 10⁻¹⁴, roughly a factor of 2 faster than WHCKL.

Figure 2.

The maximum relative energy error in simulations of the outer Solar system as a function of the runtime required to reach 1 Gyr using different symplectic integrators.

The fastest integrators of the SABA family for high accuracy simulations is SABA(10,6,4). We can see that for the same accuracy, the fastest SABACL/SABAC integrator is roughly a factor of 2–5 slower than the fastest integrators in our sample, WHCKL and SABA(10,6,4). This is consistent with the results of Wisdom (2018).

Not shown in the plots are the integrators of Wisdom et al. (1996) that use second symplectic correctors, WHCCKM, WHCCKL, and WHCCKC. We do not observe any improvement of the energy error in this test case compared to their counterparts with only first symplectic correctors applied. We also do not show the SABAC integrators because they perform similarly to their SABACL counterparts.

5 CONCLUSIONS

In this paper, we reviewed different high-order symplectic integrators for long-term direct N-body simulations of planetary systems. We implemented all of these integrators and make them freely available as part of the rebound integrator package. Some of these method have a truly remarkable performance with little to no additional cost associated when compared to the (already impressive) standard second-order WH method. For example, a typical integration of the Solar system that takes the WHCKL method 10 h to complete, would require more than a year if one were to use the standard WH method and require the same level of accuracy.

The best-performing integrators in our sample, WHCKL, SABA4CL, and SABA(10,6,4) use very different approaches to achieve a high accuracy. Reassuringly, the different approaches lead to integrators with similar performance. We find that the WHCKL integrator has a small¹⁰ advantage over the best integrators in the SABA family for moderate accuracies, being about 2–3 times faster. On the other hand, the SABA(10,6,4) integrator is faster for very high accuracy runs, ΔE/E ≲ 10⁻¹².

Aside of speed and accuracy, the integrators differentiate in other ways as well. The SABA, SABACL, WHCKC, and, WHCKL integrators only require operators that are already present in the standard WH integrator, i.e. a Kepler solver and a routine calculating the interaction terms. This makes their implementation very straightforward. Furthermore, they can be used together with forces other than Newtonian gravity, as long as these forces only depend on the particles’ positions. Astrophysically relevant forces with such a property include general relativistic corrections, quadrupole and other higher order moments of non-spherical objects, and some descriptions of tidal and radiation effects. The other integrators, SABAC and WHCKM, can also be used for these cases, but some work is needed in addition to the force implementation itself.

Our implementations do not support extended precision. For very small time-steps, finite precision of floating point numbers is therefore the limit factor in achieving even higher precision. There are currently few problems where this level of precision is required. One exception are Solar system integrations where our understanding of the physical system is now comparable to this level of precision (Laskar et al. 2011). Either a calculation in full quadruple precision or some form of compensated summation (Wisdom 2018) can be used to go beyond the limits of double precision floating point numbers.

In summary, for integrations of planetary systems in which orbits remain well separated, we recommend WHCKL, the WH method with the lazy implementer’s kernel and first symplectic correctors of order 17. To use WHCKL in rebound , one can simply set sim.integrator = "WHCKL", which configures the WHFast parameter such that it corresponds to WHCKL.¹¹ The speed of this method is very similar to the standard WH method, but the accuracy is superior in almost all cases and it can be used with a wide variety of additional forces (Tamayo et al. 2019). For very high accuracy runs, we recommend the SABA(10,6,4) integrator. To use SABA(10,6,4) in rebound , simply set sim.integrator = "SABA(10,6,4)". For systems in which close encounters occur, a different approach is needed (Rein & Spiegel 2015; Rein et al. 2019).

ACKNOWLEDGEMENTS

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).

Footnotes

1

Whereas in the Solar system the shortest orbital period is 88 d, there are many extrasolar planets with extremely short periods. Integrating a planet on a one day orbit for 10 Gyr requires roughly 10¹⁴ time-steps.

2

There is a small subtlety that originates from a permutation of the exponentials known as the Vertauschungssatz (Gröbner 1967). This has occasionally led to a sign error in the literature. The fact that there is no sign ambiguity at even orders did not help to clear up the confusion. See Lemma 5.1 in chapter III of Hairer et al. (2006).

3

As Wisdom (2018) explains, this argument is flawed because it ignores the fact that the series above might not converge everywhere in phase space. The system we end up actually solving might look nothing like the system we want to solve. The interpretation preferred by Wisdom (2018) avoids some of these problems. However, for the discussion in this paper, this is not a concern.

4

The leap-frog integrator is also a second-order symplectic splitting method, but splits the Hamiltonian into a potential and kinetic term. Both terms are equally important to describe the motion of a planet, thus ϵ ∼ 1.

5

Given that the difference between real and integrator coordinates also depends on the time-step (and vanishes in the limit dt → 0), Saha & Tremaine (1992) show that a ‘warm-up’ procedure that slowly increases the time-step from 0 to the desired finite value dt achieves essentially the same result as an explicit symplectic corrector to integrator variables. We do not discuss this method further because it has almost identical properties to the ones we describe below using symplectic correctors.

6

The jerk is the time derivative of the acceleration.

7

General relativistic corrections are in fact velocity dependent. However, the perihelion precession can be modelled with a velocity independent r⁻³ potential.

8

An autodifferentiation algorithm could do this too.

9

Except SABA1 that is equivalent to the standard WH method.

10

‘Small’ depends on the context. If one has to wait for a year for a simulation to finish, a factor 3 increase in performance might be a huge deal.

11

The settings safe_mode and keep_unsynchronized are not changed by this shortcut.

Author notes

NHFP Sagan Fellow.