General Relativity Measurements in the Field of Earth with Laser-Ranged Satellites: State of the Art and Perspectives

Abstract

Recent results of the LARASE research program in terms of model improvements and relativistic measurements are presented. In particular, the results regarding the development of new models for the non-gravitational perturbations that affect the orbit of the LAGEOS and LARES satellites are described and discussed. These are subtle and complex effects that need a deep knowledge of the structure and the physical characteristics of the satellites in order to be correctly accounted for. In the field of gravitational measurements, we present a new measurement of the relativistic Lense-Thirring precession with a $0.5\%$ precision. In this measurement, together with the relativistic effect we also estimated two even zonal harmonics coefficients. The uncertainties of the even zonal harmonics of the gravitational field of the Earth have been responsible, until now, of the larger systematic uncertainty in the error budget of this kind of measurements. For this reason, the role of the errors related to the model used for the gravitational field of the Earth in these measurements is discussed. In particular, emphasis is given to GRACE temporal models, that strongly help to reduce this kind of systematic errors.

1. Introduction

Geodetic satellites, such as the two LAGEOS (LAser GEOdynamic Satellite), represent powerful probes for the study of the Earth’s gravitational field, so to derive significant information on its internal structure, but also to test fundamental physics, such as the tiny relativistic precessions that produce secular effects in some of their orbital elements. Laser ranging to Corner Cube Retroreflectors (CCRs) is a technique used since 1964 to track both Earth’s orbiting satellites and the Moon (via the CCRs placed on its surface by the Apollo and Luna missions). It has provided a number of impressive results in the study of both Earth and Moon geophysics, as well as in gravitational physics and in General Relativity (GR). These powerful laser tracking techniques are known, respectively, as Satellite Laser Ranging (SLR) [1,2] and Lunar Laser Ranging (LLR) [3,4,5]. The observable is represented by the round-trip time of short laser pulses measured by means of very precise time devices down to a precision of about 10 ps or less in the case of the SLR, and about 100 ps in the case of the LLR1. The root-mean-square (RMS) precision for the corresponding range measurements is as small as 1 mm for SLR and about 10 mm for LLR.

It is worth mentioning that both SLR and LLR are two relevant techniques, among others, of the observational space geodesy. Indeed, Very Long Baseline Interferometry (VLBI), Global Navigation Satellite Systems (GNSS), Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS), together with SLR and LLR constitute the Global Geodetic Observing System (GGOS) [6]. Within GGOS, improvements in technology and in modelling will produce advances in geodesy and geophysics as well as in GR measurements. Therefore, space geodesy and GR are not independent, but tightly related, and improvements in either field produce benefits in the other in a positive loop. For instance, an improvement of one order of magnitude is expected, in the near future, on the accuracies of the observational components of space geodesy [7]. This would imply a tenfold advancement in a series of closely related fields, such as: modelling, measuring techniques and their accuracy, reference frames and their realization, ground stations position.

In recent years there has been a renewed interest in the geophysical and geodetic applications of the predictions of GR, thanks to increasingly precise gravimetric measurements and the ease of dissemination of precise time measurements via atomic (transportable) clocks and GNSS. In fact, relativistic geodesy [8] is becoming increasingly important and with applications more and more close to everyday life [9,10,11,12,13,14]. Consequently, from these advances, the final accuracy of fundamental physics measurements using the main space geodesy techniques will also benefit.

In this paper we deal with SLR. In particular, we will focus on the main issues regarding the GR measurements by laser tracking to geodetic and passive satellites [15,16,17,18,19,20,21,22,23,24]. These issues concern the models and the estimate of the systematic sources of error in view of new measurements and tests of gravitational physics that can be achieved in the future with this powerful technique [25]. SLR accuracy is currently limited by uncertitude on the satellite centers of mass offset, atmospheric refraction modelling and stations biases. Improvements in these fields will allow a sub-mm precision in the RMS of the SLR range measurements with significant benefits in both geophysical and GR measurements. A fundamental role in this context is played by the dynamical models used for the reconstruction of the satellite orbits. In particular, modelling both gravitational and non-gravitational perturbations in a reliable manner constitutes a delicate task, due to the subtle and complex effects they produce. Moreover, understanding the perturbations is an important prerequisite for an estimate of the systematic errors—in the case of new measurements of relativistic effects—that will need to be reliable, robust and accurate.

The activities of the LARASE (LAser RAnged Satellites Experiment) research program are here described and discussed. The final goal of LARASE is to measure relativistic effects in the weak-field, slow-motion limit and to compare the predictions of GR with those of alternative theories of gravitation in near-Earth environment [26,27,28]. In particular, we devote special attention to the precession of the orbital plane of a satellite produced by the so-called Lense-Thirring effect [29,30]. The Lense-Thirring effect is a consequence of the gravitomagnetic field of Einstein’s GR, which takes part in generating spacetime curvature, together with the gravitoelectric field produced by mass [31]. Indeed, in GR mass currents contribute, along with mass, to determining the trajectory of a test body, that is, to its motion along geodesics in spacetime. Gravitomagnetism [32,33] plays a special role in the astrophysics of compact objects [34] and it has (possible) cosmological consequences linked to Mach’s Principle [35]. To reach this objective, and to exploit at best the available precise and accurate laser range measurements, we started a deep review of the non-gravitational perturbations (NGP) that act on the orbit of the two LAGEOS and LARES (LAser RElativity Satellite) satellites, as well as an analysis of the various models that can be used for the main gravitational perturbations.

The paper is organized as follows: in Section 2, we describe the main NGP and the actions taken to improve the dynamical models used to describe their subtle effects. We focus on the spin model of the satellites, the drag from the neutral atmosphere and the Yarkovsky-Schach thermal effect. In Section 3, our recent results for the precise orbit determination (POD) of the considered satellites are shown. In Section 4, a new measurement of the Lense-Thirring effect is presented with recent results regarding the evaluation of the errors related to the gravitational field modelling. In Section 5, the state of the art of the relativistic measurements performed so far by means of laser-ranged satellites is summarized. Finally, in Section 6 our conclusions and recommendations on the topics discussed are provided.

2. Models for the Non-Gravitational Perturbations

In this Section we introduce and discuss some of the improvements achieved in the development of new models for the main NGP that perturb the orbits of the two LAGEOS [36] satellites as well as that of LARES [37]. We focus on a more detailed description of the internal structure of the satellites, in particular on their moments of inertia, on the spin dynamics and, finally, on the perturbations due to neutral drag and Yarkovsky-Schach effect. It is important to stress that the improvements in the modelling of the NGP are not only important in themselves—for a better understanding of the dynamics of orbits and the interaction of the satellites with the surrounding environment—and for the institutional applications of the considered geodetic satellites in the fields of geodynamics and geophysics. A further improvement of their modelling, in particular of the effects due to the thermal thrust perturbations, is essential to limit the overall error budget to a fraction of a percent compared to the current state of the art in the case of relativistic measurements and, in particular, in the case of the Lense-Thirring effect. This will align the accuracy of the relativistic measure, that is, the estimation of the systematic errors related to the main gravitational and non-gravitational perturbations, with the precision of the measurement itself. The key point in this regard is related to the possibility of using, as observable, the argument of pericenter $\omega$ of the satellites—in particular of LAGEOS II because of its larger eccentricity—in addition to the right ascension of their ascending node $\mathsf{\Omega }$, see Section 2.4 and Section 4.

2.1. Internal Structure

Starting from the original technical documentation of the LAGEOS and LARES satellites, we succeeded in reconstructing information about their structure, materials, and moments of inertia [38]. Through this re-analysis, we obtained an independent estimate of the moments of inertia of the satellites. The amount of deviation from spherical symmetry of a satellite is a relevant quantity that needs to be calculated (no measurements were made or released for the two LAGEOS), since it produces a gravitational torque, similar to the Hipparchus precession of the Earth, that contributes to the spin evolution of the satellite (see Section 2.2). The analysis reported in Reference [38] also impacts on the development of a improved thermal model, needed to properly account for the perturbation produced by the thermal thrust effects (see Section 2.4) which, in turn, depend on the spin evolution.

In

Table 1. Principal moments of inertia of LAGEOS, LAGEOS II and LARES in their flight arrangement.

Satellite	Moments of Inertia (kg m${}^2$)
	${\mathit{I}}_{\mathit{zz}}$	${\mathit{I}}_{\mathit{xx}}$	${\mathit{I}}_{\mathit{yy}}$
