An implicit symplectic solver for high-precision long-term integrations of the Solar System

We present FCIRK16, a 16th-order implicit symplectic integrator for long-term high-precision Solar System simulations. Our integrator takes advantage of the near-Keplerian motion of the planets around the Sun by alternating Keplerian motions with corrections accounting for the planetary interactions. Compared to other symplectic integrators (the Wisdom and Holman map and its higher-order generalizations) that also take advantage of the hierarchical nature of the motion of the planets around the central star, our methods require solving implicit equations at each time-step. We claim that, despite this disadvantage, FCIRK16 is more efficient than explicit symplectic integrators for high-precision simulations thanks to: (i) its high order of precision, (ii) its easy parallelization, and (iii) its efficient mixed-precision implementation which reduces the effect of roundoff errors. In addition, unlike typical explicit symplectic integrators for near-Keplerian problems, FCIRK16 is able to integrate problems with arbitrary perturbations (non-necessarily split as a sum of integrable parts). We present a novel analysis of the effect of close encounters in the leading term of the local discretization errors of our integrator. Based on that analysis, a mechanism to detect and refine integration steps that involve close encounters is incorporated in our code. That mechanism allows FCIRK16 to accurately resolve close encounters of arbitrary bodies. We illustrate our treatment of close encounters with the application of FCIRK16 to a point-mass Newtonian 15-body model of the Solar System (with the Sun, the eight planets, Pluto, and five main asteroids) and a 16-body model treating the Moon as a separate body. We also present some numerical comparisons of FCIRK16 with a state-of-the-art high-order explicit symplectic scheme for 16-body model that demonstrate the superiority of our integrator when very high precision is required.


Introduction
High-precision long-term dynamic simulation of planetary systems requires computationally expensive numerical integrations. Efficient computation of long-term ephemerides of the Solar System can be achieved by exploiting the near-Keplerian motion of the planets around the Sun. The main idea is to alternate Keplerian motions with appropriate corrections accounting for the planetary interactions. This is the case of symplectic splitting integrators such as the Wisdom and Holman method (Wisdom and Holman 1991) and its higher-order generalizations (Wisdom et al. 1996;Blanes et al. 2013;Farrés et al. 2013;Mikkola and Palmer 2000).
For very long-term integrations, the propagation of roundoff errors is an important limiting factor of the final accuracy, and the standard 64-bit double precision arithmetic can be insufficient in some cases. An alternative that is used in some long-term integrations of the Solar System ) is provided by the 80-bit extended precision arithmetic (one bit for the sign, 15 bits for the exponent, and 64 bits for the significand). Compared to double precision arithmetic (having 53-bit significands) 11 additional binary digits of precision are gained by making all the computations in extended precision arithmetic at the cost of approximately doubling the computing time. For further reduction of roundoff errors, one could perform all the computations in quadruple precision at the cost of a drastic increase (by a factor of twenty or more) of computing time. A more efficient approach to gain some additional digits of precision (although not providing the full accuracy of quadruple precision,) is to use a mixed-precision methodology (Baboulin et al. 2009). That is, to use quadruple precision for the most critical computations and perform the rest of the computation in a lower precision arithmetic. In the case of explicit splitting methods, a natural mixed-precision strategy might be to compute the Keplerian flows in quadruple precision, and the corrections corresponding to the planetary interactions (which are typically of smaller magnitude) in 80-bit arithmetic. This would allow us to gain a few binary digits of precision, although with a considerable increase in computing time.
Further improvement of the efficiency of numerical simulations demands the development of fast and accurate algorithms that take advantage of parallel computer architectures. Although some attempts have been made in this sense for symplectic splitting integrators for Solar System simulations (Jiménez-Pérez and Laskar 2011), the sequential nature of that kind of schemes makes it difficult to get substantial improvements from parallelization strategies.
We present FCIRK16, a code for long-term high-precision integrations of systems with Hamiltonian function of the form H (q, p) = H K (q, p) + H I (q, p), where q = (q 1 , . . . , q n ), p = ( p 1 , . . . , p n ), q i , p i ∈ R 3 , i = 1, . . . , n. The Hamiltonian H K (q, p) is the sum of n uncoupled spatial Kepler problems. For each i = 1, . . . , n, μ i (resp. k i ) is the effective mass (resp. the gravitational constant) corresponding to the ith Keplerian problem. The interaction Hamiltonian H I (q, p) is seen as a perturbation of the Keplerian Hamiltonian H K (q, p). Our code implements a 16th-order implicit symplectic method within the class of FCIRK (flow-composed implicit Runge-Kutta) schemes proposed in Antonana et al. (2018). FCIRK methods are similar to symplectic splitting schemes in that Keplerian motions are alternated with corrections of smaller magnitude. But while in symplectic splitting integrators the cor-rections are computed as the solution operator of the interaction Hamiltonian, in FCIRK methods, such corrections correspond to the application to a transformed ODE system of one step of an implicit Runge-Kutta (IRK) method based on collocation with Gauss-Legendre nodes. As their underlying IRK methods, FCIRK methods are super-convergent: the method based on s Gauss-Legendre nodes is of order 2s. Compared to symplectic splitting integrators, the implementation of FCIRK methods is more involved because an implicit system of equations has to be solved at each time-step. In return, such methods are better suited than symplectic splitting integrators (i) for exploiting parallel computer architectures (most of the computations can be done in s processors in parallel), (ii) for the mixed-precision approach mentioned above, as considerably fewer Keplerian motions in quadruple precision have to be computed for the same level of precision in a given integration interval, (iii) for the integration of equations including additional effects, such as post-Newtonian corrections, which cannot be written as the sum of integrable parts.
In our preliminary implementation of the 16th-order FCIRK scheme (Antonana et al. 2018), time-steps of constant length were applied. However, for long-term simulations of realistic models of the Solar System, the larger bodies of the asteroid belt (specially Ceres, Pallas, and Vesta) experience close encounters that have a significant impact in the chaotic behavior of the Solar System . The precision of the numerical integration performed with any integrator with constant time-steps (in physical time) is degraded during close enough encounters, a degradation that may be more pronounced for the 16th-order FCIRK method due to its higher order of precision and to the temporary loss of the hierarchy H I << H K during close approaches.
Motivated by that, we have analyzed the effect of close encounters on the local discretization errors of our FCIRK integrators. More precisely, we have obtained rigorous estimates of the leading error term by applying results and techniques presented in Antonana et al. (2020).
Based on that analysis, we propose practical mechanism to identify integration steps that involve close encounters between arbitrary bodies. Actually, our analysis could be used as a theoretical basis to tray to endow our code with an adaptive time-stepping strategy. However, it is known (Calvo and Sanz-Serna 1993) that the advantages of symplectic integrators for the long-term integration of Hamiltonian systems are lost if standard adaptive time-stepping strategies are used.
We adopt a different approach in FCIRK16. Our code is mainly intended for long-term integrations of the Solar System, where close encounters that seriously degrade the accuracy are rare. In the absence of close approaches, FCIRK16 advances in the integration with constant time-steps so that the integration enjoys the advantages of symplectic integration during long integration subintervals. In order to deal with eventual close encounters, FCIRK16 identifies and resolves in higher precision time-steps which would otherwise suffer from accuracy degradation due to the occurrence of a close approach. This paper is organized as follows. Section 2 is devoted to define FCIRK methods and state their main properties. In Sect. 3, we analyze the effect of close encounters in the local discretization errors of FCIRK methods, and provide practical functions to monitor eventual close encounters. In Sect. 4, some aspects related to the implementation of FCIRK16 are briefly discussed. Section 5 is devoted to present some numerical experiments: In Sect. 5.1, we check the results and the ideas of Sect. 3 with a 15-body point-mass Newtonian model of the Solar System (with the Sun, the eight planets, Pluto, and five main asteroids) and with a 16-body model where the Moon is treated as a separate body. In Sect. 5.2, we present some numerical comparisons of FCIRK16 with a state-of-the-art high-order explicit symplectic scheme that demonstrate the superiority of our integrator when very high precision is required. We conclude with a summary of the work in Sect. 6.

