A 1% Measurement of the Gravitomagnetic Field of the Earth with Laser-Tracked Satellites

Abstract

A new measurement of the gravitomagnetic field of the Earth is presented. The measurement has been obtained through the careful evaluation of the Lense-Thirring (LT) precession on the combined orbits of three passive geodetic satellites, LAGEOS, LAGEOS II, and LARES, tracked by the Satellite Laser Ranging (SLR) technique. This general relativity precession, also known as frame-dragging, is a manifestation of spacetime curvature generated by mass-currents, a peculiarity of Einstein’s theory of gravitation. The measurement stands out, compared to previous measurements in the same context, for its precision ($\simeq 7.4\times {10}^{-3}$, at a 95% confidence level) and accuracy ($\simeq 16\times {10}^{-3}$), i.e., for a reliable and robust evaluation of the systematic sources of error due to both gravitational and non-gravitational perturbations. To achieve this measurement, we have largely exploited the results of the GRACE (Gravity Recovery And Climate Experiment) mission in order to significantly improve the description of the Earth’s gravitational field, also modeling its dependence on time. In this way, we strongly reduced the systematic errors due to the uncertainty in the knowledge of the Earth even zonal harmonics and, at the same time, avoided a possible bias of the final result and, consequently, of the precision of the measurement, linked to a non-reliable handling of the unmodeled and mismodeled periodic effects.

1. Introduction

The theory of General Relativity (GR) [1,2,3] in over a century since its formulation has brilliantly passed many experimental and observational tests, and it presently represents our best understanding of gravitation. In fact, GR is a revolutionary, very profound, and elegant theory for the description of the gravitational interaction that interprets gravity as a manifestation of spacetime curvature. Among the many predictions of Einstein’s geometrodynamics, a very peculiar role is played by Gravitomagnetism [4,5,6,7], which shows us the greater richness of GR compared to Newton’s theory of universal gravitation. Indeed, gravitomagnetism describes gravitational phenomena that are analogous to the magnetic effects of Maxwell’s theory of electromagnetism. This analogy is quite stringent in the so-called weak-field and slow-motion (WFSM) limit of Einstein’s GR. This limit implies: $GM/r{c}^2\ll 1$, $J/Mrc\ll 1$, and $v/c\ll 1$, where G, M, and J represent, respectively, the gravitational constant, the mass of the primary and its intrinsic angular momentum, c and v are the speed of light and that of the satellite, while r represents the satellite distance from the primary. Following this approach, we have a gravitoelectric field ${\mathbf{E}}_{GE}$ produced by masses, analogous to the electric field produced by charges, and a gravitomagnetic field ${\mathbf{B}}_{GM}$ produced by mass-currents, analogous to the magnetic field produced by electric currents. This is gravitoelectromagnetism [8]. Therefore, GR predicts a formal counterpart to the Ampère’s law describing the magnetic effects that are related to the electric currents [9], which does not have an analogous counterpart in Newton’s theory, that (as is well known) has only a counterpart in Coulomb’s law of electrostatics.

Gravitomagnetism describes the gravitational effects that are produced by mass-energy currents. This points to its fundamental role in explaining high energy astrophysical phenomena, as a powerful source in creating (together with magnetohydrodynamic phenomena) the accretion disks and the jets observed in quasars and active galactic nuclei (AGN), due to the presence of a rotating supermassive black hole in their centers [10,11,12,13,14,15,16]. Currently there is not yet an observational direct evidence of the role of gravitomagnetism in the strong-field regime, as the one highlighted above. Nevertheless, in the WFSM limit of GR, two different kind of measurements have been done of the gravitomagnetic field produced by the Earth1.

One of these measurements was achieved by the NASA and Stanford University space mission denominated Gravity Probe B (GP-B) [22,23]: a gyroscope in polar orbit around the Earth is subject to two main relativistic precessions, a drift in the orbital plane, i.e., the de Sitter effect or geodetic precession [24], and a frame-dragging effect perpendicular to the orbital plane [25,26,27,28]. This frame dragging, also known as Pugh-Schiff dragging, is due to the Earth’s rotation, i.e., to gravitomagnetism. GP-B has provided a measurement of the frame-dragging effect and of the geodetic effect in agreement with the GR predictions with accuracies of about 19% and 0.28%, respectively. GP-B has been an independent and unique experiment, very beautiful, but also quite expensive, and unfortunately without any possibility of replicating it in the foreseeable future. The original goal of GP-B was a ≈1% measurement of gravitomagnetism [29].

The second kind of measurement was based on the measurement of the Lense-Thirring effect on the orbits of passive laser-ranged satellites, as the two LAGEOS [30,31,32,33,34,35], and, more recently, also including the orbit of LARES [36,37,38,39,40,41,42]. The Lense-Thirring effect is the precession of the orbit of a satellite around a primary due to the angular momentum (or spin) of the primary itself [43,44]. Contrary to the case of GP-B, these measurements may be repeated many times, allowing for replicability of the experiment from different research groups, with different software for the data analysis, modeling procedures and also using a different methodology in the evaluation of the error budget of the measurement [45,46]. Indeed, although the measurements of the Lense-Thirring effect with the analysis of the orbit of laser-tracked satellites has been subject to a number of criticisms in the past [34,47,48,49,50,51,52], their error budget, in terms of the evaluation of the main systematic sources of error, provides an error level smaller than that of GP-B.

In this paper, we present a new precise and accurate measurement of the Lense-Thirring effect on the combined orbits of the geodetic satellites LAGEOS, LAGEOS II and LARES. This work aims to complement and further extend the results of [41], which were based on a number of new features with respect to previous measurements. These mainly concern (i) the dynamical model used for the gravitational field, (ii) the way to perform the measurement of the relativistic effect from the analysis of the combined orbits of the satellites, (iii) the concomitant estimate of the low-degree coefficients of the Earth’s gravitational field with the unmodeled relativistic precession, and, finally, (iv) the accurate estimate of the systematic errors of the measurement.

This work is part of the activities of the LAser RAnged Satellites Experiment (LARASE) [46] and concludes these activities for what concerns the measurement of the Lense-Thirring effect. Other objectives of LARASE, in the field of fundamental physics measurements [53,54], will flow into the new experiment SaToR-G (Satellites Tests of Relativistic Gravity). The main goal of SaToR-G is to tests relativistic gravity in the WFSM limit and compare the predictions of GR with those of other alternative theories for gravitation to search (possibly) for new physics beyond GR [55,56].

The paper is organized as follows. In Section 2, the precessional effects that are produced by the gravitomagnetic field are briefly introduced with a discussion on the importance of a reliable measurement of gravitomagnetism. In Section 3, the main results regarding the Lense-Thirring effect are concisely mentioned with the estimate of their error budget. In Section 4, we will focus on errors deriving from imperfect knowledge of Earth’s gravitational field and how they have been handled by a careful consideration of their time dependence. In Section 5, the analyses undertaken to perform the relativistic precession measurement will be described. In Section 6, we present the results regarding the measurement of the Gravitomagnetic field of the Earth through the measurement of the Lense-Thirring parameter $\mu$ in two cases: from the residuals of $\mu$ obtained from the joint analysis of the orbits of the satellites and from their integration. In Section 7, the measurement of the relativistic parameter is presented in a statistical context, in order to decorrelate the measurement of the secular effect from the periodic terms not modeled in the data reduction and estimate (in an alternative way) their impact on the measurement itself. Finally, in Section 8, we provide our conclusions on the measurement we performed.

2. Lense-Thirring Effect and the Importance of an Accurate Measurement of the Gravitomagnetic Field

In the WFSM limit of GR, the precession of a gyroscope with respect to the asymptotic reference frame of distant stars can be written as:

$${\dot{\mathbf{\Omega }}}_{GM}=-\frac{1}{2c}{\mathbf{B}}_{GM}=\frac{G}{c^2{r}^3}\left[3(\mathbf{S}·\widehat{r})\widehat{r}-\mathbf{S}\right],$$

(1)

where $\mathbf{S}$ represents the spin of the primary. This particular precession is known as “dragging of gyroscopes” or “dragging of inertial frames”, as Einstein named it2. Indeed, local inertial frames have a peculiar role in Einstein’s theory of GR. These are the frames where locally, in space and, in time, the law of physics are those of special relativity and the Einstein Equivalence Principle (EEP) holds. EEP is at the basis of all metric theories of gravity, that is to those theories that, like GR, states that gravitation is a curved spacetime phenomenon. Consequently, in GR, the local inertial frames are not simply determined by the distribution of matter, but they are influenced by the distribution-and-flow of mass-currents and mass-energy in the Universe. Therefore, the gravitomagnetic field associated with these phenomena exerts a force (in Newtonian words) on moving bodies and changes the spacetime structure produced by the gravitoeletric field by participating in generating the curvature of spacetime.

In 1918, Thirring [43] and Lense & Thirring [44] computed the drag exerted by the spin of a central body, like a planet, on an orbiting moon3. They computed the secular rates of the longitude of pericenter and of the right ascension of the ascending node for some of the moons of the Solar System: indeed, the orbital plane of the moon can be seen as a sort of enormous gyroscope dragged by the planet’s gravitomagnetic field4. Therefore, the intrinsic angular momentum of the planet, which plays the role of mass-current, according to GR curves spacetime and induces an effect on the orbit of a satellite which does not have a counterpart in the Newtonian theory.

Consequently, an accurate and reliable measurement of the gravitomagnetic field of the Earth is not only important per se, as a further and robust test of the GR predictions in the WFSM limit. There are at least three main issues that, for their importance, require a much more precise and accurate measurement of gravitomagnetism, even in weak-field conditions as the one that is considered in this paper:

• intrinsic gravitomagnetism;
• strong fields and compact objects; and,
• Mach’s principle.

These are briefly discussed here.

2.1. Intrinsic Gravitomagnetism

Intrinsic gravitomagnetism describes the spacetime curvature effects related with mass-currents, and not simply the corresponding effects that depend on the motion of a body on a static background field. This, in a nutshell, describes the main difference between the de Sitter [24] precession, which depend on the mass, and the Lense-Thirring precession, which depend on the angular momentum [5,60].

2.2. Strong Fields and Compact Objects

