Analytically calculated post-Keplerian range and range-rate perturbations: the solar Lense-Thirring effect and BepiColombo

TL;DR

This paper analytically derives the perturbations in range and range-rate caused by the solar Lense-Thirring effect and J2, assessing their detectability with the BepiColombo mission and implications for solar system dynamics.

Contribution

It provides a general analytical framework for calculating post-Keplerian range and range-rate perturbations due to relativistic and classical effects for arbitrary orbital configurations.

Findings

Sun's gravitomagnetic field affects Earth-Mercury range by ~10 m.

BepiColombo's tracking accuracy is sufficient to detect these effects.

Neglecting Lense-Thirring shifts impacts long-term solar system models.

Abstract

We analytically calculate the time series for the perturbations $\Delta\rho(t),~\Delta\dot\rho(t)$ induced by a general disturbing acceleration $\boldsymbol{A}$ on the mutual range $\rho$ and range-rate $\dot\rho$ of two test particles $\textrm{A},~\textrm{B}$ orbiting the same spinning body. We apply it to the general relativistic Lense-Thirring effect, due to the primary's spin $\boldsymbol{S}$, and the classical perturbation arising from its quadrupole mass moment $J_2$ for arbitrary orbital geometries and orientation of the source's symmetry axis $\boldsymbol{\hat{S}}$. The Earth-Mercury range and range-rate are nominally affected by the Sun's gravitomagnetic field to the $10~\textrm{m},~10^{-3}~\textrm{cm s}^{-1}$ level, respectively, during the extended phase (2026-2028) of the forthcoming BepiColombo mission to Mercury whose expected tracking accuracy is of the order of $\simeq 0.1~\textrm{m},~2\times 10^{-4}~\textrm{cm s}^{-1}$. The competing signatures due to the solar quadrupole $J_2^\odot$, if modelled at the $\sigma_{J_2^\odot}\simeq 10^{-9}$ level of the latest planetary ephemerides INPOP17a, are nearly 10 times smaller than the relativistic gravitomagnetic effects. The position and velocity vectors $\mathbf{r},~\mathbf{v}$ of Mercury and Earth are changed by the solar Lense-Thirring effect by about $10~\textrm{m},~1.5~\textrm{m}$ and $10^{-3}~\textrm{cm s}^{-1},~10^{-5}~\textrm{cm s}^{-1}$, respectively, over 2 yr; neglecting such shifts may have an impact on long-term integrations of the inner solar system dynamics over $\sim\textrm{Gyr}$ timescales.

Equations

$$\Delta \rho =\left({\mathbf{\Delta }\text{r}}_{\text{A}}-{\mathbf{\Delta }\text{r}}_{\text{B}}\right)\cdot \hat{\rho },$$

$$\hat{\rho }=\frac{\left({\mathbf{r}}_{\text{A}}-{\mathbf{r}}_{\text{B}}\right)}{\rho },$$

$${\rho }^2=\left({\mathbf{r}}_{\text{A}}-{\mathbf{r}}_{\text{B}}\right)\cdot \left({\mathbf{r}}_{\text{A}}-{\mathbf{r}}_{\text{B}}\right).$$

$$\mathbf{\Delta }\mathrm{r}=\Delta R\text{leavevmode}\hat{\mathbf{R}}+\Delta T\text{leavevmode}\hat{\mathbf{T}}+\Delta N\text{leavevmode}\hat{\mathbf{N}},$$

$$\Delta R\left(f\right)$$

$$\Delta T\left(f\right)$$

$$\Delta N\left(f\right)$$

$$\hat{\mathbf{R}}=\hat{\mathbf{l}}\cos u+\hat{\mathbf{m}}\sin u,$$

$$\hat{\mathbf{l}}$$

$$\hat{\mathbf{m}}$$

$$\hat{\mathbf{N}}=\sin I\sin \Omega \text{leavevmode}\hat{\text{ı}}-\sin I\cos \Omega \text{leavevmode}\hat{\text{j}}+\cos I\text{leavevmode}\hat{\mathbf{k}}$$

$$\hat{\mathbf{T}}=\hat{\mathbf{N}}\times \hat{\mathbf{R}}.$$

$$r=\frac{a\left(1-e^2\right)}{1+e\cos f},$$

$$\mathbf{r}=r\text{leavevmode}\hat{\mathbf{R}},$$

$$\Delta \kappa \left(f\right)={\int }_{f_0}^f\frac{\text{d}\kappa }{\text{d}t}\frac{\text{d}t}{\text{d}{f}^{{}^{\prime }}}d{f}^{{}^{\prime }},\text{leavevmode}\kappa =a,\text{leavevmode}e,\text{leavevmode}I,\text{leavevmode}\Omega ,\text{leavevmode}\omega ,$$

$$\frac{\text{d}t}{\text{d}f}=\frac{{\left(1-e^2\right)}^{3/2}}{n_{\mathrm{b}}{\left(1+e\cos f\right)}^2},$$

$$\Delta \mathcal{M}\left(f\right)=\Delta \eta \left(f\right)+{\int }_{t_0}^t\Delta n_{\mathrm{b}}\left(t^{{}^{\prime }}\right)d{t}^{{}^{\prime }},$$

$${\int }_{t_0}^t\Delta n_{\mathrm{b}}\left(t^{{}^{\prime }}\right)d{t}^{{}^{\prime }}=-\frac{3}{2}\frac{n_{\mathrm{b}}}{a}{\int }_{f_0}^f\Delta a\left(f^{{}^{\prime }}\right)\frac{dt}{d{f}^{{}^{\prime }}}d{f}^{{}^{\prime }},$$

$$\Delta \dot{\rho }=\left({\mathbf{\Delta }\text{v}}_{\text{A}}-{\mathbf{\Delta }\text{v}}_{\text{B}}\right)\cdot \hat{\rho }+\left({\mathbf{\Delta }\text{r}}_{\text{A}}-{\mathbf{\Delta }\text{r}}_{\text{B}}\right)\cdot {\hat{\rho }}_{\text{v}},$$

