Terrestrial Exoplanet Simulator (TES): an error optimal planetary systems integrator that permits close encounters

We present TES, a new n-body integration code for the accurate and rapid propagation of planetary systems in the presence of close encounters. TES builds upon the classic Encke method and integrates only the perturbations to Keplerian trajectories to reduce both the error and runtime of simulations. Variable step size is used throughout to enable close encounters to be precisely handled. A suite of numerical improvements are presented that together makes TES optimal in terms of energy error. Lower runtimes are found in all test problems considered when compared to direct integration using ias15. TES is freely available.


INTRODUCTION
Understanding and predicting the motion of the celestial bodies has been an active field of research since the times of Newton and Kepler and is still equally as important today. Few analytical solutions exist for planetary motion and scholars have instead turned to numerical n-body techniques. N-body problems range from the most general form found in the study of plasma and star cluster dynamics (Aarseth 1999) through to systems with more inherent structure such as protoplanetary disks (Kokubo & Ida 1996;Kokubo et al. 1998) or exoplanet systems (Smith & Lissauer 2009). While integrators exist for the general case, by restricting one's self to systems that exhibit more structure one can leverage it to develop more efficient integration algorithms. In this paper we restrict ourselves to planetary systems whereby there exists a dominant central mass with any number of orbiting bodies. In this case, three key problems to address are: (i) Ensuring that solutions obtained remain accurate over the timescales required, in solar system formation and stability studies typically 10 9 dynamical periods. (ii) Ensuring that simulations can be completed within the available computing time. (iii) Ensuring that integrators can precisely model close encounters between objects.
IEEE 754 double precision floating-point numbers (IEEE 2019), typically used in scientific computing, have limited precision of about 16 decimal digits and as such errors are introduced when arithmetic operations are performed on them. All classes of numerical integrators are afflicted by this finite precision and over many time steps these round-off errors can accumulate to corrupt ★ E-mail: P.Bartram@soton.ac.uk the accuracy of solutions obtained. Furthermore, these are not the only source of error present as the choice of integrator and step size used can also allow for truncation error to be introduced, while a less than careful implementation of intrinsic functions, such as the trigonometric functions, can lead to the introduction of bias error. All of these error sources must be carefully considered when designing a numerical scheme to ensure solutions remain accurate for 10 9 dynamical periods. When round-off error can be made to dominate it has been found that specific numerical techniques can be used to ensure that the distribution of errors are symmetrical. This has led to the creation of integration schemes that are optimal in the sense that they follow Brouwer's law (Brouwer 1937) and exhibit a growth in relative energy error over integrator time, , proportional to √ . Rein & Spiegel (2015), with ias15, and Grazier et al. (2005) have both applied special numerical treatments to traditional integrators, i.e. non-symplectic integrators, used in the solution of general ordinary differential equations (ODEs) that result in integrators that follow Brouwer's law and are therefore suitable for long-term integrations. Despite being optimal in the sense of Brouwer's law, these schemes are computationally expensive and require many, typically upwards of a thousand, evaluations of the force function per orbit.
There are less computationally intensive means of ensuring invariants are preserved in long-term celestial mechanics integrations. For example, symplectic methods (Forest & Ruth 1990;Kinoshita et al. 1990;Saha & Tremaine 1992) can be used to place an upper bound on the truncation error of integrations by solving a system governed by a Hamiltonian that is slightly perturbed from that of reality. The use of symplectic methods ensures that the Poincaré invariants are conserved which in turn has favourable properties for energy and angular momentum conservation. Wisdom & Holman (1991) created a symplectic mapping, known as the Wisdom-Holman (WH) map, that revolutionised the field by reducing computational times by an order of magnitude as compared to traditional integrators. The WH map has, and continues to be, the workhorse of the field and forms the basis of many diverse propagation tools (Duncan et al. 1998;Chambers 1999;Grimm & Stadel 2014;Rein & Tamayo 2015). In planetary systems, the WH map exploits the dominant contribution to the dynamics by the star and uses the fact that secondaries, i.e. bodies in orbit around a more massive primary body such as the Sun, move on perturbed Keplerian trajectories to split the system Hamiltonian into separate Keplerian and perturbation terms that can be solved independently within a time step. This splitting allows for the WH map to make only twenty evaluations of the force per orbit for typical applications. Obtaining the solution for the Keplerian term analytically requires the solution of Kepler's equation (Battin 1987). There are many choices for solving Kepler's equation but universal variables (Battin 1987) are favoured for their versatility. Stumpff functions are typically used (Danby 1992) although other recent works shows that unbiased results can also be obtained without them (Wisdom & Hernandez 2015).
One drawback of the WH mapping, and symplectic schemes in general, is that they must use a fixed timestep to ensure that symplecticity remains unbroken. Whilst not a problem if bodies remain well separated it means integrators are unable to handle close encounters which are typically defined as encounters between bodies within one Hill radius . The ability to handle close encounters is highly important in celestial mechanics for modelling many problems. This includes: the threat of asteroids to the Earth (Giorgini et al. 2008), the behaviour of exoplanet systems after an instability event (Rice et al. 2018; Bartram et al. 2021), and the planet formation process itself (Davies et al. 2014). Options exist that enable invariants to be conserved during close encounters when using the WH map (Duncan et al. 1998;Chambers 1999;Rein et al. 2019) but they fail to obtain the true trajectories of bodies during the encounter for the typical step sizes chosen. Therefore, to study realistic trajectories of bodies during close encounters traditional integrators such as Bulirsch-Stoer (Bulirsch & Stoer 1966) or Everhart's RADAU (Everhart 1985) scheme are typically used.
In this work, we introduce our method, called the Terrestrial Exoplanet Simulator (TES), that aims to combine the accuracy and performance benefits of the analytical solution of the Keplerian motion, as found in symplectic schemes, with the flexibility of traditional integration schemes. We build upon the classic scheme of Encke, see e.g. (Wiesel 2010), to create a perturbation method that can be integrated with a traditional integrator. Importantly, we show that in this framework it is possible for close encounters to be handled precisely and with a reduction in computational cost as compared to performing a direct integration using the full nbody equation of motion. We show that through careful handling of round-off error through compensated summation (Kahan 1965;Higham 1993) the Encke method can be made optimal in the sense that it follows Brouwer's law with a relative energy error comparable to that of ias15. We offer two implementations for TES, the standard configuration where the implementation is purely in double precision floating point arithmetic and the extended configuration where 80 bit extended precision, i.e long doubles, are used in the Kepler solver. The latter implementation enables a further reduction in error of an order of magnitude for the same computational cost when compared to using ias15. Both TES implementations are capable of accurately handling close encounters and TES has already been used to study exoplanet evolution in the presence of collisions (Bartram et al. 2021).
During the preparation of this manuscript, a similar method has been published by Hernandez & Holman (2020). Both methods have been developed independently and while they make use of similar δq q q initial positions final positions reference orbits true orbits Figure 1. A three-body Encke method. For the inner planet, the position on the reference orbit, q, is shown along with the perturbation from it, q. The position on the true orbit,q, is also shown. Deviations from the reference orbits are greatly exaggerated for clarity. concepts, the end result differs, with the most significant difference being that our method can be used to model close encounters.
We begin in Section 2 with a description of the components making up TES. Section 2.2 contains the derivation of the equations of motion used. The solution of the analytical part of our model is described in Section 2.3 and the numerical part in Section 2.4. Section 3 contains details about specific numerical techniques implemented to ensure TES follows Brouwer's law. We show the impact of each of these techniques in Section 4. Beginning in Section 5, the second half of this paper contains a series of numerical experiments. In particular, Section 5.3 contains long-term integrations and Section 5.4 shows the performance in the presence of close encounters. We offer some concluding remarks in Section 6.

