New families of symplectic splitting methods for numerical integration in dynamical astronomy

We present new splitting methods designed for the numerical integration of near-integrable Hamiltonian systems, and in particular for planetary N-body problems, when one is interested in very accurate results over a large time span. We derive in a systematic way an independent set of necessary and sufficient conditions to be satisfied by the coefficients of splitting methods to achieve a prescribed order of accuracy. Splitting methods satisfying such (generalized) order conditions are appropriate in particular for the numerical simulation of the Solar System described in Jacobi coordinates. We show that, when using Poincar\'e Heliocentric coordinates, the same order of accuracy may be obtained by imposing an additional polynomial equation on the coefficients of the splitting method. We construct several splitting methods appropriate for each of the two sets of coordinates by solving the corresponding systems of polynomial equations and finding the optimal solutions. The experiments reported here indicate that the efficiency of our new schemes is clearly superior to previous integrators when high accuracy is required.


Introduction
Symplectic integrators have several features that turn out to be particularly appropriate when integrating numerically for long times evolution problems in dynamical astronomy. They preserve by construction the symplectic structure of the original Hamiltonian problem, so that the numerical solution inherits the qualitative properties of the exact one [20]. In particular, by using backward error analysis, it is possible to prove that this numerical solution is in fact exponentially close to the exact solution of a modified Hamiltonian. Moreover, although the energy is not conserved along the trajectory, the error introduced by a symplectic method of order r used with constant step size τ is of order O(τ r ) for exponentially long time intervals under rather general assumptions, whereas the error in position typically grows linearly with time [7]. Assume that, as is often the case, the Hamiltonian function is of the form H(q, p) = T (p) + U (q), where the potential energy U (q) depends on positions and the kinetic energy T (p) is a function of the conjugate momenta. Then the equations of motion corresponding to T (p) are trivially solvable, and the same happens with U (q). By composing the flows of these two special Hamiltonian systems one gets a symplectic first order approximation to the exact flow. This simple composition constitutes an example of a symplectic splitting method. Higher order approximations can be obtained by composing the flows corresponding to T (p) and U (q) with certain coefficients obtained by solving the so-called order conditions [16]. There exist in the literature a vast number of high order integrators constructed along this line (see, e.g., [1], [7], and references therein).
The non-relativistic gravitational N-body problem, in particular, belongs to this class of systems. If one considers the motion of n + 1 particles (the Sun, with mass m 0 , and n planets with masses m i , i = 1, . . . , n) only affected by their mutual gravitational interaction, the corresponding equations of motion can be derived from the Hamiltonian where q i and p i = m iqi denote the position and momenta of the n+1 bodies in a barycentric reference frame. Typically, the planets evolve around the central mass following almost Keplerian orbits, so that by an appropriate change of coordinates one can rewrite the Hamiltonian (1) as H = H K + H I , where in some sense |H I | |H K |, or equivalently, as the sum of the Keplerian motion of each planet around the central mass and a small perturbation due to the gravitational interaction between planets. Jacobi and Heliocentric coordinates constitute paradigmatic examples of canonical set of coordinates possessing this feature. Thus, the Hamiltonian (1) written as H = H K + H I is a particular example of a near-integrable Hamiltonian system, i.e, it can be expressed as H(q, p; ε) = H [a] (q, p) + εH [b] (q, p), where ε 1 and H [a] is exactly integrable. It makes sense, then, to take into account this special structure when designing integration methods to ap-proximate its dynamics. The idea consists in constructing splitting schemes as compositions of the flows corresponding to H [a] (q, p) and H [b] (q, p), assuming that they are explicitly computable or sufficiently well approximated [8,26]. In fact, since the parameter ε is small, it is possible to design methods which behave in practice as high order integrators with less severe restrictions concerning the order conditions than the usual split into kinetic and potential energy. This approach was systematically pursued by McLachlan [14], obtaining families of splitting schemes of order 2 and 4 which eliminate the most relevant error terms in ε, and further analyzed by Laskar & Robutel [12] in the context of planetary motion.
By incorporating the idea of processing, even more efficient schemes can be constructed for the Hamiltonian (2) [2]. In that case, both the kernel and the processor are taken as compositions of the flows associated with H [a] and H [b] , so that the exactly symplectic character of the integration scheme is ensured. With this approach, all terms of first order in ε in the truncation error expansion can be annihilated with the processor [15,27].
Although the symplectic methods developed in [12] and [14] for near-integrable Hamiltonian systems have proved their usefulness in long term integrations of the Solar System [11], the design of new and more efficient higher order integrators is of interest for numerical simulations of its evolution over large time spans, either by speeding up the algorithms or by providing better accuracy in the position of the different objects. Relevant examples where the new integrators could be useful include the numerical integration of the Solar System for more than 60 million years backward in time to cover the Palaeogene period to determine insolation quantities of the Earth and calibrate paleoclimatic data, studies of the planetary orbits over several billion years, etc. [13]. To this purpose, it is essential that the numerical solutions obtained are not contaminated by error accumulations along the integration and that the computations are done in a reasonable time.
These long-time numerical integrations can be combined with standard techniques of classical perturbation theory, such as the expansion of the equations of motion up to a certain order in the perturbation parameters and the use of averaging (see, e.g. [9]).
The purpose of this work is to present new families of symplectic splitting methods specifically designed for Hamiltonian systems of the form (2) appearing in many problems of dynamical astronomy, when one is interested in highly accurate results over a large time span. The schemes we propose will be useful in particular in the long time integration of the Solar System, both in Jacobi and Poincaré Heliocentric coordinates, and are more efficient than the schemes designed in [12] and [14]. Although they involve the computation of more elementary flows per step than other methods, their small error terms allow to use larger steps, which results in more efficient schemes. Obtaining these new methods requires deriving previously the necessary and sufficient order conditions to be satisfied by the coefficients (which is done here in a systematic way) and then solving these polynomial equations to get the best solutions according with some appropriately chosen optimization criteria. This is discussed in more detail in sections 2, 3 and the appendix, whereas in section 4 we consider the application of the new schemes to the integration of the Solar System. The new methods obtained in section 3 are suitable to be applied when using Jacobi coordinates and also Poincaré Heliocentric coordinates.
It is worth stressing that, while the main motivation of this work is the long time integration of Hamiltonian problems arising in dynamical astronomy, and in particular in planetary systems, the new symplectic splitting methods obtained here can also be applied to more general perturbed differential equations arising in different fields when high accuracies are required.

