High precision symplectic integrators for the Solar System

Using a Newtonian model of the Solar System with all 8 planets, we perform extensive tests on various symplectic integrators of high orders, searching for the best splitting scheme for long term studies in the Solar System. These comparisons are made in Jacobi and heliocentric coordinates and the implementation of the algorithms is fully detailed for practical use. We conclude that high order integrators should be privileged, with a preference for the new \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(10,6,4)$$\end{document} method of Blanes et al. (2013).


Introduction
Due to their simplicity and stability properties, symplectic integrators have been widely used for long-term integrations of the Solar System, starting with the work of Wisdom and Holman (1991). In many studies on the formation and evolution of the Solar System, where large numbers of particles are considered, the speed of the integrator is a major constraint and low order schemes have been often used as in the original scheme of Wisdom and Holman (1991) or Kinoshita et al. (1991) (for a review see Morbidelli 2002).
On the opposite, in the present work we are focusing on high precision symplectic integrators that are designed for the computation of long term ephemerides of the Solar System, when one searches to reduce the numerical error of the algorithm to the level of the roundoff error of the machine. These integrators will also be useful for the detailed dynamical studies of the extra solar planetary system with strong planetary interactions.
The first long term direct numerical integration of a realistic model of the Solar System, including all planets and the effects of general relativity and the Moon was made twenty years ago over 3 Myr (Quinn et al. 1991) using a high order symmetric multistep method. This solution could be compared with success with the previous averaged solutions of Laskar ( , 1990a and confirmed the existence of secular resonances in the Solar System (Laskar et al. 1992). Soon after, using a symplectic integrator with mixed variables (Wisdom and Holman 1991), Sussman and Wisdom (1992) could extend these computation to 100 Myr, confirming the chaotic behaviour of the Solar System discovered with the secular equations by Laskar ( , 1990a. As the Solar System is chaotic, the error in numerical integrations is multiplied by 10 every 10 Myr . Due to the limited accuracy of the models and initial conditions, it is thus hopeless to obtain a precise solution for the evolution of the Solar System over more than 100 Myr. The situation is even worse when the full Solar System is considered, as close encounters among the minor planets induce strong chaotic effects that will limit all possibilities of computing a precise solution for the planets to about 60 Myr (Laskar et al. 2011a,b).
Despite this limitation, there is a strong need for precise ephemerides of the planets from the paleoclimate community. Indeed, the variations of the Earth orbital elements induce some changes in the Earth climate that are reflected in the sedimentary records over million of years. This mechanism, known as Milankovitch theory (Milankovitch 1941) allows now to use the astronomical solution for the calibration of the geological time scales through the correlation of the variation of orbital and rotational elements of the Earth to geological records. This method has been successfully used for the Neogene period (Lourens et al. 2004) over 23 Myr, and a large effort is pursued at present to extend this study over the full Cenozoic era, up to about 65 Myr. This quest led to search for high order symplectic schemes that are adapted to these long time computations, where high accuracy is requested (Laskar and Robutel 2001;Laskar et al. 2004Laskar et al. , 2011a, but it should be noted that in the latest work, the integration of the Solar System model over 250 Myr, 1 including five main asteroids took more than 18 months of CPU time. Some improvements of the algorithms were thus needed, and the present paper is the outcome of the studies that we have undertaken in order to search for the best integrators for the next generations of numerical solutions. At the same time, we have compared various sets of coordinates (heliocentric and Jacobi), as the performances of these integrators depend on the choice of splitting of the Hamiltonian, and thus of the set of coordinates that correspond to these various splittings. As the integrators that are presented here are of high order, they can also be used for refined analysis of the newly discovered extra solar planetary systems, especially when the planetary interactions in the system are strong.
For the planetary case, when using an appropriate set of coordinates, the equations of motion are written as an integrable part H A , that corresponds to the Keplerian motion of each planet, and a small perturbation H B , given by the interaction of the planets between each other. Hence, the system falls into the category of Hamiltonian system of the kind Several splitting integrating schemes that take advantage of this fact to derive efficient integrators exist in the literature. McLachlan (1995) was the first to present such schemes and was followed independently by Chambers and Murison (2000) and Laskar and Robutel (2001). Most recently (Blanes et al. 2013) derived higher order schemes that present very interesting behaviour.
In this paper we describe these different splitting symplectic schemes and compare them for the case of the Solar System dynamics. We want to see which are the most efficient and accurate schemes. We will consider the gravitational N-body model and test the different integrating schemes against different planetary configurations, to be more specific: the 4 inner planets, the 4 outer planets and the 8 planets in the Solar System (Sect. 4).
The Hamiltonian of the gravitational N-body problem H = T ( p) +U (q) can be rewritten as H = H A + ε H B , using two different set of canonical coordinates: Jacobi and heliocentric coordinates (Sect. 3). The main difference between both sets of equations is that in Jacobi coordinates the small perturbation H B depends only on positions while in heliocentric coordinates the perturbation depends on both position and velocity. This is why in the literature Jacobi coordinates have been more widely used. In Sect. 5 we describe different symplectic schemes for Jacobi coordinates, and in Sect. 6 other symplectic schemes that are suitable for heliocentric coordinates. In both sections we describe and compare the different splitting schemes. Finally in Sect. 7 we compare the results for the two different set of coordinates.

Splitting symplectic integrators (general overview)
Let H (q, p) be a Hamiltonian system where (q, p) are a set of canonical coordinates (i.e. q are the positions and p the momenta). It is well known that in many mechanical problems the Hamiltonian is of the form which is separable with respect to the local canonical coordinates. Using the Lie formalism we can write the equations of motion as: where by definition L χ f := {χ, f } is the differential operator L χ , z = (q, p) and {·, ·} denotes the Poisson bracket. 2 The formal solution of Eq. 1 at time t = τ 0 + τ that starts at time t = τ 0 is given by In general the operators L T and L U do not commute, exp(τ (L T + L U )) = exp(τ L T ) exp(τ L U ), but we can find coefficients a i , b i such that for a given r , (3) Using the Baker-Campbell-Hausdorff (BCH) identity we can find relations that the coefficients a i , b i must satisfy to have a high order scheme (Koseleff 1993b(Koseleff , 1996. These are the so-called order conditions. For a given set of coefficients a i , b i satisfying Eq. 3, the composition is a symplectic map of order r . The map S(τ ) is symplectic because it is the product of elementary symplectic maps, exp(τ L T ) and exp(τ L U ), and has order r because it approximates the exact solution up to order τ r . We will refer to these kind of symplectic schemes as splitting symplectic integrators.
Some of the main advantages of these kind of integrating schemes are: (a) they are very easy to implement; (b) they preserve the symplectic character of the Hamiltonian system; and (c) in general there is no systematic drift on the conservation of the energy during the numerical integration.
These kind of symplectic schemes have been widely studied throughout the years by several authors (see Hairer et al. 2006;McLachlan and Quispel 2002 and references therein). As a matter of fact, splitting methods have been designed (often independently) and extensively used in fields as diverse as molecular dynamics, simulations of storage rings in particle accelerators, quantum chemistry and, of course, Celestial Mechanics.
There are several procedures to get the order conditions for the coefficients of the splitting scheme in Eq. 4. These are, generally speaking, large systems of polynomial equations in the coefficients that are obtained from Eq. 3. Two of the most popular are the recursive application of the BCH formula to the composition in Eq. 4, and a generalisation of the theory of rooted trees used in the analysis of Runge-Kutta methods due to Murua and Sanz-Serna (1999) (see also Hairer et al. 2006). The latter procedure, while being more systematic than the former, is however not appropriate for the case of splitting methods applied to Hamiltonians of the form A + B. In Blanes et al. (2013) a novel systematic way is proposed based on the so-called Lyndon multi-indices that is well adapted to that case.
Splitting methods of order greater than two involve necessarily some negative coefficients a i and b j (Goldman and Kaper 1996;Sheng 1989;Suzuki 1991). Although this feature does not imply in principle any special impediment for the class of systems considered in this paper, it is clear that the presence of negative coefficients may affect the numerical error and the maximal step size of the scheme. For this reason, when dealing with high order methods, minimising the size of the negative coefficients and the sum of the absolute value of all the coefficients will be a critical factor in the choice of one particular set of coefficients.
In this paper we do not intend to give the details on the derivation of the order conditions or how to find these coefficients. All these issues are analysed in detail in Blanes et al. (2013). Our aim here is to compare the performance of different splitting symplectic schemes for the specific case of the integration of the Solar System.
If we focus on the motion of the Solar System, or other planetary systems, we have a main massive body in the centre (the Sun) and the other bodies evolve around the centre mass following almost Keplerian orbits. We can take advantage of this to build more efficient schemes. Using an appropriate change of coordinates we can rewrite the Hamiltonian as, , a sum of the Keplerian motion of each planet around the central star and a small perturbation due to the interaction between the planets, where H K and H I are integrable. Wisdom and Holman (1991), Kinoshita et al. (1991) were the first to split the Hamiltonian of the N-body problem in this way for numerical simulations of the Solar System, by means of what is called a mixed variable integrator, using elliptical coordinates to integrate the Keplerian motion. Splitting the Hamiltonian as H K + H I rather than the classical T ( p)+U (q) decomposition 3 already improves the performance of the leapfrog scheme. As |H I | is small with respect to |H K |, the system falls into the class of Hamiltonian such that H = H A + ε H B for ε small. In this particular case, the truncation order of the leapfrog scheme is no longer Cτ 3 as for T ( p) + U (q), but rather C ετ 3 (McLachlan 1995;Laskar and Robutel 2001).
In Sects. 5 and 6 we will describe different families of symplectic splitting methods for Hamiltonian systems of the kind H A + ε H B and we will compare their performance for the particular case of the Solar System dynamics.

The N-body problem
Throughout this article we consider the non-relativistic gravitational N-body problem as a test model for the different integrating schemes. We are aware that to have a realistic model for the Solar System dynamics one must include effects like general relativity or tidal dissipation. Nevertheless, and for the sake of simplicity, in this paper these effects are ignored as their presence should not compromise the performance of the schemes presented here.
In a general framework, we consider the motion of n + 1 particles: the Sun and n planets, that are only affected by their mutual gravitational interaction. Let u 0 , u 1 , . . . , u n anḋ u 0 ,u 1 , . . . ,u n be the position and velocities, in a barycentric reference frame, of the n + 1 bodies and let m 0 , m 1 , . . . , m n be their respective masses. For simplicity, we consider m 0 to be the mass of the Sun and m i for i = 1, . . . , n the mass of the other planets.
Taking the conjugated momentaũ i = m iui , the equations of motion are Hamiltonian, with: In this set of coordinates the Hamiltonian naturally splits into, H = T +U , where T depends only on the momenta (ũ i ) and U depends only on the positions (u i ).
In general, when we deal with complex dynamical systems, it is important to take into account the relevant aspects of the system and use them to build efficient numerical tools to describe their dynamics. In the case of the Solar System we have a massive body in the centre and the planets evolve following Keplerian orbits around it that vary through time due to their mutual interaction.
Using an appropriate change of variables the Hamiltonian can be written as H K + H I , where |H I | is small with respect to |H K |, and both parts are integrable when we considered them on their own. There exist two canonical set of coordinates that allow us to split the Hamiltonian in this way: Jacobi and heliocentric coordinates.

Jacobi coordinates
The Jacobi set of coordinates has been widely used in Celestial Mechanics for developing analytical theories for the planetary motion. It was first used for the numerical integration of the Solar System by Wisdom and Holman (1991).
Here the position of each planet, v i for i = 1, . . . , n, is considered relative to the barycentre G i−1 of the previous i bodies, u 0 , . . . , u i−1 , and v 0 is taken as the centre of mass of the system: where η i = i j=0 m j . To have a canonical change of variables the momentaṽ i for i = 0, . . . , n, must be:ṽ In this set of coordinates the Hamiltonian in Eq. 5 takes the form (Laskar 1990b): where Δ i j = ||u i − u j || (the distance between the two bodies) can be expressed as a function of v i and v j , and r i = u i − u 0 . If we fix the centre of mass at the origin then v 0 = 0 and v 0 = 0, and we reduce by 6 the number of equations of motion.

Heliocentric coordinates
Here we consider the relative position of each planet with respect to the Sun: and to have a canonical change of variables the momenta are: In this set of coordinates the Hamiltonian in Eq. 5 takes the form (Laskar 1990b): where Δ i j = ||r i − r j || for i, j > 0. If we consider the centre of mass of the system to be fixed at the origin we have thatr 0 = 0, and we can easily recompute r 0 at all time. Hence, we have also reduced by 6 the number of equations of motion. One of the main differences between these two sets of coordinates is the size of the perturbation which in the case of Jacobi coordinates is smaller than for heliocentric coordinates (see Table 1 in Sect. 4).
Moreover, in the case of Jacobi coordinates the perturbation part (H I ) depends only on positions so it is integrable when we consider it alone. But the expressions are more cumbersome than for heliocentric coordinates (see Appendix B). While in the case of heliocentric coordinates the perturbation part depends on both position and velocities, hence it is not integrable on its own. In Sect. 6 we will show how to adapt the splitting schemes to this particular case.

Test to perform
be a splitting symplectic scheme. We say S(τ ) has s stages if it requires s evaluations of exp(τ A) exp(τ B) per step-size. The smaller the step-size τ used, the smaller is the error of the numerical solution provided by the scheme, and the larger is the computational cost, as more evaluations of exp(τ A) exp(τ B) are required to integrate over the same time period.
Usually, the higher the order of the scheme the more number s of stages it requires, increasing the computational cost to advance a given step-size τ . So a method with 4 stages will be more efficient than one with 2 stages if it can achieve a given accuracy with a stepsize which is at least two times larger than the one required for the 2-stages scheme. In this sense, we define the inverse cost of S(τ ) as τ/s, where s is the number of stages and τ is the step-size used. Thus, if one scheme achieves the same precision than another scheme with smaller inverse cost, then we can say that the former is more efficient than the latter.
It is known that, for sufficiently small step-sizes τ , the method S(τ ) integrates exactly (up to exponentially small errors that are below machine accuracy) a modified Hamiltonian system that is close to the original one. Measuring the maximum variation of the energy along a given orbit will gives us a good idea of how close is that modified Hamiltonian to the original Hamiltonian.
Motivated by that, in all our numerical tests, we measure the relative precision of a given scheme applied with a given step-size τ by computing the maximum variation (E i = max{|H (t 0 ) − H (t)|}) of the energy along a given numerical orbit obtained over 10 5 steps of the method (with the same initial conditions at the initial time t 0 ) and plot E i versus the inverse cost τ/s (both in logarithmic scale). To fix criteria we will always consider step-sizes of the form: τ i = 1/2 i for i = 0, . . . , N .
We are interested in very precise integrations of the Solar System, hence the main goal is to determine for each scheme the maximum step-size (τ i ) required to have an error in the energy variation up to machine accuracy.
Through the paper we consider three test models that we believe illustrate different particularities of the Solar System and can be extrapolated to other planetary systems. These are: a) the motion of the 4 inner planets (Mercury to Mars); b) the motion of the 4 outer planets (Jupiter to Neptune) and c) the motion of the 8 planets on the Solar System (Mercury to Neptune). The initial conditions and mass parameters have been taken from the INPOP10a Solar System ephemerides (Fienga et al. 2011) (www.imcce.fr/inpop/). Table 1 shows estimates on the size of the perturbation for these three examples for both sets of coordinates Jacobi and heliocentric. To estimate the size of the perturbation we have integrated each system over 100 years and computed the maximum values for |H I | and |H K | along this integration. Here HKep represents the size of the Keplerian part and H1max the size of the perturbation part and the estimated size of the perturbation is given by ε = H1max/HKep. We note that all the simulations in this article have been done using an extended real arithmetics and that we use the compensated summation during the intermediate evaluation of exp(a i τ A) and exp(b i τ B) (see Appendix A).

Splitting symplectic integrators for Jacobi coordinates
In Sect. 3 we have seen that with an appropriate change of variables we can rewrite the Hamiltonian of the N-body planetary system as H K + H I where |H I | |H K |. Hence, the system falls into the class of Hamiltonian that can be expressed as with |ε| 1. We can take advantage of this to build efficient high-order splitting symplectic integrators (McLachlan 1995;Laskar and Robutel 2001). In this section we summarise the main ideas behind these methods and review some of the most relevant schemes.
Using the Lie formalism the formal solution of Eq. 12 is: where to simplify notation we use We recall that H A and H B are integrable, hence we can compute explicitly exp(τ A) and exp(τ ε B).
To have a splitting symplectic integrator of order r , we need to find the coefficients a i , b i such that where H is also a Hamiltonian system and belongs to the free Lie algebra generated by A and B, L( A, B). Moreover, we can express H as a double asymptotic series in τ and ε: where p i, j are polynomials in a k and b k .
To have a symplectic scheme of order r we need: p 1,0 = a 1 + a 2 + · · · + a s = 1, , in which case all the even terms in τ in Eq. 15 vanish, having less conditions to satisfy for a scheme of a given order, r , enabling us to find high-order schemes at lower computational cost. There are two different types of symmetric compositions (Eq. 14): one in which the first and last exponentials correspond to the A part (and thus called ABA composition), and the other in which the role of exp(τ A) and exp(ετ B) is interchanged (BAB composition): All the integration schemes that we present in this paper correspond to the ABA class. For the experiments carried out, we have not found substantial differences in the efficiency with respect to methods in the BAB class. Notice that for symmetric methods, the last exponential at one step can be concatenated with the first one at the next integration step when the method is iterated, so the number of exponentials exp(τ A) and exp(ετ B) per step is s, the number of stages.
It is clear that in many cases |ε| τ (or at least ε ≈ τ ). So we can have high-order schemes by only killing the error terms with small powers of ε, and save computational cost by reducing the number of stages of the method.
Depending of the nature of the problem we can try to find the appropriate terms in ε i τ p that must vanish in order to have an optimal performance. For example, if we consider a method such that the coefficients a i , b i satisfy p 1,0 = p 1,1 = 1 and p 2,1 = p 3,1 = p 4,1 = 0, then, but as |ε| τ this method is of effective order 4. In a more general context we will have methods such that, We remark that s 1 is the order of consistency for the scheme, i.e. is the error behaviour in the limit case ε → 0. Nevertheless, in many cases for small step-sizes the method can behave as one of order s 2 . In what follows we will refer to the generalised order of a method in terms of the order in powers of ε. Hence, we will say that a method has order (s 1 , 5.1 ABA schemes of generalised order (2n, 2) McLachlan (1995) noted that as |ε| τ we can have high-order methods by only killing the terms in ετ k . Independently Chambers and Murison (2000), Laskar and Robutel (2001) dealt with this problem following similar ideas, (Laskar and Robutel 2001) providing an explicit computation of the coefficients of the remainder for all order k. One of the main advantages of only killing the terms in ετ k is that we are sure that all the coefficients a i , b i will be positive. As a consequence the coefficients a i , b i will be small and the numerical scheme will be stable.
In Table 2 we summarise the coefficients for the different ABA(2n, 2) schemes for n = 1, . . . , 4. For further details on how to find the a i , b i coefficients and the coefficients for n ≥ 4 see (McLachlan 1995;Laskar and Robutel 2001). Since all the methods we consider are symmetric, we only collect the necessary coefficients of each scheme. Thus, ABA(8, 2) corresponds to the composition We will follow this convention throughout the text.
In Fig. 1 we compare the performance of the ABA(2n, 2) for n = 1, 2, 3, 4 for the 4 inner planets (left) and the 4 outer planets (right). The x-axis corresponds to the cost of the scheme (τ/s) and the y-axis corresponds to the maximum energy variation for one integration at constant step-size τ . Laskar and Robutel (2001) already saw that the optimal schemes for this problem were those of orders (6, 2) and (8, 2) (i.e. SABA 3 and SABA 4 following their notation).
The error on the Hamiltonian approximation of these schemes is O(ετ 2n +ε 2 τ 2 ). In Fig. 1 we can see how the error in energy decreases in τ with slope 2n for large steps-sizes and slope 2 for smaller steps-sizes. We also see how for small step-sizes there is no difference between the cost of the ABA(6, 2) and ABA(8, 2) schemes. In order to improve their performance we need to kill the term in ε 2 τ 2 rather than those of order ετ 2k for k > 4, which are the limiting factor of these schemes. -4 Comparison between the ABA(2n, 2) methods for n = 1, 2, 3, 4 applied to the 4 inner planets (left) and the fou outer planets (right). The x-axis represents the cost (τ/s) and the y-axis is the maximum energy variation over one integration with constant step-size τ . Here s is the number of stages (decimal log scales) 5.2 ABA schemes of order (2n, 4) In this section we will describe three different procedures to cancel the dominant term ε 2 τ 2 in order to get methods of generalized order (2n, 4), and discuss their performance for the different test models described in Sect. 4.  Laskar and Robutel (2001) noticed that it is possible to incorporate this term into the previous compositions with a conveniently chosen constant so as to cancel the term of order ε 2 τ 2 in the asymptotic expansion Eq. 15. We note that this corrector scheme is different than the one introduced by Wisdom et al. (1996) where the corrector added at the beginning and at the end of each step-size is a change of variables. Thus, let S n (τ ) be one of the symplectic ABA schemes of order (2n, 2) described in Sect. 5.1. We can get rid of the term in ε 2 τ 2 by considering

Since in Jacobi coordinates
with the appropriate choice of the constant c. In Table 3 we show the value for the coefficient c for each of the 4 ABA(2n, 2) schemes described before. For further details see Laskar and Robutel (2001). Notice that SC n corresponds to integrating log(S n (τ )) − τ 3 ε 2 cL {{A,B},B} using the leapfrog scheme. So using Eq. 19 with any of the ABA(2n, 2) scheme in Sect. 5.1 we obtain a new integrating scheme of order (2n, 4) with no negative intermediate step. Yoshida (1990) and Suzuki (1990) independently came up with the same idea to find a symmetric scheme of order 2k + 2 from one of order 2k. They both noticed that if S(τ ) is a scheme of order 2k, then:

The composition scheme (S2 m )
is a scheme of order 2k + 2 for an appropriate choice of the constant coefficients x 0 , x 1 . One can check that x 0 , x 1 must satisfy 2x 0 + x 1 = 1 and 2x 2k+1 0 + x 2k+1 1 = 0. Notice that the second condition is used to cancel all the terms of order 2k, while the first one is only for consistency. Laskar and Robutel (2001) used this idea to turn any of the ABA schemes of order (2n, 2) into one of order (2n, 4). If S(τ ) is a symmetric ABA scheme of order (2n, 2) then the composition: is a symmetric method of order (2n, 4) if y 0 , y 1 satisfy 2my 0 + y 1 = 1 and 2my 3 0 + y 3 1 = 0 so, We have done several tests to determine the optimal value of m, and they show that this one is given by m = 2. These results are consistent with those of Suzuki (1990), McLachlan (2002) who did a similar study in a more general framework. The main advantage of this scheme is that we can use it for both heliocentric and Jacobi coordinates.

McLachlan extra stage scheme (ABA84)
McLachlan (1995) studied the possibility of adding an extra stage to the ABA(2n, 2) schemes to get rid of the ε 2 τ 2 term. To add an extra stage results in having an extra pair of coefficients a i , b i and an extra algebraic equation to satisfy. All the coefficients will no longer be positive (Suzuki 1991). In general, if the coefficients are not very large, these methods are stable. The coefficients for the ABA method of generalised order (8, 4) provided by McLachlan (1995) are collected in Table 4.  (1995) id order stg Mer -Ven -Ear -Mar (Jacobi Coord) Jup -Sat -Ura -Nep (Jacobi Coord) Mercury to Neptune (Jacobi Coord) Fig. 2 Comparison between the different schemes to kill the the ε 2 τ 2 terms in the ABA82 scheme. From left to right the 4 inner planets, the 4 outer planets and the whole Solar System. The x-axis represents the cost (τ/s) and the y-axis the maximum energy variation for one integration with constant step-size τ (decimal log scales)

Results
In Fig. 2 we compare the performance of these three different approaches to build methods of generalised order (8, 4) against the ABA(8, 2) scheme (also referred to as ABA82). In the plots we show the cost (τ/s) versus the maximum energy variation for the three test models: the 4 inner planets (left), the 4 outer planets (middle) and the 8 planets in the Solar System (right).
As we can see, the three different schemes improve considerably the performance of the ABA82 (red line). In all cases the corrector scheme SC (blue line) and the McLachlan ABA84 (purple line) show a similar quantitative behaviour. The difference between them is the cost of the extra stage in ABA84, as we are assuming that the corrector is completely free. We note that this is not entirely true if the number of bodies is large (n ≥ 4). On the other hand, the composition methods, S2 m (green line), improves the performance with respect to the ABA82 (red line) but is much more expensive than the other two schemes. 5.3 ABA schemes with generalised order (s 1 , s 2 , . . .) In the previous section we have seen that adding an extra stage to cancel the term of order ε 2 τ 2 gives good results. We can extend this idea and add more stages to kill the error terms accounting to the main limiting factor for each problem (Blanes et al. 2013). This translates into adding more constraints on the coefficients. As long as the increase in the computational cost is less than the gain in accuracy these methods will be competitive. In Fig. 2 we see that the dominant error term for the ABA84 scheme varies between the different test models. Notice that for the outer planetary system the scheme behaves as one of order 4, so the dominant term is ε 2 τ 4 . On the other hand, for the inner planetary system, the method behaves as one of order 8, now the dominant term is ετ 8 .
Hence, to improve the performance of the McLachlan's ABA84 method we need to kill different error terms depending on the problem. For the inner planets, a method of order (10, 4) should perform better than one of order (8, 6). While for the outer planets a method of order (8, 6) should give better results than one of order (10, 4).
In Blanes et al. (2013) we find details on how to solve the algebraic equations and find the set of coefficients a i , b i that provide ABA schemes for a given arbitrary order (s 1 , s 2 , . . .). We must mention that there is no unique combination of coefficients a i , b i for a given order. From all the possible solutions we have selected those that give a better approximation and whose coefficients a i , b i are small. In Table 5 we summarise the coefficients for three ABA schemes: one of order (10, 4); one of order (8, 6, 4) and one of order (10, 6, 4).

Results
In Fig. 3 we compare the performance of the three schemes summarised in Table 5 against the ABA82 and ABA84 for the three test models.
In the left-hand side of Fig. 3 we have the results for the 4 inner planets. We recall that the dominant error term in the ABA84 scheme was ετ 8 . Hence, a method of order 10 in ε should perform better than one of order 8 in ε. Nevertheless, as we can see there is no significant gain in the performance of these schemes with respect to ABA84. Apparently, for these methods the gain in precision is proportional to the computational cost in this range of accuracy.
In the middle of Fig. 3 we have the results for the 4 outer planets. We recall that here the dominant term in the ABA84 scheme was ε 2 τ 4 , hence we expect the schemes of order 6 in ε 2 to be better than the ABA84. As we can see the ABA864 and the ABA1064 schemes give better results that the ABA84. In both cases the optimal cost is around 10 −2 versus an optimal cost of around 10 −3 for the ABA84 scheme.
The main difference between the inner planets and the outer planets is the size of the perturbation. We recall that in Jacobi coordinates, for the inner planets ε ≈ 4.54 · 10 −6 , while for the outer planets ε ≈ 2.03 · 10 −4 (see Table 1). The difference of about 2 orders of magnitude should explain the difference in the performance of the different schemes, as the relevance of the terms ε i τ 2k in the error approximation will vary depending on the size of ε.
Taking this into account, one can be surprised by the performance of these schemes when we consider the whole Solar System (Fig. 3 right). Here the size of the perturbation (ε ≈ 1.96 · 10 −4 ) is of the same order of magnitude as the case of the outer planets. But as we can see in Fig. 3 the schemes behave in the same way as the case of the inner planets, where the terms of order ετ 8 dominate those of order ε 2 τ 4 . We think this is due to Mercury: its fast orbital period and large eccentricity is limiting the performance of the methods. This Table 5 Coefficients for ABA symmetric splitting methods of orders (10, 4), (8,6,4) and ( phenomena was already noticed by Wisdom et al. (1996) and re-discussed by Viswanath (2002). This is why Saha and Tremaine (1994) proposed to use independent time-steps for each planet. They used the leapfrog scheme and adapted it to take fractions of the given step-size for each planet, depending on their orbital period. It is not trivial to extend these ideas using the higher order schemes described in this section, and a second order method is not the appropriate option to achieve round-off accuracy.

Splitting symplectic integrators for Heliocentric coordinates
We recall that all the tests done in Sect. 5 have been done using Jacobi coordinates. All these integrating schemes assume that the two parts of the Hamiltonian H A , H B are integrable. This is true for Jacobi coordinates where: Mer -Ven -Ear -Mar (Jacobi Coord)

Fig. 3
Comparison between the ABA splitting schemes of arbitrary order (s 1 , s 2 , s 3 ) summarised in Table 5 against the ABA82 and ABA84. From left to right the 4 inner planets, the 4 outer planets and the whole Solar System. The x-axis represents the cost (τ/s) and the y-axis the maximum energy variation for one integration with constant step-size τ (decimal log scales) here H K (q, p) is integrable (it is a sum of independent Kepler problems) as well as H I (q) (it only depends on q). However, this is not true for heliocentric coordinates where: and H I (q, p) is not integrable, which can be a problem if we want to apply the splitting schemes presented in Sect. 5. An option to deal with the non-integrability of H I (q, p) is to use another numerical method to integrate H I (q, p) and compute the exp(b i τ B) up to machine accuracy. Unfortunately, this can drastically increase the computational cost of the schemes.
We propose to use the fact that H I (q, p) = T 1 ( p) + U 1 (q) splits naturally into two parts, one depending on positions, the other in velocities. These two parts are integrable on its own and small with respect to H K (q, p). In a general framework, the Hamiltonian splits as: where H A , H B and H C , can be integrated when we consider them separately. In the same way as in Sect. 5, we could try to find appropriate coefficients a i , b i , c i , such that approximates exp(τ L H ). As before, to simplify notation A ≡ {H A , ·}, B ≡ {H B , ·} and C ≡ {H C , ·}. Then one has to deal with the Lie Algebra generated by A, B and C. The number of order conditions as well as the complexity to solve them numerically to get the coefficients a i , b i , c i grows extraordinarily with the order (Koseleff 1993a;Blanes et al. 2013). A simple alternative way to proceed is to use the splitting symplectic schemes in Sect. 5: Here we take C as the Lie operator associated to T 1 ( p) due to its lower computational cost.
In general H B and H C do not commute ({H B , H C } = 0), so this approximation adds an extra error contribution term, ε 3 τ 2 , which will be negligible for small ε. In Fig. 4 we see the result of taking the ABA82, the ABA84 and the S2 m splitting schemes using Eq. 25 to deal with heliocentric coordinates. We can see that in general the symplectic schemes have the same behaviour as with Jacobi coordinates (Fig. 2).
We do see a difference in the case of the outer planets (Fig. 4 middle). Now the ABA84 scheme behaves as one of order 2 for small step-sizes. This is due to the extra error term ε 3 τ 2 . We recall that the main difference between the inner and the outer planets is the size of the perturbation (ε) which is smaller in the first case. Here the terms of order ε 3 τ 2 are negligible for the inner planets but not for the outer planets.
Unfortunately, when we consider high-order symplectic schemes like the ones presented in Sect. 5.3 these extra error terms will become relevant, jeopardising the performance of these schemes.
The first symplectic schemes using heliocentric coordinates (Koseleff 1993a;Touma and Wisdom 1994) used the form of heliocentric Hamiltonian where the unperturbed problem is the Keplerian motion of a planet around the Sun as in Laskar (1991), Laskar (1990b). Later on Duncan et al. (1998) proposed to rewrite the Hamiltonian in heliocentric variables for which the unperturbed problem is the Keplerian motion of a planet around a fixed Sun. For this decomposition, the unperturbed Keplerian approximation is less accurate (see Duncan et al. 1998), but H B and H C commute ({H B , H C } = 0), and thus exp(εb i τ (B + C)) = exp(εb i τ B) exp(εb i τ C). For further details see Appendix C.  We have just seen that for heliocentric coordinates we can adapt the splitting schemes described is Sect. 5 using Eqs. 24 and 25. But with this an extra term ε 3 τ 2 appears in the error approximation that will limit the performance for high-order schemes. One can check that this error term is associated to the algebraic expression b 3 1 + b 3 2 + · · · + b 3 n . We can add an extra stage to the scheme so that it also satisfies: leading to symplectic schemes with the same generalised order as before for heliocentric coordinates. for the outer planets the new ABAH844 scheme behaves better than the ABA84 for small step-sizes.
6.2 ABAH specific methods with arbitrary order (s1, s2, . . .) In the same way we can add the extra constraint b 3 1 + · · · + b 3 n = 0 to the high-order schemes in Sect. 5.3 to have high order splitting schemes specific for heliocentric coordinates. In Table 6 we show the coefficients of two ABAH schemes of orders (8, 6, 4) and (10, 6, 4). All these schemes have one more stage than the schemes presented in Table 5.
In Fig. 6 we compare the performance of the ABA82 with the other three schemes in Table 6. The behaviour of the schemes depending on its order is similar to the one presented in Jacobi coordinates. For the inner planets (Fig. 6 left) all ABAH schemes present a similar optimal cost. For the outer planets (Fig. 6 middle) the ABAH schemes of order (8, 6, 4) and (10, 6, 4) are much better than the other two schemes. We recall that here the size of the perturbation is larger and killing the terms of order ε 3 τ 4 does make a difference. Finally, if we consider the 8 planets in the Solar System (Fig. 6 right) here the ABAH schemes of order (8, 6, 4) and (10, 6, 4) do improve the performance of the schemes of order (8, 4). We recall that this was not the case in Jacobi coordinates (Fig. 3).

Jacobi versus Heliocentric coordinates
In Sects. 5 and 6 we have described different splitting schemes for both Jacobi and Heliocentric sets of coordinates. In the case of heliocentric coordinates the expressions for the Hamiltonian are less cumbersome and easier to handle (see Appendix B). But the size of the perturbation is larger than in Jacobi coordinates and an extra stage to deal with the nonintegrability of H I must be added to have high-order schemes. Here we want to compare the performance of the different schemes for both sets of coordinates.
We compare the performance of the ABA methods of order (8, 2), (8, 4) and (10, 6, 4) in both sets of coordinates, with the three test models used throughout this report. We recall that  Fig. 7 Comparison between Jacobi (continuous lines) and heliocentric (discontinous lines) coordinates using the schemes of orders (8, 2) (red), (8, 4) (blue) and (10, 6, 4) (purple). From left to right the 4 inner planets, the 4 outer planets and the whole Solar System. The x-axis represents the cost (τ/s) of the method and the y-axis the maximum energy variation for one integration with constant step-size τ (decimal log scales) the (8, 4) and (10, 6, 4) schemes have an extra stage in heliocentric coordinates. In Fig. 7 we summarise the performance of these schemes. From left to right we have the results for the inner planets, the outer planets and the whole Solar System. We distinguish the order of the schemes by the colour. The lines in red represent the schemes of order (8, 2), the blue lines those of order (8, 4) and the purple lines those of order (10, 6, 4). We use continuous lines to refer to Jacobi coordinates and discontinuous lines for heliocentric coordinates.
If we look at the results for the inner planets (Fig. 7 left), we can see there is not much difference between taking Jacobi or heliocentric coordinates (continuous vs. discontinuous lines). In both cases the size of the perturbation is small (Table 1) and there is not much difference between taking a splitting scheme of order (8, 4) or (10, 6, 4). In both cases the terms in ε in the error expansion are the ones that matter, but there is not much difference between taking order 8 or ten in ετ k . We should have to use arithmetics with higher precision to see the difference (see Appendix D). Hence the ABA84 is the best choice for this case.
If we look at the results for the outer planets (Fig. 7 middle), again there is no significant difference between Jacobi and heliocentric coordinates. But here the ABA schemes of order (10, 6, 4) perform much better that the other schemes, having an optimal step-size one order of magnitude larger than the ones for the schemes of order (8, 4).
If we look at the whole Solar System (Fig. 7 right), we see that there is a big difference between taking Jacobi or heliocentric coordinates. Looking at the ABA82 scheme (red lines) we see that the slopes are the same but that there is a difference of about one order of magnitude in accuracy for a given step-sizes. If we look at the scheme of order (8, 4) (blue lines), we see that in Jacobi coordinates the methods behave as one of order 8, while in heliocentric coordinates this one behaves as one of order 4. This can be explained by the difference in the size of the perturbation (see Table 1) in both set of coordinates. We also see that there is a big difference between the optimal step-size for both sets of coordinates, making Jacobi coordinates by far the best choice. Finally, if we compare the ABA schemes of order (10, 6, 4) (purple lines) the difference between the two set of coordinates is drastically reduced, although Jacobi coordinates still perform slightly better. However, the extra stages to go from an (8, 4) scheme to a (10, 6, 4) one do not improve in Jacobi coordinates. This is not the case of heliocentric coordinates, where the (10, 6, 4) gives the best results and the difference between the two set of coordinates is not as relevant as before.
To sum up, we recommend to use either the schemes ABA84 or ABA1064 in the case of Jacobi coordinates, while for heliocentric coordinates one should use the ABAH1064 scheme. Although Jacobi coordinates present slightly better results for most of the test models, we believe that using Jacobi or heliocentric coordinates is a matter of choice.

Conclusions
In this article we have reviewed different symplectic splitting schemes and tested their performance for the case of the planetary motion, focussing on the Solar System motion. We recall that in the case of the planetary motion, using an appropriate change of variables, the Hamiltonian of the N-body problem can be rewritten as H K + H I . A sum of independent Keplerian motions for each planet (H K ) and a small perturbation term given by the interaction between the planets (H I ).
There are two set of canonical coordinates that allow us to write the Hamiltonian in this way: Jacobi and heliocentric coordinates (Sect. 3). Although in Jacobi coordinates the size of the perturbation is smaller and H I is integrable, heliocentric coordinates seem more natural and the expressions are easier to handle (Appendix B). In this article we have compared the performance of different symplectic splitting schemes in both sets of coordinates.
In Sect. 5 we described different splitting symplectic schemes for Jacobi coordinates. In Sect. 6 we saw how to extend these schemes to use heliocentric coordinates. We note that all the splitting schemes for Jacobi coordinates can also be used in Heliocentric coordinates, but in order to have a comparable performance an extra stage to kill the terms in ε 3 τ 2 must be added (see Sect. 6).
We have seen that in Jacobi coordinates, the ABA84 scheme introduced by McLachlan (1995) and the ABA1064 scheme (Blanes et al. 2013) give the best results when we look at the motion of the whole Solar System. The high eccentricity of Mercury and its fast orbital period are the main limiting factors and taking higher order splitting schemes do not always provide significant improvements. But for different planetary configurations, as the 4 outer planets, the ABA1064 scheme has a better performance than the ABA84 method.
When we consider heliocentric coordinates, the ABAH1064 method (Blanes et al. 2013) gives the best results when we consider the whole Solar System. In this case, probably because the size of the perturbation is larger, adding extra stages to have higher order schemes does improve the results. The nominal solution La2010a of Laskar et al. (2011a) was computed with the ABA82 scheme and a stepsize of 10 −3 yr. Using the ABAH1064 scheme, the same accuracy should be reached with nearly an order of magnitude improvement in the computing time (Fig. 6). The performances of the schemes in both sets of coordinates, Jacobi or heliocentric are very similar for the scheme of order (10, 6, 4), with a slight advantage for the Jacobi coordinates. Depending on the problem, one can thus use either system of coordinates, but it is clear that using high order schemes as the ABA(H)864 and ABA(H)1064 (Blanes et al. 2013) can drastically improve the results. This should be even more the case for highly perturbed systems as some extra solar planetary systems with close planets of large masses and large eccentricities. It should be nevertheless noted that although these high order integrators behave well until very large eccentricities (up to 0.99 in some of our simulations), some specific integrators need to be used when collisional behaviors are considered (Duncan et al. 1998;Chambers 1999 Fig. 8 Comparison between the ABA schemes of order (2n, 2) for n = 1, 2, 3, 4 applied to the Sun-Jupiter-Saturn three body problem. With (CS) and without (noCS) the compensated summation. The x-axis represent the cost (τ/n) and the y-axis the maximum energy variation for one integration with constant step-size τ . Left using a double precision arithmetics; Right using an extended precision arithmetics (decimal log scales) exp(b i τ B), the increment in position and velocity is done using the compensated summation. In Fig. 8 we show the results for the ABA(2n, 2) schemes for n = 1, 2, 3, 4 using double (left) and extended precision (right). In both cases we gain almost one order of magnitude in precision when we take into account the compensated summation.

Appendix B: Integration schemes (some help on the practical coding)
In this paper we have reviewed many splitting symplectic integrating schemes, all of them of the form: where exp(τ A) and exp(τ B) can be computed explicitly. They correspond to the integrals of the two different parts of the original Hamiltonian. In this section we show how to compute explicitly exp(τ A) and exp(τ B) for the particular case of the N-body problem in Jacobi and heliocentric coordinates. We note that from now on:ũ stands for the momenta associated to u and u stands for du/dt.

B.1 Keplerian motion (H K )
We recall that in Jacobi coordinates, whereas in heliocentric coordinates, In both cases H K is a sum of independent Keplerian motions. In Jacobi coordinates each planet follows an elliptical orbit around the centre on mass of the Sun and the planets that are closer to the Sun, the mass parameter of the system is μ J = Gη i . While in heliocentric coordinates each planet follows an elliptical orbit around the planet-Sun centre of mass, and the mass parameter of the system is μ H = G(m 0 + m i ).
It is well known that Kepler problem is integrable, but the solution from time t = t 0 to t = t 0 + τ is expressed in a simple form if we consider action-angle variables. To compute exp(τ L H K ) we need to be able to compute (r(t 0 + τ ), v(t 0 + τ )) from (r(t 0 ), v(t 0 )).
An option is to change to elliptical coordinates, advance the mean anomaly and then return to Cartesian coordinates. But this can accumulate a lot of numerical errors as well as it is very expensive in terms of computational cost. Instead we use a similar idea as the Gauss f and g functions (Danby 1992), where we use an expression for the increment in position and velocities for a given step-size τ , without having to perform any change of coordinates. Let us give some details on how to derive these expressions.
In elliptical coordinates the motion of the two body problem is given by (a, e, i, Ω, ω, E), where all of the elements remain fixed except for E that varies following Kepler equation (n(t − t p ) = M = E − e sin E). Using a reference frame where the orbital plane is given by Z = 0, the X -axis is the direction of the perihelion and the Y -axis completes an orthogonal reference system on the orbital plane, the position (X, Y, 0) and velocity (X , Y , 0) are given by: where r = a(1−e cos E) and n = μ 1/2 a −3/2 . The position and velocities on a fixed reference frame are given by: where . Given that R 1 (i) = R 1 (i/2)R 1 (i/2) we have that: where p = sin i/2 sin Ω, q = sin i/2 sin Ω, and χ = 1 − p 2 − q 2 = cos i/2. From Eq. 31 we have, Hence, One can check that, where r i = a(1 − e cos E i ) for i = 0, 1. We use Kepler's equation to compute δ E = E 1 − E 0 from δt = t 1 − t 0 . Taking M i = n(t i − t p ) for i = 0, 1, we have that δ E is the solution of x − e cos E 0 sin x − e sin E 0 cos x + e sin E 0 − nδt = 0.

B.2 Jacobi coordinates
We recall that in this set of coordinates the perturbation part is given by:

B.2.1 Computing exp(L H I )
U 1 depends only on the position, hence the equations of motion are given by, As the expressions for ∂U 1 /∂v k can be a little cumbersome, we compute them separately. When we derive H I with respect to v k we must derive 3 main expressions: 1/||v i ||, 1/||r i || and 1/||r i − r j || for i < j. We first give the derivatives of these factors with respect to v k and then we will deduce ∂ H I /∂v k for k = 1, . . . , n.

B.2.2 Computing the corrector: exp(L {{A,B},B} )
In Sect. 5.2.1 we described a splitting symplectic schemed where a corrector term was added at the beginning and at the end of each step-size. The corrector term is given by, with L C = L {{A,B},B} and c a constant coefficient that depends on the order of the ABA scheme.
In Jacobi coordinates A is quadratic in p and B only depends on q so {{A, B}, B} only depends on q and {{A, B}, B} is integrable. We recall that A = H K ep = T 0 + U 0 and B = H pert = U 1 . Hence, Then the equations of motion for L C are given by: where γ k = η k η k−1 m k . As before, usingṽ i = Again the expressions for ∂ 2 U 1 ∂v i ∂v k are a little cumbersome and we first show how to derive the different parts in ∂U 1 ∂v k : v i /||v i || 3 , r i /||r i || 3 and r i − r j /||r i − r j || 3 .
From now on we call Acc(i) = γ i ∂U 1 ∂v i , and We can now give the expressions for

B.3 Heliocentric coordinates
We recall that in this set of coordinates the perturbation part is given by:

B.3.1 Computing exp(τ L T 1 )
Notice that T 1 depends only on the momenta (r). Hence, the equations of motion are given by, Finally, Notice that U 1 depends only on the positions (r). Hence, the equations of motion are given by, Given thatr k = m kṙk , we have:

Appendix C: Heliocentric coordinates (alternatives for the set of equations)
The canonical heliocentric (CH) coordinates used in Sect. 6 are canonical and the position of each body is taken with respect to the position of the Sun. The position and their associated momenta are given by: The main difference between Jacobi and heliocentric coordinates is that in the second set of coordinates the kinetic energy is not diagonal in the momenta. Instead we have: which can be rewritten as: The extra term due to the momenta of the Sun is added to the perturbation part and makes it depend on both position and velocities. In Sect. 3 we used Eq. 43 to derive the Hamiltonian expression, as in Laskar (1990bLaskar ( , 1991, Koseleff (1993a). Duncan et al. (1998) used Eq. 42 instead, and a Keplerian approximation in which the planets orbits around a fixed Sun (Eq. 48). Following Laskar (1991) we will call the first set of equations CH, and following Duncan et al. (1998) we will call the second set of equations democratic canonical heliocentric (DCH). Here we discuss the main differences between the two sets of equations and compare the performance of the integrators presented in this paper for both expressions.

Fig. 9
Comparison between the two expressions for heliocentric coordinates: the classical CH expressions (Eqs. 44-46) continuous lines and the DCH expressions (Eqs. 47-49) discontinuous lines. For the schemes ABA82 (red), ABA84 (green) and ABA1064 (blue). From left to right the 4 inner planets, the 4 outer planets and the whole Solar System. The x-axis represents the cost (τ/s) of the method and the y-axis the maximum energy variation for one integration with constant step-size τ (decimal log scales)

C.2 Comparisons between the expressions
We recall that when we use DCH variables (Eqs. 47-49) we use the splitting schemes discussed in Sect. 5 with While when we use the classical CH variables (Eqs. 44-46) we use the splitting schemes discussed in Sect. 6 with We compare the ABA schemes of orders (8, 2), (8, 4) and (10, 6, 4) for both splitting expressions. We recall that the schemes of order (8, 4) and (10, 6, 4) that use CH variables (Eqs. 44-46) have one more stage than the schemes used with DCH variables (Eqs. 47-49). Figure 9 summarise the performance of the different integrating schemes presented in Sects. 5 and 6. From left to right we have the results for the inner planets, the outer planets and the whole Solar System. The red lines show the performance of the ABA82 scheme, the green lines are for the ABA84 schemes and the blue lines are for the ABA1064 schemes. We use continuous lines when we consider CH variables and discontinuous lines for DCH variables.
As we can see, there is no significant difference between using one splitting or the other. In some cases one is better that the other. The main advantage that the splitting DCH introduced by Duncan et al. (1998) is that we do not require an extra stage for a high-order scheme. Mercury to Neptune (Helio Coord)

Fig. 11
Comparison using heliocentric coordinates between the ABAH82, ABAH84, ABAH864 and ABAH1064 schemes using quadruple precision arithmetics. From left to right the 4 inner planets, the 4 outer planets and the whole Solar System. The x-axis represents the cost (τ/s) and the y-axis the maximum energy variation for one integration with constant step-size τ (decimal log scales)

Appendix D: Comparison in quadruple precision
As we have discussed throughout the article, in many cases we have seen that despite taking higher order methods no significant improvement on the performance of the schemes was observed. This is the case of the 4 inner planets in the Solar System, where the size of the perturbation is so small that the extra stages to increase the order of the schemes are useless.
Here the round-off error dominates the terms in ε 2 and ε 4 . Similar results are also observed when we consider the whole Solar System. In order to see an improvement we need to use higher precision arithmetics. Here we have repeated the test from Sects. 5 and 6 for the different integrating schemes using quadruple precision arithmetics. We want to illustrate that the different schemes of orders (8, 6, 4) and (10, 6, 4) perform better that those of order (8, 4). In Figs. 10 and 11 we show the results for the same test models used throughout the article for Jacobi and heliocentric coordinates, respectively using quadruple precision arithmetics. For Jacobi coordinates (Fig. 10) we compare the ABA82, ABA84, ABA104, ABA864 and ABA1064 schemes. For heliocentric coordinates (Fig. 11) we compare the ABAH82, ABAH84, ABAH864 and ABAH1064. As we can see in Fig. 10 for Jacobi coordinates, the ABA864 and ABA1064 (Blanes et al. 2013) do improve the performance of the ABA84 (McLachlan 1995). Notice also that for the 4 inner planets (Fig. 10 left) and the whole Solar System (Fig. 10 right) the improvement is achieved for small step-sizes, where the energy variation is below the machines epsilon for extended arithmetics precision. In Fig. 11 similar results are observed for heliocentric coordinates.
From these experiments we see how the ABA splitting methods of orders (8, 6, 4) and (10, 6, 4) for both sets of coordinates improve the performance of the McLachlan (1995) ABA84 and the Laskar and Robutel (2001) ABA82.