A HERO for general relativity

HERO (Highly Eccentric Relativity Orbiter) is a space-based mission concept aimed to perform several tests of post-Newtonian gravity around the Earth with a preferably drag-free spacecraft moving along a highly elliptical path fixed in its plane undergoing a relatively fast secular precession. We considered two possible scenarios: a fast, 4-hr orbit with high perigee height of $1,047\,\mathrm{km}$, and a slow, 21-hr path with a low perigee height of $642\,\mathrm{km}$. HERO may detect, for the first time, the post-Newtonian orbital effects induced by the mass quadrupole moment $J_2$ of the Earth which affects the semimajor axis $a$ via a secular trend of $\simeq 4-12\,\mathrm{cm\,yr}^{-1}$, depending on the orbital configuration. Recently, the secular decay of the semimajor axis of the passive satellite LARES was measured with an error as little as $0.7\,\mathrm{cm\,yr}^{-1}$. Also the post-Newtonian spin dipole (Lense-Thirring) and mass monopole (Schwarzschild) effects could be tested to a high accuracy depending on the level of compensation of the non-gravitational perturbations. Moreover, the large eccentricity of the orbit would allow to constrain several long-range modified models of gravity and to accurately measure the gravitational red-shift as well. Each of the six Keplerian orbital elements could be individually monitored to extract the $GJ_2/c^2$ signature, or they could be suitably combined in order to disentangle the post-Newtonian effect(s) of interest from the competing mismodeled Newtonian secular precessions induced by the zonal harmonic multipoles $J_\ell$ of the geopotential. In the latter case, the systematic uncertainty due to the current formal errors $\sigma_{J_\ell}$ of a recent global Earth's gravity field model are better than $1\%$ for all the post-Newtonian effects considered, with a peak of $\simeq 10^{-7}$ for the Schwarzschild-like shifts. [Abridged]