FCIRK methods for perturbed Kepler problems
The equations of motion of the Hamiltonian function (1), written in terms of where u = (q 1 , . . . , q n , v 1 , . . . , v n ) ∈ R 6n , and with p i (i = 1, . . . , n) replaced by μ i v i in the right-hand sides of (3). This can be more compactly written in the form where k(u) corresponds to the Keplerian part and for each u ∈ R 6n In what follows, we will consider systems of the form (4) with arbitrary maps g i (u) not necessarily of the form (3) (so that (4) may include dissipative terms). We are concerned with the high-precision long-term numerical integration of (4) supplemented with the initial condition In FCIRK methods, as in the Wisdom and Holman (WH) integrator and its generalizations, one takes advantage of the integrability of the Keplerian equations of motion That is, one exploits the ability of efficiently computing for any h ∈ R (up to roundoff errors), the h-flow ϕ h : R 6n → R 6n of the unperturbed system Recall that, by definition of h-flow, ϕ h (u(t)) ≡ u(t + h) for any solution u(t) of (7). Application of the time-dependent change of variables u = ϕ t (w) transforms (4) into the non-autonomous system d dt w = (ϕ t (w)) −1 g(ϕ t (w)).
As shown in Antonana et al. (2018), any solution u(t) of (4) satisfies the identity where the map ψ h : R 6n → R 6n is defined as follows: for arbitrary u 0 ∈ R 6n , ψ h (u 0 ) = w(h), where w(t) is the solution of the initial value problem Hence, u(t + h) can be obtained from u(t) by two Keplerian motions of time-increment h/2 alternated with one application of the map ψ h that accounts for the effect of the interactions.
where the mapψ h is an approximation of the exact map ψ h accounting for the interactions.
In the WH integrator, In a FCIRK method,ψ h (u 0 ) is defined as a Runge-Kutta approximation of ψ h (u 0 ) = w(h) obtained by applying one step of length h to (9). More precisely, where the vectorsẆ i are implicitly defined bẏ Here, b i , c i , a i j (1 ≤ i, j ≤ s) are real parameters that determine the IRK scheme. The positive integer s is referred to as the number of stages of the IRK method. In Antonana et al. (2018), we considered a particular family of time-symmetric symplectic IRK methods, the collocation methods based on s Gauss-Legendre nodes. More precisely, c 1 , . . . , c s are the zeros of P s (2x −1), where P s (x) is the Legendre polynomial of degree s, and the coefficients b i (1 ≤ i ≤ s) and a i j (1 ≤ i, j ≤ s) are uniquely determined from the following equalities (Hairer et al. 2006): Notice that c 1 , . . . , c s and b 1 , . . . , b s are the nodes and weights, respectively, of the Gauss-Legendre quadrature formulas on the interval [0, 1]. The s-stage IRK method of collocation type with Gauss-Legendre nodes is of order 2s. That is,ψ as h → 0, where w (2s+1) (0) denotes the (2s +1)th order derivative of the solution w(t) of (9) evaluated at t = 0, andw (2s+1) (0) denotes the (2s + 1)th order derivative of the numerical solutionw(h) :=ψ h (u 0 ) evaluated at h = 0. In addition, if (4) are the equations of motion of a Hamiltonian function (1), then the map will be symplectic (as the map ψ h that is approximating to) (with respect to the 2-form i μ i dq i ∧ dv i ), and hence also each step u m−1 → u m of the form (10) of the FCIRK method (see Antonana et al. (2018) for further geometric properties of FCIRK methods).