TES MODEL
In this section we begin by giving an overview of the method and then describe in detail the mathematical model, coordinate system, and force function used in TES.

General Encke Method
In the Encke method, the position of a given bodyq is made up of two terms: (i) q, the two-body reference trajectory, (ii) q, the perturbation to this two-body motion, such thatq = q + q. Figure 1 illustrates this concept, where the reference trajectory, q , can be obtained analytically at any future time by solving Kepler's equation and applying the so called and functions (Battin 1987, p. 156). In contrast, the perturbation term, q, must be obtained through numerical integration. The advantages of the Encke method over, e.g., a full n-body integration stem from the fact that the ratio of | q | /|q|, henceforth called the delta ratio, is much smaller than unity meaning that a given absolute precision in q can be obtained with lower relative precision in q which means less computation by the numerical integrator. In order to keep the delta ratio small one occasionally needs to update the reference trajectory to the current true orbit, a process called rectification. With a single rectification per orbit a delta ratio of 10 −2 can reasonably be expected for simulations of the outer solar system whereas for the inner solar system this ratio remains below 10 −4 . With such delta ratios maintained, the analytical solution of the reference trajectory becomes the primary source of numerical error. In order for this method to work, it is crucial that this propagation is precise, in the case of this work this means down to machine precision.
We continue by introducing a general form for creating an Encke method through two arbitrary governing Hamiltonians which therefore requires that we are integrating a conservative system without any dissipative effects present. Later, in Section 2.2, we use this form to derive our method for two specific Hamiltonians in our chosen coordinate system. We choose to work in canonical coordinates, and throughout this work the true state vector is denoted byẑ whereq andp are conjugate position and momentum vectors. Similarly, the reference orbit state vector is where q and p represent the position and momentum components of the reference orbits. Bothẑ and z are assumed to be in the same coordinate system. Therefore, the Encke method iŝ z = z + z and the perturbation term, henceforth called simply the deltas, for which we need to derive the equations of motion, is (1) The time evolution ofẑ and z for any conservative system is governed by their respective HamiltoniansĤ (ẑ) and H (z). Making use of the canonical structure matrix where is the identity matrix of appropriate dimensions for a given problem, one can apply Hamilton's equations to yield the equations of motion as Taking the time derivative of Eq. (1) and replacing appropriate terms with those from Eq. (2) gives a formula for finding the equations of motion for the deltas themselves as

Encke Method: Democratic Heliocentric (ENCODE)
In Cartesian coordinates, we define the state vector for this system asẐ whereQ ( ) is the generalised position vector andP ( ) is the momentum vector conjugate to it. If the bodies are only interacting through mutual Newtonian gravity then the time evolution of the system in Cartesian coordinates is governed by the gravitational n-body Hamiltonian (Leimkuhler & Reich 2005) where the subscript refers to the i ℎ body in the system. Throughout this work = 0 is reserved to refer to the central body.
There are no obvious Hamiltonian splittings of Eq. (4) that allow for the dominant Keplerian motion of the secondary bodies about the primary, due to the dominant central mass, to be isolated from the general evolution of the system. For this it is necessary to introduce a different coordinate system. There are several coordinate systems available and Hernandez & Dehnen (2017) give a good overview of the canonical coordinate systems for the n-body problem. The Jacobi coordinate system was used by Roy et al. (1988) to create an Encke method to simulate the outer planets of our solar system. In this coordinate system, each body has a reference orbit taken with respect to a different moving centre of mass that depends on the position and mass of all other bodies who's orbits are smaller than its own. Therefore, as the relative size of the orbits of bodies changes in a system it becomes necessary to recalculate reference trajectories with respect to a new moving centre of mass. To avoid this complication and the numerical error it could introduce we instead opt to use democratic heliocentric (DH) coordinates (Duncan et al. 1998) which are common in celestial mechanics (Chambers 1999;Grimm & Stadel 2014). In DH coordinates the equations of motion are such that the position of each body is expressed relative to the central body, which we denote with the index zero throughout this work. The momentum of each body, however, is expressed relative to the barycentre of the system. The coordinate change from Cartesian to democratic heliocentric coordinates is given bŷ where M is the total system mass, i.e. =