Order conditions 2.1 Preliminaries
To establish the framework for the construction and analysis of the new families of integrators, we consider a generic differential equation of the form where |ε| 1 and each part is exactly solvable (or can be numerically solved up to round off accuracy) with solutions τ (x 0 ) respectively, at t = τ , the time step. If we denote by ϕ τ (x 0 ) the exact solution of (3), it is well known that ψ τ = ϕ τ provides a first-order approximation, i.e., ψ τ (x 0 ) = ϕ τ (x 0 ) + O(τ 2 ) and that higher order approximations can be obtained by taking more compositions in ψ τ , for appropriately chosen coefficients a i , b i . The splitting method ψ τ is said to be of order r if for all x ∈ R D , It is straightforward to check that the method is at least of order 1 for arbitrary problems of the form (3) if and only if the coefficients a i , b i satisfy the consistency condition We are mainly interested in symmetric methods, that is, integrators verifying ψ −τ = ψ −1 τ , or equivalently a s+2−i = a i , b s+1−i = b i (so that the composition (5) is left-right palindromic). In that case, they are automatically of even order. In particular, if a symmetric method satisfies the consistency condition (7), then it is at least of order 2.
Since the last flow ϕ τ per step is precisely s. This number is usually referred to as the number of stages in the composition.

Deriving the order conditions via the BCH formula
The conditions that the coefficients a i , b i must satisfy for a splitting method to be of order r (the so-called order conditions) can be conveniently derived by considering series of linear differential operators. We denote by A and B the Lie operators associated with f [a] and f [b] , respectively. For each smooth function g : R D → R, A g and B g are smooth functions defined as for each x ∈ R D , that is, The near-integrable Hamiltonian system (2) corresponds in this general framework to considering equation (3) with being J the canonical symplectic matrix. Therefore, for each smooth function g one has and a similar expression for B g.
It is well known that for any smooth function g, the τ -flow of (3) satisfies where e τ (A+εB) is defined as a series of linear differential operators The same is true for each part in (3): Analogously, for the integrator ψ τ in (5), one has where Ψ(τ ) is a series of linear differential operators defined as Ψ(τ ) = e a 1 τ A e b 1 τ εB · · · e asτ A e bsτ εB e a s+1 τ A .
Notice that the exponentials of Lie derivatives in (10) appear in the reverse order with respect to the maps in the integrator (5).
One of the standard ways of deriving the order conditions for splitting methods is the following. By applying repeatedly the Baker-Campbell-Hausdorff (BCH) formula [24] to the factorization (10) corresponding to a consistent (i.e., satisfying (7)) splitting method, one is able to express Ψ(τ ) as the formal exponential of only one operator: where Here the symbol [A, B] stands for the commutator of the Lie operators A and B, and p ab , p abb , p aba , p abbb , . . . are polynomials in the parameters a i , b i of the splitting scheme. In particular, and c s+1 = 1. The integrator is of order r if E(τ, ε) in (12) is of size O(τ r ), so that Ψ(τ ) agrees with the series of linear operators e τ (A+εB) of the exact flow up to terms of size O(τ r ). In consequence, the order conditions read p ab = p abb = p aba = · · · = 0 up to the order considered. For symmetric methods, Ψ −τ = Ψ −1 τ , and thus E(−τ, ε) = E(τ, ε), so that E(τ, ε) only involves even powers of τ , that is, p w = 0 for any word w with an even number of letters in the alphabet {a, b}.
In (12) we have considered the classical Hall basis associated to the Hall words a, b, ab, abb, aba, abbb, abba, abaa, . . . [19]. The coefficients p w in (12) corresponding to each Hall word w can be systematically obtained using the results in [18] in terms of rooted trees and iterated integrals. An efficient algorithm (based on the results in [18]) of the BCH formula and related calculations that allows one to obtain expression (12) up to terms of arbitrarily high degree is presented in [3].

