When the Anomalistic, Draconitic and Sidereal Orbital Periods Do Not Coincide: The Impact of Post-Keplerian Perturbing Accelerations

Abstract

In a purely Keplerian picture, the anomalistic, draconitic and sidereal orbital periods of a test particle orbiting a massive body coincide with each other. Such degeneracy is removed when post-Keplerian perturbing acceleration enters the equations of motion, yielding generally different corrections to the Keplerian period for the three aforementioned characteristic orbital timescales. They are analytically worked out in the case of the accelerations induced by the general relativistic post-Newtonian gravitoelectromagnetic fields and, to the Newtonian level, by the oblateness of the central body. The resulting expressions hold for completely general orbital configurations and spatial orientations of the spin axis of the primary. Astronomical systems characterized by extremely accurate measurements of orbital periods like transiting exoplanets and binary pulsars may offer potentially viable scenarios for measuring such post-Keplerian features of motion, at least in principle. As an example, the sidereal period of the brown dwarf WD1032 + 011 b is currently known with an uncertainty as small as ≃${10}^{-5}\mathrm{s}$, while its predicted post-Newtonian gravitoelectric correction amounts to $0.07\mathrm{s}$; however, the accuracy with which the Keplerian period can be calculated is just 572 s. For double pulsar PSR J0737–3039, the largest relativistic correction to the anomalistic period amounts to a few tenths of a second, given a measurement error of such a characteristic orbital timescale as small as $\simeq {10}^{-6}\mathrm{s}$. On the other hand, the Keplerian term can be currently calculated just to a $\simeq 9$ s accuracy. In principle, measuring at least two of the three characteristic orbital periods for the same system independently would cancel out their common Keplerian component, provided that their difference is taken into account.

1. Introducción

From a theoretical point of view, various time intervals (T) characterizing different cyclic patterns of the orbital motion of a two-body gravitationally bound system can be defined when, in addition to the dominant Newtonian inverse–square acceleration, a much smaller, post-Keplerian (Here, post-Keplerian means dynamical features arising from any acceleration—Newtonian or not—different from the simple Newtonian inverse–square one. Then, in this specific sense, the classical acceleration due to, say, the primary’s oblateness, $J_2$ has to refer to pK.) one ($\mathit{A}$) acts on a satellite. Such characteristic orbital timescales are the amounts of time elapsed between two successive passages of the latter in some directions which, in a purely Keplerian scenario, are all fixed; in this case, all such periods coincide with the Keplerian one ($T_{\mathrm{K}}$). Instead, a pK perturbation breaks such a degeneracy, and the aforementioned temporal intervals generally differ from each other.

The aim of the present work is to analytically work out the corrections ($\Delta T$ to $T_{\mathrm{K}}$) induced by the gravitoelectromagnetic accelerations arising within the General Theory of Relativity (GTR) to the first post-Newtonian (1 pN) order [1,2]. Furthermore, the impact of the oblateness of the central body [3] is worked out to the Newtonian level. In all three aforementioned cases, the anomalistic, draconitic and sidereal periods are considered. The resulting expressions turn out to be valid for completely general orbital shapes and inclinations and for arbitrary orientations of the primary’s spin axis in space. For analogous calculation restricted to some particular orbital geometries, see Iorio [4]. In addition to the traditional quantities usually adopted, like, the time-honored pericenter precession, the orbital periods, if directly measured, may offer, in principle, further ways to test the GTR and other models of gravity.

In recent years, exoplanets [5,6,7,8,9] have been attracting a growing interest as possible tools to test the GTR and modified models of gravity [10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34]. For comprehensive overviews of tests of the GTR, see, e.g., Will [35,36], Will and Yunes [37] and references therein. The results presented here may have an impact on, e.g., just exoplanetary studies, since the accuracy in measuring the orbital periods of some transiting planets is currently quite remarkable. Even better is the accuracy with which the anomalistic periods of binary pulsars are usually measured.

This paper is organized as follows. In Section 2, the treated pK accelerations are reviewed, including 1 pN gravitoelectric (Section 2.1), gravitomagnetic (Section 2.2) and Newtonian quadrupolar (Section 2.3) accelerations. Section 3 is devoted to the anomalistic period; the general calculational scheme is outlined in Section 3.1, while the 1 pN and quadrupolar corrections are worked out in Section 3.2, Section 3.3 and Section 3.4. The draconitic period is dealt with in Section 4; Section 4.1 shows how to calculate it, while the 1 pN and quadrupolar corrections are the subjects of Section 4.2, Section 4.3 and Section 4.4. The sidereal period is investigated in Section 5; the calculational approach is explained in Section 5.1, while the 1 pN and quadrupolar corrections are calculated in Section 5.2, Section 5.3 and Section 5.4. In Section 6, a numerical evaluation of the 1 pN gravitoelectric effect for a transiting exoplanet and the double pulsar whose orbital periods are accurately measured is offered. Section 7 summarizes the findings and offers conclusions.

2. The pK Accelerations

Here, a brief summary of some key concepts of celestial mechanics needed to follow the rest of the paper profitably is offered [38,39,40,41,42,43,44,45].

If $\mathit{A}$ is an arbitrary pK perturbing acceleration generally depending on the position and velocity vectors ($\mathit{r}$ y $\mathit{v}$) of the orbiter, calculating its effects on the orbital path of the latter requires knowledge of its radial, transverse and normal components ($A_r$, $A_{\tau }$ y $A_h$). They are the projections ($\mathit{A}·\widehat{\mathit{r}}$, $\mathit{A}·\widehat{\mathit{\tau }}$ y $\mathit{A}·\widehat{\mathit{h}}$) of $\mathit{A}$ onto the co-moving radial, transverse and normal unit vectors

$$\begin{array}{cc}\hfill \widehat{\mathit{r}}& =\left\{cos\mathsf{\Omega }cosu-cosIsin\mathsf{\Omega }sinu,sin\mathsf{\Omega }cosu+cosIcos\mathsf{\Omega }sinu,sinIsinu\right\},\hfill \\ \hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \widehat{\mathit{\tau }}& =\left\{-cos\mathsf{\Omega }sinu-cosIsin\mathsf{\Omega }cosu,-sin\mathsf{\Omega }sinu+cosIcos\mathsf{\Omega }cosu,sinIcosu\right\},\hfill \\ \hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \widehat{\mathit{h}}& =\left\{sinIsin\mathsf{\Omega },-sinIcos\mathsf{\Omega },cosI\right\}.\hfill \end{array}$$

(3)

In Equations (1)–(3), I is the inclination of the orbital plane relative to the reference plane of the adopted coordinate system adopted, $\mathsf{\Omega }$ is the longitude of the ascending node (an angle counted in the reference plane from the reference direction (x) to the line of nodes, which is the intersection of the orbital plane with the fundamental one) and u is the argument of latitude (a time-dependent angle reckoned in the orbital plane from the line of nodes to the instantaneous position of the test particle moving along the ellipse) defined as

$$u:=\omega +f.$$

(4)

In Equation (4), $\omega$ is the argument of the pericenter (an angle reckoned in the orbital plane from the line of nodes to the position of the pericenter on the line of apsides), and f is the true anomaly (a time-dependent angle counted in the orbital plane from the position of the pericenter on the line of apsides to the instantaneous position of the test particle moving along the ellipse; thus, it is often used as fast variable of integration when the average over one orbital period of some relevant quantity instantaneously varying during the satellite’s orbital revolution is calculated). The normal unit vector given by Equation (3) is aligned with the orbital angular momentum; Equations (1)–(3) are connected at each instant of time by the following relation:

$$\widehat{\mathit{h}}\times \widehat{\mathit{r}}=\widehat{\mathit{\tau }}.$$

(5)

Then, $A_r$, $A_{\tau }$ y $A_h$ are evaluated within a Keplerian ellipse, which is assumed to be an unperturbed reference trajectory, as expressed by

$$r=\frac{p}{1+ecosf}.$$

(6)

In Equation (6), e is the eccentricity (fixes the shape of the ellipse; $0\le e<1$, where $e=0$ corresponds to a circular orbit), and

$$p:=a\left(1-e^2\right)$$

(7)

is the orbit’s semilatus rectum, and a is the semimajor axis (determines the size of the ellipse). The position and velocity vectors generally entering the analytical expression of $\mathit{A}$ can be conveniently expressed as

$$\begin{array}{cc}\hfill \mathit{r}& =r\left(\widehat{\mathit{l}}cosu+\widehat{\mathit{m}}sinu\right),\hfill \\ \hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \mathit{v}& =\sqrt{\frac{\mathsf{\mu }}{p}}\left[-\widehat{\mathit{l}}\left(esin\omega +sinu\right)+\widehat{\mathit{m}}\left(ecos\omega +cosu\right)\right].\hfill \end{array}$$

(9)

In Equations (8) and (9), r is given by Equation (6), while

$$\begin{array}{cc}\hfill \widehat{\mathit{l}}& :=\left\{cos\mathsf{\Omega },sin\mathsf{\Omega },0\right\},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \widehat{\mathit{m}}& :=\left\{-cosIsin\mathsf{\Omega },cosIcos\mathsf{\Omega },sinI\right\}\hfill \end{array}$$

(11)

are two unit vectors lying in the orbital plane. $\widehat{\mathit{l}}$ is directed along the line of nodes, while $\widehat{\mathit{m}}$ is perpendicular to $\widehat{\mathit{l}}$ in such a way that

$$\widehat{\mathit{l}}\times \widehat{\mathit{m}}=\widehat{\mathit{h}}.$$

(12)

Finally,

$$\mathsf{\mu }:=GM$$

(13)

entering Equation (9) is the standard gravitational parameter of the primary, M is its mass and G is the Newtonian constant of gravitation.

2.1. The 1 pN Gravitoelectric Acceleration