LAGEOS	$11.42\pm 0.03$	$10.96\pm 0.03$	$10.96\pm 0.03$
LAGEOS II	$11.45\pm 0.03$	$11.00\pm 0.03$	$11.00\pm 0.03$
LARES	$4.77\pm 0.03$	$4.77\pm 0.03$	$4.77\pm 0.03$

we show the moments of inertia of the LAGEOS and LARES satellites, as computed in Reference [38]. As we can see, the two LAGEOS have almost the same oblateness, about 0.04.

A reconstructed section view of the two LAGEOS satellites is shown in Figure 1, where the two aluminum hemispheres are visible with the section of the cavities containing the CCRs together with the internal brass cylinder and the copper-beryllium shaft.

Our re-analysis was also instrumental to understand and solve some persistent disputes found in the historical literature of the older LAGEOS about some of its characteristics, such as the internal dimensions, the material of the inner cylinder and that of the shaft. We refer to Reference [38] for further details. In the case of LARES, adopting for its mass the official value of 386 kg and assuming a spherical mass distribution, we obtained moments of inertia of $4.77\pm 0.03$ kg m${}^2$. An oblateness as small as 0.002 is however possible.

2.2. Spin Dynamics

The spin of the satellites plays a major role in the modelling of thermal effects, since these depend on the satellite attitude with respect to the radiation sources (see Section 2.4). Models for the spin evolution of the LAGEOS satellites have been developed since 1991 [39,40,41], but they work successfully only in the fast-spin approximation: they provide solutions based on averaged equations, that are valid as long as the satellite rotation period is much smaller than the orbital one. A new comprehensive model, named LASSOS (LArase Satellites Spin mOdel Solutions), describes the evolution of the spin of the LAGEOS and LARES satellites [42], providing solutions for the spin evolution valid for any value of the satellite rotational period. In the fast-spin limit, LASSOS correctly reproduces the results of previous models.

In Figure 2 and Figure 3 we show, as an example for LAGEOS II, the two equatorial components ($S_x,S_y$) of the spin direction of the satellite at different times, the time evolution of the projection $S_z$ of the spin direction along the vertical axis, and the time evolution of the satellite rotational period P, respectively. The time span considered is about 20 years since the launch of LAGEOS II. In these figures we compare the predictions of the model (blue dots) with the available observations (in red those spectrally determined from the SLR data [43], and in black the collection of previous observations as reported in Reference [44]).

The agreement between our model and the experimental data is quite good over the time span where the observations of spin orientation and rate are available. Both the observations (when available) and the model prediction show no evidence of chaotic behaviour. We refer to Reference [42] for further details on the LASSOS model as well as for the results regarding LAGEOS and LARES.

2.3. Neutral Drag

Since the orbit of LARES is much lower than that of the two LAGEOS, with a height of about 1450 km vs 5900 km, a stronger impact of the drag from the atmosphere on its orbit is expected [25,45]. Therefore, an improved analysis of the effects of the neutral atmosphere on the orbit of LARES was necessary. To this purpose, we used as a benchmark the evaluation of the decay of the semi-major axis of the satellite. In this particular analysis we used both the software packages GEODYN II [46] and a modified version of SATRAP (SATellite Reentry Analysis Program) [47,48]. The main feature of SATRAP is the use of different models for the Earth’s atmosphere together with the appropriate geomagnetic and solar activity indices. In this way we directly investigated the impact of the neutral drag on the orbit of LARES, using several of the best available models for the Earth’s atmosphere [24,49].

We performed a POD of LARES over a time span of about 3.7 years (from 6 April 2012 to 25 December 2015). In the set of dynamical models of GEODYN we did not include the neutral drag, nor the thermal effects. We measured an orbital decay of the satellite semi-major axis of about 2.74 mm/d (i.e., close to 1 m/yr) with a $2\sigma$ error of about $\pm 0.02$ mm/d. Such decay corresponds to a transversal mean acceleration T of about $(-1.444\pm 0.011)\times {10}^{-11}$ m/s${}^2$. We then included the neutral drag in the dynamical model and, by means of a linear least-squares fit of the tracking data, we estimated the drag coefficient $C_{\mathrm{D}}$ of the satellite. We obtained for this coefficient a mean value of 3.9 on the time span of the analysis.

An additional analysis was then performed, using SATRAP, over the same time span with the goal of obtaining an independent estimate of the drag coefficient. Among the numerous atmospheric models available, we considered the Jacchia-Roberts 1971 model (JR-71) [50], the Mass Spectrometer and Incoherent Scatter radar 1986 model (MSIS-86) [51], the Mass Spectrometer and Incoherent Scatter radar Extended 1990 model (MSISE-90) [52], the NRLMSISE-00 model [53], developed at the US Naval Research Laboratory (NRL) and the GOST-2004 model [54], issued by the State Committee on Standardization and Metrology of the Russian Federation. Two of these models, MSIS-86 and JR-71, are also implemented in GEODYN.

In this analysis the measured decay was accounted for through the unmodeled transversal acceleration previously estimated with GEODYN; here, considering 3 independent atmospheric models we obtained $C_{\mathrm{D}}=4.0\pm 0.2$, as discussed in Reference [49]. The error provided takes into account the variability of the result according to the atmospheric model considered in the analysis. However, considering the same models implemented in GEODYN, that is, MSIS-86 and JR-71, the differences were of the order of $1\%$ or less.

This analysis proves that the current models developed to describe the neutral atmosphere can successfully reproduce the observed decay of LARES, within their intrinsic errors (around 15%) and their range of applicability. However, after modelling the neutral drag in GEODYN, a very small decay ($\approx -0.04$ mm/d) was still present in the integrated residuals (measured values minus predicted ones) of the LARES semi-major axis, corresponding to a residual, unaccounted for, transversal acceleration of about $-2\times {10}^{-13}$ m/s${}^2$ [49]. Therefore, while, for the two LAGEOS, the Earth-Yarkovsky thermal drag and charged particle drag [55,56] are known to account for $\approx 90\%$ of the observed decay, our analysis shows that the drag from the neutral atmosphere can explain $\approx 99\%$ of the observed decay of the LARES semi-major axis. In this context, a new analysis is needed to determine the contribution of the drag and of the thermal effects to the residuals of the orbital elements of LARES. We have recently undertaken this study, extending it also to the two LAGEOS satellites. In Figure 4, a new measurement for the decay of LARES semi-major axis is shown over a period of about 5.8 years, from 13 April 2012 (MJD 56030). The measured decay is about 2.44 mm/d ($\pm 0.02$ mm/d at $2\sigma$) and corresponds to an average transversal acceleration of about $(-1.289\pm 0.011)\times {10}^{-11}$ m/s${}^2$.

The value obtained for the semi-major axis decay over the longer time span of this analysis is smaller than that determined in the first analysis of about 3.7 years: this is easily explained by correlating these measurements with solar activity: indeed, the first was found during a maximum of the Sun activity, while the extended time span is characterized by a sharp decrease, that is probably approaching the deepest minimum since these measurements are performed2, see Figure 5. This indirectly confirms, once more, how sensitive these measurements, based on SLR data, are.

Since we are interested in the Lense-Thirring effect measurement, the performed analysis on the satellites orbit allows us to explicitly compute the effects of the neutral drag perturbation on the orbit of the considered satellites and, in particular, on the right ascension of the ascending node (RAAN) of each satellite. It is the out-of-plane component W of the neutral drag acceleration that affects directly this orbital element. Equation (1) provides the corresponding Gauss equation for the precession of the RAAN:

$$\dot{\mathsf{\Omega }}=\frac{W}{Hsini}rsin(\omega +f),$$

(1)

where H represents the osculating angular momentum per unit of mass, r the satellite distance from the Earth’s geocenter, while i and f are, respectively, the satellite inclination and true anomaly. In Figure 6 and Figure 7, the out-of-plane acceleration and the corresponding perturbing effects on the RAAN of LARES on a time span of about 5.8 years after the launch of the satellite are shown.

