Revisiting the 2PN Pericenter Precession in View of Possible Future Measurements

: At the second post-Newtonian (2PN) order, the secular pericenter precession ˙ ω 2PN of either a full two-body system made of well-detached non-rotating monopole masses of comparable size and a restricted two-body system composed of a point particle orbiting a ﬁxed central mass have been analytically computed so far with a variety of approaches. We offer our contribution by analytically computing ˙ ω 2PN in a perturbative way with the method of variation of elliptical elements by explicitly calculating both the direct contribution due to the 2PN acceleration A 2PN , and also an indirect part arising from the self-interaction of the 1PN acceleration A 1PN in the orbital average accounting for the instantaneous shifts induced by A 1PN itself. Explicit formulas are straightforwardly obtained for both the point particle and full two-body cases without recurring to simplifying assumptions on the eccentricity e . Two different numerical integrations of the equations of motion conﬁrm our analytical results for both the direct and indirect precessions. The values of the resulting effects for Mercury and some binary pulsars are confronted with the present-day level of experimental accuracies in measuring/constraining their pericenter precessions. The supermassive binary black hole in the BL Lac object OJ 287 is considered as well. A comparison with some of the results appeared in the literature is made.


Introduction
The problem of calculating at the second post-Newtonian (2PN) order of general relativity [1] the secular 1 precessionω 2PN of pericenter ω of a full two-body system made of a pair of detached, non-rotating masses of comparable sizes and of a restricted two-body system characterized by a test particle orbiting its massive primary has been analytically tackled several times so far with a variety of calculational approaches . Despite their formal elegance, it is not always easy to extract from them quickly understandable formulas, ready to be read and used in practical calculations in view of possible confrontation with actual data from astronomical and astrophysical scenarios of potential experimental interest. Perhaps, it is so because, e.g., of continuous references nested one inside the other to various papers pointing to a host of intermediate parameterizations, often of purely theoretical relevance that tend somehow to confuse a little bit at least some readers. Sometimes, they may wonder which numerical values of the parameters of the system under consideration out of those recorded in the literature have to be inserted in the equations. For a recent discussion on some aspects of the approaches followed in the literature so far, see Tucker & Will [19]; see also Klioner & Kopeikin [24] for a comparison of some of the parameterizations used in the literature to the 1PN level. 1 For the sake of simplicity, we will omit the brackets denoting the average over one orbital revolution here and throughout the paper.
Our aim is revisiting the issue of analytically calculating the 2PN pericenter precession by straightforwardly computing it perturbatively with the widely known method of variation of the orbital elements [25][26][27][28][29][30][31][32][33][34][35][36] in order to provide quickly understandable formulas, ready to be used in practical calculations in view of possible measurements in a not so far future, more likely in binary pulsars than in our Solar system, or to better model the dynamics of peculiar systems like, e.g., tight extrasolar planetary systems or the BL Lac object OJ 287 [37,38]. A similar strategy was adopted in Kopeikin & Potapov [12]. Because the actual data analyses of astronomical and astrophysical systems are performed by using the harmonic coordinates of PN theory, we will adopt them in our calculation (see the discussion in Section 4 of [19]). We will, first, deal with the point particle case (Section 2) by starting with the precession directly induced by the 2PN acceleration A 2PN entering the equations of motion (Section 2.1). Then, in Section 2.2, we will calculate the indirect 2PN precession arising from the fact that, to the order O c −4 , where c is the speed of light in vacuum, also the instantaneous shifts of the orbital elements occurring during an orbital revolution due to the 1PN acceleration A 1PN itself should be taken into account in the averaging procedure of the 1PN effects. Instead, neglecting such changes gives rise to the usual, time-honored Einstein-like 1PN precession. In principle, also other general relativistic precessions may be calculated, to the order O c −4 , from the mutual interaction of some 1PN accelerations induced by the bodies' mass and spin moments [39][40][41][42][43][44] entering the equations of motion; they will not be treated here because of their smallness. For some of them, see Iorio [45]. Section 2.3 contains numerical integrations of the equations of motion of some binary systems confirming our analytical result of Section 2.1 for the direct effect, and of Section 2.2 for the indirect one. It turns out that the direct 2PN perihelion precession of Mercury is smaller than the present-day observational accuracy in constraining any unmodeled perihelion precession of Mercury by about an order of magnitude or so. Currently, the 2PN equations of motion are not included in the dynamical models of the Solar system dynamics employed by the teams of astronomers producing the planetary ephemerides [46][47][48]. In Section 3 we repeat our calculation for a full two-body system by calculating both the direct (Section 3.1) and the indirect (Section 3.2) contributions to the 2PN pericenter precession in the same fashion as in Section 2. We compute them for the Hulse-Taylor binary pulsar PSR B1913+16 [49] and the double pulsar PSR J07373039A/B [50,51] by comparing the resulting predictions with the current experimental accuracy in determining their periastron precessions from timing measurements. While for PSR B1913+16 the overall 2PN periastron precession is already potentially measurable today, for PSR J07373039A/B the indirect contribution, which depends explicitly on the initial value of the orbital phase, may weaken or even cancel out the direct effect for certain values of the initial position of the pulsar along its orbit. On the other hand, for other initial positions the total 2PN periastron precession may be brought above the measurability threshold. We look also at the supermassive binary black hole in OJ 287. In Section 4, we compare our calculation with those in Kopeikin & Potapov [12] by disclosing an error in their results for the indirect effects. A comparison is made also with the results by Damour & Schafer [5]. Section 5 summarizes our findings, and offers our conclusions.

