A comment on"A test of general relativity using the LARES and LAGEOS satellites and a GRACE Earth gravity model", by I. Ciufolini et al

Recently, Ciufolini et al. reported on a test of the general relativistic gravitomagnetic Lense-Thirring effect by analyzing about 3.5 years of laser ranging data to the LAGEOS, LAGEOS II, LARES geodetic satellites orbiting the Earth. By using the GRACE-based GGM05S Earth's global gravity model and a linear combination of the nodes $\Omega$ of the three satellites designed to remove the impact of errors in the first two even zonal harmonic coefficients $J_2,~J_4$ of the multipolar expansion of the Newtonian part of the Earth's gravitational potential, they claimed an overall accuracy of $5\%$ for the Lense-Thirring caused node motion. We show that the scatter in the nominal values of the uncancelled even zonals of degree $\ell = 6,~8,~10$ from some of the most recent global gravity models does not yet allow to reach unambiguously and univocally the expected $\approx 1\%$ level, being large up to $\lesssim 15\%~(\ell=6),~6\%~(\ell=8),~36\%~(\ell=10)$ for some pairs of models.


LAGEOS
In this paper, we will quantitatively show that the claims in Ref. [5] are too optimistic despite the recent advances in modeling the Earth's gravity field.
Let us recall that the gravitomagnetic node precession of a test particle in geodesic motion around a central rotating primary is [10] Ω LT = 2GS where G and c are the Newtonian constant of gravitation and the speed of light in vacuum, respectively, S is the body's angular momentum, while a and e are the orbit's semimajor axis and eccentricity, respectively. The magnitudes of the Lense-Thirring precessions for the satellites of the LAGEOS family, along with their relevant orbital parameters, are listed in Table 1.
2 The systematic uncertainty due to the uncancelled even zonal harmonics One of the most critical issues of all the attempts planned or performed so far to detect the Lense-Thirring effect with Earth's artificial satellites is represented by the correct handling of the competing secular node precessionsΩ J ℓ induced by the even zonal harmonic coefficients 1 J ℓ = − √ 2ℓ + 1 C ℓ,0 , ℓ = 2, 4, 6, . . . of the multipolar expansion of the Newtonian terrestrial gravitational potential. Indeed, their nominal sizes are several orders of magnitude larger than the Lense-Thirring node precessions; moreover, their current level of modeling is not yet accurate enough to allow to use a single satellite's node to extract the relativistic signature. Thus, in the past years, some strategies to circumvent such an issue were devised. Contrary to what is claimed in Ref. [5], the combination presumably used in Ref. [5], not displayed there, was explicitly proposed by the present author 2 on the basis of a strategy proposed for the first time in Ref. [2]. It was suitably designed to cancel out, at least in principle, the node precessions due to the first two even zonals J 2 , J 4 , being fully impacted by all the other ones of higher degree ℓ ≥ 6. The hope is that, since the geopotential is modeled in the data reduction softwares used to process the satellites' data, the level of mismodeling in the uncancelled even zonals is accurate enough to allow for a reliable determination of the Lense-Thirring effect with a percent accuracy.
Actually, as summarized in Refs. [8,12] and references therein, there are sound doubts that such a strategy could allow to reach a ≈ 1% overall accuracy, as steadily claimed by Ciufolini et al. over the years. In this paper, we will enforce such concerns by showing that severe issues lurk even in the low-degree part of the spectrum of the geopotential.

