REBOUNDx: A Library for Adding Conservative and Dissipative Forces to Otherwise Symplectic N-body Integrations
Daniel Tamayo, Hanno Rein, Pengshuai Shi, David M. Hernandez
MMNRAS , 000–000 (0000) Preprint 10 October 2019 Compiled using MNRAS L A TEX style file v3.0
REBOUNDx: A Library for Adding Conservative and DissipativeForces To Otherwise Symplectic N-body Integrations
Daniel Tamayo (cid:63) , Hanno Rein , , Pengshuai Shi , and David M. Hernandez , , Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544 Department of Physical and Environmental Sciences, University of Toronto at Scarborough, Toronto, Ontario M1C 1A4, Canada Department of Astronomy and Astrophysics, University of Toronto, Toronto, Ontario, M5S 3H4, Canada Harvard–Smithsonian Center for Astrophysics, 60 Garden St., MS 51, Cambridge, MA 02138, USA Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, 77 Massachusetts Ave., Cambridge, MA 02139, USA RIKEN Center for Computational Science, 7-1-26 Minatojima-minami-machi, Chuo-ku, Kobe, 650-0047 Hyogo, Japan
Draft: 10 October 2019
ABSTRACT
Symplectic methods, in particular the Wisdom-Holman map, have revolutionized our abilityto model the long-term, conservative dynamics of planetary systems. However, many astro-physically important e ff ects are dissipative. The consequences of incorporating such forcesinto otherwise symplectic schemes is not always clear. We show that moving to a generalframework of non-commutative operators (dissipative or not) clarifies many of these ques-tions, and that several important properties of symplectic schemes carry over to the generalcase. In particular, we show that explicit splitting schemes generically exploit symmetries inthe applied external forces which often strongly suppress integration errors. Furthermore, wedemonstrate that so-called ‘symplectic correctors’ (which reduce energy errors by orders ofmagnitude at fixed computational cost) apply equally well to weakly dissipative systems andcan thus be more generally thought of as ‘weak splitting correctors.’ Finally, we show thatpreviously advocated approaches of incorporating additional forces into symplectic methodswork well for dissipative forces, but give qualitatively wrong answers for conservative butvelocity-dependent forces like post-Newtonian corrections. We release REBOUNDx , an open-source C library for incorporating additional e ff ects into REBOUND
N-body integrations, to-gether with a convenient
PYTHON wrapper. All e ff ects are machine-independent and we pro-vide a binary format that interfaces with the SimulationArchive class in
REBOUND to enablethe sharing and reproducibility of results. Users can add e ff ects from a list of pre-implementedastrophysical forces, or contribute new ones. Key words: methods: numerical — gravitation — planets and satellites: dynamical evolutionand stability
The long-term dynamical evolution of planetary systems remainsa rich challenge for analytical investigations. As a result, many ofthe advances in the last several decades have been driven by nu-merical studies thanks to faster computers and the development ofimproved algorithms for accurate and e ffi cient numerical solutions,e.g., the chaotic evolution of Pluto (Sussman & Wisdom 1988),chaos in the inner solar system (Laskar 1989), and the marginalinstability of Mercury over the age of the solar system Laskar &Gastineau (2009). In many cases, such chaos only manifests itselfafter millions or even billions of orbits. This imposes strong de- (cid:63) NHFP Sagan Fellow: [email protected] mands on the long-term conservation properties of numerical inte-gration methods that can be used for such studies.These concerns led to the development powerful symplec-tic 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 (Wisdom1982; Kinoshita et al. 1990; Wisdom & Holman 1991) enabledthe first direct N-body integrations of our solar system over Gyrtimescales. By respecting the Hamiltonian structure of the prob-lem, symplectic methods render the task of integration equivalent toa canonical transformation (e.g. Sanz-Serna 1992), which enforcesthe conservation of several invariants of the phase space flow. Thisavoids the secular energy error growth of many of the more generalmethods that are not restricted to Hamiltonian systems, which canlead to unphysical collisions in long-term simulations.Such symplectic integrators split the problem into multiple c (cid:13) a r X i v : . [ a s t r o - ph . E P ] O c t Daniel Tamayo, Hanno Rein, Pengshuai Shi, David Hernandez pieces that are each evolved in sequence. Apart from the robust ge-ometrical and numerical properties this provides, the modularity ofsplitting methods makes it possible to extend them for hierarchicalproblems (Fujii et al. 2007; Portegies Zwart & McMillan 2018) andto incorporate di ff erent physics in a single integration (Pelupessy &Portegies Zwart 2012; Portegies Zwart 2018). This has for examplebeen vigorously pursued in the Astrophysics Multipurpose Soft-ware Environment ( AMUSE ) package (Portegies Zwart & McMillan2018).Some e ff ects, such as bodies’ higher gravitational moments,can be expressed as simple position-dependent potentials and canbe trivially incorporated into symplectic schemes. However, manycrucial perturbations such as post-Newtonian corrections, tides, orradiation forces are velocity-dependent and even dissipative. Mal-hotra (1994) and Cordeiro et al. (1996) have proposed general-izations of symplectic schemes for incorporating weak velocity-dependent forces, but their consequences in long-term integrationsare not well understood. While some authors perform convergencetests to check the validity of their adopted stepsize using variantsof the above methods (e.g. Zhang & Hamilton 2007), the numericalrobustness of studies where this is not done is unclear.A generalization of symplectic integration theory would there-fore be valuable to understand the source of numerical error in theabove schemes, and to provide simple estimates for its magnitudeto guide the selection of timestep for a given problem.A promising direction has been o ff ered by Galley (2013) andGalley et al. (2014), who have developed a generalized theory forclassical mechanics that can incorporate dissipative processes inthe action governing the dynamics. One can use this formalism toconstruct variational integrators that generalize symplectic meth-ods to accurately track the changing momentum and energy of thesystem on long timescales (Tsang et al. 2015). Potential disadvan-tages of this elegant framework are that the equations of motion arenecessarily implicit (Tsang et al. 2015), and that adding new forcesbecomes somewhat more complicated than with traditional meth-ods, since one needs to derive a ‘nonconservative potential’ usingthe formalism in Galley et al. (2014).In this paper, we explore how the established theory of sym-plectic integration can instead be generalized in a framework ofnon-commutative operators, which can be applied to both conser-vative and dissipative systems. Such generalized ‘splitting meth-ods’ have been suggested in the field of fluid mechanics as farback as Strang (1968), and more recently for celestial mechanicsby Mikkola (1998) . We show that several strong results tradition-ally discussed in the symplectic integration literature carry over tomore general and possibly dissipative splitting methods.Readers interested in quickly using REBOUNDx to incorporateastrophysical e ff ects into N-body integrations could begin at Sec. 6for an overview of the implementation. For a more gradual intro-duction, we introduce in Sec. 2 the notation, and review symplecticintegration and the Wisdom-Holman map in an operator-centredformalism. We then show in Sec. 3 how this framework can ex-plain the energy error behaviour of such methods under weak per-turbations, both conservative and dissipative. In Sec. 4 we reviewhow the WH map can be extended to higher order without addi-tional force evaluations, and show that such ‘symplectic correc-tors’ can also be applied to dissipative systems and are thus moregeneral than typically considered. In Sec. 5, we consider the gen- See also Hernandez & Bertschinger (2018) for error analysis of one-stepmethods from the perspective of their di ff erential equations. eral case where analytic solutions to individual steps in a splitscheme are not known, and how truncation errors from numeri-cally integrating across these individual steps interact with the over-all splitting scheme errors. In particular, we show that for conser-vative, velocity-dependent forces like post-Newtonian corrections,the methods of Cordeiro et al. (1996) and Malhotra (1994) givequalitatively wrong answers, and we show how to correct them.We summarize our results and conclude in Sec. 7. We begin by reviewing the basic theory behind symplectic integra-tion. 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 discus-sion in the main text. For a more careful introduction, and to correctsome sign errors in the literature, see Appendix A.
We consider a system of N particles with positions r i and velocities v i . We define a di ff erential operator ˆ P , which acts on the currentstate of the system z = ( r , v ) to yield ˙ z , i.e., the particular set ofdi ff erential equations we are trying to solve ,ˆ P z : ˙r i = v i ˙v i = a i ( r ) , (1)where a i is the i th particle’s acceleration vector, which depends onthe positions r of all the bodies.Next, we introduce a further level of abstraction by defining asolution to Eq. 1 through a corresponding operator P ( h ). This idealintegrator exactly advances the state z by a timestep h according tothe set of di ff erential equations ˙ z = ˆ P z .In general, there is no closed form solution P (e.g., the three-body problem), but one could split the di ff erential equations intotwo pieces that could each be solved trivially in isolation,ˆ A z : ˙r i = v i ˙v i = B z : ˙r i = ˙v i = a i ( r ) , (2)In particular, the solution A ( h ) keeps the velocities constant andupdates the positions by h v i , while B ( h ) keeps the positions con-stant and similarly updates the velocities using constant acceler-ations. This splitting works for any di ff erential equations of theform in Eq. 1, and is thus widely applicable. In a Hamiltonianframework, in cases where the accelerations can be derived from aposition-dependent potential, it corresponds to splitting the Hamil-tonian into a kinetic piece containing all the momenta, and a po-tential piece containing all the positions. We therefore refer to thisscheme as a kinetic-potential splitting, which was the focus of mostearly symplectic integrators (e.g., Ruth 1983; Neri 1987; Forest &Ruth 1990; Yoshida 1990). As we will see below, better splittingsfor specialized problems are possible.The idea behind splitting methods is to alternate evolution un-der two or more operators that can be calculated exactly and e ffi -ciently. For example, we can define a first-order splitting scheme as z ( t + h ) ≈ SAB ( h ) z ( t ) ≡ A ( h ) ◦ B ( h ) z ( t ) , (3) Throughout this paper we restrict ourselves to time-independent sets ofdi ff erential equations. MNRAS , 000–000 (0000) EBOUNDx: A Library for Additional Forces in N-body simulations where B advances the state z ( t ) by h according to its half of thedi ff erential equations, and then A advances the result by h withthe remaining half. Unfortunately, doing two halves of the problemin sequence is not the same as doing the full problem all at once.Nevertheless, if we repeatedly apply this map to get from some z ( t ) to z ( t ), at least in the limit of a vanishing timestep, the splitoperations become densely interleaved and this approximation willapproach the real solution.For finite timesteps, we can also see simply that how well asplit scheme will correspond to the real solution should depend onthe degree to which the operators A ( h ) and B ( h ) commute withone another. We could always split each of the steps in Eq. 3 in halfand write it as A ( h / ◦ A ( h / ◦ B ( h / ◦ B ( h / A ( h / ◦ B ( h / ◦ A ( h / ◦ B ( h / ff ectivelybe applying both at the same time, yielding the true solution. Sowe expect the errors in splitting schemes to vanish in the limits thateither h or the commutator of A and B go to zero.The formula that quantifies the above two statements isthe Baker-Campbell-Hausdor ff (BCH) identity, which is useful inmany branches of physics involving non-commuting operators likequantum mechanics. To state it, we have to add a final layer of ab-straction and generalize beyond just how the phase space variables z change under a particular subset of the di ff erential equations (e.g.,Eq. 2), to how an arbitrary function g ( z ) changes with time. In par-ticular, for every subset of di ff erential equations, e.g., ˙z = ˆ A z , wecan define the corresponding Lie derivative L A that, when acting on g ( z ), yields dg / dt under that particular subset of di ff erential equa-tions. In the simple case of g ( z ) = z , we get back L A z = ˙ z = ˆ A z .The BCH formula says that when we apply Eq. 3, we are notsolving the original set of di ff erential equations L P z = L A z + L B z (Eqs. 1 and 2), but rather a nearby set of di ff erential equations˙ z = L SAB z = L P z + L err SAB z , (4)where the additional error term is a power series in the timestepconsisting of nested commutators [ L A , L B ] = L A L B − L A L B (e.g.Saha & Tremaine 1992), L err SAB = − h L A , L B ] + h
12 [ L A − L B , [ L A , L B ]] − O (cid:16) h (cid:17) . (5)In other words, if we integrate the set of di ff erential equations inEq. 4 numerically to high accuracy, we match the evolution gen-erated by the simple splitting scheme in Eq. 3. We will call thesedi ff erential equations that we are actually solving the ‘modified setof di ff erential equations’ for the splitting scheme. Understandingthe structure of this nearby problem gives insight into the numeri-cal errors introduced by the method .As expected from the qualitative arguments above, the modi-fied di ff erential equations approach the true ones as h approacheszero. We also see that if (and only if) [ L A , L B ] =
0, then there is noerror and A ( h ) and B ( h ) commute, i.e., A ( h ) ◦ B ( h ) = B ( h ) ◦ A ( h ).By choosing appropriate steps, one can compose A and B intohigher order schemes (e.g., Yoshida 1990). A widely used scheme It is important to point out that the BCH formula represents a formal,asymptotic series, which in particular is not guaranteed to converge every-where in phase space (Wisdom 2018). In practice, as long as the adoptedtimestep is (cid:46)
10% of the fastest timescale in the problem, this is not typi-cally a concern. is the time-symmetric ‘leapfrog’ method
SABA ( h ) ≡ A ( h / ◦ B ( h ) ◦ A ( h / . (6)Repeated application of the BCH formula shows that this particularcomposition cancels out error terms linear in h , yielding a second-order integrator with L err SABA = h
24 [2 L B + L A , [ L B , L A ]] + O (cid:16) h (cid:17) . (7)Because the scheme is time-symmetric, errors only appear at evenpowers of h . The above analysis does not depend on us solving a Hamiltoniansystem, and is general to both conservative and dissipative systems.However, the insight of symplectic integration is that if each of thesplit operators, A ( h ) and B ( h ), are exact solutions to a Hamiltoniansystem, then each of them conserve various Hamiltonian invariants.This implies that any composition of them, e.g., A ( h ) ◦ B ( h ), mustalso conserve those quantities. In particular, such schemes will con-serve the Poincar´e invariants, the linear and angular momentum,and the infinite di ff erentiability class order of the governing di ff er-ential equations (Hernandez 2019a), all of which are important forthe accuracy of solutions (Hernandez 2019b). Having these strongconservation properties built into symplectic schemes gives themexcellent long-term behaviour.The main subtlety, however, is that each split operator, e.g. A ( h ), is conserving its own Hamiltonian, e.g. H A . It is thereforenot clear what Hamiltonian an overall composition scheme like theone in Eq. 3 is conserving. Then again, a Hamiltonian H A yields theassociated set of di ff erential equations ˙ z = L A z through Hamilton’sequations, so it should not be surprising that this answer can alsobe derived from the BCH formula.One can show that one gets analogous formulae to Eqs. 5 and 7(e.g., Yoshida 1993), H err SAB = h { H A , H B } + h { H A − H B , { H A , H B }} + O ( h ) , (8)and H err SABA = h { H B + H A , { H B , H A }} + O ( h ) , (9)where curly brackets denote Poisson brackets. Thus, each of theerror terms in the modified di ff erential equations (e.g., Eq. 5) can bederived from a corresponding error Hamiltonian (e.g., Eq. 8) (e.g.,Saha & Tremaine 1992). This is another way of seeing that theintegration errors respect the symplectic geometry of the problem. While the energy errors in kinetic-potential splittings do not driftsecularly with time, they remain large throughout the integration.Physically, this is because one is making drastic deviations fromthe true trajectory each timestep, i.e., force-free motion along A ( h )alternating with order-unity kicks in the particle velocities under B ( h ).In the case of weakly perturbed systems, like planets movingon nearly Keplerian orbits, it is much more powerful to instead splitinto an integrable dominant operator and the weak perturbation.Such a scheme, today known as the Wisdom-Holman map, proveda major development for planetary integrations, enabling the first MNRAS000
24 [2 L B + L A , [ L B , L A ]] + O (cid:16) h (cid:17) . (7)Because the scheme is time-symmetric, errors only appear at evenpowers of h . The above analysis does not depend on us solving a Hamiltoniansystem, and is general to both conservative and dissipative systems.However, the insight of symplectic integration is that if each of thesplit operators, A ( h ) and B ( h ), are exact solutions to a Hamiltoniansystem, then each of them conserve various Hamiltonian invariants.This implies that any composition of them, e.g., A ( h ) ◦ B ( h ), mustalso conserve those quantities. In particular, such schemes will con-serve the Poincar´e invariants, the linear and angular momentum,and the infinite di ff erentiability class order of the governing di ff er-ential equations (Hernandez 2019a), all of which are important forthe accuracy of solutions (Hernandez 2019b). Having these strongconservation properties built into symplectic schemes gives themexcellent long-term behaviour.The main subtlety, however, is that each split operator, e.g. A ( h ), is conserving its own Hamiltonian, e.g. H A . It is thereforenot clear what Hamiltonian an overall composition scheme like theone in Eq. 3 is conserving. Then again, a Hamiltonian H A yields theassociated set of di ff erential equations ˙ z = L A z through Hamilton’sequations, so it should not be surprising that this answer can alsobe derived from the BCH formula.One can show that one gets analogous formulae to Eqs. 5 and 7(e.g., Yoshida 1993), H err SAB = h { H A , H B } + h { H A − H B , { H A , H B }} + O ( h ) , (8)and H err SABA = h { H B + H A , { H B , H A }} + O ( h ) , (9)where curly brackets denote Poisson brackets. Thus, each of theerror terms in the modified di ff erential equations (e.g., Eq. 5) can bederived from a corresponding error Hamiltonian (e.g., Eq. 8) (e.g.,Saha & Tremaine 1992). This is another way of seeing that theintegration errors respect the symplectic geometry of the problem. While the energy errors in kinetic-potential splittings do not driftsecularly with time, they remain large throughout the integration.Physically, this is because one is making drastic deviations fromthe true trajectory each timestep, i.e., force-free motion along A ( h )alternating with order-unity kicks in the particle velocities under B ( h ).In the case of weakly perturbed systems, like planets movingon nearly Keplerian orbits, it is much more powerful to instead splitinto an integrable dominant operator and the weak perturbation.Such a scheme, today known as the Wisdom-Holman map, proveda major development for planetary integrations, enabling the first MNRAS000 , 000–000 (0000)
Daniel Tamayo, Hanno Rein, Pengshuai Shi, David Hernandez integrations of the solar system over Gyr timescales (Wisdom &Holman 1991).While the development of the Wisdom-Holman map was notdriven by this operator framework (see Wisdom 2018 for a histori-cal perspective), it is instructive for our purposes to analyze it in thislight. Instead of splitting the di ff erential equations into componentsof comparable magnitudes like in Eq. 2, the Wisdom-Holman mapsplits them as (Wisdom & Holman 1991; Kinoshita et al. 1990), L K z : ˙r i = v i , ˙v i = a Kepi ( r ) (cid:15) p L I z : ˙r i = , ˙v i = (cid:15) p a inti ( r ) (10)where the a Kep i are the two-body Keplerian accelerations, and the a int i are the remaining interaction accelerations . We have also in-troduced a factor of (cid:15) p of order the characteristic planet-star massratio to keep track of the fact that the interplanetary accelera-tions (and therefore the changes in the interaction steps) are muchsmaller than the dominant Keplerian accelerations due to the cen-tral star.Like in the kinetic-potential split, the corresponding I ( h ) istrivial, since it keeps the positions (and thus the accelerations) con-stant across the timestep. The exact solution for the velocities isthen just v i ( t + h ) = v i ( t ) + (cid:15) p h a inti . (11)By contrast, the Kepler step K ( h ) updates both the positions and ve-locities of each body along their respective unperturbed two-bodyorbits. Interestingly, more accurate and e ffi cient schemes for solv-ing the Kepler problem have continually been found for over threehundred years (e.g., Newton 1687; Machin 1737; Rambaut 1890;Plummer 1896; Brown 1931; Danby 1992; Mikkola & Innanen1999; Rein & Tamayo 2015a; Wisdom & Hernandez 2015). As an example, consider a case of two Earth-mass planets withsemi-major axes of 1 and 2 AU around a solar-mass star ( (cid:15) p = × − ) with orbital eccentricities of 0.01. We specialize to secondorder schemes (Eq. 6), and compare the kinetic-potential splitting(Eq. 2) to the Wisdom-Holman splitting (Eq. 10). Figure 1 plotsboth methods’ fractional energy errors over ten thousand inner-planet orbits, using a timestep of 1% the innermost planet’s orbitalperiod.The Wisdom-Holman map conserves the energy many ordersof magnitude better than the kinetic-potential splitting, and we canmake simple order of magnitude estimates to understand why. Thisis clearest in a Hamiltonian framework for conservative systems(e.g. Wisdom & Holman 1991; Saha & Tremaine 1992), but weinstead show how it arises in the above operator-centred picturethat will later generalize to dissipative systems. While this is not a These capture both the direct gravitational interactions between planetsand their indirect e ff ects on one another as they each pull on the centralstar. There are many ways to do this split between Keplerian and interac-tion accelerations (Hernandez & Dehnen 2017; Rein & Tamayo 2019). Thesplitting above into two operators corresponds to using Jacobi coordinates(Wisdom & Holman 1991). This assumes that the distances from the planets to the central star arecomparable to their interplanetary separations. In general, (cid:15) p would capturethe ratio of the perturbation accelerations to those of the dominant Keplerianmotion. Time (Inner Planet Orbits) R e l a t i v e H a m il t o n i a n E rr o r ( h / T ) e p ( h / T ) Kinetic-PotentialWisdom-Holman
Figure 1.
Integrations of two nearly circular ( e = h of 1% the inner planet’s orbitalperiod T . We compare a kinetic-potential splitting to a Wisdom-Holmansplit. Black lines show the scheme error estimates in Sec. 2.4, except thekinetic-potential splitting error has an additional factor of e (Sec. 3.2). rigorous derivation, it is useful to be able to quickly estimate ex-pected errors and their associated scalings. While we do not plot allpossible dependencies, we have verified that the various scalingspredicted below are correct.A useful feature of the abstract Lie derivative formalism ofSec. 2 is that they describe not only how the phase space variables z evolve, but also how functions like the energy E ( z ), change withtime. In particular, just like ˙ z err = L err S z yields the additional termsin the modified di ff erential equations for the phase space variablesusing a particular split scheme S (Sec. 2), ˙ E err = L err S E ( z ) yieldsthe additional terms in the modified di ff erential equations for theenergy evolution.To approximately estimate the error ∆ E err they cause, we canmultiply the additional terms ˙ E err by the timescale over which theseerrors coherently add up. Toward this end, on the left of Fig. 1, thelogarithmic time axis makes it possible to see that the errors ac-cumulate coherently over several timesteps ( h = . ff erential equations are composed of the split opera-tors that propagate the real dynamics (Eq. 7), the errors will inheritthe timescales in the original problem. In particular, the timescaleover which errors will accumulate coherently will be the shortesttimescale in the problem. In this case of near-circular, non-resonantorbits, T would be the synodic period on which planets kick eachothers’ orbits every conjunction. In our typical example where theorbits are not too close, this is approximately the orbital period ofthe innermost planet; for a very eccentric orbit, T would be themuch shorter timescale for pericentre passage (e.g., Rauch & Hol-man 1999; Rein & Tamayo 2015b), etc. We thus have (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∆ E err E (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∼ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) T ˙ E err E (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) T L err S EE (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (12)where, for the second-order schemes used in Fig. 1, L err S is given byEq. 7.Since L A and L B each yield time derivatives according totheir own subset of the di ff erential equations, we can make a roughestimate in Eq. 12 by replacing each instance of L A and L B in L err S by their respective inverse characteristic timescale. In this case weare referring specifically to the timescale on which each operator inisolation would change the particle states by order unity, i.e., z / ˙z. MNRAS , 000–000 (0000)
EBOUNDx: A Library for Additional Forces in N-body simulations For Kepler splittings, the Kepler Lie derivative L K inducesorder-unity changes on the orbital timescale T . Perturbation op-erators like the interaction step in the WH map induce order-unity changes on timescales that are longer by roughly a factorof the perturbation strength (cid:15) (e.g., for the interaction operator,v / ˙v I = (v / ˙v K )(˙v K / ˙v I ) = (cid:15) p T , see Eq. 10).If we further assume that the two operators L A and L B do notcommute, i.e., [ L A , L B ] ∼ L A L B , then plugging Eq. 7 into Eq. 12yields a fractional energy error estimate of ∼ (cid:15) p ( h / T ) . This givesthe right scaling for a second-order method, as expected. Indeed,the same argument applied to any Kepler splitting says that if theperturbation accelerations are a factor of (cid:15) smaller than the Kep-lerian ones, then the second-order splitting scheme given by Eq. 6will yield energy errors ∼ (cid:15) ( h / T ) .By contrast, in the kinetic-potential splitting (Eq. 2), both L A and L B induce order-unity changes on the orbital timescale T , so L A ∼ L B ∼ / T . This yields a correspondingly larger error ( h / T ) .More physically, the WH method is benefiting from the factthat the dominant motion is known and solved exactly. In the limitof (cid:15) p →
0, the Wisdom-Holman map would yield an exact trajec-tory, while the kinetic-potential splitting would still be splitting theKepler problem and the errors would remain large and unchanged.These simple estimates (black arrow in Fig. 1) match well tointegrations with the Wisdom-Holman map, while our estimate forthe kinetic-potential splitting is depressed by approximately a fac-tor of e from our simple estimate. This has to do with our simpli-fication of the commutation relations, and has interesting broaderimplications that we discuss in Sec. 3.2.In summary, while the Kepler step is more expensive to com-pute than the steps in a kinetic-potential splitting, the factor of (cid:15) p gain makes reaching a given level of accuracy significantly fasterwith the Wisdom-Holman map. Perhaps more importantly, in caseswhere inter-planetary forces are very weak, the Wisdom-Holmanmethod keeps the additional terms in the modified di ff erential equa-tions smaller than the interplanetary terms. By contrast, in thekinetic-potential splitting, these additional error terms can easily belarger than the inter-planetary forces, leading to spurious results.These scalings suggest the cases where splitting methods willbe most useful, i.e. when perturbations are weak (cid:15) (cid:28) T is not too short. For exam-ple, for a planet experiencing dynamical tides on an extremely ec-centric orbit, the strength of the perturbation changes by orders ofmagnitude over an orbital period. In this case, the shortest timescalewould be the characteristic timescale of pericentre passage, whichcan be much shorter than an orbital period. Because simple sym-plectic schemes require a fixed timestep, the small h required wherethe forces are strongest acts as a bottleneck for the remaining or-bit that could otherwise be easily integrated. In such cases it istypically more e ffi cient to use a more general scheme with adap-tive timesteps like IAS15 , or in some restricted cases it is possi-ble to symplectically adapt timesteps in inverse proportion to thestrength of an external potential (Preto & Tremaine 1999; Mikkola& Tanikawa 1999; Petit et al. 2019). More precisely 2 π/ T . Considering the simple 1-D case of a circular orbitwith semimajor axis a and orbital frequency ω , ˙v = ω a , so L K ∼ v / ˙v = /ω . We are ignoring these factors of 2 π , which largely cancel with thecoe ffi cients in the BCH formula, but they can be important for higher-orderschemes that accumulate many such factors. We now move beyond point-particle gravity to consider adding ad-ditional e ff ects in N-body integrations. We will see that the aboveframework, by not specializing to Hamiltonian systems, can beused to straightforwardly understand the numerical behavior underboth conservative and dissipative forces. We begin by considering conservative forces that can be derivedfrom position-dependent potentials. Hamiltonian perturbations thatcannot be written strictly in terms of the particle positions are alsoimportant (e.g., post-Newtonian corrections), but we defer their dis-cussion to Sec. 5.We note that while any position-dependent potential can obvi-ously be trivially incorporated into a kinetic-potential split, it canalso be directly inserted into the interaction step in the Wisdom-Holman map as an additional position-dependent acceleration inEq. 10. Then Eq. 11 remains the exact solution with the accelera-tion given by the sum of all the position-dependent accelerations inthe problem.As an example, we consider the same two-planet case asabove, but now with additional perturbations from an oblate pri-mary, with the planets orbiting in the primary’s equatorial plane.The first corrections to point source gravity in a multipole expan-sion of the primary’s potential is the quadrupole term, with a po-tential in the equatorial plane of V J ( r ) = − J (cid:32) Rr (cid:33) GMr ≡ (cid:15) J GMa (cid:32) a r (cid:33) . (13)The innermost body will be most a ff ected, so we take a to be itsoriginal, reference semimajor axis, G is the gravitational constant, M and R are the primary’s mass and radius, and J is the standarddimensionless coe ffi cient for the quadrupole field. We have alsointroduced a parameter (cid:15) J = J (cid:16) Ra (cid:17) that captures the smallnessof the e ff ect relative to the dominant Keplerian potential ≈ GM / a ,since a / r is approximately unity for nearly circular orbits. We take (cid:15) J = − .We again compare second-order schemes (Eq. 6) usingkinetic-potential and Wisdom-Holman Kepler splittings in Fig. 2.Comparing with Fig. 1, we see that the kinetic-potential errors havenot changed. This is because the planetary accelerations in L B for the kinetic-potential splitting (Eq. 2) are dominated by thosedue to the star. Adding a small acceleration due to the oblatenessperturbation therefore does not noticeably change the dominant er-ror terms in Eq. 7. By contrast, for the Kepler splitting, our cho-sen oblateness perturbations are much larger than the interplane-tary ones ( (cid:15) J (cid:29) (cid:15) p ), so they change the error behavior. The sameargument following Eq. 12 suggests an energy error ∼ (cid:15) J ( h / T ) .However, both our error estimates for the kinetic-potential andKepler splitting fall short by a factor of the eccentricity (both blackarrows in Fig. 2 include this factor). This can be a significant sup-pression in the errors and is due to the geometrical properties ofsplitting methods, as we explore in the next section. Because for the remainder of the paper we focus on Kepler split-tings, 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
MNRAS000
MNRAS000 , 000–000 (0000)
Daniel Tamayo, Hanno Rein, Pengshuai Shi, David Hernandez Time (Inner Planet Orbits) R e l a t i v e H a m il t o n i a n E rr o r ( h / T ) e J ( h / T ) e Kinetic-Potential SplitKepler Split Eccentricity R e l a t i v e H a m il t o n i a n E rr o r ee Kinetic-Potential SplitKepler Split
Figure 2.
Top panel is the same as Fig. 1, except we now additionally in-clude a stellar oblateness perturbation of (cid:15) J = − . Estimated schemeerrors in the main text are labeled and denoted by black arrows. Bottompanel shows how the fractional energy error depends on the orbital eccen-tricities of the planets. It is linear in the eccentricity (solid lines) down to aneccentricity e ∼ (cid:15) J , as explained in the main text. all orbital elements, except for the mean longitude λ , which ad-vances linearly. This implies that when alternating Kepler stepswith perturbation steps, only perturbation components that changethe particles’ semi-major axes (which determine the mean motions d λ/ dt ) do not commute.For example, for a perturbation that only changed the ec-centricity vector’s magnitude and direction, it would not matterwhether the Kepler or perturbation step acted first, since the Keplerstep leaves those quantities unchanged. The ordering only mattersfor perturbations that change the semimajor axis; if the perturbationacts first, the mean motion will change, and when the Kepler stepoperates, the particle will end up at a di ff erent mean longitude thanif the Kepler step had acted first with the original value of the meanmotion.We can then see that for any purely radial force like the oblate-ness perturbation, the change in orbital energy (and thus semi-major axis) − F · v vanishes for circular motion. If we slowly in-crease the initial eccentricity, the perturbation will add and extractorbital energy as each body moves toward and away from pericen-tre, to leading order in proportion to e . So our mistake in estimatingthe energy error in Sec. 2.4 was in assuming that the commutatorsintroduced factors of order unity, i.e., that [ L A , L B ] ∼ L A L B . Ourphysical argument says that for radially perturbed, nearly circularorbits, the operators in fact approximately commute. The commu-tator [ L A , L B ] oscillates with an amplitude ∼ e L A L B and vanishesat pericentre and apocentre where F · v = This is physically what is causing the oscillationsvisible on the left of the top panel in Fig. 2, where the logarithmictime axis renders the orbital timescale visible. The planets wereboth started at pericentre, and one can see that the energy errorsremain constant at pericentre (every time unit) and apocentre (ev-ery half-time unit). We note that the flat errors for e < − in thebottom plot of Fig. 2 are an artefact of our numerical setup .The above example highlights the power of splitting methods’geometrical properties. Because all the error terms in the modifieddi ff erential equations (e.g., Eq. 7) consist of nested commutatorsof L A and L B , any symmetries (or near-symmetries) are inher-ited in the integrator’s error properties to all orders. This meansany higher-order splitting scheme that tries to correct for higherand higher order terms in the BCH expansion (e.g., Wisdom et al.1996; Laskar & Robutel 2001) will always have such a geometricalsuppression in its leading error term.This is not typically emphasized in the literature on symplecticN-body integration, presumably because the Kepler and interactionsteps do not have any such symmetries (there are always conjunc-tions between planets that change the semi-major axes at leadingorder). However, many astrophysically important e ff ects are highlysymmetric potentials, e.g. multipole gravitational potentials, radia-tion pressure, simple general relativistic corrections, etc. This fact,combined with the the highly degenerate Kepler problem, oftenlead to strongly suppressed errors for the nearly circular or nearlycoplanar orbits we often want to model.In summary, splitting methods are often powerful not only fortheir long-term conservation of important quantities, but also be-cause they can yield orders of magnitude higher accuracy for a fixedtimestep. In the above case, a specialized second-order scheme splitinto separate Kepler and perturbation steps reduced the errors by afactor of (cid:15) J e = − over short timescales compared to what onewould obtain with a generic second-order Runge-Kutta scheme,even more over long timescales. While dissipative systems are qualitatively di ff erent from the con-servative cases above, we now show that we can similarly under-stand their error behaviour using the above framework. As shouldbe clear from the operator-centred development throughout the pa-per, we could take any di ff erential equations, split them in any waywe please, and if we could find time-evolution operators for thosesubproblems, then we could compose them and find their error be-havior through the BCH formula. In particular, the error estimates Actually it would not completely vanish in this problem due to the muchweaker perturbations from the additional planet that would at these pointsbecome the dominant e ff ect. The physical argument is correct for all geometric eccentricities (i.e., thevalue measured by looking at the shape of the physical orbit in space, whichincludes the e ff ects of the perturbation) but the typical subtlety arises whendealing with osculating eccentricities (i.e. the transformation from positionsand velocities to unperturbed 2-body Kepler orbits, which ignores the e ff ectof the perturbation) smaller than the size of the perturbation (cid:15) . Because forconvenience we initialize our orbits using osculating elements, decreasingthe initial osculating eccentricities below (cid:15) J does not make them any morecircular in geometric space, because the oblateness perturbations kick themaway from Keplerian motion at order (cid:15) J . So in our integrations, the ge-ometric eccentricity that goes into the argument above reaches a floor at ∼ (cid:15) J , with a corresponding error (cid:15) J ( h / T ) e = (cid:15) J ( h / T ) , which matchesthe flat regime on the left of the bottom panel of Fig. 2.MNRAS , 000–000 (0000) EBOUNDx: A Library for Additional Forces in N-body simulations we developed in Sec. 2.4 and applied in Sec. 3.1 come from themodified di ff erential equations, without reference to conservativeor dissipative forces.On the other hand, one important distinction for symplecticsystems is that the fact that the composition of two Hamiltonianoperators must be Hamiltonian implies that all the BCH error termsare also Hamiltonian (Sec. 2.2). By contrast, when the perturbationis dissipative, the presence of L B operators in the error terms ofthe modified di ff erential equations (e.g., Eq. 7) shows that there isadditional damping (or injection of energy) introduced by the split-ting scheme. This generically leads to secular drifts, as we show inthe following example.Consider a single planet in orbit around its primary with aninitial orbital eccentricity of 0 . v , F = − m v τ a , (14)where m is the planetary mass. This parametrized prescriptionorbit-averages to yield inward migration with the semi-major axisdecaying exponentially on an e-folding timescale τ a (Papaloizou &Larwood 2000). We set τ a = τ a . In this time, the planet’s semi-major axismoves inward by a factor of exp(3) ≈
20, and the orbital perioddecreases by a factor of exp(9 / ≈ K ( h / ◦ D ( h ) ◦ K ( h / , (15)where K evolves the planet on a Kepler orbit, and D damps themotion according to Eq. 14. We adopt a timestep of 10 − times theplanet’s initial orbital period, and choose to approximate D ( h ) byintegrating across the timestep using a fourth-order Runge-Kuttascheme (see Sec. 5).With dissipation, we no longer have a conserved quantity totrack. Instead, as a proxy for the exact propagator P , we also inte-grate the system with IAS15 (Rein & Spiegel 2015), a high-orderadaptive-timestep method whose accuracy reaches machine preci-sion. We plot in orange in Fig. 3 the relative error between the en-ergy calculated in our Kepler splitting integration and the one with
IAS15 .One can see a clear secular drift. While the errors in the sym-plectic integrations in Fig. 2 oscillate and average out to yield flattime evolution, dissipative errors systematically overdamp or un-derdamp.This scaling can be understood in the same way we analyzedthe conservative case. The second-order splitting scheme error esti-mate is (Sec. 2.4) (cid:15) D ( h / T ) , with T the planet’s orbital period, and (cid:15) D = T /τ a . However, like in the case of J perturbations discussedabove, there is additionally a geometric suppression of the error bya factor of the eccentricity . We show this estimate in Fig. 3 with ablack dashed line. To 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 changethe particles’ semi-major axes don’t commute with the Kepler step. Herethe force is always pointed exactly opposite the particle’s velocity, so thetwo steps definitely do not commute. But now we additionally have to askwhether the order of operations matters for the quantity whose error we aremeasuring, i.e., the energy. For a circular orbit, it would not. Whether ornot the Kepler step moves the planet along the circle before or after the stepdoes not a ff ect the energy loss, since the problem is azimuthally symmetric.So again the energy errors are suppressed by a factor of e , though in this Time (Years) R e l a t i v e E n e r g y E rr o r D ( h / T ) e Kepler Split
Figure 3.
Second-order, Kepler-splitting integration of a single planet ini-tially on a 1-year, 0.1 eccentricity orbit, subject to a damping force thatcauses the orbital period to decay exponentially. The relative energy errorfollows the same estimate (dashed black line) as the conservative case inFig. 2, accounting for the fact that the orbital period T is changing. The secular rise then simply reflects the fact that the orbital pe-riod is changing exponentially as the planet migrates inward, whichyields a straight line on the linear-log scale.Remarkably, this suggests that exponential outward migrationwould yield errors that exponentially decrease with time. We findthat this is indeed the case, though this case is more subtle. In thecase plotted above where the errors committed each timestep in-crease with time, the most recent errors always dominate the er-ror budget. In the case where the instantaneous errors are gettingsmaller, one might expect the total error to remain at the level in-curred at the beginning of the integration, where errors are largest.But these errors are oscillatory, and as long as the errors are chang-ing adiabatically (i.e., as long as there are many orbits per migra-tion timescale), the oscillations can march toward smaller ampli-tude. Of course for exponential outward migration, it does not takemany migration timescales until the orbital period becomes com-parable to the migration timescale, at which point the applied forcein Eq. 14 (unphysically) becomes comparable to the central gravi-tational force, and we see numerical errors rise again.In summary, the error behaviour in weakly dissipative split-ting schemes can be understood in the same way as conservativecases, starting from the modified di ff erential equations yielded bythe BCH formula. Dissipative splitting schemes retain strong geo-metric properties (suppressing energy errors by a factor of e above),and conservation properties. For example, a dissipative perturba-tion that damps the radial component of particle velocities (and thusacts radially) would conserve angular momentum. Composing sucha step with conservative Kepler and interaction steps would stillconserve the total angular momentum to machine precision (sinceeach step does individually).While we have seen that dissipative splitting schemes can sys-tematically over(under)damp, we show in the following section thatthe same techniques from symplectic integration allow us to correctthese errors at no additional computational cost. case there would be no suppression of the phase errors. We consider phaseerrors more carefully in Rein et al. (2019a).MNRAS000
Second-order, Kepler-splitting integration of a single planet ini-tially on a 1-year, 0.1 eccentricity orbit, subject to a damping force thatcauses the orbital period to decay exponentially. The relative energy errorfollows the same estimate (dashed black line) as the conservative case inFig. 2, accounting for the fact that the orbital period T is changing. The secular rise then simply reflects the fact that the orbital pe-riod is changing exponentially as the planet migrates inward, whichyields a straight line on the linear-log scale.Remarkably, this suggests that exponential outward migrationwould yield errors that exponentially decrease with time. We findthat this is indeed the case, though this case is more subtle. In thecase plotted above where the errors committed each timestep in-crease with time, the most recent errors always dominate the er-ror budget. In the case where the instantaneous errors are gettingsmaller, one might expect the total error to remain at the level in-curred at the beginning of the integration, where errors are largest.But these errors are oscillatory, and as long as the errors are chang-ing adiabatically (i.e., as long as there are many orbits per migra-tion timescale), the oscillations can march toward smaller ampli-tude. Of course for exponential outward migration, it does not takemany migration timescales until the orbital period becomes com-parable to the migration timescale, at which point the applied forcein Eq. 14 (unphysically) becomes comparable to the central gravi-tational force, and we see numerical errors rise again.In summary, the error behaviour in weakly dissipative split-ting schemes can be understood in the same way as conservativecases, starting from the modified di ff erential equations yielded bythe BCH formula. Dissipative splitting schemes retain strong geo-metric properties (suppressing energy errors by a factor of e above),and conservation properties. For example, a dissipative perturba-tion that damps the radial component of particle velocities (and thusacts radially) would conserve angular momentum. Composing sucha step with conservative Kepler and interaction steps would stillconserve the total angular momentum to machine precision (sinceeach step does individually).While we have seen that dissipative splitting schemes can sys-tematically over(under)damp, we show in the following section thatthe same techniques from symplectic integration allow us to correctthese errors at no additional computational cost. case there would be no suppression of the phase errors. We consider phaseerrors more carefully in Rein et al. (2019a).MNRAS000 , 000–000 (0000) Daniel Tamayo, Hanno Rein, Pengshuai Shi, David Hernandez
Ignoring any improvements resulting from any symmetries in theproblem (Sec. 3.2), the splitting errors for a second-order Wisdom-Holman scheme should oscillate on orbital timescales with an am-plitude of ∼ (cid:15) ( h / T ) . However, these high frequency error oscilla-tions should be unimportant to the long-term evolution through theaveraging principle (Wisdom & Holman 1991).Wisdom et al. (1996) go further and show that these high-frequency oscillations of O ( (cid:15) ) can be e ffi ciently removed through anear-identity canonical transformation, which they term symplec-tic correctors, from the real action-angle variables to mapping vari-ables. In this picture, the short term oscillations intuitively arisefrom a mismatch in initial conditions. When going from the realsystem to the modified mapping, one has to correct the initial con-ditions to modified ‘mapping coordinates’. An integration is thenof the form C − ( h ) ◦ SABA ( h ) ◦ ... ◦ SABA ( h ) ◦ C ( h ) , (16)where each SABA ( h ) represents a timestep of the second-orderWisdom-Holman map (Eq. 6), and C ( h ) and C − ( h ) are the sym-plectic corrector transformations to and from mapping variables atthe beginning and end of the integration. Wisdom et al. (1996) fur-thermore prove that this canonical transformation removes all theerror terms in the modified di ff erential equations for the Wisdom-Holman map (Eq. 7) of O ( (cid:15) h n ) for all n . The leading error usingsymplectic correctors is then O ( (cid:15) h ).For generality and convenience of application, Wisdom et al.(1996) go on to show how the correctors C ( h ) for an arbitrary pairof integration steps A and (cid:15) B (hereafter we assume that A is thedominant operator and B is a perturbation) can be approximated toprogressively higher order through compositions of A and (cid:15) B withcarefully chosen timesteps forward and backward in time. Wisdom(2006) gives explicit coe ffi cients for such n th order approximations C n , which remove error terms O ( (cid:15) h k ) up to k = n in Eq. 7. Depend-ing on the problem, this leaves as the leading error the larger of thefirst uncorrected term O ( (cid:15) h n + ) and the first term at higher order inthe perturbation parameter O ( (cid:15) h ).We also briefly note that these correctors for the second-orderWisdom-Holman map SABA (Eq. 6) can also be used for the first-order scheme
SAB (Eq. 3) with a simple modification. As far asthe first-order scheme is concerned, the second-order scheme is areasonable approximation to the exact propagator P , so P ( h ) ≈ A ( h / ◦ B ( h ) ◦ A ( h / = A ( − h / ◦ A ( h ) ◦ B ( h ) ◦ A ( h / . ≡ C − ( h ) ◦ SAB ( h ) ◦ C ( h ) , (17)and thus C ( h ) = A ( h /
2) acts as a first-order corrector for the first-order splitting scheme to yield
SABA . Then, if desired, one coulduse a higher order correctors C n in Eq. 16. This trick of combininghalf steps of the second-order scheme in the middle of the integra-tion is used widely (in REBOUND this is referred to as synchroniza- We note that ordering matters once we specialize to systems where oneoperator is dominant. Wisdom et al. (1996) (with corrections in Wisdom2006) give correctors for the A ( h / ◦ (cid:15) B ( h ) ◦A ( h /
2) second-order scheme.The correctors for the alternative second-order scheme (cid:15) B ( h / ◦ A ( h ) ◦ (cid:15) B ( h /
2) would have di ff erent coe ffi cients. tion), but rarely described as a corrector that can be extended tohigher order .Finally, the fact that correctors compensate for the leading er-rors provides a straightforward way to interpret what the scheme isgetting wrong. For example, inverting Eq. 17 yields SAB ( h ) ≈ A ( h / ◦ P ( h ) ◦ A ( − h / , (18)i.e., to leading order, the first-order scheme A ( h ) ◦B ( h ) is equivalentto performing an exact step along the true solution P ( h ), exceptwith the mistake of taking additional forward and backward Keplerhalf-steps before and after the fact. While Wisdom et al. (1996) derived symplectic correctors throughphysically motivated canonical perturbation theory specific toHamiltonian systems, they point out that the formulas could alsohave been derived solely from the Lie algebra (commutators) ofthe two operators using Lie series and the BCH formula. But asmentioned above, the BCH formula is a statement about the com-position of non-commutative operators and has nothing to do withsymplecticity.‘Symplectic’ correctors therefore should correct any second-order splitting scheme involving a dominant operator and a pertur-bation. We now demonstrate explicitly that this is true even for dis-sipative perturbations. ‘Symplectic’ correctors are therefore morewidely applicable than their name implies. Specifically, they repre-sent ‘weak splitting correctors.’ In fact, like Wisdom et al. (1996),we rediscovered that similar ideas were applied as far back asButcher (1969) to general Runge-Kutta methods without restrictionto Hamiltonian systems.In particular, we take the same problem of a single planet mi-grating inward from Sec. 3.3, except we shorten the semi-majoraxis damping timescale τ a to 100 (initial) orbital periods, and in-tegrate for one damping timescale.We consider three integration schemes. First we apply the firstorder scheme K ( h ) ◦ D ( h ). From Secs. 2.4 and 3.3, this will yielderrors ∼ (cid:15) D ( h / T ) e . We then run a separate integration using first-order correctors (Eq. 17), whose error should scale as (cid:15) D ( h / T ) e ,and an integration with third-order correctors C (Eq. 16), whichshould scale as (cid:15) D ( h / T ) e for small h . The results are shown inFig. 4, with all curves following the above estimates (black dashedlines).Several features are apparent from Fig. 4. First, when thoughtof as a ‘symplectic’ corrector, it might seem surprising that sym-plectic correctors applied at the beginning and end of the integra-tion (when the energy has dissipated by a factor of 3 in between)would give a more accurate result. But as argued above, in our casethere is nothing specifically symplectic about them, they are sim-ply weak splitting correctors. Second, while all the timesteps con-sidered are small compared to the initial orbital period of 1 year,it is the final orbital period (vertical dashed line) that matters. Wesee that for timesteps (cid:38)
10% of the final orbital period, the errorrises dramatically, a consequence of higher order terms in the BCHseries becoming important.It might seem counterintuitive that symplectic correctors canimprove integrations involving dissipation. In particular, one might Jack Wisdom realized this long before we did (personal communi-cation).
MNRAS , 000–000 (0000)
EBOUNDx: A Library for Additional Forces in N-body simulations Timestep (Years) R e l a t i v e E n e r g y E rr o r h h h Roundoff dt F i n a l O r b i t a l P e r i o d T f First-Order SchemeCorrected Second-Order SchemeCorrected Fourth-Order Scheme
Figure 4.
Relative energy error in integrations analogous to Fig. 6 of a mi-grating planet, using ‘symplectic correctors’. Despite the dissipative forces,‘symplectic correctors’ fix the leading scheme errors to provide the ex-pected scaling. As discussed in the text, ‘symplectic correctors’ are reallysplitting correctors and provide a way to suppress energy errors by ordersof magnitude at fixed computational cost. The fourth-order scheme reachesroundo ff errors at small timesteps. worry that any steps backwards in time in C ( h ) would not can-cel with corresponding steps forward in time C −∞ ( (cid:104) ) when irre-versible processes are involved. For example, a step forward withfriction is not undone by a step backward in time—the second stepwould decrease the energy further. But this is not what is meantby, e.g., A ( − h ), i.e., it is not a step backward in time in the truesense where all the velocities are also flipped; rather, it merelymeans running the independent time variable run backwards, keep-ing everything else the same. Under this definition, a backwardsstep with frictional forces boosts the velocities, undoing a stepforward in time, and removing the contradiction. In the terminol-ogy of Hairer et al. (2006), the scheme is time-symmetric, butnot reversible. Hernandez & Bertschinger (2018) also find caseswhere time-symmetric, nonreversible schemes do not undergo en-ergy drifts (see their Fig. 8).The fact that ‘symplectic correctors’ can be derived directlyfrom the BCH formula implies that they must be more general,given that the BCH formula itself makes no distinction betweenconservative or dissipative operators. Finally, the fact that the greencurve in Fig. 4 reaches the theoretical round-o ff limit at ∼ − demonstrates that ‘symplectic correctors’ do improve integrationswith dissipation to high accuracy.This points out a limitation of splitting methods for dissipa-tive systems. Splitting methods require a fixed timestep (thoughsee, e.g., Preto & Tremaine 1999), and it must be shorter than thefastest timescale in the problem throughout the integration. If thereis substantial dissipation, the required timestep might be so smallthat a high-order adaptive scheme like IAS15 is more e ffi cient. Wenote that this concern also applies to the conservative N-body prob-lem. For example, if Kozai cycles or other e ff ects move planets ontovery eccentric paths at late times, the timestep for the whole inte-gration has to be chosen to match the shortest pericentre passagetimescale that occurs in the integration . There are ways to address this using time regularization.
The above discussion assumes that an analytical form for each stepcan be found, so that the only errors come from the splitting schemeitself, e.g., Eq. 7. However, for many perturbations, analytical solu-tions are not known, which naturally leads to the question of howerrors from numerically approximating the evolution across eachstep 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 theinteraction step, Eq. 10. In that case Eq. 11 remains the exact an-alytical solution, using the accumulated accelerations from all theposition-dependent e ff ects.However, many important astrophysical e ff ects like the mi-gration forces considered in the previous section, post-Newtoniancorrections, tides, etc., are velocity-dependent. Previous authorshave proposed incorporating such velocity-dependent perturbationsin either the Kepler step (Malhotra 1994), in the interaction step(Cordeiro et al. 1996), or as a separate step (Touma & Wisdom1994). 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 frameworkdiscussed above.In this case, the velocity-dependent accelerations vary acrossthe perturbation step as the velocities change, so Eq. 11 is no longerexact; it is now merely an Euler step approximation, accurate onlyto first order in h . In this section we consider the errors that suchan Euler step introduces, and compare it to higher-order numericalapproximations to the evolution across the perturbation step in thisgeneral case where the propagator cannot be found analytically. We first consider Hamiltonian velocity-dependent e ff ects, whichwe will find are qualitatively di ff erent from ones involving dissi-pation. We note that while we will find such perturbations causesevere numerical problems for these previously proposed methods,those authors were interested in dissipative perturbations, which weconsider 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 massratio. The equations of motion can be derived from a Hamiltonian(see Appendix B), whose conservation we can use to track the nu-merical accuracy.In all integrations we adopt the first-order splitting K ( h ) ◦ GR ( h ) , (19)where the corresponding di ff erential equations areˆ K z : ˙r i = v i , ˙v i = a Kepi ( r ) (cid:99) GR z : ˙r i = , ˙v i = a GRi ( r , v ) . (20)with the a GRi given in Eq. B5.We note that it is possible to find a splitting for this post-Newtonian Hamiltonian for which the evolution under each opera-tor can be solved analytically (Saha & Tremaine 1994). We never-theless choose to use the above splitting both to explore the e ff ectsof imperfect approximations across the perturbation step, and witha view toward making the REBOUNDx library a general-purpose toolfor integration. Such Hamiltonian velocity-dependent perturbationschange the relationship between the particles’ physical velocities
MNRAS , 000–000 (0000) Daniel Tamayo, Hanno Rein, Pengshuai Shi, David Hernandez and their momenta through Hamilton’s equations, so this relation-ship 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 toboth conservative and dissipative velocity-dependent forces. Thisis the same setup as in both Touma & Wisdom (1994) and Cordeiroet al. (1996); we now compare their choice of integrating across the GR step using a first-order Euler approximation (Eq. 11) to usinghigher-order methods.We first note that in this case of Hamiltonian velocity-dependent forces, Mikkola (1998) proposes using a low-order sym-plectic scheme like the implicit midpoint method to salvage thesymplecticity of the scheme. However, even a symplectic schemewill not be symplectic when applied the non-canonical position andvelocity coordinates of Eq. 20, like all the above authors consider .Therefore, one might not actually expect any advantage from usinga symplectic scheme across the perturbation step.We test this empirically on the K2-137 system (Smith et al.2017), for which post-Newtonian e ff ects are important. It consistsof a 0.89 Earth-radii ( R ⊕ ) planet on a 4.3 hour ( a ≈ .
01 AU)orbit around a 0 .
46 solar-mass primary. We assign the planet themass of the Earth. While in reality tides should have circularizedthe orbit, we inflate the initial orbital eccentricity to 0.01, in orderto see the apsidal precession induced by the relativistic e ff ects. Thepost-Newtonian perturbation strength is (cid:15) GR = GM / ac ∼ × − (e.g., Nobili & Roxburgh 1986).We compare integrating across the GR perturbation step withthe second-order, symplectic implicit midpoint method advocatedby Mikkola (1998) to three non-symplectic schemes: first-orderEuler (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 ≈ .
1% of the innermost planet’sorbital period. This implies a relative error of ∼ (cid:15) GR ( h / T ) e ∼ × − (black arrow).In the RK2, RK4 and implicit midpoint integrations in Fig. 5,the truncation errors remain much smaller than the splitting schemeerrors, and the simulations are visually indistinguishable, showingflat, oscillatory errors. As argued above, applying a symplectic al-gorithm like the implicit midpoint method to non-canonical vari-ables does not yield a symplectic scheme, so none of these inte-grations are exactly symplectic. Nevertheless, by ensuring that thetruncation errors across the perturbation step are much smaller thanthe splitting scheme errors, we only see the oscillatory symplecticbehavior.Of course, the most salient feature in Fig. 5 is that the error inthe Euler method of Touma & Wisdom (1994) and Cordeiro et al.(1996) (blue) error grows secularly, and in fact the planet is spuri-ously ejected after approximately 5.8 billion orbits, or in under 3million years. This behaviour too can be analytically estimated. As a simple analogy, consider a canonical (and thus symplectic) transfor-mation that rotates 2D Cartesian coordinates. If we instead consider a trans-formation that first converts to non-canonical polar coordinates, naively ap-plies the same symplectic rotation matrix to the polar coordinates, and thenconverts back to Cartesian coordinates, the sequence of transformations willnot be canonical. Time (Inner Planet Orbits) R e l a t i v e H a m il t o n i a n E rr o r GR ( h / T ) e EulerImplicit MidpointRK2RK4Eq. 28
Figure 5.
Integration of the ultrashort period planet K2-137b under rela-tivistic e ff ects with an inflated orbital eccentricity of 0.01. Naive incorpora-tion of the post-Newtonian accelerations into the Wisdom-Holman scheme(Euler method, blue) leads to a secular error growth, and the innermostplanet escapes the system within 6 billion orbits or less than 3 million years.Integrating across the perturbation step using higher order methods, whetherthey are symplectic (implicit midpoint) or not (RK2, RK4), all give oscilla-tory errors that remain constant over the length of the integration. For general perturbations that can not be integrated analytically,the local truncation error across the perturbation timestep must beincorporated into the error budget.The local truncation error ∆ z err LT across a single timestep usinga one-step, order- n method comes from the first neglected term atorder n + z = f ( t )(we focus on a single scalar coordinate z for simplicity), ∆ z err LT = f n + ( ξ )( n + h n + , (21)where f n + is the n + ξ is an unspecified time in the range [ t , t + h ].If we assume a solution f ( t ) = z exp( i π t / T pert ), with the char-acteristic perturbation timescale T pert a factor of (cid:15) longer than theorbital timescale T , we have ∆ z err LT z ∼ (2 π ) n + (cid:15) n + ( n + (cid:32) hT (cid:33) n + . (22)This is the error incurred every timestep. We can obtain theworst-case relative global truncation error after an integration time t by assuming the local errors add coherently, and multiplying Eq. 22by the number of steps t / h , ∆ z err LT z ∼ (2 π ) n + (cid:15) n + ( n + (cid:32) hT (cid:33) n (cid:32) tT (cid:33) (23)We overplot this estimate for the Euler method of Touma &Wisdom (1994) and Cordeiro et al. (1996) ( n =
1) in Fig. 5 and seeit gives a good match to the numerical behaviour. The truncationerrors for the second-order RK2 method would be lower by a factorof (cid:15) GR ( h / T ) ∼ − and thus never become visible.We might, however, expect di ff erent behavior under weak dis-sipation. To reach Eq. 23 we assumed the worst-case scenario thatthe one-step errors added coherently. At a more detailed level, onemust consider how such errors are propagated by the dynamicalflow itself. Consider a set of trajectories in phase space within adi ff erential error volume around the true trajectory. By the diver-gence theorem, the fractional growth of this volume element across MNRAS , 000–000 (0000)
EBOUNDx: A Library for Additional Forces in N-body simulations a timestep is given by the product of h and the divergence of thevector field of trajectories ∇· ˙z . In a Hamiltonian system like above,this divergence vanishes through Hamilton’s equations (this is Li-ouville’s theorem), so, at best, the dynamics are neutral for the ac-cumulation of errors .By contrast, under weak dissipation, the divergence of the flowis O ( (cid:15) ) and negative, so errors contract by O ( (cid:15) h ) each timestep asthe dynamics brings nearby trajectories together. One might there-fore expect that integrating across the perturbation step with even afirst-order Euler method would be good enough, given that one addsone-step errors ∼ O ( (cid:15) h ) (Eq. 22) slower than they are damped bythe dynamics. We test this empirically in the following section. Taking a somewhat realistic damping problem, we initialize twoJupiter-mass planets in a 3:2 mean-motion resonance around aSun-like star, with the innermost planet initially at 0.90 AU us-ing celmech . We then apply an exponential eccentricity damp-ing force prescription (Papaloizou & Larwood 2000) to both plan-ets with an e-folding decay timescale of τ e ≈ modify orbits forces implementation in REBOUNDx . The outer planet is then additionally acted on by an in-ward migration force (Eq. 14) with a timescale 100 times longer τ a ≈ . τ a estimate above (two equal massplanets have to be moved). We use the second-order scheme K ( h / ◦ D ( h / ◦ I ( h ) ◦ D ( h / ◦ K ( h / , (24)and integrate across each damping step using the four methodsof Sec. 5.1. We integrate for 4 × initial orbits, which trans-lates to approximately 4000 eccentricity-damping timescales, or2 semimajor-axis damping timescales. We adopt a small fixedtimestep of 10 − of the innermost planet’s initial orbital period toaccommodate the shrinking orbital periods, which decrease by afactor of ≈
25 over the integration. As above, we compare the en-ergies to the energies in a ‘perfect’ integration with
IAS15 . In reality, an initial parcel of trajectories within some error volumewould get sheared out at fixed volume celmech is a publicly available code for semi-analytic manipulations ofthe classical disturbing function in celestial mechanics. Among other things,it can perform transformations between orbital elements and resonant vari-ables: https://github.com/shadden/celmech We provide a simple argument why one can always use such atime-symmetric splitting with multiple operators to create a second-orderscheme. The timestep coe ffi cients for each operator must sum to unity inorder to match the modified di ff erential equations to the true ones at ze-roth order, and by writing an explicitly time-symmetric scheme, we ensurethat the additional terms in the modified di ff erential equations cannot de-pend on the sign of h (e.g., Eq. 7). The scheme is therefore second-order byconstruction. Time (Myr) R e l a t i v e E n e r g y E rr o r eulerimplicit_midpointrk2rk4 Figure 6.
Integration of two Jupiter-mass planets in a 3:2 mean-motion res-onance equilibrium between migration forces pushing the system deeperinto resonance (and both planets inward), and eccentricity damping push-ing the system out of resonance. Inner planet starts with an orbital periodof ≈ ≈ .
04 years by the end of the simulation. Weintegrate across the velocity-dependent damping step using four di ff erentmethods, and all give comparable results. Contrast with Fig. 5, where theEuler method fails. See text for discussion. In this case the scheme errors are dominated by the interplan-etary, resonant interactions. As in Sec. 3.3, we see in Fig. 6 thatthe energy errors grow due to the shrinking orbital periods. Like inFig. 5, the RK2, RK4 and implicit midpoint integrations are visu-ally indistinguishable. By contrast, while the integrations using thefirst-order Euler method across the damping step have higher er-rors, it is not a runaway e ff ect as in the conservative case in Fig. 5.As discussed in Sec. 5.2, this is because of the dissipative dynamicscontinually damping any accumulated energy errors.In summary, the crudest, most computationally e ffi cient Eu-ler step across a velocity-dependent force step should give reliableresults for dissipative forces, as advocated by Touma & Wisdom(1994) and Cordeiro et al. (1996). However, such a scheme wouldyield qualitatively wrong results in long-term integrations with con-servative velocity-dependent forces (Sec. 5.1). This can straightfor-wardly be solved by integrating across the perturbation step with ahigher method; for high-accuracy explicit splitting schemes that gobeyond this to correct even such higher-order integrations acrossthe perturbation step, see Blanes et al. (2013). In the previous sections we covered several technical aspects ofhow additional forces interact with splitting integration schemeslike
WHFast . However, many readers are likely interested in morepractical questions of how to quickly get something working. Wetherefore implement the integration schemes and ideas discussedabove in a new library,
REBOUNDx , which can be found at https://github.com/dtamayo/reboundx . REBOUNDx is library that provides tools and routines to ac-curately and self-consistently incorporate many di ff erent astro-physical e ff ects to N -body simulations. It seamlessly interfaceswith 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 abinary format to save
REBOUNDx configurations that interfaces with
MNRAS , 000–000 (0000) Daniel Tamayo, Hanno Rein, Pengshuai Shi, David Hernandez the
SimulationArchive class in
REBOUND , making it possible toshare and reproduce results bit by bit.We have already implemented several common forces in
REBOUNDx , some of which we used as examples in the discussionabove. We hope to increase the number of supported astrophys-ical e ff ects over time. REBOUNDx is an open-source project andwe welcome the community to contribute new e ff ects. We pro-vide many tutorials and examples in the form of Jupyter notebooksthat illustrate the library’s usage and how one can add new ef-fects with minimal e ff ort. We also share the scripts used to gen-erate 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 andlogic in the REBOUNDx package.
In a scenario where all the forces (including interplanetary forces)are always small compared to the gravitational forces from the cen-tral body, we recommended using Wisdom-Holman Kepler split-ting schemes. The implementation of this method in
REBOUND is WHFast (Rein & Tamayo 2015a). This splitting will typically befastest at reasonable levels of accuracy, and it provides strong con-servation properties in long-term integrations.Higher levels of accuracy are supported in cases where theforces can all be derived from position-dependent potentials by set-ting sim.ri whfast.corrector = N , where N is the order of thecorrector (Sec. 4.1). In velocity-dependent cases, correctors can beimplemented manually (the notebook fig4.ipynb in this paper’srepository provides an example). For even higher levels of accuracy(without significant computational cost), one can use the WHCKLmethod (Rein et al. 2019b), which is an implementation of the ker-nel method of Wisdom et al. (1996) (see also Wisdom 2018) thatalso removes errors at O ( (cid:15) ).In cases with close encounters between planets, the motionwill no longer be nearly Keplerian, and Kepler splittings will yieldlarge errors. In cases where such close encounters are rare, hybridintegrators are powerful tools (Duncan et al. 1998; Chambers 1999;Rein et al. 2019c). In REBOUND , the
MERCURIUS integrator retainse ffi ciency by using a Kepler splitting when possible, and switchesto an adaptive, high-accuracy integration with IAS15 for particlesundergoing close encounters (Rein et al. 2019c).
MERCURIUS sup-ports additional forces, but only switches to
IAS15 during closegravitational encounters. In the current implementation, any addi-tional forces must remain small throughout the integration when
MERCURIUS is used.When the above conditions are not met, the safest option (andthe default) in
REBOUND is to use the
IAS15 integrator (Rein &Spiegel 2015). Because it is an adaptive scheme,
IAS15 can handlearbitrary forces and still yield solutions that are accurate close to themachine precision.
IAS15 will typically be slower in cases where
WHFast and
MERCURIUS can be used. However,
IAS15 will alwaysbe more accurate and is applicable to a wider variety of problems . 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).
REBOUNDx
Structures
The two principal ways to incorporate additional e ff ects into an N -body simulation are through forces and operators. One way to implement an additional e ff ect is to write it explicitlyas a force. This force then contributes to the accelerations of allparticles, in addition to the standard Newtonian gravity.When the IAS15 integrator is used, the integrator automat-ically adapts the timestep such that integrations will be accurateclose to machine-precision. If one of the various splitting integra-tors available within
REBOUND (for example WH , MERCURIUS , and
LEAPFROG ) is used, then the additional forces are accumulated andapplied in the integrator’s interaction (or kick) step using Eq. 11.As discussed above, this is appropriate for position-dependent ordissipative 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 potentialsare supported.
Rather than calculating accelerations which are then integrated inthe interaction or kick step of the integrator, operator steps yieldsolutions. Operators update the positions and velocities (or evenmasses or other particle parameters) to the value they should haveat the end of a specified timestep.Operators (e.g., D ) together with a timestep (e.g., h /
4) form astep (e.g., D ( h / ff erent timesteps to construct higher-order splittingschemes (e.g., Eq. 24).For cases where the solution for the positions and velocities atthe end of the step is not known, or di ffi cult to calculate analytically,the REBOUNDx object integrate force operator can take a force,and integrate it across the timestep using any of the four integratorsdiscussed in Sec. 5. The default is RK4. This is the recommendedmethod for conservative, velocity-dependent forces (Sec. 5.1).
All forces and operator steps are stored in linked lists such thatthey are executed in the reverse order they were added. Thus, toget a scheme consisting of the steps K ( h ) ◦ D ( h ), the operator K should be added before D , i.e. left to the right, since the last stepacts first. REBOUNDx provides a simple API for adding e ff ects into N-bodyintegrations from the list of implemented options, demonstratedin the numerous examples in the documentation. However, italso provides lower-level functionality for customized tasks. Thethree 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 at-tached forces to update particle accelerations. MNRAS , 000–000 (0000)
EBOUNDx: A Library for Additional Forces in N-body simulations • pre timestep modifications . This iterates through thelinked list of attached operator steps to execute before the main REBOUND step. • post timestep modifications This iterates through thelinked list of operator steps to execute after the main
REBOUND step.The user populates these lists by either adding forces and operatorsalready implemented in
REBOUNDx , or through user-defined func-tions in either
PYTHON or C . By default, adding an operator to anintegration adds half steps before and after the main REBOUND stepto make a second-order scheme.In Fig. 7, we schematically outline a full
REBOUND timestepusing the
WHFast integrator. Di ff erent integrators switch out the reb integrator part1 and reb integrator part2 steps, butin all cases they rely on the gravitational and additional forces cal-culated in between them. In contrast to many other N-body packages,
REBOUND allows theuser to set up the simulation in an arbitrary inertial frame. The co-ordinates are returned to the user in the same inertial frame.The various splitting integrators in
REBOUND use a number ofdi ff erent coordinates internally (e.g. Jacobi coordinates). In orderto minimize logic and pitfalls, all REBOUNDx functions are alwayscalled after converting the particles’ positions and velocities backto Cartesian coordinates in the inertial frame. Force and operatorimplementations thus do not need to worry about indirect or non-inertial forces. This also implies that it is typically desirable to en-sure that the net applied forces vanish. This is discussed in detail inAppendix C.
Additional e ff ects in an N -body simulation typically come withseveral associated parameters (e.g., a migration timescale for adamping force). REBOUNDx provides an interface which allows auser to attach an arbitrary numbers of parameters to forces, oper-ators, and particles. The implementation is such that parametersremain tied to their parent structures even as objects get removedor reshu ffl ed in the simulation. ff ects A current list of implemented e ff ects can be found in the documen-tation. At the time of this writing, the following e ff ects are avail-able: • Migration forces (Papaloizou & Larwood 2000) • Eccentricity and inclination damping forces (Papaloizou & Lar-wood 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 are pre-timestepmodifications present? convert coordinatesto inertial frameapply pre-timestepmodificationsWHFast Kepler step(half timestep)calculate gravita-tional accelerationsare additionalforces present? convert coordinatesto inertial framecalculateadditional forcesWHFast interaction step(full timestep)WHFast Kepler step(half timestep)are post-timestep mod-ifications present? convert coordinatesto inertial frameapply post-timestepmodificationssearch for collisionsand ejections YesNo YesNo YesNo reb integrator part1reb integrator part2reb update acceleration
Figure 7.
Schematic outline of a full
REBOUND step when using the
WHFast integrator. Functions which are part of
REBOUNDx and called via functionpointers are shaded.
Throughout this paper, we have developed the error behaviour ofsplitting schemes in a general, non-Hamiltonian framework. Thisapproach clarifies how several properties typically ascribed to sym-plectic integrators also carry over to weakly dissipative systems.Many astrophysically relevant e ff ects have strong symmetryproperties. We showed explicitly that the remarkable degeneracyof the Kepler problem makes schemes that split the evolution intoKepler steps and perturbation steps exploit these symmetries to of-ten strongly suppress energy errors in the integration, as in the casesof radial or drag forces applied to nearly circular orbits (Sec. 3.2).We also showed that so-called ‘symplectic correctors’, which re-duce energy errors by orders of magnitude at fixed computationalcost (Wisdom et al. 1996), apply equally well to weakly dissipa- MNRAS000
Throughout this paper, we have developed the error behaviour ofsplitting schemes in a general, non-Hamiltonian framework. Thisapproach clarifies how several properties typically ascribed to sym-plectic integrators also carry over to weakly dissipative systems.Many astrophysically relevant e ff ects have strong symmetryproperties. We showed explicitly that the remarkable degeneracyof the Kepler problem makes schemes that split the evolution intoKepler steps and perturbation steps exploit these symmetries to of-ten strongly suppress energy errors in the integration, as in the casesof radial or drag forces applied to nearly circular orbits (Sec. 3.2).We also showed that so-called ‘symplectic correctors’, which re-duce energy errors by orders of magnitude at fixed computationalcost (Wisdom et al. 1996), apply equally well to weakly dissipa- MNRAS000 , 000–000 (0000) Daniel Tamayo, Hanno Rein, Pengshuai Shi, David Hernandez tive systems and can thus be more generally thought of as ‘weaksplitting correctors’ (Sec. 4.2).We also explored the general case where individual opera-tor steps can’t be solved analytically, and one is forced to inte-grate across the perturbation step. We showed that the previouslyadvocated approaches of incorporating additional forces into theWisdom-Holman map (Touma & Wisdom 1994; Malhotra 1994;Cordeiro et al. 1996) work well for dissipative e ff ects (Sec. 5.3),but give qualitatively wrong answers for conservative, velocity-dependent forces like post-Newtonian corrections of general rel-ativity (Sec. 5.1).Finally, we described REBOUNDx , an open-source C libraryfor incorporating additional e ff ects into REBOUND
N-body integra-tions, together with a convenient
PYTHON wrapper. Users can ei-ther choose to add an already implemented astrophysical forces,or easily implement new e ff ects. We hope that this modular, open-source framework for robustly incorporating new e ff ects into a va-riety of N-body integration schemes will help the community torun more accurate and realistic simulations. We encourage othersto contribute to this library of astrophysical e ff ects. ACKNOWLEDGMENTS
We would like to thank John Chambers for an insightful and con-structive review that greatly improved this manuscript. We wouldalso like to thank Jack Wisdom and Scott Tremaine for help-ful discussions. Support for this work was provided by NASAthrough the NASA Hubble Fellowship grant HST-HF2-51423.001-A awarded by the Space Telescope Science Institute, which isoperated by the Association of Universities for Research in As-tronomy, Inc., for NASA, under contract NAS5-26555. This re-search has been supported by the NSERC Discovery Grant RGPIN-2014-04553 and the Centre for Planetary Sciences at the Univer-sity of Toronto Scarborough. This research was made possible bythe open-source projects
Jupyter (Kluyver et al. 2016), iPython (P´erez & Granger 2007), and matplotlib (Hunter 2007; Droett-boom et al. 2016).
REFERENCES
Anderson J. D., Esposito P. B., Martin W., Thornton C. L., Muhleman D. O.,1975, ApJ, 200, 221Batygin K., Morbidelli A., 2012, The Astronomical Journal, 145, 1Benitez F., Gallardo T., 2008, Celestial Mechanics and Dynamical Astron-omy, 101, 289Blanes S., Casas F., Farres A., Laskar J., Makazaga J., Murua A., 2013,Applied Numerical Mathematics, 68, 58Brown E. W., 1931, Monthly Notices of the Royal Astronomical Society,92, 104Burns J. A., Lamy P. L., Soter S., 1979, Icarus, 40, 1Butcher J., 1969, in Conference on the numerical solution of di ff erentialequations. pp 133–139Chambers J. E., 1999, Monthly Notices of the Royal Astronomical Society,304, 793Cordeiro R., Gomes R., Martins R. V., 1996, Celestial Mechanics and Dy-namical Astronomy, 65, 407Danby J., 1992, Richmond: Willman-Bell,— c1992, 2nd ed.Droettboom M., et al., 2016, matplotlib: matplotlib v1.5.1,doi:10.5281 / zenodo.44579, http://dx.doi.org/10.5281/zenodo.44579 Duncan M. J., Levison H. F., Lee M. H., 1998, AJ, 116, 2067Forest E., Ruth R. D., 1990, Physica D: Nonlinear Phenomena, 43, 105 Fujii M., Iwasawa M., Funato Y., Makino J., 2007, Publ. Astron. Soc. Japan,59, 1095Galley C. R., 2013, Physical review letters, 110, 174301Galley C. R., Tsang D., Stein L. C., 2014, arXiv preprint arXiv:1412.3082Goldreich P., Schlichting H. E., 2014, AJ, 147, 32Gr¨obner W., 1967, Die lie-reihen und ihre anwendungen. Vol. 3, DeutscherVerlag der WissenschaftenHairer E., Lubich C., Wanner G., 2006, Geometric numerical integra-tion: structure-preserving algorithms for ordinary di ff erential equa-tions. Springer Science & Business MediaHernandez D. M., 2019a, arXiv preprint arXiv:1904.03364Hernandez D. M., 2019b, MNRAS, 486, 5231Hernandez D. M., Bertschinger E., 2018, Monthly Notices of the RoyalAstronomical Society, 475, 5570Hernandez D. M., Dehnen W., 2017, Monthly Notices of the Royal Astro-nomical Society, 468, 2614Hunter J. D., 2007, Computing In Science & Engineering, 9, 90Hut P., 1981, Astronomy and Astrophysics, 99, 126Kinoshita H., Yoshida H., Nakai H., 1990, Celestial Mechanics and Dynam-ical Astronomy, 50, 59Kluyver T., et al., 2016, Positioning and Power in Academic Publishing:Players, Agents and Agendas, p. 87Laskar J., 1989, Nature, 338, 237Laskar J., Gastineau M., 2009, Nature, 459, 817Laskar J., Robutel P., 2001, Celestial Mechanics and Dynamical Astron-omy, 80, 39Lithwick Y., Wu Y., 2012, The Astrophysical Journal Letters, 756, L11Machin J., 1737, Philosophical Transactions of the Royal Society of Lon-don, 40, 205Malhotra R., 1994, Celestial Mechanics and Dynamical Astronomy, 60, 373Mikkola S., 1998, Celestial Mechanics and Dynamical Astronomy, 68, 249Mikkola S., Innanen K., 1999, Celestial Mechanics and Dynamical Astron-omy, 74, 59Mikkola S., Tanikawa K., 1999, Celestial Mechanics and Dynamical As-tronomy, 74, 287Naoz S., Kocsis B., Loeb A., Yunes N., 2013, The Astrophysical Journal,773, 187Neri F., 1987, Dept. of Physics, University of MarylandNewhall X. X., Standish E. M. J., Williams J. G., 1983, A&A, 125, 150Newton I., 1687, Philosophiae Naturalis Principia Mathematica. Auctore Js.Newton, doi:10.3931 / e-rara-440.Nobili A. M., Roxburgh I. W., 1986, in Relativity in Celestial Mechanicsand Astrometry. High Precision Dynamical Theories and ObservationalVerifications. p. 105Papaloizou J., Larwood J., 2000, Monthly Notices of the Royal Astronomi-cal Society, 315, 823Pelupessy F., Portegies Zwart S., 2012, Monthly Notices of the Royal As-tronomical Society, 420, 1503P´erez F., Granger B. E., 2007, Computing in Science and Engineering, 9,21Petit A. C., Laskar J., Bou´e G., Gastineau M., 2019, arXiv preprintarXiv:1905.13240Plummer H., 1896, Monthly Notices of the Royal Astronomical Society, 56,317Portegies Zwart S., 2018, Science, 361, 979Portegies Zwart S., McMillan S., 2018, Astrophysical Recipes; The art ofAMUSE, by Portegies Zwart, Simon; McMillan, Steve. ISBN: 978-0-7503-1321-6. IOP ebooks. Bristol, UK: IOP Publishing, 2018Preto M., Tremaine S., 1999, The Astronomical Journal, 118, 2532Rambaut A. A., 1890, Monthly Notices of the Royal Astronomical Society,50, 301Rauch K. P., Holman M., 1999, The Astronomical Journal, 117, 1087Rein H., Liu S.-F., 2012, Astronomy & Astrophysics, 537, A128Rein H., Spiegel D. S., 2015, MNRAS, 446, 1424Rein H., Tamayo D., 2015a, Monthly Notices of the Royal AstronomicalSociety, 452, 376Rein H., Tamayo D., 2015b, MNRAS, 452, 376Rein H., Tamayo D., 2019, Research Notes of the AAS, 3, 16MNRAS , 000–000 (0000) EBOUNDx: A Library for Additional Forces in N-body simulations Rein H., Brown G., Tamayo D., 2019a, arXiv preprint arXiv:1908.03468Rein H., Tamayo D., Brown G., 2019b, arXiv e-prints, p. arXiv:1907.11335Rein H., et al., 2019c, Monthly Notices of the Royal Astronomical Society,485, 5490Ruth R. D., 1983, IEEE Trans. Nucl. Sci., 30, 2669Saha P., Tremaine S., 1992, The Astronomical Journal, 104, 1633Saha P., Tremaine S., 1994, arXiv preprint astro-ph / 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 thesame notation, we seek a solution to the di ff erential equations˙ z = L P z , (A1)where here we define the Lie derivative more carefully. The Liederivative with respect to a set of di ff erential equations ˙ z = ˆ F z actson functions g ( z ) to return their total derivative. Since throughoutthe paper we treat only cases that are explicitly time-independent,we have L F g ( z ) = (cid:88) i ∂ g ∂ z i ˙ z i , (A2)where the ˙ z i are given by ˆ F z .With this definition, we can Taylor expand the solution aroundthe current time , z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t + h = z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t + h d z dt (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t + h d z dt (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t + ... ≡ (cid:34) (1 + h L P + h L P + ... ) z (cid:35) t . (A3)As a concrete example, consider the simple one-dimensional dif-ferential equation ˙ z = z . We have L P z = z , L P z = L P z = z etc.. Plugging into Eq. A3 yields the Taylor expansion of the exactsolution z ( t ) = (1 / z (0) − t ) − .By analogy to the Taylor expansion for the exponential, this isoften written as z ( t + h ) ≡ e h L P z ( t ) ≡ P ( h ) z ( t ) . (A4)In words, P ( h ) is a time evolution operator, i.e., an integrator, that Throughout the paper we assume that the di ff erential equations do notdepend explicitly on time. If they did, an integral would be required in theexponential of Eq. A4. updates the phase space variables across a timestep. We then pro-ceed to split the di ff erential equations into two pieces ˙ z = L P z = ( L A + L B ) z as in the main text.One subtle notational issue that has led to some errors in theliterature is that the BCH formula is typically written as e h L A e h L B = e h L S (A5)with L S = L A + L B + h L A , L B ] + h
12 [ L A −L B , [ L A , L B ]] + O ( h ) . (A6)However, e h L A e h L B is not equivalent to the splitting scheme onemight naively implement on a computer e h L A ◦ e h L B .In e h L A e h L B z , e h L A acts on a function g ( z ) = e h L B z , and there-fore introduces partial derivatives of g through Eq. A2. By contrast,the integrator e h L A ◦ e h L B first evaluates an intermediate z (cid:48) = e h L B z ,and these updated values are used to yield a final answer of e h L A z (cid:48) .The BCH formula applies to the former sense (one can straightfor-wardly expand the three exponentials in Eq. A5 to obtain the firstterms in Eq. A6), but this interpretation is not particularly usefulfor writing integrators. Under iterated mappings, the repeated chainrules would quickly become unwieldy. By contrast, e h L A ◦ e h L B is straightforward to implement on a computer, and fortunatelycan be simply related to the scheme in the BCH formula to ana-lyze its error properties. This is sometimes referred to as the Ver-tauschungssatz (Gr¨obner 1967) (see Lemma 5.1 in chapter III ofHairer et al. (2006)); we provide a simple derivation and examplehere.Looking back at Eq. A3, one could just as well have Taylorexpanded any scalar function f ( z ), and concluded in Eq. A4 thatthe time evolution operator e h L A also evolves f ( z ) from time t totime t + h (with f ( z ) = z in Eq. A3 as a special case). Thus, sinceEq. A5 is of the form e h L A f ( z ), with f ( z ) = e h L B z , we must have e h L A f ( z ) = f ◦ ( e h L A z ) , (A7)i.e., evolving f ( z ) across the timestep (LHS) is equivalent to firstevolving z and applying the function f to the updated variables(RHS). Therefore, e h L A ◦ e h L B = e h L B e h L A . (A8)This is the reason for the sign flips at all odd orders in Eq. 5.As one might expect, Eq. A8 has led to some confusion innotation 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 orpointed out is that most integrators used in the astrophysics liter-ature are time reversible, and thus will not have any sign flips intheir corresponding BCH formulas given that they are invariant toflipping the order of operations (e.g., Yoshida 1993 and the detailedanalysis in Hernandez & Dehnen 2017).As an explicit example of the above, consider a non-dimensionalized simple harmonic oscillator with position x , mo-mentum v , and z = ( x , v ). The Hamiltonian is H = v / + x / H T = v / H V = x / ff erential equations are given by Hamilton’sequations, L T z : ˙ x T = ∂ H T ∂ v = v , ˙ v T = − ∂ H T ∂ x = , (A9) L V z : ˙ x V = ∂ H V ∂ v = , ˙ v V = − ∂ H V ∂ x = − x . (A10)(A11) MNRAS000
12 [ L A −L B , [ L A , L B ]] + O ( h ) . (A6)However, e h L A e h L B is not equivalent to the splitting scheme onemight naively implement on a computer e h L A ◦ e h L B .In e h L A e h L B z , e h L A acts on a function g ( z ) = e h L B z , and there-fore introduces partial derivatives of g through Eq. A2. By contrast,the integrator e h L A ◦ e h L B first evaluates an intermediate z (cid:48) = e h L B z ,and these updated values are used to yield a final answer of e h L A z (cid:48) .The BCH formula applies to the former sense (one can straightfor-wardly expand the three exponentials in Eq. A5 to obtain the firstterms in Eq. A6), but this interpretation is not particularly usefulfor writing integrators. Under iterated mappings, the repeated chainrules would quickly become unwieldy. By contrast, e h L A ◦ e h L B is straightforward to implement on a computer, and fortunatelycan be simply related to the scheme in the BCH formula to ana-lyze its error properties. This is sometimes referred to as the Ver-tauschungssatz (Gr¨obner 1967) (see Lemma 5.1 in chapter III ofHairer et al. (2006)); we provide a simple derivation and examplehere.Looking back at Eq. A3, one could just as well have Taylorexpanded any scalar function f ( z ), and concluded in Eq. A4 thatthe time evolution operator e h L A also evolves f ( z ) from time t totime t + h (with f ( z ) = z in Eq. A3 as a special case). Thus, sinceEq. A5 is of the form e h L A f ( z ), with f ( z ) = e h L B z , we must have e h L A f ( z ) = f ◦ ( e h L A z ) , (A7)i.e., evolving f ( z ) across the timestep (LHS) is equivalent to firstevolving z and applying the function f to the updated variables(RHS). Therefore, e h L A ◦ e h L B = e h L B e h L A . (A8)This is the reason for the sign flips at all odd orders in Eq. 5.As one might expect, Eq. A8 has led to some confusion innotation 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 orpointed out is that most integrators used in the astrophysics liter-ature are time reversible, and thus will not have any sign flips intheir corresponding BCH formulas given that they are invariant toflipping the order of operations (e.g., Yoshida 1993 and the detailedanalysis in Hernandez & Dehnen 2017).As an explicit example of the above, consider a non-dimensionalized simple harmonic oscillator with position x , mo-mentum v , and z = ( x , v ). The Hamiltonian is H = v / + x / H T = v / H V = x / ff erential equations are given by Hamilton’sequations, L T z : ˙ x T = ∂ H T ∂ v = v , ˙ v T = − ∂ H T ∂ x = , (A9) L V z : ˙ x V = ∂ H V ∂ v = , ˙ v V = − ∂ H V ∂ x = − x . (A10)(A11) MNRAS000 , 000–000 (0000) Daniel Tamayo, Hanno Rein, Pengshuai Shi, David Hernandez
The propagators can easily be found by integrating their respectivedi ff erential equations, or by noting L T n z = L V n z = n (cid:62)
2, so e h L T z = (1 + h L T ) z = (cid:32) x + hvv (cid:33) (A12) e h L V z = (1 + h L V ) z = (cid:32) xv − hx (cid:33) . (A13)We first evaluate the composition used in the SAB integrator(Eq. 3), e h L T ◦ e h L V z , where e h L T z (cid:48) is evaluated at z (cid:48) = e h L V z ( e h L T ◦ e h L V ) z = (cid:32) x + h ( v − hx ) v − hx (cid:33) . (A14)By Eq. A8, this should be equivalent to e h L V e h L T z = (1 + h L V ) (cid:32) x + hvv (cid:33) = (cid:32) x + hvv (cid:33) + h (cid:32) ˙ x V + h ˙ v V ˙ v V (cid:33) , (A15)which indeed matches Eq. A14. APPENDIX B: FIRST ORDER POST-NEWTONIANCORRECTIONS
General relativistic corrections to Newtonian dynamics are detailedthroughout the literature, but since we use these e ff ects as one ofthe main tests for the various schemes in this paper, we present theequations and details of our implementation for completeness.The first order post-Newtonian corrections (1PN) are O ( v / c ), where v is the characteristic orbital velocity and c is thespeed of light. We have implemented various levels of approxima-tion. In gr full , we include the full (1 PN) equations of motion,which are implicit and computationally expensive, but would beappropriate for equal-mass bodies.For planets around a single star, one can typically ignore cor-rections 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 e ff ect in REBOUNDx , and is the case discussed in Sec. 5.As a final level of approximation, one is typically interested inthe apsidal precession e ff ect from general relativity, in which caseit is possible to write down a potential only depending on the par-ticles’ position, which, in an orbit-averaged sense gives the correctprecession rate (Nobili & Roxburgh 1986). This is implementedin gr potential , and has the advantage that it can be integratedanalytically (Sec. 3.1), and is therefore computationally cheapest.However, this approximation introduce errors in the phases and in-stantaneous elements of O ( v / c ). Typically, gr potential is suf-ficient for planetary applications around a single star since it cap-tures the correct secular behavior, but both gr potential and gr are only appropriate for single-star systems. For higher multiplicitysystems, one must use gr full .Naoz et al. (2013) provide the full Hamiltonian for a triplesystem (see also Sch¨afer 1987). The equations of motion are givenby Newhall et al. (1983) and Benitez & Gallardo (2008). Theseaccelerations 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 e ff ect in more detail sinceit di ff ers from previous works, and is the subject of Sec. 5. The fullHamiltonian can be written as a sum of the Newtonian and 1PNHamiltonians, H = H N + H PN , where H N = (cid:88) i p i m i − (cid:88) j (cid:44) i Gm i m j r ij (B1) and (Saha & Tremaine 1994), H PN = c (cid:88) i (cid:44) (cid:32) µ m i r i − p i m i − µ p i m i r i (cid:33) (B2)Having already ignored terms of O ( m i / m ), the di ff erence betweenbarycentric, Jacobi or democratic heliocentric coordinates (and Ja-cobi or physical masses) in H PN is negligible. Interpreting all po-sitions and momenta as Jacobi coordinates considerably simplifiesthe equations of motion, since the kinetic term in H N remains a di-agonal sum of the momenta. Application of Hamilton’s equationsyields (e.g., Saha & Tremaine 1994), ˙r = ∇ p H = p m ≡ ˜v , ˙r i (cid:44) = ∇ p i H = p i m i + ∇ p i H PN ≡ ˜v i (1 + A i ) , (B3)where A i = − c (cid:32) ˜v i + µ r i (cid:33) (B4)and the ˜v i ≡ p i / m i are pseudo-velocities not equal to the physicalvelocities ˙r i . We note that Eq. B3 would have additional terms indemocratic heliocentric coordinates due to kinetic cross terms in H N . We then have ¨r = ˙˜v = , ¨r i (cid:44) = ˙˜v i (1 + A i ) + ˙ A i ˜v i (B5)where ˙˜v i (cid:44) = − m i ∇ r i H = a i + r i c (cid:32) µ ˜v i r i − µ r i (cid:33) , (B6) ˜a i are the net Newtonian accelerations on particle i , and˙ A i = − c (cid:32) ˜v i · ˙˜v i − µ r i ( r i · ˙r i ) (cid:33) . (B7)We see that the accelerations (Eq. B5) are explicitly velocity-dependent, and must be treated with care in symplectic schemes(Sec. 5). As expected, ¨r = ˜v i , but this is easily accomplished iteratively,since the 1PN approximation is only appropriate when the pertur-bation is weak, and the physical and pseudo-velocities are approxi-mately equal. Operationally, if like in REBOUND , the integrator doesnot work in Jacobi coordinates, one simply calculates the Newto-nian accelerations in inertial coordinates, transforms to Jacobi coor-dinates, calculates the Jacobi accelerations ¨r i and transforms backto inertial coordinates. Rein & Tamayo (2015a) provide unbiasedtransformation algorithms between inertial and Jacobi coordinates. APPENDIX C: CENTRE OF MASS
When inserting additional e ff ects in N-body simulations in an iner-tial frame, it is often valuable to ensure that there are no net forcesacting on the system’s centre of mass (COM). For example, theWisdom-Holman integration scheme (Kinoshita et al. 1990; Wis-dom & Holman 1991; Rein & Tamayo 2015a) assumes that theCOM moves at constant velocity and removes that degree of free-dom from the problem. This e ff ectively means that any residualnet force is incorrectly distributed among the particles’ relative co-ordinates. However, even integration schemes that track the COMdegree of freedom can be adversely a ff ected by finite forces on the MNRAS , 000–000 (0000)
EBOUNDx: A Library for Additional Forces in N-body simulations barycentre. Such a simulation’s accuracy will continually deterio-rate as the system moves away from the origin and the precision onthe inter-particle separations degrades due to subtracting ever largernearly equal numbers. These errors may be acceptable dependingon the application, but can slow down schemes like IAS15 (Rein& Spiegel 2015) that try to reach machine precision as they fail toconverge with progressively smaller timesteps.It is therefore numerically desirable to self-consistently modelall back-reactions from additional e ff ects so that the net force van-ishes. While this is trivial for forces between pairs of particles inthe simulation, one does not always wish to fully model all compo-nents.Consider a net force F D from a protoplanetary disk of mass M D on a planet. In response, the disk should feel an equal and op-posite net force, but there is no disk “particle” in the simulation onwhich to apply it. This uncompensated force will yield an acceler-ation on the COM. However, the disk is tightly coupled to the cen-tral star gravitationally, suggesting that this back-reaction shouldbe e ff ectively communicated to the central body. More quantita-tively, if we take a characteristic orbital radius for the disk r D ,then we can define a characteristic timescale for the back-reaction τ BR ∼ ( r D / a D ) / , where a D is the net acceleration of the disk F D / M D . At the same time, the disk is gravitationally coupled tothe primary on the Keplerian orbital period τ K = π ( GM / r D ) − / .If τ K (cid:28) τ BR , the central body will respond adiabatically to theaccelerations of the disk, and the back-reaction could instead beadded to the central star.In other cases it might be better motivated to add the accel-erations to the centre of mass of all or a subset of the particles,e.g., back-reactions onto a circumbinary disk. This corresponds toapplying an acceleration a = − F / M tot to the appropriate set of par-ticles, 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 ofthe Hamiltonians that govern the dynamics. Being able to preciselycalculate these Hamiltonians provids a practical check on the nu-merical accuracy of the implementation for conservative systems.A Hamiltonian that depends on the radial positions r i and r j of twoparticles, H ( r i − r j ) exerts equal and opposite forces on the twobodies. Similarly, one can show that for a Hamiltonian dependenton the distance from the system’s barycentre H ( r i − r COM ), where r COM = M − tot (cid:80) j m j r j , Hamilton’s equations dictate that if parti-cle i feels force F i = −∇ r i H ( r i − r COM ), all particles (including i )should feel an additional acceleration − a COM = − F i / M tot . Finally,for a Hamiltonian with positions referenced to the centre of mass ofa subset of particles, Hamilton’s equations require that only thoseparticles feel the reverse acceleration − F i / M sub , where M sub is themass of the subset of particles. In the case where each particle’sposition is referenced to the centre of mass of all interior bodies,this corresponds to the familiar Jacobi coordinates. MNRAS000