In the case of a binary system made of two non-rotating bodies (A and B) of finite masses ($M_{\mathrm{A}}$ y $M_{\mathrm{B}}$), the 1 pN gravitoelectric acceleration is (see, e.g., [46] (Equation (2.5), p. 111]); (Equation (A2.6), p. 166) [44]; (Equation (4.4.28), p. 154) [40]; (Equation (10.3.7), p. 381) [45]; and (Equation (10.1), p. 482) [42])

$$\begin{array}{cc}\hfill {\mathit{A}}^{1\mathrm{pN}}& =\frac{{\mathsf{\mu }}_{\mathrm{b}}}{c^2{r}^2}\left\{\left[\left(4+2\nu \right)\frac{{\mathsf{\mu }}_{\mathrm{b}}}{r}+\frac{3}{2}\nu v_r^2-\left(1-3\nu \right)v^2\right]\widehat{\mathit{r}}+\left(4-2\nu \right)v_r\mathit{v}\right\}.\hfill \end{array}$$

(14)

In Equation (14), c is the speed of light in vacuum, and

$${\mathsf{\mu }}_{\mathrm{b}}:=G{M}_{\mathrm{b}}$$

(15)

is the standard gravitational parameter of the binary whose total mass is

$$M_{\mathrm{b}}:=M_{\mathrm{A}}+M_{\mathrm{B}}.$$

(16)

Moreover,

$$\nu :=\frac{M_{\mathrm{A}}{M}_{\mathrm{B}}}{M_{\mathrm{b}}^2}$$

(17)

is the binary’s symmetric mass ratio, ranging from 0 if one of the two bodies can be considered a test particle to $1/4=0.25$ if both bodies have the same mass; $\mathit{r}$ is the position vector of one body with respect to the other one; $\widehat{\mathit{r}}$ is its unit vector; r is its magnitude, corresponding to the relative distance between the two bodies; $\mathit{v}$ is the relative velocity between them, whose magnitude is v; while

$$v_r:=\mathit{v}·\widehat{\mathit{r}}$$

(18)

is the projection of the relative velocity onto the position unit vector. By using Equations (1)–(4) and (6)–(11), the radial, transverse and normal components ($A_r^{1\mathrm{pN}}$, $A_{\tau }^{1\mathrm{pN}}$ y $A_h^{1\mathrm{pN}}$, respectively) of Equation (14) turn out to be

$$\begin{array}{cc}\hfill A_r^{1\mathrm{pN}}& =\frac{{\mathsf{\mu }}_{\mathrm{b}}^2{\left(1+ecosf\right)}^2}{4c^2{a}^3{\left(1-e^2\right)}^3}\left[e^2\left(4-13\nu \right)-4\left(-3+\nu \right)+8e\left(1-2\nu \right)cosf+e^2\left(-8+\nu \right)cos2f\right],\hfill \\ \hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill A_{\tau }^{1\mathrm{pN}}& =\frac{2e{\mathsf{\mu }}_{\mathrm{b}}^2{\left(1+ecosf\right)}^3\left(2-\nu \right)sinf}{c^2{a}^3{\left(1-e^2\right)}^3},\hfill \\ \hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill A_h^{1\mathrm{pN}}& =0.\hfill \end{array}$$

(21)

Equations (19)–(21) agree with, e.g., Equations (A2.77a)–(A2.77c), calculated with the GTR, by (Soffel [44], p. 178).

2.2. The 1 pN Gravitomagnetic Lense–Thirring Acceleration

The 1pN gravitomagnetic Lense–Thirring (LT) acceleration due to the angular momentum $\mathit{J}$ of a massive primary is, for an arbitrary orientation of the latter, [42,44,45,47,48,49]

$${\mathit{A}}^{\mathrm{LT}}=\frac{2GJ}{c^2{r}^3}\left[3\left(\widehat{\mathit{J}}·\widehat{\mathit{r}}\right)\widehat{\mathit{r}}\times \mathit{v}+\mathit{v}\times \widehat{\mathit{J}}\right].$$

(22)

In Equation (22), $\widehat{\mathit{J}}$ is the primary’s spin axis, which can be parameterized as, e.g.,

$$\begin{array}{cc}\hfill {\widehat{J}}_x& =cos{\alpha }_{J}cos{\delta }_J,\hfill \\ \hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {\widehat{J}}_y& =sin{\alpha }_{J}cos{\delta }_J,\hfill \\ \hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\widehat{J}}_z& =sin{\delta }_J,\hfill \end{array}$$

(25)

where ${\alpha }_J$ y ${\delta }_J$ are its longitude and latitude angles, respectively in some coordinate system. For a generalization of Equation (22) to a two-body system with comparable masses and spins, see, e.g., Kidder [50] (Equation (2.2.c)) and Soffel [44]. Basically, if ${\mathit{J}}_{\mathrm{A}}$ y ${\mathit{J}}_{\mathrm{B}}$ are the spin angular momenta of extended bodies A and B, $\mathit{J}$ has to be replaced in Equation (22) and in all the following equations with

$${\mathit{S}}_{\mathrm{b}}:=\left(1+\frac{3}{4}\frac{M_{\mathrm{B}}}{M_{\mathrm{A}}}\right){\mathit{J}}_{\mathrm{A}}+\left(1+\frac{3}{4}\frac{M_{\mathrm{A}}}{M_{\mathrm{B}}}\right){\mathit{J}}_{\mathrm{B}}$$

(26)

It is useful to define the following quantities:

$$\begin{array}{c}\hfill \mathtt{Jl}:=\widehat{\mathit{J}}·\widehat{\mathit{l}};\\ \hfill \end{array}$$

(27)

$$\begin{array}{c}\hfill \mathtt{Jm}:=\widehat{\mathit{J}}·\widehat{\mathit{m}};\\ \hfill \end{array}$$

(28)

$$\begin{array}{c}\hfill \mathtt{Jh}:=\widehat{\mathit{J}}·\widehat{\mathit{h}}.\end{array}$$

(29)

By means of Equations (1)–(4), (6)–(9) and (27)–(29), the radial, transverse and normal components ($A_r^{\mathrm{LT}}$, $A_{\tau }^{\mathrm{LT}}$ y $A_h^{\mathrm{LT}}$, respectively) of Equation (22) can be written as

$$\begin{array}{cc}\hfill A_r^{\mathrm{LT}}& =\frac{2n_{\mathrm{K}}GJ{\left(1+ecosf\right)}^4\mathtt{Jh}}{c^2{a}^2{\left(1-e^2\right)}^{7/2}},\hfill \\ \hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill A_{\tau }^{\mathrm{LT}}& =-\frac{2e{n}_{\mathrm{K}}GJ{\left(1+ecosf\right)}^{3}sinf\mathtt{Jh}}{c^2{a}^2{\left(1-e^2\right)}^{7/2}},\hfill \\ \hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill A_h^{\mathrm{LT}}& =-\frac{2n_{\mathrm{K}}GJ{\left(1+ecosf\right)}^3}{c^2{a}^2{\left(1-e^2\right)}^{7/2}}\left\{\left[ecos\omega -\left(2+3ecosf\right)cosu\right]\mathtt{Jl}-\frac{1}{2}\left[esin\omega +4sinu+3esin\left(2f+\omega \right)\right]\mathtt{Jm}\right\},\hfill \end{array}$$

(32)

where

$$n_{\mathrm{K}}:=\sqrt{\frac{\mathsf{\mu }}{a^3}}=\frac{2\mathsf{\pi }}{T_{\mathrm{K}}}$$

(33)

is the Keplerian mean motion.

2.3. The Newtonian Quadrupolar Acceleration

The pK acceleration induced by the first even zonal harmonic coefficient ($J_2$) of the multipolar expansion of the exterior Newtonian gravitational potential of a massive primary endowed with axial symmetry is

$${\mathit{A}}^{J_2}=\frac{3\mathsf{\mu }{J}_2{R}_{\mathrm{e}}^2}{2r^4}\left\{\left[5{\left(\widehat{\mathit{J}}·\widehat{\mathit{r}}\right)}^2-1\right]\widehat{\mathit{r}}-2\left(\widehat{\mathit{J}}·\widehat{\mathit{r}}\right)\widehat{\mathit{J}}\right\},$$

(34)

where $R_{\mathrm{e}}$ is the body’s equatorial radius.

By defining

$$\begin{array}{cc}\hfill {\widehat{T}}_1& :=1;\hfill \\ \hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill {\widehat{T}}_2& :={\mathtt{Jl}}^2+{\mathtt{Jm}}^2;\hfill \\ \hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill {\widehat{T}}_3& :={\mathtt{Jl}}^2-{\mathtt{Jm}}^2;\hfill \\ \hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill {\widehat{T}}_4& :=\mathtt{Jh}\mathtt{Jl};\hfill \\ \hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill {\widehat{T}}_5& :=\mathtt{Jh}\mathtt{Jm};\hfill \\ \hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill {\widehat{T}}_6& :=\mathtt{Jl}\mathtt{Jm}\hfill \end{array}$$

(40)

and by means of Equations (1)–(4) and (6)–(8), the radial, transverse and normal components ($A_r^{J_2}$, $A_{\tau }^{J_2}$ y $A_h^{J_2}$, respectively) of Equation (34) can be cast into the form of

$$\begin{array}{cc}\hfill A_r^{J_2}& =\frac{3\mathsf{\mu }{J}_2{R}_{\mathrm{e}}^2{\left(1+ecosf\right)}^4}{2a^4{\left(1-e^2\right)}^4}\left[-{\widehat{T}}_1+3\left(\frac{{\widehat{T}}_2}{2}+\frac{{\widehat{T}}_{3}cos2u}{2}+{\widehat{T}}_{6}sin2u\right)\right];\hfill \\ \hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill A_{\tau }^{J_2}& =\frac{3\mathsf{\mu }{J}_2{R}_{\mathrm{e}}^2{\left(1+ecosf\right)}^4}{a^4{\left(1-e^2\right)}^4}\left(\frac{{\widehat{T}}_{3}sin2u}{2}-{\widehat{T}}_{6}cos2u\right);\hfill \\ \hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill A_h^{J_2}& =-\frac{3\mathsf{\mu }{J}_2{R}_{\mathrm{e}}^2{\left(1+ecosf\right)}^4}{a^4{\left(1-e^2\right)}^4}\left({\widehat{T}}_{4}cosu+{\widehat{T}}_{5}sinu\right).\hfill \end{array}$$

