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 Ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document}, the longitude of perihelion ϖ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varpi $$\end{document}, and the mean longitude at epoch ϵ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon $$\end{document} 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 I˙1pNX=-4.3microarcsecondspercenturyμascty-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{I}_\mathrm {1pN}^\mathrm {X} = -4.3\,\mathrm {microarcseconds\,per\,century}\,\left( \upmu \mathrm {as\,cty}^{-1}\right) $$\end{document}, Ω˙1pNX=18.2μascty-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\dot{\Omega }}_\mathrm {1pN}^\mathrm {X} = 18.2\,\upmu \mathrm {as\,cty}^{-1}$$\end{document}, ϖ˙1pNX=30.4μascty-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\dot{\varpi }}_\mathrm {1pN}^\mathrm {X} = 30.4\,\upmu \mathrm {as\,cty}^{-1}$$\end{document}, ϵ˙1pNX=271.4μascty-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\dot{\epsilon }}_\mathrm {1pN}^\mathrm {X} = 271.4\,\upmu \mathrm {as\,cty}^{-1}$$\end{document}, 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 ϖ˙1pNX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\dot{\varpi }}_\mathrm {1pN}^\mathrm {X}$$\end{document} 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\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_2$$\end{document} and angular momentum S\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{S}$$\end{document}.


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 grav- 1 See, e.g., Debono & Smoot [2] and references therein for a recent overview on its status and challenges. a e-mail: lorenzo.iorio@libero.it (corresponding author) itomagnetic 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 [13]. 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 Will [13], are 2G 2 M M X c 2 r 3 X r − 6 r ·r X r X + 3 r ·r X 2r , (1) 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 2 ([9], Eq. (9.127)), G is the Newton's gravitational constant, and c is the speed of light in vacuum. Will [13] 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 [3], the Intégrateur Numérique Planétaire de l'Observatoire de Paris (INPOP) by the Institut de Mécanique Céleste et de Calcul des Éphémérides (IMCCE) at the Paris Observatory [12], and the Ephemeris of Planets and the Moon (EPM) by the Institute of Applied Astronomy (IAA) of the Russian Academy of Sciences (RAS) in Saint Petersburg [7]. Will [13] claimed that Eq. (4) would likely be detectable with the ongoing BepiColombo mission to Mercury. According to Will [13], 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 pre-cession˙ 1pN = 42.98 arcseconds per century cty −1 as little as δ˙ 1pN 0.03 mas cty −1 = 30 μas cty −1 .
Iorio [4], after having pointed out that the indirect, mixed 3 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 ([4], Table 2) X 1pN = 0.15 mas cty −1 = 150 μas cty −1 .
Iorio [4] showed also that Eq. (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 Eq. (6) is, in fact, wrong because of an error by Iorio [4] in the calculation of the precession due to Eq. (2). The correct size of the overall 1pN N −body perihelion precession of Mercury will turn out to be even smaller than Eq. (6), thus enforcing the pessimistic conclusions of Iorio [4] 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 [5] from the most recent version of the EPM ephemerides [8]. 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 4 Shapiro [11] 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.

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 1 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 HORI-ZONS 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 (1) to (3) for X ranging from Venus to Saturn over a time span 1 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 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 Eq. (6) ( [4], 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 [4] and in the consequent numerical results of the third column from the left of Table 2 in Iorio [4].

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 Eq. (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.

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 Eq. (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ė As far as the rate of change of the inclination I is concerned, we havė with The precession of the node iṡ with N A G 2 − 2 cos 2I csc I sin 2I X cos + + cos I cos 2I X (3 + cos 2 ) + 2 sin 2 .
The precession of due to Eq. (1) was correctly worked out, to the zero order in e, in Eq. (B2) of Iorio [4]; thus, we do not display it here.
The rate of change of the mean longitude at epoch iṡ where L A G 2 = −1 + 3 cos I − 3 cos 2I X + 9 cos I cos 2I X + 12 sin 2 I 2 sin 2 I X cos 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 Eq. (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 ) 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 .
For the precession of the longitude of perihelion , we havė The rate of change of the inclination I iṡ For the precession of the node we havė The precession of due to Eq. (3) was correctly calculated in Eq. (B8) of Iorio [4]; 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 [5] 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 [8]. Since, among other things, the 1pN Nbody 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 [5], 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 [8] and Iorio [5], the latter ones may be optimistic. Thus, it is difficult to deem the predicted 1pN N -body precession˙ 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 Will [13]; see the discussion in Iorio [5] 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 [5] did not calculate the uncertainty in˙ . In their Table 3, Pitjeva & Pitjev [8] 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 [6], 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 Bepi-Colombo 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 [4] 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 Eq. (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 Schwarzschild-type rates of changes was proposed by Shapiro [11]. It is also widely used in ongoing relativistic tests with geodetic satellites in the Earth's field; see, e.g., Renzetti [10], and references therein for an overview.

Summary and conclusions
Recently, Will [13] calculated a new general relativistic contribution to the Mercury's perihelion advance as large aṡ X 1pN = 220 μas cty −1 arising from an approximated form of the 1pN N -body equations of motion restricted to a hierarchical 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 [4] 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 reductions since it cannot be expressed in terms of a dedicated, solve-for parameter scaling an acceleration different from the aforementioned ones which are routinely modeled. 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 Will [13] by finding an overall precession of˙ X 1pN = 150 μas cty −1 . Iorio [4] 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 Will [13] 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 tȯ I 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 [4]. 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 [5] 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 [5]. 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.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: Because it is a theoretical study.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .