REBOUNDx: A Library for Adding Conservative and Dissipative Forces to Otherwise Symplectic N-body Integrations

TL;DR

REBOUNDx extends symplectic N-body integrators to include dissipative and conservative forces, clarifying their integration properties and enabling accurate, reproducible astrophysical simulations with an open-source library.

Contribution

It introduces a general framework for incorporating dissipative forces into symplectic schemes and demonstrates the effectiveness of symplectic correctors for weakly dissipative systems.

Findings

Symplectic schemes can be extended to dissipative forces using non-commutative operators.

Symplectic correctors are effective for weakly dissipative systems, acting as 'weak splitting correctors.'

Incorporating velocity-dependent conservative forces requires careful treatment to avoid qualitative errors.

Abstract

Symplectic methods, in particular the Wisdom-Holman map, have revolutionized our ability to model the long-term, conservative dynamics of planetary systems. However, many astrophysically important effects are dissipative. The consequences of incorporating such forces into otherwise symplectic schemes is not always clear. We show that moving to a general framework of non-commutative operators (dissipative or not) clarifies many of these questions, and that several important properties of symplectic schemes carry over to the general case. In particular, we show that explicit splitting schemes generically exploit symmetries in the applied external forces which often strongly suppress integration errors. Furthermore, we demonstrate that so-called `symplectic correctors' (which reduce energy errors by orders of magnitude at fixed computational cost) apply equally well to weakly dissipative systems and can thus be more generally thought of as `weak splitting correctors.' Finally, we show that previously advocated approaches of incorporating additional forces into symplectic methods work well for dissipative forces, but give qualitatively wrong answers for conservative but velocity-dependent forces like post-Newtonian corrections. We release REBOUNDx, an open-source C library for incorporating additional effects into REBOUND N-body integrations, together with a convenient Python wrapper. All effects are machine-independent and we provide a binary format that interfaces with the SimulationArchive class in REBOUND to enable the sharing and reproducibility of results. Users can add effects from a list of pre-implemented astrophysical forces, or contribute new ones.

Equations

$$\hat{P}\mathbf{z}$$

$$\hat{A}\mathbf{z}$$

$$\hat{B}\mathbf{z}$$

$$\mathbf{z}(t+h)\approx \mathcal{S}\mathcal{A}\mathcal{B}(h)\mathbf{z}(t)\equiv \mathcal{A}(h)\circ \mathcal{B}(h)\mathbf{z}(t),$$

$$\dot{\mathbf{z}}={\mathcal{L}}_{\mathcal{S}\mathcal{A}\mathcal{B}}\mathbf{z}={\mathcal{L}}_P\mathbf{z}+{\mathcal{L}}_{\mathcal{S}\mathcal{A}\mathcal{B}}^{\text{err}}\mathbf{z},$$

$${\mathcal{L}}_{\mathcal{S}\mathcal{A}\mathcal{B}}^{\text{err}}=-\frac{h}{2}[{\mathcal{L}}_A,{\mathcal{L}}_B]+\frac{h^2}{12}[{\mathcal{L}}_A-{\mathcal{L}}_B,[{\mathcal{L}}_A,{\mathcal{L}}_B]]-\mathcal{O}\left(h^3\right).$$

$$\mathcal{S}\mathcal{A}\mathcal{B}\mathcal{A}(h)\equiv \mathcal{A}(h/2)\circ \mathcal{B}(h)\circ \mathcal{A}(h/2).$$

$${\mathcal{L}}_{\mathcal{S}\mathcal{A}\mathcal{B}\mathcal{A}}^{\text{err}}=\frac{h^2}{24}[2{\mathcal{L}}_B+{\mathcal{L}}_A,[{\mathcal{L}}_B,{\mathcal{L}}_A]]+\mathcal{O}\left(h^4\right).$$

$$H_{\mathcal{S}\mathcal{A}\mathcal{B}}^{\text{err}}=\frac{h}{2}\{H_A,H_B\}+\frac{h^2}{12}\{H_A-H_B,\{H_A,H_B\}\}+\mathcal{O}(h^3),$$

$$H_{\mathcal{S}\mathcal{A}\mathcal{B}\mathcal{A}}^{\text{err}}=\frac{h^2}{24}\{2H_B+H_A,\{H_B,H_A\}\}+\mathcal{O}(h^4),$$

$${\mathcal{L}}_K\mathbf{z}$$

$${ϵ}_p{\mathcal{L}}_I\mathbf{z}$$

$${\mathbf{v}}_{\mathbf{i}}(t+h)={\mathbf{v}}_{\mathbf{i}}(t)+{ϵ}_{p}h{\mathbf{a}}_{\mathbf{i}}^{\mathbf{int}}.$$

$$\Bigg{|}\frac{\Delta E^{\rm err}}{E}\Bigg{|}\sim\Bigg{|}\frac{T\dot{E}^{\rm err}}{E}\Bigg{|}=\Bigg{|}\frac{T{\mathcal{L}}_{S}^{\rm err}E}{E}\Bigg{|},$$

$$V_{J2}(r)=-\frac{1}{2}J_{2}\Bigg{(}\frac{R}{r}\Bigg{)}^{2}\frac{GM}{r}\equiv\epsilon_{J2}\frac{GM}{a_{0}}\Bigg{(}\frac{a_{0}}{r}\Bigg{)}^{3}.$$

$$\mathbf{F}=-m\frac{\mathbf{v}}{2{\tau }_a},$$

$$\mathcal{K}(h/2)\circ \mathcal{D}(h)\circ \mathcal{K}(h/2),$$

$${\mathcal{C}}^{-1}(h)\circ \mathcal{S}\mathcal{A}\mathcal{B}\mathcal{A}(h)\circ ...\circ \mathcal{S}\mathcal{A}\mathcal{B}\mathcal{A}(h)\circ \mathcal{C}(h),$$

$$\mathcal{P}(h)$$

$$\mathcal{S}\mathcal{A}\mathcal{B}(h)\approx \mathcal{A}(h/2)\circ \mathcal{P}(h)\circ \mathcal{A}(-h/2),$$

$$\mathcal{K}(h)\circ \mathcal{G}\mathcal{R}(h),$$

$$\hat{K}\mathbf{z}$$

$$\widehat{GR}\mathbf{z}$$

$$\Delta z_{LT}^{\mathrm{err}}=\frac{f^{n+1}(\xi )}{(n+1)!}{h}^{n+1},$$

$$\frac{\Delta z_{LT}^{\rm{err}}}{z}\sim\frac{(2\pi)^{n+1}\epsilon^{n+1}}{(n+1)!}\Bigg{(}\frac{h}{T}\Bigg{)}^{n+1}.$$

$$\frac{\Delta z_{LT}^{\rm{err}}}{z}\sim\frac{(2\pi)^{n+1}\epsilon^{n+1}}{(n+1)!}\Bigg{(}\frac{h}{T}\Bigg{)}^{n}\Bigg{(}\frac{t}{T}\Bigg{)}$$

$$\mathcal{K}(h/2)\circ \mathcal{D}(h/2)\circ \mathcal{I}(h)\circ \mathcal{D}(h/2)\circ \mathcal{K}(h/2),$$

$$\dot{\mathbf{z}}={\mathcal{L}}_P\mathbf{z},$$

$${\mathcal{L}}_{F}g(\mathbf{z})=\sum _i\frac{\partial g}{\partial z_i}{\dot{z}}_i,$$

$${\bf z}\Bigg{|}_{t+h}={\bf z}\Bigg{|}_{t}+h\frac{d{\bf z}}{dt}\Bigg{|}_{t}+\frac{h^{2}}{2!}\frac{d^{2}{\bf z}}{dt^{2}}\Bigg{|}_{t}+...\equiv\Bigg{[}(1+h{\mathcal{L}}_{P}+\frac{h^{2}}{2!}{\mathcal{L}}_{P}^{2}+...){\bf z}\Bigg{]}_{t}.$$

$$\mathbf{z}(t+h)\equiv e^{h{\mathcal{L}}_P}\mathbf{z}(t)\equiv \mathcal{P}(h)\mathbf{z}(t).$$

$$e^{h{\mathcal{L}}_A}{e}^{h{\mathcal{L}}_B}=e^{h{\mathcal{L}}_S}$$

Full text

1 Introduction

The long-term dynamical evolution of planetary systems remains a rich challenge for analytical investigations. As a result, many of the advances in the last several decades have been driven by numerical studies thanks to faster computers and the development of improved algorithms for accurate and efficient numerical solutions, e.g., the chaotic evolution of Pluto (Sussman & Wisdom, 1988), chaos in the inner solar system (Laskar, 1989), and the marginal instability of Mercury over the age of the solar system Laskar & Gastineau (2009). In many cases, such chaos only manifests itself after millions or even billions of orbits. This imposes strong demands on the long-term conservation properties of numerical integration methods that can be used for such studies.

These concerns led to the development powerful symplectic integration techniques (e.g., Ruth, 1983; Neri, 1987; Forest & Ruth, 1990; Yoshida, 1990). Most notably, by exploiting the near-Keplerian planetary motions, the Wisdom-Holman map (Wisdom, 1982; Kinoshita et al., 1990; Wisdom & Holman, 1991) enabled the first direct N-body integrations of our solar system over Gyr timescales. By respecting the Hamiltonian structure of the problem, symplectic methods render the task of integration equivalent to a canonical transformation (e.g. Sanz-Serna, 1992), which enforces the conservation of several invariants of the phase space flow. This avoids the secular energy error growth of many of the more general methods that are not restricted to Hamiltonian systems, which can lead to unphysical collisions in long-term simulations.

Such symplectic integrators split the problem into multiple pieces that are each evolved in sequence. Apart from the robust geometrical and numerical properties this provides, the modularity of splitting methods makes it possible to extend them for hierarchical problems (Fujii et al., 2007; Portegies Zwart & McMillan, 2018) and to incorporate different physics in a single integration (Pelupessy & Portegies Zwart, 2012; Portegies Zwart, 2018). This has for example been vigorously pursued in the Astrophysics Multipurpose Software Environment (AMUSE) package (Portegies Zwart & McMillan, 2018).

Some effects, such as bodies’ higher gravitational moments, can be expressed as simple position-dependent potentials and can be trivially incorporated into symplectic schemes. However, many crucial perturbations such as post-Newtonian corrections, tides, or radiation forces are velocity-dependent and even dissipative. Malhotra (1994) and Cordeiro et al. (1996) have proposed generalizations of symplectic schemes for incorporating weak velocity-dependent forces, but their consequences in long-term integrations are not well understood. While some authors perform convergence tests to check the validity of their adopted stepsize using variants of the above methods (e.g. Zhang & Hamilton, 2007), the numerical robustness of studies where this is not done is unclear.

A generalization of symplectic integration theory would therefore be valuable to understand the source of numerical error in the above schemes, and to provide simple estimates for its magnitude to guide the selection of timestep for a given problem.