(43)

If bodies A and B are both extended and axisymmetric, Equation (34) becomes [51]

$${\mathit{A}}^{J_2}=\frac{3{\mathsf{\mu }}_{\mathrm{b}}}{2r^4}{\mathcal{F}}_{\mathrm{AB}},$$

(44)

where

$${\mathcal{F}}_{\mathrm{AB}}:=J_2^{\mathrm{A}}{R_{\mathrm{e}}^{\mathrm{A}}}^2\left\{\left[5{\left({\widehat{\mathit{J}}}^{\mathrm{A}}·\widehat{\mathit{r}}\right)}^2-1\right]\widehat{\mathit{r}}-2\left({\widehat{\mathit{J}}}^{\mathrm{A}}·\widehat{\mathit{r}}\right){\widehat{\mathit{J}}}^{\mathrm{A}}\right\}+J_2^{\mathrm{B}}{R_{\mathrm{e}}^{\mathrm{B}}}^2\left\{\left[5{\left({\widehat{\mathit{J}}}^{\mathrm{B}}·\widehat{\mathit{r}}\right)}^2-1\right]\widehat{\mathit{r}}-2\left({\widehat{\mathit{J}}}^{\mathrm{B}}·\widehat{\mathit{r}}\right){\widehat{\mathit{J}}}^{\mathrm{B}}\right\},$$

(45)

y r y $\widehat{\mathit{r}}$ entering Equations (44) and (45) refer to the relative orbit.

3. The Apsidal Period

3.1. General Calculational Scheme

The anomalistic period ($T_{\mathrm{ano}}$) is defined as the time interval between two successive instants when the real position of the test particle coincides with the pericenter position on the corresponding orbit. It can be calculated as [4,52,53]

$$T_{\mathrm{ano}}=T_{\mathrm{K}}+\Delta T_{\mathrm{ano}}={\int }_0^{2\mathsf{\pi }}\left(\frac{\mathrm{d}t}{\mathrm{d}f}\right)df,$$

(46)

where $dt/df$, when a pK acceleration ($\mathit{A}$) is present, is given by

$$\frac{\mathrm{d}t}{\mathrm{d}f}\simeq \frac{r^2}{\sqrt{\mathsf{\mu }p}}\left[1+\frac{r^2}{\sqrt{\mathsf{\mu }p}}\left(\frac{\mathrm{d}\omega }{\mathrm{d}t}+cosI\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}t}\right)\right],$$

(47)

since [38,40,42,52,54,55]

$$\frac{\mathrm{d}f}{\mathrm{d}t}=\frac{\sqrt{\mathsf{\mu }p}}{r^2}\left[1-\frac{r^2}{\sqrt{\mathsf{\mu }p}}\left(\frac{\mathrm{d}\omega }{\mathrm{d}t}+cosI\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}t}\right)\right].$$

(48)

The true anomaly (f) enters Equation (46) as a fast variable of integration just because the line of apsides is involved in the definition of the anomalistic period. In order to obtain the full correction ($\Delta T_{\mathrm{ano}}$) of the order of A to the Keplerian orbital period, the contribution of the second term of Equation (47) to Equation (46) is not enough. The partial derivatives of the Keplerian term of Equation (47), namely

$${\left.\frac{\mathrm{d}t}{\mathrm{d}f}\right|}_{\mathrm{K}}=\frac{r^2}{\sqrt{\mathsf{\mu }p}}$$

(49)

with respect to a y e multiplied by the finite variations ($\Delta a\left(f\right),\Delta e\left(f\right)$) of the same orbital elements have to be taken into consideration; in this way, one fully accounts for the fact that the Keplerian orbital elements vary instantaneously as the satellite travels along its trajectory. Thus, the following is obtained:

$$\begin{array}{cc}\hfill \Delta T_{\mathrm{ano}}& ={\int }_0^{2\mathsf{\pi }}{\left\{\frac{3}{2}\sqrt{\frac{a{\left(1-e^2\right)}^3}{\mathsf{\mu }}}\frac{\Delta a\left(f\right)}{{\left(1+ecosf\right)}^2}-\sqrt{\frac{a^3\left(1-e^2\right)}{\mathsf{\mu }}}\frac{\left[3e+\left(2+e^2\right)cosf\right]}{{\left(1+ecosf\right)}^3}\Delta e\left(f\right)+\frac{r^4}{\mathsf{\mu }p}\left(\frac{\mathrm{d}\omega }{\mathrm{d}t}+cosI\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}t}\right)\right\}}_{\mathrm{K}}df.\hfill \end{array}$$

(50)

The suffix K in Equation (50) means that the content of the curly brackets to which it is appended has to be evaluated on the unperturbed Keplerian ellipse given by Equation (6). Furthermore, the instantaneous shifts ($\Delta \mathsf{\kappa }\left(f\right)$) of $\mathsf{\kappa }=a,e$ are calculated to the first order in the perturbing pK acceleration (A) as follows:

$$\Delta \mathsf{\kappa }\left(f\right)={\int }_{f_0}^{f^{\prime }}\frac{\mathrm{d}\mathsf{\kappa }}{\mathrm{d}{f}^{\prime }}d{f}^{\prime }={\int }_{f_0}^{f^{\prime }}\frac{\mathrm{d}\mathsf{\kappa }}{\mathrm{d}t}\frac{\mathrm{d}t}{\mathrm{d}{f}^{\prime }}d{f}^{\prime },$$

(51)

where $d\mathsf{\kappa }/dt$ are the Gauss equations for the variations of $\mathsf{\kappa }=a,e$ [38,39,40,41,42,43,44,45].

$$\begin{array}{cc}\hfill \frac{\mathrm{d}a}{\mathrm{d}t}& =\frac{2}{n_{\mathrm{K}}\sqrt{1-e^2}}\left[e{A}_{r}sinf+\left(\frac{p}{r}\right)A_{\tau }\right];\hfill \\ \hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}e}{\mathrm{d}t}& =\frac{\sqrt{1-e^2}}{n_{\mathrm{K}}a}\left\{A_{r}sinf+A_{\tau }\left[cosf+\frac{1}{e}\left(1-\frac{r}{a}\right)\right]\right\}.\hfill \end{array}$$

(53)

Finally, $d\mathsf{\Omega }/dt$ y $d\omega /dt$ entering the third term of Equation (50) are the Gauss equations for the variations in the longitude of the ascending node and the argument of pericenter, respectively, given by [38,39,40,41,42,43,44,45]

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}t}& =\frac{1}{n_{\mathrm{K}}asinI\sqrt{1-e^2}}{A}_h\left(\frac{r}{a}\right)sinu,\hfill \\ \hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\omega }{\mathrm{d}t}& =\frac{\sqrt{1-e^2}}{n_{\mathrm{K}}ae}\left[-A_{r}cosf+A_{\tau }\left(1+\frac{r}{p}\right)sinf\right]-cosI\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}t}.\hfill \end{array}$$

(55)

In [4], a variant of the above calculation can be found; in Equation (49), p is adopted as an independent variable, along with the eccentricity (e), and simpler expressions for the partial derivatives of Equation (49) are obtained. The resulting expressions for calculating $\Delta T_{\mathrm{ano}}$ turn out to be

$$\begin{array}{cc}\hfill \Delta T_{\mathrm{ano}}& ={\int }_0^{2\mathsf{\pi }}{\left\{\frac{3}{2}\sqrt{\frac{p}{\mathsf{\mu }}}\frac{\Delta p\left(f\right)}{{\left(1+ecosf\right)}^2}-2\sqrt{\frac{p^3}{\mathsf{\mu }}}\frac{cosf\Delta e\left(f\right)}{{\left(1+ecosf\right)}^3}+\frac{r^4}{\mathsf{\mu }p}\left(\frac{\mathrm{d}\omega }{\mathrm{d}t}+cosI\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}t}\right)\right\}}_{\mathrm{K}}df.\hfill \end{array}$$

(56)

The first-order variation ($\Delta p\left(f\right)$) of the semilatus rectum can be calculated as [52,55]

$$\frac{\mathrm{d}p}{\mathrm{d}f}=\frac{2r^3{A}_{\tau }}{\mathsf{\mu }}.$$

(57)

In the end, both Equations (50) and (56) give the same result.

The presence or absence of the pK anomalistic correction ($\Delta T_{\mathrm{ano}}$) during the orbital period can be intuitively explained as follows. According to

$$\eta =n_{\mathrm{K}}\left(t_0-t_{\mathrm{p}}\right),$$

(58)

where $\eta$ is the mean anomaly in a given epoch, $t_0$ is the initial epoch and $t_{\mathrm{p}}$ is the time of passage at the pericenter. The rate of change of $\eta$ is proportional to the opposite of the pace of variation of $t_{\mathrm{p}}$. Thus, should $\eta$ increase, the crossing of the pericenter position would be anticipated with respect to the Keplerian case, since $t_{\mathrm{p}}$ would decrease and vice-versa. In this case, the variation of $\eta$ would result in an orbit-by-orbit advance or delay of the passages at the pericenter. As will be shown, while the 1 pN gravitoelectric acceleration of Equation (14) does induce a negative rate of $\eta$, the gravitomagnetic LT accelerationof Equation (22) leaves the mean anomaly in the given epoch unchanged. Furthermore, several modified models of gravity inducing radial pK accelerations dependent only on r secularly change both $\omega$ y $\eta$. The Newtonian acceleration raised by the primary’s oblateness ($J_2$) also affects $\eta$, among other things.

3.2. The 1pN Gravitoelectric Correction

The 1pN anomalistic period can be calculated by means of Equations (19)–(21), as explained in Section 3.1. It turns out that

