A post-Newtonian gravitomagnetic effect on the orbital motion of a test particle around its primary induced by the spin of a distant third body

We study a general relativistic gravitomagnetic 3-body effect induced by the spin angular momentum ${\boldsymbol S}_\textrm{X}$ of a rotating mass $M_\textrm{X}$ orbited at distance $r_\textrm{X}$ by a local gravitationally bound restricted two-body system $\mathcal{S}$ of size $r\ll r_\textrm{X}$ consisting of a test particle revolving around a massive body $M$. At the lowest post-Newtonian order, we analytically work out the doubly averaged rates of change of the Keplerian orbital elements of the test particle by finding non-vanishing long-term effects for the inclination $I$, the node $\Omega$ and the pericenter $\omega$. Such theoretical results are confirmed by a numerical integration of the equations of motion for a fictitious 3-body system. We numerically calculate the magnitudes of the post-Newtonian gravitomagnetic 3-body precessions for some astronomical scenarios in our solar system. For putative man-made orbiters of the natural moons Enceladus and Europa in the external fields of Saturn and Jupiter, the relativistic precessions due to the angular momenta of the gaseous giant planets can be as large as $\simeq 10-50~\textrm{milliarcseconds~per~year}~\left(\textrm{mas~yr}^{-1}\right)$. A preliminary numerical simulation shows that, for certain orbital configurations of a hypothetical Europa orbiter, its range-rate signal $\Delta\dot\rho$ can become larger than the current Doppler accuracy of the existing spacecraft Juno at Jupiter, i.e. $\sigma_{\dot\rho}=0.015~\textrm{mm~s}^{-1}$, after 1 d. The effects induced by the Sun's angular momentum on artificial probes of Mercury and the Earth are at the level of $\simeq 1-0.1~\textrm{microarcseconds~per~year}~\left(\mu\textrm{as~yr}^{-1}\right)$.


Abstract
We study a general relativistic gravitomagnetic 3-body effect induced by the spin angular momentum S X of a rotating mass M X orbited at distance r X by a local gravitationally bound restricted two-body system S of size r ≪ r X consisting of a test particle revolving around a massive body M. At the lowest post-Newtonian order, we analytically work out the doubly averaged rates of change of the Keplerian orbital elements of the test particle by finding non-vanishing long-term effects for the inclination I, the node Ω and the pericenter ω. Such theoretical results are confirmed by a numerical integration of the equations of motion for a fictitious 3-body system. We numerically calculate the magnitudes of the post-Newtonian gravitomagnetic 3-body precessions for some astronomical scenarios in our solar system. For putative man-made orbiters of the natural moons Enceladus and Europa in the external fields of Saturn and Jupiter, the relativistic precessions due to the angular momenta of the gaseous giant planets can be as large as ≃ 10 − 50 milliarcseconds per year mas yr −1 . A preliminary numerical simulation shows that, for certain orbital configurations of a hypothetical Europa orbiter, its range-rate signal ∆ρ can become larger than the current Doppler accuracy of the existing spacecraft Juno at Jupiter, i.e. σρ = 0.015 mm s −1 , after 1 d. The effects induced by the Sun's angular momentum on artificial probes of Mercury and the Earth are at the level of ≃ 1 − 0.1 microarcseconds per year µas yr −1 .
keywords General relativity and gravitation; Experimental studies of gravity; Experimental tests of gravitational theories; Satellite orbits

