Multiple timestep reversible $N$-body integrators for close encounters in planetary systems

We present new almost time-reversible integrators for solution of planetary systems consisting of"planets"and a dominant mass ("star"). The algorithms can be considered adaptive generalizations of the Wisdom--Holman method, in which all pairs of planets can be assigned timesteps. These timesteps, along with the global timestep, can be adapted time-reversibly, often at no appreciable additional compute cost, without sacrificing any of the long-term error benefits of the Wisdom--Holman method. The method can also be considered a simpler and more flexible version of the \texttt{SYMBA} symplectic code. We perform tests on several challenging problems with close encounters and find the reversible algorithms are up to $2.6$ times faster than a code based on \texttt{SYMBA}. The codes presented here are available on Github. We also find adapting a global timestep reversibly and discretely must be done in block-synchronized manner or similar.


INTRODUCTION
Orbital dynamics has been essential to the development of astronomy as a field.Newton's Universal Gravitation (Newton 1687) provides a set of second order differential equations for the rate of change of gravitational systems.Given initial conditions, if we are able to solve the set of equations, we can obtain the past or future (since the equations are time-symmetric) evolution of any gravitational system.We have been able to effectuate this procedure to gain insight into the formation of galaxies, or understanding of the longterm stability of the Solar System; indeed, entire subfields of astronomy have been created as a result.
The mathematical problem we are describing is known as the Nbody problem, in which N point masses interact through pairwise forces: it has attracted substantial attention since Newton from scientists such as Laplace, Lagrange, and Poincaré.The N-body problem has no closed-form solution in general.One particular challenge of this problem is that solving it requires a multidisciplinary effort: first, we need to develop powerful astronomical tools capable of measuring the initial conditions needed, the theory of differential equations has to be developed to solve the system, chaotic dynamics must be understood to make sense of the calculations and the randomness involved, and finally, we need the development of technology (computers) which can rapidly carry out all the calculations.
The requirement of all these ingredients has meant many of the most interesting astronomical problems have only been in reach in the last few decades.One of the least understood parts of this process of solving the N-body problem is the role chaos plays in the accuracy of solutions (Valtonen 1976;Smith 1977;Heggie 1991;Quin-lan & Tremaine 1992;Portegies Zwart & Boekholt 2014;Boekholt & Portegies Zwart 2015;Hernandez et al. 2020).We will not focus on this in this paper, but rather on the theory of the differential equations.
Dynamical astronomy made large progress due to the advent of symplectic methods (Dragt & Finn 1976;Yoshida 1990;Channell & Scovel 1990;Kinoshita et al. 1991;Wisdom & Holman 1991;Saha & Tremaine 1992;Yoshida 1993;Calvo & Sanz-Serna 1993;Wisdom et al. 1996;Chin 1997;Preto & Tremaine 1999;Mikkola & Tanikawa 1999a,b;Laskar & Robutel 2001;Yoshida 2001Yoshida , 2002;;Chambers et al. 2002;Chin & Chen 2005;Hairer et al. 2006;Rein & Tremaine 2011;Farrés et al. 2013;Blanes et al. 2013;Hernandez & Bertschinger 2015;Hernandez 2016;Petit et al. 2019).Symplectic integrators solve Hamiltonian systems like N-body systems; they conserve exactly Poincaré invariants and are particularly well-suited to study systems undergoing many orbits over long timescales, and are thus used to probe the evolution of gravitational systems like galaxies or the Solar System.Unfortunately, it is difficult to construct symplectic integrators that require changing timescales, such as in collisional systems with short relaxation times (Binney & Tremaine 2008).Conventional integration methods like the Runge-Kutta-Fehlberg method (Press et al. 2002), Bulisch-Stoer method (Press et al. 2002), or the Hermite integrator (Makino 1991) are far better suited for collisional systems with close encounters, but they lack the long-term error benefits of the symplectic methods.Ideally we could combine the long-term advantages of symplectic methods with the adaptability of conventional methods.Despite the tendency of symplectic methods to fail for collisional systems, with few good alternatives, they have remained popular.
Focusing in on the field of planetary dynamics, there have been several widely used, transformative symplectic methods.The Wisdom-Holman method (WH) (Wisdom & Holman 1991) is wellsuited for nearly Keplerian systems like the Solar System.Chambers (1999); Rein et al. (2019) develop a hybrid symplectic integrator that switches between WH and another method when near-Keplerian motion ceases.Duncan et al. (1998) (hereafter DLL98) offer another alternative, in which WH takes pairwise adaptive steps.While not perfect, these methods have significantly advanced the study of planetary dynamics.
Time-reversible integrators are an alternative to symplectic integrators, and in many cases perform just as well in long-term simulations, while being as flexible as conventional methods (Hut et al. 1995;Funato et al. 1996;Kokubo et al. 1998;Holder et al. 2001;McLachlan & Perlmutter 2004;Hairer et al. 2006;Makino et al. 2006;Hairer et al. 2009;Dehnen 2017;Hernandez & Bertschinger 2018;Boekholt et al. 2022).As long as the method recovers the initial conditions when integrating forwards and then backwards, it is time-reversible.Time-reversible methods have been largely ignored in practice, possibly due to their often large expense or the difficulty in implementing them.It has recently been shown, however, that "almost" (which we refer to simple as time-reversible) time-reversible algorithms can switch between two methods, often at no penalty to compute time or accuracy (Hernandez & Dehnen 2023).Hernandez & Dehnen (2023) used their simple time-reversible algorithm to devise a simpler, more flexible version of the symplectic hybrid integrator developed by Chambers (1999); Rein et al. (2019) (with implementation known as MERCURY/MERCURIUS).Currently, the Hernandez & Dehnen algorithm is being implemented (Lu et al., in prep) as part of the REBOUND package (Rein & Liu 2012).Among its advantages over MERCURY/MERCURIUS, for example, it can successfully integrate close encounters with the dominant star mass.
In this work, we similarly develop simple, fast, and flexible timereversible versions of the DLL98 algorithm, implemented as SYMBA.
To accomplish this, we can generalize the algorithm of Hernandez & Dehnen (2023) to allow time-reversible switching between an arbitrary number of integrators (and thus timesteps).The advantage of our approach is in the speed and simplicity of the method.We develop two time-reversible methods, MTR (Multiple Timestep Reversible) and AG (Adaptive Global), and implement a SYMBA-like symplectic integrator, MTS (Multiple Timestep Symplectic), for fair comparisons.In all our comparisons, the reversible algorithms are at least as accurate as MTS and one of the reversible methods is always faster.This all holds while maintaining the benefits of more flexible adaptive conventional integrators.For easy reference, the codes used in this work are provided at a Github link 1 .We warn the reader that they are not in a user friendly form and provided merely for reference.We have tested the reversible methods on challenging problems involving close encounters throughout the paper.Our codes can be considered as proof of concept, as they are not yet implemented in a faster compiled language like C; this is left for future work.
In Section 2, we provide background on developing integrators and on notation.In Section 3, we develop both symplectic and timereversible multiple timescale and timestep algorithms suitable for the two-body Kepler problem, with numerical demonstrations.In Section 4, we generalize these methods to systems with multiple bodies with multiple timesteps varying independently.We focus on systems with one dominant mass, the "star."We conclude in Section 5.