Generalized order
We are particularly interested in the manner in which the local error ψ τ (x) − ϕ τ (x) decreases as ε → 0. For instance, from the results in the precedent subsection, it is clear that for any consistent symmetric method the local error If in addition p aba = 0 in (12), then In that case, we say that such a method is of (generalized) order (4,2). More generally, we will say [14] that an integration method for the system (3) is of generalized order (r 1 , r 2 , . . . , r m ) (where r 1 ≥ r 2 ≥ · · · ≥ r m ) if the local error satisfies that Recall from Subsection 2.2 that for symmetric integrators, the remainder E(τ, ε) in (11) is even with respect to τ , and thus the generalized order (r 1 , r 2 , . . . , r m ) of symmetric schemes must have even r j .

Generalized order conditions
The conditions that the coefficients a i , b i must satisfy for a splitting method to be of a prescribed (generalized) order (r 1 , r 2 , . . . , r m ) can be obtained, of course, by computing the polynomials p w in expression (12) with the BCH formula and then equating each term to zero up to the considered order. Thus, in particular, a consistent symmetric scheme of order (6, 2) requires that p aba = p abaaa = 0.
There exist, however, other more systematic procedures to derive these order conditions. In what follows, we present a strategy that allows us to get in a direct way a set of necessary and sufficient independent order conditions for generic splitting methods.
As a first step, we consider Z(τ ) = e τ (A+ε B) e −τ A , which is the formal solution of the initial value problem where with On the other hand, applying repeatedly the identity e τ A e h B e −τ A = e h C(τ ) to eq. (10) and taking into account (13), we arrive at Notice that, if the splitting method is consistent, then c s+1 = 1. We thus have that a splitting method is of order (r 1 , r 2 , . . . , r m ) if and only if c s+1 = 1 and We then expand both Z(τ ) andZ(τ ) as power series of ε and compare their coefficients. First, applying Neumann iteration to (14) we get where in the last equality we have introduced explicitly the expression for C(τ ) given by (15). On the other hand, by expanding the exponentials ofZ(τ ) in (16), we have In this way a splitting method is of order (r 1 , . . . , r m ) if and only if for each k = 1, . . . , m and each multi-index (i.e., k-tuple of positive integers) (j 1 , . . . , j k ) such that j 1 + · · · + j k ≤ r k . Conditions (18) (one condition for each multi-index) have been obtained in [23] in the context of order conditions of splitting operators for unbounded operators A and B. Nevertheless, such order conditions are not all independent. For instance, it can be checked that if condition (18) holds for the multi-indices (1, 2), (2), and (1), then the condition for (2, 1) is also fulfilled. That kind of dependencies are a consequence of the fact that both Z(τ ) and Z(τ ) are exponentials of Lie series in the non-commuting indeterminates C 1 , C 2 , . . .. A set of independent order conditions can be obtained (by virtue of Theorems 3.2 and 6.1 in [19]) by considering a particular subset of multi-indices, the so-called Lyndon multi-indices. Let us consider the lexicographical order < (i.e., the order used when ordering words in the dictionary) on the set of multi-indices.
Taking into account these considerations, we finally arrive at the following result.
For illustration, in Table 1 we collect explicitly conditions (18) corresponding to some particular multi-indices, whereas in Table 2 we specify which particular Lyndon multi-indices one has to consider, or equivalently which conditions (18) must hold for each consistent symmetric splitting method of the given generalized order, according to Theorem 1. At this point some remarks must be done. The set of order conditions given by Theorem 1 is completely equivalent to the order conditions that can be obtained by following the standard approach described in subsection 2.2. On the one hand, the free Lie algebra L(C 1 , C 2 , C 3 , . . .) generated by the noncommuting indeterminates C 1 , C 2 , C 3 , . . ., admits a basis (the Lyndon basis [19]) in one-to-one correspondence with the set of Lyndon multi-indices. Clearly, if instead of directly comparing the series expansions of Z(τ ) and Z(τ ) as above, we compare the formal logarithms log(Z(τ )) and log( Z(τ )), we could obtain one order condition per element in the Lyndon basis. On the other hand, the  approach in subsection 2.2 gives one order condition per element in a basis of the free Lie algebra L(A, B) generated by the noncommuting indeterminates A and B, and Lazard elimination theorem [19] shows that, as vector spaces, the direct sum of L(C 1 , C 2 , C 3 , . . .) with the linear span of A is isomorphic to