Introduction
Let us consider a local gravitationally bound restricted two-body system S composed by a test particle completing a full orbital revolution around a planet of mass M at distance r in a time interval P b , and a distant 3rd body X with mass M X ≫ M and proper spin S X around which S revolves at distance r X ≫ r with orbital period P X b . In general, M may be endowed with its own Newtonian and post-Newtonian mass and spin multipole moments (Soffel & Frutos 2016;Frutos-Alfaro & Soffel 2018) affecting the satellite's motion with known (Capderou 2005;Poisson & Will 2014) and less known (Angélil et al. 2014;Schärer et al. 2017;Schanner & Soffel 2018) Newtonian and post-Newtonian orbital effects like the classical oblateness-driven orbital precessions, the gravitoelectric Einstein pericentre shift, the gravitomagnetic Lense-Thirring effect, etc. Let us consider a kinematically rotating and dynamically non-rotating coordinate system K (Brumberg & Kopeikin 1989;Damour, Soffel & Xu 1994;Kopeikin, Efroimsky & Kaplan 2011) attached to M in geodesic motion through the external spacetime deformed by the mass-energy currents of X, assumed stationary in a kinematically and dynamically non-rotating coordinate system K X whose axes point towards the distant quasars (Brumberg & Kopeikin 1989;Damour, Soffel & Xu 1994;Kopeikin, Efroimsky & Kaplan 2011). The planetocentric motion of the test particle referred to K is further affected by two peculiar post-Newtonian 3-body effects: the time-honored De Sitter precession due to solely the mass M X (de Sitter 1916;Schouten 1918;Fokker 1920), and a gravitomagnetic shift due to S X which, to our knowledge, has never been explicitly and clearly calculated in the literature, if ever it had been. Our purpose is to analytically work out the latter effect at the lowest post-Newtonian order without any a-priori simplifying assumptions concerning both the orbital configurations of the planetocentric satellite's motion and the trajectory of the planet-satellite system S in the external field of X, and for an arbitrary orientation of S X in space. For previous, approximate calculation restricted to the orbital angular momentum of the Moon orbiting the Earth in the field of the rotating Sun, see Gill et al. (1992, Sec. 3.3.3).
The plan of the paper is as follows. In Section 2, we analytically work out the long-term rates of change of the Keplerian orbital elements of the test particle. Section 3 is devoted to the application of the obtained results to some astronomical scenarios in our solar system. We summarize our results and offer our conclusions in Section 4. For the benefit of the reader, Appendix A contains a list of the definitions of the symbols used in the paper, while their numerical values and tables are collected in Appendix B.