BACKGROUND ON SYMPLECTIC AND TIME-REVERSIBLE INTEGRATORS
For this paper, we consider conservative Hamiltonian systems with form, where p and q are the conjugate momenta and positions, respectively.This form describes many orbital mechanics problems in astronomy.More compactly, we describe these coordinates via z(t) = (p(t), q(t)), with t the time.The equations governing them are, with formal solution, {} are Poisson brackets.For functions a(z) and b(z), the following definition holds: The Lie operator of H is Ĥ and is defined Ĥa = {a, H}.Solution (3) is not often practical because it is not obtainable in closed-form.Our task in this paper is to find more efficient and practical approximate solutions.One approach is to take advantage of the function separation in eq. ( 1) to obtain the "map," where z 1 indicates a first-order approximation in time2 .As an alternate map, switch T and V The solution of eq. ( 5) is simple: This map is symplectic because it preserves phase space volumes (Hairer et al. 2006) and, as a result, long-term orbits obtained with it are more reliable and accurate.Also, map ( 6) is derived from differential equations themselves derived from a function H1 , which is just H to 0th order in t.H1 has form, That the lowest power in t is one indicates it is a first-order method.Map (5) is not time-reversible: integrating backwards does not recover the initial conditions.In fact, any odd order method cannot be time-reversible.Time-reversible methods are better than nonreversible methods at obtaining long-term accurate orbits (Hairer et al. 2006).For this work, we will not distinguish between timesymmetric and time-reversible methods.A both symplectic and time-reversible map instead is shown as follows and known as leapfrog: Again, an alternate map switches T and V: h is now chosen as a small time interval h = t − t 0 , the "timestep," with t 0 the initial time.The associated function H2 is now, and h must be small enough to ensure H2 remains close to H. A long-term solution is obtained by repeatedly applying map (8).

Phase space dependent timesteps
Leapfrog is no longer symplectic or time-reversible if we try to adapt the stepsize as a function of phase space h = h(z(0)).The symplecticity of map ( 8) is defined by the condition, where J is the Jacobian of z(h) with respect to z(0), and Ω is a constant anti-symmetric matrix whose form depends on the ordering within z.When h is no longer a constant, the condition ( 11) is generally destroyed.That time-reversibility is lost is also seen immediately.We show how to solve these problems in Section 3.1.