New numerical schemes
There are two different types of symmetric composition schemes (5): one in which the first and last flows correspond to the A part (and thus appropriately called ABA composition), and the other in which the role of ϕ τ is interchanged (BAB composition): BAB: ϕ Notice that both types of composition are closely related: an s-stage BAB method is just an (s + 1)-stage ABA scheme with a 1 = 0, so that to construct BAB methods one has to solve the same order conditions as for ABA compositions. Although A and B are qualitatively different here, and therefore both types of composition may lead in principle to integrators with different performances, in practice, and for the examples analyzed, we have not found substantial differences, so that in what follows we only consider ABA methods for clarity in the presentation. Constructing particular methods requires solving polynomial equations (e.g. the order conditions of Table 1 for methods of Table 2), a problem whose complexity grows enormously with the number of equations and variables involved. This task can be handled by computer algebra systems when this number is relatively low. In that case one is able to get all the solutions and select the one that verifies some previously fixed optimization criterion, such as minimizing error terms at higher orders and the sum of the absolute value of the coefficients. In practice, we have followed this procedure when there are no free parameters and the number of equations to be solved is at most seven. When this number is larger than seven or there are additional parameters, another strategy based on homotopy continuation methods has been applied. In the appendix we provide a detailed treatment of the procedure for a particular method.