The Direct Pericenter Precession Due to the 2PN Acceleration
The 2PN acceleration experienced by a test particle orbiting a fixed body of mass M at distance r, written in harmonic coordinates, is (see, e.g., [23], Equation (2.3); [35], Equation (8.8.16), p. 332) In Equation (1), µ . = GM is the gravitational parameter of the primary, G is the Newtonian gravitational constant, and v r .
= v ·r is the radial velocity of the test particle, i.e., the projection of its velocity v onto the versorr of its position vector r with respect to the primary. Equation (1) can be obtained from the point particle limit of the 2PN equation of relative motion of a full two-body system treated in Section 3.1. Equation (1) can also be inferred from the equation of motion of Equation (4.4.18) of (Brumberg [27], p. 152) or Equation (1.5c) of (Damour & Schafer [5], p. 133) for the body 1 assumed as test particle orbiting the body 2 taken as its primary, i.e., for Let us analytically work out the direct long-term, i.e., averaged one orbital period P b , 2PN precession of pericenter induced solely by Equation (1) by means of the Gauss equations (e.g., [29,32,35]), valid for any additional acceleration A with respect to the Newtonian monopole A N = −µ/r 2 , where a, e, I, Ω, ω, f are the semimajor axis, eccentricity, inclination, longitude of the ascending node, argument of pericenter, and true anomaly, respectively, p . = a 1 − e 2 is the semilatus rectum, u .
= ω + f is the argument of latitude, n b . = µ/a 3 is the Keplerian mean motion, while A r , A τ , A ν are the radial, transverse and out-of-plane components of the extra-acceleration A, respectively. It is appropriate to remark that the Gauss equations are exact since the possible smallness of A with respect to A N is not assumed in their derivation ( [35], p. 108). In a perturbative calculation, which is fully adequate for the 2PN acceleration of Equation (1) in most of the situations in which a conceivable future detection could be envisaged (our Solar system, exoplanets, binary pulsars), the right-hand sides of Equations (2) and (3) have to be evaluated onto the Keplerian ellipse r = p/ (1 + e cos f ), assumed as unperturbed, reference trajectory, and averaged out over one orbital period P b . = 2π/n b by means of [26,27,32,[52][53][54][55]] In it, the derivatives of ω and Ω are given by Equations (2) and (3). To keep only terms of order O c −4 when Equation (1) is used in Equations (2) and (3), only the first term of Equation (4) has to be retained because of the presence of A itself in it through dΩ/dt, dω/dt. It is intended that in the following, the right-hand-sides of Equations (2)-(4) are evaluated onto the constant Keplerian ellipse; in order to avoid an excessively cumbersome notation, we avoid to append a subscript "K" to the orbital elements entering them.
The radial, transverse, and out-of-plane components of Equation (1), evaluated onto the reference Keplerian trajectory, turn out to be By inserting Equations (5)-(7) into Equations (2) and (3) and averaging with the first term of Equation (4) yields, to order O c −4 , the direct 2PN pericenter precessioṅ corresponding to a shift per orbit The analytical result of Equation (8) will be numerically confirmed in Section 2.3 by numerically integrating the equations of motion.

