New general relativistic contributions to Mercury's orbital elements and their measurability

We numerically and analytically work out the first-order post-Newtonian (1pN) orbital effects induced on the semimajor axis a, the eccentricity e, the inclination I, the longitude of the ascending node Ω, the longitude of perihelion ̟, and the mean longitude at epoch ǫ of a test particle orbiting its primary, assumed static and spherically symmetric, by a distant massive third body X. For Mercury, the rates of change of the linear trends found are İ 1pN = −4.3 microarcseconds per century ( μas cty−1 ) , Ω̇ X 1pN = 18.2 μas cty−1,  ̟̇ X 1pN = 30.4 μas cty−1, ǫ̇ 1pN = 271.4 μas cty−1, respectively. Such values, which are due to the added actions of the other planets from Venus to Saturn, are essentially at the same level of, or larger by one order of magnitude than, the latest formal errors in the Hermean orbital precessions calculated with the EPM2017 ephemerides. The perihelion precession  ̟̇ X1pN turns out to be smaller than some values recently appeared in the literature in view of a possible measurement with the ongoing BepiColombo mission. Linear combinations of the supplementary advances of the Keplerian orbital elements for several planets, if determined experimentally by the astronomers, could be set up in order to disentangle the 1pN N-body effects of interest from the competing larger precessions like those due to the Sun’s quadrupole moment J2 and angular momentum S. keywords gravitation − celestial mechanics − ephemerides − methods: miscellaneous


