On the 2PN Periastron Precession of the Double Pulsar PSR J0737–3039A/B

Abstract

One of the post-Keplerian (PK) parameters determined in timing analyses of several binary pulsars is the fractional periastron advance per orbit $k^{\mathrm{PK}}$. Along with other PK parameters, it is used in testing general relativity once it is translated into the periastron precession ${\dot{\omega }}^{\mathrm{PK}}$. It was recently remarked that the periastron $\omega$ of PSR J0737–3039A/B may be used to measure/constrain the moment of inertia of A through the extraction of the general relativistic Lense–Thirring precession ${\dot{\omega }}^{\mathrm{LT},\mathrm{A}}\simeq -0.{00060}^{\circ }{\mathrm{yr}}^{-1}$ from the experimentally determined periastron rate ${\dot{\omega }}_{\mathrm{obs}}$ provided that the other post-Newtonian (PN) contributions to ${\dot{\omega }}_{\exp }$ can be accurately modeled. Among them, the 2PN seems to be of the same order of magnitude of ${\dot{\omega }}^{\mathrm{LT},\mathrm{A}}$. An analytical expression of the total 2PN periastron precession ${\dot{\omega }}^{2\mathrm{PN}}$ in terms of the osculating Keplerian orbital elements, valid not only for binary pulsars, is provided, thereby elucidating the subtleties implied in correctly calculating it from $k^{1\mathrm{PN}}+k^{2\mathrm{PN}}$ and correcting some past errors by the present author. The formula for ${\dot{\omega }}^{2\mathrm{PN}}$ is demonstrated to be equivalent to that obtainable from $k^{1\mathrm{PN}}+k^{2\mathrm{PN}}$ by Damour and Schäfer expressed in the Damour–Deruelle (DD) parameterization. ${\dot{\omega }}^{2\mathrm{PN}}$ actually depends on the initial orbital phase, hidden in the DD picture, so that $-0.{00080}^{\circ }{\mathrm{yr}}^{-1}\le {\dot{\omega }}^{2\mathrm{PN}}\le -0.{00045}^{\circ }{\mathrm{yr}}^{-1}$. A recently released prediction of ${\dot{\omega }}^{2\mathrm{PN}}$ for PSR J0737–3039A/B is discussed.

1. Introduction

Recently, Hu et al. [1] performed a detailed analysis of the perspectives of measuring, or effectively constraining, the moment of inertia (MOI) ${\mathcal{I}}_{\mathrm{A}}$ of the pulsar PSR J0737–3039A [2,3] by the end of the present decade by exploiting the general relativistic spin-orbit Lense–Thirring periastron precession ${\dot{\omega }}^{\mathrm{LT},\mathrm{A}}$[4] induced by its spin angular momentum ${\mathit{S}}^{\mathrm{A}}$. Among the competing dynamical effects acting as potential sources of systematic uncertainty, Hu et al. [1] included also the periastron precession ${\dot{\omega }}^{2\mathrm{PN}}$ to the second post-Newtonian (2PN) order which, along with the much larger1 1PN precession

$${\dot{\omega }}^{1\mathrm{PN}}=\frac{3n^{\mathrm{K}}\mu }{c^{2}a\left(1-e^2\right)},$$

(1)

depends only on the masses $M_{\mathrm{A}},M_{\mathrm{B}}$ of both the neutron stars which the Double Pulsar PSR J0737–3039A/B is made of. In Equation (1), c is the speed of light in vacuum, $\mu \doteq GM$ is the gravitational parameter of the Double Pulsar given by the product of the Newtonian constant of gravitation G times the sum of the masses $M\doteq M_{\mathrm{A}}+M_{\mathrm{B}}$, a and e are the osculating numerical values of the semimajor axis and eccentricity, respectively, at the same arbitrary moment of time $t_0$ [5], while

$$n^{\mathrm{K}}\doteq \sqrt{\frac{\mu }{a^3}}$$

(2)

is the osculating Keplerian mean motion. In particular, Table 1 in Hu et al. [1] reported

$${\dot{\omega }}^{2\mathrm{PN}}=0.{000439}^{\circ }{\mathrm{yr}}^{-1}=1.58{}^{″}{\mathrm{yr}}^{-1},$$

(3)

for the 2PN periastron precession which would, thus, be prograde. In Equation (3), ${}^{″}$ stands for arcseconds. Equation (3) is to be compared with the retrograde Lense–Thirring periastron rate which, if calculated with the latest determination of ${\mathcal{I}}_{\mathrm{A}}$ by Silva et al. [6], would be of comparable magnitude 2

$${\dot{\omega }}^{\mathrm{LT},\mathrm{A}}\simeq -0.{0006}^{\circ }{\mathrm{yr}}^{-1}=-2.16{}^{″}{\mathrm{yr}}^{-1}.$$