Review of the method used so far to extract the Lense-Thirring signature
For the convenience of the reader, we will now review the aforementioned approach, which is a generalization of the one proposed for the first time in Ref. [2].
By assuming to use the nodes of, say, N different satellites, the following N linear combinations can be written down They involve the Lense-Thirring node precessions as predicted by General Relativity (see Eq. (1)), scaled by a multiplicative parameter µ LT , and errors in the computed secular node precessions due to errors in the first N − 1 even zonals J 2k , k = 1, 2, . . . N − 1, assumed as mismodeled through δJ 2k , k = 4 L. Iorio 1, 2, . . . N − 1. In the following, we will use the shorthanḋ for the partial derivative of the classical node precession with respect to the generic even zonal J ℓ of degree ℓ. Then, the N combinations of Eq. (2) are equated to the experimental residuals ∆Ω (i) , i = 1, 2, . . . N of each node of the N satellites considered. In principle, such residuals account for the purposely unmodelled Lense-Thirring effect, the mismodelling of the static and timevarying parts of the geopotential, and the non-gravitational forces. Thus, one gets If we look at the Lense-Thirring scaling parameter µ LT and the mismodeling in the even zonals δJ 2k , k = 1, 2, . . . N − 1 as unknowns, we can interpret Eq. (4) as an inhomogenous linear system of N algebraic equations in the N unknowns µ LT , δJ 2 , δJ 4 . . . δJ 2(N −1) whose coefficients arė while the constant terms are the N node residuals It turns out that, after some algebraic manipulations, the dimensionless Lense-Thirring scaling parameter, which is 1 in General Relativity, can be expressed as In Eq. (8), the combination of the N node residuals is, by construction, independent of the first N − 1 even zonals, being impacted by the other ones of degree ℓ > 2(N − 1) along with the non-gravitational perturbations and other possible orbital perturbations which cannot be reduced to the same formal expressions of the first N − 1 even zonal rates. Instead, combines the N gravitomagnetic node precessions as predicted by General Relativity. The dimensionless coefficients c j , j = 1, 2, . . . N − 1 in Eq. (9)-Eq. (10) depend only on some of the orbital parameters of the N satellites involved in such a way that, by construction, C ∆ = 0 if Eq. (9) is calculated by posing for any of the first N − 1 even zonals, independently of the value assumed for its uncertainty δJ ℓ . In our case [7], it is N = 3; the satellites employed are LAGEOS, LA-GEOS II, LARES. The resulting combination coefficients c 1 , c 2 yielding µ LT , worked out from [1] calculated for LAGEOS, LAGEOS II, LARES as in Table 1, turn out to be [7,8] In Eqs.
We point out that the numerical values of c 1 , c 2 in Eqs. (14)-(15) are quoted with nine decimal digits in order to assure a cancelation of J 2 accurate to better than 1% level. If the only concern were the evaluation of the impact of the errors in the even zonals of degree ℓ ≥ 6 discussed in Section 2.2, much less accurate values as c 1 ≈ 0.344, c 2 ≈ 0.0733 would be adequate; nonetheless, the first even zonal harmonic would nominally impact the combined Lense-Thirring trend at a 860% level. The authors of Ref. [5] did not quote any numerical values for c 1 , c 2 . In principle, the strategy reviewed above allows one to extract also δJ 2 or δJ 4 , provided the substitutionΩ The impact of the uncancelled even zonals of degree higher than ℓ = 4 and the need of using several Earth gravity models A major source of systematic uncertainty, to be reliably and realistically assessed, is represented by the impact of the mismodeling in the even zonals of degree ℓ ≥ 6 affecting the classical node precessions which are not cancelled by the adopted linear combination.
Ciufolini et al. [5] limited themselves to use only a single Earth's global gravity field solution, GGM05S [13], without any motivation for such a seemingly arbitrary choice. Moreover, in Ref. [5], it is claimed that, after arbitrarily tripling the released sigmas of the even zonal coefficients of such a model, the overall systematic uncertainty would be as little as 4%, without providing any details on either the computational strategy adopted or the maximum degree ℓ considered. Furthermore, in Ref. [5], it is also claimed, without any explicit and quantitative justifications, that different Earth gravity models would allegedly lead just to slightly different results.
Actually, in view of the nowadays ever increasing number of different Earth's global gravity models released by several independent institutions worldwide 3 , a fair practice should consist of taking into account an ensemble of various models and look into the resulting scatter of the a-priori evaluations of the overall systematic uncertainties with respect to the expected Lense-Thirring trend. Moreover, correctly assessing the realistic errors in the even zonals requires great care, following quantitative and unambiguous procedures such as those outlined in Ref. [14]; it should not be done in a hand-wavy fashion.
All such critical issues, ignored by the authors of Ref. [5], are still valid despite the most recent advances in the field of the determination of the Earth's global gravity field models. Below, we will demonstrate it by looking at some of the most recently released geopotential models.
In general, a straightforward way to analytically calculate the long-term orbital perturbations experienced by a test particle orbiting a primary of mass M due to a disturbing additional potential U pert with respect to the Newtonian monopole consists of taking its average over one orbital period P b and, then, inserting it in the Lagrange equations for the rates of change of the Keplerian orbital elements. In the case of the node, it is The temporal average of Eq. (17) has to be calculated onto the unperturbed Keplerian ellipse characterized by x = r (cos Ω cos u − cos I sin Ω sin u) , y = r (sin Ω cos u + cos I cos Ω sin u) , where ω, f, u .
= ω + f are the argument of pericenter, the true anomaly, and the argument of latitude, respectively.
In the case of an axially symmetric primary with equatorial radius R and symmetry axisk, the resulting perturbing potential of degree ℓ and order m = 0 is where J ℓ is the even zonal harmonic of degree ℓ, andr = r/r is the unit vector pointing from the primary to the test particle. If the primary's equatorial plane is assumed as reference {x, y} plane, as in the case of the Earth's artificial satellites, then In fact, also the uncancelled even zonals of degree as low as ℓ = 6, 8, 10 may have a non-negligible impact on the linear combination used because of their contribution to the nodal precession rate for LARES, as we will show in this Section. 8 L. Iorio Table 2 Estimated values C ℓ,0 and formal errors σ C ℓ,0 of the even zonal harmonics of degree ℓ = 6, 8, 10 according to the Earth's global gravity field solutions GOCO05S, ITU GRACE16, ITSG-Grace2014s, JYY GOCE04S retrievable on the Internet at http://icgem.gfz-potsdam.de/ICGEM/. Only the significant digits are displayed.
On the other hand, the remaining scatter in the estimated values of the corresponding Stokes geopotential coefficients C ℓ,0 from several global gravity models does not yet allow for an unequivocal and unambiguous mismodeling below the 1% level in the resulting node precessions of LARES.
This can be quickly shown by considering the simple differences ∆C ℓ,0 among the estimates for C ℓ,0 for various pairs of Earth's global gravity models as a naive measure of the realistic uncertainties in the even zonals of interest. In the following, we will demonstrate that also more refined calculation yield essentially the same conclusions. Table 2 displays the determined values C ℓ,0 along with their sigmas σ C ℓ,0 of some recent field solutions for the degrees ℓ = 6, 8, 10. Table 3 Absolute values of the differences ∆C 6,0 among the estimated values C 6,0 of the models of Table 2. They are very close to the uncertainties retrievable from the application of Eq. (A11) by Wagner and McAdoo [14] to the sigma of the model to be tested by contrasting it to a formally superior one according to the method of such authors.
The impact of the mismodeling in the third even zonal (ℓ = 6) It turns out that, even by restricting to ℓ = 6 alone, whose differences ∆C 6,0 for some gravity models are displayed in Table 3, the corresponding mismodelled node precession of LARES δΩ LR J6 , calculated with the value in Table 1, is as large as 104 mas yr −1 for the pair ITU GRACE16/GOCO05S, while it amounts to 77 mas yr −1 for the pair ITU GRACE16/ITSG-Grace2014s, corresponding to 15%, 11%, respectively, of the combined Lense-Thirring signature of Eq. (16). A smaller value can be obtained by considering ITU GRACE16 and JYY GOCE04S; indeed, it is δΩ LR J6 = 60 mas yr −1 , corresponding to 9% of the predicted relativistic trend. The complete list of percent uncertainty due to the errors in J 6 is quoted in Table 4. Let us, now, adopt the more advanced approach to realistically assess the uncertainty in the even zonals proposed by Wagner and McAdoo in Ref. [14]. Eq. (A11) and Eq. (A13) of Ref. [14] yield the resulting (squared) scale coef- to be applied to the sigma of the Stokes coefficient of degree ℓ of the field labeled with "test" in order to have a realistic evaluation of its uncertainty; the model labeled with "ref" has a superior formal accuracy, and is used as a reference (calibrating) gravity field [14]. By adopting ITU GRACE16 as formally superior reference model and GOCO05S as test model, it turns out that the sigma of GOCO05S, quoted in Table 2, has to be scaled by g test, 6 = 96.81.
It can be noted that Eq. (31) yields a realistic uncertainty for C 6,0 very close to the simple difference ∆C 6,0 between the estimated coefficients of ITU GRACE16 and GGM05S. Eq. (32) returns a scaling factor which is less than half the previous one; it corresponds to an overall uncertainty of 6% of the combined Lense-Thirring signature. If, instead, ITSG-Grace2014s is assumed as test model by keeping ITU GRACE16 as reference model, we have f test, 6 = 125.65, g test, 6 = 73.21.
In view of the sigma of ITSG-Grace2014s reported in Table 2, Eqs. (33)-(34) imply a mismodeled LARES node precession δΩ LR J6 as large as 78 mas yr −1 and 45 mas yr −1 , respectively, corresponding to an overall uncertainty in the combined Lense-Thirring trend of 11% and 7%, respectively. Also in this case, Eq. (29) yields essentially the same result as the simple difference between the nominal coefficients of the models considered. If the JYY GOCE04S global Table 5 Absolute values of the differences ∆C 8,0 among the estimated values C 8,0 of the models of Table 2. They are very close to the uncertainties retrievable from the application of Eq. (A11) by Wagner and McAdoo [14] to the sigma of the model to be tested by contrasting it to a formally superior one according to the method of such authors.

