On the mean anomaly and the mean longitude in tests of post-Newtonian gravity

TL;DR

This paper clarifies the definitions of orbital anomalies and longitudes in post-Newtonian gravity tests, calculates their perturbations due to relativistic and classical effects, and discusses their application in astrophysics and satellite experiments.

Contribution

It provides a detailed analysis of orbital element perturbations in post-Newtonian gravity, correcting misconceptions and proposing new observables for testing general relativity.

Findings

No net Lense-Thirring perturbation on semimajor axis and mean motion.

Quadratic effects of non-gravitational forces can be separated from relativistic signals.

Orbital anomalies and longitudes can enhance tests of gravity with satellites.

Abstract

The distinction between the mean anomaly $\mathcal{M}(t)$ and the mean anomaly at epoch $\eta$, and the mean longitude $l(t)$ and the mean longitude at epoch $\epsilon$ is clarified in the context of a possible use of such orbital elements in post-Keplerian tests of gravity, both Newtonian and post-Newtonian. In particular, the perturbations induced on $\mathcal{M}(t),\,\eta,\,l(t),\,\epsilon$ by the post-Newtonian Schwarzschild and Lense-Thirring fields, and the classical accelerations due to the atmospheric drag and the oblateness $J_2$ of the central body are calculated for an arbitrary orbital configuration of the test particle and a general orientation of the primary's spin axis $\boldsymbol{\hat{S}}$. They provide us with further observables which could be fruitfully used, e.g., in better characterizing astrophysical binary systems and in more accurate satellite-based tests around major bodies of the Solar System. Some erroneous and misleading claims by Ciufolini and Pavlis appeared in the literature are confuted. In particular, it is shown that there are no net perturbations of the Lense-Thirring acceleration on either the semimajor axis $a$ and the mean motion $n_\mathrm{b}$. Furthermore, the quadratic signatures on $\mathcal{M}(t)$ and $l(t)$ due to certain disturbing non-gravitational accelerations like the atmospheric drag can be effectively disentangled from the post-Newtonian linear trends of interest provided that a sufficiently extended temporal interval for the data analysis is assumed. A possible use of $\eta$ along with the longitudes of the ascending node $\Omega$ in tests of general relativity with the existing LAGEOS and LAGEOS II satellites is suggested.

Tables

Table 1. Table 1: Orbital and physical configuration of a fictitious terrestrial geodetic satellite. Since it is ρ LARES = 5.96 × 10 − 16 ​ kg m − 3 subscript 𝜌 LARES 5.96 superscript 10 16 superscript kg m 3 \rho_{\textrm{LARES}}=5.96\times 10^{-16}\,\textrm{kg\,m}^{-3} Pardini et al . ( 2017 ) , and ρ LAGEOS = 6.579 × 10 − 18 ​ kg m − 3 subscript 𝜌 LAGEOS 6.579 superscript 10 18 superscript kg m 3 \rho_{\textrm{LAGEOS}}=6.579\times 10^{-18}\,\textrm{kg\,m}^{-3} Lucchesi et al . ( 2015 ) , we, first, used them in ρ ​ ( h ) = ρ 0 ​ exp ⁡ [ − ( h − h 0 ) ​ Λ − 1 ] 𝜌 ℎ subscript 𝜌 0 ℎ subscript ℎ 0 superscript Λ 1 \rho(h)=\rho_{0}\exp\left[-\left(h-h_{0}\right)\Lambda^{-1}\right] , where ρ 0 subscript 𝜌 0 \rho_{0} and h 0 subscript ℎ 0 h_{0} are, in general, referred to some reference height, to determine Λ Λ \Lambda in the case h 0 = h LARES , h = h LAGEOS formulae-sequence subscript ℎ 0 subscript ℎ LARES ℎ subscript ℎ LAGEOS h_{0}=h_{\textrm{LARES}},\,h=h_{\textrm{LAGEOS}} . Then, we used the so obtained characteristic length Λ LR/L = 999.51 ​ km subscript Λ LR/L 999.51 km \Lambda_{\textrm{LR/L}}=999.51\,\mathrm{km} , valid in the range h LARES = 1 , 442.06 ​ km < h < h LAGEOS = 5 , 891.96 ​ km formulae-sequence formulae-sequence subscript ℎ LARES 1 442.06 km ℎ subscript ℎ LAGEOS 5 891.96 km h_{\textrm{LARES}}=1,442.06\,\textrm{km}<h<h_{\textrm{LAGEOS}}=5,891.96\,\textrm{km} , to calculate ρ max subscript 𝜌 max \rho_{\textrm{max}} for our orbital geometry. Instead, the value ρ min subscript 𝜌 min \rho_{\textrm{min}} is just a guess which may be even conservative. The values of the satellite’s physical parameters were taken from Pardini et al . ( 2017 ) ( m , Σ , C D 𝑚 Σ subscript 𝐶 D m,\leavevmode\nobreak\ \Sigma,\leavevmode\nobreak\ C_{\textrm{D}} ).

Orbital and physical parameter	Numerical value	Units
Mass (LARES) $m$	$386.8$	kg
Area-to-mass ratio $\Sigma$ (LARES)	$2.69\times {10}^{-4}$	${\text{m}}^2{\text{kg}}^{-1}$
Neutral drag coefficient $C_{\text{D}}$ (LARES)	$3.5$	-
Semimajor axis $a$	$12,500$	km
Orbital period $P_{\text{b}}$	$3.86$	hr
Orbital eccentricity $e$	$0.36$	-
Perigee height $h_{\text{min}}$	$1,621.86$	km
Apogee height $h_{\text{max}}$	$10,621.9$	km
Orbital inclination $I$	$63.43$	deg
Argument of perigee $\omega$	$0$	deg
Period of the node $P_{\Omega }$	$-1.76$	yr
Period of the perigee $P_{\omega }$	$-2,903.62$	yr
Neutral atmospheric density at perigee ${\rho }_{\text{max}}$	$4.71\times {10}^{-16}$	kg m-3
Neutral atmospheric density at apogee ${\rho }_{\text{min}}$	$1\times {10}^{-20}$	kg m-3
Characteristic atmospheric length scale $\Lambda$	$836.34$	km

Equations

$$\mathcal{M}(t)\doteq \eta +n_{\mathrm{b}}\left(t-t_0\right),$$

