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

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 (' 7.4× 10−3, at a 95% confidence level) and accuracy (' 16× 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.


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/rc 2 1, J/Mrc 1, and v/c 1, where G, M, and J represent, respectively, the gravitational constant, 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 µ in two cases: from the residuals of µ 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.

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: where S represents the spin of the primary. This particular precession is known as "dragging of gyroscopes" or "dragging of inertial frames", as Einstein named it 2 . 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 moon 3 . 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 field 4 . 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. 2 In a laboratory the gyroscopes are used to define the axes of a locally inertial frame. 3 See [57,58] for the description of the direct contribution of Einstein in the derivation of the so-called Lense-Thirring effect. 4 In the papers by Thirring and by Lense & Thirring, as well as in many other subsequent papers on the argument, the angular momentum of the central body is assumed aligned with its z-axis. For the general case see [59].
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; 2.
strong fields and compact objects; and, 3.
These are briefly discussed here.

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].

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.

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].

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) Ω and of the argument of pericenter ω of the considered satellites [44]: where a, e, and i represent, respectively, the semi-major axis, the eccentricity and the inclination of the orbit, and J ⊕ represents the angular momentum of the Earth. The dimensionless coefficient µ represents the Lense-Thirring effect parameter, with µ = 1 if GR is the correct theory of gravitation, and µ = 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 coefficientsC ,0 (i.e., those with degree = even and order m = 0) and of their uncertainties. In fact, these coefficients are responsible for a secular effect on Ω and ω, 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 = 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]. Tables 1 and 2 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].  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 LAREṠ
In the following, we will focus on the RAAN Ω as observable for the measurement of the relativistic precession. The pericenter ω 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.

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,C ,0 , to the geopotential, this may be written as [81,82]: where ϕ is the geocentric latitude, R ⊕ and M ⊕ are, respectively, the mean equatorial radius and the mass of the Earth, r the geocentric distance, while P ,0 (sin ϕ) 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]: whereC 2,0 andC 4,0 are the normalized Stokes coefficients of the quadrupole and octupole moments of the Earth and n = GM ⊕ /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 5 These are characterized by several periodic effects with, mainly, annual and inter-annual periodicities. multipole coefficients of higher order, up to degree = 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 solutions 7 for these coefficients provided by three independent analysis centers: CSR, GFZ, and JPL 8 . 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 = 2 and = 20 with the expression 9 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 coefficientC 20,0 , the procedure that we adopted to model the even zonal harmonics up to = 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 6 In [40], we linearly fitted the quadrupole coefficient that was obtained from GRACE monthly solutions and we used this fitted value in the data reduction of the satellites orbit. Furthermore, in that paper and for the GGM05S model, with regard to the Lense-Thirring effect measurement, we have also analyzed and compared the error related to the knowledge of the octupole coefficient with respect to the hexapole one,C 6,0 . 7 International Centre for Global Earth Models (ICGEM): Gravity Field Solutions for dedicated Time Periods: Release 05 and Release 06, http://icgem.gfz-potsdam.de/series. In the case of Release 06, no significant differences are present for the first 10 even zonal harmonics as compared to Release 05. Some small differences are present for the harmonics with degree ≥ 22, but with no significant impact in the POD of the satellites. 8 Center for Space Research (CSR), University of Texas at Austin; GeoForschungsZentrum (GFZ), German Research Centre for Geosciences, Potsdam; Jet Propulsion Laboratory (JPL), California Institute of Technology, Pasadena. 9 The even zonal harmonics between degree = 22 and degree = 30 play a minor role, with no appreciable difference in the results. For the quadrupole coefficient, we also considered a more complete non-linear fit to the time behavior outlined by GRACE, with the inclusion of two periodic terms, one with a yearly frequency and one at twice this frequency. However, the final results have not shown a noticeable difference with respect to those obtained with a simpler linear fit. 10 Usually, these trends were those suggested by IERS Conventions [91], but were not always compatible with the results from GRACE monthly solutions. In the literature of the past measurements of the Lense-Thirring effect it has never been specified whether harmonics of degree higher than = 4 were modeled according to linear trends. Within the ILRS community, and since a few years, some harmonics are modeled up to degree 6 with secular trends, and very recently up to degree 12 in the case of the JCET Analysis Center. 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 quadrupoleC 2,0 and octupoleC 4,0 coefficients, (by solving Equation (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 µ 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 µ.
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 = and with m = . 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.

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 effects 11 .
In our new analyses we used the GGM05S solution [85,99] as reference model for the Earth's gravitational field 12 . 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 data 13 . 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 used 14 . 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.