The time evolution of the perturbing acceleration has been computed with SATRAP, as previously described, using the MSIS-86 model for the neutral atmosphere. As we can see, this component of the drag perturbing acceleration is very small. It is characterized by a main oscillation at the orbital period of the satellite (about 1.9 h), further modulated by a long period oscillation of about 133 days. The maximum peak amplitude is a few pm/s${}^2$, with a mean value (on the considered time span) of about $-5.3\times {10}^{-15}$ m/s${}^2$, negligible with respect to the average value obtained for the transverse acceleration T. The effects produced by the out-of-plane acceleration on the precession of the RAAN of LARES are negligible. The oscillations at half of the orbital period of the satellite are characterized by amplitudes not larger than a few ${10}^{-3}$ mas/d, with a mean value on the considered time span of about $-1.4\times {10}^{-6}$ mas/d. As shown in Figure 7, a smoothing at 1 day (the yellow line) clearly highlights this very small impact, with a mean value practically coincident with the one computed with the non-smoothed data.

Beside this direct effect on the satellite’s node due to the out-of-plane acceleration produced by the neutral drag perturbation, it is worth of mention that there are also two indirect effects to be considered. One of these arises from the effects of the drag perturbation on the satellite inclination i. As well known in the literature, the inclination is subject to a secular effect produced by the out-of-plane acceleration due to the neutral drag [58,59,60,61]. However, this effect enters in Equation (1) by means of a trigonometric function and, consequently, does not produces any secular effect in the time behaviour of the precession of the RAAN of the satellite [62]. A second indirect effect on this observable, and in principle more significant than the one produced by the secular variation of the inclination, is that produced by the decay itself of the satellite semi-major axis. We shall return on this feature in Section 4.2.

2.4. Thermal Effects

An intricate role, among the NGPs, is played by the thermal effects produced by the radiation emitted from the satellite surface as consequence of the non-uniform distribution of its temperature. In the literature of the older LAGEOS satellite this problem has been considered since the early 1980s to explain the (apparently) anomalous behaviour of the along-track acceleration of the satellite, characterized by a complex pattern [44,56,63,64,65,66,67,68,69,70,71,72,73,74]. The solution of this dynamical problem is non trivial and requires taking into account the following main aspects:

• a complete physical characterization of the various elements that constitute the satellite, that is, emission and absorption coefficients, thermal conductivity, heat capacity, thermal inertia, …;
• knowledge of the rotational dynamics of the satellite, that is, of its spin orientation and rate;
• a reliable model for the radiation sources: Sun and (especially) Earth.

We have tackled the problem following and expanding the two approaches considered in the past: we first developed a simplified thermal model based on averaged equations, as in References [65,66,67,70], then we developed a general thermal model not restricted to time averages, as in References [44,69].

In the first, simplified case, the thermal model is mainly based on the energy balance on the satellite surface, plus a linear approximation for the distribution of the temperature over the satellite surface (restricted to its CCRs, due to their larger thermal inertia) with respect to its equilibrium (mean) temperature. Conversely, in the second (more general) case the thermal model assumes (i) a satellite (metallic structure) in thermal equilibrium, (ii) the CCRs rings at the same temperature of the satellite metallic body, and (iii) for each CCR, the thermal exchange with the satellite is considered.

In the case of the two LAGEOS, in agreement with the description of Section 2.1, for the metallic structure we considered both the external aluminum surface and the internal brass core, plus the copper-beryllium shaft. For both the simplified and general thermal model, a reliable model for the evolution of the spin vector of the satellite (as outlined in Section 2.2) is necessary. Two main thermal perturbations must be taken into account: the solar Yarkovsky-Schach effect [55,67,68,69,70] and the Earth-Yarkovsky (or Rubincam) effect [55,56,64].

The first effect is characterized by an anisotropic emission of thermal radiation that arises from the temperature gradients across the surface; such gradients are produced by the solar heating and the thermal inertia of the various elements of the satellite (mainly the CCRs). This perturbation is responsible for long-term effects on the orbit of a satellite due to modulation of thermal radiation by the eclipses. For the second effect, the temperature gradients responsible for the anisotropic emission of thermal radiation are produced by the Earth’s infrared radiation. Most of the effect is due to the CCRs and their thermal inertia. This perturbation is also responsible for long-term effects on the orbital elements of a satellite.

In particular, the Yarkovsky-Schach effect strongly influences the argument of pericenter of the LAGEOS [72] and LAGEOS II [74,75] satellites, and this is the reason why this observable is not currently used in the measurements of relativistic effects. Indeed, as briefly highlighted in the introductory part of Section 2, a significant improvement of the model for this perturbation will allow the reintroduction of the satellite(s) argument of pericenter [15,19,76,77,78,79] in the analysis of the orbit residuals that have to be performed for the measurement of the relativistic precessions [21,22,80,81,82,83]. In this context, the benefit will be twofold: (i) it will reduce the contribution of the unmodelled periodic effects in the orbital residuals, allowing for a better separation of the periodic contribution arising from the non-gravitational and the gravitational perturbations, and (ii) it will reduce the overall error budget allowing to remove one more unknown related to the non perfect knowledge of the gravitational field of the Earth, see Section 4.

In the following, we provide as an example our new results for LAGEOS II in the case of the simplified thermal model applied to the solar Yarkovsky-Schach thermal effect. The analysis was performed on a time span of 4080 days, starting from 25 October 1992 (MJD 48920). The LASSOS general (non averaged) model has been used for the satellite spin evolution [42], see Section 2.2. In Figure 8 we show the perturbing accelerations in the Gauss reference frame, that is, the radial (green), transversal (blue) and normal, or out-of-plane (red), components of the Yarkovsky-Schach acceleration. These results are in agreement with those obtained by Reference [44] and based on a more general thermal model, but with an averaged spin evolution.

In Figure 9 and Figure 10, the perturbing effects on the satellite eccentricity and argument of pericenter are shown and compared with the residuals in these orbital elements we obtained by an independent POD.

In this analysis, the EIGEN-GRACE02S solution was used as model of the gravitational field of the Earth [84]. The arc length was 7 days, no empirical accelerations3 were used and the thermal effects (Yarkovsky-Schach and Earth-Yarkovsky) were not modelled. General relativity was modelled with the exception of the Lense-Thirring effect. More details on the analysis strategy can be found in Section 3.

Figure 9 and Figure 10 show a qualitative agreement between the computed effect and the residuals: indeed, the agreement cannot be perfect since, besides the signature of the Yarkovsky-Schach effect, these residuals also contain the signature of the Earth-Yarkovsky effect, as well as that of the asymmetric reflectivity of the satellite’s hemispheres [85,86]. Details of both the simplified and the general thermal model deserves further analysis, applying them to both LAGEOS satellites and comparing the results with independent analyses performed with GEODYN4.

3. Precise Orbit Determination (POD)

The powerful SLR technique provides tracking data for the two LAGEOS and LARES satellites with a RMS precision at the mm level. These range data are known as normal points, and are obtained providing a mean value with a corresponding RMS value by means of an average of the tracking data (full rate data) on a bin size of 120 s, in the case of the LAGEOS satellites, and on a bin size of 30 s for LARES5. Consequently, a comparable precision is needed for the models included in GEODYN and for the models used, a posteriori, for the analysis of the orbital residuals, as in the case of the thermal thrust forces, see Section 2.4. A reliable POD of the two LAGEOS and LARES satellites represents an essential prerequisite for the LARASE objectives. The outcome of this procedure is a precise description of the orbit of each satellite obtained by fitting the tracking data with a suitable set of dynamical models, along with a series of estimated parameters. The orbit so determined is the starting point to compute the residuals in the Keplerian elements [57] of the satellites and then build with them the physical observables that contain the imprint of all unmodeled effects, including the GR ones [76].

Within LARASE, our strategy is to extract the relativistic effects, in particular the GR precessions, not from the range residuals of the satellites—which indeed represent the true observables that we obtain from the data reduction of the SLR data— but from the residuals in the orbit elements. This approach, which is the one followed in all previous measurements of the GR effects with the LAGEOS satellites, has the main advantage of explicitly measuring the orbital effects predicted by GR or by other theories of gravitation [26,27,28]. Traditionally, the first approach is preferred if one is more interested to estimate a (solve-for) parameter related with an effect explicitly modelled in the dynamical model [87].

We used a computation strategy, as customary in the space geodesy community, that consists in splitting the entire analysis time period in smaller intervals, called arcs; for each arc, a POD is performed and the residuals in the Keplerian elements computed as described in Reference [57]. We underline once more that, in this strategy, the relativistic part need not be included in the set of dynamical models, and the corresponding effect on the orbital elements is measured analyzing the time series of the residuals, as stated above.

The modelling setup we are currently using is shown in