The Indirect Pericenter Precession due to the 1PN Acceleration
Equation (8), although directly inferred from the 2PN acceleration of Equation (1), does not exhaust the issue of calculating the full pericenter precession to the order O c −4 . Indeed, there are also other two contributions to it, which may be dubbed as "indirect", coming from the well-known 1PN acceleration itself (e.g., [35], p. 332) Basically, they arise because during an orbital revolution of the test particle under the perturbing influence of A like Equation (10) all the orbital elements, in principle, undergo instantaneous variations changing their values from their fixed Keplerian ones referred to some reference epoch t 0 . Moreover, when the integration over f is performed in order to obtain the net change per orbit, the fact that f is reckoned from a generally varying line of apsides because of A should be taken into account as well. Such features yield additional corrections of the order of O A 2 which, in the present case, are just of the order of O c −4 . We will implement such a strategy by following Iorio [45] in which the indirect effects of order O J 2 c −2 , where J 2 is the primary's oblateness, were computed in agreement with Will [56,57].
One of the aforementioned indirect contributions to the 2PN pericenter precession, marked conventionally with the superscript (I) in the following, is obtained from the orbital average of Equations (2) and (3), calculated with Equation (10), by means of the second and third terms of Equation (4) containing just Equation (10) itself which, among other things, shifts slowly the apsidal line from which the true anomaly f is counted. By recalling that the radial, transverse, and out-of-plane components of Equation (10) the resulting indirect precessionω Please note that Equation (14) is formally singular in the limit e → 0. The second indirect contributionω 2PN (II) indir comes from the fact that in general, when an extraacceleration A like, e.g., Equation (10) enters the equations of motion, all its orbital parameters undergo instantaneous changes during an orbital period. Usually, in standard first order calculations in A, such generally slow variations are neglected by assuming the Keplerian elements as fixed to some fiducial values at a reference epoch t 0 . Instead, accounting also for such changes yield further, indirect effects of the second order in A. The resulting indirect integrated shift over one orbit of any of the orbital elements = ω, can be calculated as where the superscript (2) indicates that the calculation is to the second order in A, {. . .} K denotes that the content of the curly brackets has to be evaluated onto the unperturbed Keplerian ellipse, and ∆φ j ( f 0 , f ) (1) , j = 1, . . . 5 are the instantaneous shifts experienced by the orbital elements during the orbital revolution. The latter ones are calculated as where the superscript (1) indicates that the shifts of Equation (16) are to the first order in A.
From [26,[53][54][55]] valid to the first order in A given, in the present case, by Equation (10), and Equations (11)- (13), it turns out that in the case of pericenter, only the 1PN instantaneous shifts of a and e induced by Equation (10) are required. By recalling that the Gauss equations for such orbital elements, to the first order in A, can be written as [26,[53][54][55]] and that the radial, transverse, and out-of-plane components of Equation (10) are given by Equations (11)-(13), it is straightforward to obtain They agree with, e.g., Equation ( (15), calculated for i = 5 by means of Equation (17) Please note that also Equation (22) is formally singular in e; moreover, it depends on the initial value of the true anomaly f 0 . The indirect total 2PN precessionω 2PN indir of order O c −4 is the sum of Equation (14) and Equation (22); it readṡ It should be noticed that Equation (23) is not singular for e → 0. On the other hand, Equation (23) is not univocally determined because of the presence of f 0 . In Section 2.3, we will confirm Equation (23) by numerically integrating the equations of motion for an arbitrary fictitious system.

A Numerical Confirmation of the Direct and Indirect 2PN Pericenter Precessions
The direct 2PN precession of Equation (8) was successfully confirmed by two numerical integrations of the equations of motion of, say, Mercury in the field of the Sun over 1 century (cty).
It is worthwhile recalling that the present-day level of accuracy in constraining any anomalous perihelion precession of such a planet with the most recent ephemerides, which all model the Solar system dynamics only up to the 1PN level in harmonic coordinates, may be at the level of σω 8 microarcseconds per century µas cty −1 , or, perhaps, 10 − 50 times worse; see the discussion in Iorio [58], and references therein.
In the first run, we simultaneously integrated the Hermean equations of motion, including the Newtonian monopole and the 2PN acceleration of Equation (1), in rectangular Cartesian coordinates along with the Gauss equations for all the Keplerian orbital elements over a time span 1 cty long starting from a set of initial conditions for the state vector of Mercury retrieved from the WEB interface HORIZONS, maintained by the NASA Jet Propulsion Laboratory (JPL). The resulting time series of the solution for ω (t), in blue, is displayed in Figure 1 along with a linear fit to it, in yellow.  (1) in addition to the Newtonian monopole, in Cartesian rectangular coordinates along with the Gauss equations for all its Keplerian orbital elements. A superimposed linear fit, in yellow, to the numerically integrated time series of ω is displayed as well. Its slope of 2.6 µas cty −1 agrees with the value obtainable analytically by calculating Equation (8) with the orbital parameters of Mercury. The initial conditions were retrieved from the WEB interface HORIZONS by the NASA Jet Propulsion Laboratory (JPL) which employs the same harmonic coordinates used in obtaining Equation (1) and Equation (10) to model the dynamics of the Solar system up to the 1PN level. The same plot, not displayed here, was obtained in a second numerical integration in which the Gauss equations were not included among the differential equations to be simultaneously solved.
The same plot was obtained in a second run in which the Gauss equations were not included in the numerical integration which was limited just to the equations of motion of Mercury in rectangular Cartesian coordinates, all the rest being the same as in the first run. Then, a time series for ω (t) was straightforwardly computed from the solutions obtained for the Cartesian coordinates x (t) , y (t) , z (t) of the planet by means of the standard conversion formulas for the Keplerian orbital elements. The resulting slope of the fitted linear trend amounts to 2.6 µas cty −1 , in agreement with the first run and Equation (8) calculated with the orbital parameters of Mercury. Interestingly, such a figure is only 3 times smaller than the previously quoted value of σω which, however, as already remarked, may be optimistic by a factor of 10-50.
It should be noted that at least in principle, the direct 2PN precession of Equation (8) should be measurable since it is due to a distinct acceleration, i.e., Equation (1), which may be suitably expressed in terms of a dedicated solve-for parameter to be estimated in a least-square sense in some covariance analyses. Instead, the indirect precession of Equation (23), since it comes from the 1PN acceleration of Equation (10) which is routinely modeled in the software of all the teams currently producing the planetary ephemerides, may not be detectable as a separate effect with respect to the other 1PN features of motion. Be that as it may, Equation (23) for 0 ≤ f 0 < 360 deg. It is possible to numerically confirm our analytical findings also for the indirect 2PN precession in the following way. First of all, a straightforward numerical integration of the equations of motion of a fictitious restricted two-body system to the 1PN level, i.e., by accounting only for the 1PN acceleration of Equation (10), shows that the simple secular trend arising from the celebrated 1PN Einstein-like pericenter precessionω 1PN does not match a linear fit to the time series obtained from the numerical integration. This is clearly shown in the upper panel of Figure 2 obtained for, say, f 0 = 0. It turns out that such a feature lingers even by changing f 0 from a run to another. It is crucial to note that our analytical result for the indirect 2PN precession of Equation (23), calculated with f 0 = 0, is able to fully explain the discrepancy between the slopes of the simple analytical 1PN trend due to Equation (25) where M = n b (t − t 0 ) + M 0 is the mean anomaly, M 0 is the mean anomaly at epoch, J k (se) is the Bessel function of the first kind of order k, and s max , j max are some values of the summation indexes s, j adequate for the desired accuracy level, is the continuous brown curve for δω depicted in the lower panel of Figure 2. It can be noticed that it does not vanish, and a linear fit to it, represented by the dashed red line in the lower panel of Figure 2, returns just the same value as Equation (23). Also, in this case, it occurs by varying f 0 .
dir of the pericenter of the relative motion of a two-body system can be straightforwardly computed from Equation (28) in the same fashion as for the point particle treated in Section 2.1. The radial, transverse, and out-of-plane components of Equation (28)  +e (4 + 11 η) cos f + 3 e 2 η (3 + 2η) sin 2 f , they reduce to Equations (5)-(7) in the point particle limit, i.e., for η → 0. By averaging the right-hand sides of Equations (1) and (2), calculated with Equations (29)- (31), with the first term of Equation (4) one finally obtainsω Equation (32) The current accuracy in measuring the periastron precession of the double pulsar is [61] σω = 0.00068 deg yr −1 .
An accuracy level of the order of Equation (33) should be reached in the forthcoming years thanks to new telescopes [62]. For the historical binary pulsar PSR B1913+16, whose relevant physical and orbital parameters are [63] while the most recent determination of its periastron rate is accurate to [63] For the supermassive binary black hole in OJ 287, whose relevant orbital parameters are [38]

A Comparison with Other Works
To the knowledge of the present author, the only other work in the literature making use of the method of the variation of constants and the Gauss equations is Kopeikin & Potapov [12]. As we will show, their result is incorrect because of the treatment of what are dubbed here as indirect effects.
Equation (5.2) of Kopeikin & Potapov [12], which we reproduce here to the benefit of the reader, is their main result. It is the total 2PN pericenter shift per orbit, in units of 2π, written in terms of the constants of integration k 1 , k 2 of the solutions of the Gauss equations for the semimajor axis and the eccentricity to the 1PN level. In our notation 3 , it is, in the test particle case, Since the constants of integrations k 1 , k 2 entering Equation (50) are determined with the initial conditions at t = t 0 , they contain explicitly f 0 ; thus, Equation (5.2) of Kopeikin & Potapov [12] actually does depend on the latter one, contrary to what, at first glance, someone could argue, perhaps mislead by the notation used by Kopeikin & Potapov [12] for k 1 , k 2 . By retrieving the explicit expression of k 1 , k 2 from Eqautions (20) and (21) k 1 = a + e µ 14 + 6 e 2 cos f 0 + e (4 + 5 cos 2 f 0 ) In Kopeikin & Potapov [12], it is k 1 → a 0 , k 2 → e 0 .
where a and e entering Equations (51) and (52) which does not agree with the corresponding expression for ∆ω 2PN tot /2π obtainable from the sum of our Equation (8) and Equation (23) by taking its ratio to n b .
From what can be deduced from the description of the method followed by Kopeikin & Potapov [12], the indirect effect corresponding to ourω 2PN (II) indir arises from the replacement a → a + ∆a ( f 0 , f ) 1PN , e → e + ∆e ( f 0 , f ) 1PN in 4 Equation (17), in a series expansion of it in powers of c −1 to the order c −4 , and in an integration of the resulting expression from f 0 to f 0 + 2π. The result, not explicitly shown by Kopeikin & Potapov [12], is ∆ω which does not agree with the corresponding expression from our Equation (22) [12], agree with the corresponding shifts from our Equation (8) and Equation (14) because their sum with Equation (54)   (55) It neatly disagrees with our numerical results of Section 2.3 since, for the fictitious system treated in Figure 2, Equation (55) provides a slope as little as −0.00255 deg cty −1 .
It may be interesting to make a comparison of our results also with the seminal results by Damour & Schafer [5], despite they did not use the Gauss equations. Damour & Schafer [5], following the example by Landau & Lifshitz [65], started from the Hamiltonian of the binary system in Arnowitt-Deser-Misner (ADM) coordinates [66] and adopted the Hamilton-Jacobi method. As far as the 2PN pericenter precession of a system of two mass monopoles is concerned, their main result is Equation (3.12) where h and E are the coordinate-invariant, reduced orbital angular momentum and energy, respectively. Its translation in terms of the parameters of the Damour-Deruelle (DD) parametrization [67] is given by where n is the PN mean motion ( [67], Equation (3.6d) ) and e t is one of the DD parameters [67]. Expressing Equation (56) in terms of the osculating Keplerian orbital elements can be made in the following two steps. First, E, h are to be written in terms of the DD parameters a r , e r by inverting Equations (3.6a) and Equation (3.6b) of Damour & Deruelle [67]. The result is Then, Equations (28) to (29) of Klioner & Kopeikin [24], which in general relativity, are with the aid of Equation (14) and Equation (16) of Klioner & Kopeikin [24], whose general relativistic expressions are are used to express a r , e r as functions of the osculating Keplerian elements a, e. We obtain for a r (a, e) , e r (a, e) 8 a −1 + e 2 e r = 4 e 2 a −1 + e 2 + µ c 2 −17 + 6 η + e 2 (2 + 4 η) + Finally, an expansion of the obtained expression in powers of c −1 to the 2PN level yields, in the point particle limit, Equation (53) which, as already noted, is incorrect. On the other hand, Equation (56) and Equation (57) seem to be mutually inconsistent since their expressions in terms of a, e do not even agree each other. Indeed, by using Equations (58), (66) and (67), Equation (30) of Klioner & Kopeikin [24], which, in general relativity, reads and Equation (65) to express e t in terms of a, e 8 a −1 + e 2 e t = 4 e 2 a −1 + e 2 + µ c 2 3 (−3 + η) + e 2 (−6 + 7 η) + which disagrees even with Equation (53) itself. By expanding Equations (53) and (70) in powers of e, it turns out that their disagreement is at the order O e 2 .

Summary and Conclusions
We analytically worked out the 2PN secular pericenter precessionω 2PN of both a test particle orbiting a static central body and a full two-body system made of a pair of comparable non-rotating monopole masses with the method of variation of orbital elements.
We, first, calculated the direct precessionω 2PN dir induced by the 2PN acceleration entering the equations of motion written in harmonic coordinates. Two different numerical integrations of the equations of motion of a point particle confirmed our analytical results. For Mercury moving in the field of Sun, it isω 2PN dir = 2.6 µas cty −1 . It is just 3 times smaller than the present-day formal accuracy σω = 8 µas cty −1 in constraining any unmodeled effect in the Hermean perihelion rate with the latest planetary ephemerides, although σω may be realistically up to 10-50 times worse. In the case of the binary pulsar PSR B1913+16, the direct 2PN periastron rate isω 2PN dir = 0.000038 deg yr −1 , to be compared with the most recent determination of its periastron rate σω = 0.000005 deg yr −1 , while for the double pulsar PSR J0737-3039A/B one hasω 2PN dir = 0.00019 deg yr −1 and σω = 0.00068 deg yr −1 . The direct 2PN perinigricon precession of the supermassive binary black hole in OJ 287 amounts to 11 deg cty −1 .
Then, we computed also the indirect 2PN pericenter precessionω 2PN indir arising from the fact that the 1PN acceleration actually changes instantaneously the semimajor axis and the eccentricity, and shifts the line of apsides instant by instant during one full orbital revolution. If properly accounted for in the orbital average, such features, which are of the second order in the acceleration causing them, gives rise to a further contribution of order O c −4 to the 2PN pericenter precession which adds on the direct one. The resulting expression turns out to be dependent on the initial position f 0 along the orbit. Numerical integrations of the equations of motion confirmed also such a result. Since the orbital dynamics of our Solar system is routinely modeled up to the 1PN level in harmonic coordinates of PN theory, it is unlikely that such an indirect precession can be measured separately because it does not come from a distinct acceleration which, instead, could be suitably expressed in terms of a dedicated solve-for parameter to be estimated in specific covariance analyses. For Mercury, its nominal size amounts to 16-33 µas cty −1 , depending on f 0 . For the binary pulsars, the experimental approach is different. It implies the determination, in a phenomenological, model-independent way, of several post-Keplerian parameters, among which there is also the periastron precession, from a confrontation of an analytical timing formula with the recorded pulses. Then, model-dependent analytical expressions for the measured post-Keplerian effects are used to determine the masses of the system, and to perform one or more tests of the model of gravitation considered. In the case of PSR B1913+16, the indirect 2PN precession ranges from −0.000048 deg yr −1 to 0.001052 deg yr −1 , while for PSR J07373039A/B it is 0.00092-0.00132 deg yr −1 . This shows that the choice of f 0 may enhance or even cancel out the overall 2PN periastron precession. For OJ 287, it ranges from 20 deg cty −1 to 516 deg cty −1 ; the 1PN perinigricon precession amounts to 206.8 deg cty −1 .
We compared our formulas to some other analytical results in the literature by showing that the latter ones disagree with ours and with our numerical integrations of the equations of motion. It appears that the source of discrepancy relies in the treatment of the indirect effects arising from the inclusion of the instantaneous 1PN changes of the semimajor axis and eccentricity in the integration over one orbital revolution of the pericenter shift due to the 1PN acceleration itself.
Funding: This research received no external funding.