$$n_{\mathrm{b}}=\sqrt{\frac{\mu }{a^3}}$$

$$n_{\mathrm{b}}\to n_{\mathrm{b}}^{\text{pert}}=n_{\mathrm{b}}+\Delta n_{\mathrm{b}}(t).$$

$${\mathcal{M}}^{\text{pert}}(t)$$

$$=\eta +\Delta \eta (t)+$$

$$+n_{\mathrm{b}}\left(t-t_0\right)+{\int }_{t_0}^t\Delta n_{\mathrm{b}}\left(t^{{}^{\prime }}\right)\mathrm{d}{t}^{{}^{\prime }}.$$

$$\Delta \mathcal{M}(t)={\mathcal{M}}^{\text{pert}}(t)-\mathcal{M}(t)=\Delta \eta (t)+\Phi (t),$$

$$\Phi (t)\doteq {\int }_{t_0}^t\Delta n_{\mathrm{b}}\left(t^{{}^{\prime }}\right)\mathrm{d}{t}^{{}^{\prime }}.$$

$$\Delta n_{\mathrm{b}}(t)$$

$$=-\frac{3}{2}\frac{n_{\mathrm{b}}}{a}{\int }_{f_0}^f\frac{\text{d}a}{\text{d}t}\frac{\text{d}t}{\text{d}{f}^{{}^{\prime }}}\mathrm{d}{f}^{{}^{\prime }}=$$

$$=-\frac{3}{2}\frac{n_{\mathrm{b}}}{a}\Delta a\left(f_0,f\right),$$

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

$$\frac{\text{d}a}{\text{d}t}$$

$$\frac{\text{d}\eta }{\text{d}t}$$

$$-\frac{\left(1-e^2\right)}{n_{\mathrm{b}}ae}\left[-A_R\cos f+A_T\left(1+\frac{r}{p}\right)\sin f\right],$$

$$\frac{\text{d}t}{\text{d}f}\simeq \frac{r^2}{\sqrt{\mu p}}+\mathcal{O}\left(A\right).$$

$$⟨\stackrel{\cdot }{\Phi }(t)⟩=-\frac{3}{4\text{mbox}p}\frac{n_{\mathrm{b}}^2}{a}{\int }_{f_0}^{f_0+2\text{mbox}p}\Delta a\left(f_0,f^{{}^{\prime }}\right)\frac{\text{d}t}{\text{d}{f}^{{}^{\prime }}}\mathrm{d}{f}^{{}^{\prime }}$$

$$l(t)\doteq \varpi +\mathcal{M}(t),$$

$$\varpi \doteq \Omega +\omega$$

$$l^{\mathrm{pert}}(t)$$

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

$$\Delta l(t)=\Delta ϵ(t)+\Phi (t).$$

$$\frac{\text{d}ϵ}{\text{d}t}$$

$$+2\sqrt{1-e^2}{\sin }^2\left(\frac{I}{2}\right)\frac{\text{d}\Omega }{\text{d}t},$$

$$\frac{\text{d}\Omega }{\text{d}t}$$

$$\frac{\text{d}\varpi }{\text{d}t}$$

$$+2{\sin }^2\left(\frac{I}{2}\right)\frac{\text{d}\Omega }{\text{d}t}.$$

$$u\doteq \omega +f$$

$${\mathbf{A}}^{\text{1pN}}$$

$$-\left(1+3\zeta \right){\mathrm{v}}^2+\frac{3}{2}\zeta {\mathrm{v}}_r^2]\hat{\mathbf{r}}+$$

$$+\left(4-2\zeta \right){\mathrm{v}}_r\mathrm{v}\},$$

$${\mu }_{\text{tot}}=G\left(m_{\text{A}}+m_{\text{B}}\right)$$

$${\mathrm{v}}_r\doteq \mathrm{v}\cdot \hat{\mathbf{r}}$$

$$\zeta \doteq \frac{m_{\mathrm{A}}{m}_{\mathrm{B}}}{{\left(m_{\mathrm{A}}+m_{\mathrm{B}}\right)}^2},0\le \zeta \le \frac{1}{4}.$$

$$\Delta n_{\mathrm{b}}\left(f_0,f\right)$$