Table 2. Models currently used, within the LARASE research program, for the analysis of the orbit of the two LAGEOS and LARES satellites. The models are grouped in gravitational perturbations, non-gravitational perturbations and reference frames realizations.

Model For	Model Type	Reference
Geopotential (static)	EIGEN-GRACE02S/GGM05S	[84,90,91]
Geopotential (time-varying, tides)	Ray GOT99.2	[92]
Geopotential (time-varying, non tidal)	IERS Conventions 2010	[89]
Third–body	JPL DE-403	[93]
Relativistic corrections	Parameterized post-Newtonian	[88,94]
Direct solar radiation pressure	Cannonball	[46]
Earth albedo	Knocke-Rubincam	[63]
Earth-Yarkovsky	Rubincam	[56,64,65]
Neutral drag	JR-71/MSIS-86	[50,51]
Spin	LASSOS	[42]
Stations position	ITRF2008	[95]
Ocean loading	Schernek and GOT99.2 tides	[46,92]
Earth Rotation Parameters	IERS EOP C04	[96]
Nutation	IAU 2000	[97]
Precession	IAU 2000	[98]

. It accounts for: (i) the satellite dynamics, (ii) the measurement procedure and (iii) the reference frames transformations. In this context, our models comply, wherever possible, with the international resolutions and conventions, such as the International Astronomical Union (IAU) 2000 Resolutions [88] and the IERS Conventions (2010) [89].

In

Table 3. Mean orbital elements of LAGEOS, LAGEOS II and LARES.

$\mathbf{Element}$	Unit	Symbol	LAGEOS	LAGEOS II	LARES
Semi-major axis	[km]	a	12 270.00	12 162.07	7 820.31
Eccentricity		e	0.004433	0.013798	0.001196
Inclination	[deg]	i	109.84	52.66	69.49

, we provide our estimate of the mean orbital elements for LARES and the two LAGEOS satellites, computed from the satellites ephemerides obtained from a POD.

In Figure 11, Figure 12 and Figure 13 we show the results of the POD for the three satellites in two cases: (i) when empirical accelerations have been estimated for each arc (label 0001), and (ii) with no empirical accelerations included in the data reduction (label 0002). These are two of the standard runs that are routinely performed by the LARASE collaboration6. Neither the Lense-Thirring effect nor the thermal effects were modeled. With the colors blue, red and green are shown, respectively, the results for LAGEOS (L1), LAGEOS II (L2) and LARES (LR). The two analyses cover a time span longer than 25 years in the case of the two LAGEOS satellites, from 30 October 1992 (MJD 48925) to 31 August 2018 (MJD 58361).

The first plot of each figure shows the range residuals of the three satellites, while the second plot shows their RMS. The range residuals are computed as mean values over each arc for all the stations that perform the tracking of the satellite. The RMS represents, again for each arc, the weighted RMS (i.e., taking into account the number of observations and the degrees of freedom) of the corresponding (mean) range residual. Clearly, for the second evaluation, with no empirical accelerations, the scatter of the residuals is larger, and so their RMS is.

This comparison underlines two main aspects that we need to account for. First, empirical terms (accelerations in the case of GEODYN: in term of a constant acceleration and of once-per-revolution acceleration over each arc), must be used very carefully during a POD performed to extract geophysical parameters and/or relativistic effects. The risk is to absorb the physical quantities to which we are interested to. Second, the observed (and significant) discrepancy in both the range residuals and in their RMS clearly shows the importance of improving the dynamical model used in the code that performs the data reduction. In

Table 4. POD statistics over the time span of the analyses. The mean (M) and the standard deviation ($\sigma$) of the range residuals and of their RMS are provided for the two analyses performed with and without the use of empirical accelerations (EA).

	With EA (0001)				Without EA (0002)
	Residuals [cm]		RMS [cm]		Residuals [cm]		RMS [cm]
	M	$\sigma$	M	$\sigma$	M	$\sigma$	M	$\sigma$
LAGEOS	$-1.38$	6.13	1.58	0.56	15.25	39.09	4.40	1.87
LAGEOS II	$-0.77$	3.50	1.48	0.39	21.77	46.09	4.44	1.71
LARES	0.30	4.31	3.32	0.57	90.67	148.10	22.63	9.47

, the statistics of the two different analyses for the three satellites in terms of mean (M) and standard deviation ($\sigma$) of the (mean) range residuals and of their (weighted) RMS is given. For the two LAGEOS, the statistics refers to the entire period of the analysis, that is, to about 25 years.

The data of Table 4 clearly show the effectiveness of the empirical accelerations to absorb the effects of unmodeled processes, allowing to fit the orbits of the satellites with an RMS of about 1 or 2 cm for the two LAGEOS and $\gtrsim 3$ cm for LARES. It should also be noted that over shorter, selected time spans—and when all the GR effects are modelled—the RMS level of the POD for the two LAGEOS satellites is better by a factor of two, in the range 0.5–1 cm. We can see from Table 4 that the use of the empirical accelerations is also effective in reducing the range residuals of the PODs to values compatible with a null average.

On the other hand, when no empirical terms are used, the relative inadequacy of the modelling becomes apparent, especially in the case of LARES, whose orbit is not determined to any better than a meter. However, the use of empirical accelerations represents, as said, a workaround to handle unexplained phenomena that affect the POD. These, however, should be carefully evaluated and contextualized from the physical point of view; we believe that it is in the best interest of science (geodesy, geophysics, relativity etc.) to try to understand and model these phenomena, rather than fitting them away. Therefore, we set the goal to improve the dynamical model of the orbit of the satellites, especially for the non-gravitational forces: reliable models will allow us to reduce the use of the empirical accelerations.

4. A New Measurement of the Lense-Thirring (LT) Effect

The goals of the LARASE program in terms of relativistic measurements in the near-Earth environment have been described in Reference [25]. Among these, as mentioned in the Introduction, we devote special attention to the precession of the orbital plane of a satellite produced by the so-called Lense-Thirring (LT) effect [29,30]. This orbit precession, measurable both on the right ascension of the ascending node $\mathsf{\Omega }$ and on the argument of pericenter $\omega$, is produced by the Earth’s rotating mass distribution, i.e by Earth’s angular momentum. This effect is also known as frame-dragging effect, as originally Einstein named it [35].

Concerning the experiments proposed for the measurement of this purely relativistic precession on the orbit of a satellite around the Earth, it is worth mentioning that Yilmaz, in 1959, made the first proposal to launch a dedicated satellite in a polar orbit in such a way to cancel the classical precession of the orbit that arises from the Earth’s oblateness [99], whose uncertainty could mask the small relativistic precession. In 1976, van Patten and Everitt proposed to launch two drag-free and counter-orbiting satellites in nearly polar orbits to measure their (combined) LT precession on the two nodes [100]. Two years later, Cugusi and Proverbio suggested to use the existing laser-ranged satellites—and among them (for the first time) the LAGEOS satellite—for the measurement of the LT precession of their ascending node [101]. A few years later, in 1986, Ciufolini proposed the use of two non-polar satellites in supplementary orbital configuration (i.e., with $i_1+i_2={180}^{\circ }$) and with all the other orbital elements equal [102]. This was later on named the LAGEOS III proposal [103].

Finally, in the context of the past measurements of the terrestrial gravitomagnetic field, a special consideration deserves the Gravity Probe B (GP-B) experiment, a NASA and Stanford University space mission [104,105]. The goal of GP-B was the measurement of the frame-dragging effect on a gyroscope in orbit around the Earth. This is also known as the Schiff effect [106,107]. Indeed, a gyroscope in polar orbit around the Earth, as was the case for the GP-B experiment, is subject to two main relativistic precessions: (i) a frame-dragging effect perpendicular to the orbital plane, and (ii) a drift in the orbital plane, that is, the de Sitter (or geodetic) precession [108]. The core of this cryogenic experiment was based on four gyroscopes mounted in a superfluid helium Dewar at 1.8 K and a telescope (star-tracker). Each gyroscope was constituted by a fused quartz sphere coated with a thin layer of niobium and suspended electrically, spun up by the helium gas, and with a magnetic read-out based on the SQUID technology. GP-B has provided a measurement of the frame-dragging effect with an accuracy of about 19% after an intense and extended analysis effort to remove classical disturbances due to a patch effect. The original goal of GP-B was a $\approx 1\%$ measurement of the frame-dragging effect.

