Is there still something left that Gravity Probe B can measure?

We perform a full analytical and numerical treatment, to the first post-Newtonian (1pN) order, of the general relativistic long-term spin precession of an orbiting gyroscope due to the mass quadrupole moment $J_2$ of its primary without any restriction on either the gyro's orbital configuration and the orientation in space of the symmetry axis $\boldsymbol{\hat{k}}$ of the central body. We apply our results to the past spaceborne Gravity Probe B (GP-B) mission by finding a secular rate of its spin's declination $\delta$ which may be as large as $\lesssim 30-40\,\mathrm{milliarcseconds\,per\,year\,(\mathrm{mas\,yr}^{-1}})$, depending on the initial orbital phase $f_0$. Both our analytical calculation and our simultaneous integration of the equations for the parallel transport of the spin 4-vector and of the geodesic equations of motion of the gyroscope confirm such a finding. For GP-B, the reported mean error in measuring the spin's declination rate amounts to $\sigma^\mathrm{GP-B}_{\dot\delta}=18.3\,\mathrm{mas\,yr}^{-1}$. We also calculate the general analytical expressions of the gravitomagnetic spin precession induced by the primary's angular momentum $\boldsymbol J$. In view of their generality, our results can be extended also to other astronomical and astrophysical scenarios of interest like, e.g., stars orbiting galactic supermassive black holes, exoplanets close to their parent stars, tight binaries hosting compact stellar corpses.


Abstract
We perform a full analytical and numerical treatment, to the first post-Newtonian (1pN) order, of the general relativistic long-term spin precession of an orbiting gyroscope due to the mass quadrupole moment J 2 of its primary without any restriction on either the gyro's orbital configuration and the orientation in space of the symmetry axisk of the central body. We apply our results to the past spaceborne Gravity Probe B (GP-B) mission by finding a secular rate of its spin's declination δ which may be as large as 30 − 40 milliarcseconds per year mas yr −1 , depending on the initial orbital phase f 0 . Both our analytical calculation and our simultaneous integration of the equations for the parallel transport of the spin 4-vector S and of the geodesic equations of motion of the gyroscope confirm such a finding. For GP-B, the reported mean error in measuring the spin's declination rate amounts to σ GP−Ḃ δ = 18.3 mas yr −1 . We also calculate the general analytical expressions of the gravitomagnetic spin precession induced by the primary's angular momentum J. In view of their generality, our results can be extended also to other astronomical and astrophysical scenarios of interest like, e.g., stars orbiting galactic supermassive black holes, exoplanets close to their parent stars, tight binaries hosting compact stellar corpses.
keywords general relativity and gravitation; experimental studies of gravity; experimental tests of gravitational theories; satellite orbits; harmonics of the gravity potential field