A promising direction has been offered by Galley (2013) and Galley et al. (2014), who have developed a generalized theory for classical mechanics that can incorporate dissipative processes in the action governing the dynamics. One can use this formalism to construct variational integrators that generalize symplectic methods to accurately track the changing momentum and energy of the system on long timescales (Tsang et al., 2015). Potential disadvantages of this elegant framework are that the equations of motion are necessarily implicit (Tsang et al., 2015), and that adding new forces becomes somewhat more complicated than with traditional methods, since one needs to derive a ‘nonconservative potential’ using the formalism in Galley et al. (2014).

In this paper, we explore how the established theory of symplectic integration can instead be generalized in a framework of non-commutative operators, which can be applied to both conservative and dissipative systems. Such generalized ‘splitting methods’ have been suggested in the field of fluid mechanics as far back as Strang (1968), and more recently for celestial mechanics by Mikkola (1998)111See also Hernandez & Bertschinger (2018) for error analysis of one-step methods from the perspective of their differential equations.. We show that several strong results traditionally discussed in the symplectic integration literature carry over to more general and possibly dissipative splitting methods.

Readers interested in quickly using REBOUNDx to incorporate astrophysical effects into N-body integrations could begin at Sec. 6 for an overview of the implementation. For a more gradual introduction, we introduce in Sec. 2 the notation, and review symplectic integration and the Wisdom-Holman map in an operator-centred formalism. We then show in Sec. 3 how this framework can explain the energy error behaviour of such methods under weak perturbations, both conservative and dissipative. In Sec. 4 we review how the WH map can be extended to higher order without additional force evaluations, and show that such ‘symplectic correctors’ can also be applied to dissipative systems and are thus more general than typically considered. In Sec. 5, we consider the general case where analytic solutions to individual steps in a split scheme are not known, and how truncation errors from numerically integrating across these individual steps interact with the overall splitting scheme errors. In particular, we show that for conservative, velocity-dependent forces like post-Newtonian corrections, the methods of Cordeiro et al. (1996) and Malhotra (1994) give qualitatively wrong answers, and we show how to correct them. We summarize our results and conclude in Sec. 7.

2 Background

We begin by reviewing the basic theory behind symplectic integration. Several excellent reviews are available (e.g. Kinoshita et al., 1990; Sanz-Serna, 1992), but we choose to present it in an operator-centred framework that will generalize later to dissipative systems. We try to keep the introduction at the level required for the discussion in the main text. For a more careful introduction, and to correct some sign errors in the literature, see Appendix A.

2.1 Splitting Schemes

We consider a system of $N$ particles with positions $\bf r_{i}$ and velocities $\bf v_{i}$. We define a differential operator $\hat{P}$, which acts on the current state of the system ${\bf z}=(\bf r,\bf v)$ to yield $\dot{{\bf z}}$, i.e., the particular set of differential equations we are trying to solve222Throughout this paper we restrict ourselves to time-independent sets of differential equations.,

[TABLE]

where ${\bf a_{i}}$ is the $i$th particle’s acceleration vector, which depends on the positions ${\bf r}$ of all the bodies.

Next, we introduce a further level of abstraction by defining a solution to Eq. 1 through a corresponding operator ${\mathcal{P}}(h)$. This ideal integrator exactly advances the state ${\bf z}$ by a timestep $h$ according to the set of differential equations $\dot{\bf z}=\hat{P}{\bf z}$.

In general, there is no closed form solution $\mathcal{P}$ (e.g., the three-body problem), but one could split the differential equations into two pieces that could each be solved trivially in isolation,

[TABLE]

In particular, the solution ${\mathcal{A}}(h)$ keeps the velocities constant and updates the positions by $h{\bf v_{i}}$, while ${\mathcal{B}}(h)$ keeps the positions constant and similarly updates the velocities using constant accelerations. This splitting works for any differential equations of the form in Eq. 1, and is thus widely applicable. In a Hamiltonian framework, in cases where the accelerations can be derived from a position-dependent potential, it corresponds to splitting the Hamiltonian into a kinetic piece containing all the momenta, and a potential piece containing all the positions. We therefore refer to this scheme as a kinetic-potential splitting, which was the focus of most early symplectic integrators (e.g., Ruth, 1983; Neri, 1987; Forest & Ruth, 1990; Yoshida, 1990). As we will see below, better splittings for specialized problems are possible.

The idea behind splitting methods is to alternate evolution under two or more operators that can be calculated exactly and efficiently. For example, we can define a first-order splitting scheme as

[TABLE]

where ${\mathcal{B}}$ advances the state ${\bf z}(t)$ by $h$ according to its half of the differential equations, and then ${\mathcal{A}}$ advances the result by $h$ with the remaining half. Unfortunately, doing two halves of the problem in sequence is not the same as doing the full problem all at once. Nevertheless, if we repeatedly apply this map to get from some ${\bf z}(t_{1})$ to ${\bf z}(t_{2})$, at least in the limit of a vanishing timestep, the split operations become densely interleaved and this approximation will approach the real solution.

For finite timesteps, we can also see simply that how well a split scheme will correspond to the real solution should depend on the degree to which the operators ${\mathcal{A}}(h)$ and ${\mathcal{B}}(h)$ commute with one another. We could always split each of the steps in Eq. 3 in half and write it as $\mathcal{A}(h/2)\circ\mathcal{A}(h/2)\circ\mathcal{B}(h/2)\circ\mathcal{B}(h/2)$. If the order in which these operators are applied didn’t matter, we could switch the middle two to yield $\mathcal{A}(h/2)\circ\mathcal{B}(h/2)\circ\mathcal{A}(h/2)\circ\mathcal{B}(h/2)$. We could then continue this split-and-switch procedure ad infinitum until the operators were again densely interleaved and we would effectively be applying both at the same time, yielding the true solution. So we expect the errors in splitting schemes to vanish in the limits that either $h$ or the commutator of ${\mathcal{A}}$ and ${\mathcal{B}}$ go to zero.

The formula that quantifies the above two statements is the Baker-Campbell-Hausdorff (BCH) identity, which is useful in many branches of physics involving non-commuting operators like quantum mechanics. To state it, we have to add a final layer of abstraction and generalize beyond just how the phase space variables ${\bf z}$ change under a particular subset of the differential equations (e.g., Eq. 2), to how an arbitrary function $g({\bf z})$ changes with time. In particular, for every subset of differential equations, e.g., ${\bf\dot{z}}=\hat{A}\>{\bf z}$, we can define the corresponding Lie derivative ${\mathcal{L}}_{A}\>$ that, when acting on $g({\bf z})$, yields $dg/dt$ under that particular subset of differential equations. In the simple case of $g({\bf z})={\bf z}$, we get back ${\mathcal{L}}_{A}\>{\bf z}=\dot{{\bf z}}=\hat{A}\>{\bf z}$.

The BCH formula says that when we apply Eq. 3, we are not solving the original set of differential equations ${\mathcal{L}}_{P}\>{\bf z}={\mathcal{L}}_{A}\>{\bf z}+{\mathcal{L}}_{B}\>{\bf z}$ (Eqs. 1 and 2), but rather a nearby set of differential equations

[TABLE]

where the additional error term is a power series in the timestep consisting of nested commutators $[{\mathcal{L}}_{A},{\mathcal{L}}_{B}]={\mathcal{L}}_{A}{\mathcal{L}}_{B}-{\mathcal{L}}_{A}{\mathcal{L}}_{B}$ (e.g. Saha & Tremaine, 1992),

[TABLE]

In other words, if we integrate the set of differential equations in Eq. 4 numerically to high accuracy, we match the evolution generated by the simple splitting scheme in Eq. 3. We will call these differential equations that we are actually solving the ‘modified set of differential equations’ for the splitting scheme. Understanding the structure of this nearby problem gives insight into the numerical errors introduced by the method333It is important to point out that the BCH formula represents a formal, asymptotic series, which in particular is not guaranteed to converge everywhere in phase space (Wisdom, 2018). In practice, as long as the adopted timestep is $\lesssim 10\%$ of the fastest timescale in the problem, this is not typically a concern..

As expected from the qualitative arguments above, the modified differential equations approach the true ones as $h$ approaches zero. We also see that if (and only if) $[{\mathcal{L}}_{A},{\mathcal{L}}_{B}]=0$, then there is no error and ${\mathcal{A}}(h)$ and ${\mathcal{B}}(h)$ commute, i.e., ${\mathcal{A}}(h)\circ{\mathcal{B}}(h)={\mathcal{B}}(h)\circ{\mathcal{A}}(h)$.

By choosing appropriate steps, one can compose ${\mathcal{A}}$ and ${\mathcal{B}}$ into higher order schemes (e.g., Yoshida, 1990). A widely used scheme is the time-symmetric ‘leapfrog’ method

[TABLE]

Repeated application of the BCH formula shows that this particular composition cancels out error terms linear in $h$, yielding a second-order integrator with

[TABLE]

Because the scheme is time-symmetric, errors only appear at even powers of $h$.

2.2 Symplectic Schemes

The above analysis does not depend on us solving a Hamiltonian system, and is general to both conservative and dissipative systems. However, the insight of symplectic integration is that if each of the split operators, ${\mathcal{A}}(h)$ and ${\mathcal{B}}(h)$, are exact solutions to a Hamiltonian system, then each of them conserve various Hamiltonian invariants. This implies that any composition of them, e.g., ${\mathcal{A}}(h)\circ{\mathcal{B}}(h)$, must also conserve those quantities. In particular, such schemes will conserve the Poincaré invariants, the linear and angular momentum, and the infinite differentiability class order of the governing differential equations (Hernandez, 2019a), all of which are important for the accuracy of solutions (Hernandez, 2019b). Having these strong conservation properties built into symplectic schemes gives them excellent long-term behaviour.

The main subtlety, however, is that each split operator, e.g. ${\mathcal{A}}(h)$, is conserving its own Hamiltonian, e.g. $H_{A}$. It is therefore not clear what Hamiltonian an overall composition scheme like the one in Eq. 3 is conserving. Then again, a Hamiltonian $H_{A}$ yields the associated set of differential equations $\dot{{\bf z}}={\mathcal{L}}_{A}\>{\bf z}$ through Hamilton’s equations, so it should not be surprising that this answer can also be derived from the BCH formula.

One can show that one gets analogous formulae to Eqs. 5 and 7 (e.g., Yoshida, 1993),

[TABLE]

and

[TABLE]

where curly brackets denote Poisson brackets. Thus, each of the error terms in the modified differential equations (e.g., Eq. 5) can be derived from a corresponding error Hamiltonian (e.g., Eq. 8) (e.g., Saha & Tremaine, 1992). This is another way of seeing that the integration errors respect the symplectic geometry of the problem.

2.3 The Wisdom-Holman Map

While the energy errors in kinetic-potential splittings do not drift secularly with time, they remain large throughout the integration. Physically, this is because one is making drastic deviations from the true trajectory each timestep, i.e., force-free motion along ${\mathcal{A}}(h)$ alternating with order-unity kicks in the particle velocities under ${\mathcal{B}}(h)$.

In the case of weakly perturbed systems, like planets moving on nearly Keplerian orbits, it is much more powerful to instead split into an integrable dominant operator and the weak perturbation. Such a scheme, today known as the Wisdom-Holman map, proved a major development for planetary integrations, enabling the first integrations of the solar system over Gyr timescales (Wisdom & Holman, 1991).