In the following, we will first focus our analysis upon the knowledge of the Earth’s gravitational field and its relationship with the LT effect measurement. Then we will return to the neutral drag perturbation and discuss its impact on the LT effect measurements that involve the RAAN as observable. Finally, we will describe our last measurement of the LT effect with the new elements introduced with respect to the previous measures.

4.1. On the Role of the Background Gravitational Field

In the case of two satellites in exactly supplementary orbital configuration, the relativistic LT precession on the right ascension of their ascending node are equal, while the corresponding secular precession due to the Earth’s deviation from the spherical symmetry are equal and opposite, because they are both proportional to the cosine of the orbit inclination i, see Equations (2) and (4). Indeed, the even zonal multipolar moments of the Earth gravitational potential [109,110,111]—that measure the deviation of the Earth’s mass distribution from spherical symmetry—produce a Newtonian precession of $\mathsf{\Omega }$ with the same signature of the sought for effect [15,20,76,78,112]. These moments are described by the even zonal harmonics ${\overline{C}}_{\ell 0}$; for instance, in the case of the quadrupole moment (${\overline{C}}_{20}$) we have [109]:

$${\dot{\mathsf{\Omega }}}^{\mathrm{class}}=\frac{3}{2}n{\left(\frac{R_{\oplus }}{a}\right)}^2\frac{cosi}{{\left(1-e^2\right)}^2}\sqrt{5}{\overline{C}}_{20},$$

(2)

where $R_{\oplus }$ represents the Earth’s mean equatorial radius, n the satellite mean motion and ${\overline{C}}_{20}$ the normalized quadrupole coefficient7. Uncertainty on these moments can then obscure, as said, the LT signal. Although such experiment with two satellites with supplementary inclinations has not yet been carried out8, the approach based on the use of more than one satellite has proven to be effective, over time. Indeed, the analysis of the combined orbits of the two LAGEOS satellites and of the two LAGEOS plus the LARES—that allow a partial cancellation of the effects of the background gravitational field—have provided the currently best measurements of the relativistic LT precession [19,23,24].

These previous measurements of the LT effect with laser-ranged satellites have been proven to be precise, that is, with a small relative error of the measurement. However, an intense debate exist in the literature (see for instance References [20,112,113,114] and References [115,116]) regarding whether or not they are characterized by a robust and reliable evaluation of the main systematic error sources, in particular those of gravitational origin that ultimately strongly influence the measurement accuracy.

We believe that this was mainly due to the imperfect knowledge of the gravitational field of the Earth, as also highlighted by other authors, but with some differences with respect to previous considerations on the subject [117]. In particular, we believe that the main problem in this context is the inadequate modelling of the time evolution of the even zonal harmonics. This also strongly influence the correct comparison between different models for the gravitational field of the Earth, producing wrong conclusions on their reliability. Moreover, also the non-linear fits performed (in previous measurements of the LT effect) on the integrated residuals of the combined orbits, have contributed to absorb part of the errors due to the uncertainties of the even zonal harmonics. For instance, in Reference [24] we have shown that the sensitivity to the number of periodic effects removed by a non-linear fit was significant, between a minimum of about 3% to a maximum of about 7% of the relativistic LT precession. In what follows, we will show that this result was mainly due to the static value of the third even zonal harmonic, ${\overline{C}}_{60}$, of the GGM05S model on the the considered time span. Moreover, if the time span of the analysis of the satellites orbit (i.e., their data reduction during the POD) is long enough, it will be preferable to perform a simple linear fit of the integrated residuals of the combined right ascensions of their ascending nodes, instead of a more complex non-linear fit.

4.2. On the Role of the Neutral Drag

As briefly introduced in Section 2.3, an indirect effect on the RAAN of a satellite and, in particular, in the case of LARES, is that originates from the decay of the satellite semi-major axis produced by the neutral atmosphere. This decay couples with the classical precession produced by the even zonal harmonics coefficients, in particular by the larger quadrupole coefficient, see Equation (2), producing an additional precession of the RAAN that may mask (in principle) the relativistic LT precession. This effect was highlighted in Reference [120]. Therefore, considering the observed decay of LARES semi-major axis ($a\left(t\right)=a_0+\dot{a}t$), we must modify Equation (2) by adding a corrective term:

$$\delta {\dot{\mathsf{\Omega }}}^{\mathrm{class}}=-\frac{7}{2}\frac{\dot{a}}{a_0}{\dot{\mathsf{\Omega }}}^{\mathrm{class}}t.$$

(3)

This result for the correction to the rate of $\mathsf{\Omega }$ agrees, once integrated, with the correction term in Equation (10) of Reference [120] for the satellite right ascension. However, we derive different conclusions with respect to that paper: the effect considered has only a very small impact on the LT measurement, for the following two main reasons.

First, during the orbit determination, the data reduction is performed by means of a least-squares fit to the SLR data on (many) arcs not causally connected, with a length of 7 days in our analysis. This means that the time t in Equation (3) is equal to 7 days, and should not be considered as a single arc many years long as in Reference [120], see for instance Figure 1 of that paper. Consequently, the multiplicative coefficient to the classical precession provoked by the quadrupole coefficient in Equation (3) is as small as $\approx -8\times {10}^{-9}$ in the case of LARES (considering the entire decay, that is, without modelling the effect during the orbit determination), and smaller for the two LAGEOS because of the smaller decay rate and the larger semi-major axis. Second, since the gravitational field is (of course) modelled in our analysis, it is not the quadrupole coefficient that must be considered for the evaluation of the bias produced in the nodal rate (as implicitly assumed in Reference [120]), but its error or uncertainty. This means that the estimate of Reference [120] for the impact of the considered effect on the nodal rate of the satellites must be re-scaled by a factor of $\delta {\overline{C}}_{20}/{\overline{C}}_{20}\approx {10}^{-7}$, that is with a null impact (in practice) in the measurement of the LT secular precession with the two LAGEOS and LARES satellites.

4.3. A Precise Measurement of the LT Effect

As already anticipated, the RAAN is the observable to be carefully modelled in the POD, being less perturbed by the long-term effects produced by the thermal thrust forces [74,75,121]. Successively, after the data reduction, the residuals in this element [57] have to be computed and analyzed for measuring the relativistic precession produced by the LT effect. For this orbital element (rather, for its rate, or time derivative) the predicted secular term of the LT precession is [29,30]:

$${\dot{\mathsf{\Omega }}}^{\mathrm{LT}}=\frac{2G{J}_{\oplus }}{c^2{a}^3{(1-e^2)}^{3/2}},$$

(4)

where G is the gravitational constant, $J_{\oplus }$ the Earth’s angular momentum and c the speed of light. In

Table 5. Rate (mas/yr) for the secular Lense-Thirring precession on the RAAN of the two LAGEOS and LARES satellites.

Orbital Element	LAGEOS	LAGEOS II	LARES
${\dot{\mathsf{\Omega }}}^{\mathrm{LT}}$	30.67	31.50	118.48

, the value of the relativistic precession on the RAAN of each satellite is given. These rates have been computed using Equation (4) and the mean values of the orbital elements shown in Table 3.

From the above discussion on the errors related to the knowledge of the Earth’s gravitational field, it is clear that combining the measurements of three satellites we can improve the determination of the relativistic precession, strongly reducing the impact of these errors, if we can be confident of the overall accuracy for the background gravitational field. Therefore, following the procedure outlined in Reference [76], we constructed a linear combination of the three (node) observables in such a way to cancel out the effects of two of the lowest even zonal harmonic coefficients [122]: these ($\ell =2,4,6\dots$) carry the largest uncertainty, related both to the static and, partially, to the dynamic part of the Earth’s gravitational field. So, we will compare a derived observable, the linear combination of the rate residuals, $\delta {\dot{\mathsf{\Omega }}}_{\mathrm{comb}}^{\mathrm{res}}$:

$$\delta {\dot{\mathsf{\Omega }}}_{\mathrm{comb}}^{\mathrm{res}}\simeq \delta {\dot{\mathsf{\Omega }}}_{\mathrm{LI}}^{\mathrm{res}}+k_1\delta {\dot{\mathsf{\Omega }}}_{\mathrm{LII}}^{\mathrm{res}}+k_2\delta {\dot{\mathsf{\Omega }}}_{\mathrm{LR}}^{\mathrm{res}},$$

(5)

with the GR prediction for the combined LT precession, ${\dot{\mathsf{\Omega }}}_{\mathrm{comb}}^{\mathrm{LT}}$:

$${\dot{\mathsf{\Omega }}}_{\mathrm{comb}}^{\mathrm{LT}}={\dot{\mathsf{\Omega }}}_{\mathrm{LI}}^{\mathrm{LT}}+k_1{\dot{\mathsf{\Omega }}}_{\mathrm{LII}}^{\mathrm{LT}}+k_2{\dot{\mathsf{\Omega }}}_{\mathrm{LR}}^{\mathrm{LT}}=49.66\mathrm{mas}/\mathrm{yr}.$$

(6)

The two coefficients $k_1\simeq 0.387314$ and $k_2\simeq 0.057262$ are obtained from the solution of a linear system of three equations in three unknowns in such a way to remove from the solution the effects (and the errors) of two, low-order even zonal harmonics. Previous analyses based on the same observables have pursued this strategy by cancelling out the first ($\ell =2$, quadrupole) and second ($\ell =4$, octupole) moment. In this work we have investigated, along with the above-mentioned traditional strategy, also the alternate choice of combining these elements in such a way as to cancel, instead, the errors related to the first, ${\overline{C}}_{20}$, and third, ${\overline{C}}_{60}$ ($\ell =6$, hexadecapole moment), even zonal harmonics. Indeed, from the results of preliminary analyses on available measurements we found that the estimate of the third even zonal harmonic of GGM05S is affected by a larger uncertainty than the second one, ${\overline{C}}_{40}$, at least on the time interval of our analysis9. Figure 14 shows the time evolution of both coefficients in the considered time span, as provided by three independent Analysis Centres10: the discrepancies with respect to the static values in GGM05S (horizontal line in either plot) are of comparable amount $\approx (2.5-4)\times {10}^{-6}$, but the error for ${\overline{C}}_{60}$ is leveraged by a larger coefficient, ${\mathsf{\Omega }}_{60}$, in the solution of Equation (5).

We believe that this aspect has negatively impacted on the results (and on the precision) of previous analysis and measurements of the LT effect. At least on those based on the use of the static values for the low-degree even zonal harmonics of GGM05S and on the application of the combination that removes the uncertainties of ${\overline{C}}_{20}$ and ${\overline{C}}_{40}$, instead of those of ${\overline{C}}_{20}$ and ${\overline{C}}_{60}$ as in this work (see Equation (5)). We also believe that the non-linear fit to the combined orbital residuals, performed for example in both Reference [23] and Reference [24], and done appositely to remove (a posteriori) a few unmodelled periodic terms, has actually helped to mask the effective precision of the relativistic measurement, providing an apparently more precise result than it actually was.

It is well known in the literature that the low-degree coefficients of the gravitational field of the Earth, and in particular the quadrupole coefficient ${\overline{C}}_{20}$ because of its greater amplitude, are characterized by a significant time dependence, due to several phenomena that contribute to the variation of the Earth moments of inertia, with long-period effects of annual and inter-annual periodicity [91,123,124,125]. For this reason (and this is another element of novelty of this work), we used in our POD (unlike the current recommendations of the 2010 IERS Conventions [89]) a value for ${\overline{C}}_{20}$ that accounts for a time dependence, as described by a fit to the solutions provided by the Center for Space Research (CSR) of the University of Texas at Austin [91,125]. The CSR time series from GRACE RL05 are based on monthly estimates11.

Independent PODs for each satellite were performed over a time span of about 4.6 years, starting from 6 April 2012 (MJD 56023). The arc length was set to 7 days, no empirical accelerations were used and the thermal effects were not modeled12. For the Earth gravitational field the GGM05S model [90,91] was used. This is the model currently used by the International Laser Ranging Service (ILRS) community for their analyses, and does not rely on data from the LAGEOS satellites13, in this way we avoid undue correlations in the analysis. The residuals of each POD were combined as shown in Equation (5). In Figure 15 we show the integrated residuals for the combined right ascensions of the ascending nodes, that is, of the observable $\delta {\dot{\mathsf{\Omega }}}_{\mathrm{comb}}^{\mathrm{res}}$.

In the analysis we used the residuals in the RAAN of the satellites to evaluate the coefficients ${\overline{C}}_{20}$ and ${\overline{C}}_{60}$ together with the relativistic precession, in such a way to cancel the effects of their uncertainties on the measurement of the GR precession, as discussed above. The combined residuals shown in Figure 15 are clearly characterized by the presence of unmodelled periodic effects. A spectral analysis of the residuals in the RAAN of the three satellites refers these periodic effects to the unmodelled thermal forces (in part) and to some mismodelling of the gravitational effects, due both to the background field and to tides [126,127]. Some of the main spectral lines are summarized in

Table 6. Periodic effects on the RAAN of the LAGEOS and LARES satellites. The periods have been rounded to the closest integer number of days (the different phases were not considered). The rate $\dot{\lambda }$ refers to the time variation of ecliptic longitude of the Earth around the Sun. Some of the thermal spectral lines are also common to solid and ocean tides.

	LAGEOS	LAGEOS II	LARES
Thermal effects	Period [days]	Period [days]	Period [days]
$\dot{\mathsf{\Omega }}$	1052	570	211
$2\dot{\mathsf{\Omega }}$	526	285	105
$\dot{\lambda }$	365	365	365
$2\dot{\lambda }$	183	183	183
$2(\dot{\mathsf{\Omega }}-\dot{\lambda })$	280	111	67
$\dot{\mathsf{\Omega }}+\dot{\lambda }$	271	953	497
Solid tide 165.565	911	622	217
Ocean tide 163.555	221	138	98

.

The inclusion of some of these unmodeled periodic terms in a more general, non-linear, fit does improve the ability to reproduce the time series of the residuals. However, as long as we lack a satisfactory model for such effects, in particular of the thermal effects, this does not improve our understanding of the observed phenomena14 and, in particular, our capability to disentangle efficiently the contribution to the periodic effects due to the non-gravitational perturbations from those that have a gravitational origin. Therefore we shall focus, in this paper, on the simple linear fit of the integrated residuals.

In Figure 15, we also compare the time integral of the $\delta {\dot{\mathsf{\Omega }}}_{\mathrm{comb}}^{\mathrm{res}}$ observable for the right ascensions of the ascending nodes with its relativistic prediction: the black line in the plot represents the GR prediction for the combined nodes, that is, the integration in time of Equation (6), while the red line is the result of a linear fit to the residuals. In the scale of the plot, the two lines almost overlap. For the measured precession ${\dot{\mathsf{\Omega }}}_{\mathrm{comb}}^{\mathrm{meas}}$, that is, for the slope of the fitting line, we obtained:

$${\dot{\mathsf{\Omega }}}_{\mathrm{comb}}^{\mathrm{meas}}\simeq 49.98\pm 0.51\mathrm{mas}/\mathrm{yr},$$

(7)

where the error is at $2\sigma$, and with a fractional discrepancy from the GR prediction of Equation (6) smaller than one percent:

$$\frac{{\dot{\mathsf{\Omega }}}_{\mathrm{comb}}^{\mathrm{meas}}-{\dot{\mathsf{\Omega }}}_{\mathrm{comb}}^{\mathrm{LT}}}{{\dot{\mathsf{\Omega }}}_{\mathrm{comb}}^{\mathrm{LT}}}\simeq 6.4\times {10}^{-3}.$$

(8)

This result represents a precise measurement for the Lense-Thirring precession achieved by means of a simple linear fit to the residuals of a POD15.

The technique of combining the relativistic observables [76], as summarized by Equation (5), allows for a reduction of the impact of some of the errors related to both the static and dynamic part of the Earth gravitational field. However, applying the same combination to the errors associated with the remaining (and not cancelled) coefficients of the Earth gravitational field, is not sufficient in general to provide a robust and reliable error budget for the relativistic measurement, even if we deal with calibrated and not simply formal errors. This is true, in particular, when using a static solution for the gravitational field, whose coefficients represent an average, over a certain time interval, of the temporal trends of the same coefficients measured by GRACE. Therefore, sensitivity analyses in dedicated PODs, and with different input conditions for the background gravitational field, must be performed to have a reliable estimate of the accuracy related to the knowledge of the gravitational field of the Earth [128]. Obviously, the use of the combination of Equation (5), applied to well-calibrated errors, will have to provide a result in agreement with that of a well-executed sensitivity analysis, when the dynamic model includes the time dependency for the gravitational field coefficients. Of course, a careful and thorough analysis of the various sources of systematic errors is therefore mandatory.