Close encounters and local discretization errors
As mentioned in the introduction, the larger bodies of the asteroid belt (specially Ceres, Pallas, and Vesta) eventually experience close encounters, which degrade the accuracy of long-term simulations of realistic models of the Solar System if constant time-step implementation is used. We next try to estimate the effect of close encounters on the leading term of the local discretization error of FCIRK methods.

General analysis of local discretization errors
It is known that the local error (15) of a Runge-Kutta scheme (11)-(12) can be expanded as a B-series Hairer et al. (2006). In particular, w (2s+1) (0),w (2s+1) (0), and hence w (2s+1) (0) − w (2s+1) (0), can be written as a linear combination of elementary differentials associated to rooted trees with 2s + 1 vertices. We have checked that the coefficients of the later are substantially smaller (at least five orders of magnitude smaller) than the formers. In view of that, it is expected that the norm of the leading terms of the local error (15) will increase during close encounters in a similar way to (Here, we have collected the (6n)-vector w(t) in 2n three-dimensional sub-vectors as w(t) = (w 1 (t), . . . , w 2n (t)).) Motivated by that, we will next focus on obtaining upper bounds of (16).
In order to do that, we will apply some results and techniques (based on the analytic extension to the complex domain of N -problems) presented in Antonana et al. (2020).
Under the assumption that g : R 6n → R 6n is a real-analytic map, for each regular 1 initial state u 0 ∈ R 6n , there exists ρ = ρ(u 0 ) > 0 such that the following two conditions hold: C1 The solution u(t) = (q 1 (t), . . . , q n (t), v 1 (t), . . . , v n (t)) of (4)-(5) and g(u(t)) can be uniquely extended as analytic functions of t in the closed disk In that case, it can be proven (we provide a proof in the Appendix) that, for i = 1, . . . , n, where for i = 1, . . . , 2n,

Local error estimates for planetary systems in canonical heliocentric coordinates
Next, we will obtain suitable ρ(u 0 ) and in one particular relevant case: the case where (4) corresponds to a point-mass (n + 1)-body problem consisting of a massive body (the Sun) and n bodies (major and minor planets, and asteroids) orbiting around it, written in Poincaré's canonical heliocentric coordinates. Consider the Hamiltonian system corresponding to the (n+1)-body Hamiltonian function According to Theorem 2 in Antonana et al. (2020) (applied with λ = 1/7) the following statements S1-S2 hold for the Hamiltonian system supplemented with regular initial conditions S1 It admits a unique analytic solution ( with and The Hamiltonian H (q, p) obtained by rewriting the (n + 1)-body Hamiltonian (20) in terms of canonical heliocentric coordinates is of the form (1), with Here, we are assuming that the barycenter of the system is at rest, so that q 0 = (0, 0, 0) and v 0 = (0, 0, 0) can be ignored. The equations of motion, expressed in terms of the 3-vectors With that notation, the barycentric coordinates (Q i , V i ) (i = 0, 1, . . . , n) are recovered from the canonical heliocentric coordinates u = (q 1 , . . . , q n , v 1 , . . . , v n ) as The statements S1-S2 imply that C1-C2 hold for with where (27)-(28) is used to compute the barycentric coordinates (Q, V ) from u. We next give upper bounds for (19) required to get the estimates (17)-(18) of the (2s + 1)th derivatives w (2n+1) (0) featuring in the leading term of the local discretization error (15). The statements S1-S2 imply that, for each t ∈ D ρ(u 0 ) , From (32), one gets Lemma 1 in Antonana et al. (2020) and (31) finally imply that