=0
. Therefore,q 0 and andp 0 are the centre of mass and momentum of the system, respectively. The Hamiltonian,Ĥ (ẑ), as a function of the previously defined state vector,ẑ, in these new coordinates iŝ wherê Each of the three components ofĤ are physically identifiable: (i)Ĥ is the motion of the central body, (ii)Ĥ is the Keplerian motion of the secondary bodies about the central body, (iii)Ĥ is the gravitational interactions between secondary bodies.
An additional fourth term that only depends uponp 0 and represents the motion of the centre of mass has been excluded under the assumption of a stationary barycentre, and it is therefore unnecessary to propagateq 0 andp 0 .
Note that Equation (4) can also be used to obtain a Hamiltonian describing a system of particles that only interact with the central body by enforcing that = 0 thereby removing the gravitational interactions between secondaries. After applying the same coordinate transformation from Eqs. (5) and (6) the reference orbit Hamiltonian, H (z), as a function of the previously defined state vector, z, in democratic heliocentric coordinates reads InsertingĤ (ẑ) and H (z) into the general Encke form in Eq.
(3) yields the equations of motion for the deltas in democratic heliocentric coordinates as Here, and throughout this work, dots are used to signify the time derivative of a variable. For reasons discussed later in Section 2.4 we take an additional time derivative of q to obtain q. There are several key features to note about these equations. The summation in the trailing term in Eq. (10) starts at an index of 1 meaning that it only captures interactions between secondary bodies; the gravitational interaction with the central body, with an index of 0, is captured by the leading term instead. In this term, it can be seen that there is cancellation between similar terms that depends upon only the size of q . Therefore, so long as this term remains small then p will generally also remain small in comparison to the the reference trajectory derivative, p. The exception to this are close encounters between secondaries, where this term can grow significantly. Again, it is this reduction in the relative size of the term to be integrated that leads to the performance increases from an Encke method.
Typically, the Encke method as used in astrodynamics assumes massless secondaries, and, as such, that the centre of mass of the central body is coincident with the barycentre, i.e., the location to which reference orbits are taken is fixed. This is not the case in celestial mechanics, however, due to the relatively high mass ratio of planets in comparison to their host stars, a ratio of approximately 10 −3 for our solar system. As an example, consider a two-body problem consisting of Jupiter and the Sun. In this case, the elliptical motion of the Sun about the Sun-Jupiter barycentre is captured by the trailing term in Eq. (8) as the negative momentum of the star This shows that even in a purely two-body case there will still be a deviation from the reference trajectory around the Sun over time, and the magnitude of this deviation is inversely proportional to the ratio of the mass of the star 0 to the mass of the other bodies within the system ≡

=1
. We call this ratio the system mass ratio and define it as 0 . Practically, this means that systems with a small system mass ratio have the greatest potential performance gains. Alternatively, in systems with many bodies an axially symmetric distribution of mass reduces the motion of the central body through a cancellation of terms in . Clearly, a distribution such as this is non-physical for planetary systems, but typical accretion disks in planetary formation are symmetric enough to see cancellation.
We point out that for the previously discussed reasons it is no longer necessary to evolve the motion of the central body as an separate set of equations thereby reducing by 1. However, whereas a pure n-body integration can take advantage of a second order formulation of the equations of motion to reduce the number of equations to be integrated, especially in the case where the motion is not velocity dependent, this is no longer the case for these equations where two first order ODEs, Eqs. (8) and (10), must instead be integrated. Note that this does not double the computational cost. The additional cost in the RHS is small and scales only as O ( ) and therefore the dominant O ( 2 ) interaction term remains unchanged. Furthermore, the integration procedure itself performs identical operations for each first order ODE meaning that vectorisation, either manual or automatically performed by the compiler, can recover some of the additional computational cost. Finally, this does not have a noticeable performance penalty when calculating the reference trajectories.

Analytical Solution
In this section we describe how we solve the two-body problem with universal variables making use of the and functions (Danby 1992). If we assume that the integration of the perturbation terms can be performed in a way that ensures the truncation error is below floating point precision for z then the overall precision of the scheme depends upon the precision in the solution of the Keplerian motion. As such, we require a highly accurate solver implementation that is also non-biased to ensure that the energy growth follows Brouwer's law for long duration integrations. Several modern implementations of Kepler solvers exist that are suitably accurate, e.g. (Wisdom & Hernandez 2015;Rein & Tamayo 2015). We have chosen to follow the implementation presented in WHFAST (Rein & Tamayo 2015) but with some additional numerical improvements specific to our needs, these are discussed in Section 3.2.
We wish to solve − 1 independent two-body problems for which the time evolution is governed byĤ in Eq. (7). Typically, the f and g functions use a reduced mass parameter such that = ( 0 + ); however, in democratic heliocentric coordinates, used here, the reduced mass is given by = 0 (Grimm & Stadel 2014). In the following discussion we focus on a single two-body problem from the vectors q and p above. For the remainder of this section q and p refer to the position and momentum of a single body undergoing Keplerian motion. As a result, we drop the index referring to each body: an index of zero now refers to values at the start of an integration step. After using the WHFAST algorithm to solve for the universal anomaly and obtain the -functions (Rein & Tamayo 2015), the and functions used are where d is the time step. This allows for the solution to the Kepler problem, after a timestep of d , to be obtained from the initial values of position and momenta and a linear combination of them to be applied to each as an update term, Δq and Δp, such that When formulated in this manner, the smaller d , relative to the orbital period, the smaller the size of the bracketed terms and therefore summing them first is more numerically robust. Additionally, TES uses a value of d approximately 1 /300 th of an orbit, roughly an order of magnitude smaller than, e.g., WHFAST and this means that the relative sizes of Δq to q 0 and Δp to p 0 enables compensated summation to be used to further reduce round-off error; our algorithm for this is described in Section 3.4.