$$T_{\mathrm{ano}}^{1\mathrm{pN}}=T_{\mathrm{K}}+\Delta T_{\mathrm{ano}}^{1\mathrm{pN}},$$

(59)

with

$$\begin{array}{cc}\hfill \Delta T_{\mathrm{ano}}^{1\mathrm{pN}}& =\frac{\mathsf{\pi }\sqrt{{\mathsf{\mu }}_{\mathrm{b}}a}}{2c^2{\left(1-e^2\right)}^2}\left(36+e^2\left(42-38\nu \right)+2e^4\left(6-7\nu \right)-8\nu +3e\left\{\left[28+3e^2\left(4-5\nu \right)-12\nu \right]cos{f}_0\right.\right.\hfill \\ \hfill \\ \hfill & \left.\left.- e\left(-10+8\nu +e\nu cos{f}_0\right)cos2f_0\right\}\right).\hfill \end{array}$$

(60)

In Equation (60), $f_0$ is the true anomaly in the initial epoch. In the test particle limit ($\nu \to 0$), Equation (60) reduces to

$$\Delta T_{\mathrm{ano}}^{1\mathrm{pN}}=\frac{3\pi \sqrt{\mathsf{\mu }a}}{c^2{\left(1-e^2\right)}^2}\left[6+7e^2+2e^4+2e\left(7+3e^2\right)cos{f}_0+5e^{2}cos2f_0\right].$$

(61)

Figure 1 was obtained for generic values of the Keplerian orbital elements of a test particle revolving around a massive primary. It confirms the analytical result of Equation (61); over, say, three orbital revolutions, the satellite always reaches the precessing line of apsides after a time interval equal to $T_{\mathrm{ano}}^{1\mathrm{pN}}$. It is longer than $T_{\mathrm{K}}$, in agreement with Equation (61), which is always positive.

Furthermore, Figure 2 plots the final part of the time series of the cosine ($\widehat{\mathit{r}}·\widehat{\mathit{C}}$) of the angle between the position vector ($\mathit{r}$) and the Laplace–Runge–Lenz unit vector ($\widehat{\mathit{C}}$) versus time (t) in units of $T_{\mathrm{K}}$ for a numerically integrated fictitious test particle with and without Equation (14) starting from the moving pericenter in both cases, i.e., for ${\widehat{\mathit{r}}}_0·{\widehat{\mathit{C}}}_0=+1$. It can be seen that the orbiter comes back to the same position on the precessing line of apsides, i.e., it is $\widehat{\mathit{r}}·\widehat{\mathit{C}}=+1$ again, just after $T_{\mathrm{ano}}^{1\mathrm{pN}}=T_{\mathrm{K}}+\Delta T_{\mathrm{ano}}^{1\mathrm{pN}}$, differing from $T_{\mathrm{K}}$ by a positive amount, in agreement with Equation (61).

3.3. The 1pN Gravitomagnetic Lense–Thirring Correction

The LT anomalistic period can be calculated by means of Equations (30)–(32), as explained in Section 3.1. It turns out to be

$$T_{\mathrm{ano}}^{\mathrm{LT}}=T_{\mathrm{K}}+\Delta T_{\mathrm{ano}}^{\mathrm{LT}},$$

(62)

with

$$\Delta T_{\mathrm{ano}}^{\mathrm{LT}}=0;$$

(63)

this is an exact result that is valid for all orders in the eccentricity (e).

Figure 3, obtained for generic values of the Keplerian orbital parameters, shows just that; over three orbital revolutions, the test particle always reaches the precessing line of apsides after a time interval equal to the Keplerian orbital period after each orbit.

Furthermore, Figure 4 plots the final part of the time series of the cosine ($\widehat{\mathit{r}}·\widehat{\mathit{C}}$) of the angle between the position vector ($\mathit{r}$) and the Laplace–Runge–Lenz vector ($\mathit{C}$) versus time (t) in units of $T_{\mathrm{K}}$ for a numerically integrated fictitious test particle acted upon by Equation (22) starting from, the moving pericenter, i.e., for ${\widehat{\mathit{r}}}_0·{\widehat{\mathit{C}}}_0=+1$. It can be seen that it comes back to the same position on the precessing line of apsides, i.e., it is $\widehat{\mathit{r}}·\widehat{\mathit{C}}=+1$ again after just one Keplerian orbital period.

The fact that the gravitomagnetic apsidal period is identical to the Keplerian one can be intuitively justified, since there is no net shift per orbit of the mean anomaly at epoch $\eta$. Indeed, from Equation (58), one infers that $\eta$ is proportional to $t_{\mathrm{p}}$ through $n_{\mathrm{K}}$. Thus, since the latter stays constant because a is not secularly affected by the gravitomagnetic field, the rate of change of the mean anomaly in a given epoch is proportional to the opposite of the pace of variation of the time of passage at the pericenter according to

$$\frac{\mathrm{d}\eta }{\mathrm{d}t}=-n_{\mathrm{K}}\frac{\mathrm{d}{t}_{\mathrm{p}}}{\mathrm{d}t}.$$

(64)

Should $\eta$ increase, the crossing of the pericenter would be anticipated with respect to the Keplerian case, since $t_{\mathrm{p}}$ would decrease and vice-versa. In this case, the variation of $\eta$ would result in an orbit-by-orbit advance or delay of the passages at the pericenter, which does not occur in the present case because, in fact, ${〈d\eta /dt〉}^{\mathrm{LT}}=0$.

3.4. The Newtonian Quadrupolar Correction

The $J_2$-affected anomalistic period can be calculated by means of Equations (41)–(43), as explained in Section 3.1. It turns out to be

$$T_{\mathrm{ano}}^{J_2}=T_{\mathrm{K}}+\Delta T_{\mathrm{ano}}^{J_2},$$

(65)

with

$$\begin{array}{cc}\hfill \Delta T_{\mathrm{ano}}^{J_2}& =\frac{3\mathsf{\pi }{J}_2{R}_{\mathrm{e}}^2{\left(1+ecos{f}_0\right)}^3}{2{\left(1-e^2\right)}^3\sqrt{\mathsf{\mu }a}}\left(-2{\widehat{T}}_1+3{\widehat{T}}_2+3{\widehat{T}}_{3}cos2u_0+6{\widehat{T}}_{6}sin2u_0\right).\hfill \end{array}$$

(66)

Figure 5, obtained for generic values of the Keplerian orbital elements, confirms the analytical result of Equation (66); over three orbital revolutions, the test particle always reaches the precessing line of apsides after a time interval equal to $T_{\mathrm{ano}}^{J2}$ for each orbit. For the particular choice of the values of the primary’s spin and orbital parameters, it turns out to be longer than $T_{\mathrm{K}}$, in agreement with Equation (66).

Furthermore, Figure 6 plots the final part of the time series of the cosine ($\widehat{\mathit{r}}·\widehat{\mathit{C}}$) of the angle between the position vector ($\mathit{r}$) and the Laplace–Runge–Lenz unit vector ($\widehat{\mathit{C}}$) versus time (t) in units of $T_{\mathrm{K}}$ for a numerically integrated fictitious test particle with and without Equation (34) starting from the moving pericenter in both cases, i.e., for ${\widehat{\mathit{r}}}_0·{\widehat{\mathit{C}}}_0=+1$. It can be seen that it comes back to the same position on the precessing line of apsides, i.e., it is $\widehat{\mathit{r}}·\widehat{\mathit{C}}=+1$ again, just after $T_{\mathrm{ano}}^{J_2}=T_{\mathrm{K}}+\Delta T_{\mathrm{ano}}^{J_2}$, differing from $T_{\mathrm{K}}$ by a (positive) amount, in agreement with Equation (66) for the particular choice of the generic values of the spin and the orbital parameters adopted in the numerical integrations.

4. The Draconitic Period

4.1. General Calculational Scheme

For a perturbed trajectory, the draconitic period ($T_{\mathrm{dra}}$) is defined as the time interval between two successive instants when the real position of the test particle coincides with the ascending position on the corresponding osculating ellipse.

It can be calculated as [4,57]

$$T_{\mathrm{dra}}=T_{\mathrm{K}}+\Delta T_{\mathrm{dra}}={\int }_0^{2\mathsf{\pi }}\left(\frac{\mathrm{d}t}{\mathrm{d}u}\right)du,$$

(67)

where $dt/du$, when a pK perturbing acceleration $(\mathit{A}$) is present, can be obtained as follows.

From Equations (4) and (48), one obtains [57,58]

$$\frac{\mathrm{d}u}{\mathrm{d}t}=\frac{\sqrt{\mathsf{\mu }p}}{r^2}\left(1-\frac{r^{2}cosI}{\sqrt{\mathsf{\mu }p}}\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}t}\right).$$

(68)

Then, it is

$$\frac{\mathrm{d}t}{\mathrm{d}u}\simeq \frac{r^2}{\sqrt{\mathsf{\mu }p}}+\frac{r^{4}cosI}{\mathsf{\mu }p}\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}t}.$$

(69)

Note that $d\mathsf{\Omega }/dt$ is already expressed in terms of u, as per Equation (54). Nonsingular orbital elements can be used, which are the components of the eccentricity vector [59], an alternative formulation of the Laplace–Runge–Lenz vector [60]. In the context of pulsar astronomy, they are also known as first and second Laplace–Lagrange parameters (${ϵ}_1,{ϵ}_2$) [57,61,62].

$$\begin{array}{cc}\hfill k& :=esin\omega ,\hfill \\ \hfill \end{array}$$

(70)

$$\begin{array}{cc}\hfill q& :=ecos\omega ,\hfill \end{array}$$

(71)

Equation (6) can be cast into the following form:

$$r=\frac{p}{1+qcosu+ksinu},$$

(72)

in which p, $q$ y $k$ are treated as independent variables. By proceeding as in Section 3.1, the following can be obtained [4,57]:

