A Relativistic Motion Integrator: Numerical accuracy and illustration with BepiColombo and Mars-NEXT

for a given metric G µν where X µ = (cT, X, Y, Z) are the coordinates, τ is the proper time and Γ α µν are the Christoffel symbols of the metric considered, derived numerically. As an example, we use this integrator with the planetocentric metric advised by the IAU 2000 resolutions, described in Soffel et al. (2003) and characterised by a scalar potential W and a vector potential W i . Until now, we considered only the central body of mass M , so that we have:

2. where G is Newton’s gravitational constant, c the speed of light, C lm and S lm are related to the central gravity field, R e is the equatorial radius of the central body, while the vector S is its spin moment and R = √ X 2 + Y 2 + Z 2 . It is possible to develop analytically the equations of motion (1.1) at first Post-Newtonian (1PN) order. Doing so, one gets

, where the corrections are composed of different types of forces: a Newtonian correction coming from the harmonics (proportional to C lm or S lm ), a relativistic Schwarzschild acceleration (proportional to 1/c 2 ), a relativistic correction coming from the harmonics (proportional to C lm /c 2 or S lm /c 2 ), a relativistic coupling between harmonics (proportional 1 2 A. Hees and S. Pireaux to C lm C lm /c 2 , C lm S lm /c 2 or S lm S lm /c 2 ) and finally a relativistic Lense-Thirring acceleration (proportional to the spin momentum over c 2 ).

The analytical development described above is used to assess the order of magnitude of each separate effect and to validate the RMI method. Figures 1, 2 and 3 show the separate impact of the relativistic effects in terms of cartesian coordinates (X, Y, Z), radial distance R and L = ω + w, where ω is the argument of pericenter and w is the true anomaly. The orbital parameters of the BepiColombo mission can be found in Balog et al. (2000): a = 3389 km, e = 0.162, i = 90˚. The orbital parameters of the Mars-NEXT mission (see Chicarro et al. (2008)) are: a = 3896 km, e = 0 and i = 75˚. The numerical integration has been performed over 5 orbital periods.

The numerical derivative has to be treated carefully. In the implementation of RMI, we used a fourth-order numerical derivative

For large h, the discretization error is important, while for small h, the roundoff error increases. We use an optimal derivation step computed analytically for a function 1 r (see Figure 4), given by h opt = (45ǫa 6 c 2 /(960GM )) 1/5 (Kincaid and Cheney (2002)) where ǫ is the machine precision and a is the semi-major axis. As can be seen on Figure 4, we derive H µν = G µν -η µν , with η µν the Minkowski metric, instead of G µν since it is more stable numerically. Moreover, we use a Richardson extrapolation in order to increase the precision on the derivative (Richardson (1927)). This extrapolation uses two estimations of the derivative of order 4 with different step size (D h and D h/k with k a real factor) to construct an estimation of order 8 Typically, the value of k is often chosen as 2 or 1.5. We can see on Figure 4 that this procedure can increase the derivative precision. It is furthermore less sensitive to the choice of h opt . With such an implementation, it is possible to show that the relative precision of RMI is of the order of 10 -12 in double precision and 10 -22 in quadruple precision.

We have shown that RMI is useful to compute relativistic orbits for different missions. It includes all the relativistic effects (up to the corresponding order of the metric). It is quite easy to use, since the user only has to change the metric module if he wants to change the metric. Until now, RMI is only a prototype that is more time consuming than the usual 1PN approach. Nevertheless, it is possible to reduce drastically the integration time required by the RMI method via proper coding and using parallelization (the computation of the Christoffel symbols can easily be parallelized).