Numerical Solution
While there are many numerical integrators to choose from, a particularly efficient and accurate choice for astrophysical problems is the RADAU scheme of Everhart (1974). A useful feature of RADAU is that once a polynomial is fitted to the force, in the manner discussed below, it is possible to integrate it analytically one or multiple times allowing for both, e.g., velocity and position to be obtained from the same polynomial. Everhart (1985) found that his scheme is best suited for directly integrating second order ODEs, e.g. = ( , ), without reduction to a pair of first order equations. In fact, for the same number of evaluations of the force function he found the solution to second order ODEs could be as much 10 6 times more precise than if they were reduced to first order and integrated. In the derivation of our equations of motion we have had to reduce the system to a pair of first order equations, and integrating q and p in this manner does in fact cause a reduction in performance. We observed that integrating q and p is enough to circumvent this problem and opt to integrate these equations in TES. As discussed in Section 2.2, the effect of this choice on the computational cost is minor. Rein & Spiegel (2015) further refined the RADAU method with two new algorithms, one for convergence criteria and one for stepsize control, and additionally included compensated summation to ensure a symmetrical distribution of round-off errors. Their new implementation is 15 th order and is called ias15. This is the version we have chosen to base our implementation on to integrate q and p. Next, we introduce just enough of the ias15 algorithm to discuss our modifications to it. Only the process of integrating q is discussed but a similar process is also followed for p.
We need to simultaneously solve 3 equations of the form = (q, p, q, p, ), one for each directional component of bodies. In order to do this, RADAU expands the acceleration, , in time, , into a truncated series such that where ℎ = /d , d is the size of an integration step, and 0 and 0 are the time and acceleration at the start of an integration step. The coefficients are fitted though an iterative predictor-corrector process that performs an evaluation of at the start of an integration step, = 0 , and at seven substeps within the integration step. The sampling locations in time for each substep, where = 1...7, are chosen in accordance with Radau quadrature spacings (Radau 1880) to maximise the order of the scheme. Once the coefficients, , are obtained to a sufficient precision then Eq. (14) can be integrated analytically to obtain an estimate of q at the end of a step, 1 = 0 + d , as .
At the start of a step an analytical continuation of the curve fitted to via the accurate values calculated during the previous step are used to generate a predictor in the form of the values of to use in the current step. To ensure convergence when iterating to obtain the coefficients a convergence criterion is required. We opt to use a convergence criterion similar to that in ias15 and therefore monitor the change in the final coefficient in our truncated series, i.e 6 , from one iteration to the next, we call this change Δb 6 which is a vector containing coefficients for all 3 equations. We then compare the maximum change to the maximum magnitude of the reference orbit acceleration, q, and we terminate when We find that this criterion performs better in our use case than if q is used in the place of q in Eq. (15) which is more typical. The typical criterion would ensure that the change in the coefficients in the series expansion of the acceleration of the deltas are precise to floating point precision. By design, in TES, the deltas are much smaller than the reference orbit terms and therefore it is unnecessary to converge this far to achieve a combined relative tolerance inq+ q of 10 −16 for use within an integration step, e.g. in the predictors.
In contrast to the convergence criterion, as we wish to suppress truncation error across long duration integrations, it is a requirement that the step size, d , is chosen such that it controls the truncation error in the series expansion, Eq. (14), itself. We monitor the truncation error for all 3 equations which is estimated by . To do this, we monitor the absolute value of the 6 coefficients for all 3 equations, which we call b 6 , relative to the acceleration of the deltas at the start of an integration step to obtain an estimate of the smoothness of the acceleration of the deltas over a given step as The value of is then used to determine the next step size, d +1 , as where tol is a dimensionless tolerance parameter. We offer no analytical reasoning for choosing a value of tol; however, numerical experiments demonstrating its effect are shown in Section 5. We find that a default value of 10 −6 is suitable for maintaining Brouwer's law for 10 9 dynamical periods for simulations of the inner solar system and for handling close encounters. We have not included a step rejection algorithm as we found little benefit.

Rectification
Rectification is the name given to the process whereby a new reference trajectory is taken by adding the current reference trajectory together with the deltas. Our algorithm for doing this can be found in Section 3.4. Rectification is important and the frequency with which is it performed has two conflicting effects on the efficiency of the scheme that must be balanced. Firstly, because rectification 1 10 3 10 16 10 19 10 32 10 35 precision range of state variables q* q q* q causes the deltas to be set to zero the analytical continuation used to obtain a prediction of the values of becomes much less precise in the integrator during the subsequent step and therefore increases the number of iterations required for convergence, typically by one. Secondly, in contrast, rectifying causes the size of the deltas to be reduced to zero and therefore reduces the computational cost on the integrator in many subsequent steps. We experimented with several rectification schemes based upon the size of the deltas relative to the reference trajectories but found the optimal cutoff value depends largely upon the mass ratio of the system. Ultimately, we found that a scheme where a rectification for all particles is performed between once and twice per orbit of the shortest period body is simple to implement and provides a good balance between the two effects. We choose to rectify according to the golden ratio and perform 1.618... rectifications per orbit of the body with the smallest period in the system to help avoid any possible bias effects due to resonances. We combine this with a fall back method whereby we also rectify if the delta ratio exceeds a given value, typically 10 −3 .

IMPLEMENTATION DETAILS
In this section we introduce a collection of numerical implementation features that improve the overall energy conservation of TES for long integration times and enable it to handle close encounters.
A key component used through this work is compensated summation (Kahan 1965). Compensated summation allows for the error made during the summation of two double precision floating point numbers to be obtained and kept as an additional double precision number, known as the compensation variable. The error made is then subtracted from any future additions. This is applicable, e.g., when performing the final update of the state vector in an integration. We use this technique to ensure a symmetrical distribution of round-off error and minimise the total energy error in simulations. Compensated summation requires an additional compensation variable be kept for each variable being compensated to keep track of lost precision. To this end, we define a composite datatype, written as { , * } and comprised of two components: the variable itself, and its compensation variable which is denoted with a star; each component is stored in the computer as a double precision floatingpoint number. In infinite precision, the true value represented by a composite datatype Σ = − * where the minus sign is due to the Kahan summation algorithm. Figure 2 shows the relative ranges in magnitude that are covered by the variables used in TES. Compensated addition and subtraction operations are denoted by ⊕ and Here, a compensation variable is used to keep track of lost precision across an entire integration step, as shown by label 2. Finally, at the end of the integration step, 1 , label 3, compensated summation is used to combine the separate compensation terms. Compensated summation is also used to reduce error during rectification but this is not shown here.
respectively. We find it advantageous to use compensated summation in three ways in addition to internally within the integrator: (i) Propagating the reference trajectories.
(ii) Combining the reference trajectories and deltas.
(iii) Ensuring precision is maintained across rectifications.