$$\begin{array}{cc}\hfill \Delta T_{\mathrm{dra}}& ={\int }_0^{2\mathsf{\pi }}{\left\{\frac{3}{2}\sqrt{\frac{p}{\mathsf{\mu }}}\frac{\Delta p\left(u\right)}{{\left(1+qcosu+ksinu\right)}^2}-2\sqrt{\frac{p^3}{\mathsf{\mu }}}\frac{cosu\Delta q\left(u\right)+sinu\Delta k\left(u\right)}{{\left(1+qcosu+ksinu\right)}^3}+\frac{r^{4}cosI}{\mathsf{\mu }p}\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}t}\right\}}_{\mathrm{K}}du.\hfill \end{array}$$

(73)

The first-order variations ($\Delta p\left(u\right)$, $\Delta q\left(u\right)$ y $\Delta k\left(u\right)$) entering Equation (73) are worked out by integrating the following expressions from $u_0$ to u [57]:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}p}{\mathrm{d}u}& =\frac{2r^3{A}_{\tau }}{\mathsf{\mu }},\hfill \\ \hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}q}{\mathrm{d}u}& =\frac{r^{2}sinu{A}_r}{\mathsf{\mu }}+\frac{r^2\left[rq+\left(r+p\right)cosu\right]A_{\tau }}{\mathsf{\mu }}+\frac{cotI{r}^{3}ksinu{A}_h}{\mathsf{\mu }p},\hfill \\ \hfill \end{array}$$

(75)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}k}{\mathrm{d}u}& =-\frac{r^{2}cosu{A}_r}{\mathsf{\mu }}+\frac{r^2\left[rk+\left(r+p\right)sinu\right]A_{\tau }}{\mathsf{\mu }}-\frac{cotI{r}^{3}qsinu{A}_h}{\mathsf{\mu }p}\hfill \end{array}$$

(76)

As far as the actual measurability of the draconitic period in some astronomical scenario of interest is concerned, it was shown [63,64,65] that it is possible to measure it for an artificial Earth satellite. In their analyses, Amelin [63], Kassimenko [64] and Zhongolovich [65] used Soviet satellite 1960 $\epsilon$ 3 as the ratio of the difference in the times of passages of the sub-satellite point through a chosen parallel for two following epochs to the number of satellite revolutions corresponding to this difference. The accuracy reached at that time should be of the order of $\simeq {10}^{-4}$ s [64]; it is not unlikely that it could be improved by orders of magnitude with the most recent and currently available techniques.

4.2. The 1pN Gravitoelectric Correction

The 1pN draconitic period can be calculated by means of Equations (19)–(21), as explained in Section 4.1. It turns out to be

$$T_{\mathrm{dra}}^{1\mathrm{pN}}=T_{\mathrm{K}}+\Delta T_{\mathrm{dra}}^{1\mathrm{pN}},$$

(77)

with

$$\begin{array}{cc}\hfill \Delta T_{\mathrm{dra}}^{1\mathrm{pN}}& =\frac{\mathsf{\pi }\sqrt{{\mathsf{\mu }}_{\mathrm{b}}a}}{4c^2}\left(72+e^2\left(84-76\nu \right)+4e^4\left(6-7\nu \right)-16\nu -3e\left\{\left[8\left(-7+3\nu \right)+e^2\left(-24+31\nu \right)\right]cos{f}_0\right.\right.\hfill \\ \hfill \\ \hfill & \left.\left.+ e\left[4\left(-5+4\nu \right)cos2f_0+e\nu cos3f_0\right]\right\}-\frac{24\sqrt{1-e^2}}{{\left(1+ecos\omega \right)}^2}\right).\hfill \end{array}$$

(78)

In the test particle limit ($\nu \to 0$), Equation (78) reduces to

$$\begin{array}{cc}\hfill \Delta T_{\mathrm{dra}}^{1\mathrm{pN}}& =\frac{3\mathsf{\pi }\sqrt{\mathsf{\mu }a}}{c^2}\left[\frac{6+7e^2+2e^4+2e\left(7+3e^2\right)cos{f}_0+5e^{2}cos2f_0}{{\left(1-e^2\right)}^2}-\frac{2\sqrt{1-e^2}}{{\left(1+ecos\omega \right)}^2}\right].\hfill \end{array}$$

(79)

It can be noted that Equation (79) is always positive for all values of $e,f_0$ y $\omega$; thus, the node is reached later than in the Keplerian case.

Figure 7, obtained for generic values of the Keplerian orbital elements, confirms the analytical result of Equation (79); over three orbital revolutions, the test particle always reaches the fixed line of nodes after a time interval equal to $T_{\mathrm{dra}}^{1\mathrm{pN}}$. It is longer than $T_{\mathrm{K}}$, in agreement with Equation (79).

Furthermore, Figure 8 plots the final part of the time series of the cosine ($\widehat{\mathit{r}}·\widehat{\mathit{l}}$) of the angle between the position vector ($\mathit{r}$) and the node unit vector ($\widehat{\mathit{l}}$) versus time (t) in units of ($T_{\mathrm{K}}$) for a numerically integrated fictitious test particle with and without Equation (14) starting from the fixed ascending node in both cases, i.e., for ${\widehat{\mathit{r}}}_0·{\widehat{\mathit{l}}}_0=+1$. It can be seen that it comes back to the same position on the constant line of nodes, i.e., it is $\widehat{\mathit{r}}·\widehat{\mathit{l}}=+1$ again, just after $T_{\mathrm{dra}}^{1\mathrm{pN}}=T_{\mathrm{K}}+\Delta T_{\mathrm{dra}}^{1\mathrm{pN}}$, differing from $T_{\mathrm{K}}$ by a positive amount, in agreement with Equation (79).

4.3. The 1 pN Gravitomagnetic Lense–Thirring Correction

The LT draconitic period can be calculated by means of Equations (30)–(32), as explained in Section 4.1. It turns out to be

$$T_{\mathrm{dra}}^{\mathrm{LT}}=T_{\mathrm{K}}+\Delta T_{\mathrm{dra}}^{\mathrm{LT}},$$

(80)

with

$$\Delta T_{\mathrm{dra}}^{\mathrm{LT}}=\frac{4\mathsf{\pi }J\left(2\mathtt{Jh}+\mathtt{Jm}cotI\right)}{c^{2}M{\left(1+ecos\omega \right)}^2}.$$

(81)

The explicit form of the geometric coefficient in the numerator of Equation (81) depending on the orientation in space of both the orbital plane and the primary’s spin axis is

$$\begin{array}{cc}\hfill 2\mathtt{Jh}+\mathtt{Jm}cotI& =3cosIsin\delta +cos\delta \left(cscI-3sinI\right)sin\left({\alpha }_J-\mathsf{\Omega }\right).\hfill \end{array}$$

(82)

In general, it can be either positive and negative. For a polar orbit, i.e., for $\mathsf{\Omega }={\alpha }_J$ y $I={90}^{\circ }$, the gravitomagnetic correction vanishes, as per Equation (82). Instead, for an equatorial orbit arbitrarily oriented in space, it does not vanish, amounting to

$$\Delta T_{\mathrm{dra}}^{\mathrm{LT}}=\pm \frac{8\mathsf{\pi }J}{c^{2}M{\left(1+ecos\omega \right)}^2}.$$

(83)

Furthermore, for circular orbits, Equation (83) reduces to

$$\Delta T_{\mathrm{dra}}^{\mathrm{LT}}=\pm \frac{8\mathsf{\pi }J}{c^{2}M}.$$

(84)

If the orbital plane lies in the reference plane, i.e., for $I=0$, Equation (81) loses its meaning as expected, since in this case, the line of nodes is no longer defined.

Figure 9, obtained for generic values of the Keplerian orbital parameters, confirms the analytical result of Equation (81); over three orbital revolutions, the test particle always reaches the precessing line of nodes after a time interval equal to $T_{\mathrm{dra}}^{\mathrm{LT}}$ after each orbit. For the particular choice of the values of the primary and orbital parameters, it turns out to be longer than $T_{\mathrm{K}}$, in agreement with Equation (81).

Furthermore, Figure 10 plots the final part of the time series of the cosine ($\widehat{\mathit{r}}·\widehat{\mathit{l}}$) of the angle between the position vector ($\mathit{r}$) and the node unit vector ($\widehat{\mathit{l}}$) versus time (t) in units of $T_{\mathrm{K}}$ for a numerically integrated fictitious test particle with and without Equation (22) starting from the moving ascending node in both cases, i.e., for ${\widehat{\mathit{r}}}_0·{\widehat{\mathit{l}}}_0=+1$. It can be seen that it comes back to the same position on the precessing line of nodes, i.e., it is $\widehat{\mathit{r}}·\widehat{\mathit{l}}=+1$ again, just after $T_{\mathrm{dra}}^{\mathrm{LT}}=T_{\mathrm{K}}+\Delta T_{\mathrm{dra}}^{\mathrm{LT}}$, differing from $T_{\mathrm{K}}$ by a positive amount, in agreement with Equation (81) for the particular choice of the generic values of the spin and the orbital parameters adopted in the numerical integrations.

4.4. The Newtonian Quadrupolar Correction

The $J_2$-affected draconitic period can be calculated by means of Equations (41)–(43), as explained in Section 4.1. It turns out to be

$$T_{\mathrm{dra}}^{J_2}=T_{\mathrm{K}}+\Delta T_{\mathrm{dra}}^{J_2},$$

(85)

with

$$\begin{array}{cc}\hfill \Delta T_{\mathrm{dra}}^{J_2}& =\frac{3\mathsf{\pi }{J}_2{R}_{\mathrm{e}}^2}{2\sqrt{\mathsf{\mu }a\left(1-e^2\right)}}\left[\frac{1}{{\left(1+ecos\omega \right)}^2}\left(-2{\widehat{T}}_1+3{\widehat{T}}_2-2{\widehat{T}}_{5}cotI\right)+\frac{{\left(1+ecos{f}_0\right)}^3}{{\left(1-e^2\right)}^{5/2}}\left(-2{\widehat{T}}_1+3{\widehat{T}}_2+3{\widehat{T}}_{3}cos2u_0+6{\widehat{T}}_{6}sin2u_0\right)\right].\hfill \end{array}$$