The doubly averaged satellite's orbital precessions
In the weak-field and slow-motion approximation, the gravitomagnetic 3-body potential induced by the angular momentum S X of the external spinning object X on the planetary satellite is (1) In Equation (1), G, c are the Newtonian constant of gravitation and the speed of light in vacuum, respectively, while L is the angular momentum of the test particle's orbital motion around M. Eq. (1) comes from Eq. (2.19) of Barker & O'Connell (1979, p. 155) for the interaction potential energy V S 1,S 2 of two spins S (1) , S (2) of masses m 1 , m 2 separated by a distance r and moving with relative speed v in the limit m 2 ≡ M X ≫ m 1 ≡ M, and by assuming that the spin S (1) is the orbital angular momentum of the planetocentric satellite's motion while S (2) is the spin angular momentum S X of the distant 3rd body X. Thus, r in Eq. (2.19) of Barker & O'Connell (1979, p. 155) has to be identified with r X , and r × P is nothing but the orbital angular momentum r X ×Mv X of the motion of S around M X . It is interesting to note that, with the same identifications, V S 1 and V S 2 of Eqs. (2.17)-(2.18) in Barker & O'Connell (1979, p. 155) yield the gravitoelectric De Sitter orbital precession for the planetocentric motion of the satellite and the gravitomagnetic Lense-Thirring effect for the X-centric orbit of M, respectively.
The perturbing potential U pert to be inserted into the Lagrange planetary equations for the rates of change of the osculating Keplerian orbital elements of the test particle (Bertotti, Farinella & Vokrouhlický 2003;Kopeikin, Efroimsky & Kaplan 2011), obtained by doubly averaging Equation (1) with respect to P b , P X b for arbitrary orbital configurations of both the external body X and the test particle and for a generic orientation of S X in space, is with U = cos I 2Ŝ z − 3 sin I X Ŝ z sin I X + cos I X Ŝ y cos Ω X −Ŝ x sin Ω X + In Equations (2) and (3), a, e, I, a X , e X , I X are the semimajor axes, the eccentricities and the inclinations of the orbits of the test particle and of S, respectively, while n b is the Keplerian orbital motion of the satellite's planetary motion about M. Equations (2) and (3) were obtained in two steps. First, U GM of Equation (1) was evaluated onto the unperturbed ellipse of the planetocentric satellite motion through the standard Keplerian formulas of the restricted two-body problem (see, e.g., Equations (3.40a) to (3-41c) of Poisson & Will (2014)). Then, it was averaged over one orbital period P b to the first order in the disturbing potential by using just the Keplerian part of Equation (3.66) of Poisson & Will (2014) for d f /dt, where f is the true anomaly. Then, the resulting averaged potential U GM was, in turn, calculated onto the unperturbed X-centric Keplerian trajectory of S and averaged over P X b to the first order in the perturbation under consideration, thus finally obtaining the double average of Equations (2) and (3).
Inserting Equations (2) and (3) into the right-hand-sides of the Lagrange planetary equations allows to calculate the doubly averaged rates of change of the Keplerian orbital elements. They turn out to beȧ with We remark that Equations (4) and (11) are exact in both e and e X in the sense that the low-eccentricity approximation was not adopted in the calculation.
A more computationally cumbersome approach to obtain the same long-term rates of change of Equations (4) and (11) consists, first of all, in deriving a perturbing acceleration from Equation (1). By writing the Lagrangian per unit mass of a gravitationally bound restricted two-body system affected by a generic perturbing potential as the conjugate momentum per unit mass is, by definition, The Hamiltonian per unit mass is From the Hamilton equations of motion, it iṡ Since L pert = −H pert , by comparing Equation (14) and Equation (16), it turns out that the perturbing acceleration is just In our specific case, since H pert = U GM , we have Then, Equation (18) must be decomposed into its radial (R), transverse (T ) and out-of-plane (N) components, which are They have to be inserted into the right-hand-sides of the standard Gauss equations for the variation of the orbital elements (Bertotti, Farinella & Vokrouhlický 2003;Kopeikin, Efroimsky & Kaplan 2011;Poisson & Will 2014) which, finally, are doubly averaged with respect to P b , P X b in the same way as previously described for the disturbing potential of Equation (1).
We successfully checked our analytical results of Equations (4) to (11) as follows. We considered a fictitious system S orbiting a Jupiter-like body X along the same orbit of Callisto, whose mass was assumed for the particle's primary M, and numerically integrated its equations of motion over a time span much longer than P b , P X b with and without the post-Newtonian gravitomagnetic acceleration of Equation (18); both the integrations, which assumed a purely Keplerian motion of S about the fictitious body X, shared the same initial conditions for the test particle and its primary. For X, the same physical properties of Jupiter were assumed, including the size and the orientation of its angular momentum S. As a result, numerically produced times series of the orbital elements of the imaginary probe were produced by subtracting the purely Newtonian ones from those obtained by including also Equation (18) in the equations of motion; they are displayed in Figure 1. It turned out that the resulting numerically calculated post-Newtonian gravitomagnetic 3-body orbital shifts agree with those computed by means of the analytical formulas of Equations (4) to (11).