2.2.2
The impact of the mismodeling in the fourth even zonal (ℓ = 8) Also the mismodeling inΩ LR J8 may give a non-negligible contribution, to be added to that due to ℓ = 6. Indeed, the differences ∆C ℓ,0 of the nominal values of the even zonal with ℓ = 8 for the pairs ITU GRACE16/JYY GOCE04S and GOCO05S/JYY GOCE04S listed in Table 5, along with the precessional coefficients of Table 1, yield a mismodelled precession as large as δΩ LR J8 = 40 mas yr −1 , corresponding to δµ LT = 6% of Eq. (16). See Table 6 for the complete list of percent uncertainties δµ LT due to all the models considered for ℓ = 8. The application of the more sophisticated method by Wagner and McAdoo [14] to the models GOCO05S (reference) and JYY GOCE04S (test) yields Such scaling factors must be applied to the sigma of JYY GOCE04S, listed in Table 2, providing us with an overall uncertainty in the combined relativistic effect as large as 6% and 2%, respectively. Also in this case, the scaling factor of Eq. (29) yields a realistic uncertainty very close to the one coming from the simple difference of the estimated values; cfr. Table 5.
L. Iorio Table 6 Percent uncertainty δµ LT (%) due the errors in C 8,0 according to Table 5. Table 7 Absolute values of the differences ∆C 10,0 among the estimated values C 10,0 of the models of Table 2. They are very close to the uncertainties retrievable from the application of Eq. (A11) by Wagner and McAdoo [14] to the sigma of the model to be tested by contrasting it to a formally superior one according to the method of such authors.