New methods in the ABA class
Symmetric schemes of generalized order (2n, 2) can be obtained just by solving, in addition to consistency, the order conditions corresponding to the Lyndon multi-indices (3), (5), . . . , (2n − 1) (first line in Table 1). These equations result from approximating the integral τ 0 C(s)ds in the expression of Z(τ ) by the quadrature rule in the expansion ofZ(τ ). Equivalent order conditions were previously derived in [12,14,15], and so the same methods are obtained here. Methods in this family have all their coefficient positive and good stability properties. In the tests carried out in this paper we will take the most efficient ABA scheme of order (8,2) for comparison, which we denote by ABA82.
Generalized order (10,4). According with Table 2, there are five order conditions in addition to consistency (7), for a total number of seven equations to be satisfied by the coefficients. As a consequence, the minimum number of stages is six. A more efficient method can be obtained, however, by taking an additional stage and choosing the corresponding free parameter to reduce the error terms at a higher order. The sequence of coefficients is and their values are collected in Table 3 (method denoted by ABA104). Observe that, as expected, one of the a i and one of the b j coefficients are negative (it is known that this feature is unavoidable for any splitting method of order higher than two [6,21,22]), but they have a relatively small absolute value.
Generalized order (8,6,4). Here we have, in addition to consistency, six order conditions, for a total of eight equations, so that the minimum number of stages is seven. Hence the sequence of coefficients for the resulting methods is also as in (21). There are 30 real solutions, and the one referred to as method ABA864 in Table 3 minimizes the sum of the absolute values of its coefficients.
Generalized order (10,6,4). An additional stage is required in this case to verify the order condition associated with multi-index (9) in Table 1. Therefore, the minimum number of stages is eight, with sequence By following the construction strategy exposed in the appendix, we have obtained several solutions for the order conditions. Among those possessing reasonably small coefficients, we have selected the solution with the smallest leading terms of the local error. This corresponds to method ABA1064 in Table 3.  Table 3: Coefficients for ABA symmetric splitting methods of generalized order (10,4), (8,6,4) and (10, 6, 4).

A simple example
To illustrate the efficiency of these schemes we take the perturbed Kepler problem with Hamiltonian where r = q 2 1 + q 2 2 . This Hamiltonian is a first approximation used to describe the dynamics of a satellite moving into the gravitational field produced by a slightly oblate spherical planet and whose motion takes place in a plane containing the symmetry axis of the planet [17].
We consider this simple problem to test the relative performance of the methods obtained in this work in comparison with schemes presented in [12,14]. The following schemes are used: • ABA82: The 4-stage (8,2) ABA method given in [12,14].
For the numerical experiments we take as initial conditions q 1 = 1 − e, q 2 = 0, p 1 = 0, p 1 = (1 + e)/(1 − e), with e = 1/4, which would correspond to the eccentricity of the unperturbed Kepler problem. For this system, the strength of the perturbation depends both on the choice of the small parameter, ε, and the initial conditions. We integrate along the interval t ∈ [0, 10000] and compute the averaged error in energy as well as the averaged error in position and momenta (measured in the 2-norm) of the numerical solutions evaluated at t k = 20 · k, k = 1, 2, . . . , 500. We take as the exact solution an accurate approximation obtained using a high order method with a sufficiently small time step. This numerical test is repeated several times for each method using different time steps (changing the computational cost for the numerical integration). Finally, we plot the average errors versus the time step scaled by the number of stages per step, i.e. τ /s, in double logarithmic scale, to show how the error depends on the computational cost (the cost is inversely proportional to τ /s, and the best methods should provide a given accuracy with the largest value of τ /s). Figure 1 shows the results obtained for ε = 10 −2 , 10 −3 . In diagrams (a) and (c) we show the average error in positions and momenta, whereas in pictures (b) and (d) we measure the average error in energy. Notice how the new methods collected in Table 3 are clearly more efficient than ABA82 and ABA84.
It should be stressed that, although the coefficients in Table 3 have 40 digits of accuracy, the results displayed in Figure 1 have been obtained, for the sake of illustration, with a standard Fortran compiler in double precision. The code generating the results for the averaged error in energy is available at the website www.gicas.uji.es/software.html.