Two issues of modern astrophysics that require, for their explanation, the contribution of GR, are the mechanisms at the basis of the accretion of matter around compact objects and those of acceleration of charged particles in sources of various nature and size, with the simultaneous emission of powerful radio jets. Currently, the nature of the accretion mechanisms as well as that of the acceleration of particles and of the emission of relativistic jets is largely unknown under many aspects.

The classical physical processes presently considered as driving the interaction of the flow of matter with the central compact object, still require the correct inclusion, in the models developed to describe the details of the interaction, of the relativistic mechanisms that are produced by the gravitomagnetic field of the rotating compact object. These relativistic effects produce a dragging of the accretion disk and the alignment of the radio jets with the spin of the compact object [10,11,12,13,14,61]. Numerical relativity plays a major role in this context but, to our knowledge, papers that focus explicitly on the contribution due to the gravitomagnetic effects are not available in the literature.

2.3. Mach’s Principle

Mach’s ideas on the origin of inertia and on the relativity of rotation [62], as in the discussion regarding the well-known bucket experiment by Newton, strongly influenced Einstein in the development of GR. Mach was very precise and clear on these aspects: there is no absolute space and the Newtonian concept of acceleration has only meaning with respect to the very distant stars. Consequently, Mach is implicitly referring to a cosmological context. Again, in a nutshell, the sentence “inertia here arises from mass-energy there” [5] summarizes very well the meaning of gravitomagnetism and makes us interpret the Lense-Thirring effect (i.e., the dragging of inertial frames) as a weak GR formulation of the ideas of Mach, which is of what Einstein called Mach’s Principle [1,63,64]. For further details on the relationship between the Mach’s Principle and the Lense-Thirring effect, we refer to [65,66]. Concerning the cosmological implications of Mach’s Principle, these are complex, still debated and theory-dependent. We refer to the literature for further details, see e.g., [5,57,67,68,69].

3. On Past and Recent Measurements of the Lense-Thirring Effect

To this date, we lack an observational direct evidence of the gravitomagnetic effects in the strong field regime. On the other hand, in the WFSM limit of GR, several different works have successfully provided a measurement of the gravitomagnetic field produced by the Earth [22,23,30,33,34,36,37]. In fact, as we have already highlighted, the precession of the orbital plane of an Earth-bound satellite that is caused by the so-called Lense-Thirring effect [43,44] represents one of the most peculiar predictions of Einstein’s theory of GR related to gravitomagnetism. The orbital plane of the satellite behaves, in principle, as an external inertial frame not bound to the Earth and dragged by its rotation; such dragging is observable if the classical disturbances due to the main gravitational and non-gravitational perturbations are properly modeled during the precise orbit determination (POD), or can be removed a posteriori, or considered to be negligible on the time span of the orbital analysis.

The dragging effect manifests itself in a secular shift of the right-ascension of the ascending node (RAAN) $\Omega$ and of the argument of pericenter $\omega$ of the considered satellites [44]:

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

(2)

$${\dot{\omega }}_{\mathrm{LT}}=-\mu \frac{6G{J}_{\oplus }}{c^2{a}^3{(1-e^2)}^{3/2}}cosi=-3\mu K_{\mathrm{LT}}cosi,$$

(3)

where a, e, and i represent, respectively, the semi-major axis, the eccentricity and the inclination of the orbit, and $J_{\oplus }$ represents the angular momentum of the Earth. The dimensionless coefficient $\mu$ represents the Lense-Thirring effect parameter, with $\mu =1$ if GR is the correct theory of gravitation, and $\mu =0$ in Newtonian physics. Unfortunately, the measurements performed before 2018 had however a number of systematic errors that limited their accuracy well above 1%.

In the case of GP-B, the direct effect on an orbiting gyroscope has provided a measurement of frame-dragging with an accuracy of about 19% [22,23], after an intense and extended analysis effort needed to remove unexpected disturbances due to electric patch effects. This prevented the team to reach the original goal of 1%. Conversely, in the case of the measurements performed so far with the LAGEOS satellites [33,34,35], and also in the case of the first measurements that have included the LARES satellite [36,37], the main limitation in the accuracy of the measurement was provided by the uncertainties in the knowledge (and modeling) of the even zonal harmonic coefficients of the Earth’s gravitational potential. The accuracy of the previously cited measurements of the LT effect with laser-ranged satellites, was estimated in the 5–10% range, but gave rise to a vibrant debate in the literature [34,47,48,49,51,52].

One of LARASE’s key, albeit intermediate, objectives was to improve the modeling of the orbits of the three satellites, with the ultimate goal of refining the results that can be achieved in testing the gravitational interaction [37,40,46]. These improvements concern the modeling of both gravitational and non-gravitational perturbations [40,46,70,71,72,73,74,75,76]. In this regard, at the end of 2017, within the LARASE experiment, we started a series of analyses that aimed to improve the modeling of the gravitational field of the Earth and, in particular, the contribution of the even zonal harmonics coefficients ${\overline{C}}_{\ell ,0}$ (i.e., those with degree $\ell =\mathrm{even}$ and order $m=0$) and of their uncertainties. In fact, these coefficients are responsible for a secular effect on $\Omega$ and $\omega$, as that produced by the Earth’s gravitomagnetism, see Section 4. The analysis, at first tested on the first three even zonal harmonics [38,40,77], was later extended to the even zonal harmonics up to degree $\ell =30$ [41,72]. In these new analyses, we exploited the monthly solutions for the Earth’s gravitational field coefficients from the GRACE mission [78]. This allowed for us to considerably improve the accuracy of the measurement of the Lense-Thirring effect, strongly reducing the impact of the uncertainties of the even zonal harmonics in the final measurement of the relativistic precession. These aspects will be described in detail in Section 4. We obtained an error of about 2.3% in [40], and then reduced to about 1.6% in [41], both computed in a root-sum-square fashion. A similar result, about 2% in accuracy, has been obtained by [42].