While the development of the Wisdom-Holman map was not driven by this operator framework (see Wisdom 2018 for a historical perspective), it is instructive for our purposes to analyze it in this light. Instead of splitting the differential equations into components of comparable magnitudes like in Eq. 2, the Wisdom-Holman map splits them as (Wisdom & Holman, 1991; Kinoshita et al., 1990),

[TABLE]

where the ${\bf a_{i}^{\text{Kep}}}$ are the two-body Keplerian accelerations, and the ${\bf a_{i}^{\text{int}}}$ are the remaining interaction accelerations444These capture both the direct gravitational interactions between planets and their indirect effects on one another as they each pull on the central star. There are many ways to do this split between Keplerian and interaction accelerations (Hernandez & Dehnen, 2017; Rein & Tamayo, 2019). The splitting above into two operators corresponds to using Jacobi coordinates (Wisdom & Holman, 1991).. We have also introduced a factor of $\epsilon_{p}$ of order the characteristic planet-star mass ratio555This assumes that the distances from the planets to the central star are comparable to their interplanetary separations. In general, $\epsilon_{p}$ would capture the ratio of the perturbation accelerations to those of the dominant Keplerian motion. to keep track of the fact that the interplanetary accelerations (and therefore the changes in the interaction steps) are much smaller than the dominant Keplerian accelerations due to the central star.

Like in the kinetic-potential split, the corresponding $\mathcal{I}(h)$ is trivial, since it keeps the positions (and thus the accelerations) constant across the timestep. The exact solution for the velocities is then just

[TABLE]

By contrast, the Kepler step $\mathcal{K}(h)$ updates both the positions and velocities of each body along their respective unperturbed two-body orbits. Interestingly, more accurate and efficient schemes for solving the Kepler problem have continually been found for over three hundred years (e.g., Newton, 1687; Machin, 1737; Rambaut, 1890; Plummer, 1896; Brown, 1931; Danby, 1992; Mikkola & Innanen, 1999; Rein & Tamayo, 2015a; Wisdom & Hernandez, 2015).

2.4 Energy Error Estimates

As an example, consider a case of two Earth-mass planets with semi-major axes of 1 and 2 AU around a solar-mass star ($\epsilon_{p}=3\times 10^{-6}$) with orbital eccentricities of 0.01. We specialize to second order schemes (Eq. 6), and compare the kinetic-potential splitting (Eq. 2) to the Wisdom-Holman splitting (Eq. 10). Figure 1 plots both methods’ fractional energy errors over ten thousand inner-planet orbits, using a timestep of 1% the innermost planet’s orbital period.

The Wisdom-Holman map conserves the energy many orders of magnitude better than the kinetic-potential splitting, and we can make simple order of magnitude estimates to understand why. This is clearest in a Hamiltonian framework for conservative systems (e.g. Wisdom & Holman, 1991; Saha & Tremaine, 1992), but we instead show how it arises in the above operator-centred picture that will later generalize to dissipative systems. While this is not a rigorous derivation, it is useful to be able to quickly estimate expected errors and their associated scalings. While we do not plot all possible dependencies, we have verified that the various scalings predicted below are correct.

A useful feature of the abstract Lie derivative formalism of Sec. 2 is that they describe not only how the phase space variables ${\bf z}$ evolve, but also how functions like the energy $E({\bf z})$, change with time. In particular, just like $\dot{{\bf z}}^{\rm err}={\mathcal{L}}_{S}^{\rm err}{\bf z}$ yields the additional terms in the modified differential equations for the phase space variables using a particular split scheme ${\mathcal{S}}$ (Sec. 2), $\dot{E}^{\rm err}={\mathcal{L}}_{S}^{\rm err}E({\bf z})$ yields the additional terms in the modified differential equations for the energy evolution.

To approximately estimate the error $\Delta E^{\rm err}$ they cause, we can multiply the additional terms $\dot{E}^{\rm err}$ by the timescale over which these errors coherently add up. Toward this end, on the left of Fig. 1, the logarithmic time axis makes it possible to see that the errors accumulate coherently over several timesteps ($h=0.01$), and that on longer timescales, they oscillate. We can understand this from Sec. 2. Because the BCH formula says that all the error terms in the modified differential equations are composed of the split operators that propagate the real dynamics (Eq. 7), the errors will inherit the timescales in the original problem. In particular, the timescale over which errors will accumulate coherently will be the shortest timescale in the problem. In this case of near-circular, non-resonant orbits, $T$ would be the synodic period on which planets kick each others’ orbits every conjunction. In our typical example where the orbits are not too close, this is approximately the orbital period of the innermost planet; for a very eccentric orbit, $T$ would be the much shorter timescale for pericentre passage (e.g., Rauch & Holman, 1999; Rein & Tamayo, 2015b), etc. We thus have

[TABLE]

where, for the second-order schemes used in Fig. 1, ${\mathcal{L}}_{S}^{\rm err}$ is given by Eq. 7.

Since ${\mathcal{L}}_{A}\>$ and ${\mathcal{L}}_{B}\>$ each yield time derivatives according to their own subset of the differential equations, we can make a rough estimate in Eq. 12 by replacing each instance of ${\mathcal{L}}_{A}\>$ and ${\mathcal{L}}_{B}\>$ in ${\mathcal{L}}_{S}^{\rm err}$ by their respective inverse characteristic timescale. In this case we are referring specifically to the timescale on which each operator in isolation would change the particle states by order unity, i.e., $\text{z}/\dot{\text{z}}$.

For Kepler splittings, the Kepler Lie derivative ${\mathcal{L}}_{K}\>$ induces order-unity changes on the orbital timescale666More precisely $2\pi/T$. Considering the simple 1-D case of a circular orbit with semimajor axis $a$ and orbital frequency $\omega$, $\dot{\text{v}}=\omega^{2}a$, so ${\mathcal{L}}_{K}\>\sim\text{v}/\dot{\text{v}}=1/\omega$. We are ignoring these factors of $2\pi$, which largely cancel with the coefficients in the BCH formula, but they can be important for higher-order schemes that accumulate many such factors. $T$. Perturbation operators like the interaction step in the WH map induce order-unity changes on timescales that are longer by roughly a factor of the perturbation strength $\epsilon$ (e.g., for the interaction operator, $\text{v}/\dot{\text{v}}_{I}=(\text{v}/\dot{\text{v}}_{K})(\dot{\text{v}}_{K}/\dot{\text{v}}_{I})=\epsilon_{p}T$, see Eq. 10).

If we further assume that the two operators ${\mathcal{L}}_{A}\>$ and ${\mathcal{L}}_{B}\>$ do not commute, i.e., $[{\mathcal{L}}_{A}\>,{\mathcal{L}}_{B}\>]\sim{\mathcal{L}}_{A}\>{\mathcal{L}}_{B}\>$, then plugging Eq. 7 into Eq. 12 yields a fractional energy error estimate of $\sim\epsilon_{p}(h/T)^{2}$. This gives the right scaling for a second-order method, as expected. Indeed, the same argument applied to any Kepler splitting says that if the perturbation accelerations are a factor of $\epsilon$ smaller than the Keplerian ones, then the second-order splitting scheme given by Eq. 6 will yield energy errors $\sim\epsilon(h/T)^{2}$.

By contrast, in the kinetic-potential splitting (Eq. 2), both ${\mathcal{L}}_{A}\>$ and ${\mathcal{L}}_{B}\>$ induce order-unity changes on the orbital timescale $T$, so ${\mathcal{L}}_{A}\>\sim{\mathcal{L}}_{B}\>\sim 1/T$. This yields a correspondingly larger error $(h/T)^{2}$.

More physically, the WH method is benefiting from the fact that the dominant motion is known and solved exactly. In the limit of $\epsilon_{p}\rightarrow 0$, the Wisdom-Holman map would yield an exact trajectory, while the kinetic-potential splitting would still be splitting the Kepler problem and the errors would remain large and unchanged.

These simple estimates (black arrow in Fig. 1) match well to integrations with the Wisdom-Holman map, while our estimate for the kinetic-potential splitting is depressed by approximately a factor of $e$ from our simple estimate. This has to do with our simplification of the commutation relations, and has interesting broader implications that we discuss in Sec. 3.2.

In summary, while the Kepler step is more expensive to compute than the steps in a kinetic-potential splitting, the factor of $\epsilon_{p}$ gain makes reaching a given level of accuracy significantly faster with the Wisdom-Holman map. Perhaps more importantly, in cases where inter-planetary forces are very weak, the Wisdom-Holman method keeps the additional terms in the modified differential equations smaller than the interplanetary terms. By contrast, in the kinetic-potential splitting, these additional error terms can easily be larger than the inter-planetary forces, leading to spurious results.

These scalings suggest the cases where splitting methods will be most useful, i.e. when perturbations are weak $\epsilon\ll 1$ and when the shortest characteristic timescale $T$ is not too short. For example, for a planet experiencing dynamical tides on an extremely eccentric orbit, the strength of the perturbation changes by orders of magnitude over an orbital period. In this case, the shortest timescale would be the characteristic timescale of pericentre passage, which can be much shorter than an orbital period. Because simple symplectic schemes require a fixed timestep, the small $h$ required where the forces are strongest acts as a bottleneck for the remaining orbit that could otherwise be easily integrated. In such cases it is typically more efficient to use a more general scheme with adaptive timesteps like IAS15, or in some restricted cases it is possible to symplectically adapt timesteps in inverse proportion to the strength of an external potential (Preto & Tremaine, 1999; Mikkola & Tanikawa, 1999; Petit et al., 2019).

3 Additional Forces

We now move beyond point-particle gravity to consider adding additional effects in N-body integrations. We will see that the above framework, by not specializing to Hamiltonian systems, can be used to straightforwardly understand the numerical behavior under both conservative and dissipative forces.

3.1 Conservative Forces

We begin by considering conservative forces that can be derived from position-dependent potentials. Hamiltonian perturbations that cannot be written strictly in terms of the particle positions are also important (e.g., post-Newtonian corrections), but we defer their discussion to Sec. 5.

We note that while any position-dependent potential can obviously be trivially incorporated into a kinetic-potential split, it can also be directly inserted into the interaction step in the Wisdom-Holman map as an additional position-dependent acceleration in Eq. 10. Then Eq. 11 remains the exact solution with the acceleration given by the sum of all the position-dependent accelerations in the problem.

As an example, we consider the same two-planet case as above, but now with additional perturbations from an oblate primary, with the planets orbiting in the primary’s equatorial plane. The first corrections to point source gravity in a multipole expansion of the primary’s potential is the quadrupole term, with a potential in the equatorial plane of

[TABLE]

The innermost body will be most affected, so we take $a_{0}$ to be its original, reference semimajor axis, $G$ is the gravitational constant, $M$ and $R$ are the primary’s mass and radius, and $J_{2}$ is the standard dimensionless coefficient for the quadrupole field. We have also introduced a parameter $\epsilon_{J2}=\frac{1}{2}J_{2}\left(\frac{R}{a_{0}}\right)^{2}$ that captures the smallness of the effect relative to the dominant Keplerian potential $\approx GM/a_{0}$, since $a_{0}/r$ is approximately unity for nearly circular orbits. We take $\epsilon_{J2}=10^{-3}$.