Splitting methods with approximate flows
We have so far assumed that the exact τ -flow maps ϕ τ will no longer be available. In that case, instead of the splitting method (5), we will consider a composition of the form where ϕ τ is an approximation of ϕ τ obtained by applying some numerical integrator to the Hamiltonian εH [b] (q, p).
In what follows, we assume that ϕ τ represents one step of some 2nd-order symmetric method. In that case, the series of differential operators corresponding to ϕ instead of just e τ εB , and thus, the series Ψ h of differential operators corresponding to the method (24) can be obtained from (10) by replacing each e b j τ εB by In consequence, the leading term of the difference Ψ τ − Ψ τ of the respective series corresponding to methods (24) and (5) is It is then natural to impose, in addition to the generalized order conditions obtained in subsection 2.3, the condition with the aim of reducing the effect of replacing ϕ (5).
The order conditions of scheme (24) τ being one step of arbitrary second order symmetric integrator applied to y = εf [b] (y)) can be systematically obtained by generalizing the approach presented in section 2, just by replacing Z(τ ) in (16) by Thus, the scheme (24) will have generalized order (r 1 , r 2 , r 3 , . . . , r m ) if (5) is of generalized order (r 1 , r 2 , r 3 , . . . , r m ) (see subsection 2.3) and in addition, In particular, we have that This shows that, if the method (24) is applied with the coefficients of a standard symmetric splitting method (5) of generalized order (r, 4), the resulting method has generalized order (r, 4, 2). If the additional condition (25) holds, then (24) recovers the generalized order (r, 4) (recall that generalized order (r 1 , r 2 , . . . , r m ) of symmetric methods have even r j ). We have constructed several symmetric methods of ABA-type a 1 τ (with a 1 = 0) satisfying the additional condition (25). Notice that ABA-type compositions are more convenient than BAB-type methods, since the last stage of the method in the current step can be concatenated with the first stage at the next step. This is not possible with BAB compositions (20), because By following a similar strategy as for the methods collected in section 3, we have constructed symmetric schemes within this family of generalized order (8,4), (8,6,4) and (10,6,4), all of them involving the minimum number of stages.
The corresponding coefficients are collected in Table 4. For schemes (8,4) we have found all the real solutions and selected the solution that minimizes the sum of the absolute values of the coefficients (method ABAH844). This method has the following structure The procedure for constructing method ABAH1064 is detailed in the appendix, and a similar strategy has been used to build scheme ABAH864.  Table 4: Coefficients for ABA symmetric splitting methods of generalized order (8,4), (8,6,4) and (10,6,4) especially adapted to be used when the flow ϕ τ is approximated by a symmetric 2nd-order method ϕ τ . This happens, in particular, when Heliocentric coordinates are used for the integration of the Solar System, as it is shown in section 4.