(86)

It can be noted that Equation (86) is not defined for $I\to 0$ because of the term

$$\begin{array}{cc}\hfill {\widehat{T}}_{5}cotI& =cotI\left[sinIsin{\delta }_J+cosIcos{\delta }_{J}sin\left({\alpha }_J-\mathsf{\Omega }\right)\right]\left[cosIsin{\delta }_J-cos{\delta }_{J}sinIsin\left({\alpha }_J-\mathsf{\Omega }\right)\right],\hfill \end{array}$$

as it is expected, since in this case, the line of nodes is no longer defined.

Figure 11, obtained for generic values of the Keplerian orbital elements, confirms the analytical result of Equation (86); over three orbital revolutions, the test particle always reaches the precessing line of nodes after a time interval equal to $T_{\mathrm{dra}}^{J_2}$ after each orbit. For the particular choice of the values of the primary’s spin and orbital parameters, it is shorter than $T_{\mathrm{K}}$, in agreement with Equation (86).

Furthermore, Figure 12 plots the final part of the time series of the cosine ($\widehat{\mathit{r}}·\widehat{\mathit{l}}$) of the angle between the position vector ($\mathit{r}$) and the node unit vector ($\widehat{\mathit{l}}$) versus time (t) in units of $T_{\mathrm{K}}$ for a numerically integrated fictitious test particle with and without Equation (34) starting from the moving ascending node in both cases, i.e., for ${\widehat{\mathit{r}}}_0·{\widehat{\mathit{l}}}_0=+1$. It can be seen that it comes back to the same position on the precessing line of nodes, i.e., it is $\widehat{\mathit{r}}·\widehat{\mathit{l}}=+1$ again, just after $T_{\mathrm{dra}}^{J_2}=T_{\mathrm{K}}+\Delta T_{\mathrm{dra}}^{J_2}$, differing from $T_{\mathrm{K}}$ by a positive amount, in agreement with Equation (86) for the particular choice of the generic values of the spin and the orbital parameters adopted in the numerical integrations.

5. The Sidereal Period

5.1. General Calculational Scheme

In general, both the line of nodes and the line of apsides vary over time because of one or more pK accelerations. Thus, it may be useful to look at a characteristic orbital timescale involving the crossing of some fixed reference direction in space; the sidereal period ($T_{\mathrm{sid}}$), defined as the time interval between two successive instants when the real position of the test particle lies in a given reference direction, plays such a role.

For an arbitrarily inclined orbit, the sidereal period can be calculated as

$$T_{\mathrm{sid}}=T_{\mathrm{K}}+\Delta T_{\mathrm{sid}}={\int }_0^{2\mathsf{\pi }}\left(\frac{\mathrm{d}t}{\mathrm{d}\varphi }\right)d\varphi ,$$

(87)

where $\varphi \left(t\right)$ is the azimuthal angle reckoned from the reference x-axis in the fundamental plane; when the latter is assumed to be coincident with the Earth’s equatorial plane in some reference epoch, $\varphi \left(t\right)$ is the right ascension ($\alpha \left(t\right)$) of the celestial body of interest. From

$$\begin{array}{cc}\hfill x\left(t\right)& =r\left(t\right)\left[cos\mathsf{\Omega }cosu\left(t\right)-cosIsin\mathsf{\Omega }sinu\left(t\right)\right],\hfill \\ \hfill \end{array}$$

(88)

$$\begin{array}{cc}\hfill y\left(t\right)& =r\left(t\right)\left[sin\mathsf{\Omega }cosu\left(t\right)+cosIcos\mathsf{\Omega }sinu\left(t\right)\right],\hfill \end{array}$$

(89)

one obtains $\varphi \left(t\right)$ as

$$\varphi \left(t\right)=arctan\left[\frac{y\left(t\right)}{x\left(t\right)}\right];$$

(90)

it is a function of the generally varying $I\left(t\right)$, $\mathsf{\Omega }\left(t\right)$, and of $u\left(t\right)$, i.e., $\varphi \left(t\right)=\varphi \left(I\left(t\right),\mathsf{\Omega }\left(t\right),u\left(t\right)\right)$. Since the ongoing calculation is to the first order in the pK acceleration, the differential $d\varphi$ in Equation (87) can be written as

$$d\varphi \simeq \left(\frac{\partial \varphi }{\partial u}\right)du.$$

(91)

The integrand of Equation (87) can be obtained as

$$\begin{array}{cc}\hfill \frac{\mathrm{d}t}{\mathrm{d}\varphi }& =\frac{1}{\frac{\mathrm{d}\varphi }{\mathrm{d}t}}=\frac{1}{\frac{\partial \varphi }{\partial I}\frac{\mathrm{d}I}{\mathrm{d}t}+\frac{\partial \varphi }{\partial \mathsf{\Omega }}\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}t}+\frac{\partial \varphi }{\partial u}\frac{\mathrm{d}u}{\mathrm{d}t}}=\frac{1}{\frac{\mathrm{d}u}{\mathrm{d}t}\frac{\partial \varphi }{\partial u}\left[1+\frac{\partial u}{\partial \varphi }\left(\frac{\partial \varphi }{\partial I}\frac{\mathrm{d}I}{\mathrm{d}u}+\frac{\partial \varphi }{\partial \mathsf{\Omega }}\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}u}\right)\right]}.\hfill \end{array}$$

(92)

Thus, to the first order in the pK acceleration, the integral of Equation (87) can be approximated as

$$\begin{array}{cc}\hfill T_{\mathrm{sid}}& \simeq {\int }_0^{2\mathsf{\pi }}\frac{\mathrm{d}t}{\mathrm{d}u}\left[1-\frac{\partial u}{\partial \varphi }\left(\frac{\partial \varphi }{\partial I}\frac{\mathrm{d}I}{\mathrm{d}u}+\frac{\partial \varphi }{\partial \mathsf{\Omega }}\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}u}\right)\right]du={\int }_0^{2\mathsf{\pi }}\left(\frac{\mathrm{d}t}{\mathrm{d}u}\right)du-{\int }_0^{2\mathsf{\pi }}\frac{1}{\frac{\partial \varphi }{\partial u}}\left(\frac{\partial \varphi }{\partial I}\frac{\mathrm{d}I}{\mathrm{d}u}+\frac{\partial \varphi }{\partial \mathsf{\Omega }}\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}u}\right)\left(\frac{\mathrm{d}t}{\mathrm{d}u}\right)du.\hfill \end{array}$$

(93)

The first term in Equation (93) is nothing but the draconitic period and can be calculated to the order $\mathcal{O}\left(A\right)$, as outlined in Section 4.1. The second term in Equation (93) is a correction to the former

$$\Delta T_{\mathrm{sid}\mathrm{II}}:=-{\int }_0^{2\mathsf{\pi }}\frac{1}{\frac{\partial \varphi }{\partial u}}\left(\frac{\partial \varphi }{\partial I}\frac{\mathrm{d}I}{\mathrm{d}u}+\frac{\partial \varphi }{\partial \mathsf{\Omega }}\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}u}\right){\left(\frac{\mathrm{d}t}{\mathrm{d}u}\right)}_{\mathrm{K}}du,$$

(94)

taking into account the fact that, in general, the orbital plane is displaced by the pK acceleration; indeed, the rates of I y $\mathsf{\Omega }$ enter it. In Equation (94), the suffix K appended to $dt/du$ implies that it has to be calculated on the unperturbed Keplerian ellipse in order to keep the calculation to the first order in A.

If the orbital plane coincides with the fundamental one, the previously outlined calculational strategy may lead to analytical expressions for $T_{\mathrm{sid}}$ that, for some pK accelerations, are singular in $I=0$. In such cases, the sidereal period can be straightforwardly calculated by means of the true longitude.

$$l:=\varpi +f,$$

(95)

where

$$\varpi :=\mathsf{\Omega }+\omega$$

(96)

is the longitude of the pericenter. It is a dogleg angle [66,67], since $\mathsf{\Omega }$ y $\omega$ generally lie in different planes [4] as

$$T_{\mathrm{sid}}=T_{\mathrm{K}}+\Delta T_{\mathrm{sid}}={\int }_0^{2\mathsf{\pi }}\left(\frac{\mathrm{d}t}{\mathrm{d}l}\right)dl,$$

(97)

in close analogy with Section 3.1 y Section 4.1. It should be recalled that l is generally a dogleg angle, since $\mathsf{\Omega }$ y u are located in different planes; it is the true longitude of the test particle actually moving along its real orbit only if $I=0$. When a pK perturbing acceleration ($\mathit{A}$) enters the equations of motion, $dt/dl$ can be obtained in the following way.

According to Equations (47) and (111),

$$\frac{\mathrm{d}l}{\mathrm{d}t}=\frac{\sqrt{\mathsf{\mu }p}}{r^2}\left[1+\frac{2r^2{sin}^2\left(I/2\right)}{\sqrt{\mathsf{\mu }p}}\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}t}\right].$$

(98)

Then, it can be written as

$$\frac{\mathrm{d}t}{\mathrm{d}l}\simeq \frac{r^2}{\sqrt{\mathsf{\mu }p}}-\frac{2r^4{sin}^2\left(I/2\right)}{\mathsf{\mu }p}\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}t}.$$

(99)

The sine of the argument of latitude (u) entering Equation (54) for $d\mathsf{\Omega }/dt$ can be written in terms of l as $sin\left(l-\mathsf{\Omega }\right)$. By introducing nonsingular equinoctial elements [68]

$$\begin{array}{cc}\hfill Q& :=ecos\varpi ,\hfill \\ \hfill \end{array}$$

(100)

$$\begin{array}{cc}\hfill K& :=esin\varpi ,\hfill \end{array}$$

(101)

Equation (6) can be rewritten as

$$r=\frac{p}{1+Qcosl+Ksinl},$$

(102)

in which p, $Q$ y $K$ are treated as independent variables. By proceeding as in Section 3.1 y Section 4.1, one obtains [4]