(4)

An accurate prediction of the 2PN periastron precession, or of the experimental quantity related to it, which is actually determined in real data analyses, is of the utmost importance since, according to Equation (3), it may cancel Equation (4) to a large extent. To this aim, it is important to stress that, although seemingly unnoticed so far in the literature, a certain amount of uncertainty should be deemed as still lingering on that matter because, perhaps, of how $\dot{\omega }$ is routinely expressed in most of the papers devoted to binary pulsars. Indeed, as it will be shown here, the way usually adopted in the literature to write the total 2PN periastron precession hides its dependence on the initial conditions which, indeed, is buried in the parameterization used. Such a distinctive feature does not occur at the 1PN level whose averaged orbital precessions such as Equation (1) and the Lense–Thirring one ${\dot{\omega }}^{\mathrm{LT}}$ are independent of the orbital phase at a reference epoch. Thus, while the predictions of the 1PN precessions are valid for any starting time, it is not so for the 2PN ones, despite their—purely formal—independence of the initial conditions in certain parameterizations. Moreover, there is some confusion about the periastron precession and how to correctly calculate it from the fractional periastron advance per orbit. Finally, in numerically calculating ${\dot{\omega }}^{2\mathrm{PN}}$, the fact that also the formal 1PN term contributes it in a subtle way is often overlooked yielding incomplete results.

The paper is organized as follows. In Section 2, the total 2PN3 periastron rate is calculated (see Equation (18)) by using the osculating Keplerian orbital elements from existing expressions in the literature for the fractional PN periastron shift per orbit $k^{1\mathrm{PN}}+k^{2\mathrm{PN}}$. In particular, Equation (21) by Iorio [11] is used as starting point in Section 2.1, where an error by Iorio [11] in obtaining the true total 2PN periastron rate is disclosed and corrected. In Section 2.2, Equation (18) is obtained starting from Equation (5.18) by Damour and Schäfer [4], expressed in the DD parameterization, after a proper conversion from the latter to the osculating Keplerian orbital elements. In Section 3, Equation (18) is confirmed by numerically integrating the PN equations of motion up to the 2PN order for a fictitious binary system. The results of Section 2 are applied to other astrophysical and astronomical systems of interest in Section 4. The case of PSR J0737–3039A/B is dealt with in Section 4.1, where Equation (3) is discussed as well. Section 4.2 treats Mercury, the spacecraft Juno orbiting Jupiter, the Earth’s artificial satellites LAGEOS II, and the S-star S4711 around Sgr A${}^{*}$. Section 5 summarizes the findings obtained and offers concluding remarks.

2. How to Correctly Calculate the 2PN Periastron Precession from the Fractional 1PN + 2PN Periastron Shift per Orbit Using the Osculating Keplerian Orbital Elements

In pulsar timing analyses, one of the so-called post-Keplerian (PK) parameters that are determined is the fractional periastron shift per orbit $k^{\mathrm{PK}}$ [4,12,13] defined as

$$k^{\mathrm{PK}}\doteq \frac{〈\Delta {\omega }^{\mathrm{PK}}〉}{2\pi }.$$

(5)

In Equation (5), $\omega$ is the argument of periastron, $\Delta {\omega }^{\mathrm{PK}}$ is the time-dependent shift of periastron induced by some PK dynamical extra-acceleration with respect to the Newtonian inverse-square one, and the angular brackets $〈\cdots 〉$ denote the average over the orbital period $P^{\mathrm{PK}}$ which, in presence of PK accelerations, has to be meant as the anomalistic period $P_{\mathrm{ano}}^{\mathrm{PK}}$, i.e., the time span between two successive crossings of the (moving) periastron position. Nonetheless, it is common practice to deal with the averaged4 periastron precession ${\dot{\omega }}^{\mathrm{PK}}$ which is connected to $k^{\mathrm{PK}}$ through

$$k^{\mathrm{PK}}=\frac{{\dot{\omega }}^{\mathrm{PK}}}{n^{\mathrm{PK}}},$$

(6)

where

$$n^{\mathrm{PK}}\doteq \frac{2\pi }{P^{\mathrm{PK}}}$$

(7)

is the PK mean motion. In general, Equation (7) differs from Equation (2).

2.1. Starting from the Formula by Iorio in Osculating Keplerian Orbital Elements