We again compare second-order schemes (Eq. 6) using kinetic-potential and Wisdom-Holman Kepler splittings in Fig. 2. Comparing with Fig. 1, we see that the kinetic-potential errors have not changed. This is because the planetary accelerations in ${\mathcal{L}}_{B}\>$ for the kinetic-potential splitting (Eq. 2) are dominated by those due to the star. Adding a small acceleration due to the oblateness perturbation therefore does not noticeably change the dominant error terms in Eq. 7. By contrast, for the Kepler splitting, our chosen oblateness perturbations are much larger than the interplanetary ones ($\epsilon_{J2}\gg\epsilon_{p}$), so they change the error behavior. The same argument following Eq. 12 suggests an energy error $\sim\epsilon_{J2}(h/T)^{2}$.

However, both our error estimates for the kinetic-potential and Kepler splitting fall short by a factor of the eccentricity (both black arrows in Fig. 2 include this factor). This can be a significant suppression in the errors and is due to the geometrical properties of splitting methods, as we explore in the next section.

3.2 Splitting Methods Exploit Symmetries

Because for the remainder of the paper we focus on Kepler splittings, we consider in detail the eccentricity factor in the Wisdom-Holman errors seen in Fig. 2.

The Kepler problem is special in its degeneracies; it conserves all orbital elements, except for the mean longitude $\lambda$, which advances linearly. This implies that when alternating Kepler steps with perturbation steps, only perturbation components that change the particles’ semi-major axes (which determine the mean motions $d\lambda/dt$) do not commute.

For example, for a perturbation that only changed the eccentricity vector’s magnitude and direction, it would not matter whether the Kepler or perturbation step acted first, since the Kepler step leaves those quantities unchanged. The ordering only matters for perturbations that change the semimajor axis; if the perturbation acts first, the mean motion will change, and when the Kepler step operates, the particle will end up at a different mean longitude than if the Kepler step had acted first with the original value of the mean motion.

We can then see that for any purely radial force like the oblateness perturbation, the change in orbital energy (and thus semi-major axis) $-{\bf F}\cdot{\bf v}$ vanishes for circular motion. If we slowly increase the initial eccentricity, the perturbation will add and extract orbital energy as each body moves toward and away from pericentre, to leading order in proportion to $e$. So our mistake in estimating the energy error in Sec. 2.4 was in assuming that the commutators introduced factors of order unity, i.e., that $[{\mathcal{L}}_{A}\>,{\mathcal{L}}_{B}\>]\sim{\mathcal{L}}_{A}\>{\mathcal{L}}_{B}\>$. Our physical argument says that for radially perturbed, nearly circular orbits, the operators in fact approximately commute. The commutator $[{\mathcal{L}}_{A}\>,{\mathcal{L}}_{B}\>]$ oscillates with an amplitude $\sim e{\mathcal{L}}_{A}\>{\mathcal{L}}_{B}\>$ and vanishes at pericentre and apocentre where ${\bf F}\cdot{\bf v}=0$ and the semimajor axis is unchanged.777Actually it would not completely vanish in this problem due to the much weaker perturbations from the additional planet that would at these points become the dominant effect. This is physically what is causing the oscillations visible on the left of the top panel in Fig. 2, where the logarithmic time axis renders the orbital timescale visible. The planets were both started at pericentre, and one can see that the energy errors remain constant at pericentre (every time unit) and apocentre (every half-time unit). We note that the flat errors for $e<10^{-3}$ in the bottom plot of Fig. 2 are an artefact of our numerical setup888The physical argument is correct for all geometric eccentricities (i.e., the value measured by looking at the shape of the physical orbit in space, which includes the effects of the perturbation) but the typical subtlety arises when dealing with osculating eccentricities (i.e. the transformation from positions and velocities to unperturbed 2-body Kepler orbits, which ignores the effect of the perturbation) smaller than the size of the perturbation $\epsilon$. Because for convenience we initialize our orbits using osculating elements, decreasing the initial osculating eccentricities below $\epsilon_{J2}$ does not make them any more circular in geometric space, because the oblateness perturbations kick them away from Keplerian motion at order $\epsilon_{J2}$. So in our integrations, the geometric eccentricity that goes into the argument above reaches a floor at $\sim\epsilon_{J2}$, with a corresponding error $\epsilon_{J2}(h/T)^{2}e=\epsilon_{J2}^{2}(h/T)^{2}$, which matches the flat regime on the left of the bottom panel of Fig. 2..

The above example highlights the power of splitting methods’ geometrical properties. Because all the error terms in the modified differential equations (e.g., Eq. 7) consist of nested commutators of ${\mathcal{L}}_{A}\>$ and ${\mathcal{L}}_{B}\>$, any symmetries (or near-symmetries) are inherited in the integrator’s error properties to all orders. This means any higher-order splitting scheme that tries to correct for higher and higher order terms in the BCH expansion (e.g., Wisdom et al., 1996; Laskar & Robutel, 2001) will always have such a geometrical suppression in its leading error term.

This is not typically emphasized in the literature on symplectic N-body integration, presumably because the Kepler and interaction steps do not have any such symmetries (there are always conjunctions between planets that change the semi-major axes at leading order). However, many astrophysically important effects are highly symmetric potentials, e.g. multipole gravitational potentials, radiation pressure, simple general relativistic corrections, etc. This fact, combined with the the highly degenerate Kepler problem, often lead to strongly suppressed errors for the nearly circular or nearly coplanar orbits we often want to model.

In summary, splitting methods are often powerful not only for their long-term conservation of important quantities, but also because they can yield orders of magnitude higher accuracy for a fixed timestep. In the above case, a specialized second-order scheme split into separate Kepler and perturbation steps reduced the errors by a factor of $\epsilon_{J2}e=10^{-5}$ over short timescales compared to what one would obtain with a generic second-order Runge-Kutta scheme, even more over long timescales.

3.3 Dissipative Forces

While dissipative systems are qualitatively different from the conservative cases above, we now show that we can similarly understand their error behaviour using the above framework. As should be clear from the operator-centred development throughout the paper, we could take any differential equations, split them in any way we please, and if we could find time-evolution operators for those subproblems, then we could compose them and find their error behavior through the BCH formula. In particular, the error estimates we developed in Sec. 2.4 and applied in Sec. 3.1 come from the modified differential equations, without reference to conservative or dissipative forces.

On the other hand, one important distinction for symplectic systems is that the fact that the composition of two Hamiltonian operators must be Hamiltonian implies that all the BCH error terms are also Hamiltonian (Sec. 2.2). By contrast, when the perturbation is dissipative, the presence of ${\mathcal{L}}_{B}\>$ operators in the error terms of the modified differential equations (e.g., Eq. 7) shows that there is additional damping (or injection of energy) introduced by the splitting scheme. This generically leads to secular drifts, as we show in the following example.

Consider a single planet in orbit around its primary with an initial orbital eccentricity of $0.1$ and an orbital period of 1 year, subject to a simple damping force directed opposite to the planet’s velocity vector ${\bf v}$,

[TABLE]

where $m$ is the planetary mass. This parametrized prescription orbit-averages to yield inward migration with the semi-major axis decaying exponentially on an e-folding timescale $\tau_{a}$ (Papaloizou & Larwood, 2000). We set $\tau_{a}=1000$ years, and integrate for three damping timescales $\tau_{a}$. In this time, the planet’s semi-major axis moves inward by a factor of $\text{exp}(3)\approx 20$, and the orbital period decreases by a factor of $\text{exp}(9/2)\approx 90$.

We integrate the system using a second-order Kepler splitting

[TABLE]

where $\mathcal{K}$ evolves the planet on a Kepler orbit, and $\mathcal{D}$ damps the motion according to Eq. 14. We adopt a timestep of $10^{-3}$ times the planet’s initial orbital period, and choose to approximate $\mathcal{D}(h)$ by integrating across the timestep using a fourth-order Runge-Kutta scheme (see Sec. 5).

With dissipation, we no longer have a conserved quantity to track. Instead, as a proxy for the exact propagator $\mathcal{P}$, we also integrate the system with IAS15 (Rein & Spiegel, 2015), a high-order adaptive-timestep method whose accuracy reaches machine precision. We plot in orange in Fig. 3 the relative error between the energy calculated in our Kepler splitting integration and the one with IAS15.

One can see a clear secular drift. While the errors in the symplectic integrations in Fig. 2 oscillate and average out to yield flat time evolution, dissipative errors systematically overdamp or underdamp.

This scaling can be understood in the same way we analyzed the conservative case. The second-order splitting scheme error estimate is (Sec. 2.4) $\epsilon_{D}(h/T)^{2}$, with $T$ the planet’s orbital period, and $\epsilon_{D}=T/\tau_{a}$. However, like in the case of $J_{2}$ perturbations discussed above, there is additionally a geometric suppression of the error by a factor of the eccentricity999To see this we have to go one step further than the arguments in Sec. 3.2. There we argued that only the components of the perturbation that change the particles’ semi-major axes don’t commute with the Kepler step. Here the force is always pointed exactly opposite the particle’s velocity, so the two steps definitely do not commute. But now we additionally have to ask whether the order of operations matters for the quantity whose error we are measuring, i.e., the energy. For a circular orbit, it would not. Whether or not the Kepler step moves the planet along the circle before or after the step does not affect the energy loss, since the problem is azimuthally symmetric. So again the energy errors are suppressed by a factor of $e$, though in this case there would be no suppression of the phase errors. We consider phase errors more carefully in Rein et al. (2019a).. We show this estimate in Fig. 3 with a black dashed line.

The secular rise then simply reflects the fact that the orbital period is changing exponentially as the planet migrates inward, which yields a straight line on the linear-log scale.

Remarkably, this suggests that exponential outward migration would yield errors that exponentially decrease with time. We find that this is indeed the case, though this case is more subtle. In the case plotted above where the errors committed each timestep increase with time, the most recent errors always dominate the error budget. In the case where the instantaneous errors are getting smaller, one might expect the total error to remain at the level incurred at the beginning of the integration, where errors are largest. But these errors are oscillatory, and as long as the errors are changing adiabatically (i.e., as long as there are many orbits per migration timescale), the oscillations can march toward smaller amplitude. Of course for exponential outward migration, it does not take many migration timescales until the orbital period becomes comparable to the migration timescale, at which point the applied force in Eq. 14 (unphysically) becomes comparable to the central gravitational force, and we see numerical errors rise again.

In summary, the error behaviour in weakly dissipative splitting schemes can be understood in the same way as conservative cases, starting from the modified differential equations yielded by the BCH formula. Dissipative splitting schemes retain strong geometric properties (suppressing energy errors by a factor of $e$ above), and conservation properties. For example, a dissipative perturbation that damps the radial component of particle velocities (and thus acts radially) would conserve angular momentum. Composing such a step with conservative Kepler and interaction steps would still conserve the total angular momentum to machine precision (since each step does individually).

While we have seen that dissipative splitting schemes can systematically over(under)damp, we show in the following section that the same techniques from symplectic integration allow us to correct these errors at no additional computational cost.

4 Higher Accuracy Splitting Methods At Fixed Computational Cost

4.1 Symplectic Correctors

