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

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.


Introduction
In regard to possible tests of post-Newtonian (pN) features of general relativity with, e.g., Earth's artificial satellites and solar system's planets, there is a considerable confusion in the literature about the possible use of the mean anomaly M as potential observable in addition to the widely inspected argument of pericentre ω and longitude of the ascending node Ω. Indeed, it is as if some researchers, including myself, 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 η, or they systematically neglected a potentially non-negligible contribution to the change of the mean anomaly induced indirectly by the semimajor axis a through the mean motion n 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 & 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 × 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 & 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 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 satellites, we will show that 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 & 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 ǫ 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 & 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 disturbing 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 M(t) is one of the three 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, it is defined as where 1 η is the mean anomaly at the reference epoch t 0 , and is the Keplerian mean motion. In Equation (2), µ 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 µ = const. The mean anomaly at epoch η is one of the six Keplerian orbital elements parameterizing the orbit of the test particle in space. In the unperturbed case, M(t) is a linear function of time because both a and η are constants of motion. If a relatively small perturbing acceleration A is present, both a and η are, in general, affected by it, becoming time-dependent. As a result, also the mean motion is, in general, modified so that Thus, the perturbed mean anomaly is the sum of the now time-dependent mean anomaly at epoch η(t) and a function of time ̺(t) whose derivative is equal to the (perturbed) mean motion, i.e., The resulting change ∆M of the mean anomaly with respect to the unperturbed case is, thus, where we defined 1 The symbol η is used for the mean anomaly at epoch by Milani, Nobili & Farinella (1987). In the notation by Brumberg (1991), the mean anomaly is l, while the mean anomaly at epoch is l 0 . Kopeikin, Efroimsky & Kaplan (2011) denote η as M 0 , while Bertotti, Farinella & Vokrouhlický (2003a) 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 varying semimajor axis a is so that The shifts ∆a(t) and ∆η(t) can be perturbatively calculated by evaluating the right-hand-sides of the Gauss equations for their rates of change (Bertotti, Farinella & Vokrouhlický 2003a) onto the unperturbed Keplerian ellipse. In Equations (9) to (10), e is the eccentricity, p a 1 − e 2 is the semilatus rectum, r = p/ (1 + e cos f ) is the (unperturbed) distance of the test particle from the primary, and A R , A T are the projections of the perturbing acceleration 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, Depending on the disturbing acceleration, Φ(t) is linear in time if the average over an orbital period P b of its rate of change 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. In general, the calculation of Φ(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 M may not be a wise choice because of the disturbances introduced by Φ(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 data reductions. On the contrary, the mean anomaly at the epoch η, 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, in principle, as an additional piece of information to improve some tests of pN gravity on the same foot of ω and Ω. This fact seems to have gone unnoticed so far in the literature, as in the case of Ciufolini & Pavlis (2005).

The mean longitude
Similar considerations hold for the mean longitude l defined as where is the longitude of pericentre. If a disturbing acceleration A is present, it can be expressed in terms of the mean longitude at epoch 2 ǫ (Milani, Nobili & Farinella 1987;Soffel 1989;Brumberg 1991;Bertotti, Farinella & Vokrouhlický 2003b) as so that its shift is ∆l(t) = ∆ǫ(t) + Φ(t).
The perturbation ∆ǫ(t) of the mean longitude at epoch can be straightforwardly worked out by means of the Gauss equation for its variation (Bertotti, Farinella & Vokrouhlický 2003a) In Equation (18), A N is the projection of the perturbing acceleration A onto the out-of-plane direction, while is the argument of latitude.
3. The secular rates of change of η(t), ǫ(t), Φ(t) for some pN accelerations Here, we will look at the effects due to the standard 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Ŝ to assume any orientation in space.

The 1pN gravitoelectric Schwarzschild-like acceleration
To the first pN order, the relative acceleration for two pointlike bodies of masses m A , m B separated by a distance r and orbiting with relative velocity v is (Damour & Deruelle 1985;Soffel 1989) where c is the speed of light in vacuum, is the total gravitational parameter of the binary system, is the the radial velocity of the relative orbital motion, and In Sections 3.1.1 to 3.1.2 we will work out the effect of Equation (21) on Φ and η, ǫ, respectively.
From it, the rate of change of Φ(t) averaged over one orbital period P b can be straightforwardly worked out as where From the analytical expression of the right-hand-side of Equation (27), 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, Φ(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 (27) on f 0 . We are not able to analytically calculate Equation (26) unless a power expansion in e of Equation (27) 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, whose relevant 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 (27). Figure 1 displays the plot of Equation (27), in milliarcseconds per year mas yr −1 , for the orbital parameters of Table 1, and the numerically produced time series of Φ(t), in mas, over 1 yr for the same orbital configuration; the agreement between the slope of Φ(t) and the area under the curve of Equation (27) is remarkable. From Figure 1, it can be noted that, as expected, the 1pN Schwarzschild-like acceleration induces a secular variation on Φ(t) which has to be added to that affecting η and ǫ displayed in Section 3.1.2. Table 1: Orbital and physical configuration of a fictitious terrestrial geodetic satellite characterized by the frozen perigee configuration. The period P ω of the perigee is mainly determined by its secular precession due to J 3 , J 4 because of the critical inclination which makes the secular precession due to J 2 equal to zero. Since the values ρ LARES = 5.96 × 10 −16 kg m −3 (Pardini et al. 2017), and ρ LAGEOS = 6.579 × 10 −18 kg m −3 (Lucchesi et al. 2015) of the neutral atmospheric density at h LARES and h LAGEOS are known, we, first, used them in where ρ 0 and h 0 are, in general, referred to some reference height, to determine Λ in the case  (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 · Φ(t) may represent an important contribution to the overall long-term rate of change of M(t). Suffice it to say that, in the case of Mercury and the Sun, it is 3.1.2. The mean anomaly at epoch η and the mean longitude at epoch ǫ The Gauss equations for the variation of η and ǫ (Equation (10) and Equation (17)) allow to straightforwardly work out their secular rates of change which turn out to be 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 (32) to (30). In the case of Mercury and the Sun, Equations (32) to (30)

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 S, on a test particle orbiting it with velocity v it turns out that for an arbitrary orientation of the body's spin axisŜ in space. Thus, it is It implies that the claims by 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 b of a satellite are, in fact, erroneous for any spacecraft.

Moreover, it is also
for anyŜ as well. In Equation (37), is the unit vector directed along the orbital angular momentum along the out-of-plane direction, whilem = {− cos I sin Ω, cos I cos Ω, sin I} 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Ŝ, Equation (37) reduces to Incidentally, Equation (35) and Equation (40) 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 ≃ 1 − 0.5 cm level.

The secular rates of change of η(t), ǫ(t), Φ(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Ŝ 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.

The quadrupole mass moment J 2
To the Newtonian level, the external potential of an oblate body at the position r is where J 2 is the first even zonal harmonic coefficient of the multipolar expansion of its classical gravitational potential, ξ Ŝ ·r is the cosine of the angle between the primary's spin axis and the particle's position, and 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 In Sections 4.1.1 to 4.1.2 we will work out its impact on Φ(t) and η, ǫ, respectively.

The shift Φ(t) due to the variation of the mean motion
It turns out that so that Φ(t), which depends on f 0 , is linear in time. It is not possible to explicitly display the analytical expression which we obtained for (1/P b ) dΦ/d f in the case of an arbitrary orientation ofŜ 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 the sky 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Ŝ, we have + 24e 3 cos 3 f cos 2u sin 2 I + 3 cos f e 4 + e 2 (1 + 3 cos 2I) + 24e cos 2u sin 2 I + + 3 −4 2 + 3e 2 + 3e 2 cos 2 f sin 2 f + 8 (1 + e cos f 0 ) 3 sin 2 f 0 sin 2 I sin 2ω.
The numerical value of the area under the plot of Equation (46), depicted in the upper panel of Figure 2, is confirmed by the time series for Φ(t) produced by numerically integrating the equations of motion of the fictitious satellite of Table 1, and displayed in the lower panel of Figure 2.  Table 1 and with f 0 = 228 deg, over a full orbital cycle of the true anomaly f . Its area, giving · Φ in mas yr −1 , turns out to be equal to 3.8 × 10 7 mas yr −1 . Lower panel: Numerically produced time series, in mas, of Φ(t) over 1 yr obtained by integrating the equations of motion in rectangular Cartesian coordinates for the fictitious Earth's satellite of Table 1. The Newtonian acceleration of Equation (44) due to J 2 was added to the Newtonian monopole. As initial value for the true anomaly, f 0 = 228 deg was adopted. The slope of the linear trend amounts just to the area under the curve in the upper panel.

The mean anomaly at epoch η and the mean longitude at epoch ǫ
The Gauss equations for the variations of η and ǫ (Equation (10) and Equation (17)) allow to straightforwardly obtain wherel = {cos Ω, sin Ω, 0} is the unit vector directed along the line of the nodes such that l ×m =ĥ. Also Equations (A7) to (49) 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 (A7) to (49) reduce to · η = 3 n b R 2 e J 2 (1 + 3 cos 2I) 8 a 2 1 − e 2 3/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 b (t) and Φ(t).
For a cannonball geodetic satellite, the drag acceleration can be expressed as In Equation (52), C D , Σ, ρ, 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, V is where Ψ is the Earth's angular velocity. We will model the atmospheric density as where ρ 0 refers to some reference distance r 0 , while Λ is the characteristic scale length. By assuming Λ can be determined as where ρ min = ρ(r max ), 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 calculate the averaged rates of change of its Keplerian orbital elements by keeping ρ 0 fixed during one orbital period P b . An exact analytical calculation without recurring to any approximation in both e and ν Ψ/n b is difficult.
In Sections 4.2.1 to 4.2.2 we will calculate the impact of Equation (52) on Φ(t) and η, ǫ, respectively.

The shift Φ(t) due to the variation of the mean motion
Let us, now, start to look at ∆n b (t) by means of Equation (7). We will show that it is linear in time because . ∆n b 0. The analytical expression of 1/P b d∆n b /d f is where V 2 ( f ) = 1 − ν 2 1 − e 2 3/2 cos I 1 + e 2 + 2 e cos f + ν 2 1 − e 2 3 3 + cos 2I + 2 sin 2 I cos 2u 4 (1 + e cos f ) 2 1 + e 2 + 2 e cos f .
Since it is not possible to analytically integrate Equation (59) with Equation (60) in the most general case without recurring to approximations in e and ν, 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 Figure 3 depicts Equation (59), while the lower panel displays the time series for ∆n b (t) calculated from a numerical integration of the satellite's equations of motion in rectangular Cartesian coordinates over 1 yr.  Table 1 and with f 0 = 228 deg, over a full orbital cycle of the true anomaly f . Its area, giving . ∆n b in mas yr −2 , turns out to be equal to 107, 217 mas yr −2 . Lower panel: Numerically produced time series, in mas yr −1 , of ∆n b (t) over 1 yr obtained by integrating the equations of motion in rectangular Cartesian coordinates for the fictitious Earth's satellite of Table 1. The drag acceleration of Equation (52) was added to the Newtonian monopole. As initial value for the true anomaly, f 0 = 228 deg was adopted. The linear trend is apparent, and its slope amounts just to the area under the curve in the upper panel.
The fact that . ∆n b 0 implies that ∆n b is linear 3 in time and, thus, Φ(t) is quadratic. It is explicitly shown in Figure 4 by the time series calculated for Equation (8) Table 1. The drag acceleration of Equation (52) was added to the Newtonian monopole. As initial value for the true anomaly, f 0 = 228 deg was adopted. The quadratic signature is apparent, and its final value is in agreement with what expected from Figure 3.
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 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 Φ(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 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 & Pavlis (2005). The dependence of Φ(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 · Φ to be added to the further contribution due to · η , · ǫ .

The mean anomaly at epoch η and the mean longitude at epoch ǫ
About the secular rates of η and ǫ, the Gauss equations for their variations allow to obtain 4π e (1 + e cos f ) 4 2 + 3e 2 + 2e 2 + e 2 cos f + e 2 cos 2 f − −ν 1 − e 2 3/2 (2 + e cos f ) cos I , Since it is not possible to analytically integrate Equations (61) to (62) in an exact form, we, first, plot them as functions of f over a full orbital cycle in Figure 5 for the orbital configuration of Table 1, and, then, numerically calculate the areas under their curves in order to obtain · η , · ǫ .  Table 1 and with f 0 = 228 deg, over a full orbital cycle of the true anomaly f . Their areas give · η , · ǫ in mas yr −1 . In this case, they vanish, as confirmed also by a numerical integration of the satellite's equations of motion for the same physical and orbital parameters. Also in this case, a numerical integration of the satellite's equations of motion turns out to confirm such results.

Some possible uses with the LAGEOS and LAGEOS II satellites
Here, we will look at the possibility of using the nodes Ω and the mean anomalies at epoch η of 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 η in addition to Ω 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. 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 solutions 4 releasing the experimentally estimated values of the geopotential's parameters, do not 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 (2013) and references therein. To the benefit of the reader, we review the linear combination approach, which is a generalization of that proposed for the first time by Ciufolini et al. (1996) to test the gravitomagnetic field of the Earth with artificial satellites of the LAGEOS family. It should be noted that, actually, it is quite general, being not necessarily limited just to the Lense-Thirring case. By looking at N orbital elements 5 κ (i) experiencing classical secular precessions due to the even zonals of the geopotential, the following N linear combinations can be written down They involve the pN averaged precessions · κ (i) pN as predicted by General Relativity and scaled by a multiplicative parameter µ pN , and the errors in the computed secular node precessions due to the uncertainties in the first N − 1 even zonals J 2s , s = 1, 2, . . . N − 1, assumed as mismodeled through δJ 2s , s = 1, 2, . . . N − 1. In the following, we will use the shorthand for the partial derivative of the classical averaged precession · κ J ℓ with respect to the generic even zonal J ℓ of degree ℓ. Then, the N combinations of Equation (63)  , i = 1, 2, . . . N of each of the N orbital elements considered. In principle, such residuals account for the purposely unmodelled pN effect, the mismodelling of the static and time-varying parts of the geopotential, and the non-gravitational forces. Thus, one gets If we look at the pN scaling parameter 6 µ pN and the mismodeling in the even zonals δJ 2s , s = 1, 2, . . . N − 1 as unknowns, we can interpret Equation (65) as an inhomogenous linear system of N algebraic equations in the N unknowns while the constant terms are the N orbital residuals It turns out that, after some algebraic manipulations, the dimensionless pN scaling parameter, which is 1 in General Relativity, can be expressed as In Equation (69), the combination of the N orbital residuals 6 In general, it is not necessarily one of the parameters of the parameterized post-Newtonian (PPN) formalism, being possibly a combination of some of them.
is, by construction, independent of the first N − 1 even zonals, being impacted by the other ones of degree ℓ > 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. Instead, combines the N pN orbital precessions as predicted by General Relativity. The dimensionless coefficients c j , j = 1, 2, . . . 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, C δ = 0 if Equation (70) is calculated by posing for any of the first N − 1 even zonals, independently of the value assumed for its uncertainty δJ ℓ .
As far as the Lense-Thirring effect and the satellites LAGEOS and LAGEOS II are concerned, the linear combination of the four experimental residuals δΩ L , δΩ L II , δη L , δη 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 C δ = δΩ L + c 1 δΩ L II + c 2 δη L + c 3 δη L II (73) whose coefficients c 1 , c 2 , c 3 are purposely constructed with the results of Section A.1. They turn out to be -26 -Their numerical values, computed with the satellites' orbital elements inserted in Equations (A6) to (A11), are Thus, the predicted combined Lense-Thirring signature is The combination of Equation (73) is affected by the orbital precessions induced by the fourth even zonal harmonic of the geopotential. The resulting mismodeled combined signal can be evaluated by means of Equations (A12) to (A13) 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 (80) would be much smaller than 1%. Indeed, from, e.g., the zero-tide model Tongji-Grace02s (Chen et al. 2018), it is 7 σ C 8,0 = 1.3 × 10 −14 . It implies a combined mismodeled precessions as little as 0.01 mas yr −1 , corresponding to 0.01% of the combined Lense-Thirring effect. If, instead, the difference ∆C 8,0 between the values of 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 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 η 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(Lucchesi , 2002(Lucchesi , 2003Pardini 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.

Summary and overview
The Newtonian and pN shifts ∆M(t), ∆l(t) experienced by the mean anomaly M(t) and the mean longitude l(t) with respect to their Keplerian temporal linear trends when a post-Keplerian 7 The zonal harmonics J ℓ of the geopotential are connected with its fully normalized Stokes coefficients C ℓ,0 by the relation J ℓ = − √ 2ℓ + 1 C ℓ,0 , ℓ = 2, 3, 4, . . . disturbing acceleration acts upon an orbiting test particle are, in general, due to the perturbations ∆η(t), ∆ǫ(t) of the mean anomaly at epoch η and mean longitude at epoch ǫ, and the change ∆n b (t) in the mean motion n b which, in some cases, can induce a quadratic shift Φ(t) in 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 Φ(t) quadratically; nonetheless, the pN linear trends of interest can be effectively separated from such a competing aliasing effect if a sufficiently long time span for the data analysis is adopted. Thus, also M(t) and l(t) can, in principle, be employed in tests of pN gravity even with passive satellites, not to mention, then, the use of drag-free apparatuses. If, instead, η and ǫ are adopted, such an issue is a-priori circumvented because they are not impacted by the possible change in the mean motion n b . Since they undergo secular precessions due to the even zonal harmonics J ℓ , ℓ = 2, 4, . . . of the geopotential, it is possible, in principle, to use them in combination with, say, the nodes Ω 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 η and ǫ by all the most relevant non-gravitational accelerations should be performed. There are no net Lense-Thirring rates of change of a and n b .
In astronomical and astrophysical binary systems, not plagued by non-gravitational perturbations, using η 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 pN effects on η are often larger than the corresponding pericenter rates.