As far as the 2PN periastron advance is concerned, Equation (21) in Iorio [11] correctly calculated $k^{2\mathrm{PN}}$, up to the scaling factor $n^{\mathrm{K}}$entering Equation (21) of Iorio [11], with the Gauss perturbing equations in terms of the osculating Keplerian orbital elements (see Equation (9)) by showing that his expression agrees with those obtained by Kopeikin and Potapov [14] with the same perturbative technique but a different calculational strategy, and by Damour and Schäfer [4] who, instead, used the Hamilton-Jacobi method and the Damour–Deruelle (DD) parameterization [15] which is nowadays routinely used in standard pulsar timing analyses [12,13]. Iorio [11], after having scaled $k^{2\mathrm{PN}}$ by Equation (2), erroneously claimed that the resulting expression for $n^{\mathrm{K}}{k}^{2\mathrm{PN}}$ is the total 2PN pericentre precession, which is not the case, as it will be shown below. Here, the explicit expressions of $k^{1\mathrm{PN}},k^{2\mathrm{PN}}$ in terms of the osculating Keplerian orbital elements are reported. They are

$$\begin{array}{cc}\hfill k^{1\mathrm{PN}}& =\frac{3\mu }{c^{2}a\left(1-e^2\right)},\hfill \\ \hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill k^{2\mathrm{PN}}& =\frac{3{\mu }^2\left[2-4\nu +e^2\left(1+10\nu \right)+16e\left(-2+\nu \right)cos{f}_0\right]}{4c^4{a}^2{\left(1-e^2\right)}^2},\hfill \end{array}$$

(9)

where $f_0$ is the osculating numerical value of the true anomaly f at some arbitrary moment of time $t_0$, and

$$\nu \doteq \frac{M_{\mathrm{A}}{M}_{\mathrm{B}}}{M^2}.$$

(10)

To correctly calculate the total 2PN pericentre precession ${\dot{\omega }}^{2\mathrm{PN}}$, some characteristic time interval playing the role of “orbital period” has to be worked out to the 1PN order. In the present case, the anomalistic period, i.e., the time interval between two successive crossings of the (moving) pericentre position, fulfils such a requirement. To the 1PN order, it can be written as

$$P_{\mathrm{ano}}^{1\mathrm{PN}}=P^{\mathrm{K}}+\Delta P_{\mathrm{ano}}^{1\mathrm{PN}},$$

(11)

where the osculating Keplerian period is

$$P^{\mathrm{K}}\doteq \frac{2\pi }{n^{\mathrm{K}}}=2\pi \sqrt{\frac{a^3}{\mu }},$$

(12)

and the 1PN correction, calculated according to the strategy followed by Iorio [16], turns out to be

$$\Delta P_{\mathrm{ano}}^{1\mathrm{PN}}=\frac{\pi \sqrt{a\mu }}{2c^2{\left(1-e^2\right)}^2}{\mathcal{T}}_{\mathrm{ano}}^{1\mathrm{PN}},$$

(13)

with

$$\begin{array}{cc}\hfill {\mathcal{T}}_{\mathrm{ano}}^{1\mathrm{PN}}& \doteq 36+e^2\left(42-38\nu \right)+2e^4\left(6-7\nu \right)-8\nu +\hfill \\ \hfill \\ \hfill & +3e\left\{\left[28+3e^2\left(4-5\nu \right)-12\nu \right]cos{f}_0-e\left(-10+8\nu +e\nu cos{f}_0\right)cos2f_0\right\}.\hfill \end{array}$$

(14)

In the point particle limit corresponding to $\nu \to 0$, Equation (14) reduces to Equation (72) of Iorio [16]. The 1PN mean motion is, thus,

$$n^{1\mathrm{PN}}\doteq \frac{2\pi }{P_{\mathrm{ano}}^{1\mathrm{PN}}}=\frac{n^{\mathrm{K}}}{1+\frac{\mu }{4c^{2}a{\left(1-e^2\right)}^2}{\mathcal{T}}_{\mathrm{ano}}^{1\mathrm{PN}}}.$$

(15)

Note that Equation (21) of Iorio [11], i.e., the product of Equation (2) times Equation (9)

$$n^{\mathrm{K}}{k}^{2\mathrm{PN}}=\frac{3{\mu }^{5/2}\left[2-4\nu +e^2\left(1+10\nu \right)+16e\left(-2+\nu \right)cos{f}_0\right]}{4c^4{a}^{7/2}{\left(1-e^2\right)}^2},$$

(16)

has just formally the dimensions of a pericentre precession of the order of $\mathcal{O}\left(c^{-4}\right)$, but, contrary to what mistakenly claimed by Iorio [11], it is not the total 2PN pericentre rate ${\dot{\omega }}^{2\mathrm{PN}}$. Indeed, the correct analytical expression for it can only be obtained by retaining the term of the order of $\mathcal{O}\left(c^{-4}\right)$ in the expansion in powers of $c^{-1}$ of the product of Equation (15) times the sum of Equations (8) and (9). If, on the one hand, replacing $n^{\mathrm{K}}$ with Equation (15) does not affect Equation (16) in the power expansion to the 2PN order, on the other hand, it does matter when it is Equation (8) that is multiplied by Equation (15) and power-expanded to the order of $\mathcal{O}\left(c^{-4}\right)$. Indeed, from Equations (8) and (15) one has