Ignoring any improvements resulting from any symmetries in the problem (Sec. 3.2), the splitting errors for a second-order Wisdom-Holman scheme should oscillate on orbital timescales with an amplitude of $\sim\epsilon(h/T)^{2}$. However, these high frequency error oscillations should be unimportant to the long-term evolution through the averaging principle (Wisdom & Holman, 1991).

Wisdom et al. (1996) go further and show that these high-frequency oscillations of $\mathcal{O}(\epsilon)$ can be efficiently removed through a near-identity canonical transformation, which they term symplectic correctors, from the real action-angle variables to mapping variables. In this picture, the short term oscillations intuitively arise from a mismatch in initial conditions. When going from the real system to the modified mapping, one has to correct the initial conditions to modified ‘mapping coordinates’. An integration is then of the form

[TABLE]

where each ${\mathcal{S}}{\mathcal{A}}{\mathcal{B}}{\mathcal{A}}(h)$ represents a timestep of the second-order Wisdom-Holman map (Eq. 6), and $\mathcal{C}(h)$ and $\mathcal{C}^{-1}(h)$ are the symplectic corrector transformations to and from mapping variables at the beginning and end of the integration. Wisdom et al. (1996) furthermore prove that this canonical transformation removes all the error terms in the modified differential equations101010We note that ordering matters once we specialize to systems where one operator is dominant. Wisdom et al. (1996) (with corrections in Wisdom 2006) give correctors for the $\mathcal{A}(h/2)\circ\epsilon\mathcal{B}(h)\circ\mathcal{A}(h/2)$ second-order scheme. The correctors for the alternative second-order scheme $\epsilon\mathcal{B}(h/2)\circ\mathcal{A}(h)\circ\epsilon\mathcal{B}(h/2)$ would have different coefficients. for the Wisdom-Holman map (Eq. 7) of $\mathcal{O}(\epsilon h^{n})$ for all $n$. The leading error using symplectic correctors is then $\mathcal{O}(\epsilon^{2}h^{2})$.

For generality and convenience of application, Wisdom et al. (1996) go on to show how the correctors ${\mathcal{C}}(h)$ for an arbitrary pair of integration steps $\mathcal{A}$ and $\epsilon\mathcal{B}$ (hereafter we assume that $\mathcal{A}$ is the dominant operator and $\mathcal{B}$ is a perturbation) can be approximated to progressively higher order through compositions of $\mathcal{A}$ and $\epsilon\mathcal{B}$ with carefully chosen timesteps forward and backward in time. Wisdom (2006) gives explicit coefficients for such $n$th order approximations $\mathcal{C}_{n}$, which remove error terms $\mathcal{O}(\epsilon h^{k})$ up to $k=n$ in Eq. 7. Depending on the problem, this leaves as the leading error the larger of the first uncorrected term $\mathcal{O}(\epsilon h^{n+2})$ and the first term at higher order in the perturbation parameter $\mathcal{O}(\epsilon^{2}h^{2})$.

We also briefly note that these correctors for the second-order Wisdom-Holman map ${\mathcal{S}}{\mathcal{A}}{\mathcal{B}}{\mathcal{A}}$ (Eq. 6) can also be used for the first-order scheme ${\mathcal{S}}{\mathcal{A}}{\mathcal{B}}$ (Eq. 3) with a simple modification. As far as the first-order scheme is concerned, the second-order scheme is a reasonable approximation to the exact propagator $\mathcal{P}$, so

[TABLE]

and thus ${\mathcal{C}}_{1}(h)=\mathcal{A}(h/2)$ acts as a first-order corrector for the first-order splitting scheme to yield ${\mathcal{S}}{\mathcal{A}}{\mathcal{B}}{\mathcal{A}}$. Then, if desired, one could use a higher order correctors ${\mathcal{C}}_{n}$ in Eq. 16. This trick of combining half steps of the second-order scheme in the middle of the integration is used widely (in REBOUND this is referred to as synchronization), but rarely described as a corrector that can be extended to higher order111111Jack Wisdom realized this long before we did (personal communication)..

Finally, the fact that correctors compensate for the leading errors provides a straightforward way to interpret what the scheme is getting wrong. For example, inverting Eq. 17 yields

[TABLE]

i.e., to leading order, the first-order scheme ${\mathcal{A}}(h)\circ{\mathcal{B}}(h)$ is equivalent to performing an exact step along the true solution $\mathcal{P}(h)$, except with the mistake of taking additional forward and backward Kepler half-steps before and after the fact.

4.2 Symplectic Correctors Are Splitting Correctors

While Wisdom et al. (1996) derived symplectic correctors through physically motivated canonical perturbation theory specific to Hamiltonian systems, they point out that the formulas could also have been derived solely from the Lie algebra (commutators) of the two operators using Lie series and the BCH formula. But as mentioned above, the BCH formula is a statement about the composition of non-commutative operators and has nothing to do with symplecticity.

‘Symplectic’ correctors therefore should correct any second-order splitting scheme involving a dominant operator and a perturbation. We now demonstrate explicitly that this is true even for dissipative perturbations. ‘Symplectic’ correctors are therefore more widely applicable than their name implies. Specifically, they represent ‘weak splitting correctors.’ In fact, like Wisdom et al. (1996), we rediscovered that similar ideas were applied as far back as Butcher (1969) to general Runge-Kutta methods without restriction to Hamiltonian systems.

In particular, we take the same problem of a single planet migrating inward from Sec. 3.3, except we shorten the semi-major axis damping timescale $\tau_{a}$ to 100 (initial) orbital periods, and integrate for one damping timescale.

We consider three integration schemes. First we apply the first order scheme $\mathcal{K}(h)\circ\mathcal{D}(h)$. From Secs. 2.4 and 3.3, this will yield errors $\sim\epsilon_{D}(h/T)e$. We then run a separate integration using first-order correctors (Eq. 17), whose error should scale as $\epsilon_{D}(h/T)^{2}e$, and an integration with third-order correctors ${\mathcal{C}}_{3}$ (Eq. 16), which should scale as $\epsilon_{D}(h/T)^{4}e$ for small $h$. The results are shown in Fig. 4, with all curves following the above estimates (black dashed lines).

Several features are apparent from Fig. 4. First, when thought of as a ‘symplectic’ corrector, it might seem surprising that symplectic correctors applied at the beginning and end of the integration (when the energy has dissipated by a factor of 3 in between) would give a more accurate result. But as argued above, in our case there is nothing specifically symplectic about them, they are simply weak splitting correctors. Second, while all the timesteps considered are small compared to the initial orbital period of 1 year, it is the final orbital period (vertical dashed line) that matters. We see that for timesteps $\gtrsim 10\%$ of the final orbital period, the error rises dramatically, a consequence of higher order terms in the BCH series becoming important.

It might seem counterintuitive that symplectic correctors can improve integrations involving dissipation. In particular, one might worry that any steps backwards in time in $\mathcal{C}(h)$ would not cancel with corresponding steps forward in time $\mathcal{C^{-1}(h)}$ when irreversible processes are involved. For example, a step forward with friction is not undone by a step backward in time—the second step would decrease the energy further. But this is not what is meant by, e.g., $\mathcal{A}(-h)$, i.e., it is not a step backward in time in the true sense where all the velocities are also flipped; rather, it merely means running the independent time variable run backwards, keeping everything else the same. Under this definition, a backwards step with frictional forces boosts the velocities, undoing a step forward in time, and removing the contradiction. In the terminology of Hairer et al. (2006), the scheme is time-symmetric, but not reversible. Hernandez & Bertschinger (2018) also find cases where time-symmetric, nonreversible schemes do not undergo energy drifts (see their Fig. 8).

The fact that ‘symplectic correctors’ can be derived directly from the BCH formula implies that they must be more general, given that the BCH formula itself makes no distinction between conservative or dissipative operators. Finally, the fact that the green curve in Fig. 4 reaches the theoretical round-off limit at $\sim 10^{-14}$ demonstrates that ‘symplectic correctors’ do improve integrations with dissipation to high accuracy.

This points out a limitation of splitting methods for dissipative systems. Splitting methods require a fixed timestep (though see, e.g., Preto & Tremaine 1999), and it must be shorter than the fastest timescale in the problem throughout the integration. If there is substantial dissipation, the required timestep might be so small that a high-order adaptive scheme like IAS15 is more efficient. We note that this concern also applies to the conservative N-body problem. For example, if Kozai cycles or other effects move planets onto very eccentric paths at late times, the timestep for the whole integration has to be chosen to match the shortest pericentre passage timescale that occurs in the integration121212There are ways to address this using time regularization..

5 Imperfect Operators

The above discussion assumes that an analytical form for each step can be found, so that the only errors come from the splitting scheme itself, e.g., Eq. 7. However, for many perturbations, analytical solutions are not known, which naturally leads to the question of how errors from numerically approximating the evolution across each step interact with the splitting errors discussed above.

As mentioned above, this is not an issue for position-dependent forces, since they can be trivially incorporated into the interaction step, Eq. 10. In that case Eq. 11 remains the exact analytical solution, using the accumulated accelerations from all the position-dependent effects.

However, many important astrophysical effects like the migration forces considered in the previous section, post-Newtonian corrections, tides, etc., are velocity-dependent. Previous authors have proposed incorporating such velocity-dependent perturbations in either the Kepler step (Malhotra, 1994), in the interaction step (Cordeiro et al., 1996), or as a separate step (Touma & Wisdom, 1994). The methods achieve comparable accuracy (Cordeiro et al., 1996), so we focus on the methods of Touma & Wisdom (1994) and Cordeiro et al. (1996), which fit directly into the framework discussed above.

In this case, the velocity-dependent accelerations vary across the perturbation step as the velocities change, so Eq. 11 is no longer exact; it is now merely an Euler step approximation, accurate only to first order in $h$. In this section we consider the errors that such an Euler step introduces, and compare it to higher-order numerical approximations to the evolution across the perturbation step in this general case where the propagator cannot be found analytically.

5.1 Hamiltonian Velocity-Dependent Forces

We first consider Hamiltonian velocity-dependent effects, which we will find are qualitatively different from ones involving dissipation. We note that while we will find such perturbations cause severe numerical problems for these previously proposed methods, those authors were interested in dissipative perturbations, which we consider in Sec. 3.3.

A good test case is the velocity-dependent first-order post-Newtonian correction for general relativity (Anderson et al., 1975), where, given the dominant central mass in planetary applications, we ignore second-order corrections of order the planet-star mass ratio. The equations of motion can be derived from a Hamiltonian (see Appendix B), whose conservation we can use to track the numerical accuracy.

In all integrations we adopt the first-order splitting

[TABLE]

where the corresponding differential equations are

[TABLE]

with the ${\bf a_{i}^{GR}}$ given in Eq. B.

We note that it is possible to find a splitting for this post-Newtonian Hamiltonian for which the evolution under each operator can be solved analytically (Saha & Tremaine, 1994). We nevertheless choose to use the above splitting both to explore the effects of imperfect approximations across the perturbation step, and with a view toward making the REBOUNDx library a general-purpose tool for integration. Such Hamiltonian velocity-dependent perturbations change the relationship between the particles’ physical velocities and their momenta through Hamilton’s equations, so this relationship varies depending on what forces are added to an integration. To minimize the necessary logic and possible pitfalls, we choose in REBOUNDx to always use the splitting 20, which can be applied to both conservative and dissipative velocity-dependent forces. This is the same setup as in both Touma & Wisdom (1994) and Cordeiro et al. (1996); we now compare their choice of integrating across the ${\mathcal{GR}}$ step using a first-order Euler approximation (Eq. 11) to using higher-order methods.

