On the 2PN periastron precession of the Double Pulsar PSR J0737-3039A/B

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,\,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_\mathrm{exp}$ can be accurately modeled. Among them, the 2PN one seems to be of the same order of magnitude of $\dot\omega^\mathrm{LT,\,A}$. An analytical expression of the total 2PN periastron precession $\dot\omega^\mathrm{2PN}$ in terms of the osculating Keplerian orbital elements, valid not only for binary pulsars, is provided elucidating the subtleties implied in correctly calculating it from $k^\mathrm{1PN}+k^\mathrm{2PN}$ and correcting some past errors by the present author. The formula for $\dot\omega^\mathrm{2PN}$ is demonstrated to be equivalent to that obtainable from $k^\mathrm{1PN}+k^\mathrm{2PN}$ by Damour and Sch\"afer expressed in the Damour-Deruelle (DD) parameterization. $\dot\omega^\mathrm{2PN}$ actually depends on the initial orbital phase, hidden in the DD picture, so that $-0.00080^\circ\,\mathrm{yr}^{-1} \leq\dot\omega^\mathrm{2PN}\leq -0.00045^\circ\,\mathrm{\,yr}^{-1}$. A recently released prediction of $\dot\omega^\mathrm{2PN}$ for \psrab\, is discussed.


Introduction
Recently, Hu et al. (2020) performed a detailed analysis of the perspectives of measuring, or effectively constraining, the moment of inertia (MOI) I A of the pulsar PSR J0737-3039A (Burgay et al. 2003;Lyne et al. 2004) by the end of the present decade by exploiting the general relativistic spin-orbit Lense-Thirring periastron precessionω LT, A (Damour & Schäfer 1988) induced by its spin angular momentum S A . Among the competing dynamical effects acting as potential sources of systematic uncertainty, Hu et al. (2020) included also the periastron precessionω 2PN to the second post-Newtonian (2PN) order which, along with the much larger 1 1PN precessioṅ depends only on the masses M A , M 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, µ GM is the gravitational parameter of the Double Pulsar made of the product of the Newtonian constant of gravitation G times the sum of the masses M M A + M 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 (Klioner & Kopeikin 1994), while is the osculating Keplerian mean motion. In particular, Hu et al. (2020, Table 1) reporteḋ 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 I A by Silva et al. (2021), would be of comparable magnitudeω It is clear that 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ω 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ω 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ω 2PN , the fact that also the formal 1PN term contributes it in a subtle way is often overlooked yielding wrong results.
The paper is organized as follows. In Section 2, the total 2PN 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 1PN + k 2PN . In particular, Equation (21) by Iorio (2021) is used as starting point in Section 2.1, where an error by Iorio (2021) 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 & Schäfer (1988), 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 which are determined for several binary pulsars is the fractional periastron shift per orbit k PK defined as .
In Equation (5), ω is the argument of periastron, ∆ω 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 · · · denote the average over the orbital period P PK which, in presence of PK accelerations, has to be meant as the anomalistic period P PK ano , i.e. the time span between two successive crossings of the (moving) periastron position. Nonetheless, it is common practice to deal with the averaged 2 periastron precessionω PK which is connected with k PK through where n PK 2π P PK (7) is the PK mean motion. In general, Equation (7) differs from Equation (2).