Introduction
The (slow) motion of a test particle moving in the spacetime (weakly) deformed by the mass-energy content of an isolated, axially symmetric rotating body of mass M, angular momentum S, polar and equatorial radii R p , R e , ellipticity ε = 1 − R 2 p /R 2 e , dimensionless quadrupole mass moment J 2 exhibits several post-Newtonian (pN) features. Some of them have never been put to the test so far because of their smallness; they are the gravitoelectric and gravitomagnetic effects associated with the asphericity of the central body induced by its mass quadrupole and spin octupole moments, respectively (Soffel et al. 1987;Soffel 1989;Brumberg 1991;Panhans & Soffel 2014;Meichsner & Soffel 2015).
Instead, the pN orbital effects which have been extensively tested so far in several terrestrial and astronomical scenarios are the gravitoelectric and gravitomagnetic precessions due to the mass monopole and spin dipole moments, respectively. The former is responsible of the time-honored, previously anomalous perihelion precession of Mercury (Le Verrier 1859), whose explanation by (Einstein 1915) was the first empirical success of his newly born theory of gravitation. It was later repeatedly measured with radar measurements of Mercury itself (Shapiro et al. 1972;Shapiro 1990), of other inner planets (Anderson et al. 1978(Anderson et al. , 1992, and of the asteroid Icarus (Shapiro, Ash & Smith 1968;Shapiro et al. 1971) as well. Also binary pulsars (Kramer et al. 2006) and Earth's artificial satellites (Lucchesi & Peron 2010 have been used so far. The latter is the so-called Lense-Thirring effect (Lense & Thirring 1918;Pfister 2007Pfister , 2008Pfister , 2012Pfister , 2014 which is currently under scrutiny in the Earth's surrounding with the geodetic satellites of the LAGEOS family; see, e.g., Renzetti (2013), and Lucchesi et al. (2019), and references therein. Another gravitomagnetic effect-the Pugh-Schiff rates of change of orbiting gyroscopes (Pugh 1959;Schiff 1960)-was successfully tested in the field of the Earth with the dedicated Gravity Probe B (GP-B) spaceborne mission a few years ago (Everitt et al. 2011(Everitt et al. , 2015 to a 19% accuracy level, despite the originally expected one was of the order of 1% (Everitt et al. 2001).
By assuming the validity of general relativity, its mass quadrupole and spin octupole accelerations are, to the first post-Newtonian (pN) order, A pNS J 2 = 3GS R 2 e ε 2 7c 2 r 5 v × 5ξ 7ξ 2 − 3 r + 3 1 − 5ξ 2 Ŝ , where G is the Newtonian constant of gravitation, µ GM is the gravitational parameter of the primary, c is the speed of light in vacuum,Ŝ is the unit vector of the rotational axis, is the cosine of the angle between the directions of the body's angular momentum and the orbiter's position vector r, and v r v ·r, are the components of the particle's velocity v along the radial direction and the primary's spin respectively. The averaged rates of change of the semimajor axis a, the eccentricity e, the inclination I, the longitude of the ascending node Ω and the argument of pericenter ω induced by Equations (1) to (2) were analytically calculated for a general orientation ofŜ in space by Iorio (2015Iorio ( , 2019; previous derivations of the gravitoelectric mass quadrupole effects in the particular case of an equatorial coordinate system with its reference z axis aligned withŜ can be found in Soffel et al. (1987); Soffel (1989); Brumberg (1991); .
In this paper, we will preliminarily explore the perspectives of measuring, for the first time, some consequences of Equations (1) to (2) by suitably designing a dedicated drag-free satellite-based mission around the Earth encompassing a highly eccentric geocentric orbit exploiting the frozen perigee configuration; we provisionally name it as HERO (Highly Eccentric Relativity Orbiter). For some embryonal thoughts about the possibility of using an Earth's spacecraft to measure the pN gravitoelectric effects proportional to GJ 2 /c 2 , see Iorio (2013Iorio ( , 2015; for deeper investigations concerning a possible probe around Jupiter to measure them and the pN gravitomagnetic signature proportional to GS ε 2 /c 2 , see Iorio (2013Iorio ( , 2019. About the propagation of the electromagnetic waves in the deformed spacetime of an oblate body and the perspectives of measuring the resulting deflection due to Jupiter with astrometric techniques, see, e.g., Crosta & Mignard (2006); Kopeikin & Makarov (2007); Le Poncin-Lafitte & Teyssandier (2008); Abbas, Bucciarelli & Lattanzi (2019), and references therein. We will show that the size of the secular rate of a predicted by Equation (1) falls within the recently reached experimental accuracy in measuring phenomenologically such a kind of an effect with the existing passive LARES satellite (Lucchesi et al. 2019). Be that as it may, we will show that, as a by-product, also other general relativistic features of motion could be measured with high accuracy, at least as far as the systematic error due to the current formal level of mismodeling in the competing classical precessions due to the zonal harmonic coefficients J ℓ , ℓ = 2, 3, 4, . . . of the multipolar expansion of the Earth's gravity potential is concerned. To this aim, it is crucial to assess the level of possible cancelation of the non-gravitational perturbations by the drag-free technology. Its evaluation is outside the scopes of the present paper. However, we will look in detail at the atmospheric drag, which is one of the major disturbing non-conservative accelerations inducing relevant competing signatures, especially on a. For the sake of simplicity, we will assume a spherical shape for a passive spacecraft.
The high eccentricity of the suggested orbit of HERO would allow also for accurate tests of the gravitational redshift provided that the spacecraft is endowed with accurate atomic clocks; see Table A2 and Table A6 for the expected sizes of it. For recent tests of such an effect performed with the H-maser clocks carried onboard the satellites GSAT0201 (5-Doresa) and GSAT0202 (6-Milena) of the Galileo constellation by exploiting their fortuitous rather high eccentricity due to their erroneous orbital injection, see Herrmann et al. (2018); Delva et al. (2018).
Finally, several long-range modified models of gravity (Clifton et al. 2012) imply spherically symmetric modifications of the Newtonian inverse-square law which induce net secular precessions of the pericenter and the mean anomaly at epoch. They would represent further valuable goals for HERO.
The paper is organized as follows. Section 2 is the main body of the paper; it discusses two possible orbital configurations along with the magnitude of the various pN effects and the size of the corresponding systematic errors due to the current level of mismodeling in the multipolar zonal coefficients of the geopotential. It deals also with the linear combination approach which could be implemented in order to reduce the bias due to the latter ones. Section 3 summarizes our findings and offer our conclusions. Appendix A contains the tables and the figures. Appendix B displays the analytically calculated orbital precessions due to the zonal harmonics of the geopotential up to degree ℓ = 8. Appendix C deals in detail with the impact of the atmospheric drag on all the orbital elements of a spherical, passive geodetic satellite in a highly elliptical orbit.

Two different orbital configurations for HERO
In Table A2 and Table A6, two different orbital configurations are proposed. They imply highly eccentric orbits, characterized by values of the eccentricity as large as e = 0.45 and e = 0.82, respectively, and the critical inclination I crit = arcsin 2/ √ 5 which allows to keep the argument of perigee ω essentially fixed over any reasonable time span for an actual data analysis and the longitude of the ascending node Ω circulating at a relatively high pace. Their orbital periods are P b = 4.3 hr and P b = 21.3 hr, with perigee heights of h min = 1, 047 km and h min = 642 km, respectively.
One of the most interesting relativistic features of motion is, perhaps, the relatively large value of the expected semimajor axis increase ȧ induced by the pN gravitoelectric quadrupolar acceleration of Equation (1); let us recall that it is (Soffel et al. 1987;Brumberg 1991;Iorio 2013) da dt = 9 a n 3 b R 2 e J 2 e 2 6 + e 2 sin 2 I sin 2ω 8 c 2 1 − e 2 4 , where n b = µ/a 3 = 2π/P b is the Keplerian mean motion. Indeed, according to Table A3 and  Table A7, its predicted rates for the orbital geometries considered in Table A2 and Table A6 are ȧ = 3.8 cm yr −1 and ȧ = 11.6 cm yr −1 , respectively. Suffice it to say that for the existing passive satellites of the LAGEOS family, whose orbits are essentially circular, secular decay rates have been measured over the last decades with an experimental accuracy of σ ȧ ≃ 3 cm yr −1 for LAGEOS (Rubincam 1982;Sośnica et al. 2014;Sośnica 2014), and σ ȧ = 0.7 cm yr −1 for LARES (Lucchesi et al. 2019). It is arguable that an active mechanism of compensation of the non-gravitational accelerations would allow to increase such accuracies, allowing, perhaps, to measure the pN quadrupolar effect at a ≃ 1 − 10% level, depending on the orbital configuration adopted. As far as possible competing effects of gravitational origin are concerned, neither the static and time-dependent parts of the geopotential nor 3rd-body lunisolar attractions induce nonvanishing averaged perturbations on a. Thus, it is of the utmost importance the reduction of the non-conservative accelerations. Among them, a prominent role is played the atmospheric drag, treated in detail in Appendix C, because it causes a net long-term averaged decay rate of a. By modeling the spacecraft as a LARES-like cannonball geodetic satellite for the sake of simplicity, it can be shown that, in the case of the orbital configuration of Table A2, the average acceleration due to the neutral drag only amounts to A drag = (2.3 − 0.8) × 10 −11 m s −2 over one orbital period, so that the resulting effect on a would be as large as ȧ drag = −(5.1 − 2.0) m yr −1 ; see Table A5, and Figures A1 to A2. For the more eccentric orbital configuration of Table A6, we have A drag = (1.22 − 7.3) × 10 −11 m s −2 and ȧ drag = −(27.6 − 164.6) m yr −1 , as shown by Table A9, and Figures A3 to A4. An inspection of Figure A1 and Figure A3 reveals that, as expected, most of the disturbing effect occurs around the perigee passage, i.e. for f ≃ 0. This may help in suitably calibrating the counteracting action of the drag-free mechanism. Table A3 and Table A6 show that also all the other Keplerian orbital elements exhibit non-zero secular rates of change due to Equation (1). This is a potentially important feature since they could be linearly combined, as suggested by 1 Shapiro (1990) in a different context, in order to decouple the pN effect(s) of interest from the disturbing mismodeled Newtonian secular precessions induced by the Earth's zonal multipoles. Indeed, contrary to a, the other Keplerian orbital elements do exhibit long-term averaged precessions due to the classical zonal harmonics J ℓ , ℓ = 2, 3, 4, . . . of the geopotential; they are analytically calculated in Appendix B up to degree ℓ = 8. Depending on the specific orbital geometry, they can be both secular and harmonic, or entirely 2 secular. To this aim, we complete the set of the pN orbital effects by analytically calculating the averaged rate of the mean anomaly at epoch due to the pN accelerations considered. It turns out that Equations (1) to (2) induce η = µ n b R 2 J 2 − 80 + 73 e 2 (1 + 3 cos 2I) − 84 1 + 2 e 2 sin 2 I cos 2ω 32 c 2 a 3 1 − e 2 5/2 , η = 9 G S R 2 ε 2 5 cos 3I + cos I 3 + 10 sin 2 I cos 2ω 56 a 5 c 2 1 − e 2 2 , respectively, while the gravitoelectric mass monopole acceleration yields where ζ M m / (M + m) 2 , and m is the satellite's mass. Instead, it turns out that there is no net Lense-Thirring effect on η. To the benefit of the reader, we review the linear combination approach, which is a generalization of that proposed explicitly for the first time by Ciufolini et al. (1996) to test the pN Lense-Thirring effect in the gravitomagnetic field of the Earth with the artificial satellites of the LAGEOS family. It should be noted that, actually, it is quite general, being not necessarily limited just to the pN spin dipole case. By looking at N orbital elements 3 κ (i) experiencing classical long-term precessions due to the zonals of the geopotential, the following N linear combinations can be written down They involve the pN averaged precessions κ (i) pN as predicted by general relativity and scaled by a multiplicative parameter µ pN , and the errors in the computed secular node precessions due to the uncertainties in the first N − 1 zonals J s , s = 2, 3, . . . N, assumed as mismodeled through δJ s , s = 2, 3, . . . N. In the following, we will use the shorthanḋ for the partial derivative of the classical averaged precession κ J ℓ with respect to the generic even zonal J ℓ of degree ℓ; see Appendix B. Then, the N combinations of Equation (10) are posed equal to the experimental residuals δκ (i) , i = 1, 2, . . . N of each of the N orbital elements considered. In principle, such residuals account for the purposely unmodelled pN effect, the mismodelling of the static and time-varying parts of the geopotential, and the non-gravitational forces. Thus, one gets If we look at the pN scaling parameter 4 µ pN and the mismodeling in the first N − 1 zonals δJ s , s = 2, 3, . . . N as unknowns, we can interpret Equation (12) as an inhomogenous linear system of N algebraic equations in the N unknowns whose coefficients are while the constant terms are the N orbital residuals It turns out that, after some algebraic manipulations, the dimensionless pN scaling parameter, which is 1 in general relativity, can be expressed as In Equation (16), the combination of the N orbital residuals is, by construction, independent of the first N − 1 zonals, being impacted by the other ones of degree ℓ > N 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 zonal rates. Instead, combines the N pN orbital precessions as predicted by general relativity. The dimensionless coefficients c j , j = 1, 2, . . . N − 1 in Equation (17)-Equation (18) depend only on some of the orbital parameters of the satellite(s) involved in such a way that, by construction, C δ = 0 if Equation (17) is calculated by posing for any of the first N − 1 zonals, independently of the value assumed for its uncertainty δJ ℓ .
As far as HERO is concerned, the linear combination of the four experimental residuals δΩ, δη, δe, δω of the satellite's node, mean anomaly at epoch, eccentricity and perigee suitably designed to cancel out the secular precessions due to the first three zonal harmonics J 2 , J 3 , J 4 of the geopotential is The coefficients c 1 , c 2 , c 3 turn out to be Their numerical values, computed with the formulas of Appendix B for the orbital configurations of Table A2 and Table A6, are listed in Table A4 and Table A8, respectively. In them, the combined mismodeled classical precessions due to the uncancelled zonals, calculated by assuming the formal, statistical sigmas σ J ℓ , ℓ = 5, 6, . . . of the recent global gravity field solution Tongji-Grace02s (Chen et al. 2018) as a measure of their uncertainties δJ ℓ , ℓ = 5, 6, . . ., are reported as well. However, caution is in order since the realistic level of mismodeling in the geopotential's coefficients is usually larger than the mere formal errors released in the models produced by various institutions and publicly available on the Internet at http://icgem.gfz-potsdam.de/tom longtime. A correct evaluation of the actual uncertainties in the zonal harmonics require great care by suitably comparing different global gravity field models; we will not deal with such a task here. From an inspection of Table A4 and Table A8, it can be noted that the (formal) impact of the uncancelled zonals on the combined pN mass quadrupole effect (GJ 2 /c 2 ) is at the ≃ 0.6 − 0.07% level for the proposed orbital configurations of Table A2 and Table A6. If the pN spin dipole Lense-Thirring effect (GS /c 2 ) is considered, the systematic error due to the mismodeling in J ℓ , ℓ > 4 is about ≃ 0.1 − 0.03%. The pN mass monopole combined precessions (GM/c 2 ) are affected at the ≃ (20 − 5) × 10 −7 relative level. Instead, it turns out that the the combined mismodelled classical precessions are at the same level of the pN spin octupole trends (GS ε 2 /c 2 ). It may not be unrealistic to expect that, when the forthcoming global gravity field models based on the analysis of the entire long data records of the dedicated GRACE and GOCE missions will be finally available, the current merely formal level of uncertainties in the geopotential's zonal harmonics may be considered as realistic. Moreover, in the next years, the mission GRACE-FO (GFO) (Sheard et al. 2012), launched in May 2018, will also contribute to the production of new global Earth's gravity field models of increased quality.
On the other hand, the size of the coefficients c 1 , c 2 , c 3 amplifies the impact of any non-gravitational perturbations that may affect the spacecraft; thus, they should be effectively counteracted by some active drag-free apparatus. In particular, the coefficient c 2 of the eccentricity is ≃ 30 − 40; Table A5 and Table A9 show that the expected secular decrease rate of e due to the atmospheric drag is rather large. Thus, some trade-off may be required among the need of reducing the systematic error of gravitational origin and the actual performance of the drag-free mechanism by looking, e.g., at different linear combinations. It may be interesting to note the case of the mean anomaly at epoch. Indeed, in the case of the high perigee orbital configuration of Table A2, Table A3 and Table A5 tell us that the neutral atmospheric drag would represent just ≃ 1 − 2% of the predicted pN GJ 2 /c 2 precession on η. On the other hand, the present-day formal mismodeling in the classical J 2 -induced rate is about 19% of it. If the pN Schwarzschild-like effect is considered, the formal bias due to J 2 is at the ≃ 10 −5 level, while the impact of the atmospheric drag is as little as ≃ 10 −6 . The neutral atmospheric drag has a larger impact on the pN precessions of η in the case of the low perigee configuration of Table A6, as shown by Table A7  and Table A9.

Summary and overview
The HERO concept, meant as a hopefully drag-free spacecraft moving in a highly eccentric orbit in a frozen perigee configuration aimed to perform several tests of relativistic gravity in the Earth's spacetime, represents, in principle, a promising opportunity to measure, for the first time, a general relativistic effect which has never received the same attention of the more known Schwarzschild and Lense-Thirring precessions so far: the post-Newtonian gravitoelectric orbital shifts due to the mass quadrupole moment of the Earth. Indeed, the systematic uncertainty in the combined satellite's precessions due to the formal, statistical errors in the competing Newtonian mass multipoles of the geopotential, as per one of the most recent global gravity field models, is currently below the per cent level for both the orbital configurations proposed. A unique feature of such a post-Newtonian effect is that also the semimajor axis a undergoes a long-term variation which, for a frozen perigee configuration, resembles a secular trend of the order of ≃ 4 − 11 cm yr −1 , depending on the orbital geometry chosen. At present, the secular decay of the semimajor axis of the existing passive geodetic satellite LARES has been measured to an accuracy better than 1 cm yr −1 at 2σ level. As far as the traditional Lense-Thirring and Schwarzschild-like post-Newtonian precessions, the formal systematic bias due to the present-day mismodeling in the classical Earth's zonal harmonics is currently 0.1% and 0.0002%, respectively if a suitable linear combination of some of the orbital elements of HERO is adopted. However, it must be stressed that the actual uncertainties in the zonal multipoles of the terrestrial gravity field may usually be (much) worse than the sigmas released in the various global gravity solutions. Nonetheless, it cannot be ruled out that, if and when HERO will fly, our knowledge of the Earth's gravity field will have reached such levels that today's only formal uncertainties can finally be considered as truly realistic. In addition to the post-Newtonian accelerations, HERO may perform an accurate test of the gravitational red-shift in view of its high eccentricity. Also several models of modified gravity, which generally affect the perigee and the mean anomaly at epoch with secular precessions, could be fruitfully put to the test. A crucial aspect is represented by the level of compensation of the non-gravitational perturbation which will be practically attainable with some drag-free apparatus; suffice it to say that the nominal size of the competing secular decrease of the semimajor axis due to, e.g., the neutral atmospheric drag can reach the ≃ 5 −160 m yr −1 level if a passive, cannonball satellite is considered. Its investigation deserves a dedicated publication. Table A1: Relevant physical parameters of the Earth (Sehnal 1981;Petit, Luzum & et al. 2010;Durand-Manterola 2009). The zonal harmonics of the geopotential of degree ℓ are given by J ℓ = − √ 2ℓ + 1 C ℓ,0 , ℓ = 2, 3, 4, . . ., where C ℓ,0 , ℓ = 2, 3, 4, . . . are the fully normalized Stokes coefficients of degree ℓ and order m = 0 of the multipolar expansion of the Newtonian part of the Earth's gravity field. The formal, statistical errors in the first seven Stokes coefficients of the geopotential, along with their nominal values, were retrieved from the global gravity field solution Tongji-Grace02s (Chen et al. 2018 9.571989759740 × 10 −7 -Normalized Stokes coefficient C 4,0 5.399893295930 × 10 −7 -Normalized Stokes coefficient C 5,0 6.86499810446677 × 10 −8 -Normalized Stokes coefficient C 6,0 −1.49976729587105 × 10 −7 -Normalized Stokes coefficient C 7,0 9.05017773295824 × 10 −8 -Normalized Stokes coefficient C 8,0 4.94794369681244 × 10 −8 -Formal error σ C 2,0 2.98340899705584 × 10 −13 -Formal error σ C 3,0 8.39284383652709 × 10 −14 -Formal error σ C 4,0 4.07426781903578 × 10 −14 -Formal error σ C 5,0 2.57688174349872 × 10 −14 -Formal error σ C 6,0 1.89009491873398 × 10 −14 -Formal error σ C 7,0 1.50081719867797 × 10 −14 -Formal error σ C 8,0 1.27528335995664 × 10 −14 - Table A2: Orbital configuration of the proposed satellite HERO: high perigee case. The orbital motion is rather fast since the orbital period P b is as short as ≃ 4 hr. The period P ω of the perigee is mainly determined by its secular precession due to J 3 , J 4 because of the critical inclination which makes the secular precession due to J 2 nominally vanishing. Note the relatively short period P Ω of the node, amounting to less than 2 yr. 3.7 × 10 −10 − Table A3: Nominal pN (first four rows from the top) and mismodeled Newtonian (first seven rows from the bottom) rates of change, averaged over one orbital revolution, of the semimajor axis a, the eccentricity e, the inclination I, the longitude of the ascending node Ω, the argument of pericenter ω, and the mean anomaly at epoch η for the ideal (no orbital injection error on I assumed) orbital configuration of Table A2. The units are cm yr −1 for ȧ , and mas yr −1 for ė , İ , Ω , ω , η . The uncertainties in the classical rates of change due to the geopotential are the formal, statistical errors σ J ℓ in J ℓ , ℓ = 2, 3, . . . 8 of the model Tongji-Grace02s (Chen et al. 2018) quoted in Table A1.   (20) canceling out the classical precessions induced by the first three zonal harmonics J 2 , J 3 , J 4 for the orbital configuration of Table A2. Middle rows: Uncancelled mismodeled precessions due to the zonal harmonics J 5 , J 6 , J 7 , J 8 , in mas yr −1 , linearly combined according to Equation (20); the formal, statistical errors σ J ℓ in J ℓ , ℓ = 5, 6, 7, 8 of the model Tongji-Grace02s (Chen et al. 2018), quoted in Table A1 Table A5: Numerically integrated nominal rates of change, averaged over one orbital revolution, of the semimajor axis a, the eccentricity e, the inclination I, the longitude of the ascending node Ω, the argument of pericenter ω, and the mean anomaly at epoch η induced by the neutral atmospheric drag for the orbital configuration of Table A2. The units are m yr −1 for ȧ , and mas yr −1 for ė , İ , Ω , ω , η . For the satellite, assumed spherical in shape and passive, we adopted C D = 3.5, Σ = 2.69 × 10 −4 as for the existing LARES (Pardini et al. 2017). In regard to the Earth's neutral atmospheric density, we adopted r 0 = r min = a (1 − e) = 1, 046.86 km, ρ 0 = ρ max = (7.3 − 2.8) × 10 −15 kg m −3 , ρ min = 0.001 ρ L , λ = 872.87 km − 938.49 km; the neutral atmospheric density at the height of LAGEOS is ρ L = 6.579 × 10 −18 kg m −3 (Lucchesi et al. 2015). See Appendix C for details. Neither approximations in e nor in ν Ψ/n b were used. The value of ρ 0 = ρ max was kept fixed over one orbital revolution. Cfr. with the analytical plots in Figure A1 and the numerically produced time series in Figure A2.  Table A6: Orbital configuration of the proposed satellite HERO: low perigee case. The orbital motion is relatively slow since the orbital period P b is as long as more than 21 hr. The period P ω of the perigee is mainly determined by its secular precession due to J 3 , J 4 because of the critical inclination which makes the secular precession due to J 2 nominally vanishing. The period P Ω of the node is rather long, amounting to more than 13 yr.  Table A7: Nominal pN (first four rows from the top) and mismodeled Newtonian (first seven rows from the bottom) rates of change, averaged over one orbital revolution, of the semimajor axis a, the eccentricity e, the inclination I, the longitude of the ascending node Ω, the argument of pericenter ω, and the mean anomaly at epoch η for the ideal (no orbital injection error on I assumed) orbital configuration of Table A6. The units are cm yr −1 for ȧ , and mas yr −1 for ė , İ , Ω , ω , η . The uncertainties in the classical rates of change due to the geopotential are the formal, statistical errors σ J ℓ in J ℓ , ℓ = 2, 3, . . . 8 of the model Tongji-Grace02s (Chen et al. 2018) quoted in Table A1.   (20) canceling out the classical precessions induced by the first three zonal harmonics J 2 , J 3 , J 4 for the orbital configuration of Table A6. Middle four rows: Uncancelled mismodeled precessions due to the zonal harmonics J 5 , J 6 , J 7 , J 8 , in mas yr −1 , linearly combined according to Equation (20); the formal, statistical errors σ J ℓ in J ℓ , ℓ = 5, 6, 7, 8 of the model Tongji-Grace02s (Chen et al. 2018), quoted in Table A1, were used. Lower four rows: pN precessions, in mas yr −1 , linearly combined according to Equation (20).

Appendix A Tables and figures
mas yr −1 M c −2 6, 033.6 mas yr −1 Table A9: Numerically integrated nominal rates of change, averaged over one orbital revolution, of the semimajor axis a, the eccentricity e, the inclination I, the longitude of the ascending node Ω, the argument of pericenter ω, and the mean anomaly at epoch η induced by the neutral atmospheric drag for the orbital configuration of Table A6. The units are m yr −1 for ȧ , and mas yr −1 for ė , İ , Ω , ω , η . For the satellite, assumed spherical in shape and passive, we adopted C D = 3.5, Σ = 2.69 × 10 −4 as for the existing LARES (Pardini et al. 2017). In regard to the Earth's neutral atmospheric density, we adopted r 0 = r min = a (1 − e) = 641.86 km, ρ 0 = ρ max = (6.9 − 1.11) × 10 −14 kg m −3 , ρ min = 0.0001 ρ L , λ = 3, 463.23 km − 3, 843.48 km; the neutral atmospheric density at the height of LAGEOS is ρ L = 6.579 × 10 −18 kg m −3 (Lucchesi et al. 2015). See Appendix C for details. Neither approximations in e nor in ν Ψ/n b were used. The value of ρ 0 = ρ max was kept fixed over one orbital revolution. Cfr. with the analytical plots in Figure A3 and the numerically produced time series in Figure A4.  Table A2. For the satellite, assumed spherical in shape and passive, we adopted C D = 3.5, Σ = 2.69 × 10 −4 as for the existing LARES (Pardini et al. 2017). In regard to the Earth's atmospheric density, we adopted r 0 = r min = a (1 − e) = 1, 046.86 km, ρ 0 = ρ max = 7.3 × 10 −15 kg m −3 , λ = 872.87 km. Neither approximations in e nor in ν Ψ/n b were used. The areas of the regions delimited by the curves and the f axis are the rates of change of the orbital elements averaged over one orbital period P b ; they are numerically calculated and displayed in Table A5. The value of ρ 0 = ρ max was kept fixed over one orbital revolution.  Table A2. For the satellite, assumed spherical in shape and passive, we adopted C D = 3.5, Σ = 2.69 × 10 −4 as for the existing LARES (Pardini et al. 2017). In regard to the Earth's atmospheric density, we adopted r 0 = r min = a (1 − e) = 1, 046.86 km, ρ 0 = ρ max = 7.3×10 −15 kg m −3 , λ = 872.87 km. Cfr. with the semi-analytical results of Figure A1 and Table A5.  Table A6. For the satellite, assumed spherical in shape and passive, we adopted C D = 3.5, Σ = 2.69 × 10 −4 as for the existing LARES (Pardini et al. 2017). In regard to the Earth's atmospheric density, we adopted r 0 = r min = a (1 − e) = 641.86 km, ρ 0 = ρ max = 6.9 × 10 −14 kg m −3 , λ = 3, 463.23 km. Neither approximations in e nor in ν Ψ/n b were used. The areas of the regions delimited by the curves and the f axis are the rates of change of the orbital elements averaged over one orbital period P b ; they are numerically calculated and displayed in Table A9. The value of ρ 0 = ρ max was kept fixed over one orbital revolution.  Table A6. For the satellite, assumed spherical in shape and passive, we adopted C D = 3.5, Σ = 2.69 × 10 −4 as for the existing LARES (Pardini et al. 2017). In regard to the Earth's atmospheric density, we adopted r 0 = r min = a (1 − e) = 641.86 km, ρ 0 = ρ max = 6.9 × 10 −14 kg m −3 , λ = 3, 463.23 km. Cfr. with the semi-analytical results of Figure A3 and Table A9.
Appendix B Classical long-term rates of change of the Keplerian orbital elements up to degree ℓ = 8 Here, we will analytically work out the coefficientsκ .ℓ ∂ κ J ℓ /∂J ℓ , κ = e, I, Ω, ω, η of the long-term rates of change of the eccentricity e, the inclination I, the longitude of the ascending node Ω, the argument of perigee ω, and the mean anomaly at epoch η induced by the first seven zonal harmonics of the geopotential up to degree ℓ = 8; the averaged rates of change of the semimajor axis a are all vanishing. We will adopt the Lagrange planetary equations (Bertotti, Farinella & Vokrouhlický 2003;Capderou 2005;Xu 2008;Kopeikin, Efroimsky & Kaplan 2011) applied to the perturbing potential of degree ℓ where P ℓ (ξ) is the Legendre polynomial of degree ℓ, averaged over one orbital period as disturbing function. We will not make any a-priory assumption on the orbital configuration of the satellite. As far as the orientation of the Earth's spin axis, we will align it to the reference z axis of an equatorial coordinate system. The following formulas include both the genuine secular and the harmonic components with the frequency of the perigee and its integer multiples.
In the case of the orbital configuration of Table A2, the perigee height is h min = 1, 046.86 km; we will determine the corresponding neutral atmospheric density ρ 0 = ρ max as follows. The neutral atmospheric density at the LARES height, which is h LR = 1, 450 km, amounts to ρ LR = 5.644 × 10 −16 kg m −3 (Pardini et al. 2017). According to Brito, Celestino & Moraes (2015), the neutral atmospheric density at h 700 = 700 km is 6.9 × 10 −14 kg m −3 (TD-88 model), or 1.11 × 10 −14 kg m −3 (NASA model). Thus, it is possible to calculate the characteristic scale length λ valid for the range 700 km < h < 1, 450 km as depending on the value adopted for ρ at 700 km. Since for the orbital configuration of HERO of Table A2 it is h 700 < h min < h LR , one can determine ρ max by using just Equation (C8) in depending on the value of λ 700/LR adopted. As far as ρ min is concerned, since h max = 13, 169.9 km is much larger than the height of, say, the LAGEOS satellite (h L = 5, 891.87 km), for which it is ρ L = 6.579 × 10 −18 kg m −3 (Lucchesi et al. 2015), it does not seem unreasonable to assume or so. Thus, Equation (C5), applied to ρ max , given by Equation (C9), and ρ min , given by Equation (C10), allows to infer the value for λ which must be used in the calculation of the drag effect for HERO. It is λ = 872.87 km − 938.49 km, depending on the value of ρ max adopted. As far as the more eccentric orbital configuration of Table A6 is concerned, by adopting for ρ 0 the two possible values of ρ 700 and, say, ρ min = 0.0001 ρ L , it turns out that λ = 3, 463.23 km − 3, 843.48 km.
Actually, even the density at a given height may not be regarded as truly constant because of a variety of geophysical phenomena characterized by quite different time scales. Anyway, in order to have an order-of-magnitude evaluation of the perturbing action of Equation (C1) on the motion of HERO, we calculate the averaged rates of change of its Keplerian orbital elements by keeping ρ 0 = ρ max in Equation (C3) fixed during one orbital period P b ; given its short duration, at least in the case of Table A2, it is not an unreasonable assumption. The large value of the eccentricity makes most of the existing results in the literature unsuitable to the present case; moreover, an exact analytical calculation without recurring to any approximation in both e and ν Ψ/n b is difficult. Thus, we follow two complementary approaches. In one of them, we, first, plot in Figure A1 and Figure A3 the analytical expressions of the ratios of the right-hand-sides of the Gauss perturbing equations for the rates of change of the Keplerian orbital elements, evaluated onto the unperturbed Keplerian ellipse, to P b as functions of the true anomaly f . In the most general case, by assuming the co-rotation of the atmosphere, they are n b 2 π da d f = − C D Σ ρ a n b √ 1 − e 2 V 1 + 2 e cos f + e 2 − ν 1 − e 2 3/2 cos I 2 π (1 + e cos f ) 2 , (C13) 4 π e (1 + e cos f ) 4 2 + 3 e 2 + 2 e 2 + e 2 cos f + e 2 cos 2 f − −ν 1 − e 2 3/2 (2 + e cos f ) cos I , where V 2 = 1 − ν 2 1 − e 2 3/2 cos I 1 + e 2 + 2 e cos f + ν 2 1 − e 2 3 3 + cos 2I + 2 sin 2 I cos 2u 4 (1 + e cos f ) 2 1 + e 2 + 2 e cos f , and u ω + f (C20) is the argument of latitude. It is worthwhile noticing that, in general, being even possible that V 2 − 1 1 (C22) for some values of f , thus preventing from expanding it in powers of ν. Then, we numerically calculate the areas under the resulting curves, i.e. we numerically integrate the averaged rates of change of the orbital elements for the given orbital configuration: see Table A5 and Table A9.
In the second approach, we numerically integrate the equations of motion of the satellite in rectangular Cartesian coordinates, referred to a geocentric equatorial coordinate system, with and without Equation (C1) over 1 yr; both the runs share the same initial conditions. Then, we subtract the resulting time series of the orbital elements in order to single out their shifts due to the disturbing acceleration. Finally, we perform a linear fit of them, and look at their slopes; we plot the fitted trends as functions of time t in Figure A2 and Figure A4. Both the methods reciprocally agree well, as shown by Figures A1 to A2 and Figures A3 to A4 for the neutral drag. A slight reduction turns out to occur if the decrease of the atmospheric co-rotation with height is taken into account as modeled by Membrado & Pacheco (2010) by assuming the simpler case of a constant viscosity. We tested our approach by checking that it was able to reproduce the observed features of the semimajor axis decay of LARES recently determined in Pardini et al. (2017). However, it must be stressed that such findings should be deemed just as indicative of the limitations of the scenario considered if a passive spacecraft were to be adopted. Indeed, they were computed preliminarily by assuming the same physical properties of the existing LARES satellite which, in principle, could well be superseded by a new, specifically manufactured spacecraft able to reduce both the drag coefficient C D and the area-to-mass ratio Σ. Moreover, also the actual temporal variability of the atmospheric density over timescales larger than the satellite's orbital period P b should be taken into account, especially if data will be collected during temporal intervals several years long. To this aim, it is important to note that an inspection of Equations (C13) to (C18) shows that no other sources of long-term modulation are present. Indeed, the circulating node Ω does not enter them, contrary to the perigee ω which, however, is held fixed by the adopted value of the inclination I.