Indeed, within the LARASE activities, we have undertaken a number of actions to provide a correct estimate of these systematics that will result in a robust and reliable error budget [24,25,129]. For instance, we started an activity based on a wider use of GRACE data to account for the time dependence of the main even zonal harmonics. In particular, we are fitting the GRACE monthly solutions for these coefficients and use these fitted values in our PODs16. Preliminary results in this sense show a good agreement between the results of the sensitivity analyses and the error estimated by means of the calibrated errors of GGM05S [128]. Under these circumstances, we are confident that the application of Equation (5) to the calibrated errors of GGM05S provides a reliable result for the accuracy of the errors related to the background gravitational field for the measurement of the LT effect described above. By applying Equation (5) to the calibrated errors of GGM05S up to the degree $\ell =20$ (computed conservatively as sum of their absolute values) we obtain, for the static errors:

$$\delta {\dot{\mathsf{\Omega }}}_{\mathrm{comb}}^{\mathrm{LT}-\mathrm{stat}}\simeq 1.8\%{\dot{\mathsf{\Omega }}}_{\mathrm{comb}}^{\mathrm{LT}},$$

(9)

when the third harmonic ${\overline{C}}_{60}$ is removed along with the first one ${\overline{C}}_{20}$.

Finally, in the post-fit analysis we are performing (after the SLR data reduction) to measure the LT precession, we are also estimating the corrections to the low degree even zonal harmonics coefficients removed by the adopted combination together with their degree of correlation with the measured LT effect [128]. In Figure 16 and Figure 17 we show, for the analysis and measurement described in this paper, the correlations plots among the three estimated parameters, that is, a LT parameter $\mu$ (with $\mu =1$ if GR is correct and $\mu =0$ in Newtonian physics), and the corrections $\delta {\overline{C}}_{20}$ and $\delta {\overline{C}}_{60}$ to the quadrupole and hexadecapole coefficients17.

In

Table 7. Correlations among the three estimated quantities obtained from the solution of a system of three equations in three unknowns: the correction $\delta {\overline{C}}_{20}$ to the quadrupole coefficient, the correction $\delta {\overline{C}}_{60}$ to the hexadecapole coefficient and, finally, the Lense-Thirring parameter $\mu$.

	$\mathit{\delta }{\overline{\mathit{C}}}_{20}$	$\mathit{\delta }{\overline{\mathit{C}}}_{60}$	$\mathit{\mu }$
$\delta {\overline{C}}_{20}$	+1.000	+0.354	$-0.024$
$\delta {\overline{C}}_{60}$	+0.354	+1.000	$+0.278$
$\mu$	$-0.024$	$+0.278$	+1.000

, we conversely show the values obtained for the correlations coefficients among the three estimated quantities.

As we can see from these results, while the correlation between the LT parameter and the correction to the quadrupole coefficient is small, neither the correlation between the two corrections to the even zonal harmonics coefficients nor the correlation between the LT parameter and the correction to the hexadecapole coefficient is negligible. However, the values obtained for the degree of correlation among the three parameters represent an improvement if we compare them with those we obtained in previous analyses [24,45,129], and these correlations are further reduced in the analyses we are currently performing with a greater number (in the a-priori dynamical model) of time dependent coefficients for the even zonal harmonics [128].

5. State-Of-The-Art of Relativistic Measurements with Laser-Ranged Satellites

In this section we summarize the main results obtained so far for the relativistic measurements performed by the use of the SLR technique18. This technique has been successfully employed since the 1990s to pursue relativistic tests in the weak-field and slow-motion limit of metric theories [26]. The main achievements have been the measurements of the precession of the orbit of satellites produced by both the gravitoelectric and gravitomagnetic field of the planet. Both effects, that involve two of the three Euler angles defining the orbit orientation with respect to an inertial frame, have been successfully measured: the former is the gravitoelectric, or Einstein, or Schwarzschild precession [21,22,80], while the second is the gravitomagnetic, or Lense-Thirring precession, the main topic of the current paper, measured since 1996 [15].

From the historical point of view, it is worth of mention that Rubincam, in 1977, was the first who proposed the analysis of LAGEOS orbit to provide a direct measurement of the Schwarzschild precession on an Earth-orbiting satellite [132]. Unfortunately, at that time, the overall accuracy (due to observations plus models) was not enough to extract from the data analysis the sought for precession. Indeed, a long time elapsed with, in the meantime, several works focused upon this argument [80,81,82], before the first reliable measurement was finally achieved in 2010 [21].

The results of these measurements, compatible with the predictions of GR, have also allowed to constrain a number of alternative theories for gravitation. These alternative scenarios involve possible violations of the inverse-square law [133,134,135], that is, the existence of a new long-range interaction at the scale of the orbit of the satellites [83], but also the existence of non-symmetric or torsional contributions to the gravitational interaction [136,137,138].

In

Table 8. Measurements of the Lense-Thirring secular precession due to the Earth gravitomagnetic field and of the Schwarzschild secular precession due to the Earth gravitoelectric field. The normalized Lense-Thirring precession $\mu$ was obtained: combining the nodes of the two LAGEOS with the pericenter of LAGEOS II in the case of the measurement (1); combining the nodes of the two LAGEOS for the measurements (2) and (3); combining the nodes of the two LAGEOS with that of LARES for the other measurements. In the case of the Schwarzschild precession $ϵ$, the pericenter of LAGEOS II was used as the relativistic observable.

		Lense-Thirring Precession				Schwarzschild Precession
	$\mathit{\mu }-\mathbf{1}$	$\mathit{\delta }\mathit{\mu }$	Year	Model	Ref.	$\mathit{ϵ}-\mathbf{1}$	$\mathit{\delta }\mathit{ϵ}$	Year	Ref.
(1)	$+0.10$	$\pm 0.03$	1998	EGM96	[18]	$+0.28\times {10}^{-3}$	$\pm 2.14\times {10}^{-3}$	2010	[21]
(2)	$-0.01$	$\pm 0.02$	2004	EIGEN-GRACE02S	[19]	$-0.12\times {10}^{-3}$	$\pm 2.10\times {10}^{-3}$	2014	[22]
(3)	$-0.01$	$\pm 0.01$	2004	EIGEN2S	[77,79]
(4)	$-0.006$	$\pm 0.002$	2016	GGM05S	[23]
(5)	$-0.001$	$\pm 0.001$	2017	GGM05S	[24]
(6)	$+0.006$	$\pm 0.005$	2019	GGM05S	This paper

the best measurements obtained so far for both the LT and Schwarzschild precession are shown. The results are expressed in terms of normalized parameters: $\mu$ for the Lense-Thirring precession and $ϵ$ in the case of the Schwarzschild precession. The Table also indicates the background gravitational field used for the considered measurements of the LT effect. In the case of the two measurements of the Schwarzschild precession, the EIGEN-GRACE02S model was used.

In the case of the LT effect measurements, the reported errors $\delta \mu$ are derived by the authors from a non-linear fit to the residuals and have to be considered as 1$\sigma$ formal errors. The only exceptions are the measurements (3) and (6), where a simple linear fit has been performed. In the case of the present measurement (6), the error $\delta \mu$ has been scaled at 1$\sigma$ to better compare with previous measurements. All the above measurements have been also supplemented with an estimate of their systematic errors19. These estimates are however not completely reliable, as discussed above, and were not included in Table 8.

For this reason, and on the basis of previous discussion in Section 4.3, we believe that the precisions of the last two measures, respectively $0.2\%$ for measurement (4) and $0.1\%$ for measurement (5), are to be considered optimistic. In fact, in both measurements the ${\overline{C}}_{60}$ was not cancelled, and the bias introduced in the slope of the integrated residuals was accidentally absorbed by the amplitude of the periodic terms removed by the non-linear fit. We are particularly confident in the correctness of this analysis on the precision of previous measurements in the case of the measurement (5), that was performed by our team 20.

In the case of the measurements of the Schwarzschild precession, the error $\delta ϵ$ is based on a sensitivity analysis of all the parameters involved in the non-linear fit performed, that is, it is not a simple formal error. For these measurements, a detailed error budget of the systematic errors due to the main gravitational and non-gravitational disturbing effects has also been provided [22], at the level of $2.5\%$.

In

Table 9. Summary of the constraints on fundamental physics obtained in Reference [22] with their corresponding errors. The constraints are compared with previous results [80,82,139,140].