The effect of close encounters in the local discretization errors
Consider, as in previous Subsection, an (n + 1)-body system described in canonical heliocentric coordinates, and assume that the mass ratios i = m i /m 0 (i = 1, . . . , n) are small. We want to study the effect of close encounters on the local discretization errors ) the 3-vector that collects the components of δ(u 0 , h) corresponding to the positions (resp. velocities) of the ith body.
While integrating with FCIRK, whenever all the bodies orbiting around the central massive body are well separated in w m (we next set u 0 := w m ), -The upper bounds (33)-(34) will remain small, and -Typically, ρ(u 0 ) = ρ 0,1 (u 0 ), where i = 0 corresponds to the central massive body (say, the Sun) and body i = 1 is nearest the central body (say, Mercury).
In that case, ρ(u 0 ) will oscillate with approximately the orbital period of Mercury, with lowest values of ρ(u 0 ) at Mercury's perihelion. In view of (15) and (17)- (18), the components of the local discretization errors are expected to oscillate likewise. However, if body i = k and j = are close enough to each other in u 0 , then ρ k, (u 0 ) will be smaller than the value of ρ 0,1 (u 0 ) at Mercury's perihelion, and thus ρ(u 0 ) = ρ k, (u 0 ). (In addition, among the upper bounds (33)-(34), those for i = k and i = will dominate over the rest for close enough encounters between bodies k and . ) We thus expect that the local error will grow roughly as ∼ (h/ρ(w m )) 2s+1 when ρ(w m ) decreases due to a close encounter. This suggests that -Close encounters that may degrade the accuracy of the FCIRK integration (10) can be identified by monitoring the value of ρ(w m ), and -In such cases, the local discretization could be controlled by decreasing the step size h of the IRK scheme as ρ(w m ) decreases to keep h/ρ(w m ) roughly constant.

The case of models of the Solar System with the moon as a separate body
The discussion in previous two subsections was applicable to models of the Solar System where the Earth-Moon system is treated as a point mass located at their barycenter. We now focus on the case of (n + 1)-body models of the Solar System with the Earth and the Moon treated as two separate point masses, say, bodies i = n − 1 and i = n, respectively. We consider canonical Heliocentric coordinates for bodies i = 1, 2, . . . , n − 2 and for the Earth-Moon barycenter, and geocentric positions and corresponding conjugate momenta for the Moon. The (n + 1)-body problem with Hamiltonian (20) described in such canonical coordinates is of the form (1), where and The barycentric coordinates Q i , . . , n, as follows: Q n−1 = Q 0 + q n−1 − m n m n−1 q n , Q n = Q 0 + q n−1 + q n , It is not difficult to prove, as in Sect. 3.2, that C1-C2 also hold in that case, with ρ(u) given by (29)-(30), where L i, j (Q, V ) is given by (23) and (35) is applied to compute the barycentric coordinates (Q, V ) from u. Upper bounds of M i (u) similar to those obtained in Sect. 3.2 can also be obtained. The local error (15) can thus be estimated with the help of (33)-(34).
We eventually arrive at the conclusion that the local error will grow roughly as ∼ (h/ρ(w m )) 2s+1 when ρ(w m ) decreases due to some close encounter. During most of the simulation, ρ(w m ) = ρ n−1,n (w m ) due to the highly oscillatory nature of the orbit of the Moon, and only for very extreme close encounters between body k and does ρ(w m ) decrease noticeably (as ρ k, (w m ) becomes smaller than ρ n−1,n (w m ) and consequently ρ(w m ) < ρ n−1,n (w m )). However, we have observed in practice that some close encounters that degrade the accuracy of the simulation are not close enough to have ρ(w m ) < ρ n−1,n (w m ), and hence cannot be detected by monitoring the value of ρ(w m ). We believe that this is due to the fact that the contribution of the high oscillation of the Moon on the local discretization errors of the rest of the bodies is typically small. Moreover, the contribution of such highly oscillatory local errors over accumulated errors is further reduced because they are partially averaged out.
Motivated by that, we propose to modify the function ρ(u) for monitoring close encounters when the Moon is treated as a separated body by excluding ρ n−1,n (u) from the definition of ρ(u).

Implementation of the IRK approximation
In our implementation of FCIRK methods, we compute the vectorsẆ i from (12) by using the fixed-point implementation of IRK schemes presented in Antonana et al. (2017). This consists of a robust and efficient implementation that avoids the systematic accumulation of errors in energy (due to roundoff errors and iteration errors) observed in standard fixed-point implementations of symplectic implicit RK schemes (Hairer et al. 2008;Antonana et al. 2017).
We compute the initial guesses to start the fixed-point iteration with an extrapolation procedure that allows us to substantially reduce the required number of fixed-point iterations for small enough time-steps h. However, for FCIRK methods with large number s of implicit stages, this is only true for extremely small values of the step size. Furthermore, the extrapolation formulae become numerically unstable (due to cancellation errors caused by large coefficients of opposite sign) for more large values of s. These two factors limit, for practical computations in 64-bit or 80-bit floating point arithmetic, the efficiency of FCIRK methods for large values of s, say s ≥ 12. After some preliminary numerical experiments, we have concluded that s = 6, 7, 8, 9 are good choices in this respect. Similar considerations apply to implicit Runge-Kutta schemes of collocation type: in Hairer et al. (2006), the choice s = 6 is favored for 64-bit floating point arithmetic in the case of Gauss-Legendre nodes, and s = 8 is chosen in Everhart (1985), Rein and Spiegel (2014) in the case of Gauss-Radau nodes. For FCIRK methods, we finally choose the scheme with s = 8 stages (hence of order 2s = 16) because it allows a more efficient parallelization with four or eight threads.
We refer to our implementation of the 16th-order FCIRK method as FCIRK16. In FCIRK16, the evaluation of the right-hand side of (12) for each i ∈ {1, . . . , 8} is optionally performed in parallel (with openMP) with a prescribed number of threads.