MULTIPLE TIMESTEP INTEGRATORS FOR THE KEPLER PROBLEM
3.1 Symplectic algorithm Skeel & Biesiadecki (1994) showed a strategy to essentially adapt leapfrog's timestep in a symplectic way, by decomposing the interaction potential into a sum of components, each treated with a different timestep.We describe the implementation of this idea by Lee et al. (1997).We consider a 2D Kepler Hamiltonian.Using Hamiltonian (1), choose the functions, and with q = q 2 x + q 2 y .Then we select a set of cutoff radii r 1 > r 2 > ... > r n with constant ratio r i /r i+1 = R.Each cutoff radius range, and by extension, potential range, has a corresponding timestep, such that if q > r 1 , h = h 0 ; if r 1 > q > r 2 , h = h 1 ; and so on until h = h n−1 .Also, h i /h i+1 = M and the global timestep is h 0 .Now we decompose the potential into, where . ., with associated forces, Fi = −∂ Ṽi /∂q.When an operator Vi appears, we substitute in terms of Vj .We focus on the form of these forces adopted in the tests of DLL98: where f (x) = 2x 3 − 3x 2 + 1. Fi transitions smoothly from 0 to −q/q 3 .This polynomial has the properties f (1) = d f dx 0 = d f dx 1 = 0 and f (0) = 1.This form ensures V(q) has continuous derivatives, and thus the symplecticity of methods containing operators Vi is ensured (Hernandez 2019a).Fig. 1 of DLL98 and Fig. 3 of Lee et al. (1997) plot the force decompositions F i for forces in a single dimension.Now, define, We approximate the φi using leapfrog methods φi , such that φi = φi + O(h 2 ), but the two sets of operators in eq. ( 8) are now φi+1 and Vi or Vim and T : φ0 is the Lee et al. (1997) map.i m is a maximum recursion level which may vary during evaluation of φ0 .More details are provided in Section 3.1.1.Note each φi is M leapfrogs with stepsize h i = h 0 /M i .In principle, the recursion of map ( 17) is infinite, but we can truncate it at the i m if we can be sure q > r im+1 during global step h 0 .
At best, we can only estimate the maximum h i during a global step; we outline such an estimate in Section 3.1.1.That we can only estimate the maximum h i is a disadvantage of this method that will be avoided with the reversible scheme derived in Section 3.2.Map (17) allows us to take adaptive stepsizes while maintaining symplecticity and forms the basis of the SYMBA code (DLL98).We have independently implemented map (17) and name it MTS for Multiple Timestep Symplectic.DLL98 derive a function H from which the differential equations of the map can be derived, and is close to H: H Lee err constitutes the error Hamiltonian, which must be kept small.Define F i = −∂V i /∂q = Fi − Fi−1 .Assume at some point in time q ≈ r i+1 .Then, from eq. ( 15), F i = −q/q 3 , while We have q ∝ p −1/2 , which implies, from eq. ( 19), that H Lee err = O(h 2 i q −4 ).Thus, choosing R = √ M makes H 2 independent of the close encounter strength and H Lee err = O(h 2 0 ).Another logical choice is to set M = 3 and R = 2, which yields h ∝ q 3/2 , a timestep proportional to the free-fall time, so that the phase is resolved evenly during free-fall.This gives H Lee err = O(h 2 0 q −1 ), adequate if the encounter is not extremely close.In addition, the smoothness of f ensures the h 2 i term is defined.Making f even smoother (in the sense that higher order derivatives of f are 0 at x = 0 and x = 1) can increase the accuracy of the map (Hernandez 2019b), and increases the number of terms at higher powers of h i well-defined in H.

Estimating the maximum recursion level
The code SYMBA estimates the maximum h i required during a global step (17), as explained above.We describe this implementation here, and also implement it in MTS.For full details, refer to our code 1 .If any of the following conditions are satisfied before evaluating φi , solve this operator and require recursion at φi+1 .Otherwise, set i = i m .
If min(r min , q) < r i " is the condition.
Our estimates only use positions and velocities/momenta, and no higher time derivatives, but work well enough in practice.One inefficiency of our implementation MTS is that once we have entered a recursion map φi , the maximum recursion level for the remainder of the global timestep h 0 must be at least i and cannot decrease.A more efficient algorithm would allow i to decrease.Our goal is to implement a simple and accurate MTS method, and the focus of this paper is on reversible, not symplectic, methods.

Adapting pairwise timesteps
We now describe an algorithm which is the first novelty of this paper.We wish to have a simpler multiple timestepping algorithm which avoids smoothing functions and can adapt the time-step not only based on the position.The algorithm will be time-reversible but non-symplectic, which is just as good in many cases.We call it MTR for Multiple Timestep Reversible.This algorithm is a generalization of reversible maps from Hernandez & Dehnen (2023) which only switch between two global timesteps.MTS will handle multiple adaptive pairwise steps.First, we substitute the timestep criteria r i for timestep criteria g i .g(q, p) is a useful possibly momentumdependent function such as a free-fall time or a semi-major axis.The timestep level i is then determined by the condition g i+1 < g(q, p) < g i .As long as we set g 0 , the g i follow from the adopted R value.Simply setting g i = r i works for a variety of problems and, indeed, in this paper we assume this choice unless otherwise stated.Define the map, which is composed of M i substeps and evolves for global step h 0 .A timestep level j k is collected after each substep k.Then, the condition for using Ĉi and We collect timestep levels at three points.For Ĉ2 , we would collect nine timestep levels.On occasion, global timesteps must be repeated once or even more when one of the h j collected at the substeps is less than h i .An algorithm demonstrating MTR is Listing 1.The while loop rarely requires more than one pass; it makes sure the maximum level has not changed after redoing a step.The initial timestep level at time 0 is calculated from the initial conditions.Every time MTR is deciding between two Ĉi maps, it is subject to rare reversibility errors described in Hernandez & Dehnen (2023); the so called "ambiguous" and "inconsistent" steps.
Listing 2 shows an alternate implementation of MTR, which gives identical results.We denote it MTRa.zeta is given by eq. ( 9) and is a function of the state z 1 and substep h 0 /M i .The level function calculates the timestep level.We present MTRa to aid in understanding MTR.
Listing 1: Python script for the MTR time-reversible algorithm, the first main result of this paper.MTR redoes timesteps as needed to enforce time-reversibility.