Table 1. Mean orbital elements of the satellites.

Element	Unit	Symbol	LAGEOS	LAGEOS II	LARES
semi-major axis	[km]	a	12,270.00	12,162.07	7820.31
eccentricity		e	0.0044	0.0138	0.0012
inclination	[deg]	i	109.84	52.66	69.49

and

Table 2. Rate, in millisecond of arc per year (mas/yr), for the secular Lense–Thirring precession on the right-ascension of the ascending node and on the argument of pericenter of LAGEOS, LAGEOS II, and LARES satellites.

Rate in the Element	LAGEOS	LAGEOS II	LARES
${\dot{\Omega }}_{\mathrm{LT}}$	+30.67	+31.50	+118.48
${\dot{\omega }}_{\mathrm{LT}}$	+31.23	−57.31	−334.68

show, respectively, the mean orbital elements of the satellites and the expected relativistic Lense-Thirring precession on their orbits. These mean Keplerian elements have been determined by a dedicated analysis of the orbits of the satellites [46].

In the following, we will focus on the RAAN $\Omega$ as observable for the measurement of the relativistic precession. The pericenter $\omega$ of the satellites is influenced, especially for the two LAGEOS, by the effects of thermal thrust. A significant activity for the development of new models for the thermal forces—that allow us to recover the pericenter of the satellites and, in particular, that of LAGEOS II—has been carried out in recent years [73,75,79,80]. Preliminary results are encouraging and we plan to adopt them in future analyses.

4. On the Accuracy of the Even Zonal Harmonics

The gravitational potential V of the Earth is usually expanded in series of spherical harmonics in order to account for the non-uniform distribution of the mass of our planet. Concerning the contribution of the even zonal harmonics coefficients, ${\overline{C}}_{\ell ,0}$, to the geopotential, this may be written as [81,82]:

$$V_{\ell ,0}=-\frac{G{M}_{\oplus }}{r}\sum _{\ell =2}^{\infty }{\left(\frac{R_{\oplus }}{r}\right)}^{\ell }{P}_{\ell ,0}(\sin \phi ){\overline{C}}_{\ell ,0},$$

(4)

where $\phi$ is the geocentric latitude, $R_{\oplus }$ and $M_{\oplus }$ are, respectively, the mean equatorial radius and the mass of the Earth, r the geocentric distance, while $P_{\ell ,0}(\sin \phi )$ are the Legendre polynomials. The deviation from the spherical symmetry for the mass distribution of the Earth, described by these harmonics, is responsible for a classical precession of the RAAN of the satellite orbit, and with much larger amplitudes than Lense-Thirring precession. For instance, the secular effect in the rate of the RAAN due to the first two even zonal harmonics is [81]:

$${\dot{\Omega }}_{\mathrm{class}}=\frac{3}{2}n{\left(\frac{R_{\oplus }}{a}\right)}^2\frac{\cos i}{{\left(1-e^2\right)}^2}\left\{\sqrt{5}{\overline{C}}_{2,0}+\frac{15}{8}{\left(\frac{R_{\oplus }}{a}\right)}^3(7{\sin }^{2}i-4)\frac{(1+\frac{3}{2}{e}^2)}{{(1-e^2)}^2}{\overline{C}}_{4,0}\right\}=K_{2,0}{\overline{C}}_{2,0}+K_{4,0}{\overline{C}}_{4,0},$$

(5)

where ${\overline{C}}_{2,0}$ and ${\overline{C}}_{4,0}$ are the normalized Stokes coefficients of the quadrupole and octupole moments of the Earth and $n=\sqrt{G{M}_{\oplus }/a^3}$ is the satellite mean motion. The order-of-magnitude of this classical precession is about +126 deg/yr for LAGEOS, −231 deg/yr for LAGEOS II, and about −624 deg/yr in the case of LARES. Thus, in order to avoid a large systematic error in the measurement of the Lense-Thirring precession, we started a careful analysis of the Earth’s multipolar moments and of their dependence on time [83,84,85,86]5. The first steps in this direction were taken in [40] in the case of the quadrupole coefficient that is responsible for the largest effect; in [41] we extended those actions to multipole coefficients of higher order, up to degree $\ell =30$6. As mentioned above, in this work, we take up, extend, and complete these considerations.

Part of these improvements were made possible by the new solutions of the Earth’s gravitational field provided, especially, by the GRACE space mission. The twin GRACE satellites, launched on March 2002 by NASA and DLR, have provided a much improved determination of the Earth’s gravitational field, both in its static and dynamical components [78,87,88,89], with respect to previous results based on multi-satellite data. Moreover, the correlations among various coefficients, in particular the even zonal ones, have been reduced with respect to previous models, such as EGM96 [90].