$$\times \{4\left[-7+3\zeta +e^2\left(-3+4\zeta \right)\right]+$$

$$+e[e\zeta \cos 2f+4\left(-5+4\zeta \right)\cos f_0+$$

Full text

keywords General relativity and gravitation; Experimental studies of gravity; Experimental tests of gravitational theories; Satellite orbits

1 Introduction

In regard to possible tests of post-Newtonian (pN) features of general relativity111For a recent overview of the current status and challenges of the Einsteinian theory of gravitation, see, e.g., Debono and Smoot Debono & Smoot (2016), and references therein. and of alternative models of gravity with, e.g., Earth’s artificial satellites, Solar system’s planets and other astrophysical binaries, there is a considerable confusion in the literature about the possible use of the mean anomaly $\mathcal{M}\left(t\right)$ as potential observable in addition to the widely inspected argument of pericentre $\omega$ and, to a lesser extent, longitude of the ascending node $\Omega$. Indeed, it is as if some researchers, including the present author, who have tried to compute perturbatively the mean rate of change of the mean anomaly in excess with respect to the Keplerian case due to some pN accelerations were either unaware of the fact that what they, actually, calculated was the secular precession of the mean anomaly at the epoch $\eta$, or they systematically neglected a potentially non-negligible contribution to the overall change of the mean anomaly induced indirectly by the semimajor axis $a$ through the mean motion $n_{\rm b}$. Such a confusion has produced so far some misunderstanding which led, e.g., to unfounded criticisms about alleged proposals of using the mean anomaly, especially in the case of man-made spacecraft orbiting the Earth, or even uncorrect evaluations of the total pN effects sought. An example that sums up well the aforementioned confusion and misunderstanding, even in the peer-reviewed literature, is the following one. Ciufolini and Pavlis Ciufolini & Pavlis (2005) wrote “[…] one of the most profound mistakes and misunderstandings of Iorio (2005) is the proposed use of the mean anomaly of a satellite to measure the Lense-Thirring effect […] This is simply a nonsense statement: let us, for example, consider a satellite at the LAGEOS altitude, the Lense-Thirring effect on its mean longitude is of the order of 2 m/y, however, the mean longitude change is about $1.8\times 10^{11}$ m/y. Thus, from Kepler’s law, the Lense-Thirring effect corresponds to a change of the LAGEOS semi-major of less than $0.009$ cm! Since, even a high altitude satellite such as LAGEOS showed a semimajor axis change of the order of 1 mm/day, due to atmospheric drag and to the Yarkoski-Rubincam effect (because of atmospheric drag, the change of semimajor axis and mean motion is obviously much larger for lower altitude satellites), and since the present day precision of satellite laser ranging is, even in the case of the best SLR stations, of several millimeters, it is a clear nonsense to propose a test of the Lense-Thirring effect based on using the mean anomaly of any satellite, mean anomaly largely affected by non-conservative forces.” It is difficult to understand what is the target of the arrows by Ciufolini and Pavlis Ciufolini & Pavlis (2005) since the mean anomaly is not even mentioned in the published version of the criticized paper by the present author, not to mention any explicitly detailed proposal to use it. Be that as it may, in the following, we will show that, actually, using the mean anomaly, or the mean longitude $l\left(t\right)$, in pN tests with artificial Earth’s satellites may be feasible, provided that certain non-gravitational perturbations are compensated by some active drag-free mechanism. However, even in case of passive, geodetic satellites, we will show that, under certain conditions, it is possible to separate the relativistic linear trends of interest from the unwanted parabolic signatures of non-conservative origin. Furthermore, the arguments provided by Ciufolini and Pavlis Ciufolini & Pavlis (2005) about the Lense-Thirring effect and the mean longitude are erroneous. Finally, the use of the mean anomaly at epoch or of the mean longitude at epoch $\epsilon$ is, in principle, possible even with passive, geodetic spacecraft like those of the LAGEOS family because they are, by construction, free from the aforementioned potential drawbacks exhibited by the mean anomaly and the mean longitude themselves, which was completely ignored or unrecognized by Ciufolini and Pavlis Ciufolini & Pavlis (2005).

The paper is organized as follows. In Section 2, we will review the basics of the mean anomaly, the mean anomaly at epoch, the mean longitude, and the mean longitude at epoch along with the calculation of their perturbations with respect to the purely Keplerian case in presence of a generic disturbing post-Keplerian (pK) acceleration. Section 3 is devoted to the calculation of the effects of some well-known pN accelerations (Schwarzschild and Lense-Thirring), while the impact of the atmospheric drag and the oblateness of the primary are treated in Section 4. The potential of a possible use of the mean anomaly at epoch in the ongoing tests with the satellites LAGEOS and LAGEOS II is discussed in Section 5. Section 6 summarizes our findings, and offers our conclusions.

2 The mean anomaly and the mean longitude

2.1 The mean anomaly

In the restricted two-body problem, the mean anomaly $\mathcal{M}(t)$ is one of the three time-dependent fast angular variables which, in celestial mechanics, can be used to characterize the instantaneous position of a test particle along its Keplerian ellipse, being the eccentric anomaly $E$ and the true anomaly $f$ the other two anomalies. In the unperturbed Keplerian case, the mean anomaly is defined as

[TABLE]

where222The symbol $\eta$ is used for the mean anomaly at epoch by Milani et al. Milani, Nobili & Farinella (1987). In the notation by Brumberg Brumberg (1991), the mean anomaly is $l$, while the mean anomaly at epoch is $l_{0}$. Kopeikin et al. Kopeikin, Efroimsky & Kaplan (2011) denote $\eta$ as $\mathcal{M}_{0}$, while Bertotti et al. Bertotti, Farinella & Vokrouhlický (2003) adopt $\epsilon^{{}^{\prime}}$. $\eta$ is the mean anomaly at the reference epoch $t_{0}$, and

[TABLE]

is the Keplerian mean motion. In Equation (2), $\mu\doteq GM$ is the gravitational parameter of the primary having mass $M$, while $G$ is the Newtonian constant of gravitation; in the following, we will assume $\mu=\textrm{const}$. The mean anomaly at epoch $\eta$ is one of the six Keplerian orbital elements parameterizing the orbit of a test particle in space. In the unperturbed case, $\mathcal{M}(t)$ is a linear function of time $t$ because both $a$ and $\eta$ are constants of motion. If a relatively small perturbing pK acceleration $\boldsymbol{A}$ is present, both $a$ and $\eta$ are, in general, affected by it, becoming time-dependent. As a result, also the mean motion is, in general, modified so that

[TABLE]

Thus, the perturbed mean anomaly is the sum of the now time-dependent mean anomaly at epoch $\eta(t)$ and a function of time $\varrho(t)$ whose derivative is equal to the (perturbed) mean motion, i.e.,

[TABLE]

The resulting change $\Delta\mathcal{M}$ of the mean anomaly with respect to the unperturbed case is, thus,

[TABLE]

where we defined

[TABLE]

as a function whose derivative yields the perturbation of the mean motion. In Equation (6), the instantaneous shift of the mean motion due to the time-varying semimajor axis333It should be recalled that we kept $\mu$ constant. $a$ is

[TABLE]

so that

[TABLE]

The shifts $\Delta a(t)$ and $\Delta\eta(t)$ can be perturbatively calculated by evaluating the right-hand-sides of the Gauss equations for their rates of change Bertotti, Farinella & Vokrouhlický (2003)

[TABLE]

onto the unperturbed Keplerian ellipse. In Equations (9) to (2.1), $e$ is the eccentricity, $p\doteq a\left(1-e^{2}\right)$ is the semilatus rectum, $r=p/\left(1+e\cos f\right)$ is the (unperturbed) distance of the test particle from the primary, and $A_{R}$, $A_{T}$ are the projections of the perturbing pK acceleration $\boldsymbol{A}$ onto the radial and transverse directions, respectively. The derivative of $t$ with respect to $f$ entering Equations (7) to (8) is, up to terms of the first order in the perturbing acceleration $A$,

[TABLE]

Depending on the disturbing acceleration, $\Phi(t)$ is linear in time if the average over an orbital period $P_{\rm b}$ of its rate of change

[TABLE]

is constant. Otherwise, it may exhibit a more complex temporal pattern, as when the semimajor axis $a$ undergoes a secular change due to, e.g., some non-gravitational perturbing accelerations as in artificial satellites’ dynamics. In general, the calculation of $\Phi(t)$ is rather cumbersome since it involves two integrations. Moreover, it depends on $f_{0}$.

From such considerations it follows that, at first sight, using the mean anomaly $\mathcal{M}(t)$ may not be a wise choice because of the disturbances introduced by $\Phi(t)$, especially in non-trivial scenarios in which several perturbing accelerations of different nature act simultaneously on the test particle inducing non-vanishing long-term effects on the semimajor axis $a$. Actually, we will show that it may not be the case in practical satellite data reductions if certain conditions are fulfilled. On the contrary, the mean anomaly at the epoch $\eta$, which is one of the six osculating Keplerian orbital elements in the perturbed restricted two-body problem, is not affected by such drawbacks. As such, it can be safely used, at least in principle, as an additional piece of information to improve some tests of pN gravity on the same foot of $\omega$ and $\Omega$. This fact seems to have gone unnoticed so far in the literature, as in the case of Ciufolini and Pavlis Ciufolini & Pavlis (2005).

2.2 The mean longitude

Similar considerations hold for the mean longitude $l$ defined as

[TABLE]

where

[TABLE]

is the longitude of pericentre. If a disturbing acceleration $\boldsymbol{A}$ is present, it can be expressed in terms of the mean longitude at epoch444It is more suited than $\eta$ at low orbital inclinations Bertotti, Farinella & Vokrouhlický (2003). $\epsilon$ Milani, Nobili & Farinella (1987); Soffel (1989); Brumberg (1991); Bertotti, Farinella & Vokrouhlický (2003) as

[TABLE]

so that its shift is

[TABLE]

The shift $\Delta\epsilon(t)$ of the mean longitude at epoch can be worked out by means of the Gauss equation for its variation Bertotti, Farinella & Vokrouhlický (2003)

[TABLE]

where Bertotti, Farinella & Vokrouhlický (2003)

[TABLE]

In Equation (18), $A_{N}$ is the projection of the perturbing acceleration $\boldsymbol{A}$ onto the out-of-plane direction, while

[TABLE]

is the argument of latitude.

3 The secular rates of change of $\eta(t),\,\epsilon(t),\,\Phi(t)$ for some pN accelerations

Here, we will preliminarily look at the effects due to the standard general relativistic pN accelerations induced by the static, gravitoelectric (Schwarzschild, Section 3.1) and stationary, gravitomagnetic (Lense-Thirring, Section 3.2) components of the spacetime of an isolated rotating body. We will not restrict to almost circular orbits; furthermore, we will allow the primary’s spin axis $\boldsymbol{\hat{S}}$, entering the Lense-Thirring acceleration, to assume any orientation in space.

3.1 The 1pN gravitoelectric Schwarzschild-like acceleration

To the first pN order (1pN), the relative acceleration for two pointlike bodies of masses $m_{\mathrm{A}},\,m_{\mathrm{B}}$ separated by a distance $r$ and moving with relative velocity $\mathbf{v}$ is Damour & Deruelle (1985); Soffel (1989)

[TABLE]

where $c$ is the speed of light in vacuum,

[TABLE]

is the total gravitational parameter of the binary system,

[TABLE]

is the the radial velocity of the relative orbital motion, and

[TABLE]

In Sections 3.1.1 to 3.1.2, we will work out the effect of Equation (21) on $\Phi$, and $\eta$ and $\epsilon$, respectively.

3.1.1 The shift $\Phi(t)$ due to the variation of the mean motion

By using Equation (21) in Equation (7) yields

[TABLE]

From it, the rate of change of $\Phi(t)$ averaged over one orbital period $P_{\rm b}$ can be straightforwardly worked out as

[TABLE]

where

[TABLE]

From the analytical expression of the right-hand-side of Equation (3.1.1), it turns out that the true anomaly $f$, and, thus, also the time $t$, appears only in trigonometric functions. This implies that, in this case, $\Phi(t)$ does not exhibit a polynomial temporal pattern, being, at most, linear in $t$ provided that Equation (26) is not vanishing. Note also the dependence of Equation (3.1.1) on $f_{0}$. We are not able to analytically calculate Equation (26) unless a power expansion in $e$ of Equation (3.1.1) is made. Nonetheless, it is possible to perform a numerical integration of Equation (26) for given values of the physical and orbital parameters entering it without any restriction on $e$. We successfully tested it for a fictitious cannonball geodetic satellite moving along an eccentric orbit, whose arbitrarily chosen physical and orbital parameters are displayed in Table 1,

by numerically integrating its equations of motion in rectangular Cartesian coordinates, and by numerically performing the integral of Equation (26) with Equation (3.1.1). Fig. 1 displays the plot of Equation (3.1.1), in milliarcseconds per year $\left(\mathrm{mas\,yr}^{-1}\right)$, for the orbital parameters of Table 1, and the numerically produced time series of $\Phi(t)$, in mas, over 1 yr for the same orbital configuration; the agreement between the slope of $\Phi(t)$ and the area under the curve of Equation (3.1.1) is remarkable. From Fig. 1, it can be noted that, as expected, the 1pN Schwarzschild-like acceleration induces a secular variation on $\Phi(t)$ which has to be added to those affecting $\eta$ and $\epsilon$ displayed in Section 3.1.2.

The opportunity offered by the exact expression of Equation (3.1.1) to calculate $\left\langle\overset{\cdot}{\Phi}(t)\right\rangle$ as per Equation (26) is important also in astronomical and astrophysical scenarios, like the almost circular orbital motions of the major bodies of our solar system and the much more eccentric ones of various types of binary systems (extrasolar planets, binary stars, binary pulsars hosting at least one emitting neutron star, stellar systems revolving around supermassive galactic black holes, etc.), in which secular variations of the semimajor axis $a$-or even of the masses involved-are absent or negligible with respect to either the duration of the typical data analyses or to the observational accuracy. Indeed, in all such cases, the perturbed evolution of the mean anomaly can, in principle, be monitored as well, and $\left\langle\overset{\cdot}{\Phi}(t)\right\rangle$ may represent an important contribution to the overall long-term rate of change of $\mathcal{M}(t)$. Suffice it to say that, in the case of Mercury and the Sun, it is

[TABLE]

3.1.2 The mean anomaly at epoch $\eta$ and the mean longitude at epoch $\epsilon$

The Gauss equations for the variation of $\eta$ and $\epsilon$ (Equation (2.1) and Equation (2.2)) allow to straightforwardly work out their secular rates of change which turn out to be

[TABLE]

They were confirmed by a numerical integration of the equations of motion in the case of the satellite’s orbital configuration of Table 1 which returned linear times series whose slopes agree with Equations (29) to (30). In the case of Mercury and the Sun, Equations (29) to (30) yield

[TABLE]

3.2 The 1pN gravitomagnetic Lense-Thirring acceleration

In the case of the 1pN gravitomagnetic Lense-Thirring acceleration Soffel (1989) induced by the spin dipole moment of the central mass, i.e. its proper angular momentum $\boldsymbol{S}$, on a test particle orbiting it with velocity $\mathbf{v}$

[TABLE]

it turns out that

[TABLE]

for an arbitrary orientation of the body’s spin axis $\boldsymbol{\hat{S}}$ in space. Thus, it is

[TABLE]

It implies that the claims by Ciufolini and Pavlis Ciufolini & Pavlis (2005) about an alleged non-vanishing perturbing effect of the gravitomagnetic field of the Earth on both the semimajor axis $a$ and the mean motion $n_{\rm b}$ of a satellite are, in fact, erroneous for any spacecraft.

Moreover, it is also

[TABLE]

for any $\boldsymbol{\hat{S}}$ as well. In Equation (37),

[TABLE]

is the unit vector directed along the orbital angular momentum along the out-of-plane direction, while

[TABLE]

is the unit vector directed transversely to the line of the nodes in the orbital plane. In the case of an Earth’s satellite, by assuming, as usual, an equatorial coordinate system with its reference $z$ axis directed along $\boldsymbol{\hat{S}}$, Equation (37) reduces to

[TABLE]

Equation (35) and Equation (40) show that the claim by Ciufolini and Pavlis Ciufolini & Pavlis (2005) “[…] let us, for example, consider a satellite at the LAGEOS altitude, the Lense-Thirring effect on its mean longitude is of the order of 2 m/y, […]” is wrong. Indeed, the gravitomagnetic linear shift corresponding to Equation (40) amounts to $3.68\,\mathrm{m\,yr}^{-1}$ for LAGEOS; it is an enormous discrepancy with respect to the statement by Ciufolini and Pavlis Ciufolini & Pavlis (2005) since the present-day accuracy in reconstructing the orbits of the laser-ranged satellites of the LAGEOS type is notoriously at the $\simeq 1-0.5\,\mathrm{cm}$ level.

4 The secular rates of change of $\eta(t),\,\epsilon(t),\,\Phi(t)$ for some Newtonian perturbing accelerations

Here, we will deal with the impact of the oblateness of the primary (Section 4.1), whose spin axis $\boldsymbol{\hat{S}}$ is assumed arbitrarily oriented in space, and of the atmospheric drag (Section 4.2). The small eccentricity approximation for the satellite’s orbit will not be adopted. Such classical accelerations represent two of the most important sources of systematic errors in accurate tests of pN gravity with artificial satellites. On the other hand, they can be considered interesting in themselves if one is interested in better characterizing the shape and the inner mass distribution of the primary like, e.g., a star, at hand, and the properties of the atmosphere of the orbited planet.

4.1 The quadrupole mass moment $J_{2}$

To the Newtonian level, the external potential of an oblate body at the outside position $\boldsymbol{r}$ is

[TABLE]

where $J_{2}$ is the first even zonal harmonic coefficient of the multipolar expansion of its classical gravitational potential,

[TABLE]

is the cosine of the angle between the primary’s spin axis and the particle’s position, and

[TABLE]

is the Legendre polynomial of degree 2. The Newtonian acceleration due to $J_{2}$ experienced by a test particle orbiting the distorted axisymmetric primary is

[TABLE]

In Sections 4.1.1 to 4.1.2, we will work out its impact on $\Phi(t)$, and $\eta$ and $\epsilon$, respectively.

4.1.1 The shift $\Phi(t)$ due to the variation of the mean motion

It turns out that

[TABLE]

so that $\Phi(t)$, which depends on $f_{0}$, is linear in time. It is not possible to explicitly display the analytical expression which we obtained for $\left(1/P_{\rm b}\right)\,\mathrm{d}\Phi/\mathrm{d}f$ in the case of an arbitrary orientation of $\boldsymbol{\hat{S}}$ in space because of its cumbersomeness. However, it can be fruitfully used with, e.g., any astronomical binary systems since, in general, their spin axes are not aligned with the line of sight which, usually, is assumed as reference $z$ axis of the coordinate systems adopted. In regard to an Earth’s satellite, whose motion is customarily studied in an equatorial coordinate system whose reference $z$ axis is aligned with $\boldsymbol{\hat{S}}$, we have

[TABLE]

with

[TABLE]

The numerical value of the area under the plot of Equation (46), depicted in the upper panel of Fig. 2, is confirmed by the time series for $\Phi(t)$ produced by numerically integrating the equations of motion of the fictitious satellite of Table 1, and displayed in the lower panel of Fig. 2.

4.1.2 The mean anomaly at epoch $\eta$ and the mean longitude at epoch $\epsilon$

The Gauss equations for the variations of $\eta$ and $\epsilon$ (Equation (2.1) and Equation (2.2)) allow to straightforwardly obtain

[TABLE]

where $\boldsymbol{\hat{l}}=\left\{\cos\Omega,\,\sin\Omega,\,0\right\}$ is the unit vector directed along the line of the nodes such that $\boldsymbol{\hat{l}}\boldsymbol{\boldsymbol{\times}}\boldsymbol{\hat{m}}=\boldsymbol{\hat{h}}$. Also Equations (48) to (4.1.2) can be used with any astronomical binary system in view of their generality. In the case of a coordinate system with its reference $z$ axis aligned with the body’s spin axis, as in the case of an Earth’s satellite referred to an equatorial coordinate system, Equations (48) to (4.1.2) reduce to

[TABLE]

4.2 The atmospheric drag

The atmospheric drag induces, among other things, a secular decrease of the semimajor axis $a$ which, in turn, has an impact on $n_{\rm b}(t)$ and $\Phi(t)$.

For a cannonball geodetic satellite, the drag acceleration can be expressed as

[TABLE]

In Equation (52), $C_{\textrm{D}},\leavevmode\nobreak\ \Sigma,\leavevmode\nobreak\ \rho,\leavevmode\nobreak\ \boldsymbol{V}$ are the dimensionless drag coefficient of the satellite, its area-to-mass ratio, the atmospheric density at its height, and its velocity with respect to the atmosphere, respectively. In the following, we will assume that the atmosphere co-rotates with the Earth. Thus, $\boldsymbol{V}$ is

[TABLE]

where $\boldsymbol{\Psi}$ is the Earth’s angular velocity. We will model the atmospheric density as

[TABLE]

where $\rho_{0}$ refers to some reference distance $r_{0}$, while $\Lambda$ is the characteristic scale length. By assuming

[TABLE]

$\Lambda$ can be determined as

[TABLE]

where

[TABLE]

are the values of the atmospheric density at the apogee and perigee heights, respectively. Table 1 shows the neutral atmospheric density at the perigee height chosen as inferred from existing data on LAGEOS and LARES. On the other hand, the values reported for the apogee are purely speculative and should be regarded as subjected to huge uncertainties. Actually, even the density at a given height may not be regarded as truly constant because of a variety of geophysical phenomena characterized by quite different time scales. Anyway, in order to have an order-of-magnitude evaluation of the perturbing action of Equation (52) on the motion of the fictitious satellite of Table 1, we will make our calculation by keeping $\rho_{0}$ fixed during one orbital period $P_{\textrm{b}}$. An exact analytical calculation without recurring to any approximation in both $e$ and $\nu\doteq\Psi/n_{\rm b}$ is difficult.

In Sections 4.2.1 to 4.2.2, we will calculate the impact of Equation (52) on $\Phi(t)$, and $\eta$ and $\epsilon$, respectively.

4.2.1 The shift $\Phi(t)$ due to the variation of the mean motion

Let us, now, start to look at $\Delta n_{\rm b}(t)$ by means of Equation (7). We will show that it is linear in time because $\left\langle{\overset{.}{\Delta n_{\rm b}}}\right\rangle\neq 0$. The analytical expression of $1/P_{\rm b}\,\mathrm{d}\Delta n_{\rm b}/\mathrm{d}f$ is

[TABLE]

where

[TABLE]

Since it is not possible to analytically integrate Equation (4.2.1) with Equation (4.2.1) in the most general case without recurring to approximations in $e$ and $\nu$, we will plot it as a function of $f$ over a full orbital cycle and integrate it numerically for the physical and orbital parameters of Table 1. The upper panel of Fig. 3 depicts Equation (4.2.1), while the lower panel displays the time series for $\Delta n_{\rm b}(t)$ calculated from a numerical integration of the satellite’s equations of motion in rectangular Cartesian coordinates over 1 yr.

The fact that $\left\langle\overset{.}{\Delta n_{\rm b}}\right\rangle\neq 0$ implies that $\Delta n_{\rm b}$ is linear555Strictly speaking, it is, in general, true only for fast satellites orbiting in much less than a day, so that the term proportional to $\nu^{2}$ in Equation (4.2.1), which contains $\omega$, can be neglected. However, in the particular case of the fictitious satellite of Table 1, $\omega$ stays essentially constant because of the frozen perigee configuration. in time and, thus, $\Phi(t)$ is quadratic. It is explicitly shown in Fig. 4 by the time series calculated for Equation (8) from the same integration of the satellite’s equations of motion.

It is an important feature because it allows to accurately separate the unwanted parabolic signature due to the atmospheric drag from the relativistic trend of interest affecting the time series of $\mathcal{M}(t)$ or $l(t)$, provided that a sufficiently long time span is chosen for the data analysis. The same holds, in principle, also for any other perturbing acceleration of non-gravitational origin inducing a secular trend in the satellite’s semimajor axis like, e.g., the Yarkovsky-Rubincam thermal effect. We numerically confirmed that by integrating the equations of motion of the fictitious satellite of Table 1 including the 1pN Schwarzschild-like and the atmospheric drag accelerations, and fitting a linear plus quadratic model to the resulting time series of $\Phi(t)$ over, say, 5 yr for a given value of $f_{0}$. As a result, we were able to accurately recover the slope of the relativistic secular signal. We successfully repeated it for different values of $f_{0}$ as well. It turns out that the longer the data span is, the more accurate the recovery of the linear signal. This suggests that, actually, also the mean anomaly $\mathcal{M}(t)$ and the mean longitude $l(t)$ may be fruitfully used in tests of pN gravity in the field of the Earth even with passive artificial satellites, contrary to the claims by Ciufolini and Pavlis Ciufolini & Pavlis (2005). The dependence of $\Phi(t)$ on $f_{0}$ may even represent an advantage to enhance the signal-to-noise ratio since, in principle, one can choose $f_{0}$ in order to maximize the relativistic rate for $\left\langle\overset{\cdot}{\Phi}\right\rangle$ to be added to the further contribution due to $\left\langle\overset{\cdot}{\eta}\right\rangle,\,\left\langle\overset{\cdot}{\epsilon}\right\rangle$.

4.2.2 The mean anomaly at epoch $\eta$ and the mean longitude at epoch $\epsilon$

About the secular rates of $\eta$ and $\epsilon$, the Gauss equations for their variations allow to obtain

[TABLE]

Since it is not possible to analytically integrate Equation (4.2.2) and Equation (4.2.2) in an exact form, we, first, plot them as functions of $f$ over a full orbital cycle in Fig. 5 for the orbital configuration of Table 1, and, then, numerically calculate the areas under their curves in order to obtain $\left\langle\overset{\cdot}{\eta}\right\rangle,\,\left\langle\overset{\cdot}{\epsilon}\right\rangle$.

Also in this case, a numerical integration of the satellite’s equations of motion turns out to confirm such results.

5 Some possible uses with the LAGEOS and LAGEOS II satellites

As an illustrative example, here we will look at the possibility of using the nodes $\Omega$ and the mean anomalies at epoch $\eta$ of, say, the existing satellites LAGEOS and LAGEOS II in order to propose an accurate test of the 1pN Lense-Thirring effect exploiting their multidecadal data records.

The availability of $\eta$ in addition to $\Omega$ may be particularly important in view of the fact that the competing classical secular precessions due to the even zonals of low degree, which have just the same time signature of the gravitomagnetic ones of interest, are nominally several orders of magnitude larger than them; thus, the signal-to-noise ratio must be somehow enhanced. The present-day level of actual mismodeling in the geopotential coefficients, which should be considered as (much) worse than the mere formal, statistical sigmas of the various global gravity field solutions666They are freely available on the Internet at the webpage of the International Centre for Global Earth Models (ICGEM), currently located at http://icgem.gfz-potsdam.de/tom$\_$longtime. releasing the experimentally estimated values of the geopotential’s parameters, does not yet allow to use the residuals of a single orbital element separately. To circumvent such an issue, some strategies involving the simultaneous use of more than one orbital element have been devised so far over the years: for a general overview, see, e.g., Renzetti Renzetti (2013), and references therein. To the benefit of the reader, we review here the linear combination approach, which is an extension of the one proposed by Ciufolini Ciufolini (1996) to test the gravitomagnetic field of the Earth with artificial satellites of the LAGEOS family. In turn, it is a generalization of the strategy put forth, for the first time, by I.I. Shapiro Shapiro (1990) who, at that time, wanted to separate the Sun-induced 1pN gravitoelectric perihelion precession from that due to the solar quadrupole mass moment $J_{2}$ by using other planets or highly eccentric asteroids.

By looking at $N$ orbital elements777At least one of them must be affected also by the 1pN effect one is looking for. The $N$ orbital elements $\kappa^{(i)}$ may be different from one another belonging to the same satellite, or some of them may be identical belonging to different spacecraft (e.g., the nodes of two different vehicles). $\kappa^{(i)},\,i=1,2,\ldots N$ experiencing, among other things, classical secular precessions due to the even zonals of the geopotential, the following $N$ linear combinations can be written down

[TABLE]

They involve the 1pN averaged precessions $\left\langle\overset{\cdot}{\kappa}\right\rangle_{\mathrm{1pN}}^{(i)}$ as predicted by General Relativity and scaled by a multiplicative parameter888It is equal to 1 in the Einstein’s theory of gravitation, and 0 in the Newtonian one. In general, ${\mathchoice{\mbox{m}}{\mbox{m}}{\mbox{m}}{\mbox{m}}}_{\mathrm{1pN}}$ is not necessarily one of the parameters of the parameterized post-Newtonian (PPN) formalism, being possibly a combination of some of them. ${\mathchoice{\mbox{m}}{\mbox{m}}{\mbox{m}}{\mbox{m}}}_{\mathrm{1pN}}$, and the errors in the computed secular node precessions due to the uncertainties in the first $N-1$ even zonals $J_{2s},\leavevmode\nobreak\ s=1,2,\ldots N-1$, assumed as mismodeled through $\delta J_{2s},\leavevmode\nobreak\ s=1,2,\ldots N-1$. In the following and in Appendix A, we will use the shorthand

[TABLE]

for the partial derivative of the classical averaged precession $\left\langle\overset{\cdot}{\kappa}\right\rangle_{J_{\ell}}$ with respect to the generic even zonal $J_{\ell}$ of degree $\ell$. Then, the $N$ combinations of Equation (63) are posed equal to the experimental residuals $\delta\overset{\cdot}{\kappa}^{(i)},\leavevmode\nobreak\ i=1,2,\ldots N$ of each of the $N$ orbital elements considered getting

[TABLE]

It should be recalled that, in principle, the residuals $\delta\overset{\cdot}{\kappa}^{(i)}$ account for the purposely unmodelled 1pN effect, the mismodelling of the static and time-varying parts of the geopotential, and the non-gravitational forces. If we look at the 1pN scaling parameter ${\mathchoice{\mbox{m}}{\mbox{m}}{\mbox{m}}{\mbox{m}}}_{\mathrm{1pN}}$ and the mismodeling in the even zonals $\delta J_{2s},\leavevmode\nobreak\ s=1,2,\ldots N-1$ as unknowns, we can interpret Equation (65) as an inhomogenous linear system of $N$ algebraic equations in the $N$ unknowns

[TABLE]

whose coefficients are

[TABLE]

while the constant terms are the $N$ orbital residuals

[TABLE]

It turns out that, after some algebraic manipulations, the dimensionless 1pN scaling parameter can be expressed as

[TABLE]

In Equation (69), the combination of the $N$ orbital residuals

[TABLE]

is, by construction, independent of the first $N-1$ even zonals, being, instead, impacted by the other ones of degree $\ell>2(N-1)$ along with the non-gravitational perturbations and other possible orbital perturbations which cannot be reduced to the same formal expressions of the first $N-1$ even zonal rates. On the other hand,

[TABLE]

combines the $N$ 1pN orbital precessions as predicted by General Relativity. The dimensionless coefficients $c_{j},\ j=1,2,\ldots N-1$ in Equation (70)-Equation (71) depend only on some of the orbital parameters of the satellite(s) involved in such a way that, by construction, ${\mathcal{C}}_{\delta}=0$ if Equation (70) is calculated by posing

[TABLE]

for any of the first $N-1$ even zonals, independently of the value assumed for its uncertainty $\delta J_{\ell}$.

As far as the Lense-Thirring effect and the satellites LAGEOS and LAGEOS II are concerned, the linear combination of the four experimental residuals $\delta\Omega^{\mathrm{L}},\,\delta\Omega^{\mathrm{L\,II}},\,\delta\eta^{\mathrm{L}},\,\delta\eta^{\mathrm{L\,II}}$ of the satellites’s nodes and mean anomalies at epoch suitably designed to cancel out the secular precessions due to the first three even zonal harmonics $J_{2},\,J_{4},\,J_{6}$ of the geopotential is

[TABLE]

whose coefficients $c_{1},\leavevmode\nobreak\ c_{2},\leavevmode\nobreak\ c_{3}$ are purposely constructed with the results of Section A.1. They turn out to be

[TABLE]

where the common denominator is

[TABLE]

Their numerical values, computed with the satellites’ orbital elements inserted in Equations (A6) to (A.1), are

[TABLE]

Thus, the predicted combined Lense-Thirring signature is

[TABLE]

The combination of Equation (73) is mainly affected by the orbital precessions induced by the fourth even zonal harmonic $J_{8}$ of the geopotential. The resulting mismodeled combined signal can be evaluated by means of Equations (A.1) to (A.1) along with some measure of the uncertainty in $J_{8}$. If one were to rely upon on the formal sigmas of the latest global Earth’s gravity field models by the dedicated GRACE and GOCE missions, the resulting impact on Equation (81) would be much smaller than $1\%$. Indeed, from, e.g., the zero-tide model Tongji-Grace02s Chen et al. (2018), it is999The zonal harmonics $J_{\ell}$ of the geopotential are connected with its fully normalized Stokes coefficients ${\overline{C}}_{\ell,0}$ by the relation $J_{\ell}=-\sqrt{2\ell+1}\,{\overline{C}}_{\ell,0},\,\ell=2,\,3,\,4,\ldots$ ${\mathchoice{\mbox{s}}{\mbox{s}}{\mbox{s}}{\mbox{s}}}_{{\overline{{C}}_{8,0}}}=1.3\times 10^{-14}$. It implies a combined mismodeled precessions as little as $0.01\,\mathrm{mas\,yr}^{-1}$, corresponding to $0.01\%$ of the combined Lense-Thirring effect. If, instead, the difference $\Delta{\overline{C}}_{8,0}$ between the values of ${\overline{C}}_{8,0}$ from Tongji-Grace02s and the zero-tide model ITU$\_$GRACE16 Akyilmaz et al. (2016), whose formal errors are comparable, is adopted as a measure of the actual uncertainty in the even zonal of degree 8, the resulting mismodeled signal amounts to $2.1\,\mathrm{mas\,yr}^{-1}$ corresponding to a percent error in the Lense-Thirring combined signature of $1.8\%$.

In fact, an accurate investigation, both analytical and numerical, of the perturbations on $\eta$ induced by the main non-gravitational accelerations acting on the LAGEOS-type satellites like, e.g., the direct solar radiation pressure, the Earth’s albedo, the Earth’s direct infrared radiation pressure, the Earth’s Yarkovsky-Rubincam and Solar Yarkovsky-Schach thermal effects, possible anisotropic reflectivity, etc. Lucchesi (2001, 2002, 2003); Pardini et al. (2017); Visco & Lucchesi (2018); Lucchesi et al. (2019) is required to realistically assess the overall error budget of the promising combination of Equation (73). This is outside the scopes of the present paper.

6 Summary and overview

In presence of Newtonian, general relativistic 1pN or modified gravity-induced disturbing accelerations, the shifts $\Delta\mathcal{M}(t)$ and $\Delta l(t)$ of the mean anomaly $\mathcal{M}(t)$ and the mean longitude $l(t)$ with respect to their Keplerian linear trends are, in general, due to the perturbations $\Delta\eta(t)$ and $\Delta\epsilon(t)$ of the mean anomaly at epoch $\eta$ and mean longitude at epoch $\epsilon$, and the change $\Delta n_{\rm b}(t)$ in the mean motion $n_{\rm b}$ which, in some cases, can induce a quadratic shift $\Phi(t)$ in $\mathcal{M}(t)$ and $l(t)$ depending on the true anomaly at epoch $f_{0}$.

In the case of an Earth’s artificial satellite, the atmospheric drag affects $\Phi(t)$ quadratically; nonetheless, the non-Newtonian linear trends of interest may be effectively separated from such a potentially competing aliasing effect if a sufficiently long time span for the data analysis is adopted. Thus, also $\mathcal{M}(t)$ and $l(t)$ can, in principle, be employed in gravity tests even with passive geodetic satellites, not to mention the use of drag-free apparatuses. If, instead, $\eta$ and $\epsilon$ are adopted, such an issue is a-priori circumvented because they are not impacted by the possible change in the mean motion $n_{\rm b}$. Since $\eta$ and $\epsilon$ undergo secular precessions due to the even zonal harmonics $J_{\ell},\,\ell=2,\,4,\ldots$ of the geopotential, it is possible, in principle, to use them in combination with, say, the nodes $\Omega$ to reduce the impact of the mismodeled even zonals in experiments of fundamental physics with existing satellites. In an actual test, a detailed analysis of the perturbations affecting $\eta$ and $\epsilon$ by all the most relevant non-gravitational accelerations should be performed. There are no net Lense-Thirring rates of change of the semimajor axis $a$ and $n_{\rm b}$.

In astronomical binary systems, not affected by non-gravitational perturbations, using $\eta$ may provide a further valuable observable in addition to the usual periastron precession to put to the test general relativity and, say, modified models of gravity, or to better characterize the physical properties of the bodies like, e.g., their oblateness $J_{2}$ and their orbital configurations as well. Indeed, the 1pN effects on $\eta$ are often larger than the corresponding pericenter rates.

Appendix A Mean orbital precessions of $\Omega$ and $\eta$ due to the even zonal harmonics of the geopotential

Here, we analytically calculate the coefficients

[TABLE]

of the precessions

[TABLE]

of the node $\Omega$ and of the mean anomaly at epoch $\eta$ averaged over one full orbital period $P_{\rm b}$, induced by the first four even zonal harmonics $J_{\ell}$. To this aim, we use the standard Lagrange planetary equations Bertotti, Farinella & Vokrouhlický (2003)

[TABLE]

In them, the correction of degree $\ell$

[TABLE]

to the Newtonian monopole is straightforwardly averaged over one full orbital revolution by using the Keplerian ellipse as reference unperturbed orbit. In Equation (A5), $\mathcal{P}_{\ell}\left(\xi\right)$ is the Legendre polynomial of degree $\ell$. As a result, two kind of averaged, long-term effects occur: secular precessions, explicitly displayed in Section A.1 and labelled with a superscript “s”, and long-periodic signatures, not shown here, having a harmonic pattern characterized by a frequency which is an integer multiple of that of perigee $\omega$. In the calculation, the Earth’s symmetry axis $\boldsymbol{\hat{S}}$ is assumed to be aligned with the reference $z$ axis; moreover, no a-priori simplifying assumptions concerning the orbital geometry of the satellite were made at all.

A.1 Secular effects

[TABLE]