Starting from the formula by Iorio in osculating Keplerian orbital elements
As far as the 2PN periastron advance is concerned, Iorio (2021, Equation (21)) correctly calculated k 2PN , up to the scaling factor n K , 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 & Potapov (1994) with the same perturbative technique but a different calculational strategy, and by Damour & Schäfer (1988) who, instead, used the Hamilton-Jacobi method and the Damour-Deruelle (DD) parameterization (Damour & Deruelle 1985) which is nowadays routinely used in standard pulsar timing analyses Damour & Deruelle (1986); Damour & Taylor (1992). Iorio (2021), after having scaled k 2PN by Equation (2), erroneously claimed that the resulting expression for n K k 2PN is the total 2PN pericentre precession, which is not the case, as it will be shown below. Here, the explicit expressions of k 1PN , k 2PN in terms of the osculating Keplerian orbital elements are reported. They are where f 0 is the osculating numerical value of the true anomaly f at some arbitrary moment of time In order to correctly calculate the total 2PN pericentre precessionω 2PN , 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 where the osculating Keplerian period is and the 1PN correction, calculated according to the strategy followed by Iorio (2016), turns out to be with In the point particle limit corresponding to ν → 0, Equation (14) reduces to Equation (72) of Iorio (2016). The 1PN mean motion is, thus, Note that Equation (21) of Iorio (2021), i.e. the product of Equation (2) times Equation (9) n K k 2PN = 3 µ 5/2 2 − 4 ν + e 2 (1 + 10 ν) + 16 e (−2 + ν) cos f 0 4 c 4 a 7/2 1 − e 2 2 , has formally the dimensions of a pericentre precession of the order of O c −4 , but, contrary to what mistakenly claimed by Iorio (2021), it is not the total 2PN pericentre rateω 2PN . Indeed, the correct analytical expression for it can only be obtained by retaining the term of the order of O c −4 in the expansion in powers of c −1 of the product of Equation (15) times the sum of Equation (8) and Equation (9). If, on the one hand, replacing n 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 O c −4 . Indeed, from Equation (8) and Equation (15) which, added to Equation (16) This is the right analytical expression for the full 2PN pericentre precession expressed in terms of the osculating Keplerian orbital elements.
Recapitulating, on the one hand, Iorio (2021, Equation (21)) correctly worked out the 2PN fractional pericentre shift per orbit k 2PN up to n K as scaling factor. On the other hand, Iorio (2021), after having multiplied it by the osculating Keplerian mean motion n K , mistakenly claimed that the resulting expression for n K k 2PN was the total 2PN pericentre precession, missing a further contribution from the power expansion to the 2PN order of the product n 1PN k 1PN .

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 (Damour & Schäfer 1988, Equation (5.18 In Equation (19), x while e T and n DD are members of the Damour-Deruelle (DD) formalism (Damour & Deruelle 1985) which, in the limit c → ∞, reduce to the Keplerian eccentricity e and mean motion n K , as it will be shown below.
are numerically integrated over 1 yr with and without Equations (37)-(38) in each run, and time series of ω (t) are correspondingly calculated. In Equations (37)-(39), r is the relative distance between A and B,r is the versor of the relative position vector, v is the velocity vector of the relative motion, and v r v ·r is the radial velocity. Then, the difference between the Newtonian and the PN time series for ω (t) is taken to extract the time-dependent PN shift ∆ω PN (t). 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. In order to single out just the total 2PN effect (up to other PN contributions of higher order) ∆ω 2PN (t), the 1PN linear trend, analytically calculated by multiplying k 1PN of Equation (8) times n K t, is subtracted from the time series ∆ω PN (t). 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 ∆ω 2PN (t) 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 • yr −1 , can be compared with the lower panel of Figure 1 displaying the plot of Equation (18) as shown by Figure 2, or, equivalently,

Summary and conclusions
Calculating correctly the total 2PN periastron precessionω 2PN from the fractional periastron advance per orbit k PN = k 1PN + k 2PN requires to multiply the latter one by the 1PN mean motion n 1PN instead of the osculating Keplerian one n K , as incorrectly done by Iorio (2021), and to expand the resulting expression to the order of O c −4 . It remains true independently of the parameterization used. Adopting the osculating Keplerian orbital elements allows to obtain Equation (18) forω 2PN . 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 1PN + k 2PN by Damour & Schäfer (1988), 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.
The valueω 2PN = 0.000439 • yr −1 by Hu et al. (2020) comes from having neglected to multiply k 1PN by the 1PN mean motion and to expand the resulting product to the order of O c −4 .
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 yr −1 ≤ω 2PN ≤ 0.074 yr −1 .