We first note that in this case of Hamiltonian velocity-dependent forces, Mikkola (1998) proposes using a low-order symplectic scheme like the implicit midpoint method to salvage the symplecticity of the scheme. However, even a symplectic scheme will not be symplectic when applied the non-canonical position and velocity coordinates of Eq. 20, like all the above authors consider131313As a simple analogy, consider a canonical (and thus symplectic) transformation that rotates 2D Cartesian coordinates. If we instead consider a transformation that first converts to non-canonical polar coordinates, naively applies the same symplectic rotation matrix to the polar coordinates, and then converts back to Cartesian coordinates, the sequence of transformations will not be canonical.. Therefore, one might not actually expect any advantage from using a symplectic scheme across the perturbation step.

We test this empirically on the K2-137 system (Smith et al., 2017), for which post-Newtonian effects are important. It consists of a 0.89 Earth-radii ($R_{\oplus}$) planet on a 4.3 hour ($a_{1}\approx 0.01$ AU) orbit around a $0.46$ solar-mass primary. We assign the planet the mass of the Earth. While in reality tides should have circularized the orbit, we inflate the initial orbital eccentricity to 0.01, in order to see the apsidal precession induced by the relativistic effects. The post-Newtonian perturbation strength is $\epsilon_{GR}=3GM/ac^{2}\sim 3\times 10^{-6}$ (e.g., Nobili & Roxburgh, 1986).

We compare integrating across the GR perturbation step with the second-order, symplectic implicit midpoint method advocated by Mikkola (1998) to three non-symplectic schemes: first-order Euler (Touma & Wisdom, 1994; Cordeiro et al., 1996), second-order Ralston’s Runge Kutta (RK2), and fourth-order Runge Kutta (RK4). In all cases the timestep is $\approx 8.1\%$ of the innermost planet’s orbital period. This implies a relative error of $\sim\epsilon_{GR}(h/T)e\sim 2\times 10^{-9}$ (black arrow).

In the RK2, RK4 and implicit midpoint integrations in Fig. 5, the truncation errors remain much smaller than the splitting scheme errors, and the simulations are visually indistinguishable, showing flat, oscillatory errors. As argued above, applying a symplectic algorithm like the implicit midpoint method to non-canonical variables does not yield a symplectic scheme, so none of these integrations are exactly symplectic. Nevertheless, by ensuring that the truncation errors across the perturbation step are much smaller than the splitting scheme errors, we only see the oscillatory symplectic behavior.

Of course, the most salient feature in Fig. 5 is that the error in the Euler method of Touma & Wisdom (1994) and Cordeiro et al. (1996) (blue) error grows secularly, and in fact the planet is spuriously ejected after approximately 5.8 billion orbits, or in under 3 million years. This behaviour too can be analytically estimated.

5.2 Truncation Errors

For general perturbations that can not be integrated analytically, the local truncation error across the perturbation timestep must be incorporated into the error budget.

The local truncation error $\Delta z_{LT}^{\rm{err}}$ across a single timestep using a one-step, order-$n$ method comes from the first neglected term at order $n+1$ in a Taylor series approximation to the solution $z=f(t)$ (we focus on a single scalar coordinate $z$ for simplicity),

[TABLE]

where $f^{n+1}$ is the $n+1$th derivative of the true solution across the perturbation step, and $\xi$ is an unspecified time in the range $[t,t+h]$. If we assume a solution $f(t)=z_{0}\exp(i2\pi t/T_{pert})$, with the characteristic perturbation timescale $T_{pert}$ a factor of $\epsilon$ longer than the orbital timescale $T$, we have

[TABLE]

This is the error incurred every timestep. We can obtain the worst-case relative global truncation error after an integration time $t$ by assuming the local errors add coherently, and multiplying Eq. 22 by the number of steps $t/h$,

[TABLE]

We overplot this estimate for the Euler method of Touma & Wisdom (1994) and Cordeiro et al. (1996) ($n=1$) in Fig. 5 and see it gives a good match to the numerical behaviour. The truncation errors for the second-order RK2 method would be lower by a factor of $\epsilon_{GR}(h/T)\sim 10^{-7}$ and thus never become visible.

We might, however, expect different behavior under weak dissipation. To reach Eq. 23 we assumed the worst-case scenario that the one-step errors added coherently. At a more detailed level, one must consider how such errors are propagated by the dynamical flow itself. Consider a set of trajectories in phase space within a differential error volume around the true trajectory. By the divergence theorem, the fractional growth of this volume element across a timestep is given by the product of $h$ and the divergence of the vector field of trajectories $\nabla\cdot{\bf\dot{z}}$. In a Hamiltonian system like above, this divergence vanishes through Hamilton’s equations (this is Liouville’s theorem), so, at best, the dynamics are neutral for the accumulation of errors141414In reality, an initial parcel of trajectories within some error volume would get sheared out at fixed volume.

By contrast, under weak dissipation, the divergence of the flow is $\mathcal{O}(\epsilon)$ and negative, so errors contract by $\mathcal{O}(\epsilon h)$ each timestep as the dynamics brings nearby trajectories together. One might therefore expect that integrating across the perturbation step with even a first-order Euler method would be good enough, given that one adds one-step errors $\sim\mathcal{O}(\epsilon^{2}h^{2})$ (Eq. 22) slower than they are damped by the dynamics. We test this empirically in the following section.

5.3 Dissipation

Taking a somewhat realistic damping problem, we initialize two Jupiter-mass planets in a 3:2 mean-motion resonance around a Sun-like star, with the innermost planet initially at 0.90 AU using celmech 151515celmech is a publicly available code for semi-analytic manipulations of the classical disturbing function in celestial mechanics. Among other things, it can perform transformations between orbital elements and resonant variables: https://github.com/shadden/celmech. We then apply an exponential eccentricity damping force prescription (Papaloizou & Larwood, 2000) to both planets with an e-folding decay timescale of $\tau_{e}\approx 1000$ initial inner-planet orbits using the modify_orbits_forces implementation in REBOUNDx. The outer planet is then additionally acted on by an inward migration force (Eq. 14) with a timescale 100 times longer $\tau_{a}\approx 0.1$ Myr. The eccentricities grow as the system evolves deeper into resonance until the migration forces pushing the system into resonance are balanced against the dissipative divergence of the orbits caused by the eccentricity damping (e.g., Batygin & Morbidelli, 2012; Lithwick & Wu, 2012; Goldreich & Schlichting, 2014), reaching equilibrium eccentricities of 0.022 and 0.025 for the inner and outer planet, respectively.

In Fig. 6 we then plot the energy errors during the subsequent integration, where the eccentricities remain in a quasi-steady state as the two plants migrate inward together, on a timescale approximately two times longer than the $\tau_{a}$ estimate above (two equal mass planets have to be moved). We use the second-order scheme161616We provide a simple argument why one can always use such a time-symmetric splitting with multiple operators to create a second-order scheme. The timestep coefficients for each operator must sum to unity in order to match the modified differential equations to the true ones at zeroth order, and by writing an explicitly time-symmetric scheme, we ensure that the additional terms in the modified differential equations cannot depend on the sign of $h$ (e.g., Eq. 7). The scheme is therefore second-order by construction.

[TABLE]

and integrate across each damping step using the four methods of Sec. 5.1. We integrate for $4\times 10^{5}$ initial orbits, which translates to approximately 4000 eccentricity-damping timescales, or 2 semimajor-axis damping timescales. We adopt a small fixed timestep of $10^{-4}$ of the innermost planet’s initial orbital period to accommodate the shrinking orbital periods, which decrease by a factor of $\approx 25$ over the integration. As above, we compare the energies to the energies in a ‘perfect’ integration with IAS15.

In this case the scheme errors are dominated by the interplanetary, resonant interactions. As in Sec. 3.3, we see in Fig. 6 that the energy errors grow due to the shrinking orbital periods. Like in Fig. 5, the RK2, RK4 and implicit midpoint integrations are visually indistinguishable. By contrast, while the integrations using the first-order Euler method across the damping step have higher errors, it is not a runaway effect as in the conservative case in Fig. 5. As discussed in Sec. 5.2, this is because of the dissipative dynamics continually damping any accumulated energy errors.

In summary, the crudest, most computationally efficient Euler step across a velocity-dependent force step should give reliable results for dissipative forces, as advocated by Touma & Wisdom (1994) and Cordeiro et al. (1996). However, such a scheme would yield qualitatively wrong results in long-term integrations with conservative velocity-dependent forces (Sec. 5.1). This can straightforwardly be solved by integrating across the perturbation step with a higher method; for high-accuracy explicit splitting schemes that go beyond this to correct even such higher-order integrations across the perturbation step, see Blanes et al. (2013).

6 REBOUNDx

In the previous sections we covered several technical aspects of how additional forces interact with splitting integration schemes like WHFast. However, many readers are likely interested in more practical questions of how to quickly get something working. We therefore implement the integration schemes and ideas discussed above in a new library, REBOUNDx, which can be found at https://github.com/dtamayo/reboundx.

REBOUNDx is library that provides tools and routines to accurately and self-consistently incorporate many different astrophysical effects to $N$-body simulations. It seamlessly interfaces with the REBOUND N-body integration package (Rein & Liu, 2012). Like REBOUND, REBOUNDx is written in C, and provides a PYTHON wrapper to interface with other libraries. Also, like REBOUND, the REBOUNDx source code is machine independent. We implement a binary format to save REBOUNDx configurations that interfaces with the SimulationArchive class in REBOUND, making it possible to share and reproduce results bit by bit.

We have already implemented several common forces in REBOUNDx, some of which we used as examples in the discussion above. We hope to increase the number of supported astrophysical effects over time. REBOUNDx is an open-source project and we welcome the community to contribute new effects. We provide many tutorials and examples in the form of Jupyter notebooks that illustrate the library’s usage and how one can add new effects with minimal effort. We also share the scripts used to generate all the plots in this paper in a separate repository at https://github.com/dtamayo/reboundxpaper.

We now give an overview of the main concepts, structures and logic in the REBOUNDx package.

6.1 Which Integrator To Use?

In a scenario where all the forces (including interplanetary forces) are always small compared to the gravitational forces from the central body, we recommended using Wisdom-Holman Kepler splitting schemes. The implementation of this method in REBOUND is WHFast (Rein & Tamayo, 2015a). This splitting will typically be fastest at reasonable levels of accuracy, and it provides strong conservation properties in long-term integrations.

Higher levels of accuracy are supported in cases where the forces can all be derived from position-dependent potentials by setting sim.ri_whfast.corrector $=N$, where $N$ is the order of the corrector (Sec. 4.1). In velocity-dependent cases, correctors can be implemented manually (the notebook fig4.ipynb in this paper’s repository provides an example). For even higher levels of accuracy (without significant computational cost), one can use the WHCKL method (Rein et al., 2019b), which is an implementation of the kernel method of Wisdom et al. (1996) (see also Wisdom, 2018) that also removes errors at $\mathcal{O}(\epsilon^{2})$.