$${\hat{\rho }}_{\text{v}}=\frac{\left({\text{v}}_{\text{A}}-{\text{v}}_{\text{B}}\right)-\dot{\rho }\hat{\rho }}{\rho },$$

$$\dot{\rho }=\left({\text{v}}_{\text{A}}-{\text{v}}_{\text{B}}\right)\cdot \hat{\rho }.$$

$$\mathbf{\Delta }\mathrm{v}=\Delta {\text{v}}_R\text{leavevmode}\hat{\mathbf{R}}+\Delta {\text{v}}_T\text{leavevmode}\hat{\mathbf{T}}+\Delta {\text{v}}_N\text{leavevmode}\hat{\mathbf{N}}.$$

$$\Delta {\mathrm{v}}_R\left(f\right)$$

$$-\frac{n_{\mathrm{b}}{a}^2}{r\left(f\right)}\sqrt{1-e^2}\left[\cos I\Delta \Omega \left(f\right)+\Delta \omega \left(f\right)\right],$$

$$\Delta {\mathrm{v}}_T\left(f\right)$$

$$\Delta {\mathrm{v}}_N\left(f\right)$$

$$\mathbf{v}=\sqrt{\frac{GM}{a\left(1-e^2\right)}}\left[-\hat{\mathbf{P}}\sin f+\hat{\mathbf{Q}}\left(e+\cos f\right)\right]$$

$$\hat{\mathbf{P}}=\hat{\mathbf{l}}\cos \omega +\hat{\mathbf{m}}\sin \omega ,$$

$$\hat{\mathbf{Q}}=-\hat{\mathbf{l}}\sin \omega +\hat{\mathbf{m}}\cos \omega ;$$

$$\mathbf{r}=r\left(\hat{\mathbf{P}}\cos f+\hat{\mathbf{Q}}\sin f\right).$$

$$\Delta {\mathcal{M}}^{\text{GE}}\left(f\right)=C_1\left\{\arctan \left[\frac{\left(-1+e\right)\tan \left(\frac{f}{2}\right)}{\sqrt{1-e^2}}\right]-\arctan \left[\frac{\left(-1+e\right)\tan \left(\frac{f_0}{2}\right)}{\sqrt{1-e^2}}\right]\right\}+C_2,$$

$$C_1$$

$$+3e\left(\left(56+e^2\left(24-31\nu \right)-24\nu \right)\cos f_0+e\left(4\left(5-4\nu \right)\cos 2f_0-e\nu \cos 3f_0\right)\right)),$$

$$C_2$$