$$\begin{array}{cc}\hfill {\left.n^{1\mathrm{PN}}{k}^{1\mathrm{PN}}\right|}_{2\mathrm{PN}}& =\frac{3{\mu }^{5/2}}{4c^4{a}^{7/2}{\left(1-e^2\right)}^3}\left(-36+8\nu +2e^4\left(-6+7\nu \right)+e^2\left(-42+38\nu \right)+\right.\hfill \\ \hfill \\ \hfill & +3e\left\{\left[4\left(-7+3\nu \right)+3e^2\left(-4+5\nu \right)\right]cos{f}_0+\right.\hfill \\ \hfill \\ \hfill & +\left.\left.e\left(-10+8\nu +e\nu cos{f}_0\right)cos2f_0\right\}\right),\hfill \end{array}$$

(17)

which, added to Equation (16), yields

$$\begin{array}{cc}\hfill {\dot{\omega }}^{2\mathrm{PN}}& =\frac{3{\mu }^{5/2}}{8c^4{a}^{7/2}{\left(1-e^2\right)}^3}\left\{-68+8\nu +e^4\left(-26+8\nu \right)+2e^2\left(-43+52\nu \right)+\right.\hfill \\ \hfill \\ \hfill & +e\left[8\left(-29+13\nu \right)+e^2\left(-8+61\nu \right)\right]cos{f}_0+\hfill \\ \hfill \\ \hfill & +\left.3e^2\left[4\left(-5+4\nu \right)cos2f_0+e\nu cos3f_0\right]\right\}.\hfill \end{array}$$

(18)

This is the right analytical expression for the full 2PN pericentre precession expressed in terms of the osculating Keplerian orbital elements. In the limit $e\to 0$, Equation (18) is independent of $f_0$.

Recapitulating, on the one hand, Equation (21) in Iorio [11] correctly worked out the 2PN fractional pericentre shift per orbit $k^{2\mathrm{PN}}$ up to $n^{\mathrm{K}}$ as scaling factor. On the other hand, Iorio [11] mistakenly claimed that the resulting expression for $n^{\mathrm{K}}{k}^{2\mathrm{PN}}$ was the total 2PN pericentre precession, missing a further contribution from the power expansion to the 2PN order of the product $n^{1\mathrm{PN}}{k}^{1\mathrm{PN}}$.

2.2. Starting from the Formula by Damour and Schäfer in the Damour–Deruelle Parameterization

Equation (18) is in agreement also with the expression for the total 2PN pericentre precession, written in terms of the osculating Keplerian orbital elements, which can be extracted from Equation (5.18) in [4]

$$\begin{array}{cc}\hfill k^{1\mathrm{PN}}+k^{2\mathrm{PN}}& =\frac{3{\left(\mu n_{\mathrm{DD}}\right)}^{2/3}}{c^2\left(1-e_{\mathrm{T}}^2\right)}\left[1+\frac{{\left(\mu n_{\mathrm{DD}}\right)}^{2/3}}{c^2\left(1-e_{\mathrm{T}}^2\right)}\left(\frac{39}{4}{x}_{\mathrm{A}}^2+\frac{27}{4}{x}_{\mathrm{B}}^2+15x_{\mathrm{A}}{x}_{\mathrm{B}}\right)-\right.\hfill \\ \hfill \\ \hfill & \left.-\frac{{\left(\mu n_{\mathrm{DD}}\right)}^{2/3}}{c^2}\left(\frac{13}{4}{x}_{\mathrm{A}}^2+\frac{1}{4}{x}_{\mathrm{B}}^2+\frac{13}{3}{x}_{\mathrm{A}}{x}_{\mathrm{B}}\right)\right].\hfill \end{array}$$

(19)

In Equation (19),

$$\begin{array}{cc}\hfill x_{\mathrm{A}}& \doteq \frac{M_{\mathrm{A}}}{M},\hfill \\ \hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill x_{\mathrm{B}}& \doteq \frac{M_{\mathrm{B}}}{M}=1-x_{\mathrm{A}},\hfill \end{array}$$

(21)

while $e_{\mathrm{T}}$ and $n_{\mathrm{DD}}$ are members of the Damour–Deruelle (DD) formalism [15] which, in the limit $c\to \infty$, reduce to the Keplerian eccentricity e and mean motion $n^{\mathrm{K}}$, as it will be shown below.

The “proper time” eccentricity $e_{\mathrm{T}}$ reads p. 272 in [12]