In cases with close encounters between planets, the motion will no longer be nearly Keplerian, and Kepler splittings will yield large errors. In cases where such close encounters are rare, hybrid integrators are powerful tools (Duncan et al., 1998; Chambers, 1999; Rein et al., 2019c). In REBOUND, the MERCURIUS integrator retains efficiency by using a Kepler splitting when possible, and switches to an adaptive, high-accuracy integration with IAS15 for particles undergoing close encounters (Rein et al., 2019c). MERCURIUS supports additional forces, but only switches to IAS15 during close gravitational encounters. In the current implementation, any additional forces must remain small throughout the integration when MERCURIUS is used.

When the above conditions are not met, the safest option (and the default) in REBOUND is to use the IAS15 integrator (Rein & Spiegel, 2015). Because it is an adaptive scheme, IAS15 can handle arbitrary forces and still yield solutions that are accurate close to the machine precision. IAS15 will typically be slower in cases where WHFast and MERCURIUS can be used. However, IAS15 will always be more accurate and is applicable to a wider variety of problems171717 By using adaptive timesteps, IAS15 can even be faster in cases, e.g., very eccentric orbits, where one has to choose a small timestep with WHFast or MERCURIUS in order to resolve a short timescale in the problem (e.g., pericenter passage)..

6.2 REBOUNDx Structures

The two principal ways to incorporate additional effects into an $N$-body simulation are through forces and operators.

6.2.1 Forces

One way to implement an additional effect is to write it explicitly as a force. This force then contributes to the accelerations of all particles, in addition to the standard Newtonian gravity.

When the IAS15 integrator is used, the integrator automatically adapts the timestep such that integrations will be accurate close to machine-precision. If one of the various splitting integrators available within REBOUND (for example WH, MERCURIUS, and LEAPFROG) is used, then the additional forces are accumulated and applied in the integrator’s interaction (or kick) step using Eq. 11. As discussed above, this is appropriate for position-dependent or dissipative forces (Secs. 3.1 and 5.3). For conservative, velocity-dependent forces, one should add a separate operator (Sec. 6.2.2).

For the hybrid symplectic MERCURIUS integrator (Rein et al., 2019c), only forces derivable from position-dependent potentials are supported.

6.2.2 Operators

Rather than calculating accelerations which are then integrated in the interaction or kick step of the integrator, operator steps yield solutions. Operators update the positions and velocities (or even masses or other particle parameters) to the value they should have at the end of a specified timestep.

Operators (e.g., ${\mathcal{D}}$) together with a timestep (e.g., $h/4$) form a step (e.g., ${\mathcal{D}}(h/4)$). The same operator can thus be used in multiple steps with different timesteps to construct higher-order splitting schemes (e.g., Eq. 24).

For cases where the solution for the positions and velocities at the end of the step is not known, or difficult to calculate analytically, the REBOUNDx object integrate_force operator can take a force, and integrate it across the timestep using any of the four integrators discussed in Sec. 5. The default is RK4. This is the recommended method for conservative, velocity-dependent forces (Sec. 5.1).

6.3 Ordering

All forces and operator steps are stored in linked lists such that they are executed in the reverse order they were added. Thus, to get a scheme consisting of the steps $\mathcal{K}(h)\circ\mathcal{D}(h)$, the operator $\mathcal{K}$ should be added before $\mathcal{D}$, i.e. left to the right, since the last step acts first.

6.4 Timestep logic

REBOUNDx provides a simple API for adding effects into N-body integrations from the list of implemented options, demonstrated in the numerous examples in the documentation. However, it also provides lower-level functionality for customized tasks. The three main function pointers in a REBOUNDSimulation object that REBOUNDx sets and calls behind the scenes every timestep are

• •

additional_forces. This iterates through the linked list of attached forces to update particle accelerations.

• •

pre_timestep_modifications. This iterates through the linked list of attached operator steps to execute before the main REBOUND step.

• •

post_timestep_modifications This iterates through the linked list of operator steps to execute after the main REBOUND step.

The user populates these lists by either adding forces and operators already implemented in REBOUNDx, or through user-defined functions in either PYTHON or C. By default, adding an operator to an integration adds half steps before and after the main REBOUND step to make a second-order scheme.

In Fig. 7, we schematically outline a full REBOUND timestep using the WHFast integrator. Different integrators switch out the reb_integrator_part1 and reb_integrator_part2 steps, but in all cases they rely on the gravitational and additional forces calculated in between them.

6.5 Coordinate systems

In contrast to many other N-body packages, REBOUND allows the user to set up the simulation in an arbitrary inertial frame. The coordinates are returned to the user in the same inertial frame.

The various splitting integrators in REBOUND use a number of different coordinates internally (e.g. Jacobi coordinates). In order to minimize logic and pitfalls, all REBOUNDx functions are always called after converting the particles’ positions and velocities back to Cartesian coordinates in the inertial frame. Force and operator implementations thus do not need to worry about indirect or non-inertial forces. This also implies that it is typically desirable to ensure that the net applied forces vanish. This is discussed in detail in Appendix C.

6.6 Parameters

Additional effects in an $N$-body simulation typically come with several associated parameters (e.g., a migration timescale for a damping force). REBOUNDx provides an interface which allows a user to attach an arbitrary numbers of parameters to forces, operators, and particles. The implementation is such that parameters remain tied to their parent structures even as objects get removed or reshuffled in the simulation.

6.7 Implemented effects

A current list of implemented effects can be found in the documentation. At the time of this writing, the following effects are available:

• •

Migration forces (Papaloizou & Larwood, 2000)

• •

Eccentricity and inclination damping forces (Papaloizou & Larwood, 2000)

• •

General relativistic 1 PN corrections (Newhall et al., 1983)

• •

General relativistic 1 PN corrections for a dominant central mass (Anderson et al., 1975)

• •

General relativistic 1 PN corrections with a simplified potential (Nobili & Roxburgh, 1986)

• •

Radiation Forces (Burns et al., 1979)

• •

Mass loss

• •

Precession from conservative equilibrium tides (Hut, 1981)

• •

General central forces

• •

Higher order gravitational moments, J2 and J4

7 Conclusion

Throughout this paper, we have developed the error behaviour of splitting schemes in a general, non-Hamiltonian framework. This approach clarifies how several properties typically ascribed to symplectic integrators also carry over to weakly dissipative systems.

Many astrophysically relevant effects have strong symmetry properties. We showed explicitly that the remarkable degeneracy of the Kepler problem makes schemes that split the evolution into Kepler steps and perturbation steps exploit these symmetries to often strongly suppress energy errors in the integration, as in the cases of radial or drag forces applied to nearly circular orbits (Sec. 3.2). We also showed that so-called ‘symplectic correctors’, which reduce energy errors by orders of magnitude at fixed computational cost (Wisdom et al., 1996), apply equally well to weakly dissipative systems and can thus be more generally thought of as ‘weak splitting correctors’ (Sec. 4.2).

We also explored the general case where individual operator steps can’t be solved analytically, and one is forced to integrate across the perturbation step. We showed that the previously advocated approaches of incorporating additional forces into the Wisdom-Holman map (Touma & Wisdom, 1994; Malhotra, 1994; Cordeiro et al., 1996) work well for dissipative effects (Sec. 5.3), but give qualitatively wrong answers for conservative, velocity-dependent forces like post-Newtonian corrections of general relativity (Sec. 5.1).

Finally, we described REBOUNDx, an open-source C library for incorporating additional effects into REBOUND N-body integrations, together with a convenient PYTHON wrapper. Users can either choose to add an already implemented astrophysical forces, or easily implement new effects. We hope that this modular, open-source framework for robustly incorporating new effects into a variety of N-body integration schemes will help the community to run more accurate and realistic simulations. We encourage others to contribute to this library of astrophysical effects.

Acknowledgments

We would like to thank John Chambers for an insightful and constructive review that greatly improved this manuscript. We would also like to thank Jack Wisdom and Scott Tremaine for helpful discussions. 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 has been supported by the NSERC Discovery Grant RGPIN-2014-04553 and the Centre for Planetary Sciences at the University of Toronto Scarborough. 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 Operator theory

While not strictly necessary for the development in the main text, we review here some of the gaps glossed over in Sec. 2. Using the same notation, we seek a solution to the differential equations

[TABLE]

where here we define the Lie derivative more carefully. The Lie derivative with respect to a set of differential equations $\dot{{\bf z}}=\hat{F}{\bf z}$ acts on functions $g({\bf z})$ to return their total derivative. Since throughout the paper we treat only cases that are explicitly time-independent, we have

[TABLE]

where the $\dot{z}_{i}$ are given by $\hat{F}{\bf z}$.

With this definition, we can Taylor expand the solution around the current time ,

[TABLE]

As a concrete example, consider the simple one-dimensional differential equation $\dot{z}=z^{2}$. We have ${\mathcal{L}}_{P}z=z^{2}$, ${\mathcal{L}}_{P}^{2}z={\mathcal{L}}_{P}z^{2}=2z^{3}$ etc.. Plugging into Eq. 27 yields the Taylor expansion of the exact solution $z(t)=(1/z(0)-t)^{-1}$.

By analogy to the Taylor expansion for the exponential, this is often written as181818Throughout the paper we assume that the differential equations do not depend explicitly on time. If they did, an integral would be required in the exponential of Eq. 28.

[TABLE]

In words, $\mathcal{P}(h)$ is a time evolution operator, i.e., an integrator, that updates the phase space variables across a timestep. We then proceed to split the differential equations into two pieces $\dot{{\bf z}}={\mathcal{L}}_{P}\>{\bf z}=({\mathcal{L}}_{A}\>+{\mathcal{L}}_{B}\>){\bf z}$ as in the main text.

One subtle notational issue that has led to some errors in the literature is that the BCH formula is typically written as

[TABLE]

with

[TABLE]

However, $e^{h{\mathcal{L}}_{A}\>}e^{h{\mathcal{L}}_{B}\>}$ is not equivalent to the splitting scheme one might naively implement on a computer $e^{h{\mathcal{L}}_{A}\>}\circ e^{h{\mathcal{L}}_{B}\>}$.

In $e^{h{\mathcal{L}}_{A}\>}e^{h{\mathcal{L}}_{B}\>}{\bf z}$, $e^{h{\mathcal{L}}_{A}\>}$ acts on a function $g({\bf z})=e^{h{\mathcal{L}}_{B}\>}{\bf z}$, and therefore introduces partial derivatives of $g$ through Eq. 26. By contrast, the integrator $e^{h{\mathcal{L}}_{A}\>}\circ e^{h{\mathcal{L}}_{B}\>}$ first evaluates an intermediate ${\bf z}^{\prime}=e^{h{\mathcal{L}}_{B}\>}{\bf z}$, and these updated values are used to yield a final answer of $e^{h{\mathcal{L}}_{A}\>}{\bf z^{\prime}}$. The BCH formula applies to the former sense (one can straightforwardly expand the three exponentials in Eq. 29 to obtain the first terms in Eq. 30), but this interpretation is not particularly useful for writing integrators. Under iterated mappings, the repeated chain rules would quickly become unwieldy. By contrast, $e^{h{\mathcal{L}}_{A}\>}\circ e^{h{\mathcal{L}}_{B}\>}$ is straightforward to implement on a computer, and fortunately can be simply related to the scheme in the BCH formula to analyze its error properties. This is sometimes referred to as the Vertauschungssatz (Gröbner, 1967) (see Lemma 5.1 in chapter III of Hairer et al. (2006)); we provide a simple derivation and example here.