Application to the integration of the Solar System
In this section we illustrate how the new families of methods proposed here behave when they are used in the numerical integration of the simplest model of the Solar System, i.e., a main massive body (the Sun) and a set of particles (the planets) orbiting the Sun following almost Keplerian trajectories. It is not our intention to carry out a detailed treatment of this problem, but rather to check the performance of the new methods and compare them with other well established schemes designed for near-integrable Hamiltonian systems such as those presented in [12] and [14]. We instead refer the reader to reference [4], where this issue is handled in much more detail.
As stated in the Introduction, integrating numerically the gravitational Nbody problem requires first to choose a convenient set of canonical coordinates. Two widely used coordinate systems where the corresponding Hamiltonian (1) adopts the form (2), suitable to the application of the integration schemes developed in this work, are Jacobi and Heliocentric coordinates. In the former, the position of each planet is taken relative to the barycenter of the previous i bodies, whereas in the later the position of each planet is taken with respect to the Sun.
Methods of Table 3 are particularly appropriate for long time integrations of the N-body problem in Jacobi coordinates, and extensive numerical experiments with different planetary configurations have been carried out in [4]. Here we will restrict ourselves to Heliocentric coordinates.
In this set the coordinates r i are the relative positions of each planet with respect to the Sun: whereas the conjugate momenta read r 0 = p 0 + · · · + p n ,r i = p i (27) and the Hamiltonian (1) is given by [10] where ∆ ij = r i − r j for i, j > 0. Heliocentric coordinates have the advantage, compared with Jacobi coordinates, that adding a new body to the model does not change the origin and the new Hamiltonian is easily updated. On the other hand, the perturbation H [b] depends on both positions and momenta and is not integrable by itself. Hence, when considering the splitting (3), the equation τ is approximated by a 2nd-order method, we can use the new splitting methods of Table 4. Notice that H [b] is the sum of two terms, one depending only on positions and the other depending only on momenta τ . Therefore, we can approximate ϕ τ by the second order symmetric scheme ϕ Another possibility consists in taking splitting methods of the form a 1 τ and obtaining the appropriate coefficients a i , b i , c i satisfying the required order conditions. These can be derived, for instance, by analyzing the free Lie algebra generated by the three Lie derivatives corresponding to each piece of the Hamiltonian. The number of order conditions grows rapidly with the order, however, in comparison with splitting schemes involving only two parts. In particular, (10,6,4) methods require, in addition to consistency, 22 order conditions, while methods of the form (24) with (29) only need to satisfy 8 order conditions. This being the case, in what follows we consider splitting methods of the form (24), with the approximate flow ϕ τ given by the leapfrog composition (29). In Figure 2 we can see the results achieved by the new ABA splitting schemes of Table 4 on the N-body problem in Heliocentric coordinates for different planetary configurations. Method ABA82 refers here to the composition (24) with the coefficients of the ABA scheme of generalized order (8,2) given in [12,14], but with the flow ϕ We have integrated the same initial conditions for each scheme using different step sizes τ . For each step size (τ i = 1/2 i for i = 1, 15) we have computed the numerical trajectory over niter = 10 5 evaluations of the integration scheme (i.e. if τ = 0.5, then the final time is t f = 50000 years). For each trajectory we plot the maximum variation in energy along the trajectory versus the inverse of the computational cost, τ /s, both in logarithmic scale.
All the simulations have been done in Fortran using extended double precision and compensated summation during the evaluation of the inner stages of each scheme.
Notice that in the case of the four inner planets (Figure 2 (a)) the performance for the different ABA schemes in Table 4 are better than method ABA82, but there is not much difference between the ABA schemes ABAH844, ABAH864 and ABAH1064. Nevertheless, if we look at the results for the four outer planets and the whole Solar System (Figure 2 (b) and (c) respectively) we can see that ABAH864 and ABAH1064 show better results than ABAH844.
Despite the fact that the size of the perturbation for the different planetary configurations presented can be very different [4], we believe that the difference between the performance of the schemes for the inner and the outer planets in the Solar System is mainly due to Mercury. Its fast orbital period and relatively high eccentricity are the main limiting factors when one tries to improve the efficiency of the higher order schemes. Notice that for the inner planets the orbital period is much shorter than for the outer planets, and so, according with [25], this imposes a restriction on the step size to be used and the number of evaluations per orbital period.