$$e_{\mathrm{T}}=e_t\left(1+\delta \right)+e_{\theta }-e_r,$$

(22)

where Equation (3.8b) in [15]

$$e_t=\frac{e_R}{1+\frac{\mu }{c^2{a}_R}\left(4-\frac{3}{2}\nu \right)},$$

(23)

Equation (4.13) in [15]

$$e_{\theta }=e_R\left(1+\frac{\mu }{2c^2{a}_R}\right),$$

(24)

Equation (6.3b) in [15]

$$e_r=e_R\left[1-\frac{\mu }{2c^2{a}_R}\left(x_{\mathrm{A}}^2-\nu \right)\right],$$

(25)

and Equation (20) in [12]

$$\delta =\frac{\mu }{c^2{a}_R}\left(x_{\mathrm{A}}{x}_{\mathrm{B}}+2x_{\mathrm{B}}^2\right).$$

(26)

In Equations (23)–(26), $a_R$ is another member of the DD parameterization. According to Equations (23)–(26), Equation (22) can be expressed in terms of only $a_R,e_R$ as

$$\frac{e_{\mathrm{T}}}{e_R}=\frac{1+\frac{\mu }{2c^2{a}_R}\left[4+3\left(x_{\mathrm{A}}-2\right)x_{\mathrm{A}}\right]+\frac{{\mu }^2}{4c^4{a}_R^2}\left(8-3\nu \right)x_{\mathrm{A}}^2}{1+\frac{\mu }{2c^2{a}_R}\left(8-3\nu \right)}.$$

(27)

The DD mean motion is Equation (3.7) in [15]

$$n_{\mathrm{DD}}\doteq \sqrt{\frac{\mu }{a_R^3}}\left[1+\frac{\mu }{2c^2{a}_R}\left(-9+\nu \right)\right].$$

(28)

Equations (27) and (28) are both functions of $a_R,e_R$ which, in turn, can be expressed in terms of the osculating Keplerian semimajor axis a and eccentricity e by means of Equations (28) and (29) in [5]

$$\begin{array}{cc}\hfill a_R& =a-d{a}_0-\frac{\mu }{c^2{\left(1-e^2\right)}^2}\left[-3+\nu +e^2\left(-13+e^2+7\nu +2e^2\nu \right)\right],\hfill \\ \hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill e_R& =e-d{e}_0-\frac{e\mu }{2c^{2}a\left(1-e^2\right)}\left[-17+6\nu +e^2\left(2+4\nu \right)\right],\hfill \end{array}$$

(30)

with Equation (14) in [5]

$$\begin{array}{cc}\hfill d{a}_0& =\frac{e\mu }{4c^2{\left(1-e^2\right)}^2}\left\{\left[8\left(-7+3\nu \right)+e^2\left(-24+31\nu \right)\right]cos{f}_0+\right.\hfill \\ \hfill \\ \hfill & \left.+e\left[4\left(-5+4\nu \right)cos2f_0+e\nu cos3f_0\right]\right\},\hfill \end{array}$$

(31)

and Equation (16) in [5]

$$\begin{array}{cc}\hfill d{e}_0& =\frac{\mu }{8c^{2}a\left(1-e^2\right)}\left\{\left[8\left(-3+\nu \right)+e^2\left(-56+47\nu \right)\right]cos{f}_0+\right.\hfill \\ \hfill \\ \hfill & \left.+e\left[4\left(-5+4\nu \right)cos2f_0+e\nu cos3f_0\right]\right\}.\hfill \end{array}$$

(32)

Note that Equations (29)–(32) are written for general relativity; their general expressions for a given class of alternative theories of gravitation can be found in Klioner and Kopeikin [5]. The final expressions for $a_R,e_R$ are

$$\begin{array}{cc}\hfill \frac{a_R}{a}& =1-\frac{\mu }{c^{2}a{\left(1-e^2\right)}^2}\left[-3+\nu +e^4\left(1+2\nu \right)+e^2\left(-13+7\nu \right)\right]+\hfill \\ \hfill \\ \hfill & +e\frac{\mu }{4c^{2}a{\left(1-e^2\right)}^2}\left\{\left[56+e^2\left(24-31\nu \right)-24\nu \right]cos{f}_0+\right.\hfill \\ \hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill & \left.+e\left[4\left(5-4\nu \right)cos2f_0-e\nu cos3f_0\right]\right\},\hfill \\ \hfill \\ \hfill \frac{e_R}{e}& =1-\frac{\mu }{2c^{2}a\left(1-e^2\right)}\left[-17+6\nu +e^2\left(2+4\nu \right)\right]+\hfill \\ \hfill \\ \hfill & -\frac{\mu }{8c^{2}ae\left(1-e^2\right)}\left\{\left[8\left(-3+\nu \right)+e^2\left(-56+47\nu \right)\right]cos{f}_0+\right.\hfill \\ \hfill & +\left.e\left[4\left(-5+4\nu \right)cos2f_0+e\nu cos3f_0\right]\right\}.\hfill \end{array}$$

