On testing frame-dragging with LAGEOS and a recently announced geodetic satellite

Recently, Ciufolini and coworkers announced the forthcoming launch of a new cannonball geodetic satellite in 2019. It should be injected in an essentially circular path with the same semimajor axis $a$ of LAGEOS, in orbit since 1976, and an inclination $I$ of its orbital plane supplementary with respect to that of its existing cousin. According to their proponents, the sum of the satellites' precessions of the longitudes of the ascending nodes $\Omega$ should allow one to test the general relativistic Lense-Thirring effect to a $\simeq 0.2\%$ accuracy level, with a contribution of the mismodeling in the even zonal harmonics $J_\ell,~\ell=2,4,6,\ldots$ of the geopotential to the total error budget as little as $0.1\%$. Actually, such an ambitious goal seems to be hardly attainable because of the direct and indirect impact of, at least, the first even zonal $J_2$. On the one hand, the lingering scatter of the estimated values of such a key geophysical parameter from different recent GRACE/GOCE-based global gravity field solutions is representative of an uncertainty which may directly impact the summed Lense-Thirring node precessions at a $\simeq 70-80\%$ in the worst scenarios, and to a $\simeq 3-10\%$ level in other, more favorable cases. On the other hand, the phenomenologically measured secular decay $\dot a$ of the semimajor axis of LAGEOS (and, presumably, of the other satellite as well), currently known at a $\sigma_{\dot a}\simeq 0.03~\textrm{m~yr}^{-1}$ level after more than 30 yr, will couple with the sum of the $J_2$-induced node precessions yielding an overall bias as large as $\simeq 20-40\%$ after $5-10$ yr. A further systematic error of the order of $\simeq 2-14\%$ may arise from an analogous interplay of the secular decay of the inclination $\dot I$ with the oblateness-driven node precessions.