Introduction
In its weak-field and slow-motion approximation, general relativity 1 predicts that, in addition to the time-honored first-order post-Newtonian (1pN) gravitoelectric and gravitomagnetic precessions induced by the mass monopole M (Schwarzschild) and the spin dipole S (Lense-Thirring) moments of the central body acting as source of the gravitational field, further 1pN orbital effects due to the presence of other interacting masses arise as well . Let us consider a nonrotating primary of mass M, assumed as origin of a locally inertial coordinate system, orbited by a test particle located at r and moving with velocity v. If a distant, pointlike body X of mass M X is present at r X and moves with velocity v X with respect to M, the test particle experiences certain 1pN accelerations which, from Eq. (4) of , are 2G 2 MM X c 2 r 3 X r − 6 (r ·r X )r X + 3 (r ·r X ) 2r , -3 - In Eqs.
(1) to (3), which are a particular case of the full 1pN equations of motion for a system of N pontlike, massive bodies mutually interacting through gravitation (Poisson & Will 2014, Eq. (9.127)), G is the Newton's gravitational constant, and c is the speed of light in vacuum.  looked at the longitude of perihelion ̟ of Mercury finding an additional contribution to its 1pN secular precession of abouṫ ̟ X 1pN = 0.22 milliarcseconds per century mas cty −1 = = 220 microarcseconds per century µas cty −1 .
Eq. (4) was obtained by making some simplifying assumptions about the orbital geometries of both the perturbed and the perturbing bodies, and includes the combined actions of Venus, Earth, Mars, Jupiter and Saturn. It should be a direct effect of the accelerations of Eqs.
(1) to (3), and an indirect consequence of the interplay between the usual Newtonian N−body pull by the other planets and the Sun-only 1pN gravitoelectric acceleration. Eqs.
(1) to (3) and all the standard Newtonian and 1pN N-body dynamics is routinely modeled in the data reduction softwares of the teams of astronomers producing the planetary ephemerides like the Development Ephemeris (DE) by the NASA Jet Propulsion Laboratory (JPL) in Pasadena (Folkner et al. 2014  , it would be so because the expected ≃ 10 −6 accuracy with which the parameterized Post-Newtonian (PPN) parameters β, γ should be measured by such a spacecraft would correspond to an uncertainty in the main contribution to the Mercury's 1pN perihelion precession̟ 1pN = 42.98 arcseconds per century ′′ cty −1 as little as δ̟ 1pN ≃ 0.03 mas cty −1 = 30 µas cty −1 .
Iorio (2018), after having pointed out that the indirect, mixed 2 effects should likely be not measurable in practical planetary data reductions, analytically worked out the direct perihelion precessions due to Eqs. (1) to (3) for arbitrary orbital configurations of both the test particle and the perturbing body X. The total 1pN rate of change induced on the perihelion of Mercury by all the other planets of the solar system from Venus to Saturn would amount to (Iorio 2018, Table 2) ̟ X 1pN = 0.15 mas cty −1 = 150 µas cty −1 .
Iorio (2018) showed also that Equation (6) would likely be overwhelmed by the larger systematic errors due to the mismodeling in the competing secular precessions due to the Sun's oblateness J 2 and angular momentum S (1pN Lense-Thirring effect).
In this paper, we will show that the value reported in Equation (6) is, in fact, wrong because of an error by Iorio (2018) in the calculation of the precession due to Equation (2). The correct size of the overall 1pN N−body perihelion precession of Mercury will turn out to be even smaller than Equation (6), thus enforcing the pessimistic conclusions of Iorio (2018) about its possible measurability. As such, we will further explore the consequences of Eqs. (1) to (3) by numerically working out the secular shifts induced by them on all the other orbital elements, i.e. the semimajor axis a, the eccentricity e, the inclination I, the longitude of the ascending node Ω, and the mean longitude at epoch ǫ, and will compare them with the uncertainties in the planetary orbital motions inferred by Iorio (2019) from the most recent version of the EPM ephemerides (Pitjeva & Pitjev 2018). Indeed, if and when the astronomers will observationally produce the supplementary rates of change ∆ȧ obs , ∆ė obs , ∆İ obs , ∆Ω obs , ∆̟ obs , and ∆ǫ obs of as many planets as possible, it will be possible to generalize the approach proposed by 3 Shapiro (1990) by suitably combining them in order to disentangle the effects of Eqs. (1) to (3) in from the other competing precessions due to, e.g., the Sun's J 2 and S.
2. The 1pN N−body secular changes of the orbital elements

Numerical integration of the equations of motion
We simultaneously integrate the equations of motion of Mercury in Cartesian rectangular coordinates and the Gauss equations for each orbital element with and without the fifteen terms of the sum of Eqs. (1) to (3) calculated for Venus, Earth, Mars, Jupiter and Saturn over a time span as long as 100 cty in order to clearly single out the sought features of motion: both runs share the same initial conditions retrieved on the Internet from the WEB interface HORIZONS maintained by the JPL. For consistency reasons with the planetary data reductions available in the literature, we use the equatorial coordinates of the International Celestial Reference System (ICRS). Then, for each orbital element, we plot in Fig. 1 the time series (blue curve) resulting from the difference between the runs with and without the 1pN N−body accelerations. Finally, we fit a linear model (yellow line) to its numerically produced signal, and estimate its slope: the outcome is collected in the caption of Fig. 1. From Fig. 1, the secular trends of I, Ω, ̟, ǫ are apparent, while a and e seem to experience long-term harmonic variations. The size of the slopes of the precessions of the angular rates of change vary in the range ≃ 1 − 100 µas cty −1 = 0.001 − 0.1 mas cty −1 . In particular, it turns out that the secular precession of the perihelion is about five times smaller than Equation (6) (Iorio 2018, Table 2), being as little aṡ ̟ X 1pN = 30 µas cty −1 = 0.03 mas cty −1 .
Numerical tests conducted by switching off from time to time each of Eqs.
(1) to (3) for every single perturbing planet X showed that the issue resides in the analytical calculation of Eq. (B5) in Iorio (2018) and in the consequent numerical results of the third column from the left of Table 2 in Iorio (2018).

Analytical calculation
It is also possible to analytically work out the long-term rates of change of the Keplerian orbital elements of the test particle with the Gauss perturbative equations applied to Eqs. (1) to (3) by doubly averaging their right-hand-sides over the orbital periods P b and P X of the perturbed body and the perturber X, respectively. The resulting expressions, especially those due to Eqs. (1) to (2), are very cumbersome. Thus, we display just approximate formulas for them to their leading order in e. The shifts due to Equation (3), which are relatively less involved, are displayed in full. In the next Sects., we use the shorthand ∆Ω Ω − Ω X .
It turns out that there is an excellent agreement among the numerical results of Sect. 2.1 and the analytical results shown below.
2.2.1. The doubly averaged rates of change of the orbital elements due to A G 2 Here, we analytically calculate the doubly averaged rates of change of the Keplerian orbital elements of the test particle, to their leading order in e, due to Equation (1). No further approximations in the orbital configurations of both the perturbed body and X are made. They are as follows.
The semimajor axis a stays constant sincė The rate of change of the eccentricity e turns out to bė with As far as the rate of change of the inclination I is concerned, we havė with I A G 2 sin I X (cos I cos I X + sin I sin I X cos ∆Ω) sin ∆Ω.
The precession of ̟ due to Equation (1) was correctly worked out, to the zero order in e, in Eq. (B2) of Iorio (2018); thus, we do not display it here.
2.2.2. The doubly averaged rates of change of the orbital elements due to A G Here, we analytically work out the doubly averaged rates of change of the Keplerian orbital elements of the test particle, to their leading order in e, induced by Equation (2). No further approximations in the orbital configurations of both the perturbed body and X are made. We list them below.
For the semimajor axis a, we havė with A A G = sin I X (− sin I cos I X + cos I sin I X cos ∆Ω) sin ∆Ω.
The rate of change of the eccentricity e iṡ with E A G = sin I X (− sin I cos I X + cos I sin I X cos ∆Ω) sin ∆Ω.
The rate of change of the inclination I turns out to bė with I A G sin I X (cos I cos I X + sin I sin I X cos ∆Ω) sin ∆Ω.
Eq. (25)-eq. (26), which correct Eq. (B5) of Iorio (2018), allow to calculate the same values for Mercury which are obtained with our numerical integrations of Sect. 2.1, limited to Equation (2) only, for each of the perturbing planets at a time.
2.2.3. The doubly averaged rates of change of the orbital elements due to A v X Here, we analytically calculate the doubly averaged rates of change of the Keplerian orbital elements of the test particle caused by Equation (3). No approximations in the orbital configurations of both the perturbed body and X are made; the following expressions are exact.
The semimajor axis a and the eccentricity e are constant sincė The rate of change of the inclination I iṡ For the precession of the node Ω we havė The precession of ̟ due to Equation (3) was correctly calculated in Eq. (B8) of Iorio (2018); as such, it is not shown here.
The rate of change of the mean longitude at epoch ǫ does depend on e. It turns out to bė where

Confrontation with the observations
Iorio (2019) attempted to calculate the formal uncertainties in the secular rates of change of a, e, I, Ω, and ̟ of the planets of the solar system from the recently released formal errors in a and the nonsingular orbital elements e sin ̟, e cos ̟, sin I sin Ω, and sin I cos Ω estimated for the same bodies with the EPM2017 ephemerides by Pitjeva & Pitjev (2018). Since, among other things, the 1pN N-body equations of motion are routinely included in the EPM software dynamics, such errors should be overall regarded as representative of the current level of modeling the solar system dynamics along with measurement errors. As such, they may be viewed as the uncertainties that would affect a putative measurement of the effects worked out in Sect. 2 if they were explicitly measured in some dedicated data analysis. From the column dedicated to Mercury in Table 1 of Iorio (2019), it can be noted that the 1 − σ error inȧ amounts to δȧ obs = 0.003 m cty −1 , while for the other Keplerian orbital elements we have δė obs = 0.6 µas cty −1 , δİ obs = 3 µas cty −1 , δΩ obs = 24 µas cty −1 , and δ̟ obs = 8 µas cty −1 . From a comparison with the expected 1pN rates of change of Fig. 1, it turns out that, with the possible exception of the perihelion, they are about of the same order of magnitude of the aforementioned uncertainties. Moreover, as discussed in Pitjeva & Pitjev (2018) and Iorio (2019), the latter ones may be optimistic. Thus, it is difficult to deem the predicted 1pN N-body precessioṅ ̟ X 1pN = 30 µas cty −1 as realistically measurable compared to a merely formal uncertainty δ̟ obs = 8 µas cty −1 . It is worth noticing that such a tiny error would correspond to current bounds in the PPN parameters β, γ as little as ≃ 10 −7 , which are better than the expected accuracy from the ongoing BepiColombo mission quoted by ; see the discussion in Iorio (2019) about the reliability of such an evaluation. The mean longitude at epoch ǫ seem, at first sight, more interesting since its 1pN N-body rate is as large asǫ X 1pN = 270 µas cty −1 = 0.27 mas cty −1 . Iorio (2019) did not calculate the uncertainty inǫ. In their Table 3, Pitjeva & Pitjev (2018) released the formal uncertainty in the planetary mean longitudes, dubbed there as λ; for Mercury, it is as little as δλ obs = 3.3 µas. This implies that, in order to retrieve the uncertainty inǫ, the errors in the mean motion n b due to the mismodeling of the Sun's gravitational parameter µ and of the planet's semimajor axis are required as well. Since δµ obs = 1 × 10 10 m 3 s −2 (Pitjeva 2015a), the resulting error in the Hermean mean motion is as large as δn obs b = 20 mas cty −1 . It vanishes the possibility of measuring the 1pN N-body effect on ǫ. As such, only a dramatic improvement in the determination of the Hermean orbit, which might be obtained when all the data from BepiColombo will be collected and processed, may bring the 1pN N-body precessions due to the direct effect of Eqs. (1) to (3) in the measurability domain.
On the other hand, even should this finally be the case, the concerns raised by Iorio (2018) about the systematic errors caused by the competing Sun's quadrupole and Lense-Thirring rates of change are even reinforced by the present analysis since the actual size of the 1pN N-body perihelion precession of Mercury turned out to be smaller than the incorrect value of Equation (6). Thus, it is hopeful that the astronomers will finally provide the community with the supplementary advances of all the other Keplerian orbital elements in addition to the perihelion. Indeed, if and when it will happen, it would, then, be possible to set up linear combinations of them suitably designed to cancel out, by construction, the other unwanted precessions. An analogous approach, originally limited just to the perihelia of other planets and asteroids in order to separate the disturbing Sun's J 2 action from the Schwarschild-type rates of changes was proposed by Shapiro (1990). It is also widely used in ongoing relativistic tests with geodetic satellites in the Earth's field; see, e.g., Renzetti (2013), and references therein for an overview.

Summary and conclusions
Recently, Will (2018) calculated a new general relativistic contribution to the Mercury's perihelion advance as large as̟ X 1pN = 220 µas cty −1 arising from an approximated form of the 1pN N-body equations of motion restricted to a hierarchic three body system. He claimed that it may be measured in the next future by the ongoing BepiColombo mission to Mercury if it will reach a ≃ 10 −6 accuracy level in constraining the PPN parameters β, γ. Later, the present author first remarked in Iorio (2018) that the indirect precession due to the interplay of the Newtonian N-body and the 1pN Sun's Schwarzschild-like accelerations in the equations of motion is likely undetectable in actual data reduction since it cannot be expressed in terms of a dedicated, solve-for parameter. Then, he calculated analytically the individual contributions to the perihelion advance induced directly by each of the approximated 1pN N-body accelerations put forth by  by finding an overall precession of̟ X 1pN = 150 µas cty −1 . Iorio (2018) discussed also the impact of the systematic aliasing due to the competing perihelion rates induced by the Sun's quadrupole mass moment J 2 and angular momentum A via the Lense-Thirring effect by noting that their mismodeling would likely compromise a clean recovery of the 1pN effect of interest.
Here, the secular rates of change of all the other Keplerian orbital elements a, e, I, Ω, ̟, and ǫ caused by the same approximated 1pN N-body accelerations by  were analytically worked out. A numerical integration of the equations of motion confirmed such findings in the case of Mercury acted upon by the other planets from Venus to Saturn. The resulting rates of change amount toİ X 1pN = −4.3 µas cty −1 ,Ω X 1pN = 18.2 µas cty −1 ,̟ X 1pN = 30.4 µas cty −1 , ǫ X 1pN = 271.4 µas cty −1 . As a result, the Hermean 1pN N-body perihelion precession turned out to be smaller than the previously reported values because of an error explicitly disclosed, at least in the calculation by Iorio (2018). This makes even more difficult than before its possible present and future measurement. A comparison with the merely formal uncertainties in some of the orbital secular rates of Mercury, recently obtained by Iorio (2019) from the EPM2017 ephemerides, showed that the sizes of the predicted 1pN N-body precessions are just at the same level or even below them if, more realistically, they are rescaled by a factor of ≃ 10 − 50 (Iorio 2019). If our future knowledge of the orbit of the closest planet to the Sun will be adequately improved, the systematic bias caused by other competing precessions could be removed by suitably designing linear combinations of the other Keplerian orbital elements of Mercury, provided that the astronomers will determine also their supplementary advances in addition to the perihelion's one. (1) to (3) for X ranging from Venus to Saturn over a time span 100 cty long. The units are m for a and microarcseconds (µas) for all the other orbital elements. They were obtained for each orbital element as differences between two time series calculated by numerically integrating the barycentric equations of motion of all the planets from Mercury to Saturn in Cartesian rectangular coordinates with and without the aforementioned 1pN N-body accelerations. The initial conditions, referred to the Celestial Equator at the reference epoch J2000, were retrieved from the WEB interface HORIZONS by NASA JPL; they were the same for both the integrations. The slopes of the secular trends, in yellow, fitted to the blue time series of ∆I(t), ∆Ω(t), ∆̟(t), and ∆ǫ(t) areİ X 1pN = −4.3 µas cty −1 ,Ω X 1pN = 18.2 µas cty −1 ,̟ X 1pN = 30.4 µas cty −1 ,ǫ X 1pN = 271.4 µas cty −1 , respectively.