(34)

By using Equations (33) and (34), Equations (27) and (28) can be finally expressed, to the order of $\mathcal{O}\left(c^{-2}\right)$, as

$$\begin{array}{cc}\hfill \frac{8c^{2}a\left(e-e_{\mathrm{T}}\right)\left(1-e^2\right)}{\mu }& =\left[8\left(-3+\nu \right)+e^2\left(-56+47\nu \right)\right]cos{f}_0+\hfill \\ \hfill \\ \hfill & +e\left(4\left\{-13+3\nu -3\left(-2+x_{\mathrm{A}}\right)x_{\mathrm{A}}+\right.\right.\hfill \\ \hfill \\ \hfill & \left.+e^2\left[-2+7\nu +3\left(-2+x_{\mathrm{A}}\right)x_{\mathrm{A}}\right]\right\}+\hfill \\ \hfill & \left.+4\left(-5+4\nu \right)cos2f_0+e\nu cos3f_0\right),\hfill \end{array}$$

(35)

$$\begin{array}{c}\hfill \\ \hfill \left(\frac{n_{\mathrm{DD}}}{n^{\mathrm{K}}}-1\right)\frac{8c^{2}a{\left(1-e^2\right)}^2}{\mu }& =8\left(-9+2\nu \right)+4e^4\left(-6+7\nu \right)+e^2\left(-84+76\nu \right)+\hfill \\ \hfill \\ \hfill & +3e\left\{\left[8\left(-7+3\nu \right)+e^2\left(-24+31\nu \right)\right]cos{f}_0+\right.\hfill \\ \hfill & \left.+e\left[4\left(-5+4\nu \right)cos2f_0+e\nu cos3f_0\right]\right\}.\hfill \end{array}$$

(36)

A power expansion to the order of $\mathcal{O}\left(c^{-4}\right)$ of the product of Equation (28) by Equation (19), calculated with Equations (35) and (36), yields just Equation (18).

About the fractional periastron shift per orbit, Iorio [11] demonstrated the equivalence of his Equation (21), up to the scaling factor $n^{\mathrm{K}}$, with the total explicit components of the order of $\mathcal{O}\left(c^{-4}\right)$ of Equation (19) and of an analogous formula by Kopeikin and Potapov [14], once the proper translation of the latter ones into the osculating Keplerian orbital elements was appropriately carried out.

3. Numerically Integrating the 1PN + 2PN Equations of Motion

The correctness of Equation (18), and also its general applicability to whatsoever binary system for which the PN approximation is deemed applicable, can be numerically demonstrated in the following way. For the sake of clarity, a fictitious two-body system made of, say, two supermassive black holes with $M_{\mathrm{A}}=1\times {10}^6{M}_{\odot },M_{\mathrm{B}}=2\times {10}^6{M}_{\odot }$ orbiting along a highly eccentric ($e=0.75$) orbit in $0.05\mathrm{yr}$ is considered. For a given set of initial conditions, parameterized in terms of the5 Keplerian orbital elements, the equations of motion, in rectangular Cartesian harmonic coordinates, including the PN accelerations [17] (see, e.g., Equation (4.4.28), p. 154); [18] (Equation (A2.6), p. 166); [19] (Equation (10.3.7), p. 381)

$$\begin{array}{cc}\hfill {\mathbf{A}}^{1\mathrm{PN}}& =\frac{\mu }{c^2{r}^2}\left\{\left[\left(4+2\nu \right)\frac{\mu }{r}+\frac{3}{2}\nu {\mathrm{v}}_r^2-\left(1+3\nu \right){\mathrm{v}}^2\right]\widehat{\mathit{r}}+\left(4-2\nu \right){\mathrm{v}}_r\mathbf{v}\right\},\hfill \end{array}$$

(37)

and [17] (see, e.g., Equation (4.4.29), p. 154); [20] (Equation (2.2d), p. 825); [21] (Equation (B11), p. 10)