Some potentially interesting astronomical scenarios
For the sake of simplicity, we will consider a circular orbit (e = 0) for the test particle motion around M.
In the case of a hypothetical orbiter of the Kronian natural satellite Enceladus in the external field of Saturn, Equations (6) and (7) and Equations (9) and (10), referred to the mean Earth's equator at the reference epoch J2000.0 as reference {x, y} plane, yielḋ Ω = −49.9 mas yr −1 + cot I A eq cos Ω + ϕ eq , with A eq = −5.7 mas yr −1 , ϕ eq = 49.4 deg .
Instead, if the mean ecliptic at the reference epoch J2000.0 is adopted as reference {x, y} plane, we haveİ with A ecl = 23.9 mas yr −1 , By looking at a putative orbiter of the Jovian natural satellite Europa in the external field of Jupiter, we haveİ Ω = −9.9 mas yr −1 + cot I A eq cos Ω + ϕ eq , with A eq = 4.8 mas yr −1 , ϕ eq = 2.9 deg, with ϕ ecl = 31.0 deg .
We considered just Enceladus and Europa because they are of great planetological interest in view of the possible habitability of their oceans beneath their icy crusts (Lunine 2017).
As a consequence, they are the natural targets of several concept studies and proposals for dedicated missions to them, including also orbiters (Razzaghi et al. 2008;Spencer & Niebur 2010;MacKenzie et al. 2016;Verma & Margot 2018;Sherwood et al. 2018). Since, at present, sending a spacecraft to Europa seems more likely than to Enceladus, as it can be learnt at https://europa.nasa.gov/about-clipper/overview/ and http://sci.esa.int/juice/ on the Internet, we investigated in a little more detail this potentially appealing Jovian scenario, even if it is not said that the actually approved missions will finally involve the use of an orbiter. In such kind of endeavours, the observable quantity is typically the Earth-probe range-rateρ, whose accuracy for, e.g., the ongoing mission Juno (Bolton et al. 2017) around Jupiter is σρ ≃ 0.015 mm s −1 (Iess et al. 2018). Figure 2 shows the numerically simulated Earth-spacecraft range-rate signature due to the post-Newtonian gravitomagnetic 3-body acceleration of Equation (18) for a generic orbital configuration of the hypothesized orbiter. To produce it, we numerically integrated the equations of motion in Cartesian rectangular coordinates of the Earth, Jupiter, its Galilean moons and a fictitious test particle orbiting Europa over 1 d. In both runs, sharing the same initial conditions retrieved from the database JPL HORIZONS (https://ssd.jpl.nasa.gov/?horizons) at the arbitrary epoch of midnight of 1st January 2030, we modeled the mutual attractions among all the bodies involved to the Newtonian level, with the exception of Equation (18) which was added to the other classical gravitational pulls felt by the probe in one of the runs. Then, we numerically calculated two range-rate time series, and subtracted the purely Newtonian one from that including also the post-Newtonian gravitomagnetic acceleration. It can be noted that, for the orbital configuration chosen, the range-rate relativistic signature ∆ρ reaches the 0.05 mm s −1 level after just 1 d. Thus, the scenario considered seems worth of further, dedicated analyses investigating the actual measurability of Equation (18) in a realistic error budget analysis and mission proposal. It should take into account several concurring perturbations of gravitational and non-gravitational nature, and also several technological and engineering issues.
In the case of an artificial satellite orbiting a planet in the field of the Sun, the effects are much smaller. For an Earth's spacecraft, we havė with ϕ eq = 9.13 deg, with For a probe orbiting Mercury one getṡ with with