2.2.3
The impact of the mismodeling in the fifth even zonal (ℓ = 10) From Table 1 and Table 7 , it turns out that also the even zonal of degree ℓ = 10 may act as a source of non-negligible systematic bias. Indeed, Table 8 It turns out that the uncertainty resulting from the scaling of the sigma σ C 10,0 of JYY GOCE04S, quoted in Table 1, by Eq. (39) is essentially identical to the difference ∆C 10,0 between the estimated values for ITU GRACE16 and JYY GOCE04S. On the other hand, using the smaller scaling factor of Eq. (40) returns a bias δµ LT as large as 13%. Table 9 Largest values of the percent uncertainties δµ LT (%) by degree ℓ according to Tables 4, Table 6, and Table 8. As pointed out at the end of Section 2.1, the linear combination approach allows, in principle, to determine also the corrections δJ 2 , δJ 4 to the nominal values of the specific background reference model adopted for the geopotential in the LAGEOSs' data reduction, or even J 2 , J 4 themselves if a truncated model not including them at all is used. Such a test would be of great relevance since the results obtained in this way could be compared to those inferred in the GRACE/GOCE-based global solutions. At present, it has not been implemented.
Remarkably, Ciufolini et al. [4] actually estimated some geopotential coefficients as solve-for parameters in a past data reduction of the LAGEOS/LARES observations. If, on the one hand, this fact further strengthens our concerns about the continuing lack of an explicit determination of a dedicated Lense-Thirring parameter itself, on the other hand, as remarked in Ref. [11], the results by Ciufolini et al. [4] for the estimated Stokes coefficients were up to two orders of magnitude less accurate than in the latest GRACE/GOCE-based global solutions.

Summary and conclusions
The recently published test of frame-dragging with the terrestrial geodetic satellites of the LAGEOS family by Ciufolini et al. is flawed by the same fundamental issues as many of the other previous works by the same authors.
As in previous occasions, the linear combination of the nodes of the three satellites used as observable, explicitly published by the present author in several papers since 2005 on the basis of a method proposed by Ciufolini in 1996, was attributed also this time by Ciufolini et al. to themselves.
Our analysis shows that the uncertainties in the uncancelled even zonals of degree ℓ = 6, 8, 10, computed both by taking the simple differences in the estimated geopotential's coefficients in satellite-based global models for such degrees and with a method published in the literature by Wagner and McAdoo a few years ago to quantitatively assess their realistic accuracy, may induce a systematic bias up to 15% (ℓ = 6), 6% (ℓ = 8), 36% (ℓ = 10) for certain models, as resumed by Table 9, because of the scatter in the determined values of the corresponding geopotential coefficients from several recent global gravity fields produced by different institutions worldwide.
As a further consistency test of the approach followed so far, we suggest to use suitably designed linear combinations to determine also the corrections δJ 2 , δJ 4 to the reference background geopotential model used in the reduction of the LAGEOSs data.
Finally, it is remarkable that, after about twenty years since the first reported tests with LAGEOS and LAGEOS II and four years since the launch of LARES, nobody has yet published any genuinely independent test of the Lense-Thirring effect with such geodetic satellites in the peer-reviewed literature, especially in view of how many researchers around the world constitute the global satellite laser ranging community.