Introduction
To the first post-Newtonian (1pN) level of order O c −2 , where c is the speed of light in vacuum, the geodesic motion of a test particle through the deformed spacetime outside an axially symmetric oblate body of mass M, equatorial radius R and dimensionless mass quadrupole moment J 2 is characterized by certain secular orbital precessions (Soffel et al. 1987;Soffel 1989;Brumberg 1991;Iorio 2015). They have recently gained attention, being possibly detectable in some proposed space-based experiments like, e.g., HERO (Iorio 2019).
Here, we will look at the long-term 1pN rate of change, proportional to J 2 c −2 , of the spin S of a pointlike gyroscope freely moving with velocity v around an oblate primary. The analogous 1pN gyro's precessional effects due to only the mass monopole (the mass M) and the spin dipole (the proper angular momentum J) moments of the central body acting as source of the gravitational field are the time-honored de Sitter-Fokker (or geodetic) (de Sitter 1916;Fokker 1920) and Pugh-Schiff (Pugh 1959;Schiff 1960) precessions, respectively. They were recently measured by the spaceborne mission Gravity Probe B (GP-B) in the field of Earth to ≃ 0.3% and ≃ 19%, respectively, (Everitt et al. 2011(Everitt et al. , 2015, despite a higher accuracy had been originally expected (Everitt 1974;Everitt et al. 2001). We will not restrict ourselves to any particular orbital configuration of the moving gyroscope, and the symmetry axis of the oblate primary will retain an arbitrary orientation in space. We will calculate the sought effect both numerically and analytically by finding, among other things, that it depends on the initial position of the gyro along its orbit. In the case of GP-B, it turns out that the rate of change of the spin's declination (DEC) δ, averaged over an orbital revolution, may be as large as ≃ 30 − 40 milliarcseconds per year mas yr −1 . Thus, it may be potentially measurable in a future data reanalysis since the reported average experimental accuracy in measuring the temporal evolution of δ is (Everitt et al. 2011(Everitt et al. , 2015 σ GP−Ḃ δ = 18.3 mas yr −1 . For previous analytical calculations, relying upon various simplifying assumptions concerning the gyro's orbit and different computational approaches, see O'Connell (1969); Barker & O'Connell (1970);Breakwell (1988); Adler & Silbergleit (2003). Even putting aside the issue of the particular orbital configurations adopted, they are, at least, incomplete since they neglect an important feature in the averaging procedure yielding to the dependence on the gyro's initial conditions which, instead, we will take into account. Our simultaneous numerical integrations of the equations of motion of the gyro and of its spin will display it, by supporting our analytical findings. Moreover, it seems that the aforementioned works return incorrect results even for the part which is independent of the initial conditions, being also in mutual disagreement. In the following, we will not deal too much with the spin's right ascension (RA) α since it turns out that, for GP-B, its total rate of change of the order of O J 2 c −2 is negligible.
For the sake of completeness, we will analytically derive also the generalization of the Pugh-Schiff gravitomagnetic spin precession valid for an arbitrary orientation of the primary's angular momentum J and for a generic orbital configuration of the gyroscope.
The generality of our approach allows our results to be extended also to other astronomical and astrophysical scenarios of interest like, e.g., other planets of our solar system, exoplanets, binaries with compact stellar corpses, supermassive black holes orbited by planets and stars. To this aim, it may be interesting to recall that Haas & Ross (1975) investigated the possibility of using spacecraft-based missions to measure the angular momenta of Jupiter an the Sun by means of the gravitomagnetic Pugh-Schiff spin precession.
The outline of the paper is as follows. In Section 2, we numerically calculate the total spin precession of the order of O J 2 c −2 by simultaneously integrating the equations for the parallel transport of the gyro's spin 4-vector and the geodesic equations of motion of the gyroscope. The spin and orbital configurations of GP-B are used. Section 3 is devoted to the analytical calculation. It, first, includes the direct effects (Section 3.1), which are the de Sitter precession for an arbitrary orbital configuration (Section 3.1.1), and the component of the spin rate of change of the order of O J 2 c −2 arising from using a fixed Keplerian ellipse for the orbital average (Section 3.1.2). Then, in Section 3.2, we calculate the indirect, or mixed, components of the sought precession. They are those arising from averaging the instantaneous 1pN de Sitter-like spin rate over the orbital period of a J 2 -driven precessing ellipse (Section 3.2.1), and those coming from the inclusion of the instantaneous orbital shifts caused by J 2 in the averaging procedure (Section 3.2.2). The total analytical spin precession of the order of O J 2 c −2 is discussed in Section 3.3, where the GP-B case is illustrated and compared with the numerical results of Section 2. The general expression of the gravitomagnetic spin precession is analytically calculated in Section 4. Section 5 summarizes our finding and offers our conclusions.
2. Numerical simulations: simultaneously integrating the equations for the motion of the gyroscope and of its spin The equations for the parallel transport of the spin 4-vector S of a pointlike gyroscope freely moving in the deformed spacetime of a central body are (Zee 2013;Ohanian & Ruffini 2013;Misner, Thorne & Wheeler 2017;Will 2018) where τ is the gyro's proper time, are the spacetime's Christoffel symbols, g γλ , g γλ , γ, λ = 0, 1, 2, 3 are the components of the spacetime metric tensor and of its inverse, respectively, and dx β /dτ, β = 0, 1, 2, 3 are the components of the gyro's 4-velocity u. The space-like components of S are the components of the gyro's spin vector S, i.e. S i = S i , i = 1, 2, 3. The time-like component S 0 of S is determined by the constraint 1 The geodesic equations of motion of the pointlike gyroscope are In standard pN isotropic coordinates, the components of the metric tensor of the spacetime of an isolated body are (Soffel 1989;Poisson & Will 2014) where is the Kronecker delta. In the following, we will use cartesian coordinates, so that x 1 = x, x 2 = y, x 3 = z, S 1 = S x , S 2 = S y , S 3 = S z . In Equations (5)-(7), the potential U (r) of the oblate mass is In Equation (9), µ GM is the gravitational parameter of the central body, G is the Newtonian constant of gravitation, is the Legendre polynomial of degree 2, while is the cosine of the angle between the body's symmetry axisk and the unit position vectorr. In the case of a diagonal metric, as for Equations (5)-(7), Equation (3) yields We set up a numerical code to simultaneously integrate both Equation (1) and Equation (4) for an arbitrary orientation ofk in space and unrestricted orbital configurations for the moving gyroscope. The space-like components of S are parameterized in terms of two spherical angles α, δ as which, in the case of Earth and an equatorial coordinate system, are the spin's right ascension and declination, respectively. As initial conditions for both the gyroscope orbit and its spin, we adopt those of GP-B (Kahn 2007), summarized in Table 1. As far as the initial value of S 0 is concerned, it can be retrieved from the condition of Equation (3). The initial values of the space-like components of the 4-velocity u can be obtained from where v i , i = 1, 2, 3 are the components of the velocity v (see Equation (46)), and We, first, test our routine by successfully reproducing the de Sitter precession, shown in Figure 1. The time series in it were obtained by switching off J 2 in both Equation (1) and Equation (4). They correspond to the orbital average over a Keplerian ellipse 2 of the 1pN components of the right-hand-sides of Equation (1) for i = 1, 2, 3 and J 2 = 0. As expected, all the signatures in Figure 1 are independent of f 0 . Figure 2 displays the "direct" part of the spin precession of the order of O J 2 c −2 obtained by restoring J 2 in Equation (1), but not in Equation (4), and subtracting from the resulting signatures the purely de Sitter ones. It essentially corresponds to the orbital average of the 1pN components of the right-hand-sides of Equation (1) for i = 1, 2, 3 and J 2 0 over an actually non-existent Keplerian ellipse 3 . Clearly, it is an unphysical situation which is just an intermediate check of our analytical calculation, to be displayed in Section 3.1.2, and of the results existing in the literature. Its slope amounts to 5.8 mas yr −1 , and is independent of f 0 . As we will see in Section 3.1.2, our  (1) and Equation (4) with J 2 = 0 in Equations (5)-(9), and calculating arcsin S z (τ) for the resulting solution S z (τ) of each run. It turned out that essentially the same outcome can also be obtained by replacing Equation (4) with the 3-dimensional Newtonian acceleration A N = − µ/r 2 r for the gyroscope. The initial conditions adopted, common to all of the integrations, were those of GP-B (Kahn 2007), summarized in Table1. All the shifts are independent of f 0 , and their slopes amount just toδ dS = −6603.8 mas yr −1 , as expected.
analytical outcome for the direct precession of the order of O J 2 c −2 agrees with Figure 2 to within ≃ 0.6 mas yr −1 . Instead, the part of Equation (53) of Barker & O'Connell (1970) containing J 2 allows to obtain∆δ = 4.1 mas yr −1 , while the J 2 -dependent part of Ω G in Adler & Silbergleit (2003, pag. 153) corresponds to∆δ = 6.6 mas yr −1 . As it will be demonstrated in Section 3.1.2, both of them disagree with our analytical calculation.
The total spin precessions of the order of O J 2 c −2 , obtained by simultaneously integrating both Equation (1) and Equation (4) with J 2 0 for different values of f 0 and subtracting the purely de Sitter trends from the resulting signatures, are displayed in Figure 3. They can be thought as the sum of the direct precession of Figure 2 and of the "indirect", or "mixed", ones arising from the fact that, in this case, the trajectory of the gyroscope is, more realistically, a (slowly) Fig. 2.-Numerically produced direct 1pN J 2 -induced yearly shifts ∆δ (τ), in mas, of the declination δ of the spin axis of a gyroscope orbiting the Earth along a fictitious Keplerian ellipse for some different initial values of the true anomaly f 0 . Each of the times series ∆δ (τ) was obtained by simultaneously integrating Equation (1) with J 2 0 and Equation (4) with J 2 = 0, and subtracting from each of them the corresponding time series obtained by integrating both Equation (1) and Equation (4) with J 2 = 0 (the de Sitter trends). Then, arcsin S z (τ) was calculated for each of the resulting solutions S z (τ). The initial conditions adopted, common to all of the runs, were those of GP-B (Kahn 2007), summarized in Table 1. All the shifts are independent of f 0 , and agree with the result calculated analytically in Section 3.1.2. The slope amounts to∆δ = 5.8 mas yr −1 . precessing ellipse mainly driven by 4 J 2 . It can be thought as if, in addition to the Keplerian average of the J 2 -dependent parts of the space-like components of Equation (1), the de Sitter-like 1pN components of the right-hand-sides of Equation (1) for i = 1, 2, 3 and J 2 = 0 were averaged over one orbital revolution by taking now into account also the J 2 -induced instantaneous changes of the osculating Keplerian orbital elements parameterizing the varying ellipse, and the fact that the orbital period is the time interval between two successive passages at a changing perigee. The same, in principle, would hold also for the 1pN orbital changes which, however, would affect the spin precession to the 1/c 4 level. Effects of the order of O J 2 2 c −2 would arise by repeating the same average for the 1pN components of the right-hand-sides of Equation (1) for i = 1, 2, 3 and J 2 0. Our numerical integration accounts simultaneously for all such negligible effects of higher order as well. A striking feature of Figure 3 is that the indirect effects induce a neat  (1) and Equation (4), both with J 2 0 in Equations (5)-(9) and subtracting the corresponding de Sitter trends from each of them, and calculating arcsin S z (τ) for the resulting solution S z (τ) of each run. The initial conditions adopted, common to all of the integrations, were those of GP-B (Kahn 2007), summarized in Table 1. All the numerically integrated shifts agree with those calculated analytically in Section 3.3 to within 5 − 8 mas yr −1 . Such a discrepancy is not statistically significative since it is smaller than σ GP−Ḃ δ = 18.3 mas yr −1 (Everitt et al. 2011(Everitt et al. , 2015. dependence on f 0 which can yield spin precessions as large as ≃ 30 − 40 mas yr −1 . It is a quite important finding since the reported mean error in measuring the spin's declination precession of GP-B is σ GP−Ḃ δ = 18.3 mas yr −1 (Everitt et al. 2011(Everitt et al. , 2015, and it may prompt some reanalysis of the mission data. Such a dependence on f 0 induced by the mixed effects is captured and reproduced by our analytical calculation of the overall precession in Sections 3.1.2 to 3.2.2 to within 5 − 8 mas yr −1 ; cfr. with Figure 4 in Section 3.3. Instead, it is missing in the literature. Indeed, if, on the one hand, Adler & Silbergleit (2003) seemingly dealt only with the direct J 2 -induced precession, on the other hand, Barker & O'Connell (1970) were aware of such an issue, but they somehow treated it only partly since their Equation (52) does not contain any dependence on the initial orbital phase. Should it ever be related to the aforementioned issue of the orbital period in a precessing orbit, it is in disagreement with our analytical results for it, as we will show in Section 3.2.1.

Analytical calculation
calculated with Equation (9) in Equations (5)-(7), to the order of O c −2 , one obtains the instantaneous rates of change of the gyro's spin components as where the coefficients of the matrices T dS , T J 2 are, in general, time-dependent. They are and As far as the rates of change of the spin's spherical angles α, δ are concerned, from Equations (13)-(15) one gets Since we are interested in the long-term rate of change of S, we must properly average the right-hand-sides of Equations (19)-(21) over one orbital period P b . It requires care, especially for the effects of the order of O J 2 c −2 . Indeed, the actual orbital path of the gyroscope around its distorted primary is a generally slowly precessing ellipse (Capderou 2005), not a fixed Keplerian one as it would be if it were 5 J 2 = 0. This implies that, during an orbital revolution, all the Keplerian orbital elements characterizing the shape, the size and the orientation of the ellipse undergo instantaneous variations due to J 2 which should be taken into account in the averaging procedure since they give rise to effects which are just of the order of 6 O J 2 c 2 . Moreover, the fact that the line of the apsides, from which the time-dependent true anomaly 7 f is reckoned, does vary during the orbital motion because of J 2 has to be taken into account as well, yielding further contributions of the order of O J 2 c 2 . Such "indirect", or "mixed", features are to be added to the direct ones arising from a straightforward average of Equations (31)-(39) over an unperturbed Keplerian ellipse assumed as reference trajectory.
From a computational point of view, we can split the calculation of the averaged 1pN gyro's spin precession in two parts.

The direct effects
The first one deals with what one may define as the "direct" effects, denoted in the following with the superscript "dir", arising from averaging Equations (22) an unchanging 8 Keplerian ellipse. The latter is characterized by (Brumberg 1991) In Equations (45)- (46), it isP =l cos ω +m sin ω, In Equations (42)-(50), p, a , e , I , Ω , ω are the semilatus rectum, the semimajor axis, the eccentricity, the inclination, the longitude of the ascending node, and the argument of pericentre, respectively, of the Keplerian ellipse. The size and the shape of the latter are fixed by a and e, respectively. The inclination and the position of the orbital plane with respect to the reference {x, y} plane are determined by I and Ω, respectively; the line of the nodes is the intersection of the orbital plane with the reference {x, y} plane. The orientation of the ellipse within the orbital plane itself is characterized by ω. The unit vectorl is directed along the line of the nodes toward the ascending node, whilem lies in the orbital plane perpendicularly tol. The unit vectorP is directed along the line of the apsides toward the pericentre in the orbital plane whereQ stays transversely toP itself. Finally, we mention also the unit vector directed along the orbital angular momentum perpendicularly to the orbital plane 9 .
The resulting direct effects consist of the usual de Sitter precession, and of one part of the 1pN spin's rate of change due to J 2 . In Section 3.1.1 and Section 3.1.2, we will display the explicit expressions of the averaged matrix elements of Equations (22)-(39). For the sake of simplicity, we will omit the brackets . . . denoting the average over one orbital period throughout the paper.

The de Sitter precession
Let us introduce the following dimensional amplitude having the dimension of reciprocal time where R s 2 µ/c 2 is the primary's Schwarzschild radius. The analytical expressions of the average of Equations (22)-(30) yield the geodetic precession for an arbitrary orbital configuration of the moving gyroscope. We have T dS T dS T dS 9 It turns out thatl,m,ĥ are a right-handed triad of unit vectors.
T dS T dS From Equations (40)-(41) and Equations (53)-(61), it is possible to obtain  (61) show that the 1pN spin rate due to the mass monopole of the primary can be written as with The vectorial expression of Equation (64) agrees with, e.g., (10.146a) of Poisson & Will (2014) in the limit e → 0.

The indirect effects
This part treats what one may call the "indirect", or "mixed", effects arising from the precession of the orbit of the gyro caused by the oblateness of the primary. When applied to Equations (22)-(30), they give rise to further components of the gyro's spin rate of change of the order of J 2 c −2 which are to be added to the direct ones of Section 3.1.2 in order to have the total expression of the 1pN spin rate due to J 2 . In turn, the calculation of the mixed effects can be split into two parts.
The first one, tagged in the following with the superscript "mix I", consists of averaging Equations (22)-(30), to be evaluated onto the unperturbed Keplerian ellipse, by means of (Brumberg 1991;Poisson & Will 2014) It accounts for the instantaneous change of the line of the apsides; indeed, the orbital period P b is just the time required by the test particle to return at the (moving) pericentre position along its path. In Equation (80), are the radial and transverse components, respectively, of the perturbing acceleration A inducing the slow variation of the otherwise fixed Keplerian ellipse. In the present case, it is The second part, labeled in the following with the superscript "mix II", takes into account the J 2 -driven instantaneous changes experienced by the osculating Keplerian elements during an orbital revolution. The mean variation of any of the spin components' rates dS i /dt, i = 1, 2, 3 over an orbital period occurring due to the aforementioned shifts can be worked out as where f 0 is the true anomaly at a referenced epoch t 0 , and φ 1 a, φ 2 e, φ 3 I, φ 4 Ω, φ 5 ω.
The instantaneous shifts of the Keplerian orbital elements
dα dt 2 = 9 256 From Equations (108)-(109) it can be noted that, since the pericentre of a polar orbit, in general, does undergo a secular precession due to J 2 (Capderou 2005), the shift of the spin's right ascension is, actually, a harmonic signal with half the period 10 P ω of the pericentre, while the spin's declination experiences a genuine secular trend superimposed to a harmonic pattern with P ω /2.
In the case of GP-B, we plot its spin's declination precession as a function of f 0 in Figure 4. It can be noted that the predicted rate is larger than σ GP−Ḃ analytically computed 1pN J 2 -induced rate of changeδ (yellow curve, in mas yr −1 ) of the declination δ of the spin axis of a gyroscope orbiting the oblate Earth as a function of the initial value f 0 of the true anomaly of the gyro's orbit. We adopted the GP-B's orbital and spin configuration (Kahn 2007) (Everitt et al. 2011(Everitt et al. , 2015 in measuring the long-term rates of change of δ. Cfr. with Figure 3. 250 • , 325 • f 0 ≤ 360 • , with peaks of more than 30 mas yr −1 . A comparison with Figure 3 shows agreement between our analytical and numerical results up to a few mas yr −1 .
We do not display the total GP-B's 1pN right ascension rate due to J 2 since it turned out to be smaller than 2 mas yr −1 , while the reported experimental accuracy in measuringα is as large as σ GP−Ḃ α = 7.2 mas yr −1 (Everitt et al. 2011(Everitt et al. , 2015.

The gravitomagnetic spin precession
The long-term gravitomagnetic spin precession induced by the proper angular momentum J of the primary can be analytically worked out by including (Rindler 2001;Ruggiero & Tartaglia 2002) where is the 3-dimensional Levi-Civita symbol (Tyldesley 1975), in the spacetime metric tensor of Equations (5)- (7), and averaging the resulting J-dependent part of the expansion of Equation (18) to the order of O c −2 over a Keplerian ellipse. Smaller terms of the order of O J 2 J c −2 , arising from using a J 2 -driven precessing ellipse for the orbital average, will be neglected.

Summary and conclusions
The quadrupole mass moment J 2 of a body affects, among other things, also the general relativistic precession of the spin of an orbiting gyroscope. We worked out it, to the 1pN level, both numerically and analytically by taking into account also the effect that the J 2 -driven change of the gyro's orbit has on the the long-term spin rate itself. Indeed, limiting to averaging out the instantaneous J 2 -dependent part of the spin precession onto a Keplerian orbit is not sufficient to correctly reproduce the total spin rate of change to the order of O J 2 c −2 . Also the instantaneous Newtonian orbital shifts due to J 2 have to be taken into account when the average of the 1pN de Sitter-like instantaneous part of the spin precession is performed. The latter contribution introduces a dependence of the total averaged spin rate of the order of O J 2 c −2 on the initial orbital phase f 0 . Such a feature was confirmed, among other things, also by the simultaneous numerical integrations of the equations for the parallel transport of the spin and of the geodesic equations of the gyro's motion that we performed by varying f 0 .
We applied our results to the past GP-B mission in the field of Earth by finding a net precession of the declination of the spin axis which may be as large as ≃ 30 − 40 mas yr −1 . Since the reported error in measuring the GP-B's declination rate amounts to 18.3 mas yr −1 , our result may prompt a reanalysis of the data in order to see if the effect we predicted could be detected.
For the sake of completeness, we analytically worked out, to the 1pN level, also the general expression of the gravitomagnetic spin precession induced by the proper angular momentum J of the central body.
Both our numerical and analytical methods hold for an arbitrary orientation of the body's symmetry axis and for a general orbital configuration of the gyro. As such, they can be extended also to other astronomical and astrophysical scenarios of interest like, e.g., other planets of our solar system, exoplanets close to their parent stars, stars orbiting galactic supermassive black holes, tight binaries hosting compact stellar corpses. It is hardly necessary to mention that, years ago, spacecraft-based missions were proposed to measure the angular momenta of Jupiter and the Sun by means of the gravitomagnetic Pugh-Schiff spin precessions.