$$\begin{array}{cc}\hfill \Delta T_{\mathrm{sid}}& ={\int }_0^{2\mathsf{\pi }}{\left\{\frac{3}{2}\sqrt{\frac{p}{\mathsf{\mu }}}\frac{\Delta p\left(l\right)}{{\left(1+Qcosl+Ksinl\right)}^2}-2\sqrt{\frac{p^3}{\mathsf{\mu }}}\frac{cosl\Delta Q\left(l\right)+sinl\Delta K\left(l\right)}{{\left(1+Qcosl+Ksinl\right)}^3}-\frac{2r^4{sin}^2\left(I/2\right)}{\mathsf{\mu }p}\frac{\mathrm{d}\mathsf{\Omega }}{\mathrm{d}t}\right\}}_{\mathrm{K}}dl.\hfill \end{array}$$

(103)

The first-order variations $(\Delta p\left(l\right)$, $\Delta Q\left(l\right)$ y $\Delta K\left(l\right)$) entering Equation (103) can be worked out by integrating the following expressions from $l_0$ to l [4]:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}p}{\mathrm{d}l}& =\frac{2r^3{A}_{\tau }}{\mathsf{\mu }},\hfill \\ \hfill \end{array}$$

(104)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}Q}{\mathrm{d}l}& =\frac{r^{2}sinl{A}_r}{\mathsf{\mu }}+\frac{r^2\left[rQ+\left(r+p\right)cosl\right]A_{\tau }}{\mathsf{\mu }}-\frac{tan\left(I/2\right)r^{3}Ksin\left(l-\mathsf{\Omega }\right)A_h}{\mathsf{\mu }p},\hfill \\ \hfill \end{array}$$

(105)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}K}{\mathrm{d}l}& =-\frac{r^{2}cosl{A}_r}{\mathsf{\mu }}+\frac{r^2\left[rK+\left(r+p\right)sinl\right]A_{\tau }}{\mathsf{\mu }}+\frac{tan\left(I/2\right)r^{3}Qsin\left(l-\mathsf{\Omega }\right)A_h}{\mathsf{\mu }p}\hfill \end{array}$$

(106)

If the orbital plane is aligned with the fundamental one, Equations (103), (105) and (106) have to be calculated with $I=0$.

It is generally expected that if the orbital plane stays constant in space, i.e, if neither the nodes, when defined, nor the orbit’s projection onto the fundamental plane change over time, the sidereal period coincides with the draconitic one, since the line of nodes is oriented in a fixed direction in space.

5.2. The 1 pN Gravitoelectric Correction

As shown in Section 5.1, the sidereal period for a generic perturbed orbit is the sum of the draconitic period, calculated as explained in Section 4.1, and the term given by Equation (94). For Equations (14) and (94), it turns out to be

$$\Delta T_{\mathrm{sid}\mathrm{II}}^{1\mathrm{pN}}=0.$$

(107)

Thus, in this case, the sidereal period coincides with the draconitic one.

This is shown in Figure 13, which plots the final part of the time series of the cosine of the angle ($\varphi$) normalized to its initial value ($cos{\varphi }_0$) versus time (t) in units of $T_{\mathrm{K}}$ for a numerically integrated fictitious test particle with and without Equation (14) starting from the same generic initial position. It can be seen that it comes back to the same position in the fixed direction chosen in the reference plane, i.e., it is $cos\varphi /cos{\varphi }_0=+1$ again, just after $T_{\mathrm{sid}}^{1\mathrm{pN}}=T_{\mathrm{dra}}^{1\mathrm{pN}}$, differing from $T_{\mathrm{K}}$ by a positive amount, in agreement with Equation (79).

5.3. The 1 pN Gravitomagnetic Lense–Thirring Correction

As shown in Section 5.1, the sidereal period for a generic perturbed orbit is the sum of the draconitic period, calculated as explained in Section 4.1, and the term given by Equation (94). Equations (22) and (94) yield

$$\begin{array}{cc}\hfill \Delta T_{\mathrm{sid}\mathrm{II}}^{\mathrm{LT}}& =\frac{4\mathsf{\pi }\mathit{J}cotI}{c^{2}M{e}^2\sqrt{1-e^2}}·\left\{\widehat{\mathit{m}}\left[-e^2+2\left(2-e^2-2\sqrt{1-e^2}\right)cos2\omega \right]+2\widehat{\mathit{l}}\left(-2+e^2+2\sqrt{1-e^2}\right)sin2\omega \right\}.\hfill \end{array}$$

(108)

In the equatorial case, the orbital plane stays constant in space, Equation (108) vanishes and the sidereal period coincides with the draconitic one, as expected, since neither the line of nodes nor the orbit’s projection onto the reference plane change. By taking the sum of Equations (81) and (108) into consideration, the full expression of the gravitomagnetic correction of the sidereal period ($\Delta T_{\mathrm{sid}}^{\mathrm{LT}}$) is obtained. It can be noted that for a generic eccentric orbit, $\Delta T_{\mathrm{sid}}^{\mathrm{LT}}$ is not defined if the orbital plane lies in the fundamental one. Nonetheless, for $e=0$, it reduces to

$$\Delta T_{\mathrm{sid}}^{\mathrm{LT}}=\frac{8\mathsf{\pi }J}{c^{2}M}\left[cosIsin{\delta }_J-cos{\delta }_{J}sinIsin\left({\alpha }_J-\mathsf{\Omega }\right)\right],$$

(109)

which is not singular when $I=0$. By using the true longitude (l) in the case of $I=0$, it turns out that

$$\Delta T_{\mathrm{sid}}^{\mathrm{LT}}=\frac{8\mathsf{\pi }Jsin{\delta }_J}{c^{2}M{\left(1+ecos\varpi \right)}^2}.$$

(110)

In the limit ($e\to 0$), it reduces to

$$\Delta T_{\mathrm{sid}}^{\mathrm{LT}}=\frac{8\mathsf{\pi }Jsin{\delta }_J}{c^{2}M},$$

(111)

which agrees with Equation (109) calculated with $I=0$. In turn, if ${\delta }_J=\pm {90}^{\circ }$, corresponding to the case of an equatorial orbit whose orbital plane coincides with the reference plane, Equation (111) becomes

$$\Delta T_{\mathrm{sid}}^{\mathrm{LT}}=\pm \frac{8\mathsf{\pi }J}{c^{2}M},$$

(112)

in agreement with Equation (84).

Figure 14 confirms the analytical results of Equations (81) and (108). Indeed, over three orbital revolutions, the projection of a generic LT perturbed orbit in the fundamental plane ($\left\{x,y\right\}$) crosses a fixed direction in the latter set by a certain value (${\varphi }_0$), always after a time interval equal to $T_{\mathrm{sid}}^{\mathrm{LT}}=T_{\mathrm{dra}}^{\mathrm{LT}}+\Delta T_{\mathrm{sid}\mathrm{II}}^{\mathrm{LT}}$ for each orbit. With the particular choice of the primary’s spin and the orbital parameters used in the picture, $T_{\mathrm{sid}}^{\mathrm{LT}}$ turns out to be shorter than $T_{\mathrm{K}}$, in agreement with Equations (81) and (108).

Furthermore, Figure 15 plots the final part of the time series of the cosine of the $\varphi$ angle normalized to its initial value ($cos{\varphi }_0$) versus time (t) in units of $T_{\mathrm{K}}$ for a numerically integrated fictitious test particle with and without Equation (22) starting from the same generic initial position. It can be seen that it comes back to the same position in the fixed direction chosen in the reference plane, i.e., it is $cos\varphi /cos{\varphi }_0=+1$ again, just after $T_{\mathrm{sid}}^{\mathrm{LT}}=T_{\mathrm{dra}}^{\mathrm{LT}}+\Delta T_{\mathrm{sid}\mathrm{II}}^{\mathrm{LT}}$, differing from $T_{\mathrm{K}}$ by a positive amount, in agreement with Equations (81) and (108) for the particular choice of the generic values of the spin and the orbital parameters adopted in the numerical integrations.

5.4. The Newtonian Quadrupolar Correction

As shown in Section 5.1, the sidereal period for a generic perturbed orbit is the sum of the draconitic period, calculated as explained in Section 4.1, and the term given by Equation (94). For Equations (34) and (94), it turns out to be

$$\begin{array}{cc}\hfill \Delta T_{\mathrm{sid}\mathrm{II}}^{J_2}& =\frac{3\mathsf{\pi }{J}_2{R}_{\mathrm{e}}^{2}cotI}{e^2\sqrt{\mathsf{\mu }a\left(1-e^2\right)}}\left\{{\widehat{T}}_5\left[e^2+2\left(-2+e^2+2\sqrt{1-e^2}\right)cos2\omega \right]-2{\widehat{T}}_4\left(-2+e^2+2\sqrt{1-e^2}\right)sin2\omega \right\}.\hfill \end{array}$$

(113)

For equatorial orbits, Equation (113) vanishes, and the sidereal period reduces to the draconitic one. The oblateness of the sidereal period ($\Delta T_{\mathrm{sid}}^{J_2}$) can be obtained by summing Equations (86) and (113); for an elliptic orbit, it turns out to be singular when $I=0$. Instead, in the limit ($e\to 0$), it reduces to