Introduction
The cannonball geodetic satellites of the LAGEOS family, i.e. LAGEOS (L), LAGEOS II (L II) and LARES (LR), entirely covered by passive retroreflectors and tracked on a continuous basis from several ground stations scattered around the world with the Satellite Laser Ranging (SLR) technique (Combrinck 2010), are currently used, among other things, to put to the test some predictions of the Einstein's General Theory of Relativity (GTR) (Combrinck 2011(Combrinck , 2013Peron 2014;Lucchesi et al. 2015;Ciufolini et al. 2015). One of them is known as 1 Lense-Thirring (LT) effect (Lense & Thirring 1918), and its measurement is one of the current goals of the LAGEOS-type satellites in fundamental physics. It consists of small secular precessions of some orbital elements of a test particle in geodesic motion around a rotating primary which, in the case of the aforementioned spacecraft, amount to some dozens milliarcseconds per year (mas yr −1 ). Such long-term orbital rates of change are induced by the post-Newtonian (pN) gravitomagnetic field (Thorne, MacDonald & Price 1986;Thorne 1988) of the central body. It is generated by the mass-energy currents of its angular momentum S, and it is believed to play an important role in several dynamical effects taking place around spinning Kerr black holes (Penrose & Floyd 1971;Thorne 1988;Williams 2005). For an overview of the so-far performed attempts to measure it with artificial satellites in the field of the Earth, see, e.g., Iorio et al. (2011); Ciufolini et al. (2013); Renzetti (2013), and references therein.
The space environment of an astronomical body like the Earth is characterized by several other competing accelerations, both of non-gravitational (Sehnal 1975;Milani, Nobili & Farinella 1987b;de Moraes 1994) and gravitational (Lambeck, Cazenave & Balmino 1974;Rosborough 1986;Kaula 2000) origin, some of which have just the same temporal signature of the relativistic ones of interest. In view of their relatively small magnitudes with respect to the Newtonian gravitational monopole of the Earth, they can be treated perturbatively with the standard methods of perturbation theory and celestial mechanics; see, e.g., ; Bertotti, Farinella & Vokrouhlický (2003a); Kopeikin, Efroimsky & Kaplan (2011); Poisson & Will (2014). All of them contribute in determining the actual path followed by a test particle which, thus, deviates more or less notably from the Keplerian ellipse (Capderou 2014). As a consequence, extracting the gravitomagnetic signature and assessing realistically the uncertainty with which such a task can be implemented is not easy, and requires a careful analysis of all such other biasing disturbances. Since the gravitomagnetic effects are linear trends (see Equation (1) below), this is particularly true for those perturbations which are either secular rates too, like those due to the asphericity of the Earth's gravitational potential (Rosborough 1986;Kaula 2000) (see Equation (2) below) or to some non-gravitational accelerations (Lucchesi 2002), or are harmonic with so small characteristic frequencies that they may mimic the action of biasing linear trends over observational time spans which cover just a fraction of their period of variation like certain tidal perturbations (Pavlis & Iorio 2002). In 1976, van Patten & Everitt (1976a noticed that the nodes Ω of two counter-revolving satellites sharing the same orbital parameters undergo identical secular LT precessionṡ c 2 a 3 1 − e 2 3/2 (1) and opposite secular Newtonian rates of change due to the even zonal harmonic coefficients J ℓ , ℓ = 2, 4, 6 . . . of the multipolar expansion of the Earth's gravitational potential. In Equation (1), G, c are the Newtonian constant of gravitation and the speed of light in vacuum, respectively, while a, e are the semimajor axis and the eccentricity of the satellite's orbit, respectively. The largest one is induced by the first even zonal harmonic J 2 : it iṡ and its nominal value is usually several orders of magnitude larger than the gravitomagnetic one of Equation (1); see Table 1 for the currently orbiting satellites of the LAGEOS family, their relevant orbital parameters and precessions. In Equation (2), R is the Earth's mean equatorial radius, while n b , I are the Keplerian mean motion and the inclination to the Earth's equator of the satellite's orbit, respectively. At that time, the state-of-the-art of modeling the Earth's geopotential would not have allowed one to use the node of a single satellite because of the still too large uncertainty in J 2 . It would have yielded a mismodelled node precession which would have completely overwhelmed the LT one in view of its much larger size and identical temporal pattern. Such a state of affairs still persists today, despite the steady efforts in producing global gravity field models of ever increasing accuracy by several institutions throughout the world. Unfortunately, there are no other Keplerian orbital elements affected by both the gravitomagnetic field of the Earth and the geopotential with different temporal patterns, so that they could be used to decouple their mutual impact. Indeed, only the perigee undergoes a secular Lense-Thirring rate; actually, apart from the fact that it is much more heavily impacted by the non-gravitational perturbations than the node, also the even zonal harmonics J ℓ , ℓ = 2, 4, 6 . . . induce just the same kind of linear trend on it (Kaula 2000). Thus, van Patten & Everitt (1976a,b) considered the sum of the nodes of both their counter-orbiting spacecraft, which should have been endowed with drag-free apparatus to counteract the non-gravitational perturbations. Indeed, such an arrangement would allow, at least in principle, to exactly cancel out the classical perturbations due to the even zonals and add up the LT ones. In 1986, Ciufolini (1986) put forth an essentially equivalent version of the scenario by van Patten & Everitt (1976a,b) suggesting to launch a new passive geodetic satellite X with the same orbital configuration of LAGEOS, launched in 1976, apart from the inclination I X which should have been displaced by 180 deg from that of the already existing spacecraft. Such a "butterfly" orbital geometry has the same main features of that by van Patten & Everitt (1976a,b) as far as the classical and relativistic node precessions are concerned. Thus, also Ciufolini (1986) proposed to monitor the sum of the nodes of LAGEOS and of the proposed new spacecraft X. The non-gravitational perturbations affecting the nodes would not have posed a too severe threatening to the LT measurement because of their relatively accurate modeling for cannonball satellites like the LAGEOS-type ones, being at the 1% level of the LT. Such a conclusion is still valid today; see, e.g., Sehnal (1981); Lucchesi (2001Lucchesi ( , 2002; Lucchesi et al. (2015); Lhotka, Celletti & Galeş (2016); Pardini et al. (2017), and references therein.
In the following years, LAGEOS II and LARES were actually launched, but none of them in the orbital configuration desired by Ciufolini (1986) for his satellite which assumed various provisional names over the years like Lageos-3, LARES/WEBER-SAT. Various LT tests conducted by combining data of LAGEOS and LAGEOS II (Ciufolini et al. 2013), and more recently also of LARES (Ciufolini et al. 2016), have been reported so far by Ciufolini and collaborators; there is currently a lingering debate on several issues which would plague them like their realistic overall accuracy (Iorio et al. 2011;Ciufolini et al. 2013;Renzetti 2013).
Recently, Ciufolini et al. (2017a) announced that a further LAGEOS-type satellite, which we will denote as CiufoLares (CL) in honor of its proponent instead of the rather anodyne LARES 2, is planned for launch in 2019 with the new VEGA C rocket. This time, the orbit of the forthcoming SLR target seems to be the right one: indeed, it should match the originally proposed butterfly configuration with LAGEOS, up to allowed offsets in a and I as little as ∆a CL = 20 km, ∆I CL = 0.15 deg. Ciufolini et al. (2017a) claimed an overall accuracy in the LT measurement of the order of 0.2%, with a total contribution from the uncertainties in the even zonal harmonics as little as ≃ 0.1%.
In our preliminary analysis, we will show that, actually, the overall impact of just the first even zonal harmonic of the geopotential, including both its direct effect due to the mismodeling in J 2 and the indirect one due to the interplay with the measured secular decay of the semimajor axis, and, perhaps, of the inclination as well, of LAGEOS and, likely, of CL as well, may be up to 200 − 800 times larger over a data analysis 5 − 10 yr long. The paper is organized as follows. In Section 2, we will deal with the currently existing scatter in the estimated values of the Earth's oblateness from the latest global gravity field solutions produced by several independent institutions. Instead of taking into account the more or less realistically re-scaled sigmas of just a single, favorable Earth's gravity field, which is also likely a-priori imprinted by the LT itself, by comparing the estimated values of the coefficient C 2,0 of several pairs global gravity models of similar accuracy, we will show that the total bias in the LT summed node rates can be as large as up to ≃ 10 − 80% for certain geopotential models. Here, C 2,0 is the fully normalized Stokes coefficient of degree ℓ = 2 and order m = 0 of the multipolar expansion of Earth's gravitational potential (Petit, Luzum & et al. 2010). It is related to the first even zonal harmonic of the geopotential by J 2 = − √ 5 C 2,0 . Section 3 is devoted to quantitatively assessing the indirect effect of the secular orbital decays measured so far for all the existing satellites of the LAGEOS family on their node precessions due to J 2 showing that the resulting systematic bias can reach ≃ 20 − 40% of the combined LT node signal after 5 − 10 yr. An analogous effect due to a possible secular decay of the inclination, yielding a bias of the order of ≃ 2 − 14%, is dealt with in Section 4. Section 5 summarizes our findings and offers our conclusions. A list of definitions and conventions used throughout the text is displayed in Appendix A for the benefit of the reader. Appendix B contains the tables and the figures.
2. The direct impact of the mismodeling of the even zonal harmonics of the geopotential By introducing the coefficientΩ the mismodelled part of the sum of the node precessions of L and CL due to our imperfect knowledge of the Earth's oblateness can be calculated as where δJ 2 represents some quantitative measure of the actual uncertainty in the first even zonal harmonic of the geopotential. Note that, since J 2 is an overall external parameter not depending on the satellite, the orbital configuration of CL is crucial in the assessment of Equation (4) since the coefficient proportional to δJ 2 is made just of the sum of the J 2 −induced node rates of L and CL. In principle, if the orbital parameters of CL were exactly equal to their nominal values, Equation (4) would vanish sinceΩ L .2 = −Ω CL .2 . The percent error ψ in the sum of the Lense-Thirring node rates can be straightforwardly evaluated as the ratio of Equation (4) to the sum of the nominal gravitomagnetic node precessionsΩ i.e. ψ δ Ω tot The key factors in Equation (4) will be the actual departures ∆a CL , ∆I CL of the CL orbit with respect to the ideal one, and a realistic evaluation of the lingering uncertainty in J 2 . About the orbital configuration of CL, it should be noted that the inclination of the existing LARES satellite, which should have been inserted in an orbital plane supplementary to that of LAGEOS, exhibits an offset of ∆I LR = 0.7 deg with respect to the ideal case; thus, in the following we will prudentially assume for CL ∆I CL = 0.5 deg with respect to the smaller value quoted in Table 1.
As far as the evaluation of the uncertainty in J 2 is concerned, an ever increasing number of global gravity field solutions produced by several independent institutions throughout the world chiefly from data of the dedicated GRACE and GOCE space-based missions is nowadays available. They are made freely available by the International Center for Global Earth Models (ICGEM) which collects them in its webpage at http://icgem.gfz-potsdam.de/tom longtime on the Internet. Thus, it is neither realistic to assume the mere statistical, formal errors σ J 2 released by the various models as a trustable measure of δJ 2 nor to pick up just that model which gives the smallest overall uncertainty in the LT test by discarding other ones if less favorable. Furthermore, it must also be stressed that some models return determinations of C 2,0 obtained directly from the constellation of the existing SLR satellites among which the LAGEOS family plays an important role; such values must be, of course, discarded to avoid any possible a priori "imprint" of GTR itself (Ciufolini 1996(Ciufolini , pag. 1718. Finally, in order to meaningfully compare the values of the first even zonal retrieved from different models, it is important that they share the same tide system (zero-tide or tide-free). For an explanation of such definitions, see Petit, Luzum & et al. (2010, Sect. 1.1 and 6.2.2). Actually, Ciufolini et al. (2017b) did not take ito account any of such issues. Indeed, instead of considering several different global gravity field solutions, they came to their claimed 0.1% error due to the even zonals by using only a single Earth's gravity model, i.e. GOCO05s (Mayer-Gürr & the GOCO consortium 2015), by using its 1 − σ formal errors, apart from C 2,0 whose uncertainty was specifically evaluated in a quite hand-waving, confusing and arbitrary fashion. Suffice it to say that Ciufolini et al. (2017b) resorted to a historical time series of its SLR-based determinations covering a temporal interval  which will necessarily have nothing to do with L and CL. As a further drawback, such a model relies upon data records from CHAMP, GRACE, GOCE and six SLR satellites (LAGEOS, LAGEOS II, Starlette, Stella, Ajisai, Larets). The data from the geodetic satellites are routinely used in some global gravity field models to determine more accurately just the even zonals of low degree. Thus, the values for C 2,0 from GOCO05s are unavoidably plagued by an a-priori "imprint" of the LT effect itself which mainly concentrates just in the even zonals of the lowest degree (Ciufolini 1996(Ciufolini , pag. 1719. As such, GOCO05s should not be considered as a viable background model to be used in dedicated data reductions to measure the gravitomagnetic field involving LAGEOS itself. Furthermore, also the evaluation of the error in C 2,0 proposed by Ciufolini et al. (2017b) should be deemed as a priori strongly imprinted by the LT itself since it entirely relies upon SLR data. It would have been much more useful and appropriate if Ciufolini et al. (2017b) had simulated a realistic data reduction of LAGEOS and the new satellite by simultaneously estimating, among other things, both a LT parameter, or, even better, the Earth's angular momentum S itself, and C 2,0 on some predetermined temporal basis (say, weekly, monthly, etc.), and the resulting correlations between S and C 2,0 had been inspected. Such a procedure should have been repeated for several different Earth's global gravity field solutions not including previous SLR data from LAGEOS itself as background models. Inexplicably, it has never been yet implemented, at least publicly, by any group, not even with the real data of the existing satellites of the LAGEOS family, after more than 20 yr since the first published tests . About the issue of the a priori "imprint" of the LT in the existing estimated values of C 2,0 , it should be noted that, although indirectly, it may somewhat affect also some models based solely on GRACE/GOCE data. Indeed, many of them use previously obtained global solutions which rely upon just SLR data for the low-degree even zonals as background gravity models.
In Tables 2 to 3, we adopt the differences ∆C 2,0 between the nominal values of the estimated coefficients C i/ j 2.0 of several pairs of geopotential models i, j, not directly relying upon the SLR data of the LAGEOS family itself, for the uncertainty δJ 2 entering Equation (4); (e.g. Ciufolini 1996Ciufolini , pag. 1713. From Table 2, which compares some models of the zero-tide system built from GRACE data records differing by their temporal extensions and type of measurements, it turns out that the maximum uncertainty in the systematic bias due to the first even zonal can reach the ≃ 10% level by using the GRACE-based HUST-Grace2016s (Zhou et al. 2017) and ITU GRACE16 (Akyilmaz et al. 2016) models; cfr. with Figure 1. Interestingly, the formal errors of such solutions are comparable (see Ciufolini et al. 2009, pag. 91 so that σ ITU GRACE16 HUST-Grace2016s is a static global gravity field model obtained by using approximately 13 yr (spanning from January 2003 to April 2015) of GRACE-only data (Zhou et al. 2017). ITU GRACE16 is a static global gravity field model computed from GRACE Satellite-to-Satellite-Tracking (SST) data of 50 months collected between April 2009 to October 2013 (Akyilmaz et al. 2016). For the details of the other models used in Table 2, see the dedicated entries at the ICGEM Webpage http://icgem.gfz-potsdam.de/tom longtime. Table 3, showing the impact of some tide-free models from various GOCE data records analyzed with different approaches, tells us that the J 2 -driven overall uncertainty in the LT test can be as large as ≃ 85% if the GOCE-based IGGT R1 (Lu et al. 2018) and IfE GOCE05s (Wu, Müller & Brieden 2017) models are considered; see also Figure 2. Moreover, the (formal) released sigmas of IGGT R1 and IfE GOCE05s are comparable since so that σ IGGT R1 IGGT R1 (Lu et al. 2018) relies upom the use of the three invariants of the gravitational gradient tensor (IGGT) to process about 1 yr of GOCE data (2009)(2010), while in IfE GOCE05s (Wu, Müller & Brieden 2017) almost the same GOCE-only data set is analyzed with either the SST and the Satellite Gravity Gradient (SGG) techniques. For the details of the other tide-free models used in Table 3, see their dedicated entries at http://icgem.gfz-potsdam.de/tom longtime.
It is difficult to properly argue about the differences between Table 2 and Table 3, also because, after all, new global gravity field solutions will be available if and when the LT-dedicated L-CL analysis will be performed. Our results should likely be looked just as illustrative of the current state-of-affairs and of the need of not limiting just to one particular model, given the sound possibility that the increasing number of values of J 2 will not finally converge to the desired level of accuracy.
A further source of potentially non-negligible systematic uncertainty is as follow. The precessions of Equations (1) to (2), on which Equation (6) is based, hold in a coordinate system whose reference z axis is aligned with the Earth's spin axis. Actually, real data reductions are performed with respect to the International Celestial Reference Frame (ICRF) (Petit, Luzum & et al. 2010) which adopts the mean Earth's equator at epoch J2000.0 as fundamental reference {x, y} plane. In fact,Ŝ does vary in time because of a variety of physical phenomena among which the lunisolar torques inducing the precession of the equinoxes is a major one. Thus, in correctly assessing the systematic bias due to J 2 one has to properly account for the fact that, at the time of the data reduction,Ŝ will be displaced by a certain amount with respect to the case of Equations (1) to (2), and adequate formulas have to be used. To this aim, the general gravitomagnetic and J 2 -driven node precessions are (Iorio 2017) In Equations (11) to (12),m is a unit vector in the orbital plane directed transversely to the line of the nodes, whilen is a unit vector perpendiculrar to the orbital plane along the orbital angular momentum. AboutŜ, its temporal variation, which introduces a time-dependence in the LT and J 2 rates of Equations (11) to (12), can be expressed as (Seidelmann et al. 2007) δ ≃ 90 deg − 0.557 deg cty −1 τ.
In Equations (13) to (14), α, δ are the right ascension (RA) and the declination (DEC) of the Earth's spin axis, respectively, while τ is the interval in Julian centuries (of 36525 days) from the standard epoch JD 2451545.0. The Earth's spin axis can be expressed in terms of α, δ aŝ S = {cos α cos δ, sin α cos δ, sin δ} .
Thus, if one calculates Equation (6) by means of Equations (11) to (14), a dependence on t and the initial values Ω L 0 , Ω CL 0 of the nodes of L and CL at the starting epoch of the data reduction which may have a non-negligible impact is introduced. By expanding Equation (15) according to Equations (13) to (14) and neglecting terms quadratic inα,δ, a linear dependence on t occurs.
It implies a quadratic signature in the integrated node residuals giving rise to a mild parabolic residual signal. Over temporal intervals just some years long, it would likely superimpose to the linear LT trend by potentially corrupting its recovery. A similar issue occurs also for other potential sources of systematic errors (see Sections 3 to 4). Figure 3 depicts the scatter in the values of ψ over 10 yr for δC 2,0 = 2.5 × 10 −10 (see Table 2) and ∆a CL = 20 km, ∆I CL = 0.5 deg by allowing Ω L 0 , Ω CL 0 to vary independently of each other within a 360 deg range. It can be noted that ψ finally vary from 5% to about 15%.

The indirect impact of the Earth's oblateness through the secular decay of the satellites' semimajor axes
It has been known for several decades that the semimajor axis a of the existing LAGEOS-type satellites experiences a secular decayȧ < 0 as large as (Rubincam 1982;Sośnica 2014;Sośnica et al. 2014;Pardini et al. 2017) |ȧ| ≃ 0.2 − 0.9 m yr −1 experimentally determined with an accuracy of the order of (Rubincam 1982;Sośnica 2014;Sośnica et al. 2014) σȧ ≃ 0.03 − 0.1 m yr −1 ; see Table 2 of Iorio (2016). About LAGEOS, it should be remarked that the experimental accuracy in determining its orbital decay rate has not significantly improved over the past few decades. Indeed, Rubincam (1982) reported σȧL = 0.1 mm d −1 = 0.036 m yr −1 , while σȧL = 0.035 m yr −1 is quoted in Sośnica (2014); Sośnica et al. (2014). Note that the quoted σȧ are not due to a priori uncertainties in the acceleration models used; instead, they are a posteriori errors generated in the data reduction procedure.
Neither the pN gravitomagnetic field of the Earth nor its Newtonian oblateness J 2 directly affect the semimajor axis a with secular perturbations. Nonetheless, as shown by Iorio (2016), the secular decayȧ < 0 of such an orbital element has a long-term indirect effect on the node Ω through its interplay with its classical precessionΩ J 2 driven by J 2 . The resulting change over a time span T is (Iorio 2016) In Equation (19), a 0 is the semimajor axis at some reference epoch assumed as initial value at the beginning of the observational time span T . To the first order inȧ, Equation (19) can be written as where It is quite plausible and reasonable to assume that also CL will finally experience such a subtle orbital decay, although much will likely depend on the actual manufacturing of CL and its surface properties. Thus, it is important to calculate its impact on the proposed use of the sum of the nodes of LAGEOS and CL. It is worthwhile remarking that we will exclusively rely upon the so-far phenomenologically determined values of the secular decays of the semimajor axes of the satellites of the LAGEOS family from real data reductions of long observational records. Stated differently, we will not try to model such an observed orbital feature in terms of some exotic or mundane physical phenomena; for several attempts towards this goal in terms of standard non-gravitational effects, see, e.g., Rubincam (1982Rubincam ( , 1987Rubincam ( , 1988; Scharroo et al. (1991); Lucchesi et al. (2015); Pardini et al. (2017). Since each satellite experiences its own semimajor axis secular rate,ȧ must be treated as an independent variable for both satellites in calculating the associated error in the sum of their node shifts which, thus, reads Note that it was assumed that the errors σȧL, σȧCL stay constant during the data analysis. Should they change, things would be even more complicated. In this case, contrary to the bias due to the Earth's even zonal J 2 , the peculiar orbital configuration of CL does not provide benefits in calculating the systematic error due to σȧ. Indeed, ifȧ were, say, an external parameter common to both satellites as it occurs for J 2 , the coefficient proportional to σȧ would be K L J 2 + K CL J 2 . As a result, the corresponding percent error ξ δ ∆Ω toṫ a ∆Ω tot LT (23) in the total Lense-Thirring node shift, for given values of ∆a CL , ∆I CL and σȧL, σȧCL, is linear in T . Figure 4 depicts it over a time shift 10 yr long for ∆a CL = 10 km, ∆I CL = 0.15 deg, σȧL = 0.035 m yr −1 . For CL, we hypothesize an improvement of, say, a factor of 3 in the determination of its possible orbital rate decay with respect to LAGEOS by assuming σȧCL = 0.01 m yr −1 . It can be noted that ξ is as large as 20% after 5 yr, reaching 40% after 10 yr. It turns out that, even if it were σȧCL = 0, the situation would not change too much, with a maximum bias of roughly ≃ 35% after 10 yr. A breakthrough in the accuracy of the determination of the orbital decay rate of LAGEOS by a factor of 10 would be needed to bring such source of systematic uncertainty down to the percent level; in view of what happened so far in the last decades, it does not seem plausible to occur in the foreseeable future. At least as far as the existing LAGEOS is concerned, the experimental error σȧL has shown so far a very weak time dependence, staying essentially constant over the decades.
For the sake of completeness, we also note that the impact of the uncertainty in J 2 in the sum of the corresponding node shifts due toȧ of Equation (19) is negligible with respect to the added Lense-Thirring node rates.
Finally, it is worth noting that the effect considered in this Section would act as some sort of parabolic signature in the integrated node residuals, being quadratic in time. Nonetheless, it would be premature and unjustified to argue that its mismodelled component would necessarily decouple from the LT linear signal. Indeed, the resulting parabola would be rather flat, especially over not too long temporal intervals T . Thus, it is likely that the recovery of the gravitomagnetic slope would be impacted by the aforementionded mismodeled competing quadratic effect. A quantitative assessment of the level of such a potential bias is beyond the scope of the present paper since it would require ad-hoc data simulations and reductions.

The indirect impact of the Earth's oblateness through the secular decay of the satellites' inclination
Despite it seems that, at present, no experimental determinations of a possible secular decaẏ I of the inclination of the LAGEOS-type satellites are publicly available in the literature, there are several standard physical mechanisms of non-gravitational origin which, in principle, are able to induce such an effect. Thus, despite the content of the present Section may be deemed as still hypothetical, we will prefer to accurately consider also its possible indirect impact on the classical node precessions driven by the Earth's oblateness as done for the semimajor axis decay in Section 3.

By inserting
in Equation (2), the total node shift including also the part due to the secular decay of the inclination, integrated over a time span T , is In Equations (24) to (25), I 0 is the inclination at some reference epoch assumed as initial value at the beginning of the observational time span T . Note that, in the limitİ → 0, Equation (25) reduces just to Equation (2). By posing as it is the case for the LAGEOS-type satellites over even multidecadal time spans T (see Equation (47) and Equation (49) below), let us expand the trigonometric functions sin α, cos α in sin (I 0 + α) entering Equation (25) as Thus, from Equation (25) one can approximately obtain for theİ-induced node shift From Equation (29), it turns out By posing the total systematic bias in the sum of the nodes is, thus, Not only the experimental accuracy σ˙I in measuring a possible decay of the orbital plane, but also its nominal valueİ enter Equation (32) through Equation (31). Thus, at least in principle, it is important to try to give a plausible order of magnitude of such a subtle phenomenon. Below, we will consider only the effect of the neutral and charged drag, whose disturbing acceleration can be modeled as (King-Hele 1987;Milani, Nobili & Farinella 1987b) In Equation (33), C D , Σ, ρ, V are the satellite's drag coefficient and area-to-mass ratio, the atmospheric density at the satellite's height, and the satellite's velocity with respect to the Earth's atmosphere, respectively. It can be shown that, among other things, Equation (33) also induces a secular variation of the inclination I given by (Milani, Nobili & Farinella 1987b;Iorio 2010) In Equation (34), ω atm is the angular velocity of Earth's atmosphere. We will, first, confirm the validity of Equation (34) in the case of LARES by comparing its predicted theoretical rate with a numerical integration of its equations of motion. For such a satellite it is (Pardini et al. 2017) ρ neutral = 5.9 × 10 −16 kg m −3 , so that Equation (34) allows to predict a nominal decay of the order of (Iorio 2010) In Equation (38), ω ⊕ is the angular velocity of Earth. The analytical result of Equation (39) is supported by a numerical integration of the equations of motion of LARES in rectangular Cartesian coordinates with and without Equation (33). Figure 5 displays the resulting time series for the drag-induced inclination shift ∆I drag (t) of LARES over a time span 1 yr long; a negative trend with just the same size of Equation (39) is clearly apparent. It can be shown that our approach is also able to reproduce, both analytically and numerically, the observed secular decay of the semimajor axis of LARES, in agreement with the finding by Pardini et al. (2017) who found that the neutral atmospheric drag is able to explain almost entirely such a phenomenon. Thus, confident of our strategy, we can, now, safely apply it to LAGEOS and CL. At the altitude of 5900 km, the total neutral atmospheric density experienced by LAGEOS is of the order of ρ L neutral ≃ 8.4 × 10 −18 kg m −3 (40) due to neutral hydrogen H, corresponding to a hydrogen number density of 5 × 10 9 m −3 (Rubincam 1990). As far as the charged particles in the plasmasphere are concerned, comes from protons H + for a proton number density of 3 × 10 9 m −3 (Rubincam 1990), refers to Helium ions He + , while is attributable to Oxygen ions O + . However, as noted in Rubincam (1990), Equations (41) to (43) have to be divided by a factor of 2 since LAGEOS spends just half its time in the plasmasphere. Thus, the total charged density actually felt by LAGEOS is about In calculating the acceleration due to charged drag with Equation (33), the charged drag coefficient of LAGEOS can be as large as (Milani, Nobili & Farinella 1987b) C ch D ≃ 20, while its area-to-mass ratio amounts to Σ = 7 × 10 −4 m 2 kg −1 .
Thus, both Equation (34) and a numerical integration of the equations of motion, whose outcome is displayed in Figure 6, return a secular decay of the orbital plane of LAGEOS due to neutral and charged atmospheric drag of the order oḟ In the case of CL, since from the satellite's physical parameteres released by Ciufolini et al. (2017a) it can be inferred Σ = 4 × 10 −4 m 2 kg −1 , the drag-induced inclination decay would amount tȯ by assuming the same value of Equation (45) for its C ch D . Incidentally, Equation (47) and Equation (49) confirm the validity of the assumption made in Equation (26). Also thermal effects connected with the orientation and magnitude of the satellite's spin axis (Rubincam 1987) can induce a secular decay on I (Lucchesi 2002); thus, the actual overall size of such a phenomenon may be larger than Equation (47) and Equation (49).
For some reasons, onlyȧ has been experimentally investigated so far by satellite geodesists; it is expected that, sooner or later, they will look also at the inclination by releasing a best estimate for the predicted inclination decayİ along with its associated error σ˙I. For the sake of convenience, we may tentatively assume σ˙I ≃ (5 − 10%)İ = 0.035 − 0.07 mas yr −1 ; it is the same percent error of the actually measuredȧ. The present assumption that the expected secular decay of the inclination will be measured with the same fractional accuracy as the decay of the semimajor axis should be regarded just as a, hopefully, plausible guess in view of the lingering lacking of actual experimental determinations of it. Note that Ciufolini et al. (2009) claim to be able to determine the inclination of LAGEOS within an accuracy of just 30 µas = 0.03 mas. About the overall temporal pattern of the bias due to such a potential source of systematic error, the same considerations as in Sections 2 to 3 about the temporal variation of the Earth's spin axis and the cross correlation between J 2 andȧ hold. Figure 7 depicts for given values of ∆a CL , ∆I CL ,İ L ,İ CL , σ˙IL, σ˙ICL over a time shift 10 yr long. For the orbital shifts of CL we assumed ∆a CL = 10 km, ∆I CL = 0.15 deg. As far as the inclination decays and their putative errors, we took σ˙IL/CL = 0.03 mas yr −1 along with Equation (47) and Equation (49). It turns out that the nominal valueİ of the expected inclination decay does not has a relevant impact on the total bias χ; even by rescalingİ by a factor of 10 for both satellites do not alter it. Thus, accurately predicting all the possible contributions toİ does not appear, actually, so important in the present context. On the contrary, σ˙I is of crucial importance. Indeed, while for their previously listed values the total bias can reach 2% of the summed LT shifts, as shown by Figure 7, by increasing them up to, say, σ˙I = 0.2 mas yr −1 would push χ up to to 14%.

Summary and conclusions
According to Ciufolini and coworkers, the sum of the precessions of the nodes of the existing LAGEOS satellite and the future CiufoLares (also known-in a more anodyne fashion-as LARES 2), to be launched in 2019 with a VEGA C rocket, should provide us with a ≃ 0.2% test of frame-dragging in the field of the Earth whose imperfectly known even zonal harmonics should contribute just 0.1% the total error budget. Actually, such an ambitious goal seems difficult to be reached because of the direct and indirect impact of the competing classical node precessions mainly induced by the first even zonal of the geopotential.
On the one hand, according to some of the latest GRACE/GOCE-based global gravity field solutions released by various international institutions in 2016-2017, the lingering scatter of their determinations of such a fundamental geophysical parameter would directly affect the sum of the nodes with an uncancelled mismodelled total precession which, in some cases, may reach even several dozens percentage points of the total Lense-Thirring effect for departures of the orbital elements of CiufoLares as little as 20 km and 0.5 deg. Other more favorable scenarios point towards uncertainties of the order of ≃ 3 − 10%. A further source of systematic uncertainty is given by the dependence on the initial values of the satellites' nodes introduced by the unavoidable displacement of the Earth's spin axis with respect to its direction at the reference epoch J2000.0 when the data reduction will be finally performed in the next years by using the ICRF. The insertion of the new SLR target in an orbit as closest as possible to the ideal butterfly configuration with LAGEOS will be crucial for effectively controlling such an insidious source of systematic error in view of the persisting difficulties in effectively constraining, at least to the desired level of accuracy for a successfull test of the Lense-Thirring effect, the first even zonal of the geopotential.
On the other hand, the classical node precession due to the oblateness of the Earth exerts a further, even subtler, but not less potentially insidious, aliasing effect on the relativistic signature of interest in an indirect way through its coupling with the secular decay which has been phenomenologically measured so far for the semimajor axis of all the existing members of the LAGEOS family. Indeed, it turned out that the associated bias on the sum of the nodes is unaffected by the peculiar orbital geometry of CiufoLares, being impacted crucially by the accuracy in determining the satellites' orbital decay. Unfortunately, in the case of LAGEOS, after more than 30 yr, no substantial improvements have occurred so far in improving our knowledge of the rate of decreasing of its semimajor axis. Even by assuming that the forthcoming satellite will experience a more accurately known analogous effect, the resulting systematic uncertainty may reach ≃ 20 − 40% of the Lense-Thirring signal after 5 − 10 yr. An analogous indirect bias should occur due to the expected secular decay of the inclination of the orbital planes, although no experimental records for such an effect on the LAGEOS-type spacecraft are yet available in the literature. Thus, such a potential further source of systematic error should be currently deemed as hypothetical, although plausible and well rooted in estabilished standard non-gravitational physics. Depending on its actual experimental accuracy, its impact on the combined Lense-Thirring signature may be in the range ≃ 2 − 14% after 10 yr. We note that the indirect effects of the Earth's rotation pole position, and the secular decays of the semimajor axis and the inclination would affect the time series of the integrated node residuals with a quadratic temporal signature. This does not necessarily mean that it would neatly decouple from the linear Lense-Thirring effect since the residual parabolic curve would be likely rather smooth and poorly distinguishable from a secular trend, especially over relatively short observational time spans. In order to quantitatively assess how much the recovery of the gravitomagnetic slope would be impacted by the aforementionded mismodeled competing quadratic effects, a full simulation of the time series of the sum of the nodes along with a best fit with a parabolic curve would be required; it is beyond the scopes of the present preliminary analysis.
In conclusion, the desired goal of a ≃ 0.1% test of the gravitomagnetic orbital precessions with LAGEOS and CiufoLares seems, at least at present, rather unlikely to be met. In addition to a very strict adherence of the actual orbit of the new spacecraft to its ideal one, a breakthrough of one order of magnitude in determining both the terrestrial first even zonal harmonic and the secular decrease of the semimajor axis of LAGEOS (and, in perspectives, of CiufoLares as well) would be required to yield a test of some percent by using the sum of the nodes.

Appendix A Notations and definitions
Here, some basic notations and definitions used in the text are presented (Brumberg 1991;Milani, Nobili & Farinella 1987a;Bertotti, Farinella & Vokrouhlický 2003b) G 2,0 of the differences between the estimated normalized Stokes coefficients C 2,0 of degree ℓ = 2 and order m = 0 for the recent GRACE/GOCE-based global gravity field solutions (tide system: zerotide) i, j =Tongji-Grace02s (Chen et al. 2016), ITU GRACE16 (Akyilmaz et al. 2016), HUST-Grace2016s (Zhou et al. 2017), XGM2016 (Pail et al. 2018) retrieved from the section Static Models of the WEB page of the International Center for Global Earth Models (ICGEM) at http://icgem.gfz-potsdam.de/tom longtime. Strictly lower triangular matrix: maximum values of the percent error ψ computed from Equation (6) by assuming the figures in the strictly upper triangular matrix for the uncertainty in C 2,0 . The semimajor axis a and the inclination I of CL were allowed to vary within a range of a CL = 12270 ± 20 km, I CL = 70.16 ± 0.5 deg, respectively. Cfr. with Figure 1.