Concluding remarks
In reference [12], symplectic splitting methods of generalized order (2n, 2) up to n = 5 (first described in [14]) were systematically derived and tested on the Sun-Jupiter-Saturn system over 25000 years in Jacobi coordinates, observing an improvement in the accuracy with respect to the leapfrog integrator by several orders of magnitude at the same computational cost. Methods in this family have all the coefficients positive and good stability properties. Scheme (8,2) in particular has been used in several long term simulations of the whole Solar System (e.g., [11,13]) and corresponds to method ABA82 in the examples reported here. All the tests carried out in [12] showed that the error term ε 2 τ 3 was the main limiting factor in the performance of the integrators, so the natural question was whether schemes of higher order (and thus already involving some negative coefficients) could be useful for integrating planetary N-body problems.
As a matter of fact, methods of order (8,4) obtained in [14] do improve the performance of ABA82 for this problem in Jacobi coordinates, as the experiments reported in [4] show. In this work we have pursued this line of research and constructed new families of higher order splitting methods specifically oriented to the numerical integration of near-integrable Hamiltonian systems, and in particular for planetary N-body problems, both in Jacobi and Poincaré Heliocentric coordinates. For this purpose, first we have derived explicitly the set of independent necessary and sufficient order conditions that splitting methods must verify to achieve a certain order of accuracy and then we have solved these equations. A non-trivial task that requires the use of homotopy continuation techniques and optimization criteria to select the most appropriate solution.
Although the new methods involve some negative coefficients, and thus one could think that their numerical stability might be compromised, they have been selected to minimize the error terms at higher orders and the sum of the absolute values of their coefficients. As a result, the size of the negative coefficients of our new methods is relatively small. In any case, the experiments reported here clearly indicate that the new methods of order (8,6,4) and (10,6,4) achieve accuracy up to round off error with larger step sizes than 2nd-order schemes.
There are near-integrable systems of the form (3) where the exact flow ϕ corresponding to the perturbation is not available. In that case, we have constructed splitting methods of the form (24), where ϕ τ is a 2nd-order symmetric approximation of the exact flow. This class of schemes has shown to be particularly efficient for long time integrations of N-body planetary systems in Poincaré Heliocentric coordinates when the leapfrog approximation (29) is considered.
Numerical simulations show that the efficiency of the new integrators presented here is essentially similar in both Jacobi and Heliocentric coordinates. We believe this result is worth remarking, since canonical Heliocentric variables provide very often a more convenient formulation of the problem. The improvement of the new integrators presented here (in particular, methods of order (8,6,4) and (10,6,4)) with respect to previous schemes is most notably exhibited when they are applied for the numerical integration of the outer planets. When the whole Solar System is considered, although methods of order (8,6,4) and (10,6,4) still provide the best results, it is Mercury with its relatively high eccentricity and fast orbital period which constitutes the main limiting factor in all simulations.
When designing the new methods of generalized order (8,6,4) and (10, 6, 4), we have only considered compositions with the minimum number of stages to solve all the order conditions. It might be the case, as in other contexts, that introducing more stages with additional free parameters could lead to more efficient schemes. We intend to explore this possibility and eventually collect the new methods obtained, both in the ABA and BAB classes, in our website (www.gicas.uji.es/software.html).
It is not difficult to check that the sequence of coefficients corresponds precisely to the 5-stage symmetric ABA method of generalized order (10,2) with positive coefficients considered in [12,14,15]. Thus, it is our starting point in the search of an efficient (10,6,4) method.
• We then choose an arbitrary orthogonal matrix M ∈ R 3×10 , and for a randomly chosen complex number γ ∈ C, consider the following oneparameter family of systems of polynomial equations: For a generic M and γ, there exists a unique continuous curve x = ρ(t) ∈ C 10 of solutions of this family of polynomial systems such that ρ(0) = x 0 .
Here t ∈ [0, 1) denotes the continuation parameter. If x = ρ(1) ∈ R 10 , then x is a real solution of the original system (31). We have applied a numerical continuation algorithm to compute such a solution for several values of γ ∈ C, and found two real solutions. The solution x with smaller norm x gives the method ABAH1064 of generalized order (10,6,4) displayed in Table 4. Notice the small absolute value of both a 3 and a 4 in the resulting scheme (recall that this solution has been obtained starting with x 0 such that a 0 3 = a 0 4 = 0).
We have followed a similar procedure to compute solutions starting with x 0 such that a 0 i = a 0 j = 0 for indices (i, j) = (3,4). Such procedure leads to nonreal solutions x for some of the choices of (i, j), and for other choices gives real solutions with larger norm x and larger error terms than method ABAH1064.