$$\begin{array}{cc}\hfill {\mathbf{A}}^{2\mathrm{PN}}& =\frac{\mu }{c^4{r}^2}\left\{\left[\nu \left(-3+4\nu \right){\mathrm{v}}^4+\frac{15}{8}\nu \left(-1+3\nu \right){\mathrm{v}}_r^4+\nu \left(\frac{9}{2}-6\nu \right){\mathrm{v}}^2{\mathrm{v}}_r^2+\nu \left(\frac{13}{2}-2\nu \right)\frac{\mu }{r}{\mathrm{v}}^2+\right.\right.\hfill \\ \hfill \\ \hfill & \left.+\left(2+25\nu +2{\nu }^2\right)\frac{\mu }{r}{\mathrm{v}}_r^2-\left(9+\frac{87}{4}\nu \right)\frac{{\mu }^2}{r^2}\right]\widehat{\mathit{r}}+\left[\nu \left(\frac{15}{2}+2\nu \right){\mathrm{v}}^2-\nu \left(\frac{9}{2}+3\nu \right){\mathrm{v}}_r^2-\right.\hfill \\ \hfill \\ \hfill & \left.\left.-\left(2+\frac{41}{2}\nu +4{\nu }^2\right)\frac{\mu }{r}\right]{\mathrm{v}}_r\mathbf{v}\right\}\hfill \end{array}$$

(38)

in addition to the Newtonian monopole

$${\mathit{A}}^{\mathrm{N}}=-\frac{\mu }{r^2}\widehat{\mathit{r}},$$

(39)

are numerically integrated over 1 yr with and without Equations (37) and (38) in each run, and time series of $\omega \left(t\right)$ are correspondingly calculated. In Equations (37)–(39), r is the relative distance between A and B, $\widehat{\mathit{r}}$ is the versor of the relative position vector, $\mathrm{v}$ is the velocity vector of the relative motion, and ${\mathrm{v}}_r\doteq \mathrm{v}·\widehat{\mathit{r}}$ is the radial velocity. Then, the difference between the Newtonian and the PN time series for $\omega \left(t\right)$ is taken to extract the time-dependent PN shift $\Delta {\omega }^{\mathrm{PN}}\left(t\right)$. By construction, it includes both the full 1PN and 2PN contributions along with other terms of higher order due to the interplay between the 1PN and 2PN accelerations, not of interest here. To single out just the total 2PN effect (up to other PN contributions of higher order) $\Delta {\omega }^{2\mathrm{PN}}\left(t\right)$, the 1PN linear trend, analytically calculated by multiplying $k^{1\mathrm{PN}}$ of Equation (8) times $n^{\mathrm{K}}t$, is subtracted from the time series $\Delta {\omega }^{\mathrm{PN}}\left(t\right)$. The same procedure is repeated by varying the true anomaly at epoch $f_0$ leading to a change of the initial conditions. The resulting time-dependent signatures for $\Delta {\omega }^{2\mathrm{PN}}\left(t\right)$ are displayed in the upper panel of Figure 1. A linear fit to each of them is performed, and the resulting straight lines are superimposed. Their slopes, in ${}^{\circ }{\mathrm{yr}}^{-1}$, can be compared with the lower panel of Figure 1 displaying the plot of Equation (18) as a function of $f_0$; the agreement is neat.

4. Application to Double Pulsar and Other Systems

4.1. The Case of PSR J0737–3039A/B

In the case of PSR J0737–3039A/B, Equation (18) yields

$$-0.{00080}^{\circ }{\mathrm{yr}}^{-1}\le {\dot{\omega }}^{2\mathrm{PN}}\le -0.{00045}^{\circ }{\mathrm{yr}}^{-1},$$

(40)

as shown by Figure 2, or, equivalently,

$$-2.8{}^{″}{\mathrm{yr}}^{-1}\le {\dot{\omega }}^{2\mathrm{PN}}\le -1.6{}^{″}{\mathrm{yr}}^{-1}.$$

(41)

Equations (40) and (41) correct the wrong range for ${\dot{\omega }}^{2\mathrm{PN}}$ which one would obtain for the Double Pulsar by summing the values in Equation (20) and Equation (21) of6 Iorio [11]. From Figure 2, the 2PN periastron precession of PSR J0737–3039A/B is always retrograde and does not vanish for any value of $f_0$.

About Equation (3) quoted by Hu et al. [1] (Table 1), it was obtained as follows. By taking the product of Equation (28), or, in this case, also of Equation (2), times the sum of only the second and the third term in7 Equation (19) and expanding the resulting expression to the order of $\mathcal{O}\left(c^{-4}\right)$, one gets the following quantity which is dimensionally a 2PN precession

$${\dot{\psi }}^{2\mathrm{PN}}=\frac{n^{\mathrm{K}}{\mu }^2\left\{78-28\nu +e^2\left[3+2\left(23-5x_{\mathrm{A}}\right)x_{\mathrm{A}}\right]\right\}}{4c^4{a}^2{\left(1-e^2\right)}^2}=0.{000439}^{\circ }{\mathrm{yr}}^{-1},$$

(42)