In our analyses, we considered the monthly solutions7 for these coefficients provided by three independent analysis centers: CSR, GFZ, and JPL8. We then fitted these coefficients with a linear trend on the time span of our analysis, about 6.5 years, modeling each even zonal harmonic between $\ell =2$ and $\ell =20$ with the expression9:

$${\overline{C}}_{\ell ,0}(t)={\overline{C}}_{\ell ,0}(t_0)+{\dot{\overline{C}}}_{\ell ,0}(t-t_0),$$

(6)

where $t_0$ represents a reference epoch. In the Lense-Thirring effect measurements prior to 2019, only the quadrupole and octupole coefficients where modeled, by means of a linear trend [91]10. The linear trend “follows”, on average, the more complex time dependence of these coefficients, which are also characterized by the mentioned annual and semi-annual periodicities, plus other subtler periodic effects [85,92]. This approach strongly reduces the discrepancies between the effective time behavior of these coefficients and the models that use static values for the same coefficients.

It is important to stress that the static models provide values for the coefficients of the Earth’s gravity field that are averages on the time span of GRACE data. However, in different subintervals of the entire period of GRACE data, as in the case of our analysis, which starts after the launch of the LARES satellite, these average values may be quite different compared with their effective time behavior provided by GRACE monthly solutions.

In Figure 1, we exemplified, for the case of the coefficient ${\overline{C}}_{20,0}$, the procedure that we adopted to model the even zonal harmonics up to $\ell =30$. First, we compared the time behavior of each coefficient obtained from GRACE data with the corresponding static value provided in GGM05S. This is our main reference model for the gravitational field of the Earth, see Section 5. In general, the value provided by GGM05S represents (as highlighted) an average—in the time interval over which the gravitational field has been determined—of the time-dependent values obtained from the GRACE mission. However, this overall average value is, in general, far from the average value valid in the time span of our analysis, which starts after the launch of LARES, in February 2012. Moreover, as we can see from the figure, a simple average for the coefficient over the time span of the analysis, is not generally able to capture the time evolution of the coefficient, which is better described by a linear fit to the GRACE-derived coefficients.

Hereinafter, as usual practice in this kind of analysis, we estimated, together with the relativistic precession, the corrections to the quadrupole ${\overline{C}}_{2,0}$ and octupole ${\overline{C}}_{4,0}$ coefficients, (by solving Equations (7) below); these are the coefficients that most contribute to the systematic errors, if not properly accounted for in the POD of the satellites; practically, the value for $\mu$ is independent of the errors related to the imperfect knowledge of these coefficients (see Section 6). We also explicitly computed the correlation of these coefficients with the relativistic parameter $\mu$.

In the next section, we apply our considerations and analyses to several models for the background gravitational field of the Earth, in order to extract a reliable and robust measurement of the Lense-Thirring effect, as in [41]. This allows us, once the first 10 even zonal harmonics for each field were set in agreement with GRACE monthly solutions, to highlight the variability of the systematic errors among the considered fields. Consequently, this variability is due to even zonal harmonics of higher order, as well as to the non-zonal ones (tesseral and sectorial), respectively with $m\ne \ell$ and with $m=\ell$. Below, we first describe the characteristics of the new analyses that we have performed, and then we present the results that were obtained for the measurement of the relativistic precession of the orbital plane of terrestrial satellites due to the gravitomagnetic field.

5. The New Analyses

We analyzed the tracking data of LAGEOS, LAGEOS II, and LARES on a time span of 2359 days (about 6.5 years), starting from 6 April, 2012 (MJD 56023). The GEODYN II [93,94] software was used for the data reduction of the SLR observations of the three satellites in the normal point format [95]. These PODs have been performed by adjusting the state-vector (i.e., position and velocity) of each satellite on an arc length of seven days, to best fit the SLR data, for an overall number of 337 arcs not causally connected. The relativistic Lense-Thirring acceleration was not included in the dynamical model of GEODYN [96]. Also the subtle thermal thrust effects due to Earth’s infrared radiation (Earth-Yarkovsky effect) and to the Sun visible radiation, modulated by the eclipses (solar Yarkovsky-Schach effect) [97], were not included in the dynamical model, because, unfortunately, the routines included in GEODYN for their modelling can only handle a fast rotation regime of the satellites, which is not longer applicable for the two LAGEOS [75,98], averaging out many relevant thermal effects11.

In our new analyses we used the GGM05S solution [85,99] as reference model for the Earth’s gravitational field12. Other models considered in this work are: EIGEN-GRACE02S (2004) [89], ITU_GRACE16 (2016) [101], and Tonji-Grace02s (2017) [102,103]. All of these models were obtained from GRACE data13. For the considered gravitational models, the first 10 even zonal harmonics were modeled with the time behavior described above, while, for the remaining harmonics, the values of each single solution were used14. For solid tides we used the Colombo’s model [104], while ocean tides were described by the GOT99.2 model [105]. Furthermore, we adhered to the ILRS recommendations (those relevant for the time span of our analyses) regarding the on-ground tracking stations [106]. These recommendations regard mainly an a posteriori correction of the SLR data to avoid (or reduce) possible systematic errors that were related to several station-dependent features such as calibration issues, a possible change of the laser frequency, or of the coordinates of the station itself (within the ITRF) because of an earthquake, for instance. Among these corrections, the station biases play a special role, in particular the so-called range bias, which is correlated with the satellite center-of-mass correction [107]. We remark that the considered time span, about 6.5 years, is close to twice the period of the node of LAGEOS (1052 days), four times the period of the node of LAGEOS II (570 days) and 11 times that of LARES (211 days). This allows us to reduce the impact on the relativistic measurement of all the unmodeled (or poorly modeled) effects related to the non-gravitational and gravitational perturbations that are characterized by these periodicities. For instance, this is the case for some thermal thrust perturbations [97], or of the long-period ocean tide $K_1$ [71]. In the following sections, the most significant results obtained from the analysis of the residuals in the RAAN of the three satellites are described.