Encke Method: Democratic Heliocentric (ENCODE)
Equation (10) contains a subtraction between two terms of similar size and, owing to the finite precision of floating point numbers, this causes cancellation of significant digits which leads to a large decrease in the relative precision of p. The loss of relative precision here depends on how small the difference between the terms is. In particular, after a rectification the difference can be very small, e.g. 10 −12 has been observed. In addition to the risk of introducing numerical error when calculating the acceleration this also poses a problem for the step size control algorithm in Eq. (16) and (17). The algorithm works by ensuring that the step size is chosen such that the acceleration approximated by the expansion in Eq. (14) is smooth to a precision of tol. Oftentimes, the numerical cancellation in Eq. (10) means that the required degree of smoothness cannot be met. This can lead to a situation where the step size will shrink uncontrollably as the algorithm tries to shrink the step size further to reach an unattainable smoothness. To remedy this situation, it is possible to reformulate the problematic term to avoid the subtraction of like terms and rewrite the term p as (Battin 1987, p. 449) where p is now obtained without loss of significance. This is the equation that we have implemented and it decreases numerical error as well as preventing a step size lockup from occurring.

Analytical Solution
Due to the relative size of the reference trajectories q compared to the deltas q, the dominant contribution to error growth stems from roundoff errors in q. To minimise this potential we have used compensated summation in the final update step of the and functions, Eq. (12) and (13). With compensated summation they become Due to the relatively small update terms, Δq and Δp, this allows for the value of q and p to be maintained to machine precision across Kepler solver steps. Owing to non-linearity when calculating thefunctions it is more precise to take steps of size than one step of size . We therefore go further and also minimise the value of d in Eq. (11) (the simulation time passed since the and function basis vectors were calculated) which in turn decreases the size of Δq and Δp. To do this, we only ever take a single step in the Kepler solver before we calculate new basis vectors, i.e., we calculate new basis vectors at the start of each step as well as at each substep required by the integrator. The locations in time that we perform both the compensation in Eq. (18) and a recalculation of basis vectors are illustrated in Fig. 3 and are marked by the label 1. Here, it is shown that the universal variables compensation is used at the start of a step, 0 , at the end of a step 1 , and also at all substeps required by the integrator, . Whilst the standard TES configuration only makes use of double precision floating-point arithmetic throughout, there is also a build configuration that allows for the selected use of extended precision arithmetic through the C long double datatype. We only use extended precision to perform the entirety of the reference trajectory calculations, and we achieve an improvement of energy conservation of an order of magnitude for long duration integrations. Using extended precision sparingly like this allows the compiler to optimise the majority of the code to make use of single instruction multiple dispatch (SIMD) operations which are not generally available for extended precision variables on Intel or AMD64 hardware. Figure 2 shows an example of the relative scales of the various state variables used within TES: the state vector variables are in blue and the compensation variables are in orange. The cross-hatched area shows the extra precision required in the reference trajectory, q, such that the extra relative precision in q can be used to create a scheme with a round-off error, per step, of 10 −19 , assuming that q remains suitably small. The cross-hatched area is also exactly the extra precision that can be obtained through the use of long doubles in the calculation of the reference trajectories, and, as will be shown in Section 4, allows for a reduction in energy violation of over an order of magnitude.

Numerical Solution
In a similar fashion to Eq. (18), the numerical integrator also obtains a compensation variable at the end of a step for each variable being integrated, as is shown in the top panel of Fig. 3. In our case, this means that at the end of an integration step we obtain { q, q * } and { p, p * }. We therefore have two separate sets of compensation variables: one for the reference trajectories, q * and p * , and one for the deltas, q * and p * ; however, these two sets must be combined at the end of a step to ensure the correct error is used in the subsequent step. Combination of compensation variables is performed at the end of each step as can be seen by label 3 in Fig. 3. It is achieved through the following algorithm: * ← 0, where * is a temporary summation variable for each body. This algorithm begins by combining the two compensation variables q * and q * into a new compensation variable * . This is done as a compensated subtraction into the reference trajectory q in case the subtraction causes the range of the double precision variable * to overlap with q. After this, the range of * now overlaps with the least significant region of q and a third compensated subtraction is therefore used to update q accordingly. In this final subtraction the reference trajectory compensation variable, q * , is used to store the final summation error such that it can be used immediately in the subsequent Kepler step. The same process is also used for the momenta terms, p and p. Note that this process differs from rectification, described next, as only the reference trajectories are used here.

Rectification
Compensated summation can also be applied to maintain precision across a rectification by following a very similar algorithm: * ← 0, where again * is a temporary summation variable for each body. This algorithm begins by performing the rectification process by summing the reference trajectories and deltas together and using temporary summation variables to capture any lost precision from the addition. The reference trajectories, q, are then refined using the compensation variables q * and q * with any lost precision being captured again by * . Due to the sign convention used in the compensated summation this means that − * now contains the rectified value of the deltas which is then placed into q. This process is also used to rectify the momenta, p.