in agreement with Equation (3). This demonstrates that the prediction for the 2PN periastron precession by Hu et al. [1] (Table 1) is, in fact, incomplete since it neglected the contribution to ${\dot{\omega }}^{2\mathrm{PN}}$ of the product of Equation (28) times the first term in Equation (19), which is formally of the order of $\mathcal{O}\left(c^{-2}\right)$.

4.2. Other Astronomical Systems

For the Sun and Mercury, Equation (18) yields

$$-18\mathsf{\mu }\mathrm{as}{\mathrm{cty}}^{-1}\le {\dot{\omega }}^{2\mathrm{PN}}\le -4\mathsf{\mu }\mathrm{as}{\mathrm{cty}}^{-1},$$

(43)

while for the spacecraft Juno currently orbiting Jupiter, it is

$$-4\mathsf{\mu }\mathrm{as}{\mathrm{yr}}^{-1}\le {\dot{\omega }}^{2\mathrm{PN}}\le 0\mathsf{\mu }\mathrm{as}{\mathrm{yr}}^{-1}.$$

(44)

The 2PN perigee precession of the Earth’s artificial satellite LAGEOS II is as little as

$$-0.0108\mathsf{\mu }\mathrm{as}{\mathrm{yr}}^{-1}\le {\dot{\omega }}^{2\mathrm{PN}}\le -0.0100\mathsf{\mu }\mathrm{as}{\mathrm{yr}}^{-1}.$$

(45)

In Equations (43)–(45), $\mathsf{\mu }\mathrm{as}$ stands for microarcseconds.

Larger values occur, e.g., for the recently discovered S-star S4711 [22] orbiting the supermassive black hole in the Galactic Center at Sgr A${}^{*}$; it revolves around its primary in $7.6\mathrm{yr}$ along an orbit with an eccentricity as large as $e=0.768$. Its 2PN periastron precession range turns out to be

$$-1.4{}^{″}{\mathrm{yr}}^{-1}\le {\dot{\omega }}^{2\mathrm{PN}}\le 0.074{}^{″}{\mathrm{yr}}^{-1}.$$

(46)

5. Summary and Conclusions

Calculating correctly the total 2PN periastron precession ${\dot{\omega }}^{2\mathrm{PN}}$ from the fractional periastron advance per orbit$k^{\mathrm{PN}}=k^{1\mathrm{PN}}+k^{2\mathrm{PN}}$ requires to multiply the latter one by the 1PN mean motion $n^{1\mathrm{PN}}$ instead of the osculating Keplerian one $n^{\mathrm{K}}$, as incorrectly done by Iorio [11], and to expand the resulting expression to the order of $\mathcal{O}\left(c^{-4}\right)$. Adopting the osculating Keplerian orbital elements allows to obtain Equation (18) for ${\dot{\omega }}^{2\mathrm{PN}}$. It has a general validity, being straightforwardly applicable to whatsoever two-body system whose data are not analyzed within the DD framework, and clearly shows that the total 2PN periastron precession does depend on the initial conditions, as confirmed also by the numerical integration of the 1PN + 2PN equations of motion for a fictitious binary displayed in Figure 1. Also, the formula for $k^{1\mathrm{PN}}+k^{2\mathrm{PN}}$ by Damour and Schäfer [4], written in terms of the DD parameters, yields Equation (18) if properly multiplied by the DD version of the 1PN mean motion and after appropriate conversion from the DD parameters to the osculating Keplerian ones.

For PSR J0737–3039A/B, ${\dot{\omega }}^{2\mathrm{PN}}$ is retrograde for any value of the initial orbital phase, as shown by the plot of Equation (18) in Figure 2. Since it turns out that $-{0.00080}^{\circ }{\mathrm{yr}}^{-1}\le {\dot{\omega }}^{2\mathrm{PN}}\le -{0.00045}^{\circ }{\mathrm{yr}}^{-1}$, it adds up to the spin-orbit Lense–Thirring precession ${\dot{\omega }}^{\mathrm{LT},\mathrm{A}}\simeq -{0.0006}^{\circ }{\mathrm{yr}}^{-1}$.

The value ${\dot{\omega }}^{2\mathrm{PN}}={0.000439}^{\circ }{\mathrm{yr}}^{-1}$ by Hu et al. [1] comes from having neglected to multiply $k^{1\mathrm{PN}}$ by the 1PN mean motion and to expand the resulting product to the order of $\mathcal{O}\left(c^{-4}\right)$.

For some astronomical systems in the Solar System of potential interest, the 2PN pericenter precession is negligible, while for the S-star S4711 orbiting Sgr A${}^{*}$ it amounts to $-1.4{}^{″}{\mathrm{yr}}^{-1}\le {\dot{\omega }}^{2\mathrm{PN}}\le 0.074{}^{″}{\mathrm{yr}}^{-1}$.