Summary and overview
In the weak-field and slow-motion approximation of general relativity, we analytically worked out the post-Newtonian gravitomagnetic long-term rates of change of the relevant Keplerian orbital elements of a test particle orbiting a primary M at distance r from it which, in turn, moves in the external spacetime deformed by the mass-energy currents of the spin angular momentum S X of a distant (r X ≫ r) 3rd body X with mass M X ≫ M. We did not assume any preferred orientation for the spin axisŜ X of the external body; moreover, we did not make simplifying assumptions pertaining the orbital configurations of both the M's satellite and of M itself in its motion around M X . Thus, our calculation have a general validity, being applicable to arbitrary astronomical systems of potential interest. It turns out that, by doubly averaging the perturbing potential employed in the calculation with respect to the orbital periods P b , P X b of both M and M X , the semimajor axis a and the eccentricity e do not experience long-term variations, contrary to the inclination I of the orbital plane, the longitude of the ascending node Ω and the argument of pericenter ω. While the gravitomagnetic ratesİ andω are harmonic signatures characterized by the frequency of the possible variation of the node Ω, induced by other dominant perturbations like, e.g., the Newtonian quadrupole mass moment of the satellite's primary M, the gravitomagnetic node rateΩ exhibits also a secular trend in addition to a harmonic component with the frequency of the node itself. A numerical integration of the equations of motion of a fictitious 3-body system made of a distant spinning body with the same physical properties of Jupiter, a primary with the same orbital and physical characteristics of Callisto and a test particle orbiting it confirms our analytical results.
The Sun's angular momentum exerts very small effects on spacecraft orbiting Mercury (≃ 1 µas yr −1 ) and the Earth (∼ 0.1 µas yr −1 ). Instead, the angular momenta of the gaseous giant planets like Jupiter and Saturn may induce much larger perturbations of the orbital motions of hypothetical anthropogenic orbiters of some of their major natural moons like, e.g., Europa ( 10 mas yr −1 ) and Enceladus ( 50 mas yr −1 ). Such natural satellites have preeminent interest in planetology making them ideal targets for future, dedicated spacecraft-based missions which may be opportunistically exploited to attempt to measure such relativistic effects as well. In the case of Europa, for whose exploration there are already approved missions by NASA and ESA, a preliminary numerical simulation of the signature induced by the post-Newtonian gravitomagnetic 3-body effect of interest on the range-rate of a putative orbiter shows that, for certain orbital configurations, its magnitude can become larger than the present-day accuracy σρ = 0.015 mm s −1 of the current Juno mission around Jupiter after 1 d.
S eq y = cos δ X sin α X : component of the 3rd-body's spin axis w.r.t. the reference y axis of an equatorial coordinate system S eq z = sin δ X : component of the 3rd-body's spin axis w.r.t. the reference z axis of an equatorial coordinate system r X : position vector towards the 3rd-body X r X : distance of S to the 3rd-body X r X r X /r X : versor of the position vector towards the 3rd-body X a X : semimajor axis of the orbit about the 3rd-body X n X b µ X /a 3 X : mean motion of the orbit about the 3rd-body X P X b 2π/n X b : orbital period of the orbit about the 3rd-body X e X : eccentricity of the orbit about the 3rd-body X I X : inclination of the orbital plane of orbit about the 3rd-body X to the reference {x, y} plane of some coordinate system Ω X : longitude of the ascending node of the orbit about the 3rd-body X referred to the reference {x, y} plane of some coordinate system M : mass of the primary (planet or planetary natural satellite) orbited by the test particle and moving in the external field of the 3rd-body X µ GM : gravitational parameter of the primary orbited by the test particle and moving in the external field of the 3rd-body X R : radius of the primary (planet or planetary natural satellite) orbited by the test particle and moving in the external field of the 3rd-body X S : angular momentum of the primary J ℓ , ℓ = 2, 3, 4 : zonal multipole moments of the classical gravitational potential of the primary r : position vector of the test particle with respect to its primary r : magnitude of the position vector of the test particle v : velocity vector of the test particle L r × v : orbital angular momentum per unit mass of the test particle a : semimajor axis of the test particle's orbit n b µ/a 3 : Keplerian mean motion of the test particle's orbit P b 2π/n b : orbital period of the test particle's orbit e : eccentricity of the test particle's orbit f : true anomaly of the test particle's orbit I : inclination of the orbital plane of the test particle's orbit to the reference {x, y} plane of some coordinate system Ω : longitude of the ascending node of the test particle's orbit referred to the reference {x, y} plane of some coordinate system   1.-Numerically computed time series of the post-Newtonian gravitomagnetic 3rd-body shifts experienced by the inclination I, node Ω and pericentre ω of a fictitious test particle moving around a Callisto-like primary M which, in turn, orbits a Jupiter-type 3rd-body X. They were obtained by numerically integrating the equations of motion of the orbiter in Cartesian rectangular coordinates referred to the Earth's mean equator at the epoch J2000.0 with and without Equation (18). Both runs shared the same set of arbitrary initial conditions for the probe P b = 10.07 d, e 0 = 0.3, I 0 = 80 deg, Ω 0 = 230 deg, ω 0 = 40 deg, f 0 = 50 deg and the primary; as far as the motion of M with respect to X is concerned, the initial state vector of the Callisto-Jupiter relative motion was adopted from the database JPL HORIZONS (https://ssd.jpl.nasa.gov/?horizons). For each Keplerian orbital element, its time series calculated from the purely Newtonian run was subtracted from that obtained from the post-Newtonian integration in order to obtain the signatures displayed here. The resulting shifts, in mas yr −1 , agree with the analytically computed ones in Equations (6) to (8).  (18). We numerically integrated the solar system barycentric equations of motion in Cartesian rectangular coordinates of the Earth, Jupiter, its Galilean moons and a fictitious test particle orbiting Europa over 1 d. In both runs, sharing the same initial conditions for all the existing natural bodies retrieved from the database JPL HORIZONS (https://ssd.jpl.nasa.gov/?horizons) at the arbitrary epoch of midnight of 1st January 2030, we modeled the mutual attractions among all the planets and the satellites involved to the Newtonian level, with the exception of Equation (18) which was added to the other classical gravitational pulls felt by the orbiter in one of the runs. Then, we numerically calculated two range-rate time series, and subtracted the purely Newtonian one from that including also Equation (18). The orbital configuration adopted for the spacecraft, referred to Europa, was a 0 = 3.55 R, e 0 = 0.69, I 0 = 100 deg, Ω 0 = 90 deg, ω 0 = 40 deg, f 0 = 50 deg, where R is the radius of Jovian moon.