Mixed-precision implementation of FCIRK16
A mixed-precision approach proposed in Antonana et al. (2018) to reduce the effect of using finite precision arithmetic is implemented in FCIRK16. The main idea consists of applying most of the computations in a basic floating point arithmetic and performing a few critical operations in a higher precision arithmetic.
In our current implementation of FCIRK16,ŵ m =ψ h (w m ) in (11) is evaluated as where Φ(u m , h) is evaluated in 80-bit extended precision, while the rest of the computations in (10), that is, the evaluations of the map ϕ h/2 and the sum (36) are performed in quadruple precision (128-bit floating point arithmetic). This is motivated by the following: -The evaluation of Φ(w m , h) represents the bulk of the computations of each step (10). Note that this is true regardless of the complexity of the evaluation of the interaction term g(u) in (4). Indeed, if m fixed-point iterations are needed at each step (typically, m = 5, 6), sm = 8m evaluations of ϕ h and its partial derivatives (in addition to 8m evaluations of g(u)) are required to evaluate Φ(w m , h), compared to two evaluations of ϕ h (and 6n additions) in quadruple precision. Hence, the increase in computing time of one step (10) due to the mixed-precision implementation compared to one step fully implemented in 80-bit extended precision is always relatively small. -The components of the vector Δ m in (36) are typically of considerably smaller magnitude than the corresponding component of w m . If their magnitude is, say, 2 −k smaller, then the overall precision of each step (10) will correspond to a significand of 64 + k bits compared to a precision of 64-bit significand one would have if all the computations in (36) were performed in 80-bit floating point arithmetic.
It is worth mentioning that, for planetary systems, canonical heliocentric coordinates are better suited than Jacobi coordinates for applying the mixed-precision implementation described above. This is due to the fact that in the case of Jacobi coordinates, cancellation errors will typically occur in the evaluation of Δ m , which may increase considerably the effect of roundoff errors. Motivated by that, we will focus on the application of planetary systems described in canonical heliocentric.

Dealing with close encounters
We aim at enjoying the favorable behavior of constant time-step symplectic integration during long integration subintervals, and only occasionally (during eventual close encounters) perform a more dedicated (computationally more expensive) numerical integration to compute u m ≈ u(mh) from u m−1 . We thus have endowed FCIRK16 with a mechanism to detect the time-steps (from now on, critical time-steps) involving a close encounter that would substantially degrade the accuracy of the numerical integration in constant step size mode.
In such critical time-steps,ψ h (w m ) in (10) is replaced by a better approximation of ψ h (w m ). While in ordinary stepsû m =ψ h (w m ) is obtained by applying one step of length h of the underlying IRK scheme applied to the transformed system, in critical steps,û m will be computed by applying k steps of constant length h/k of the IRK scheme to the transformed system. In addition, since the mixed-precision approach becomes less effective during close encounters (as the dominance of the Keplerian Hamiltonian H K over the interaction Hamiltonian is reduced or lost during such events), we implement critical time-steps fully in quadruple precision 128-bit floating point arithmetic.
We next specify the criterion to identify critical time-steps that we have implemented in FCIRK16, and the rule used to choose the number k of substeps of the IRK schemes to be applied in each critical step. Both tasks are done with the help of the monitoring function ρ(u) defined in Sect. 3.2 for (n + 1)-body problems described in heliocentric coordinates, and the modified monitoring function (that we will also refer to as ρ(u)) defined at the end of Sect. 3.4 for the case where the Moon is treated as a separate body. For systems of the form (4) obtained by describing an (n + 1)-body problem in other coordinates the monitoring function ρ(u) could be defined similarly.
In FCIRK16, the mth step (10) will be considered critical if, once where μ (resp. σ ) is the arithmetic mean (resp. the standard deviation) of the values of ρ(w k ) for all previous ordinary steps, and ν is a prescribed positive number. We set ν = 1.6 as default value in FCIRK16. Motivated by second item at the end of Sect. 3.3, we choose the number k of substeps of constant length h/k of the IRK scheme to be performed when computingŵ m from w m , as the positive integer k such that

