Symplectic propagators for the Kepler problem with time-dependent mass

New numerical integrators specifically designed for solving the two-body gravitational problem with a time-varying mass are presented. They can be seen as a generalization of commutator-free quasi-Magnus exponential integrators and are based on the compositions of symplectic flows. As a consequence, in their implementation they use the mapping that solves the autonomous problem with averaged masses at intermediate stages. Methods up to order eight are constructed and shown to be more efficient than other symplectic schemes on numerical examples.


Introduction
The problem of determining the motion of two bodies under their mutual gravitational attraction (the so-called Kepler problem) is one of the most studied dynamical systems in classical mechanics. It not only constitutes a paradigmatic example of an integrable system (Arnold 1989), but in addition it is also used as one of the first test bench for illustrating the main features of whatever numerical integration methods for differential equations are proposed, ranging from classical embedded Runge-Kutta schemes with adaptive step size, to multistep methods (Hairer et al. 1993). Moreover, given all its geometric properties, it has received great attention for checking the recent class of structure-preserving numerical methods, such as energy-preserving and symplectic integrators (Hairer et al. 2006;Sanz-Serna and Calvo 1994).
Even in the long-time integration of planetary systems involving several bodies, the Kepler problem plays a fundamental role. This is so because, using an appropriate set of coordinates (e.g., heliocentric or Jacobi coordinates Laskar 1991), the equations of motion can be shown to derive from a Hamiltonian function of the form where H K corresponds to the Keplerian motion of each planet, H I is a (small) perturbation given by the interaction between the planets, and ε denotes a small parameter. An early reference in this respect is Wisdom and Holman (1991), where a symmetric second-order symplectic scheme was successfully used to study the chaotic behavior of the solar system. Since then, many other highly efficient symplectic integrators have been designed taking into account this near-integrable structure of the problem (see, e.g., Blanes et al. 2000;Farrés et al. 2013;Laskar and Robutel 2001;McLachlan 1995, and references therein). All these schemes require, at each intermediate stage, solving one or several times a transcendental equation to determine with great accuracy the position in phase space of each body subjected to the Hamiltonian H K . This is done in practice by numerical iteration as follows. If we write the Hamiltonian function describing the Kepler problem as where μ = G M, in which G is the gravitational constant and M is the reduced mass, q, p ∈ R 3 , and r = q = q T q; the map advancing the solution in time from t 0 to t can be expressed as (q(t), p(t)) ≡ (q 0 , p 0 ; t, μ) = ( f q 0 + g p 0 , f p q 0 + g p p 0 ) in terms of functions f , f p , g, g p that are determined through (see, e.g., Blanes and Casas 2016;Danby 1988) where x is evaluated by numerical iteration. As a matter of fact, this can be done very efficiently with only 2-5 iterations, depending on the time step and the values of q 0 , p 0 to reach round-off accuracy.
There are other astronomical problems that can be also modeled as (5), but, as they frequently involve some loss of mass, the corresponding Hamiltonian system depends explicitly on time. Examples are the evolution of planetary systems with time-dependent stellar loss of mass (Adams et al. 2013), evolution of exoplanets around binary star systems (Rahoma 2016) with stellar mass loss (Veras et al. 2013), and the two-body problem with varying mass (El-Saftawy and El-Salam 2017; Rahoma et al. 2009), among others (Li 2013). The Hamilton systems to be solved are still near integrable, but non-autonomous where now H K corresponds to the Keplerian motion of each planet (with time-dependent mass). If one takes t as a new coordinate, say q t = t, and its associated momentum, p t , the equations obtained from the autonomous Hamiltonian have the same solution for q, p as the original non-autonomous problem. In those situations, one has to solve in an accurate and efficient way the dynamics of the Hamiltonian H K (t, q, p) or, equivalently, each of the independent Hamiltonians, i.e., where now μ(t) is a time-dependent function, and the perturbed part ε H I (t, q, p) has to be solved with the time frozen since this part does not depend on p t . Then, splitting symplectic integrators for near-integrable problems can be used. The natural question arising here is whether the procedure (3, 4) can be conveniently adapted to deal with this problem. In this work, we show that this is actually the case and present several numerical algorithms with different orders of accuracy that use the map in (3) with properly averaged values of the time-dependent mass μ(t) that are more efficient than other numerical integrators in shortand long-time integrations. The new algorithms can be seen as a specially tuned class of schemes called commutator-free quasi-Magnus (CFQM) exponential integrators, originally intended for linear non-autonomous problems (Alvermann and Fehske 2011;Blanes et al. 2017Blanes et al. , 2018. They are based on compositions of the symplectic flows corresponding to certain linear combinations of the time-dependent Hamiltonian (7) evaluated at appropriately chosen times. For illustration, the map given by the composition where 3 , provides a fourth-order approximation to the exact solution of (7) at time t 1 = t 0 + h, in the sense that (q 1 , p 1 ) = (q(t 1 ), p(t 1 )) + O(h 5 ).
In the following, we show how scheme (8) is obtained and also construct similar approximations of higher orders. Then, we illustrate the advantages of the new procedures in comparison with other numerical algorithms on some examples.

Hamiltonian vector fields and Poisson brackets
The derivation process of the new schemes can be carried out in a suitable way by introducing the formalism of Lie derivatives and Lie transformations (Abraham and Marsden 1978;Arnold 1989). Thus, starting from the equations of motion for the autonomous case (2) and the corresponding Hamiltonian vector X H = (∇ p H , −∇ q H ) T , one introduces the Lie derivative L X H , whose action on a differentiable function G(q, p) is Here, x = (q, p) T , J is the basic canonical symplectic matrix and {H , G} denotes the Poisson bracket of the two scalar functions H (q, p) and G(q, p), so that, in terms of the Poisson bracket, Eq. (9) can be written simply aṡ In case G is a vector function, L X H in (10) acts on any of its components. If ϕ t denotes the flow corresponding to (9), for each infinitely differentiable map G, the Taylor series of G(ϕ t (x 0 )) at t = t 0 is given by where exp(t L X H ) is the so-called Lie transformation. If we introduce the operator t acting on differentiable functions as G(ϕ t (x)) = t [G](x), then we can write so that the solution of (9) is obtained by replacing G(x) in (12) by the identity map Id(x) = x: Lie transformations obey the following important property for any pair of arbitrary vector fields. Given the flows ϕ [1] t 1 and ϕ [2] t 2 corresponding to the differential equationsẋ . Notice where the indices 1 and 2 appear depending on whether one is dealing with maps or with exponentials of operators. This relation generalizes by induction to any number of flows (Hairer et al. 2006).
In the non-autonomous problem, one can still formally write the operator t associated with the exact flow as the Lie transformation corresponding to an unknown function (t, t 0 ), which can be approximated for sufficiently small time intervals [t 0 , t] by truncating the corresponding Magnus expansion (Blanes 2018;Oteo and Ros 1991). This is done by the infinite series which involves multiple integrals of nested Poisson brackets: Recursive procedures to generate the terms in the Magnus series are found, e.g., in Blanes et al. (2009).

Derivation of the new methods
The methods we propose for approximating the operator t up to order r in the general non-autonomous case (7) have the form Alternatively, in terms of maps, the approximation reads where ϕ [ j] h denotes the flow corresponding to the Lie transformation exp(B j (h)). Eventually, in the setting of the non-autonomous Kepler problem, the particularly simple algebraic structure of the problem will allow us to include additional Poisson brackets in the linear combinations B j to reduce the number of flows. In this way, the computational cost is reduced and the overall efficiency is improved. Note that, by construction, these schemes preserve the symplectic character of the exact flow.
The requirement that scheme (15) provides an approximation to the exact solution up to a given order r implies that the coefficients a jk , c k have to satisfy a certain number of conditions. These can be derived by reproducing the exact solution provided by the Magnus expansion (13) with (15) up to order r .
For simplicity in the presentation, we take t 0 = 0 and the Gauss-Legendre quadrature rule, c i , i = 1, . . . , s, of order 2s, although the schemes can be easily adapted to any other quadrature rule. Then, Theorem 2.1 in Blanes (2018) applied to a non-autonomous Hamiltonian system (11) leads to the following Theorem 1 If c i , i = 1, . . . , s, are the s nodes of the Gauss-Legendre quadrature rule of order 2s, L i (t) denotes the Lagrange polynomial and x(h) is the exact solution of (11) at t = t 0 + h, then the solution of the differential equationẏ with Since H (2s) (t, y, h) is a polynomial in t of degree s − 1, we can compute analytically all the terms in the Magnus expansion corresponding to the initial value problem (18) and get an approximation to the exact solution x(h) of (11) up to a given order r by computing the first terms of the series (14). Specifically, for the non-autonomous Kepler problem one has In terms of the map , the scheme is nothing but the well-known second-order midpoint rule We next apply the same strategy to construct higher-order approximations.
Order four. The quadrature rule with s = 2 has nodes and weights respectively, and the interpolating polynomial reads . Now the Lie transformation associated with [4] (and therefore the exact flow) can be correctly reproduced up to order O(h 4 ) by Here and in the sequel, for simplicity, we denote byα i the Lie derivative associated with the vector field corresponding to α i , i.e.,α i = L X α i . Finally, expressing α 1 , α 2 in terms of the evaluations H 1 , H 2 results in or in a more compact way The map corresponding to this operator is precisely the scheme given in (8).
Order six. In this case, one has to use the nodes and weights of the sixth-order Gauss-Legendre quadrature, Then, or, equivalently, Notice that, since only α 1 depends on momenta (through T ), {α 2 , α 3 } = 0, and thus, the number of conditions required to approximate exp(L X [6] ) is six instead of seven. Moreover, a simple calculation shows that As already illustrated by the previous fourth-order approximation, a linear combination of α i , i = 1, 2, . . ., gives rise to a Hamiltonian function which corresponds to a scaled autonomous Kepler problem with a modified (but constant) massμ, i.e., with x 1 = 0, whose exact flow can be determined with the map given in (3). On the other hand, and according to (23), the flow of can also be trivially obtained, since it only depends on coordinates. In consequence, we propose the following composition to approximate exp(L X [6] ) up to order six. Here, α 2 ,α 1 ,α 2 ≡ α 2 , [α 1 ,α 2 ] is the Lie bracket corresponding to the Hamiltonian vector field of the function {α 1 , α 2 , α 1 }. It turns out that there is only one solution for the coefficients x i , namely In consequence, the scheme can be expressed in a compact way as and The step (q n , p n ) → (q n+1 , p n+1 ) with this scheme can be obtained with Algorithm 1. It requires basically the same computational effort as the fourth-order method (8) and (22), since the evaluation of r in the last map, r = Q 2 , can be reused in the first map in the next step as well as in the next map . Accordingly, the computational cost of evaluating the first and fourth flows in (25) can be neglected.
We have analyzed several time-symmetric compositions with enough parameters to satisfy the 16 conditions required to achieve order eight with a reduced number of exponentials associated with the Kepler map, since this is the most expensive part. It turns out that at least five maps for the scaled autonomous Kepler problem are then necessary.
The following composition is first considered with  (29) Algorithm 2 shows how to advance one step using this scheme.
Algorithm 2: Eighth-order method [8] 5a for one step t n → t n+1 = t n + h . . , 7, j = 1, 2, 3, 4; 3 a i = 4 j=1 a i j , i = 2, 3, 4, 5, 6; 4 M i = 4 j=1 a i j μ j , i = 1, 2, . . . , 7; 5 μ i = 4 j=1 D i j μ j , i = 1, 2, 3; 6 Q 0 = q n , P 0 = p n ; 7 r = Q 0 ; We have analyzed other compositions with additional exponentials and free parameters that do not increase the overall computational cost. Thus, in particular, the following scheme only contains five exponentials involvingα 1 (corresponding to five Kepler maps): with 2, . . . , 9, j = 1, 2, 3, 4, so that x 52 = x 54 = 0. This time-symmetric composition has 19 variables, and thus, one can choose three of them as free parameters (in particular, x 42 , x 43 , and x 44 ). We have observed in practice that the performance of the resulting methods obtained according to different criteria depends on the particular problem to be solved. For example, if we take x 42 = x 43 = x 44 = 0, there are four real solutions, and the best results correspond to x 11 = 0.67021911442375565293, Alternatively, if in the composition (28) we take x 31 = x 51 = 0, the cost of the scheme is reduced from 5 to 3 maps, but we end up with only 15 parameters to solve the 16 order conditions. The composition is given by With this, obviously, not all eighth-order conditions can be satisfied. We have considered two possibilities: Either we discard the equation corresponding to [111112] or This scheme can be considered as an optimized sixth-order method that uses an eighth-order quadrature rule and satisfies many order conditions at order eight.

Other symplectic integrators to compare
It is worth comparing the mentioned time-symmetric symplectic integration algorithms specifically designed for the non-autonomous Kepler problem with other well established symplectic schemes such as splitting and composition methods. It is not our purpose here to provide a full treatment of this class of methods, but just to indicate how they can be applied to the particular problem we are considering here. For further details, we refer, e.g., to Blanes (2018), Blanes and Casas (2016), Hairer et al. (2006), andMcLachlan andQuispel (2002). As indicated in "Introduction," if one takes t as a new coordinate, say q t = t, and its associated momentum, p t , the equations obtained from the autonomous Hamiltonian have the same solution for q, p as the original non-autonomous problem. Since the Hamiltonian K 1 = 1 2 p T p − μ(q t ) 1 r can be considered as one describing an autonomous Kepler problem (q t does not change), and the dynamics of K 2 = p t is trivial (it only advances q t , i.e., time), we can use splitting methods defined by a set of coefficients for the numerical integration of (34). This results in the following algorithm for advancing from (q n , p n ) to (q n+1 , p n+1 ): where μ i = μ(t n + d i h) and d i = i j=1 a j . On the other hand, we can use the second-order midpoint rule (21) as a basic integrator of m-stage symmetric composition methods of order r (Blanes and Casas 2016;Hairer et al. 2006;McLachlan and Quispel 2002), with α m+1−i = α i . It turns out that the same algorithm (35) can be applied to implement scheme (36) with Methods of orders r = 4, 6, 8, and 10 involving up to m = 35 basic integrators ϕ [2] α i h have been obtained in the literature (see Hairer et al. 2006 and references therein).
Splitting and composition methods of order higher than two necessarily involve some negative coefficients and therefore a backward integration in time at some inner stages. If, in the specific problem, we are considering μ(t) is a decreasing function, this represents a non-physical effect of an increment of the mass.
Commutator-free exponential integrators, originally proposed for the numerical integration of linear problems (Alvermann and Fehske 2011;Blanes et al. 2017Blanes et al. , 2018, can also be adapted to the nonlinear setting, eventually resulting in schemes of the form (15). Timesymmetric methods of order four with J = 2, 3, of order six with J = 5, 6, and of order eight with J = 11 can be found in Alvermann and Fehske (2011). We should notice that, whereas scheme (22) belongs to this class, this is not the case of the new integrators of orders six and eight, which have been specifically designed taking into account the special structure of the non-autonomous Kepler problem.

A pair of numerical examples
Next we compare the performance of the methods proposed in this work with some of the schemes mentioned before. Specifically, we consider m-stage symmetric compositions of order r of the form with the midpoint rule (21) as basic integrator ϕ [2] h and m = 2k + 1. In particular, we take -SS [4] 5 : The 5-stage fourth-order Suzuki composition, -SS [6] 9 : The 9-stage sixth-order composition, -SS [8] 17 : The 17-stage eighth-order composition, whose coefficients can be found, e.g., in Hairer et al. (2006). Within the class of commutator-free (quasi-)Magnus exponential integrators, we select, in addition to CF [4] 2 given by (8), the following schemes presented in Alvermann and Fehske (2011): : The 3-stage optimized fourth-order composition, -CF [6] 5 : The 5-stage optimized sixth-order composition, -CF [8] 11 : The 11-stage eighth-order composition.
In all cases, the computational cost is counted as the number of evaluations of the map . Finally, for the sake of completeness, we also compare with a pair of ODE solvers provided by MATLAB, namely -RK [4(5)] : Ode45 Runge-Kutta solver, based on the embedded Dormand-Prince 4(5) pair.
-AB [13] : Variable-step and variable-order ode113 Adams-Bashforth-Moulton solver, designed to be more efficient than ode45 at problems with stringent error tolerances.
In this last case, it is less obvious how to estimate the computational cost in comparison with the previous algorithms. We have counted the number of times the solver calls to the vector field. However, to take into account the relative cost with respect to the cost of the map , one should know the computational cost for estimating the local and global errors and to change either the time step or the order of the method as well as the cost of evaluating the mass μ(t). (For some problems, this could be the most costly part, depending on the model used.) For these reasons, we have decided to count two evaluations of the vector field (and the cost to change the order and/or time step) as the cost to evaluate one map . Obviously, different choices for the relative cost would result in a (small) shift of the corresponding curves, although the overall conclusions remain valid.

Example 1
We illustrate the performance of the new schemes on the classical Eddington-Jeans law for the secular evolution of mass in binary systems (Hadjidemetriou 1967), or equivalently Note that, for non-integer values of δ, the computational cost of evaluating μ(t) cannot be neglected. We consider the 2-dimensional case with δ = 1.4, μ 0 = 1, γ = 10 −2 , initial conditions with e = 0.2 and e = 0.8 and final time t f = 20. We compute the solution at the final time, (q(t f ), p(t f )), numerically to high precision and plot the two-norm error of (q, p) at the final time versus the computational cost (measured as the number of times the Kepler map φ is called) for different choices of the time step. Figure 1 shows the results obtained for methods of each order: fourth-order in the top, sixth-order in the second, and eighth-order in the third row. The graphs for the smaller value of the eccentricity are given in the left column. The best method of each order, along with the classic reference methods, ode45 and ode113 of MATLAB, are summarized in the last row.
From the graphs we conclude that the symplectic methods have better general performance than the classical adaptive methods. Moreover, the adaptation of methods by the inclusion of cheap Poisson brackets leads to a notable performance gain compared to the general-purpose CF methods.
Among the new schemes and for this type of problems, [6] 2 is the method of choice for low and medium accuracies, while [8] 5b should be chosen for high accuracies.
Example 2 As a second illustration, we consider the two-dimensional case with a decreasing function given by μ(t) = 1 + exp − 1 5 (t + 1 4 sin 2 (4t)) , and the same initial conditions as in the previous example for the time interval, t ∈ [0, 20]. Figure 2 supports the results of Example 1. In this case, with a faster decaying mass, the symplectic methods perform better than the classical methods. The main change we see is that in this example the difference between the sixth-order methods, [6] 2 and , is more notable, depending on the value of e. In addition, the second eighth-order methods are marginally better in this problem.

Conclusions
In this work, we have considered the numerical integration of the two-body gravitational problem with a time-varying mass. The exact flow corresponds to a symplectic transformation, and different symplectic integrators from the literature can be adapted to solve the non-autonomous systems. However, none of these symplectic methods were designed to  solve Hamiltonian systems with this particular structure. This is a relevant problem, and new specifically designed symplectic integrators have been built. These new schemes can be seen as a generalization of the commutator-free quasi-Magnus exponential integrators and are based on compositions of symplectic flows corresponding to linear combinations of the Hamiltonian function and certain Poisson brackets. The implementation makes use of the mapping that solves the autonomous problem for averaged masses at each intermediate stage.
In the autonomous case, the schemes provide the exact solution, so they also show a high performance in the adiabatic limit.
We have built time-symmetric methods of order four, six and eight that can be used with any quadrature rule of the order of the method or higher. Some of the proposed methods are optimized using a quadrature rule of higher order than the order of the method as well as by adding free parameters into the scheme in order to satisfy certain order conditions at higher orders.
Since the proposed methods provide the exact solution in the limit when the mass is constant, it is clear that the performance of the algorithm will improve if a time step control is used. For instance, one could take into account the time derivative of the mass in applications that involve very slowly varying mass.
The new methods have shown to be more efficient than other symplectic schemes to all desired accuracies on several numerical examples.

Compliances with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.