Adapting the global step itself
So far, we have kept h 0 fixed, but now we show how to adapt it reversibly.In a system with a hierarchy of timesteps (see Section 3), perhaps we could set h 0 as the lowest timestep being used by any pair.In a first example, we vary h 0 and do not use any other timestep, which is sufficient for the Kepler problem.We need only consider the leapfrog method, eq. ( 8), φ0 , which we write as mapphi0 in our algorithm.We still increase or decrease h 0 via a factor M. h 0 is easily reduced now when necessary, but we have found that increasing h 0 must be done with care.h 0 is only increased if the number of h 0 steps taken modulo M is 0. This algorithm does not exactly change steps at block synchronized points.We have found this requirement via experimentation.Dehnen (2017), Section 6.1.1,makes brief remarks about a similar finding, and provides a possible explanation.
Listing 3 shows explicitly how this algorithm works and is our second main result, labeled AG for Adaptive Global.lev is an array of the number of steps taken each timestep level, modulo M.
Listing 3: Python script for the AG time-reversible algorithm, the second main result of this paper.In contrast to MTR, AG varies the global timestep reversibly, and does not allow an increase of the stepsize at every point in time, which is required to maintain long-term accuracy.

Numerical comparisons among algorithms
3.3.1 0.9 eccentricity We now compare the performance of the symplectic MTS method and reversible MTR and AG methods.It is most advantageous to use MTR for systems with more than two planets, as presented in Section 4.2 , but we test it here for completeness.We consider the 2D Kepler problem, equations ( 1), ( 12), and (13).We set R = √ 2 and M = 2 a suggested in Section 3.1, which is a steep scaling.Also, we set the semi-major axis as a = 1, r 1 = √ 2, the eccentricity as e = 0.9, and h 0 = P/20000, with P the period (= 2π).R, M, e, and h 0 are chosen to match the parameters of Lee et al. (1997), Fig. 4. The A Kepler orbit with e = 0.9 is integrated for 1000 periods P, and the energy error is recorded over time.We compare results to the SYMBA-like MTS method.The errors are all stable over time and the reversible methods are the fastest algorithms.
initial conditions are at apocenter.The energy error as a function of time is plotted in Fig. 1.Here E is the energy (eq.( 1)) so ∆E/E is the relative energy error.The evolution is calculated up to time t = 1000P at 10000 linearly spaced outputs.
The algorithms are implemented in the MATLAB language and the compute times for MTS, MTR, and AG are 23, 13, and 9 s, respectively.The relative compute times are a useful guide for the result in faster compiled languages.Of course, these numbers depend on the processor being used but the important quantity to note are relative compute times, which are a fair comparison.The reversible methods have a speed advantage.
All methods keep a small energy error over time which is not drifting secularly.The errors of MTR and AG are not identical, although they appear close.The tendency of the error to be more positive or negative is a function of the phase of the orbit of the initial conditions, which we verified.MTR redid a fraction 4 × 10 −4 of the timesteps (800 out of 2000001 steps).The fraction for AG was 6 × 10 −3 (8000 out of 13309460 steps).The smallest timestep used by any algorithm was h 9 = P/1024000, near pericenters.For MTR, the timestep level between adjacent global timestep jumps by ±1 always.The global timestep level for AG also jumps by ±1 always.

0.999 eccentricity
We next want to test extreme eccentricity orbits, which will require extremely small steps near pericenter.The motivation for this extreme orbit is that Lee et al. (1997) also test extreme eccentricity orbits.Except for setting e = 0.999, we use the same parameters and algorithms as in Section 3.3.1.This test will require timestep level changes by more than one for MTR.This happens because timestep levels near pericenter can vary significantly during a global timestep.The MTS maximum recursion algorithm of Section 3.1.1becomes important for efficiency and catching rapid timestep level changes.We expect the efficiency of MTR to be worse for e = 0.999 compared to e = 0.9 since the timestep levels cannot be varied within a global step.The runtimes are 2.7 × 10 3 , 170, and 270 s for MTR, AG, and MTS, respectively, so the runtime of AG is 60% that of MTS, while MTR is about 10 times slower than MTS.The fastest method is a reversible method.We plot the energy error over time for the methods in Fig. 2. The errors over time for all methods again remain small.In this case, the error amplitudes for AG and MTR are slightly larger in magnitude than those of MTS.The median errors are × 10 −8 , −2.0 × 10 −7 , and 2.0 × 10 −7 for MTS, MTR, and AG, respectively: the reversible methods have absolute median errors 3.7 times that of MTS.The smallest timestep used is h 0 /2 22 (the 22nd timestep level).For MTR, timestep level changes by up to ±10 between global timesteps occur.