Numerical experiments
We next present some numerical experiments: The experiments in the first subsection are aimed to assess the effectiveness of the technique for monitoring close encounters introduced in Sect. 3. In the second subsection, we demonstrate the performance of FCIRK16 for highprecision long-term integrations of the Solar System by comparing it with a state-of-the-art explicit symplectic integrator. We consider two Newtonian point-mass (n + 1)-body models of the Solar System: From one hand, a 15-problem with the Sun, the eight planets, Pluto, and five main bodies of the asteroid belt (Ceres, Pallas, Vesta, Iris, and Bamberga) described in Poincare's canonical heliocentric coordinates. On the other hand, we consider the 16-problem obtained from the former by considering the Moon as a separate body, with canonical coordinates described in Sect. 3.4.
We consider the initial values at Julian time (TDB) 2440400.5 (the 28th of June of 1969), from the DE430 Ephemerides (Folkner et al. 2014).
All the numerical tests presented here are run in a HP-Z6-G4-Workstation with 8 cores of Intel©Xeon©Silver 4110 CPU @ 2.10 GHz and 32GiB (16 +16 DIMM DDR4) of RAM memory. The code is written in C and has been compiled under Ubuntu 18.04.3 LTS operating system with gcc version 7.5.0 with options -O2 -std=c99 -fno-common -mfma. The parallelized implementation of FCIRK16 is based on the OpenMP library.

Monitoring close encounters
We first conduct some experiments to check the techniques of monitoring close encounters proposed and discussed in Sect. 3. For that purpose, we have first made two integrations of the 15-body problem with FCIRK16 forward in time for 50000 years with a constant time-step of h = 1.5 days and h = 3 days, respectively.
The evolution of relative errors in energy is displayed in the upper subplot of Fig. 1, while the middle (resp. lower) subplot shows the evolution of the errors in position 2 for each of the bodies orbiting around the Sun for the non-adaptive (adaptive) integration with constant step size h = 1.5.
Some jumps in the energy errors are clearly observed in the integration with h = 3. (The largest jump coincides with a close encounter between Pallas and Vesta.) On the contrary, there is no clear sign of accuracy degradation in the plot of the evolution of relative errors in energy for h = 1.5. However, a substantial loss of precision in the positions of some of the bodies can clearly be seen in the middle subplot of Fig. 1. We have also made an analogous Fig. 1 Evolution of relative errors in energy for FCIRK16 with h = 1.5 and h = 3 (top), and evolution of errors in position of each of the body for the non-adaptive (resp. adaptive with ν = 1.6) implementation of FCIRK16 with h = 1.5 in the middle (resp. lower) subplot experiment for the 16-body problem, and obtained very similar plots. We conclude that monitoring the energy error is not enough to identify critical close encounters between minor bodies in the asteroid belt. This is due to the marginal contribution of the asteroids to the total energy of the system.
The lower subplot of Fig. 1 shows that the technique to reduce the accuracy loss due to close encounters proposed in Sect. 4.3 has been effective in this example. We have used the default value ν = 1.6, which has resulted in 203 critical steps out of 122×10 5 ordinary steps. It is worth stressing that, due to the low percentage of critical steps, applying our adaptive technique does not have any practical impact in the execution time of the integration.
We next take a closer look to the effectiveness of the monitoring technique proposed in Sect. 3. Figure  The sharpest spikes of local error (corresponding to close encounters between asteroids) nicely coincide in time with the sharpest local minima of the monitoring function. We conclude that the proposed monitoring technique works very well in this case. In Fig. 3, the results of a similar experiment performed for the 16-body problem are reported. This time, the three curves (μ, μ − σ , and min ρ) corresponding to two different monitoring functions are displayed: (i) the original monitoring function in the upper subplot, (ii) and the modified monitoring function proposed at last item at the end of Sect. 3.4 (which excludes the Earth-Moon interaction) in the middle subplot. In the lower subplot, the evolution of the estimated local errors in positions of each of the bodies is displayed. The upper nearly horizontal curve corresponds to the local error in position of the Moon. (The local errors in position of Mercury and the Earth are five orders of magnitude smaller.) The sharp local minima of the original monitoring function identify correctly the spikes of local error that arise well above the local error of the Moon, but some of the less pronounced spikes are not reflected in the original monitoring function. The modified monitoring function, displayed in the middle subplot, is able to identify nicely the spikes of local error that arise above Mercury's local error.