VALIDATION OF IMPLEMENTATION DETAILS
In the previous section, we described four key numerical features present in TES, in summary they are: (i) Kepler compensation described in Section 3.2.
(ii) Kepler long doubles described in Section 3.2.
(iii) Combining compensation variables described in Section 3.3. (iv) Rectification compensation described in Section 3.4.
With the exception of Kepler long doubles, the default configuration of TES enables all of these features. Figure 4 shows the effect of disabling compensation at various points in TES has upon the energy conservation in a simulation of the inner planets of our solar system over a period of 10 8 Mercury orbits. The data plotted is the RMS of twenty realisations of the initial conditions randomly perturbed on the order of 10 −15 . The modification to the force function in Section 3.1 to avoid numerical issues due to cancellation of similar size terms is enabled for all configurations as without it a step size lockup can occur. As a baseline for the performance without any numerical improvements a naive Encke implementation is included that has no numerical improvements other than the reformulation of the acceleration. The default configuration of TES using only double precision with all compensated summation schemes enabled is shown under TES default settings. Here, the conservation of energy for TES is two orders of magnitude better than for the naive Encke scheme. If the use of extended precision floating point variables is permitted in the Kepler solver then the energy conservation in TES can be further improved by an order of magnitude as compared to the default configuration. Disabling individual compensation schemes results in energy conservation at least ten times worse than that of TES with default settings. In the worst case, when compensated summation is not used in the final update step of the and functions, the conservation in energy can be as poor as just using the naive Encke method which highlights how important a precise solution to the dominant Keplerian motion is to schemes of this nature.
In order to ensure that TES can handle close encounters we ran a simulation of three Earth mass planets orbiting a solar mass star at 1 AU as per Bartram et al. (2021). Planets are tightly packed and over time the system becomes unstable causing close encounters between them. Figure 5 shows one of these encounters where the planets passed closer to one another than the Moon is to the Earth. Firstly, in the left hand panel, a close encounter after approximately five orbits leads to a much closer encounter an orbit later. Next, in the central panel the step size controller shrinks and subsequently 10 −14 with 20 steps per orbit expands the step size appropriately. Finally, in the right hand panel, the relative energy error can be seen to be maintained below 10 −15 throughout. Interestingly, the relative energy error values in this plot take six distinct values indicating that the error is confined to the three least significant bits of our double precision representation of the energy. These examples validate that the TES model derived in Section 2 and implemented as described in Section 3 is indeed capable of performing highly accurate long term integration as well as handle close encounters between bodies effectively.

NUMERICAL EXPERIMENTS
To further investigate the performance of TES in a variety of settings, in this section we perform a series of numerical experiments. We also provide comparisons against a number of other integrators. The schemes used are: TES (double); TES (long double), which makes use of long doubles in the Kepler solver; naive Encke; ias15 (Rein & Spiegel 2015) from the Rebound package (Rein & Liu 2012); and Bulirsch-Stoer (Bulirsch & Stoer 1966) as well as the hybrid (Chambers 1999) integrators from the MERCURY package. Table 1 contains the default tolerances used throughout these experiments unless otherwise specified. In the case of the TES, naive Encke and ias15 schemes the tolerances used are the recommended defaults. All runtime measurements were performed on an Intel Core i7-6700 CPU running at 3.4 GHz.

Efficiency Mass Dependence
As discussed previously, the magnitude of the deltas in comparison to the magnitude of the reference trajectory, the delta ratio, must be kept small in order to maximize the efficiency of an Encke method. Simultaneously, one must not rectify too frequently as rectifications degrade the precision of the predictor in the subsequent step and thus incur a performance penalty. In the absence of close encounters, the dominant contribution to the acceleration of the deltas is related to the motion of the central body which in turn is dependent on the system mass ratio. Therefore, this experiment is designed to understand in which regions of system mass ratio TES is effective.
We perform integrations of twenty-one two-body problems for our full selection of integration packages. We have opted to examine the range of system mass ratios that can be found in our own solar system if each planet is taken in isolation with the Sun. Therefore, this experiment ranges from a secondary of Jupiter mass, , down to a mass of 10 −4 , approximately equal to that of Mercury. The samples across the range of masses are logarithmically spaced, and the primary is always a solar mass star. The secondary body is placed on a circular, co-planar orbit at 1 AU and integrations are performed The primary is a solar mass star and the mass of the secondary is varied across a range coincident with that of our solar system. The Encke based methods must account for the motion of the central body and the two-body problem is therefore still an appropriate test case. Each data point is the average of twenty identical integrations.
for 10 4 orbits. Runtime is calculated as the mean of twenty identical integrations for each two-body problem for each integrator. Figure 6 shows the relative energy error achieved in these experiments. Here, with the exception of the hybrid scheme, we find no dependence between the relative energy error and the system mass ratio, simply meaning that the error control algorithms within each integrator are performing as expected. However, we do find large differences in the precision of the various integration schemes. In particular, the Bulirsch-Stoer, hybrid and naive Encke schemes all fail to reach the regions of highest energy conservation. In contrast, the schemes that are floating point arithmetic aware, i.e TES and ias15, are much more precise. TES (double) and ias15 have almost identical performance across the entire range of system mass ratios examined. TES (long double) is the best performer overall and outperforms TES (double) and ias15 by up to two orders of magnitude. Figure 7 shows the runtime for the same experiments where a cluster of curves can be seen as well as the hybrid scheme which is between six and eight times faster than the non-symplectic schemes. We find that the standard deviation in runtime across all TES realisations is 146 ms. The Bulirsch-Stoer, hybrid and ias15 schemes can be seen to exhibit no dependency of the runtime on the system mass ratio. However, all of the Encke based schemes, i.e. TES and naive Encke, show a positive correlation between the runtime and the system mass ratio, as predicted. An interesting comparison is that of TES (double) and ias15 where for systems with a smaller mass ratio TES is able to achieve the same levels of precision with only 75% of the computational cost. TES (double) remains more efficient until the mass of the secondary is roughly 10 −1 . Therefore, for maximum benefit, TES should be applied to systems with a system mass ratio below this value. Finally, TES (long double) can also be seen to perform well with a runtime below ias15 despite having a better conservation in energy of up to two orders of magnitude.

