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.

💡 Research Summary

The paper addresses the long‑standing challenge of constructing highly accurate symplectic integrators for near‑integrable Hamiltonian systems, with a focus on planetary N‑body problems where very small numerical errors must be maintained over millions of years. The authors begin by formulating the problem in the standard split form H = A + εB, where A represents the exactly integrable Keplerian part and εB the small perturbative interactions among planets. While the classic Wisdom–Holman method and its later refinements (e.g., Laskar–Robutel schemes) are widely used, they are limited to relatively low order and become inefficient when the user demands sub‑machine‑precision accuracy.

To overcome this limitation, the authors develop a systematic theory of order conditions for splitting methods. Using the Baker–Campbell–Hausdorff expansion, they express the overall propagator of a general s‑stage splitting scheme as a formal series in the step size τ. By collecting terms of equal τ‑power, they obtain a set of polynomial equations that the coefficients a_i (weights for the A‑flow) and b_i (weights for the B‑flow) must satisfy in order to achieve a prescribed order k. Crucially, they prove that these equations are both necessary and sufficient, and they identify an independent subset that eliminates redundancies present in earlier works. This “generalized order condition” framework provides a clear map of the degrees of freedom that remain after fixing the desired order.

The paper then distinguishes two coordinate representations commonly employed in planetary dynamics: Jacobi coordinates and Poincaré heliocentric coordinates. In Jacobi coordinates the perturbation B retains a simple pairwise form, so the general order conditions alone guarantee the required accuracy. In contrast, Poincaré coordinates introduce additional nonlinear terms in B, which generate extra Poisson‑bracket contributions in the BCH expansion. Consequently, achieving the same order of accuracy demands an extra polynomial constraint on the splitting coefficients (for example, Σ a_i³ = 0 or Σ a_i b_i² = 0). The authors derive these extra constraints explicitly and show how they modify the feasible coefficient space.

Armed with the full set of algebraic conditions, the authors solve for optimal coefficients. They employ Gröbner‑basis techniques to reduce the symbolic system and then apply numerical optimization (non‑linear least squares combined with evolutionary algorithms) to explore the remaining free parameters. For each coordinate system they construct 4th‑, 6th‑, and 8th‑order schemes, deliberately leaving 3–5 free parameters to be tuned for minimal energy‑error growth and maximal computational efficiency.

Extensive numerical experiments validate the new methods. The authors integrate the full Solar System (Mercury to Neptune) over 10⁷ years using step sizes ranging from 0.5 to 2 days. In Jacobi coordinates, the 6th‑order scheme reduces the global energy error by roughly a factor of three compared with the best existing 4‑step Laskar–Robutel integrator at the same computational cost. In Poincaré coordinates, the newly derived 8th‑order scheme keeps the total energy drift below 10⁻¹⁴ over the same interval, a level comparable to machine precision. Additional tests on synthetic exoplanetary systems confirm that the methods retain their superior accuracy and stability across a variety of mass ratios and orbital configurations.

Performance profiling shows that, for a given target accuracy, the new schemes require 30–45 % less CPU time than the traditional Wisdom–Holman or higher‑order composition methods. This efficiency gain stems from the higher order per stage and the careful balancing of the free parameters to minimize the leading error coefficients.

In summary, the paper delivers a complete, mathematically rigorous framework for deriving high‑order symplectic splitting integrators tailored to specific coordinate choices in celestial mechanics. By providing both the generalized order conditions and the additional constraints required for Poincaré heliocentric coordinates, the authors enable the construction of integrators that are demonstrably more accurate and faster than existing alternatives. The work has immediate implications for long‑term stability studies of the Solar System, high‑precision ephemeris generation, and the dynamical modeling of exoplanetary systems where sub‑millisecond timing accuracy over astronomical timescales is essential.