Comparison with a state-of-the art explicit symplectic integrator
The efficiency of several generalizations of the classical Wisdom and Holman integrators (in Jacobi coordinates) for high-precision long-term integration of the Solar System are compared in Rein et al. (2019). It is there concluded that, for high-precision integrations, the symplectic integrator ABA(10,6,4) presented in Blanes et al. (2013), Farrés et al. (2013) is the most efficient among all the considered methods. The integrator ABA(10,6,4) requires that the interaction Hamiltonian H I only depend on the positions (or at least, that H I can be written as the sum of two commuting integrable Hamiltonians). This is for instance the case if Jacobi coordinates are used to rewrite the N -body Hamiltonian (20) in the form (1). If Poincare's canonical Heliocentric coordinates are used instead, the interaction Hamiltonian is of the form In Blanes et al. (2013), the integrator ABAH(10,6,4) is constructed as an alternative to ABA(10,6,4) when Poincaré's Heliocentric coordinates are used instead of Jacobi coordinates (more generally, for interaction Hamiltonians of the form (39)). In the numerical integrations of the Solar System performed in Farrés et al. (2013), it is concluded that ABAH(10,6,4) applied with Heliocentric coordinates is as efficient as ABA(10,6,4) applied with Jacobi coordinates.
The three positive integers in the name of the integrators ABA(10,6,4) and ABAH(10,6,4) refer to their global error estimates, which are in both cases of the form O( h 10 + 2 h 6 + 3 h 4 ), where is the maximum of the mass ratios m i /m 0 . In contrast, the global error of FCIRK16 is of the form O( h 16 ).
In the case of the two models of the Solar System that we are considering in our numerical experiments (that is, the interaction Hamiltonians described in Sects. 3.2 and 3.4,) are also of the form (39). We will thus compare FCIRK16 with the symplectic integrator ABAH(10,6,4) presented in Blanes et al. (2013), Farrés et al. (2013).
The splitting symplectic integrator ABAH(10,6,4) is implemented by alternating evaluations of the Keplerian motions with nine evaluations of the gradients of A(q) and B( p) in (39). No symplectic correction is applied. We apply our own optimized implementation of the later in 80-bit extended precision arithmetic (with Kahan's compensated summation (Kahan 1965)), which makes use of the same routines for the evaluation of ϕ h and the gradient of the interaction Hamiltonian as in FCIRK16. We do not attempt the implementation of parallelization techniques for ABAH(10,6,4): splitting integrators, unlike the FCIRK methods, require the sequential evaluation of the gradients of the interaction Hamiltonians. Consequently, parallelization of ABAH(10,6,4) is only possible inside the evaluation of the gradient, hence with a lower granularity than that of FCIRK16, with would have a considerably higher parallelization overhead.
The technique proposed in Sect. 4.3 can easily be adapted to splitting symplectic integrators and in particular to the method ABAH(10,6,4). For a fair comparison with FCIRK16, we have implemented ABAH(10,6,4) with a similar mechanism to deal with close encounters. In order to illustrate the advantage of incorporating that technique into our implementation of ABAH(10,6,4), we have made two integrations of the 15-body problem considered in the previous subsection with step size h = 1, one with fully constant step size, and the other one with our technique for detecting and coping with close encounters. In the upper plot of Fig. 4, we display the evolution of the position errors of each of the orbiting bodies for the integration with fully constant step size h = 1. One can observe a clear degradation of the precision of the positions of some of the bodies due to close encounters. In the lower plot of Fig. 4, we can see the evolution of the errors in the integrations corresponding to ν = 1.6 (304 critical steps out of 183 × 10 5 ordinary steps have been required in that case).
In Antonana et al. (2018), the efficiency of a preliminary implementation of our 16thorder FCIRK method was compared to that of ABAH(10,6,4) (and also to the 16th-order implicit Runge-Kutta implemented in Antonana et al. (2017)). A simple 10-body model of the Solar System (including no asteroids, nor the Moon as a separate body) was considered. In particular, efficiency diagrams displaying the error in energy versus execution time for different values of the constant time-step were presented for that purpose. We have repeated such experiments for the 15-body model and 16-body model of the Solar System described above for our current implementation of FCIRK16 (with the default value ν = 1.6 in (37)). We have obtained very similar efficiency diagrams (not shown here), the main difference (compared to those presented in Antonana et al. (2017)) being that, due to the higher complexity of the evaluation of the gradient of the interaction Hamiltonian, higher parallelization speedup is obtained for the execution of FCIRK16. The parallel execution of FCIRK16 with four threads is (both in the 15-body problem and the 16-body problem) approximately twice as fast as the sequential execution.
As in Antonana et al. (2018), we conclude that ABAH(10,6,4) is more efficient than FCIRK16 for lower precisions. However, if higher precision is required, ABAH(10,6,4) is limited by the precision of the 80-bit floating point arithmetic. At that limit, thanks to its mixed-precision implementation, FCIRK16 is able to reduce errors by two orders of magnitude with similar computing time. We have also tested a mixed-precision implementation ABAH(10,6,4) that computes all the Keplerian flows in quadruple precision, but its execution time increases by a factor between 6 and 9 in our examples. This is an expensive toll to be paid for a little of extra precision, while in FCIRK16, this is achieved essentially for free.
In order to illustrate the superiority of FCIRK16 over ABAH(10,6,4) near the limit of precision of the 80-bit floating point arithmetic, we next compare the propagated errors We will next show that the ABAH(10,6,4) integration with h = 1 (resp. h = 0.5), which is comparable in computing time to the sequential (resp. parallel) execution of FCIRK, gives less precise results than FCIRK with h = 3. We first display the evolution of relative errors in energy for each of the three integrations in Fig. 5. Figure 6 (resp. Fig. 7) shows the evolution of the errors in position (resp. in eccentricities) of each of the bodies: for FCIRK16 with h = 3 (top), for ABAH(10,6,4) with h = 0.5 (middle), for ABAH(10,6,4) with h = 1 (bottom). The evolution of errors in velocity and in semi-major axis (not included here,) show a behavior similar to that of positions and eccentricities, respectively.
Comparison of the evolution of errors for ABAH(10,6,4) with h = 1 and h = 0.5, suggests that no substantial increase in the precision could be expected for ABAH(10,6,4) with step sizes smaller than h = 0.5. We conclude that for that accuracy regime (limited by the precision of the 80-bit floating point arithmetic), FCIRK is able to produce (about two orders of magnitude) more precise results than ABAH(10,6,4) for similar execution times. We stress that, as illustrated in Fig. 4, our implementation of ABAH(10,6,4) also benefits greatly from the ideas and techniques introduced in Sects. 3 and 4.3. The evolution of errors for the fully constant step size implementation of ABAH(10,6,4) would be considerably less favorable due to the occurrence of close encounters between asteroids.

Summary of the work
We have presented FCIRK16, an integrator for high-precision long-term integration of the Solar System. It consists of an optimized and robust implementation of a 16th-order scheme within the class of FCIRK methods presented in Antonana et al. (2018). Such methods are implicit schemes that take advantage of the near-Keplerian motion of the planets around the Sun by alternating Keplerian motions with corrections accounting for the planetary interactions.
We have presented a novel analysis of the local errors of FCIRK methods intended to understand how discretization errors increase during close encounters. That allows us to endow FCIRK16 with practical mechanisms to monitor and accurately resolve eventual close encounters.
We have compared our integrator with a state-of-the-art high-precision long-term integrator (which is an explicit scheme that also takes advantage of the near-Keplerian motion of the planets), the symplectic splitting method ABAH(10,6,4) Farrés et al. 2013). Compared to ABAH(10,6,4) (or any available explicit symplectic integrator), our new integrator FCIRK16 -Is better suited to parallel implementations, -Admits the application of a very efficient mixed-precision technique to reduce the accumulation of roundoff errors, -Admits any kind of interaction Hamiltonian, not necessarily the sum of integrable parts.
-Includes a mechanism to detect and effectively cope with eventual close encounters.
We present numerical experiments on a point-mass 16-body Newtonian model (including the Moon as a separate body) of the Solar System to compare our integrator to ABAH(10,6,4) for an integration backward in time of 350 thousand years. We compare our own optimized implementation in 80-bit extended precision of the ABAH(10,6,4) with a mixed-precision implementation (that combines 80-bit extended precision with 128-bit quadruple precision) of FCIRK16. Actually, we have compared FCIRK16 with an improved implementation of ABAH(10,6,4) endowed with the same mechanism to detect and accurately resolve close encounters.
In our numerical experiments with that 16-body of the Solar System, we conclude that: (i) ABAH(10,6,4) is more efficient than FCIRK16 for the regime where local truncation errors of ABAH(10,6,4) clearly dominate over the contribution of roundoff errors, (ii) for more precise integrations, FCIRK16 is more efficient than ABAH(10,6,4), (iii) the accuracy of ABAH(10,6,4) is limited by the precision of the 80-bit extended precision floating point arithmetic (and we do not see how to gain some digits of precision without compromising seriously its efficiency), (iv) in that limit, thanks to application of the mixed-precision technique, FCIRK16 is able to improve the precision with similar computing time. We stress that our implementation of ABAH(10,6,4) also benefits from the analysis of Sect. 3 and the technique for resolving close encounters proposed in Sect. 4.3. The superiority for high precisions of FCIRK16 over the standard implementation (with fully constant time-steps) of ABAH(10,6,4) would be considerably higher.
We believe that, for very long integrations such as , the abovementioned extra accuracy may be critical in some cases. For such high precision, it may certainly make sense to use a more realistic model of the Solar System instead of the simple Newtonian point-mass model we have considered in our numerical tests. That scenario actually favors FCIRK16 over ABAH(10,6,4). Indeed, due to the higher computational cost of the evaluation of g(u) of more realistic models of the Solar System (including expensive relativistic corrections), (i) the relative parallelization overhead of FCIRK16 should be greatly reduced (hence getting a higher parallelization speedup) compared to the simple model of the Solar System we have considered in our numerical experiments; (ii) in addition, the fixed-point implementation could be adapted to perform fewer evaluations of the relativistic corrections per step. These two factors would potentially improve the relative efficiency of FCIRK16 compared to ABAH(10,6,4). In addition, the inclusion of physical effects that cannot be written as the sum of integrable parts may hinder the efficient implementation of ABAH(10,6,4), while this poses no problem to FCIRK16.
FCIRK16 is intended to be applied for planetary systems, like the Solar System, having low eccentricity orbits, and where close encounters between the modeled bodies are relatively rare. For higher eccentricity planetary systems or systems with potentially frequent close encounters, a fully adaptive version of FCIRK integrators should be implemented. We believe that the analysis and ideas presented in Sect. 3 could be helpful to pursue that goal.
The current C language implementation of FCIRK16 can be downloaded from our Github software repository: https://github.com/mikelehu/FCIRK/. User documentation is provided and many examples are included, ensuring that the results in the present work are reproducible (LeVeque et al. 2012).