Convergence and Runtime Comparisons
Next, we study the convergence and runtime of TES in comparison to the wider field of integrators over a period of ten thousand orbits of the innermost planet. In the previous section we showed that TES is most effective in systems with a low system mass ratio, and we have therefore chosen the inner planets of our solar system as a test problem using initial conditions taken from the NASA Horizons database (NASA 2021). The tolerance of each of the nonsymplectic integrators is varied over a range of values such that the relative energy error no longer converges. TES and ias15 use a range ending at the recommended operating tolerances in Table 1.
The hybrid scheme uses a fixed tolerance of 10 −14 throughout but the step size, which must be kept constant within an integration, is varied until the relative energy error no longer converges. Figure 8 shows the relative energy error against the number of steps per orbit. Once again we can see a divide between the optimal and non optimal schemes, TES and ias15 clearly conserve energy more precisely than the other schemes with TES (long double) being the most precise by roughly an order of magnitude once the round-off error dominated regime is entered at roughly fifteen steps per orbit. The power of the Encke method can be seen in two places in this plot. Firstly, in the truncation error dominated region, i.e. the region below roughly fifteen steps, where the relative change in energy for a given step size is approximately three orders of magnitude smaller than any of the direct integrations. Secondly, in the furthest right data point for TES and ias15 where the recommended default tolerances for the Encke based methods yield a reduction in the number of steps taken per orbit when compared to a direct integration. Figure 9 shows how these benefits manifest themselves in the runtime. Immediately, the hybrid scheme can be seen to stand apart from the others and is indeed much faster than any of the nonsymplectic integrators; however, as we will show in Section 5.4 the relatively low precision of the hybrid scheme is not entirely suitable to modelling exoplanet evolution in the presence of close encounters. The Bulirsch-Stoer scheme has a poor runtime in comparison to the other integrators and does not reach the highest levels of precision either. The naive Encke method has a favourable runtime in the truncation error dominated regime but is not capable of the energy conservation of the optimal, floating point implementation aware methods. Of the three remaining integrators, TES (double), TES (long double) and ias15, the performance is similar. However, for the recommended default tolerances, the furthest right data point for each integrator, TES (double) is the fastest and is approximately 20% faster than the slowest scheme. Interestingly, TES (long double) has the best energy conservation by up to an order of magnitude and has very comparable runtime to ias15 despite the disadvantage of non being able to use vectorisation in the Kepler solver.

Long-term Integrations of the Inner Solar System
In the field of exoplanet modelling it is typical for simulations to span a billion dynamical periods. It is therefore of great importance to ensure that propagation schemes follow the optimal error growth of Brouwer's Law, i.e. ∝ √ for the relative energy error, where is the number of steps taken. We have therefore chosen to perform long-term simulations of the inner solar system to ensure TES conforms to this requirement. We use the tolerance in Table 1 for TES and ias15. In keeping with the original ias15 experiments (Rein & Spiegel 2015) and to generate a statistical sample we perform twenty integrations for TES and ias15 with a perturbation in the initial conditions on the order of 10 −15 . The RMS of these twenty realisations is plotted in Fig. 10. Additionally, results of three integrations highlighting the performance of the Bulirsch-Stoer integrator are also shown for tolerances of 10 −13 , 10 −14 and 10 −15 . Integrations performed with TES (double) span the full 10 9 Mercury orbital periods whereas all other schemes are terminated after 10 8 Mercury orbital periods to save on computation. Finally, two slopes are included: one in grey marking the linear error growth typically associated with truncation dominated regimes, and another, in brown, showing the optimal error growth associated with the symmetrical distribution of round-off error required for Brouwer's Law. Figure 10 highlights that TES (double) and TES (long double) both follow Brouwer's law for the full integration duration and therefore show that these schemes are well suited to the long duration integrations required in exoplanet modelling. Note that despite the larger steps taken by TES (double) the use of the Encke method has enabled the truncation error growth to be suppressed for the entirety of integrations. When performed by an integrator that follows Brouwer's law, the RMS relative energy error of a suitably large number of realisations ≈ √ where is a constant approximately at floating point precision, i.e. ≈ 10 −16 , and is the number of steps taken. In double precision, TES performs marginally better than ias15 and we believe this is due to the larger step size taken by the Encke method reducing the √ term. TES (long double) performs only the solution of the Keplerian motion in extended precision, i.e. the integrator and force models use only double precision. However, even with this sparing use of extended precision TES is able to attain a relative energy error of over an order of magnitude better than if only double precision is used throughout, and importantly, it does this without excessive computational cost. Most notably, TES (long double) is able to integrate for up to 10 5 orbits before there is any noticeable growth in relative energy error above the floating point floor. Both TES (double) and ias15 already start to show error growth after just hundreds of orbits.