Looking back at Eq. 27, one could just as well have Taylor expanded any scalar function $f({\bf z})$, and concluded in Eq. 28 that the time evolution operator $e^{h{\mathcal{L}}_{A}\>}$ also evolves $f({\bf z})$ from time $t$ to time $t+h$ (with $f({\bf z})={\bf z}$ in Eq. 27 as a special case). Thus, since Eq. 29 is of the form $e^{h{\mathcal{L}}_{A}\>}f({\bf z})$, with $f({\bf z})=e^{h{\mathcal{L}}_{B}\>}{\bf z}$, we must have

[TABLE]

i.e., evolving $f({\bf z})$ across the timestep (LHS) is equivalent to first evolving ${\bf z}$ and applying the function $f$ to the updated variables (RHS). Therefore,

[TABLE]

This is the reason for the sign flips at all odd orders in Eq. 5.

As one might expect, Eq. 32 has led to some confusion in notation and sign errors in the literature (e.g., footnote 3 in Saha & Tremaine 1992 and Eq. 10 in Hernandez & Dehnen 2017). However, one reason why this is probably not typically noticed or pointed out is that most integrators used in the astrophysics literature are time reversible, and thus will not have any sign flips in their corresponding BCH formulas given that they are invariant to flipping the order of operations (e.g., Yoshida 1993 and the detailed analysis in Hernandez & Dehnen 2017).

As an explicit example of the above, consider a non-dimensionalized simple harmonic oscillator with position $x$, momentum $v$, and ${\bf z}=\left(x,v\right)$. The Hamiltonian is $H=v^{2}/2+x^{2}/2$, and we use a kinetic-potential splitting $H_{T}=v^{2}/2$ and $H_{V}=x^{2}/2$. The corresponding differential equations are given by Hamilton’s equations,

[TABLE]

The propagators can easily be found by integrating their respective differential equations, or by noting ${\mathcal{L}}_{T}\>^{n}{\bf z}={\mathcal{L}}_{V}\>^{n}{\bf z}=0$ for $n\geq 2$, so

[TABLE]

We first evaluate the composition used in the ${\mathcal{S}}{\mathcal{A}}{\mathcal{B}}$ integrator (Eq. 3), $e^{h{\mathcal{L}}_{T}\>}\circ e^{h{\mathcal{L}}_{V}\>}{\bf z}$, where $e^{h{\mathcal{L}}_{T}\>}{\bf z}^{\prime}$ is evaluated at ${\bf z}^{\prime}=e^{h{\mathcal{L}}_{V}\>}{\bf z}$

[TABLE]

By Eq. 32, this should be equivalent to

[TABLE]

which indeed matches Eq. 38.

Appendix B First Order Post-Newtonian Corrections

General relativistic corrections to Newtonian dynamics are detailed throughout the literature, but since we use these effects as one of the main tests for the various schemes in this paper, we present the equations and details of our implementation for completeness.

The first order post-Newtonian corrections (1PN) are $\mathcal{O}(v^{2}/c^{2})$, where $v$ is the characteristic orbital velocity and $c$ is the speed of light. We have implemented various levels of approximation. In gr_full, we include the full (1 PN) equations of motion, which are implicit and computationally expensive, but would be appropriate for equal-mass bodies.

For planets around a single star, one can typically ignore corrections to the GR perturbation of order the planet-star mass ratio. This significantly simplifies the equations, but involves velocity-dependent accelerations. This approximation is implemented in the gr effect in REBOUNDx, and is the case discussed in Sec. 5.

As a final level of approximation, one is typically interested in the apsidal precession effect from general relativity, in which case it is possible to write down a potential only depending on the particles’ position, which, in an orbit-averaged sense gives the correct precession rate (Nobili & Roxburgh, 1986). This is implemented in gr_potential, and has the advantage that it can be integrated analytically (Sec. 3.1), and is therefore computationally cheapest. However, this approximation introduce errors in the phases and instantaneous elements of $\mathcal{O}(v^{2}/c^{2})$. Typically, gr_potential is sufficient for planetary applications around a single star since it captures the correct secular behavior, but both gr_potential and gr are only appropriate for single-star systems. For higher multiplicity systems, one must use gr_full.

Naoz et al. (2013) provide the full Hamiltonian for a triple system (see also Schäfer, 1987). The equations of motion are given by Newhall et al. (1983) and Benitez & Gallardo (2008). These accelerations and Hamiltonian are implemented in gr_full.

In the gr implementation, we assume a dominant central mass, and ignore additional corrections of order the planet-star mass ratio. We outline our implementation for this effect in more detail since it differs from previous works, and is the subject of Sec. 5. The full Hamiltonian can be written as a sum of the Newtonian and 1PN Hamiltonians, $\mathcal{H}=\mathcal{H}_{N}+\mathcal{H}_{PN}$, where

[TABLE]

and (Saha & Tremaine, 1994),

[TABLE]

Having already ignored terms of $\mathcal{O}(m_{i}/m_{0})$, the difference between barycentric, Jacobi or democratic heliocentric coordinates (and Jacobi or physical masses) in $\mathcal{H}_{PN}$ is negligible. Interpreting all positions and momenta as Jacobi coordinates considerably simplifies the equations of motion, since the kinetic term in $\mathcal{H}_{N}$ remains a diagonal sum of the momenta. Application of Hamilton’s equations yields (e.g., Saha & Tremaine, 1994),

[TABLE]

where

[TABLE]

and the $\mathbf{\tilde{v}_{i}}\equiv\mathbf{p_{i}}/m_{i}$ are pseudo-velocities not equal to the physical velocities $\mathbf{\dot{r}_{i}}$. We note that Eq. B would have additional terms in democratic heliocentric coordinates due to kinetic cross terms in $\mathcal{H}_{N}$. We then have

[TABLE]

where

[TABLE]

$\mathbf{\tilde{a}_{i}}$ are the net Newtonian accelerations on particle $i$, and

[TABLE]

We see that the accelerations (Eq. B) are explicitly velocity-dependent, and must be treated with care in symplectic schemes (Sec. 5). As expected, $\mathbf{\ddot{r}_{0}}=0$ and the centre of mass moves at constant speed. Equation B needs to be solved implicitly for the pseudo-velocities $\mathbf{\tilde{v}_{i}}$, but this is easily accomplished iteratively, since the 1PN approximation is only appropriate when the perturbation is weak, and the physical and pseudo-velocities are approximately equal. Operationally, if like in REBOUND, the integrator does not work in Jacobi coordinates, one simply calculates the Newtonian accelerations in inertial coordinates, transforms to Jacobi coordinates, calculates the Jacobi accelerations $\mathbf{\ddot{r}_{i}}$ and transforms back to inertial coordinates. Rein & Tamayo (2015a) provide unbiased transformation algorithms between inertial and Jacobi coordinates.

Appendix C Centre of Mass

When inserting additional effects in N-body simulations in an inertial frame, it is often valuable to ensure that there are no net forces acting on the system’s centre of mass (COM). For example, the Wisdom-Holman integration scheme (Kinoshita et al., 1990; Wisdom & Holman, 1991; Rein & Tamayo, 2015a) assumes that the COM moves at constant velocity and removes that degree of freedom from the problem. This effectively means that any residual net force is incorrectly distributed among the particles’ relative coordinates. However, even integration schemes that track the COM degree of freedom can be adversely affected by finite forces on the barycentre. Such a simulation’s accuracy will continually deteriorate as the system moves away from the origin and the precision on the inter-particle separations degrades due to subtracting ever larger nearly equal numbers. These errors may be acceptable depending on the application, but can slow down schemes like IAS15 (Rein & Spiegel, 2015) that try to reach machine precision as they fail to converge with progressively smaller timesteps.

It is therefore numerically desirable to self-consistently model all back-reactions from additional effects so that the net force vanishes. While this is trivial for forces between pairs of particles in the simulation, one does not always wish to fully model all components.

Consider a net force $\mathbf{F_{D}}$ from a protoplanetary disk of mass $M_{D}$ on a planet. In response, the disk should feel an equal and opposite net force, but there is no disk “particle” in the simulation on which to apply it. This uncompensated force will yield an acceleration on the COM. However, the disk is tightly coupled to the central star gravitationally, suggesting that this back-reaction should be effectively communicated to the central body. More quantitatively, if we take a characteristic orbital radius for the disk $r_{D}$, then we can define a characteristic timescale for the back-reaction $\tau_{BR}\sim(r_{D}/a_{D})^{1/2}$, where $a_{D}$ is the net acceleration of the disk $F_{D}/M_{D}$. At the same time, the disk is gravitationally coupled to the primary on the Keplerian orbital period $\tau_{K}=2\pi(GM/r_{D}^{3})^{-1/2}$. If $\tau_{K}\ll\tau{BR}$, the central body will respond adiabatically to the accelerations of the disk, and the back-reaction could instead be added to the central star.

In other cases it might be better motivated to add the accelerations to the centre of mass of all or a subset of the particles, e.g., back-reactions onto a circumbinary disk. This corresponds to applying an acceleration $\mathbf{a}=\mathbf{-F}/M_{tot}$ to the appropriate set of particles, where $M_{tot}$ is their total mass.

In the case of forces that only depend on particle positions, the above choices are straightforwardly reflected in the structure of the Hamiltonians that govern the dynamics. Being able to precisely calculate these Hamiltonians provids a practical check on the numerical accuracy of the implementation for conservative systems. A Hamiltonian that depends on the radial positions $\mathbf{r_{i}}$ and $\mathbf{r_{j}}$ of two particles, $\mathcal{H}(\mathbf{r_{i}}-\mathbf{r_{j}})$ exerts equal and opposite forces on the two bodies. Similarly, one can show that for a Hamiltonian dependent on the distance from the system’s barycentre $\mathcal{H}(\mathbf{r_{i}-r_{COM}})$, where $\mathbf{r_{COM}}=M_{tot}^{-1}\sum_{j}m_{j}\mathbf{r_{j}}$, Hamilton’s equations dictate that if particle $i$ feels force $\mathbf{F_{i}}=-\nabla_{r_{i}}\mathcal{H}(\mathbf{r_{i}-r_{COM}})$, all particles (including $i$) should feel an additional acceleration $-\mathbf{a_{COM}}=-\mathbf{F_{i}}/M_{tot}$. Finally, for a Hamiltonian with positions referenced to the centre of mass of a subset of particles, Hamilton’s equations require that only those particles feel the reverse acceleration $-\mathbf{F_{i}}/M_{sub}$, where $M_{sub}$ is the mass of the subset of particles. In the case where each particle’s position is referenced to the centre of mass of all interior bodies, this corresponds to the familiar Jacobi coordinates.