6. The Measurement of $\mathbf{\mu }$

Our observable for the measurement of the Lense-Thirring effect is represented by the RAAN $\Omega$ of the three satellites, as anticipated at the end of Section 3. Furthermore, as discussed with some details in Section 4, since the main effects that limit the measurement of the relativistic precession are the uncertainties in the knowledge of even zonal harmonics coefficients ${\overline{C}}_{\ell ,0}$, in particular of the first two coefficients, we have solved (separately for each arc) a system of three equations in the three unknowns $\delta {\overline{C}}_{2,0}$, $\delta {\overline{C}}_{4,0}$, and $\mu$ (following the original proposal of [108]15):

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{K}_{2,0}^{L1}\delta {\overline{C}}_{2,0}+K_{4,0}^{L1}\delta {\overline{C}}_{4,0}+K_{LT}^{L1}\mu =\delta {\dot{\Omega }}_{res}^{L1}\\ K_{2,0}^{L2}\delta {\overline{C}}_{2,0}+K_{4,0}^{L2}\delta {\overline{C}}_{4,0}+K_{LT}^{L2}\mu =\delta {\dot{\Omega }}_{res}^{L2}\\ K_{2,0}^{LR}\delta {\overline{C}}_{2,0}+K_{4,0}^{LR}\delta {\overline{C}}_{4,0}+K_{LT}^{LR}\mu =\delta {\dot{\Omega }}_{res}^{LR},\end{array}\right.\end{array}$$

(7)

In these equations16, the constant terms $\delta {\dot{\Omega }}_{res}^{L1}$, $\delta {\dot{\Omega }}_{res}^{L2}$ and $\delta {\dot{\Omega }}_{res}^{LR}$ represent the residuals in the rate of the RAAN of each satellite determined, after the POD of their orbit, with the method described in [110]. In the left hand side of each equation, the coefficients $K_{2,0}^{Lj}$ and $K_{4,0}^{Lj}$ (with $Lj=(L1)\mathrm{LAGEOS},(L2)\mathrm{LAGEOS}\mathrm{II},(LR)\mathrm{LARES}$) are function of the orbital parameters a, e and i of each satellite, as explicitly shown in Equation (5)17. The three coefficients $K_{LT}^{Lj}$ represent the relativistic Lense-Thirring precession for the three satellites, see Equation (2) and Table 2. Finally, concerning the three unknowns of the system, the quantities $\delta {\overline{C}}_{\ell ,0}$ (with $\ell =2$ and $\ell =4$) represent the mismodeling of the even zonal harmonics of the background model of the Earth’s gravitational field, while $\mu$ represents the dimensionless Lense-Thirring parameter introduced in Equation (2)18. Therefore, with this solution, we can remove, from the estimate of $\mu$, the uncertainties that are related to the knowledge of the first two even zonal harmonics19. Figure 2 shows the results for the corrections to the quadrupole and octupole coefficients, while, in Figure 3, we show the results for $\mu$, evaluated, for each arc, as a solution of the system of Equation (7).

Figure 4, plots $\mu$ as a cumulative sum of the values shown in Figure 3. The set shown was computed using the GGM05S Earth model. The other models considered provide similar results: see below for a comparison.

To extract the Lense-Thirring parameter, we compute the mean value in the case of Figure 3, and the slope of a best fitting straight line in the case of Figure 4. It can be seen that periodic effects perturb, especially in the case of Figure 3, the measurement of the relativistic secular effect embedded in the time series of the residuals. However, we underline that the unmodeled periodic effects that characterize the residuals in Figure 3 behave—in relation of their peculiar nature that alternatively increases or decreases the unmodeled precession on various arcs—as a Gaussian (with a Skewness of $-8.4\times {10}^{-3}$ and Kurtosis of +3.097) superimposed to the relativistic precession.

The periodic effects that characterize the residuals in $\mu$ are due to both gravitational and non-gravitational perturbations20. However, the separation of the different effects is not simple, being influenced by the time dependence of the even zonal harmonics and, above all, by the variation of the spin of the satellites. Among the main spectral lines, we have the periodicities of the eclipses of the satellites $2(\dot{\Omega }-\dot{\lambda })$21, especially that of LAGEOS (about 280 days), and the periodicities, already mentioned, related to the time variation of the RAAN of the satellites $\dot{\Omega }$.

Because the cumulative sum acts as a low-pass filter, the residuals shown in Figure 4 are less influenced by the highest frequencies present in Figure 3, and allow for a cleaner, and consequently more precise, measurement of the Lense-Thirring effect. In