Apophis 2029 Encounter
The last example we present studies the performance of TES in the presence of close encounters. When the dynamics of planetary systems are such that bodies do not remain well separated indefinitely it is important that close encounters are handled accurately, e.g. when modelling the behaviour of systems after an instability event (Rice et al. 2018;Bartram et al. 2021). Closer to home, another example is the accurate modelling of near Earth objects (NEOs) to understand potential hazards to humanity. Of the known NEOs, one of particular interest is Apophis owing to the fact that it will pass within 100, 000 km of the Earth in 2029 bringing it closer than the Moon. The simplified model we use (Amato et al. 2017) evolves only the Sun, Earth, and Apophis and makes use of purely Newtonian dynamics throughout. While not a realistic model, it allows us to showcase the behaviour of different integrators in the presence of strongly perturbing close encounters. To obtain a set of initial con- ditions and a high precision reference solution trajectory, first the position of the three bodies were obtained from the NASA Horizons system at the point of closest approach during the 2029 encounter. Next, a Dormand-Prince integrator (Dormand & Prince 1986) operating in quadruple precision was used to integrate backwards for fifty years to obtain the initial conditions. Finally, the same integrator was used to propagate the initial conditions forward for one hundred years to obtain the final conditions. We consider this high precision integration the truth, such that the final positions can also be used as an error metric in addition to energy conservation. The trajectories followed by the three bodies are plotted in Fig. 11 and the transition of Apophis from the Apollo to the Atens group after its flyby with Earth can clearly be seen. Figure 12 shows the separation between Apophis and the Earth over the simulated period of 100 years. Here, we see that Apophis makes many close approaches within 10 Earth Hill radii, or 0.1 AU, but only one very close approach which is in 2029. This approach distance is approximately 2.5 × 10 −4 AU or roughly 17, 000 km, well within the geosynchronous orbital altitude at 35, 786 km. We perform integrations at 61 different tolerance settings for each of the integrators in Table 1. The tolerances used are chosen in the same way as in Section 5.2. At each tolerance we perform twenty identical integrations and take the mean value for the runtime. Figure 13 shows the relative energy error for a given runtime. First, consider the behaviour of TES (double) in this plot. We find that in the truncation dominated regime it is amongst the fastest schemes competing only with the naive Encke and TES (long double). TES (double) can also be seen to be fastest at the default tolerance (Table 1) which is represented by the furthest right blue data point; this is the recommended default setting for users. Once TES (double) has converged such that only round-off error is present it is clear that the numerical implementation described in Section 3 has improved the performance by between one and two orders of magnitude in comparison to a naive Encke method (red). Due to the very short integration timescale of only one hundred orbits TES (long double) shows no advantage over TES (double) and the runtime for the default tolerance setting is comparable to that of ias15. The MERCURY integrators, Bulirsch-Stoer and hybrid, fail to conserve energy as precisely as the schemes based upon Everhart's Radau and reach final values of relative energy error of 10 −13 and 10 −12 , respectively. Interestingly, in contrast to Fig. 9, we can see that in the presence of repeated close encounters the hybrid scheme runtime increases to be comparable to that of the non-symplectic integrators.
While the conservation of energy is a common metric for the accuracy, its value is usually dominated by the most massive objects in a given system. Given the low mass of Apophis in comparison to the Earth we also look at the final position of Apophis as output by each integrator in comparison to our quadruple precision numerical solution. Figure 14 shows the error in the position of Apophis, , for a given runtime. Again starting with TES (double), there are two regions of performance: one where the runtime is low, below 32 ms, and it performs poorly, and one where the runtime is higher, above 32 ms, where the scheme performs well. What is interesting is the lack of a transition period between these two regions which is present in, e.g., ias15; instead, all of the Encke based methods are either highly precise or highly inaccurate. Fortunately, the recommended default tolerance for both implementations of TES is well within the highly precise solution region. Given a suitable tolerance we can see that both implementations of TES and ias15 are comparable in precision, with TES being slightly faster for the default tolerances (Table 1).
For this problem, TES (long double) shows interesting behaviour once it has converged to roughly = 10 −7 AU. At this point, its performance becomes highly consistent regardless of decreases in tolerance. This is in contrast to all of the other Everhart based schemes which exhibit variance in the precision of the final position of Apophis with changing tolerance. We find that the precisions seen here are highly consistent with Amato et al. (2017) who find their implementation of Everhart's Radau scheme integrating using Cowell's formulation converging to roughly 10 −7 AU. In order to obtain more precise solutions authors had to employ regularisation techniques on the equations of motion.
The positional precision obtained after a close encounter is highly important if one wishes to model repeated encounters. Deviations from a true trajectory are greatly amplified during a close encounter and any inaccuracies in modelling the first encounter in a series are therefore increased in all subsequent encounters. We find that the precision achieved by TES in these experiments equates to a positional error at the time of closest approach of approximately 10 cm whereas the naive Encke method achieves an error of 10 m. Given that the size of keyholes, i.e. regions of space on the b-plane formed by the separation vector at the point of closest approach, found by Farnocchia et al. (2013) that can lead to a resonant return trajectory for Apophis are between 6 cm and 600 m, the precision in the naive Encke makes it unsuitable for use in this challenging application domain. Therefore, the numerical improvements in Section 3 that comprise TES are sufficient to improve a naive Encke method to the point where it can now be used in the study of NEO asteroid dynamics.
We find that the hybrid scheme is not precise enough at the step sizes presented to accurately model the trajectory of Apophis during its flyby with Earth. However, we mirror the findings of Amato et al. (2017) and find that if the hybrid scheme step size is reduced by roughly a factor of one hundred then it is possible to obtain positional errors as low as 10 −5 AU although the computational cost for doing this means this option is of little practical importance.

CONCLUSIONS
We introduced TES, a new integrator for planetary systems that follows Brouwer's law and permits close encounters between massive bodies. TES builds upon the classical Encke method and takes advantage of the dominant nature of the star in planetary systems. We show that TES is effective across a wide range of planet to star mass ratios but find that the more dominant the central body the more effective the scheme is, with excellent improvements in speed being seen in simulations of the inner solar system.
In Section 2.2 we derived a new version of the Encke method in democratic heliocentric coordinates (ENCODE) and presented the equations of motion in this coordinate system. Additionally, we implemented a series of numerical improvements that reduced the round-off error by two orders of magnitude as compared to a naive Encke method. TES is optimal in that it follows Brouwer's law and has an RMS energy error slightly below that of ias15.
We performed extensive comparisons against ias15 in Rebound, and the Bulirsch-Stoer and symplectic hybrid schemes within the MERCURY package. We found that for well separated systems TES is the fastest non-symplectic scheme for a given precision. In the presence of close encounters, we found that TES is able to reach precision much greater than either of the MERCURY schemes, in terms of both conservation of energy and final position, with a precision comparable to ias15. TES is open source and accessible at https://github.com/ PeterBartram/TES. It is available in a "double" version using only double precision floating point arithmetic and a "long" version using extended, 80 bit, floating point arithmetic in the Kepler solver. In double precision, we found that TES is only 15% faster than the extended precision implementation which is two orders of magnitude more precise, but for portability we recommend the double precision version as the default. Regardless of the version used, we found that TES is faster than ias15 for the same energy conservation in the problems examined although some of this performance may be lost for systems with more massive secondaries. We also found that TES can handle close encounters such as the Apophis flyby in 2029 efficiently and accurately.