Parameter	Values or Uncertainties	Previous Values or Uncertainties
$\left|\alpha \right|$	$\simeq |(0.5\pm 8.0)\pm 101|\times {10}^{-12}$	$\pm 1\times {10}^{-8}$
${\mathcal{C}}_{\oplus \mathrm{LII}}$	$\le {\left(0.003\mathrm{km}\right)}^4\pm {\left(0.036\mathrm{km}\right)}^4\pm {\left(0.092\mathrm{km}\right)}^4$	$\pm {\left(0.16\mathrm{km}\right)}^4$; $\pm {\left(0.087\mathrm{km}\right)}^4$
$|2t_2+t_3|$	$\simeq 3.5\times {10}^{-4}\pm 6.2\times {10}^{-3}\pm 7.49\times {10}^{-2}$	$3\times {10}^{-3}$

, the constraints obtained in Reference [22] for possible alternative theories of the gravitational interaction are shown and compared with previous results. The physical parameters included in this Table are responsible, if they exist in nature, for a secular effect on the argument of pericenter of a satellite. For the measurement of the LAGEOS II pericenter advance, it was possible to carry on an accurate analysis of the systematic errors [22]: we can therefore provide (see Table 9), besides the errors obtained from a sensitivity analysis of the parameters of the fit, as previously described, also the error budget for each constraint: the last of the three values given for each constraint.

The parameter $\alpha$ refers to the strength of a Yukawa-like potential predicted by some (possible) new long-range force21 for a characteristic range $\lambda$ of the order of half the semi-major axis of LAGEOS II, that is, close to the radius of the Earth. The parameter ${\mathcal{C}}_{\oplus \mathrm{LII}}$ refers to a possible additional interaction between the Earth and the satellite that comes to play in the case of a non-symmetric gravitation theory [136]. If a non-vanishing torsional tensor is present, because of non-symmetric connection coefficients [141,142], a modification of spacetime will be produced22. The parameters $t_2$ and $t_3$ are two of a set of parameters that allow to test for torsion and which appear in the metric [138]. We finally remark that the measurement of the Schwarzschild precession, reported in Table 8, corresponds to a measurement of a combination of the post-Newtonian parameters $\beta$ and $\gamma$ [26]: $\epsilon =\frac{|2+2\gamma -\beta |}{3}$.

6. Conclusions and Perspectives

Since its final formulation, in 1916, Einstein’s GR [143] has passed an incredible number of observational and experimental tests [28], and with still impressive results in recent years. In 2016, the first century of GR [144] was celebrated with the first direct detection of gravitational waves by the LIGO interferometers [145], that is, with the last of the many predictions of the theory and providing insight into the dynamical regime of gravitation. This result opens the way to the gravitational astronomy, thus enriching more and more the fields of cosmology and astrophysics. In 2017, a new era of multi-messenger astronomy arose with the first observation of gravitational waves from the merger of a binary system of two neutron stars [146]. Therefore, although several open problems are still there, we can continue to state that GR still represents the standard model for gravitation [147]. The Lense-Thirring effect we have considered in this paper, and more generally gravitomagnetism, plays a special role among the many predictions of Einstein’s theory [35], since it describes the gravitational effects produced by mass-currents that have no classical counterpart.

We have presented some of the results of the LARASE research program which are developed with the final goal to provide precise and accurate measurements of the gravitational interaction in the field of the Earth by means of laser measurements to passive satellites. In this framework, we have undertaken an activity aimed at improving the models of non-gravitational forces.

We completed a new general model (named LASSOS) for the spin evolution of the two LAGEOS satellites and of LARES (Section 2.2) that is used in a new, detailed thermal model of the satellites to properly account for the thermal thrust effects. These perturbing effects are mainly due to the Sun visible radiation (the Yarkovsky-Schach effect) and to the Earth’s infrared radiation (the Earth-Yarkovsky effect) and critically depend on the spin evolution of each satellite.

We presented the results of a simplified thermal model, based on averaged equations: it provides interesting results when its predictions for the evolution of the orbital elements of a satellite are compared with the results for the residuals in the same elements independently obtained by a POD. In Section 2.4, the results for LAGEOS II in the case of the solar Yarkovsky-Schach effect have been shown. Further work is required for a detailed and comprehensive description of both the simplified and the general thermal models.

As highlighted in previous Sections, an improved model for the thermal thrust forces, especially for the solar Yarkovsky-Schach effect, is mandatory to "clean" the LAGEOS II argument of pericenter from the disturbing effects related with these NGP. This result will allow the introduction of one more observable, whose residuals will help in a further reduction of the impact, in the relativistic measurements, of the errors related with the gravitational field of the Earth.

Among the NGP, special attention was devoted to the effects of the atmospheric neutral drag. This perturbation is particularly important for LARES, because of its lower height with respect to the LAGEOS satellites. In Section 2.3, our results for the decay of the semi-major axis of LARES have been shown. The correlation between the decay rate and the varying solar activity has been higlighted. We have shown that almost all of the observed decay of the satellite semi-major axis can be explained by the drag of the neutral atmosphere. Anyway, it seems that a very small residual decay $\lesssim 1\%$ is still unaccounted for. This may be due to the thermal effects [148,149] and/or to the effects of the Coulombian interaction between the satellite potential and the charged particles in its surrounding environment. These issues are currently under investigation. We have also shown that the neutral drag perturbation has no significant impact in the measurement of the Lense-Thirring effect, see Section 4.2.

We presented in Section 4.3 a preliminary result for our last measurement of the Lense-Thirring effect. This new measurement is based on a simple linear fit to the integrated and combined residuals of the two LAGEOS and LARES satellites. A measurement with a discrepancy of 0.6% with respect to the prediction of GR was obtained. In terms of the Lense-Thirring parameter $\mu$, our result is:

$$\mu -1=0.006\pm 0.005\pm \delta {\mu }_{sys},$$

(10)

where the $2\sigma$ error of Equation (7) has been re-scaled to $1\sigma$ for a direct comparison with previous measurements, see Section 5. On the basis of the discussion presented in Section 4, we retain this $0.5\%$ precision reliable.

With regard to the systematic sources of error, $\delta {\mu }_{sys}$, Section 4.1 and Section 4.3 were also devoted to a discussion of the systematic error in the Lense-Thirring effect measurement related to the knowledge of the gravitational field of the Earth. This represents, for relativistic measurements with laser-ranged satellites, the main source of error and also one of the most complex to manage in a reliable way. In this regard, we have shown the importance of the inclusion in the dynamical model of GRACE temporal solutions for the even zonal coefficients, starting from the quadrupole one. The error was set to about $2\%$ of the relativistic precession predicted by GR for the analysed combined orbits, see Equation (9). The GGM05S model was used to model the gravitational field of the Earth.

In this context, a final robust and reliable error budget also needs a correct evaluation of the contribution of the periodic effects that are currently not modelled in the POD and, probably, some improvements in the dynamical model related with the gravitational effects that are responsible of periodic signatures in the right ascension of the ascending node of the two LAGEOS satellites and in that of the LARES one. Consequently, we have prepared a line-up of improvements in the analysis techniques, such as introducing time-dependent models for the zonal harmonics, performing sensitivity analysis runs and improving the thermal models, that will further improve the quality of the measurement for the Lense-Thirring precession, and reduce its error bar. One more important aspect is the estimate of the coefficients of the gravitational field together with the Lense-Thirring effect parameter $\mu$, see Section 4.3. Other sources of systematic errors to be included in the error budget are the uncertainties in the tidal model and in the de Sitter precession.

Concerning tides, both solid and ocean, we have considered their effects and errors in the case of the Lense-Thirring effect measurement in Reference [127]. Solid tides have a negligible impact in this relativistic measurement, while ocean tides may impact (conservatively) at a level of about $\lesssim 0.5\%$ of the relativistic precession. The de Sitter precession was modelled in our POD following [94]. Assuming conservatively an upper bound error of $0.64\%$ [150], the impact in the Lense-Thirring effect measurement is about $0.3\%$ of the relativistic precession.

We conclude underlying that thanks to the efforts of the ILRS, the precise measurements offered by the SLR technique cover a very important role for fundamental physics measurements, beside their “institutional” role in space geodesy, geodynamics and geophysics. Consequently, it is very important that the wealth of the SLR ranging observations be guaranteed and, possibly, improved in the future.