Table 3. Correlations among the three estimated quantities evaluated as solutions of Equation (7) for the GGM05S model.

	$\mathit{\delta }{\overline{\mathit{C}}}_{2,0}$	$\mathit{\delta }{\overline{\mathit{C}}}_{4,0}$	$\mathit{\mu }$
$\delta {\overline{C}}_{2,0}$	+1.000	+0.082	+0.071
$\delta {\overline{C}}_{4,0}$		+1.000	−0.179
$\mu$			+1.000

, we show, still in the case of the GGM05S model, the correlation matrix for the three solutions of the system of Equation (7).

The correlation is particularly small between $\mu$ and the correction $\delta {\overline{C}}_{2,0}$, as well as between the two corrections, $\delta {\overline{C}}_{2,0}$ and $\delta {\overline{C}}_{4,0}$. Conversely, the correlation is moderately larger between the Lense-Thirring parameter and the correction $\delta {\overline{C}}_{4,0}$. The relative low degree of correlation among the relativistic parameter $\mu$ with the corrections to the first two even zonal harmonics can be represented in the colorbar plot that is shown in Figure 5, where the values for $\mu$ are distributed in relation to those assumed by $\delta {\overline{C}}_{2,0}$ and $\delta {\overline{C}}_{4,0}$.

This analysis provides better (smaller) correlations with respect to a previous measurement [40], where the Lense-Thirring effect was estimated together with the corrections to the quadrupole and hexapole, $\delta {\overline{C}}_{6,0}$, coefficients22. Similar results are found when using all other gravitational fields. This is due to the common use we made, among the considered gravitational fields, of the first 10 even zonal harmonics estimated from GRACE monthly solutions. This is well represented in Figure 6, where the comparison among the results that were obtained for the Lense-Thirring parameter for the various models considered is shown.

That is to say, when the even zonal harmonics are well modeled in the POD of the different satellites, in the sense of considering for their dependence on time a model as close as possible to the temporal trend of these coefficients as measured by the monthly solutions of GRACE, the differences between the different gravitational field models tend to considerably decrease as regards the measurement of the Lense-Thirring precession.

In

Table 4. Comparison of the results for the measurement of the relativistic Lense-Thirring precession among different models for the background gravitational field of the Earth, computed with the method of cumulative residuals. Errors appear to be equal due to rounding to the most significant digit. These results for the measured parameter $\mu$ are compatible (within their errors) with the prediction of GR for the Lense-Thirring precession: ${\mu }_{GR}=1$.

Model	$\mathit{\mu }\pm \mathit{\delta }\mathit{\mu }$
GGM05S	$1.0053\pm 0.0074$
EIGEN-GRACE02S	$1.0002\pm 0.0074$
ITU_GRACE16	$0.9996\pm 0.0074$
Tonji-Grace02s	$1.0008\pm 0.0074$

, we show the resulting value for the Lense-Thirring parameter $\mu$, estimated from the slope of the linear fit to the cumulative residuals, for each of the gravitational models considered in our analyses.

We notice that, in all these analyses, the agreement with the GR prediction is at a fraction of percent. Additionally, the precision $\delta \mu$ (at a 95% confidence level) is at a fraction of percent of the relativistic precession. As already highlighted, the variability of the results, actually very small, is mainly related to possible systematic differences among the different solutions for the gravitational field, due to even zonal harmonic coefficients with degree $\ell >20$. Assuming, conservatively, for this measurement of the Lense-Thirring effect, the mean of the measured values of $\mu$ reported in Table 4, we obtain:

$${\mu }_{\mathrm{meas}}-1=(1.5\pm 7.4)\times {10}^{-3}\pm 16\times {10}^{-3},$$

(8)

where the first terms within parentheses represent the measurement and its formal error, while the last term represents our estimate for the error budget, 1.6%. This error budget derives from a root-sum-square (RSS) of the main systematic effects that are related to the gravitational and non-gravitational perturbations that act on the orbits of the satellites. It accounts for the errors related to the static field (≃1.0%), to ocean tides (≲0.6%) 23, to other periodic effects (≃1.0%) and to the error related to the knowledge of the de Sitter precession (≃0.3%). Conversely, if we pessimistically consider the sum of the absolute values of these errors, the error budget increases to about 2.9%.

7. A statistical Approach to the Measurement of $\mathbf{\mu }$

Hereinafter, we consider two aspects introduced in the previous sections: the time interval of the analysis and the unmodeled periodic effects. In Section 5, we have pointed out that the time interval of our analysis has the advantage to be a integer multiple of the periods of the RAAN of the three satellites. In this case, as discussed in Section 6, the mismodeling or poorly modeling (or even the non-modeling) of some perturbation linked with these main periods has no impact in the slope of the best fitting straight line through the integrated residuals of the relativistic parameter $\mu$. Consequently, the gravitational and non-gravitational perturbations that are characterized by these three main periods cannot be considered dangerous neither for the measurement of the relativistic effect nor for their impact in the final error budget of the measurement, i.e., for the estimate of the systematic errors. We underlined in Section 6 that the unmodeled periodic effects in the residuals of Figure 3 represent a sort of noise, with Gaussian distribution, masking the relativistic effect that represents the common (constant) signal along each arc. Of course, this periodic signal represents, in the first place, what is not modeled and which we must therefore find in the residuals. The choice of the total time interval of the analysis was made precisely to exploit this correlation, so that the main effects of some long-term perturbations mediate at zero, as explained, both in the case of the analysis of the residuals of Figure 3 and in the case of the analysis of the integrated residuals of Figure 4.