MULTIPLE TIMESTEP ALGORITHMS FOR PLANETARY SYSTEMS WITH A DOMINANT MASS
In Section 3 we showed that we can reversibly adapt the timestep in a single orbit with two different approaches, MTR and AG.AG was faster than a symplectic SYMBA-like algorithm with adaptive steps.
In this section, we generalize the algorithms to work for systems with multiple bodies which have a dominant mass which we call "Sun."The other bodies are called "planets."For simplicity, we allow close encounters between planets only, and we don't allow close Sun-planet encounters.For problems close encounters with the Sun, solutions are presented elsewhere (Levison & Duncan 2000, Lu et al., in prep).Each pair of planets will be evolved on different timescales, so that we can think of each pair of planets having a different timestep.We begin with our 3D Hamiltonian of form (1): Here, m i are masses, G is the gravitational constant, and q i j = |q i − q j |. (q, p) are the inertial Cartesian canonical coordinates.Although we could proceed using Hamiltonian (22), it has the disadvantage that N-body maps based on it perform poorly in frames that are not at the center of mass.Thus, we recast Hamiltonian ( 22) in Democratic Heliocentric Coordinates (DLL98), denoted by (Q, P).Q 0 is the position of the center of mass and P 0 is the total momentum.For i > 0, Q i are heliocentric positions and P i are baryentric momenta.In these coordinates, Q 0 is cyclic and we can ignore a bulk kinetic term P 2 0 /(2M) in the Hamiltonian, with M = i m i .The total Hamiltonian is H = H Kep + H Sun + V, where , and Here, H Sun is the Solar kinetic energy.We are now ready to construct a second order map that solves this Hamiltonian, known as the Wisdom-Holman map (Wisdom & Holman 1991).In analogy to eq. ( 8) for leapfrog, we have, Because {H Sun , V} = 0, the order we apply ĤSun and V in is irrelevant.The Wisdom-Holman map ( 24) is efficient for studying planetary systems in which the planetary orbits remain well separated and thus Solving H Kep involves employing an incremental implicit Kepler equation solver.We use the approach in Universal Variables, such that the equations of motion no longer have singularities at Q i = 0. Our implementation is described in Wisdom & Hernandez (2015).

Symplectic integrator (DLL98)
We now allow all the pairs of planetary bodies in eq. ( 24) to be evolved on different timescales.While this was first described in DLL98, here we provide more details and fill in gaps in the derivation.Each pairwise distance will now be classified into a radius range (and thus, timestep level), which transitions smoothly into other shells using eq.( 15).First, rewrite V in analogy to eq. ( 14) as, Here, k is a shell index and (i, j) labels the planetary pair.Now, we update map (17) using the fact that, for any indices i, j, k, l, m, n: k m is the maximum recursion level which, as in Section 3.1.1,may vary during a global timestep h 0 .W k is the vector of indices i such that for > k, j i V li j = 0. Here, (i, j) pairs indicates a sum over all pairs (i, j) in some order which does not change the integrator due to eq. ( 26).As in eq. ( 16), we define, ψi = exp where ψi = ψi + O(h 2 ).Similar to the situation in Section 3.1.1,the maximum recursion levels within ψ0 can only increase in our implementation.We have implemented eq. ( 27), which is the SYMBA map, only for three-body systems, and call this implementation MTS again.We do not implement this algorithm for systems with more particles as the estimates for the pairwise maximum recursion levels (Section 4.1.1)have complications and our focus here is not on symplectic algorithms.

Estimating the maximum recursion level
For three bodies, we can apply the algorithm of Section 3.1.1 to determine the maximum h i for the planetary pair.Eq. ( 25) simplifies to V k12 , and We also make the substitutions For systems with more bodies, it's unclear a pairwise maximum recursion level can be adapted within the global step while maintaining symplecticity and time-reversibility, so we do not attempt it.The SYMBA code may have achieved this although we have not studied their implementation in detail.

Time-reversible integrators
As in Section 3.2, as a next novelty, we now develop a simpler, more flexible time-reversible version of the algorithm of Section 4.1.First, given z(0), identify the set of all pairs of particles S k at each timestep h k (based on some criteria like separation).The set at the smallest timestep is S km .Next, define maps, Ṽ is defined from eq. ( 25).The sum in Âk is over all pairs in S k .The same index can be present in different S k .S ′ k is a 1D array of all the indices in S k , but if any indices are present in S j for j > k, they have been removed.In this way, Bk consists of Kepler operators for particles not treated at smaller timesteps h j (for j > k).It also follows that, (for j > k) We now write down an evolution map, D0 generalizes map (20).Due to property (30), D0 is timesymmetric.We can now write generalizations of the reversible algorithms MTR and AG, Listings 1 and 3, that work with more than two bodies.For MTR, mapC is given by eq. ( 31).Before ( 31) is applied, we construct a timestep level matrix iN0, indicating each pair's timestep level.iN0 has size n(n − 1)/2, with n the number of planets.For problems with large n, not considered in this paper, manipulating an array of this size may become expensive, and a more efficient computational approach could be devised.We apply (31) and timestep levels are computed at each substep for each planet pair.The maximum level for each pair is recorded as iN, substituting Line 07 of Listing 1.If any level has increased compared to iN0, the step is redone.If close encounters occur frequently, one can expect the compute effort to, at most, double compared to an algorithm that does not redo steps.This is because redoing steps more than once is a rare case in our tests.After all maps are applied, iNt is calculated at the final coordinates and momenta, to be used at the beginning of the next global timestep.
To generate AG, in Line 07 of Listing 3, j must be the minimum over all planetary pairs.The leapfrog method, mapphi0, is replaced by the Wisdom-Holman map, eq. ( 24).Adapting the global timestep alone with AG will only be efficient for problems in which one pair of planets has close encounters at a time (more precisely, we consider it a close encounter if the pair timestep is smaller than h 0 ).For problems in which different pairs of planets have simultaneous close encounters, for efficiency, we should use MTR.An even more efficient approach would be to simultaneously adapt the pairwise steps with MTR and adapt the global step with AG, as needed.For many problems which are reasonably stable, a global timestep can be kept fixed, but if the system settles into a state with slower or faster timescales, adapting h 0 may make sense.