$$+e^2\left(-84+76\nu \right)+3e(\left(8\left(-7+3\nu \right)+e^2\left(-24+31\nu \right)\right)\cos f_0+$$

Full text

keywords gravitation–celestial mechanics–Sun: rotation–space veichles

1 Introduction

Let us consider a pair of test particles $\textrm{A},\leavevmode\nobreak\ \textrm{B}$ following two generally different elliptical paths around the same spinning, oblate massive body of mass $M$, equatorial radius $\mathcal{R}$, quadrupole mass moment $J_{2}$ and angular momentum $\boldsymbol{S}$, whose symmetry axis is arbitrarily oriented in space; they could typically be two planets orbiting a star, two man-made spacecraft circling a planet or a satellite (Wolff, 1969) like the GRACE (Tapley et al., 2004) and GRAIL (Zuber et al., 2013) space missions, or a subset of constellations of more than two spacecraft like, e.g., the planned eLISA (Amaro-Seoane et al., 2013). The two-body range $\rho$, i.e. the mutual distance $\rho$ between A and B, and its temporal rate of change $\dot{\rho}$, i.e. the range-rate, contain a wealth of information about the evolution of the system itself and the properties of the common primary’s gravitational field in which $\textrm{A},\leavevmode\nobreak\ \textrm{B}$ move. Think about, e.g., the past and future long-term evolution of our solar system which has been-and will perhaps be-studded with collisions among planetesimals and planets themselves (Canup & Asphaug, 2001; Laskar & Gastineau, 2009), or, to much shorter timescales, the accurate mapping of the Earth’s gravity field performed by GRACE and its planned follow-on GRACE-FO (Loomis, Nerem & Luthcke, 2012). At intermediate time-scales, geocentric range and range-rate measurements to spacecraft orbiting some major bodies of our solar system are able to improve, among other things, their orbit determination (Konopliv et al., 2011; Verma et al., 2014; Folkner, Jacobson & Jones, 2015). As a consequence, intersatellite and interplanetary tracking through electromagnetic waves of different frequencies have reached an impressive accuracy level (Bender et al., 2003; Kim & Lee, 2009; Sheard et al., 2012; Iess, Asmar & Tortora, 2009). One of the main application of accurate ranging is testing the general theory of relativity, whose current status was recently overviewed by, e.g., Debono & Smoot (2016). In this paper, we analytically calculate the mutual two-body range and range-rate perturbations $\Delta\rho,\leavevmode\nobreak\ \Delta\dot{\rho}$ induced by the general relativistic gravitomagnetic111It is generated by mass-energy currents encoded in the $g_{0j},\leavevmode\nobreak\ j=1,2,3$ components of the spacetime metric tensor $g_{\mu\nu},\leavevmode\nobreak\ \mu,\nu=0,1,2,3$ (Thorne, 1986). It is believed to play important roles in astrophysical scenarios involving spinning black holes (Thorne, 1988). At present, its main experimental confirmation relies upon the measurement performed by the GP-B probe in the Earth’s gravitational field, accurate to $19$ per cent (Everitt et al., 2011); for other performed, ongoing and proposed attempts involving the analysis of the orbital motions of natural and artificial test particles around spinning planets and stars, see, for example, (Iorio et al., 2011; Ciufolini et al., 2013; Renzetti, 2013). field of the common primary for an arbitrary orientation of its angular momentum and without making a-priori simplifying assumptions on the eccentricity of both the orbiting test particles. This choice is motivated by reasons of generality and in view of the fact that the spatial orientation of the angular momentum of astronomical bodies which could be considered as potentially interesting sources of testable gravitomagnetic fields is not always known with sufficient accuracy. Furthermore, the coordinate systems used in several practical astronomical data analyses do not generally adopt the primary’s equator as their fundamental plane. Moreover, over multiyear data analyses, the angular momentum of the source of the gravitomagnetic field often changes its position with respect to that originally occupied at the reference epoch of the coordinate system used because of a variety of more or less accurately known physical phenomena; also such time-dependent evolution is necessarily known with a given level of uncertainty. Our general approach is fully able to cope with such issues allowing for suitably designing dedicated scenarios and performing sensitivity analyses in view of, e.g., forthcoming space-based missions like BepiColombo at Mercury (Balogh et al., 2007; Schettino & Tommei, 2016). This is important either if one is considering to put directly to the test the Lense-Thirring effect via range and range-rate measurements, and if, conversely, one trusts general relativity also as far as the gravitomagnetic field is concerned and uses it as a potential tool to measure the angular momentum of an astronomical body. Even if one is interested in other dynamical effects, like, e.g., characterizing the Newtonian multipolar structure of the gravitational potential of its source, the Lense-Thirring effect must be taken into account to avoid a priori biased results, as it would be the case for the Sun’s oblateness (Folkner et al., 2014). For the sake of completeness, we deal also with the Schwarzschild-type gravitoelectric range and range-rate perturbations. Almost all the results required to build them can be already found in the literature except the short-term variation of the mean anomaly, which is correctly calculated in the present paper. It is important to point out that, in view of its generality and lack of a priori restrictions on the orbital geometries involved, our strategy can be straightforward extended to any other perturbing acceleration induced, e.g., by modified models of gravity as well.

The plan of the paper is as follows. In Section 2, we outline our general approach to analytically calculate the perturbations induced by any small post-Keplerian disturbing acceleration on the mutual range and range-rate of a pair of test particles orbiting a common primary. It is applied to the Schwarzschild-type acceleration in Section 3.1 and to the Lense-Thirring effect in Section 3.2, where exact analytical expressions for the building blocks of the range and range-rate gravitomagnetic shifts are explicitly displayed. Section 4 deals with the range and range-rate perturbations due to the first even zonal harmonic $J_{2}$ of the multipolar expansion of the primary’s gravitational field since it is a major source of systematic bias of dynamical origin to be accounted for in any realistic preliminary sensitivity analysis. In view of their extreme cumbersomeness in the case of arbitrary spatial orientation of the primary’s spin axis, it is not possible to explicitly show the analytical results obtained by applying the general method of Section 2 to this specific disturbance. The shifts treated in Sections 3.2 to 4 are due only to the orbital motions of $\textrm{A},\leavevmode\nobreak\ \textrm{B}$ in the field of their primary. We will show that the effects of the propagation of the electromagnetic waves in its gravitomagnetically deformed spacetime are negligible in the specific scenario treated in Section 5 which deals with the geocentric Hermean222From ‛Ερμῆς (‘Hermes’), corresponding to the ancient Roman deity Mercury range and range-rate during the planned two-years extended mission of BepiColombo orbiting Mercury. In Section 6, we apply the results of Section 3.2 to evaluate the corrections to the Hermean and terrestrial state vectors in view of possible consequences of neglecting the Sun’s gravitomagnetic field in the so-far performed long-term numerical integrations of the solar system dynamics over future eons. Our finding and conclusions are summarized in Section 7.

2 General calculational method of the range and range-rate perturbations

According to (Cheng, 2002), the perturbation $\Delta\rho$ of the mutual range $\rho$ of two test particles $\textrm{A},\leavevmode\nobreak\ \textrm{B}$ orbiting the same central body is

[TABLE]

where

[TABLE]

and

[TABLE]

In turn, the perturbation $\mathbf{\Delta\mathrm{r}}$ experienced by the position vector $\mathbf{r}$ of any of the two bodies $\textrm{A},\leavevmode\nobreak\ \textrm{B}$ is

[TABLE]

In Equation (4), the instantaneous radial, transverse and normal perturbations $\Delta R,\leavevmode\nobreak\ \Delta T,\leavevmode\nobreak\ \Delta N$ of the position vector $\mathbf{r}$, expressed in terms of the osculating semimajor axis $a$, eccentricity $e$, inclination $I$, longitude of the ascending node $\Omega$, argument of pericenter $\omega$, mean anomaly $\mathcal{M}$, are (Casotto, 1993)

[TABLE]

In Equations (5) to (7), $f$ is the true anomaly, assumed as fast variable encoding the time dependence (see Equation (47)), and $u=\omega+f$ is the argument of latitude. The radial unit vector $\mathbf{\hat{R}}$ entering Equation (4) can be expressed as (Brumberg, 1991)

[TABLE]

where the unit vectors ${\boldsymbol{\hat{l}}},\leavevmode\nobreak\ {\boldsymbol{\hat{m}}}$ are defined as (Brumberg, 1991)

[TABLE]

The vector ${\boldsymbol{\hat{l}}}$ is directed along the line of the nodes toward the ascending node, while ${\boldsymbol{\hat{m}}}$ is directed transversely to the line of the nodes in the orbital plane; $\boldsymbol{\hat{\textbf{\i}}},\leavevmode\nobreak\ \boldsymbol{\hat{\j}},\leavevmode\nobreak\ \mathbf{\hat{k}}$ are the usual unit vectors spanning the reference $x,\leavevmode\nobreak\ y,\leavevmode\nobreak\ z$ axes of the coordinate system adopted. The normal unit vector $\mathbf{\hat{N}}$ in Equation (4) is (Brumberg, 1991)

[TABLE]

Thus, the transverse unit vector $\mathbf{\hat{T}}$ appearing in Equation (4) can straightforwardly be obtained as (Brumberg, 1991)

[TABLE]

Both Equations (1) to (3) and Equations (5) to (7) must be evaluated onto the unperturbed Keplerian ellipse

[TABLE]

assumed as reference orbit; thus the position vector of any of $\textrm{A},\leavevmode\nobreak\ \textrm{B}$, entering Equations (2) to (3), is

[TABLE]

where $r$ is given by Equation (13). Furthermore, the instantaneous shifts $\Delta\kappa\left(f\right),\leavevmode\nobreak\ \kappa=a,\leavevmode\nobreak\ e,\leavevmode\nobreak\ I,\leavevmode\nobreak\ \Omega,\leavevmode\nobreak\ \omega$ entering Equation (4) are to be calculated as

[TABLE]

by taking $d\kappa/dt,\leavevmode\nobreak\ \kappa=a,\leavevmode\nobreak\ e,\leavevmode\nobreak\ I,\leavevmode\nobreak\ \Omega,\leavevmode\nobreak\ \omega$ from the right-hand-sides of the usual Gauss equations for the variation of the elements (Brumberg, 1991; Milani, Nobili & Farinella, 1987; Soffel, 1989; Bertotti, Farinella & Vokrouhlický, 2003), evaluated onto Equation (13), and with

[TABLE]

in which $n_{\rm b}=\sqrt{GMa^{-3}}$ is the Keplerian mean motion; $G$ is the Newtonian constant of gravitation. The instantaneous change $\Delta\mathcal{M}\left(f\right)$ of the mean anomaly $\mathcal{M}=n_{\rm b}\left(t-t_{p}\right)$, where $t_{p}$ is the time of passage at pericenter, must be evaluated as

[TABLE]

where $\Delta\eta\left(f\right)$ is the instantaneous change of the mean anomaly at epoch $\eta$, worked out with the corresponding Gauss equation (Brumberg, 1991; Milani, Nobili & Farinella, 1987; Soffel, 1989; Bertotti, Farinella & Vokrouhlický, 2003), and

[TABLE]

Depending on the specific perturbation at hand, the explicit calculation of Equation (18) may turn out rather cumbersome.

It is possible to analytically calculate also the shift $\Delta\dot{\rho}$ experienced by the range-rate $\dot{\rho}$ as follows (Cheng, 2002)

[TABLE]

where

[TABLE]

with

[TABLE]

The perturbation $\mathbf{\Delta\mathrm{v}}$ of the velocity vector $\mathbf{v}$ of any of the two bodies $\textrm{A},\leavevmode\nobreak\ \textrm{B}$ can be written as

[TABLE]

In Equation (22), the instantaneous radial, transverse and normal perturbations $\Delta\textrm{v}_{R},\leavevmode\nobreak\ \Delta\textrm{v}_{T},\leavevmode\nobreak\ \Delta\textrm{v}_{N}$ of the position vector $\mathbf{v}$, expressed in terms of the osculating orbital elements, are (Casotto, 1993)

[TABLE]

Since, also in this case, the Keplerian ellipse has to be adopted as unperturbed reference trajectory,

[TABLE]

must be used in Equations (20) to (21). The unit vectors $\boldsymbol{\hat{P}},\leavevmode\nobreak\ \boldsymbol{\hat{Q}}$ in Equation (26) are defined as (Brumberg, 1991)

[TABLE]

while $\boldsymbol{\hat{P}}$ is directed along the line of the apsides toward the pericenter, $\boldsymbol{\hat{Q}}$ lies transversely to the line of the apsides in the orbital plane. Incidentally, it turns out that also the position vector $\mathbf{r}$ can be expressed in terms of $\boldsymbol{\hat{P}},\leavevmode\nobreak\ \boldsymbol{\hat{Q}}$ as (Brumberg, 1991)

[TABLE]

3 The general relativistic range and range-rate perturbations

3.1 The gravitoelectric Schwarschild-like effect

As far as the post-Newtonian gravitoelectric333In our solar system, it induces the formerly anomalous, time-honored perihelion precession of Mercury of $\dot{\omega}^{\textrm{GE}}_{\mercury}=42.98\leavevmode\nobreak\ \textrm{arcsec\leavevmode\nobreak\ cty}^{-1}$ (Nobili & Will, 1986). acceleration due to the static part of the primary’s field is concerned, most of the building blocks required to calculate its range and range-rate perturbations are available in the literature. Indeed, the instantaneous variations $\Delta\kappa^{\textrm{GE}}\left(f\right),\leavevmode\nobreak\ \kappa=a,\leavevmode\nobreak\ e,\leavevmode\nobreak\ I,\leavevmode\nobreak\ \Omega,\leavevmode\nobreak\ \omega$ can be found in Equations (A2.78b) to (A2.78d) of Soffel (1989, p. 178). The correct calculation of the variation of the mean anomaly requires more care. According to Equations (17) to (18), its instantaneous change turns out to be

[TABLE]

with

[TABLE]

In Equations (3.1) to (3.1), the dimensionless mass parameter $\nu=mM\left(M+m\right)^{-2}$, where $m$ is the mass of any of the two objects A, B, vanishes in the test particle limit. Equations (A2.78e) to (A2.78f) of Soffel (1989, p. 178) allow to obtain the instantaneous shift of the mean anomaly in terms of the three anomalies $f,\leavevmode\nobreak\ E,\leavevmode\nobreak\ \mathcal{M}$; instead, only the true anomaly $f$ enters our Equation (30). See also Equations (3.1.102) to (3.1.107) of Brumberg (1991, p. 93). Let us remark that Equation (30) is exact in $e$.

3.2 The gravitomagnetic Lense-Thirring effect

The Lense-Thirring acceleration ${\boldsymbol{A}}_{\textrm{LT}}$, written in harmonic, post-Newtonian coordinates, can be found in several references; see, e.g., Brumberg (1991); Soffel (1989); Petit & Luzum (2010). Its radial, transverse and normal components, for a generic orientation of the spin axis unit vector ${\boldsymbol{\hat{S}}}$, i.e. for ${\boldsymbol{\hat{S}}}\neq\mathbf{\hat{k}}$, can be retrieved, e.g., in Iorio (2017); they are

[TABLE]

where $c$ is the speed of light. Thus, the shifts $\Delta R_{\textrm{LT}},\leavevmode\nobreak\ \Delta T_{\textrm{LT}},\leavevmode\nobreak\ \Delta N_{\textrm{LT}}$ of the radial, transverse and normal components of the position vector $\mathbf{r}$ turn out to be

[TABLE]

where

[TABLE]

From Equation (39), it can be noted that $\boldsymbol{w}$ is not a unit vector. The shifts $\Delta\textrm{v}^{\textrm{LT}}_{R},\leavevmode\nobreak\ \Delta\textrm{v}^{\textrm{LT}}_{T},\leavevmode\nobreak\ \Delta\textrm{v}^{\textrm{LT}}_{N}$ of the radial, transverse and normal components of the velocity vector $\mathbf{v}$ are

[TABLE]

where

[TABLE]

According to Equation (43), $\boldsymbol{d}$ is not a unit vector.

4 The oblateness-induced range and range-rate perturbations

If the main target is testing some predictions of general relativity, one of the major sources of systematic uncertainties is represented by the departures from spherical symmetry of the primary. They are usually modeled as an expansion in multipoles of its gravitational potential. A class of such multipolar coefficients cause long-term Newtonian orbital perturbations whose nominal size is usually orders of magnitude larger than the non-Newtonian effects of interest. Conversely, the determination of the multipoles themselves can be one of the main scientific goals of a given mission; in this case, it is the Lense-Thirring effect that, if not explicitly modeled, can bias the outcome of the experiment, as recently recognized for the Sun by Folkner et al. (2014) on the basis of Iorio et al. (2011).

As far as the disturbing acceleration induced by the quadrupole mass moment $J_{2}$ of the primary is concerned, its radial, transverse and normal components for an arbitrary orientation of the unit vector ${\boldsymbol{\hat{S}}}$ of its symmetry axis are Iorio (2017)

[TABLE]

where $\mu=GM$ is the gravitational parameter of the primary.

Without recurring to some a priori approximations about the orbital configuration of the test particle and the orientation of the primary’s spin axis, very cumbersome analytical expressions are obtained for $\Delta R_{J_{2}},\leavevmode\nobreak\ \Delta T_{J_{2}},\leavevmode\nobreak\ \Delta N_{J_{2}},\leavevmode\nobreak\ \Delta\textrm{v}_{R}^{J_{2}},\leavevmode\nobreak\ \Delta\textrm{v}_{T}^{J_{2}},\leavevmode\nobreak\ \Delta\textrm{v}_{N}^{J_{2}}$ from Equations (44) to (46); thus, they cannot be explicitly displayed. In particular, the calculation of Equation (18) entering the total shift of the mean anomaly in Equation (17) turns out to be particularly unwieldy. They can be conveniently simplified depending on the specific scenario at hand by expanding them in powers of $e$ to the desired level of accuracy.

5 The BepiColombo range and range-rate perturbations

As a concrete application of our results in Sections 3.2 to 4, let us consider the geocentric Hermean range and range-rate during the expected extended phase (2026 March 14-2028 May 1) of the BepiColombo mission to Mercury (Benkhoff et al., 2010, 2017); according to http://sci.esa.int/bepicolombo/47346-fact-sheet/, its launch is currently scheduled in October 2018. It should greatly improve, among other things, the accuracy of the orbital determination of the small rocky planet to the $\sigma_{\rho}\simeq 0.1\leavevmode\nobreak\ \textrm{m},\leavevmode\nobreak\ \sigma_{\dot{\rho}}\simeq 2\times 10^{-4}\leavevmode\nobreak\ \textrm{cm\leavevmode\nobreak\ s}^{-1}$ level for the range and the range-rate, respectively (Milani et al., 2010). Here, we will not deal in detail with the gravitoelectric effect of Section 3.1 since it has been already treated numerically in the literature (Milani et al., 2002; Ashby, Bender & Wahr, 2007; Milani et al., 2010; Schettino & Tommei, 2016; Imperi, Iess & Mariani, 2018). Recently, more and more extended portions of the data record of the NASA MESSENGER mission to Mercury, ended on April 30, 2015, have started to be analyzed in order to look at the Lense-Thirring effect explicitly by modeling it (Park, Folkner & Konopliv, 2015; Park et al., 2017; Genova et al., 2018); the first preliminary results evidenced a strong correlation of the gravitomagnetic signature with the solar oblateness, thus limiting the accuracy of a possible direct detection of the relativistic effect to about $\sim 20-25\%$. The perspectives of a direct detection of the solar gravitomagnetic field opened up by the recent advances in the field of the planetary ephemerides were pointed out by the present author more than a decade ago; see, e.g., Iorio et al. (2011) and references therein.

Before proceeding further, it is advisable to briefly review some basic features of the confrontation between theory and observations in Relativistic Celestial Mechanics (RCM) (Brumberg, 2010b, a; Kopeikin, Efroimsky & Kaplan, 2011; Brumberg, 2013). Such a task implies solving not only the dynamics of the specific problem at hand, i.e. the equations of motion of the massive bodies involved, but also the equations of propagation of the electromagnetic waves and the description of the observational procedures (the kinematical part of RCM). Both parts should be investigated in the same coordinates to exclude unphysical coordinate–dependent spurious effects and to present the results in terms of measurable quantities. Indeed, contrary to the Newtonian case, the issue of coordinate–dependent quantities is the main qualitative new feature of RCM. Among the main possibilities to overcome such a problem emerged in the literature over the past decades, in 1991 it was pragmatically adopted by the International Astronomical Union (IAU) the approach of forgetting about the general relativistic arbitrariness in the coordinate conditions, and to use one specific type of coordinates-the harmonic ones-for both the dynamical and the kinematical parts of RCM once for all; see Soffel et al. (2003) for the more recent IAU 2000 resolutions updating the earlier IAU 1991 ones. In the case of interplanetary ranging experiments in the Solar System, it is customary to adopt some realization of a Barycentric Reference System (BRS) along with some suitable Barycentric Coordinate Time (TCB) as global coordinate system. Suitable and well established spacetime coordinate transformations from the specific local Planetocentric Reference Systems (PRSs) of the major bodies involved in the experiment under examination to the BRS allow to avoid the insurgence of spurious terms like in the nowadays outdated problem of the unphysical effects plaguing the Earth-Moon range when calculated in the BRS instead of some realization of a Geocentric Reference System (GRS)444In analyzing the data provided by the Earth-based Lunar Laser Ranging (LLR) technique, a GRS is more suited than a BRS for the reasons explained in, e.g., Brumberg (2010b). (Brumberg, 2010b). For details of the computation of observables for BepiColombo and the required coordinate transformations, see, e.g., Section 3.1 of Schettino & Tommei (2016).

In Figure 1, the nominal Earth-Mercury range and range-rate signatures induced by the Sun’s angular momentum via the Lense-Thirring effect (red) and its quadrupole mass moment (blue) are depicted according to our analytical results in Sections 3.2 to 4.

In order to produce time series, we adopted the following expansion of the planetary true anomaly $f$ in terms of the mean anomaly $\mathcal{M}$ (Brouwer & Clemence, 1961, p. 77)

[TABLE]

where $J_{k}\left(se\right)$ is the Bessel function of the first kind of order $k$ and $s_{\textrm{max}},\leavevmode\nobreak\ j_{\textrm{max}}$ are some values of the summation indexes $s,\leavevmode\nobreak\ j$ set by the desired accuracy level. It can be noted that the $J_{2}^{\odot}$-induced signals are much larger than the gravitomagnetic ones since their amplitudes can be as large as $300\leavevmode\nobreak\ \textrm{m}$ and $0.03\leavevmode\nobreak\ \textrm{cm\leavevmode\nobreak\ s}^{-1}$, while the general relativistic ones are as little as $10\leavevmode\nobreak\ \textrm{m}$ and $0.0010\leavevmode\nobreak\ \textrm{cm\leavevmode\nobreak\ s}^{-1}$. Such figures for the predicted Lense-Thirring range and range-rate shifts fall well within the previously mentioned improvements in the orbit determination of Mercury (Milani et al., 2010). However, the Sun’s quadrupole field is routinely included in the dynamical models of the current planetary ephemerides. Figure 2 shows the mismodeled $J_{2}^{\odot}$ signature according to (Park et al., 2017) $\sigma_{J_{2}^{\odot}}=9\times 10^{-9}$ recently inferred by processing the ranging data of the MESSENGER mission.

Now, the size of the quadrupolar signatures is at the same level of the gravitomagnetic ones; the picture highlights the different temporal patterns characterizing the two effects under consideration. Finally, Figure 3 depicts the case in which a much smaller uncertainty in the Sun’s first even zonal is assumed: $\sigma_{J_{2}^{\odot}}=1\times 10^{-9}$ obtained with the recent INPOP17a ephemerides Viswanathan et al. (2017) in a global fit to an almost centennial record of data of several types including also, among other things, ranging to MESSENGER.

It turns out that the Lense-Thirring signatures are neatly predominant with respect to the residual traces left by the mismodeling in $J_{2}^{\odot}$ which are about 10 times smaller than the gravitomagnetic effects. Recently, Genova et al. (2018) obtained, among other things, $\sigma_{J_{2}^{\odot}}=2.2\times 10^{-9}$ by analyzing the MESSENGER data. It is important to remark that simulations of the BepiColombo mission, performed so far without modeling the Lense-Thirring effect, point towards an accuracy level in determining the Sun’s quadrupole of the order of $\sigma_{J_{2}^{\odot}}\simeq 4.1-5.5\times 10^{-10}$ (Schettino & Tommei, 2016; Imperi, Iess & Mariani, 2018). The corresponding plots are displayed in Figure 4; the maximum values of the mismodeled quadrupolar signals are now of the order of just $0.6\leavevmode\nobreak\ \textrm{m},\leavevmode\nobreak\ 4\times 10^{-5}\textrm{cm\leavevmode\nobreak\ s}^{-1}$.

Our analytical calculation is based on a perturbative approach in terms of the osculating Keplerian orbital elements, which are used in celestial mechanics, applied to the Lense-Thirring acceleration ${\boldsymbol{A}}_{\textrm{LT}}$ written in the standard harmonic post-Newtonian coordinates555According to http://iaaras.ru/en/dept/ephemeris/epm/2017/, ${\boldsymbol{A}}_{\textrm{LT}}$ was explicitly included in the latest version EPM2017 of the EPM ephemerides produced by the team led by E.V. Pitjeva. (Brumberg, 1991; Soffel, 1989; Petit & Luzum, 2010). In order to check its validity, we numerically integrated the equations of motion of all the major bodies of the Solar System with and without the solar gravitomagnetic acceleration ${\boldsymbol{A}}_{\textrm{LT}}$ over the extended BepiColombo mission starting from the same initial conditions, and produced numerical time series $\Delta\rho_{\textrm{num}},\leavevmode\nobreak\ \Delta\dot{\rho}_{\textrm{num}}$ for the gravitomagnetic Earth-Mercury range and range-rate shifts of orbital origin. Then, we compared them to the corresponding analytically worked out time series $\Delta\rho_{\textrm{anal}},\leavevmode\nobreak\ \Delta\dot{\rho}_{\textrm{anal}}$ displayed in Figures 1 to 4. Figure 5 depicts their differences; they are well below the expected experimental accuracy level since they amount to about $\left|\Delta\rho_{\textrm{anal}}-\Delta\rho_{\textrm{num}}\right|\lesssim 5\times 10^{-5}\leavevmode\nobreak\ \textrm{m},\leavevmode\nobreak\ \left|\Delta\dot{\rho}_{\textrm{anal}}-\Delta\dot{\rho}_{\textrm{num}}\right|\lesssim 1\times 10^{-4}\leavevmode\nobreak\ \textrm{cm\leavevmode\nobreak\ s}^{-1}$, respectively. Furthermore, Figure 5 demonstrates pragmatically and effectively that using the osculating Keplerian orbital elements to perform our analytical gravitomagnetic calculation as in Section 3.2 did not introduce any spurious, unphysical harmonics, contrary to what could, in principle, be argued on the basis of what happened in the literature for the post-Newtonian gravitoelectric field in the case of the two-body problem and its different orbital parameterizations (Klioner & Kopeikin, 1994). Moreover, to the best of the author’s knowledge, no other orbital parameterizations than the osculating elements have been used so far in the case of the Lense-Thirring effect; see, e.g., Damour & Schafer (1988) and references therein.

As far as the gravitomagnetic propagation delay $\Delta t_{\textrm{LT}}$ is concerned, it turns out to be negligible in the present scenario. Indeed, it can be shown that is proportional to (Kopeikin, 1997; Wex & Kopeikin, 1999; Ciufolini et al., 2003)

[TABLE]

where $\mathcal{F}$ is a geometric factor depending on the mutual orientation of the primary’s spin axis $\boldsymbol{\hat{S}}$ and the position vectors ${\mathbf{r}}_{\textrm{A}},\leavevmode\nobreak\ {\mathbf{r}}_{\textrm{B}}$ of the orbiting bodies. For the Sun and Mercury, Equation (48) yields a time delay as little as $\left|\Delta t_{\textrm{LT}}\right|\lesssim 5\times 10^{-14}\leavevmode\nobreak\ \textrm{s}$, corresponding to a range shift of the order of just $\simeq 1.6\times 10^{-5}\leavevmode\nobreak\ \textrm{m}$. Furthermore, while the range perturbation of orbital origin is cumulative in time, the propagation delay of Equation (48) is periodic. The usual Shapiro time delay caused by the post-Newtonian gravitoelectric field of the Sun

[TABLE]

where $\gamma$ is the parameter of the Parameterized Post-Newtonian (PPN) formalism accounting for the spatial curvature, yields an Earth-Mercury range signature which is nominally much larger than the gravitomagnetic one of orbital origin (see the upper panel of Figure 6). On the other hand, it is routinely included in the data reduction softwares, and its mismodeling due to the uncertainty in $\gamma$ is small enough not to represent a problem, as shown by the lower panel of Figure 6 obtained for $\sigma_{\gamma}=2.3\times 10^{-5}$ (Bertotti, Iess & Tortora, 2003); indeed, it is apparent that it could not cancel out the Lense-Thirring Earth-Mercury range perturbation also because its temporal pattern is quite different.

In conclusion, it is not possible that the Lense-Thirring range shift, arising from the equations of motion of massive bodies, can be canceled by the range perturbations induced by the post-Newtonian delays in the propagation of the electromagnetic waves. Indeed, the latter ones have different temporal patterns, and their nominal (gravitomagnetic) and mismodelled (gravitoelectric) magnitudes are much smaller than the gravitomagnetic signal of interest.

6 Gravitomagnetic corrections to the state vector of Mercury

In Figures 7 to 8, we plot the Lense-Thirring corrections $\Delta x_{\textrm{LT}},\leavevmode\nobreak\ \Delta y_{\textrm{LT}},\leavevmode\nobreak\ \Delta z_{\textrm{LT}},\leavevmode\nobreak\ \Delta\dot{x}_{\textrm{LT}},\leavevmode\nobreak\ \Delta\dot{y}_{\textrm{LT}},\leavevmode\nobreak\ \Delta\dot{z}_{\textrm{LT}}$ to the position and velocity vectors $\mathbf{r},\leavevmode\nobreak\ \mathbf{v}$ of Mercury and Earth over a 2-yr time span; for the sake of definiteness, we adopt the timeframe of the extended mission of BepiColombo.

It turns out that the components of the position vectors are shifted by about $10\leavevmode\nobreak\ \textrm{m}$ (Mercury) and $1.5\leavevmode\nobreak\ \textrm{m}$ (Earth), while the velocity’s components are changed by about $0.001\leavevmode\nobreak\ \textrm{cm\leavevmode\nobreak\ s}^{-1}$ (Mercury) and $0.00003\leavevmode\nobreak\ \textrm{cm\leavevmode\nobreak\ s}^{-1}$ (Earth).

Such results may be used as possible suggestions to consider also the solar gravitomagnetic field in future accurate long-term integrations of the solar system aimed to test its stability over the eons. Indeed, it would be interesting to quantitatively check if neglecting the Lense-Thirring effect may have some appreciable impact on the outcomes of those existing studies predicting non-zero probabilities of interplanetary collisions or ejections. To the post-Newtonian level, only the gravitostatic, Schwarzschild-like acceleration was included so far in the existing literature; see, e.g., the Appendix in Zeebe (2015b). (Laskar & Gastineau, 2009; Laskar et al., 2011; Lithwick & Wu, 2011; Boué, Laskar & Farago, 2012; Batygin, Morbidelli & Holman, 2015; Zeebe, 2015a, b, 2017). In view of Figures 7 to 8, it does not seem unreasonable since, as pointed out in Laskar (1989); Zeebe (2015b), a difference in initial position of $1\leavevmode\nobreak\ \textrm{cm}$ grows to about $1\leavevmode\nobreak\ \textrm{au}$ after $90-150\leavevmode\nobreak\ \textrm{Myr}$. Furthermore, in the analysis by Zeebe (2015b), the Hermean initial radial distance was offset by $1.75\leavevmode\nobreak\ \textrm{mm}$ between every two adjacent orbits while the largest overall offset was $2.80\leavevmode\nobreak\ \textrm{m}$. Such a figure is smaller than the typical Lense-Thirring shift after just two years, as per Figure 7.

7 Summary and conclusions

In view of their wide application in practical investigations about several astronomical and astrophysical scenarios of interest for tests of fundamental physics, we looked at the mutual range $\rho$ and range-rate $\dot{\rho}$ of two test particles $\textrm{A},\leavevmode\nobreak\ \textrm{B}$ orbiting a common rotating primary. We set up an approach to explicitly calculate the perturbations $\Delta\rho,\leavevmode\nobreak\ \Delta\dot{\rho}$ due to a generic disturbing acceleration $\boldsymbol{A}$ for which an explicit analytical expression is available. We applied our computational strategy to the post-Newtonian and Newtonian effects of lowest order induced by the rotation of the central body: the general relativistic Lense-Thirring field, generated by the angular momentum $\boldsymbol{S}$ of the primary, and the classical perturbation caused by the first even zonal harmonic $J_{2}$ of its non-spherical potential. Our analytical results, which are particularly cumbersome in the case of the source’s oblateness, are completely general since they do not rely upon any a priori simplifying assumptions pertaining both the particles’ orbital configurations and the orientation of the primary’s symmetry axis. Thus, they can be applied to any system of intercommunicating probes designing dedicated missions, like, e.g., a GRACE-type tandem orbiting some giant planet of our solar system, performing sensitivity analyses, reinterpreting existing data, and looking also at the long-term dynamics of, say, the inner solar system. We looked also at the general relativistic Schwarzschild-like effect by calculating exactly the corresponding full instantaneous shift of the mean anomaly; the variations of the other orbital elements can be found in the existing literature.

As a practical application of our analytical calculation, we considered the geocentric range and range-rate of Mercury during the planned extended phase of the forthcoming BepiColombo mission, to be launched in late 2018, which, among other things, should notably improve the Hermean ephemerides. It turned out that the expected nominal Lense-Thirring perturbations can reach the $10\leavevmode\nobreak\ \textrm{m},\leavevmode\nobreak\ 1\times 10^{-3}\leavevmode\nobreak\ \textrm{cm\leavevmode\nobreak\ s}^{-1}$ level, well within the tracking accuracy of BepiColombo which is of the order of $\sigma_{\rho}\simeq 0.1\leavevmode\nobreak\ \textrm{m},\leavevmode\nobreak\ \sigma_{\dot{\rho}}\simeq 2\times 10^{-4}\leavevmode\nobreak\ \textrm{cm\leavevmode\nobreak\ s}^{-1}$. The competing signatures induced by the Sun’s quadrupole moment $J_{2}^{\odot}$, if modeled at the level of accuracy reached by the recent INPOP17a ephemerides, i.e. $\sigma_{J_{2}^{\odot}}=1\times 10^{-9}$, would be about 10 times smaller than the relativistic signals of interest. Furthermore, BepiColombo should be able to constrain $J_{2}^{\odot}$ down to the $\sigma_{J_{2}\odot}\simeq 5\times 10^{-10}$ level. We successfully checked our analytical results for the Lense-Thirring range and range-rate shifts by comparing them to numerically produced ones by integrating the equations of motion; indeed, their differences ar as little as $\left|\Delta\rho_{\textrm{anal}}-\Delta\rho_{\textrm{num}}\right|\lesssim 5\times 10^{-5}\leavevmode\nobreak\ \textrm{m},\leavevmode\nobreak\ \left|\Delta\dot{\rho}_{\textrm{anal}}-\Delta\dot{\rho}_{\textrm{num}}\right|\lesssim 1\times 10^{-4}\leavevmode\nobreak\ \textrm{cm\leavevmode\nobreak\ s}^{-1}$. It turned out that the gravitomagnetic time delay due to the propagation of the electromagnetic waves is negligible since it leads to a range shift of the order of $\lesssim 10^{-5}\leavevmode\nobreak\ \textrm{m}$. The gravitomagnetic field of the Sun has been always neglected so far in all the existing studies dedicated to the scientific return of BepiColombo. Our results show that it is time to explicitly account for the Lense-Thirring effect in future analyses in order to investigate its actual detectability.

We also looked at the Lense-Thirring corrections $\Delta x_{\textrm{LT}},\leavevmode\nobreak\ \Delta y_{\textrm{LT}},\leavevmode\nobreak\ \Delta z_{\textrm{LT}},\leavevmode\nobreak\ \Delta\dot{x}_{\textrm{LT}},\leavevmode\nobreak\ \Delta\dot{y}_{\textrm{LT}},\leavevmode\nobreak\ \Delta\dot{z}_{\textrm{LT}}$ to the state vectors of Mercury and Earth as a preliminary insight for future, accurate investigations of their impact on long-term integrations of the solar system dynamics over the past and future $\sim\textrm{Gyr}$ scale. We found that, after just two years, the position $\mathbf{r}$ and the velocity $\mathbf{v}$ of Mercury and Earth are changed by $10\leavevmode\nobreak\ \textrm{m},\leavevmode\nobreak\ 1.5\leavevmode\nobreak\ \textrm{m}$ and $10^{-3}\leavevmode\nobreak\ \textrm{cm\leavevmode\nobreak\ s}^{-1},\leavevmode\nobreak\ 10^{-5}\leavevmode\nobreak\ \textrm{cm\leavevmode\nobreak\ s}^{-1}$, respectively. In light of the existing studies, all neglecting the general relativistic gravitomagnetic field of the Sun, such shifts may not be negligible over the eons; suffice it to say that it has been demonstrated in the literature that an error as little as $1\leavevmode\nobreak\ \textrm{cm}$ grows to about $1\leavevmode\nobreak\ \textrm{au}$ after about $100\leavevmode\nobreak\ \textrm{Myr}$.

Acknowledgements

I am grateful to an anonymous referee for his insightful comments and critical remarks