The Measurement of µ
Our observable for the measurement of the Lense-Thirring effect is represented by the RAAN Ω 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 coefficientsC ,0 , in particular of the first two coefficients, we have solved (separately for each arc) a system of three equations in the three unknowns δC 2,0 , δC 4,0 , and µ (following the original proposal of [108] 15 ): 11 For the same reason, we could not include the pericenter in our analysis. Our ongoing effort to improve the model for the thermal thrust perturbations [40,73,75,97] will hopefully allow us to include, in forthcoming measurements, also this element in the analysis, especially for LAGEOS II. 12 This is the solution currently used by the International Laser Ranging Service (ILRS) to realize the International Terrestrial Reference Frame (ITRF2014) [100], i.e., the practical realization of the International Terrestrial Reference System (ITRS) [91]. 13 For the maximum degree of the field, we used an expansion up to degree and order 30 for the two LAGEOS and up to degree and order 90 for LARES, because of its smaller semi-major axis with respect to that of the two LAGEOS. In the case of LAGEOS satellites, the IERS Conventions [91] suggest at least an expansion up to degree and order 20. We have verified that for the two LAGEOS-considering an expansion up to degree and order 20, or 30, or 90-no significant difference in the POD of the satellites is produced. 14 From here on, when we refer generically to one of these solutions for the Earth's gravitational field, we will always imply the modified solution for the first 10 even zonal harmonics, as described in the text. 15 More precisely, the solution of these equations for the combined RAAN of the three satelltes was explicitly computed in [48] following the original proposal of [108]. This approach was proposed, for the first time, by [109] in a Sun-planet context.
In these equations 16 , the constant terms δΩ L1 res , δΩ L2 res and δΩ LR res 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 Lj 2,0 and K Lj 4,0 (with Lj = (L1)LAGEOS, (L2)LAGEOS II, (LR)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 Lj LT 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 δC ,0 (with = 2 and = 4) represent the mismodeling of the even zonal harmonics of the background model of the Earth's gravitational field, while µ represents the dimensionless Lense-Thirring parameter introduced in Equation (2) 18 . Therefore, with this solution, we can remove, from the estimate of µ, the uncertainties that are related to the knowledge of the first two even zonal harmonics 19 . Figure 2 shows the results for the corrections to the quadrupole and octupole coefficients, while, in Figure 3, we show the results for µ, evaluated, for each arc, as a solution of the system of Equation (7).   16 In this system of Equations we are neglecting the contributions from the mismodeling of the higher harmonics. Their contribution, however, is explicitly considered in the evaluation of the systematic errors, i.e., in the overall error budget of the measurement. 17 The orbital parameters are known with a very small relative uncertainty, such that the K coefficients can be considered, for our purposes, as error-free. 18 We directly provide the result for µ and not, as done in previous measurements of the Lense-Thirring effect, for the combined rate of the RAAN of the three satellites. This combination corresponds to a precession of 50.17 mas/yr. 19 We underline that the solution obtained for µ is independent of a possible imprint of the Lense-Thirring effect in the gravitational field coefficients. In fact, if this imprint were present, it would be mainly contained in the quadrupole coefficient, which is however eliminated by the combination that provides the Lense-Thirring parameter. However, from the tests we have performed in the past, there is no evidence of such a possibility.      Figure 3. The horizontal axis shows the number of arcs of the satellites POD. The slope of the linear fit (red line) provides a value of µ = 1.0053 ± 0.0074 (at a 95% confidence level), see Table 4. The impact of the unmodeled periodic effects is reduced in the case of the cumulative residuals. The GGM05S model was considered for the background gravitational field of the Earth.
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 × 10 −3 and Kurtosis of +3.097) superimposed to the relativistic precession.
The periodic effects that characterize the residuals in µ are due to both gravitational and non-gravitational perturbations 20 . 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(Ω −λ) 21 , especially that of LAGEOS (about 280 days), and the periodicities, already mentioned, related to the time variation of the RAAN of the satellitesΩ.
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, we show, still in the case of the GGM05S model, the correlation matrix for the three solutions of the system of Equation (7). Table 3. Correlations among the three estimated quantities evaluated as solutions of Equations (7) for the GGM05S model. The correlation is particularly small between µ and the correction δC 2,0 , as well as between the two corrections, δC 2,0 and δC 4,0 . Conversely, the correlation is moderately larger between the Lense-Thirring parameter and the correction δC 4,0 . The relative low degree of correlation among the relativistic parameter µ 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 µ are distributed in relation to those assumed by δC 2,0 and δ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, δC 6,0 , coefficients 22 . 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.   Table 4.
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, we show the resulting value for the Lense-Thirring parameter µ, estimated from the slope of the linear fit to the cumulative residuals, for each of the gravitational models considered in our analyses. 22 We reiterate that the Lense-Thirring parameter µ obtained from the solution of the system of Equation (7) is independent of any static and dynamic error that characterizes the first two even zonal harmonics. The parameter µ obtained is however influenced by any error related to the even zonal harmonics with degree ≥ 6. The correlations described in the text are strictly those between the relativistic parameter and the corrections to the first two even zonal harmonics with respect to their a priori values used in the POD. We notice that, in all these analyses, the agreement with the GR prediction is at a fraction of percent. Additionally, the precision δµ (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 > 20. Assuming, conservatively, for this measurement of the Lense-Thirring effect, the mean of the measured values of µ reported in Table 4, we obtain: 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%.

A statistical Approach to the Measurement of µ
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 µ. 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 µ (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) distribution 24 . In Figure 7, we show the results obtained with the GGM05S model after 50,000 random permutations.  For the relativistic parameter we obtained: whereμ meas = 0.9889 represents the mean of the Gaussian and the error represents the standard deviation σ µ 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 measurement 25 . Finally, in Figure 8 and Table 5, 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 µ with this statistical approach is approximately 1.9% when compared with the average result that is presented in the previous Section 6.  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 σ µ represents the standard deviation of each Gaussian fit at a 95% confidence level.

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 µ, the corrections to the quadrupoleC 2,0 and octupoleC 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: 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 µ 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.