Comparing the error analysis for two algorithmic
approaches for close encounters

Hybrid symplectic integrators
In this work, we have focused on decreasing pairwise timesteps during planetary close encounters.An alternate symplectic strategy was presented by Chambers (1999), which we describe here, and we study its error.Chambers (1999) describes a hybrid symplectic integrator to resolve close encounters.The idea is to use a cheap map M1 when there is no close encounter and switch to a more expensive map M2 during close encounters.M1 can be written, where Â1 = ĤKep and B1 For M1 , we find, and it follows H M 1 err = O(ϵh 2 ).During a very close encounter, this error scaling becomes only O(h 2 ) so a hybrid method switches to map M2 , where, and Â2 = ĤKep + (i, j)∈W Vij and B2 = ĤSun + (i, j) W Vij .W is the set of pairs in a close encounter.Switching to M2 allows hybrid symplectic integrators to maintain excellent error behavior.A disadvantage is that A 2 needs to be solved via expensive methods like Bulirsch-Stoer.

Symplectic multiple timestep algorithms
It is straightforward to write the error Hamiltonian of ψ0 ; this was also derived in eq. ( 35) of DLL98: In the absence of close encounters, H DLL98 err = O(ϵh 2 0 ).We now further simplify eq. ( 35), as we did in eq. ( 19).The force experienced by i due to j is found by, During a close (i, j) encounter, the dominant contribution to the error Hamiltonian is, Note the difference with eq. ( 19).We can see that ), adequate if the close encounter is not extreme.

Discussion
The hybrid and multiple timestep integrators are two different algorithmic approaches to ensuring the error of the Wisdom-Holman method remains small during planetary close encounters.The approach of hybrid methods is more straightforward, but requires use of expensive conventional integrators like Bulirsch-Stoer during close encounters.Conceptually, it is unappealing to use a conventional integrator within a symplectic integrator.Multiple timestep algorithms can rely entirely on Kepler solvers and need not invoke conventional integrators.Conceptually, it is also unappealing to take the Wisdom-Holman method, which is designed under the assumption of no close encounters occurring, and to use brute force and decrease its timestep, to make it able to handle close encounters.Clearly, there is no perfect solution to the problem of planetary close encounters in N-body symplectic simulations.

Numerical comparisons among algorithms
We now present numerical experiments that demonstrate the performance of the reversible methods for systems with more than one planet (with a dominant mass).