After clearing these aspects, we can estimate the impact of the periodic terms in the evaluation of the relativistic effect $\mu$ (the slope of the best fitting straight line of the integrated residuals of Figure 4 in the following way: we perform a series of random permutations of the unintegrated residuals of Figure 3, in order to decorrelate the relativistic effect present in them from the periodic effects that characterize the temporal sequence considered. For each permutation, we determine the new slope from a linear fit over the entire number of arcs. Finally, we check how the new solutions are distributed. Of course, for the Central Limit Theorem, we expect that these results are distributed very close to a Normal (i.e., Gaussian) distribution24. In Figure 7, we show the results obtained with the GGM05S model after 50,000 random permutations.

For the relativistic parameter we obtained:

$${\overline{\mu }}_{\mathrm{meas}}-1=(-1.11\pm 2.94)\times {10}^{-2},$$

(9)

where ${\overline{\mu }}_{\mathrm{meas}}=0.9889$ represents the mean of the Gaussian and the error represents the standard deviation ${\sigma }_{\mu }$ of the Gaussian with a 95% confidence level. This result has a discrepancy of about 1% with respect to the result that was obtained in previous Section 6; consequently, this value confirms a 1% contribution of the unmodeled periodic terms in the overall error budget of the measurement25.

Finally, in Figure 8 and

Table 5. Results for the Lense-Thirring parameter after the Gaussian fit described in the text and based on 50,000 random permutations of the residuals obtained for each 7-day arc from the POD described in Section 6. The error ${\sigma }_{\mu }$ represents the standard deviation of each Gaussian fit at a 95% confidence level.

Model	$\overline{\mathit{\mu }}\pm {\mathit{\sigma }}_{\mathit{\mu }}$
GGM05S	$0.9889\pm 0.0294$
EIGEN-GRACE02S	$0.9778\pm 0.0292$
ITU_GRACE16	$0.9772\pm 0.0291$
Tonji-Grace02s	$0.9788\pm 0.0292$

, the results for the Lense-Thirring parameter after the 50,000 random permutations are shown for the other gravitational fields considered in the present work. On average, the discrepancy in measuring $\mu$ with this statistical approach is approximately 1.9% when compared with the average result that is presented in the previous Section 6.

8. Conclusions

We presented the results of a precise and accurate measurement of the gravitomagnetic field of the Earth through the measurement of the Lense-Thirring effect on the combined orbits of the geodetic satellites LAGEOS, LAGEOS II, and LARES. These kinds of measurements raised criticisms [34,47,48,49,51,52] in the past regarding the estimate of the systematic sources of error due to the background gravitational field, which is responsible for a classical orbital precession that completely mask, if not adequately considered, the precession envisaged by Einstein’s theory of GR. We believe that the analysis reported here well overcomes these critical issues. Indeed, several new improvements have been introduced. Among these, the most significant is the modeling—in the POD of the satellites—of the even zonal harmonics of the Earth’s gravitational field through linear trends that fit the monthly solutions of the corresponding GRACE coefficients in the time interval of our analyses, about 6.5 years. We also estimated, together with the relativistic parameter $\mu$, the corrections to the quadrupole ${\overline{C}}_{2,0}$ and octupole ${\overline{C}}_{4,0}$ coefficients and the correlation values among these three quantities.

These improvements have allowed us to perform the measurement by means of a simple linear fit to the orbit residuals. Indeed, whenever some of the coefficients of the gravitational field are not correctly modeled in the POD, the simple linear fit preserves the analysis from the risk of forcing the result to the desired value, using an excessive number of fitting parameters: because these could absorb some of the periodic effects beyond our ability of understanding their origin. Our final result is:

$${\mu }_{\mathrm{meas}}-1=(1.5\pm 7.4)\times {10}^{-3}\pm 16\times {10}^{-3},$$

(10)

where the formal error from the fit (0.74%) has been conservatively given with a 95% confidence level. Usually, in the literature of the Lense-Thirring effect measurements, this formal error is provided at a 68% confidence level; consequently, for comparison with other recent works [40,42], our error should be rescaled to 0.38%. This result is compatible (within the errors) with the prediction of GR for the relativistic precession produced by the gravitomagnetic field of the Earth, i.e., with ${\mu }_{GR}=1$. The discrepancy with GR is only 0.15%, and the accuracy, calculated in RSS fashion, is 1.6%, robustly improved with respect to previous estimates.

Concerning the contribution of the unmodeled periodic effects to the error budget, the realization of a single measurement over an interval of 6.5 years has turned into 50,000 different realizations by means of the random permutation of the residuals obtained for each of the 337 orbital arcs of seven days that make up the entire period of our analysis. This allowed us to obtain a further, and independent, estimate of the errors that are related to the unmodeled periodic terms, of the order of 1%. It is clear that future work in this context will require a deeper insight and comprehension of the nature of these periodic effects in order to disentangle the contribution of the gravitational perturbations from the subtler, non-gravitational ones.