$$\begin{array}{cc}\hfill \Delta T_{\mathrm{sid}}^{J_2}& =\frac{3\mathsf{\pi }{J}_2{R}_{\mathrm{e}}^2}{2\sqrt{\mathsf{\mu }a}}\left\{-4+6{cos}^2{\delta }_J{cos}^2\left({\alpha }_J-\mathsf{\Omega }\right)+6cos{\delta }_{J}cos\left({\alpha }_J-\mathsf{\Omega }\right)sin2u_0\left[sinIsin{\delta }_J+cosIcos{\delta }_{J}sin\left({\alpha }_J-\mathsf{\Omega }\right)\right]\right.\hfill \\ \hfill & + 6{\left[sinIsin{\delta }_J+cosIcos{\delta }_{J}sin\left({\alpha }_J-\mathsf{\Omega }\right)\right]}^2+3cos2u_0\left[cos{\delta }_{J}cos\left({\alpha }_J-\mathsf{\Omega }\right)-sinIsin{\delta }_J-cosIcos{\delta }_{J}sin\left({\alpha }_J-\mathsf{\Omega }\right)\right]\hfill \\ & \left.\times \left[cos{\delta }_{J}cos\left({\alpha }_J-\mathsf{\Omega }\right)+sinIsin{\delta }_J+cosIcos{\delta }_{J}sin\left({\alpha }_J-\mathsf{\Omega }\right)\right]\right\},\hfill \end{array}$$

(114)

which is also defined for that value of the inclination. In such a case, using the true longitude (l) yields

$$\begin{array}{cc}\hfill \Delta T_{\mathrm{sid}}^{J_2}& =-\frac{3\mathsf{\pi }{J}_2{R}_{\mathrm{e}}^2}{4{\left(1-e^2\right)}^2\sqrt{\mathsf{\mu }a}}\left[\frac{\left(-2+3{cos}^2{\delta }_J\right)}{{\left(1+ecos\varpi \right)}^2}\left[2+e^2-2{\left(1-e^2\right)}^{3/2}+4ecos\varpi +e^{2}cos2\varpi \right]\right.\hfill \\ & + \frac{1}{2\left(1-e^2\right)}\left(\left(4+e^2\right)\left(1-3cos2{\delta }_J\right)-e\left[-1+6{cos}^2\left(l_0-{\alpha }_J\right)cos2{\delta }_J\right]\left[3ecos\left(2l_0-2\varpi \right)+6cos\left(l_0-\varpi \right)\right.\right.\hfill \\ & \left.\left.\left.+ 2e^2{cos}^3\left(l_0-\varpi \right)\right]-3cos\left(2l_0-2{\alpha }_J\right)\left\{\left(2+3e^2\right)cos2{\delta }_J+2{\left[1+ecos\left(l_0-\varpi \right)\right]}^3\right\}\right)\right].\hfill \end{array}$$

(115)

In the limit ($e\to 0$), Equation (115) agrees with Equation (114) calculated for $I=0$.

Figure 16 confirms the analytical results of Equations (86) and (113). Indeed, over three orbital revolutions, the projection of a generic $J_2$-perturbed orbit in the fundamental plane ($\left\{x,y\right\}$) crosses a fixed direction in the latter set by a certain value (${\varphi }_0$), always after a time interval equal to $T_{\mathrm{sid}}^{J_2}=T_{\mathrm{dra}}^{J_2}+\Delta T_{\mathrm{sid}\mathrm{II}}^{J_2}$ after each orbit. For the particular choice of the primary’s spin and the orbital parameters used in the picture, $T_{\mathrm{sid}}^{J_2}$ turns out to be shorter than $T_{\mathrm{K}}$, in agreement with Equations (86) and (113).

Furthermore, Figure 17 plots the final part of the time series of the cosine of the angle ($\varphi$) normalized to its initial value ($cos{\varphi }_0$) versus time (t) in units of $T_{\mathrm{K}}$ for a numerically integrated fictitious test particle with and without Equation (34) starting from the same generic initial position. It can be seen that it comes back to the same position in the fixed direction chosen in the reference plane, i.e., it is $cos\varphi /cos{\varphi }_0=+1$ again, just after $T_{\mathrm{sid}}^{J_2}=T_{\mathrm{dra}}^{J_2}+\Delta T_{\mathrm{sid}\mathrm{II}}^{J_2}$, differing from $T_{\mathrm{K}}$ by a positive amount, in agreement with Equations (86) and (113) for the particular choice of the generic values of the spin and the orbital parameters adopted in the numerical integrations.

6. Some Numerical Evaluations

The accuracy in measuring the orbital period of several transiting exoplanets [69] is nowadays at the

$${\sigma }_T\simeq {10}^{-7}-{10}^{-8}\mathrm{d}\simeq 9\times {10}^{-3}-{10}^{-4}\mathrm{s},$$

(116)

level (see https://exoplanet.eu/home/, accessed on 1 July 2024 and https://exoplanetarchive.ipac.caltech.edu/, accessed on 1 July 2024). Even better accuracies are available for other classes of objects; suffice it to say that for the period of WD1032 + 011 b, an inflated brown dwarf in an old eclipsing binary with a white dwarf with

$$\begin{array}{c}\hfill \nu =0.11,\end{array}$$

(117)

is known with an uncertainty as small as [70]

$${\sigma }_T\simeq 4.5\times {10}^{-10}\mathrm{d}=3.8\times {10}^{-5}\mathrm{s}.$$

(118)

The relevant orbital and physical parameters of WD1032 + 011 b are [70]

$$\begin{array}{cc}\hfill M_{\mathrm{p}}/M_{\odot }& =0.0665\pm 0.0061\hfill \\ \hfill \end{array}$$

(119)

$$\begin{array}{cc}\hfill M_{\mathrm{s}}/M_{\odot }& =0.4502\pm 0.0500\hfill \\ \hfill \end{array}$$

(120)

$$\begin{array}{cc}\hfill a/R_{\odot }& =0.6854\pm 0.0244,\hfill \end{array}$$

(121)

where p and s designate the planet and the star, respectively, and $M_{\odot }$ y $R_{\odot }$ are the mass and the radius of the Sun, respectively.

By assuming that the measured orbital period is the sidereal, a circular orbit is assumed for WD1032 + 011 b [70]. Equation (78), calculated with Equations (119)–(121), yields

$$\Delta T_{\mathrm{dra}}^{1\mathrm{pN}}=\Delta T_{\mathrm{sid}}^{1\mathrm{pN}}=0.07\pm 0.004\mathrm{s}.$$

(122)

Equation (118) shows that, in principle, the 1 pN gravitoelectric correction to the Keplerian orbital period given by Equation (122) falls within the measurability regime. On the other hand, any excessive optimism should be tempered, since the Keplerian term should be ultimately subtracted from the measured period in order to extract the 1 pN component. This implies that the values of the parameters entering the former should be known to a sufficiently high level of accuracy, which is not yet the case. Indeed, from the errors in Equations (119)–(121), one can calculate that the resulting uncertainty in the Keplerian period is as large as

$$\begin{array}{c}\hfill {\sigma }_{T_{\mathrm{K}}}\simeq 572\mathrm{s}.\end{array}$$

(123)

An evaluation of the other pK corrections is not possible, since they depend on $\mathit{J}$, $J_2$ and the mutual spin–orbit orientation about which no information was provided by Casewell et al. [70].

The importance of the issue of the mismodeling in the Keplerian period can be clearly understood in the case of double pulsar PSR J0737–3039 [71,72], characterized by

$$\nu =0.2497,$$

(124)

and whose (anomalisticA. Possenti, personal communication, May 2024.) orbital period is measured with an accuracy as good as [73]

$${\sigma }_T=5\times {10}^{-11}\mathrm{d}\simeq 4.32\times {10}^{-6}\mathrm{s}.$$

(125)

In principle, the exquisite accuracy of Equation (125) would allow for an accurate determination of the 1 pN gravitoelectric correction to the anomalistic period of PSR J0737–3039. Indeed, according to Equation (60), it nominally falls in the range of

$$0.27\mathrm{s}\lesssim \Delta T_{\mathrm{ano}}^{1\mathrm{pN}}\lesssim 0.40\mathrm{s},$$

(126)

depending on $f_0$. Unfortunately, the uncertainty in the Keplerian period, calculated by propagating the errors in the relevant parameters entering it [73], turns out to beIncidentally, the difference between the measured orbital period and the Keplerian one is not statistically significant, since it can be calculated to be as small as $1.9\mathrm{s}$.

$${\sigma }_{T_{\mathrm{K}}}\simeq 9\mathrm{s}.$$

(127)

Should it be possible to independently measure at least two of the three characteristic orbital periods for the same system, their common Keplerian component would be automatically canceled by determining their difference.

7. Summary and Conclusions

It was shown that post-Keplerian accelerations perturbing a two-body gravitationally bound system, like those arising from the oblateness of the central body to the Newtonian level and from the post-Newtonian gravitoelectromagnetic mass and spin-dependent components of its external gravitational potential, breaks the degeneracy between the otherwise coincident anomalistic, draconitic and sidereal orbital periods.

The resulting corrections to the Keplerian orbital period are generally different for the aforementioned characteristic timescales. The sidereal period still coincides with the draconitic one when the 1 pN gravitoelectric acceleration is taken into account, both being different from the anomalistic period. The 1 pN gravitomagnetic LT acceleration leaves the anomalistic period unaffected with respect to the Keplerian case, while it differently modifies the draconitic and sidereal periods. Finally, the oblateness of the central body alters all three orbital periods in different ways. In general, all the non-vanishing corrections to the Keplerian orbital period depend on the true anomaly in the given epoch.

The resulting analytical expressions are completely general, since they hold for arbitrary values of the orbital eccentricity and inclination. They are also valid for generic orientations of the primary’s symmetry axis in space.

For transiting brown dwarf WD1032 + 011 b ($\nu =0.11$), the predicted 1 pN gravitoelectric correction to the (sidereal) orbital period amounts to 0.07 s, while the current uncertainty in measuring it is as small as $\simeq {10}^{-5}$ sand the mismodeling in the Keplerian part is as large as 572 s. For double pulsar PSR J0737–3039 ($\nu =0.2497$), the 1 pN gravitoelectric correction to the anomalistic period, which is the measured one for this astrophysical system, is as large as a few tenths of a second; although the experimental accuracy in measuring the apsidal period is $\simeq {10}^{-6}\mathrm{s}$, the present-day uncertainty in the calculated value of the Keplerian component is still too large, amounting to about 9 s. If it were possible to independently measure at least two of the three characteristic orbital timescales for the same system, the difference between them would allow the Keplerian term to be canceled a priori.