The restricted three-body problem
As one of the simplest chaotic systems with more than one planet, we consider the planar, restricted three-body problem, in which there is one conserved quantity, the Jacobi integral.Our test problem is taken from Wisdom (2017), and consists of the "Sun," "Jupiter," on a circular orbit, and a test particle with initial conditions such that it can visit the other two masses.The test particle has multiple close encounters with Jupiter.For example, in one integration, there were 34 close encounters in the first 1000 years.Here, a close encounter is defined by the (only) timestep level jumping from 1 to a higher level, and back down to 1.A timestep level of 1 is the global timestep h 0 .As there is only one planetary pair in this problem, we can also apply AG to it.We set h 0 = 8 days, run for 10000 yrs, set r 1 = 5r Hill , with r Hill the Hill radius of Jupiter, and M = 4 and R = 2 Figure 3.A comparison of the error performance of the symplectic and reversible methods for a planar restricted three-body problem, in which the test particle can visit primary and secondary and has many close encounters with secondary.We record the error in the Jacobi constant over 843 secondary periods.The errors of all methods are comparable and the fastest method is the reversible AG method.
for all methods.Output is generated at one year intervals (Jupiter's period is 11.86 yrs).We compare the error in the Jacobi integral over time for MTS, MTR, and AG in Fig. 3.The errors of all methods do not drift significantly and have similar median errors, respectively: −2.3 × 10 −7 , −9.5 × 10 −7 , and 1 × 10 −6 .The runtimes are 256, 280 and 145 s, respectively: the fastest method is again AG with compute time about 57% that of MTS.Note we are deliberately avoiding comparing compute times with codes like MERCURIUS or SYMBA, since our implementations are not written in fast languages.

Violent Solar System
Next, we will test a problem demonstrating the remarkable difference a few irreversible timesteps can have on the accuracy of the evolution of a dynamical system.We consider a violent outer Solar System, introduced by DLL98.We reuse the following parameters: the system configuration of Sun plus outer giant planets, with planetary masses scaled up by a factor 50; h 0 = 0.03 yrs for MTR, and a runtime of 3000 yrs.To match the DLL98 experiment as closely as possible, we use position dependent shells with r 1 = 1.52 au, about 3 initial Saturn Hill radii.r 1 is used for all pairs.We set R = 2 as in DLL98.However, in contrast to DLL98, we use M = 4 rather than their M = 3.We were unable to achieve good accuracy with their choice for both MTR and AG.This can have many explanations; in particular as we note below, our dynamical evolution of this chaotic system differs from that of DLL98; note our initial conditions are different and the Lyapunov time in this problem can be short; e.g., 40 yrs (Hernandez 2016).We have also discarded the possibility that higher accuracy can be attributed to the maximum recursion level algorithm of Section 3.1.1.To check this, we implemented a higher accuracy version of MTR in which we not only collect the separations between substeps, but also the minimum separation within substeps, to determine timestep levels.There was no improvement in any result.We also tested a Kepler problem and found that checking minimum separations within substeps in no case altered recursion levels for the symplectic method.So checking separations within timesteps leave both the symplectic and reversible .Error over time of a time integration of a violent outer Solar System, consisting of the Sun and outer giant planets, whose masses have been scaled up by a factor 50.We compare the MTR method ("Reversible"), in which 8 steps have been repeated (a fraction 8 × 10 −5 ), to the same method without any repetition ("Naive").While "Naive" reaches an error of order unity, MTR maintains an error to about 2 parts in a million from redoing the 8 steps.
results unaltered in our tests.From previous experience, we do not expect that SYMBA's smooth switching function, which contrasts to our discrete switching function, is responsible for differences in accuracy.
In Fig. 4, we compare the performance of MTR (labeled "Reversible") and a naive, non-reversible method (labeled "Naive"), which is simply MTR with no repeated steps.MTR only redoes 8 steps in the entire simulation, a fraction of 8×10 −5 , and is almost like the naive method in the sense that it rarely repeats steps.Although the number of close encounters is small for this tests, we know from our long-term extreme eccentricity Kepler orbit tests of MTR (Section 3.3.2) that good accuracy is maintained even with more frequent encounters.The MTR error is similar to that obtained with SYMBA as reported by Fig. 6 of DLL98.MTR's error oscillates mostly between 1 and −2 × 10 −6 , while theirs oscillates between about 2 and −4 × 10 −6 .For the naive method, the error blows up to order unity by the end of the integration.The naive and reversible curves are identical until t = 281 yrs.
To more carefully understand the dynamical behavior, Fig. 5 plots the six planetary pairwise distances as a function of time for MTR.Saturn's initial Hill radius is indicated with a dashed line.The Jupiter-Saturn separation reaches a minimum of < 0.077 au at 282 yrs (because the separation is measured at output, 0.077 au is an upper bound).This is 160 Jupiter radii.MTR attempts a timestep level of 6, corresponding to h = 15.9 min.Saturn is ejected at ≈ 430 yrs, which contrasts to the result in DLL98, where it is ejected at 1950 yrs.Uranus and Neptune have a series of three close encounters in which the timestep level jumps up to 2.Besides these four events, no other timestep levels ever jump past 1.In contrast to the result of DLL98, Uranus is never ejected.As discussed above, local Lyapunov times for this system can be small, and our initial conditions likely differ from those of DLL98, so the different dynamical evolution is not surprising.

Planetary system with hierarchical binaries
As a final test, we wish to simulate a complicated system in which different pairs of particles are at different varying timestep levels simultaneously.We use MTR to study this system.AG alone would not give us the flexibility to adapt all pairwise steps.To find such a system, we build on the binary planet problem introduced by DLL98, Section 6.1.We have a Solar mass star with one binary planet placed at 1 au and another binary planet at 3 au, for a total of five bodies.The first binary has a = 0.0125 au and e = 0.6, while the second binary has a = 0.013 au and e = 0.2.All planets have mass 10 −3 M ⊙ .As in DLL98, Section 6.1, we set h 0 = 0.01 yr and integrate for 100 yrs, so that the total number of periods in the system for the different hierarchical binaries are 19, 100, 3017, and 3200.Although we could construct velocity-dependent timestep levels for this problem, in contrast to the case of SYMBA, we find instead that timestep levels based on the relative free-fall times t ff work well.For this system, we prefer this to timestep levels based on Hill radii as proposed by DLL98, because the Hill radius is less well motivated.For pair of planets (i, j), t ff = [Q 3 i j /(G(m i + m j ))] 1/2 .To get the timestep levels, use shells of a nondimensional function g = t ff /h 0 , and let g i /g i+1 = R, in analogy to the radius ratios explored in Section 3.1.We set g 1 = 30, which was determined numerically to work well.We set R = 2 and M = 3. Fig. 6 plots the energy error over time again for a reversible and a naive method.For the naive method, the error grows linearly in time by about four orders of magnitude.The error is well preserved for the reversible method.Fig. 7 displays a and e for the two planetary binaries in one year using the reversible method.The elements do not drift significantly.We also show the timestep levels of the two binaries.Between global timesteps, the inner binary's attempted step level varies by ±3, between 5 and 8.The inner binary actually uses timestep levels between 6 and 8 that change by at most ±2.For the outer binary, the variation is ±1, between 6 and 7, both for attempted and used timestep levels.All other planetary pairs remain at level 1 Figure 6.Error over time for an integration of a hierarchical five-body problem, with the two methods of Fig. 4. The system consists of a Solar mass star, a binary planet at 1 au, and another binary planet at 3 au, with both pairs in eccentric orbits.The pairwise timestep levels vary substantially, as shown in Fig. 7.The integration time is 100 yrs, over which the binary planets take about 3200 orbits.The reversible method conserves energy to better than a part in a million, while the naive method's error is worse than 0.1%.(= h 0 ).67% of steps were redone, while a 10 −4 steps were redone No steps were redone more twice.Thus, showed reversible method can handle timestep level changes by more than ±1 in a global timestep, and pairwise timestep levels can vary independently, all without incurring secular energy error drifts.Finally, we comment that the reversible method we used for this system might be made more efficient still by combining Listings 1 and 4, to adapt the pairwise steps locally in time; exploration of this is left for future work.

CONCLUSIONS
We presented two time-reversible algorithms, AG (Adaptive Global) and MTR (Multiple Timesteps Reversible).MTR can be understood as a Wisdom-Holman method (Wisdom & Holman 1991) (WH) in which the planetary pairs can have their own independently varying timestep, increasing the WH efficiency.When adapting the timesteps, a global timestep is simply redone if it is determined to not be reversible; the repeated timesteps are often a small fraction of all steps.AG is WH in which its global timestep can be adapted at block-synchronized times, without losing any of its beneficial longterm error properties.Indeed, AG can simply replace WH in most circumstances to increase its efficiency by making it more flexible and adaptable.These methods combined can be considered a simple, flexible, and adaptive version of the multiple timescale algorithm SYMBA, presented in DLL98.We implemented a symplectic algorithm based on SYMBA which we call MTS (Multiple Timestep Symplectic).We tested the symplectic and reversible algorithms on various challenging planetary N-body problems exhibiting close encounters.The accuracy of the reversible algorithms was as good as MTS and SYMBA (as seen in the figures of DLL98), and AG was up to 2.6 times faster than MTS.In contrast to SYMBA, the reversible methods presented here require no switching functions and are able to adapt timestep levels based on not only the pairwise separations between bodies.
To construct the reversible methods, we generalized the results of Hernandez & Dehnen (2023), who show how to switch reversibly between two timesteps.Our new algorithms instead can switch between an arbitrary number of maps, each with different pairwise planetary timesteps or global timesteps.The methods simply change timesteps when needed and then we check if the global timestep was reversible.If not (often in a small fraction of cases), the global timestep is redone, resulting in an "almost" reversible method (Hernandez & Dehnen 2023).For maximum efficiency, MTR and Listing 4 (not explored in this work) can be combined, so that steps are adapted locally in time, and we hope to present this implementation in future work.The codes studied in this paper are available for reference on Github 1 , and are written in an interpreted language.We leave a user friendly implemented in a compiled language for future work.
Our tests were limited to challenging systems with at most five bodies, including hierarchical binaries and a violent outer Solar System, but we expect the algorithms to work well in large particle systems, which can be studied in a language like C and is, again, left for future work.We remark that we have assumed conservative Hamiltonian systems when deriving our methods, but we can apply the methods here to conservative systems with nonconservative perturbations.
Taking a broader view, we see that traditional symplectic methods can be made simpler and more flexible simply by translating them to a reversible framework.While we have limited ourselves to planetary dynamics problems, there is no reason not to believe we can similarly transform other symplectic methods in other domains of astronomy.
We also remark on a significant finding in the course of this work, even if it is not the main result.Adapting the global timestep of any time-symmetric (including symplectic) algorithm via the simple prescription here, which builds on Hernandez & Dehnen (2023), can be immediately applied to many codes.However, it is curious that we found that, while the timestep can be arbitrarily decreased, it can only be increased at block synchronized times, or similar prescriptions.This should not be viewed as a limitation of AG because it should not really affect performance, but rather a feature of algorithms attempting to discretely vary the stepsize.Without taking care of this, the long-term error performance will be destroyed.This finding contrasts with the results of Hut et al. (1995), who vary timesteps reversibly and continuously but find no similar requirement.We note the Hut et al. (1995) algorithm is implicit, rendering it more expensive to use.

Figure 1 .
Figure1.A first test of the near explicit time-reversible methods MTR and AG.A Kepler orbit with e = 0.9 is integrated for 1000 periods P, and the energy error is recorded over time.We compare results to the SYMBA-like MTS method.The errors are all stable over time and the reversible methods are the fastest algorithms.

Figure 2 .
Figure2.The same test as in Fig.1, except now with e = 0.999.The errors over time for all methods are again comparable, and the fastest method is AG.

Figure 4
Figure 4. Error over time of a time integration of a violent outer Solar System, consisting of the Sun and outer giant planets, whose masses have been scaled up by a factor 50.We compare the MTR method ("Reversible"), in which 8 steps have been repeated (a fraction 8 × 10 −5 ), to the same method without any repetition ("Naive").While "Naive" reaches an error of order unity, MTR maintains an error to about 2 parts in a million from redoing the 8 steps.

Figure 5 .
Figure5.Pairwise distances over time for the system of Fig.4, obtained with the reversible method.Saturn's initial Hill radius is indicated by a dashed line.There is a close encounter between Jupiter and Saturn at 282 yrs.Then there are three additional close encounters between Uranus and Neptune.Saturn is ejected at ≈ 430 yrs.

Figure 7 .
Figure7.Attempted timestep levels, semi-major axes, and eccentricities for the pair of binary planets of the system of Fig.6.The results of MTR are used to make this pot which shows 1% of the time of Fig.6.The two pairs have timestep levels varying independently, which MTR handles well.