Detecting the General Relativistic Orbital Precession of the Exoplanet HD 80606b

We investigate the relativistic effects in the orbital motion of the exoplanet HD80606b with a high eccentricity of $e\simeq0.93$ .We propose a method to detect these effects (notably the orbital precession) based on measuring the successive eclipse and transit times of the exoplanet. In the case of HD80606b, we find that in ten years (after approximately 33 periods) the instants of transits and eclipses are delayed with respect to the Newtonian prediction by about three minutes due to relativistic effects. These effects can be detected by comparing at different epochs the time difference between a transit and the preceding eclipse, and should be measurable by comparing events already observed on HD80606 in 2010 with the Spitzer satellite together with those to be observed in the future with the James Webb Space Telescope.


Introduction
Whereas more than 4000 exoplanets have been discovered so far, the planet orbiting the G5 star HD 80606 remains a remarkable and unique case. Its eccentricity is extremely high: e 0.93. HD 20782b is the only planet reported to have a higher eccentricity (O'Toole et al. 2009) but its high e-value stands on one single measurement and has not been confirmed up to now. The high eccentricity of HD 80606b was well established as soon as the planet was discovered through radial velocity measurements (Naef et al. 2001) and has been largely confirmed by subsequent observations. Today the most accurate system parameters are those determined by Hébrard et al. (2011), who provide e = 0.9330 ± 0.0005.
The long orbital period of the planet (P = 111.4367 ± 0.0004 days) implied a low probability of the orbital plan to be aligned with the line of sight to the Earth. Still, HD 80606b was successively discovered to pass behind its parent star (planetary eclipses) by Laughlin et al. (2009), then in front of it (planetary transits) by Moutou et al. (2009) (see also Garcia-Melendo & McCullough 2009;Fossey et al. 2009;Winn et al. 2009;Hidas et al. 2010). That configuration allows the planetary radius to be measured (R p = 0.981±0.023 R Jup ); the inclination of the orbit is also known (i = 89.269 • ±0.018 • ), which provides the mass M p = 4.08 ± 0.14 M Jup from the sky-projected mass derived with radial velocities. The obliquity of the system could also be measured from the Rossiter-McLaughlin anomaly observed during transits, which revealed that the orbit of HD 80606b is prograde but inclined (Moutou et al. 2009;Pont et al. 2009;Winn et al. 2009);Hébrard et al. (2011) reported an obliquity of λ = 42 • ±8 • . Such a peculiar orbit may be explained by the influence of the companion star HD 80607, located 1200 AU further away, through a Kozai-Lidov mechanism (Pont et al. 2009;Correia et al. 2011).
The particularly high eccentricity of the planetary orbit raises the question of the feasibility of detecting the general relativistic (GR) precession of its periastron. Indeed it should be enhanced to ∆ GR = 6πGM ac 2 (1 − e 2 ) 215 arcsec/century , as compared to 43 arcsec/century for Mercury (see also Correia et al. 2011). Besides Mercury, the relativistic precession has been measured in the solar system for some asteroids like Icarus (Gilvarry 1953). The tests done in the solar system agree with GR to within 10 −4 (Will 1993). Outside the solar system the effect is well known in binary pulsars like the Hulse-Taylor pulsar (Hulse & Taylor 1975) where it has been measured with high precision (see Taylor & Weisberg 1982), and in other binary pulsars such as the double pulsar (Kramer et al. 2006). To a further distance, there are now attempts to measure the relativistic precession of the star S2 close to the galactic centre (Parsa et al. 2017). The measure of the relativistic periastron advance in stellar systems was proposed by Gimenez (1985), but it is quite complicated to discriminate it from the tidal precession (see e.g. Wolf et al. (2010)). Figueira et al. (2016) reported that the relativistic precession was not detectable in the HD 80606 system in spite of the large dataset secured over several years. They measured the variation of the longitude of the periastron using the whole available datasets:ω = 9 720 ± 11 160 arcsec/century, thus in a non-significant way. Indeed, the accuracy reached on ω mainly from radial velocities is ± 540 arcsec (Hébrard et al. 2011), which is not accurate enough.
The transiting nature of HD 80606b allows a better accuracy to be reached thanks to the precise timing that can be measured on the mid-point of each transit or eclipse. Typical accuracies of ± 85 sec (Hébrard et al. 2011) and ± 260 sec (Laughlin et al. 2009) were obtained on the mid-point of transits and eclipses, respectively, using the Spitzer Space Telescope and its IRAC instruments. Indeed, only a satellite on an Earth-trailing heliocentric orbit such as Spitzer could allow the whole 12-hour-long transit of HD 80606b to be continuously observed together with sufficiently long off-transit references immediately before and after the event.
In this paper we propose to detect the GR precession of the periastron (and the relativistic orbital motion) for the exoplanet HD 80606b by using a method based on the transit times, computing a small drift in the successive instants of transits due to relativity. We find that in ten years (after 33 periods), due to the relativistic effects, the instants of transits are delayed with respect to the Newtonian predicted ones by about three minutes. Whereas the transit mid-time could be measured with such an accuracy, the absolute value of that shift could not be detected because the system's orbital parameters are not known today with a high enough accuracy.
However, we argue that the effect could actually be measurable from the observation of the change of the elapsed time t tr−ec = t tr − t ec between an eclipse and the next transit. 1 The time between an eclipse and a transit was measured in January 2010 using Spitzer to be t tr−ec (Jan. 2010) 5.9 days , with a good accuracy of ± 275 sec. Due to the GR precession of the periastron, and the modification of the period due to relativity, we find that this time difference will be reduced by 183 sec in 2020 (after 33 orbits), and by 271 sec in 2024 (49 orbits). Such shift would be difficult to detect in 2020 using IRAC on Spitzer, but should be detectable a few years later using NIRCam or MIRI on James Webb Space Telescope (JWST). Most of the paper is devoted to theoretical aspects, and it ends with a conclusion and discussion in Sect. 6 where we give further comments on previous observations and on the Newtonian perturbing effects. In Sect. 2 we introduce the geometrical conventions used to compute the times of transits and eclipses. Those are computed in two different ways: a post-Keplerian parametrisation of the orbit in Sect. 3, and a Hamiltonian method using Delaunay-Poincaré canonical variables in Sect. 4. Those methods give the exact solution of the first post-Newtonian equations of motion. Our final results for the transit times are given in Tables 1, 2, and 3 below. Furthermore, as a check of the latter methods, we also computed the transit times with a method relying on celestial-mechanics perturbation theory in Sect. 5. As we used it, this procedure is valid on average over one orbit, and takes only into account the secular effects. We find that it is in good agreement with the post-Keplerian and the Hamiltonian methods. A quick but rough estimate of the GR shift of transit times is also given in Appendix A.

Geometry of the planetary transits
We introduce a frame {x, y, z} with origin at the centre of the star. The observer is in the direction x (Earth's direction), while the projected angular momentum of the star is by convention along the direction z. Hence the plane {y, z} is the plane of the sky. The orbital plane of the planet is defined by its inclination I with respect to the plane {x, y}, and by its longitude Ω with respect to the plane {x, z}, following our conventions described in Fig. 1a. As depicted in Fig. 1c, the motion of the planet in the orbital plane is parametrised by polar coordinates (r, ϕ), where ϕ denotes the sum of the true anomaly and of ω 0 , the latter denoting a constant angle, which can be viewed as defining the initial argument of the periastron at a reference time. Then the coordinates (x, y, z) of the planet in this frame read (see Figs. 1a and 1c) x = r − cos ϕ cos Ω + sin ϕ cos I sin Ω , Since x is the direction of the observer, the condition for the (partial or complete) transit of the planet is that y 2 + z 2 (R + R p ) 2 , where R and R p are the radius of the star and planet, hence sin 2 I sin 2 ϕ + cos I cos Ω sin ϕ + sin Ω cos ϕ Since R +R p r, it is easy to see that this condition can be satisfied when sin I 1 or sin ϕ 1, and | cos I cos Ω sin ϕ+ sin Ω cos ϕ| 1. With the conventions of Fig. 1a we are interested in a transit for which ϕ π and Ω 0. 2 On the other hand, for an eclipse during which the planet is passing beyond the star as seen from the observer, we have ϕ 0 and Ω 0. In an expansion to the first order in R +Rp r 1, we find that the coordinates of the planet during and around the transit are parameter is given by (see Fig. 1b where r(π) and r(0) are the values of the radial coordinate of the planet at ϕ = π and ϕ = 0, when the transit and eclipse occur. Since the eclipse will be seen from behind the star and so much farther away, we shall have to add to the instants of eclipse a Roemer time delay accounting for the propagation of light at finite velocity c, and which is a 0.5PN effect ∝ 1/c. Thust i −→t i + ∆t R where ∆t R = r(0) + r(π) c 2.7 min (for HD 80606b) .
We define the positions of the planet {T 1 , T 2 , T 3 , T 4 } to be respectively the entry of the transit, the start of the full transit, and the exits of the full and partial transits. Similarly the points {T 1 ,T 2 ,T 3 ,T 4 } corresponding to the eclipse are shown in Fig. 1c. The transit and eclipse conditions for the polar coordinates (r, ϕ) at the passage of these points read Given a model r(ϕ) for the planetary orbit, this determines the true anomaly ϕ for each of these points, where the vertical coordinate Y (b) is given in terms of the impact parameter b by (see Figs. 1b and 1c) The cases for which (R ± R p ) 2 = b 2 correspond to the planet just grazing tangentially the star, either from the exterior or from the interior of the stellar disc. In addition to Eqs. (8)-(9) we compute also the points corresponding to the minimal approach to the center of the star, T m andT m , whose transit conditions are obviously (see Fig. 1b) Having determined the polar coordinates of the planet we shall deduce the instants t i andt i of passage at each of these points by using the law of motion of the planet on the trajectory.

Post-Newtonian motion of the planet
In this section we compute the times of transits and eclipses, that is, the instants of the successive passages at the positions T i andT i defined by the geometric conditions in Eqs. (8)- (10), in a relativistic model, including the first post-Newtonian (1PN) corrections beyond a Keplerian orbit. More precisely, we compute the difference between the transit times predicted by the relativistic 1PN model and those predicted by the Newtonian model. This will permit us to assess whether the latter relativistic corrections could be detectable, and we propose an observable quantity, directly measurable in future observations of the planet HD 80606b.

Keplerian parametrisation of the orbit
In the Keplerian model the orbit of the planet and its motion on the orbit are given by the usual six orbital elements {a, e, I, , ω 0 , Ω}, where I and Ω specify the orbital plane (see Fig. 1a), a is the semi-major axis, e is the eccentricity, ω 0 denotes the constant angular position of the periastron, and the mean anomaly is defined by = n 0 (t − t 0,P ), where t 0,P is the instant of passage at the periastron and n 0 = 2π/P 0 is the mean motion, with P 0 the period. The Keplerian relative orbit in polar coordinates (r, ϕ), is conveniently parametrised by the eccentric anomaly ψ as The relative orbit in polar coordinates reads from which we deduce the impact parameter in the Newtonian model, say b 0 , to be In the Newtonian model we have Kepler's third law n 0 = (GM/a 3 ) 1/2 , where M = M + M p is the total mass of the star plus planet system. Furthermore the energy and the angular momentum of the orbit are given by the usual Keplerian formulas where The conditions for the transits and eclipses are given by Eqs. (8) with (9), and determine the values of the eccentric anomaly at these points. In the Newtonian analysis the eccentric anomaly ψ 0 is obtained by solving where Y 0 = Y (b 0 ) denotes the vertical coordinate (Eq. (9)) computed with the Newtonian impact parameter given by Eq. (13).

Relativistic corrections to the Keplerian motion
Next we consider the relativistic 1PN model (see e.g. Wagoner & Will 1976). For this we shall use an explicit solution of the equations of motion that is particularly elegant, called the quasi-Keplerian representation of the 1PN motion (Damour & Deruelle 1985). In this representation the orbit and motion are parametrised by the eccentric anomaly ψ in a Keplerianlooking form as where = n(t − t P ) and n = 2π/P , with now P being the period of the relativistic model. The motion is parametrised by six constants: the mean motion n, a particular definition of semi-major axis a r differing from a in the Newtonian model by a small 1PN term, the precession of the orbit K including the relativistic precession, and three types of eccentricities e r , e t , and e ϕ differing from e and from each other by small 1PN terms. We shall pose The constant k denotes the relativistic advance of the periastron per orbital revolution, which is k = ∆ 2π where ∆ is the angle of return to the periastron, and whose GR value has been given in Eq. (1). The effect of precession can be seen clearly with the expression of the relative orbit in polar coordinates, which takes the more involved form (Damour & Deruelle 1985) Besides the usual contribution of the precession ∝ ϕ−ω0 K , we note the presence of a term vanishing in the small mass ratio limit, when ν = µ/M → 0.
Crucial to the post-Keplerian formalism are the explicit expressions of all these constants in terms of the conserved energy E and angular momentum J of the orbit. Since we want to compare the relativistic model of Eq. (16) to the Newtonian one of Eq. (11) we will compare, for convenience, orbits that have the same energy E and angular momentum J. Thus we need to express the relations in terms of the "Newtonian" semi-major axis a and eccentricity e defined in terms of E and J by the Newtonian formulas in Eq. (14). We have (Damour & Deruelle 1985) (see also Blanchet 2014) Article number, page 5 of 15 Among the effects described by these constants, the main ones are those associated with the mean motion n and the periastron precession K. The expressions for n and K in terms of E and J are invariant, meaning they do not depend on the coordinate system. Since we expect the effect of the precession of the orbit to be dominant over the shift in the orbital period. The other relations, linking the semi-major axis and eccentricities to E and J, depend on the coordinate system; they are given here in harmonic coordinates. In this paper we compute an observable quantity (defined by Eq. (30) below), which thus does not depend on the choice of coordinate system. We now investigate the difference between the analysis using the relativistic model and the Newtonian one concerning the planetary transits. The conditions for the transits and eclipses are still given by Eqs. (8) with (9), with the impact parameter given by Eq. (6). However, the impact parameter will differ from the Newtonian value by an amount δb. We readily find this amount thanks to the explicit formula Eq. (18), using a first order expansion in the small parameters given by Eq. (19), where ϕ N = (2N + 1)π − ω 0 andφ N = 2N π − ω 0 . Here N is the number of orbits since the reference time, with N = 0 corresponding to the transit observed in January 2010. Then the coordinate Y (b) is modified by δY = dY db | b0 δb and we find (see Eq. (9) Naturally for T m andT m , Eq. (10), we just have δY = δb cot I.
We are now in a position to solve the transit conditions and find the modification of the eccentric anomaly δψ with respect to the Newtonian model, that is, the solution of Eq. (15). Among all the relativistic effects contributing to δψ, only the one associated with the relativistic precession k is "secular", and therefore grows when ψ 0 −→ ψ 0 + 2πN proportionally to the number of orbits N . The other effects, associated with the relativistic parameters ξ, ε r , and ε ϕ , are periodic and therefore average to zero in the long term. Accordingly we split Similarly, we distinguish δY = δY secular + δY periodic , where the secular part is given by the N -dependent terms in Eq. (21). Perturbing Eq. (8) around the solution ψ 0 of the Newtonian condition of Eq. (15), we explicitly find Finally, having computed δψ with the transit condition we readily obtain the time difference between the instants of transits (and eclipses) in the relativistic model and in the Newtonian model by using the equation for the mean anomaly = n(t − t P ) = ψ − e t sin ψ. We notice that there is another secular term associated with the correction in the period, n = n 0 (1 + ζ), where ζ is the second invariant of the motion besides the relativistic precession k. We thus write similarly δt = δt secular + δt periodic and find When computing δt, the 1.5PN shift in the Roemer effect also has to be taken in account, although we find that this effect is very small. With the help of Eq. (7), it comes after N orbits Finally let us comment on the approximation we have made here, namely the expansion to first order in the relativistic parameters k, ζ, ξ, and so on. This approximation is valid as long as the number of orbits N over which we integrate is much smaller than the inverse of these relativistic effects, for example as long as N 1/k. In the case of HD 80606b we have k 5 · 10 −7 and we are considering N = 33 orbits between the transit of January 2010 and that of 2020 (see the Introduction), so our approximation is justified.
We present our results 5 in Tables 1, 2, and 3. From the above analysis we compute for each orbit N the quantities (with i ∈ {1, 2, m, 3, 4}) where t 0,i (N ) andt 0,i (N ) denote the Newtonian values. Obviously we have t 0,i (N ) = t 0,i (0) + N P 0 andt 0,i (N ) = t 0,i (0) + N P 0 since the Newtonian orbit is fixed, and the motion periodic with period P 0 = 2π/n 0 . Together with the shift of instants of transit, we compute also the relativistic effect on the total duration and the entrance of the N -th transit, as defined by together with similar quantities δt 14 (N ) and δt 12 (N ) for the N -th eclipse. The results are given in Tables 1 and 2. number N of transit δt 1 (N ) δt 2 (N ) δt m (N ) δt 3 (N ) δt 4 (N ) δt 14 (N ) δt 12 (N ) 0 (January 2010 As we see in Table 1, the relativistic model predicts that after ten years the planet will begin the transit about three minutes earlier than the Newtonian prediction, and about 4.5 minutes earlier after 15 years. Given that the precision on the measured instants of transits can be as good as 85 seconds (Hébrard et al. 2011), this is already an indication number N of eclipse δt 1 (N ) δt 2 (N ) δt m (N ) δt 3 (N ) δt 4 (N ) δt 14 (N ) δt 12 (N ) 0 (January 2010) 7.3 · 10 −3 7.3 · 10 −3 8.6 · 10 −3 8.2 · 10 −3 8.3 · 10 −3 1.0 · 10 −3 6.0 · 10  that the relativistic effects might be measurable if the orbital parameters were particularly well known. However, the ten-years-long relativistic effect on the transit duration δt 14 and transit entrance δt 12 are on the order of seconds, and seem a priori hardly detectable. Concerning the eclipse in Table 2, the effect is much smaller, of the order of ten seconds after ten years with respect to the Newtonian model. As the effect on the transit and the eclipse timings is different, this provides a way to directly detect it.  Table 3: Differences (in seconds) between the time interval t tr−ec (N ) between the passage at the minimum approach point during the N -th eclipse and the passage at the minimum approach point during the N -th transit, and the time interval t tr−ec (0) at the reference period, when the eclipse and transit where measured in January 2010 with the Spitzer satellite (Hébrard et al. 2011). The values for the next five years are also given. The 30th transit occurred in March 2019 and the 49th will occur in December 2024.
Next we define t tr−ec (N ) to be the time interval between the passage at the minimum approach point during the N -th eclipse and the passage at the minimum approach point during the N -th transit: By a straightforward calculation, using the fact that the Newtonian orbit is fixed, we obtain The important point about this result is that the right-hand side of the equation represents the relativistic prediction, which is tabulated in Tables 1 and 2, while the left-hand side is directly measurable: namely t tr−ec (0) has been measured with good precision between the eclipse and the transit of 2010, while future observation campaigns could allow t tr−ec (N ) to be measured ten to fifteen years later. From Table 3 we predict that the observations in 2020 will measure that the time difference between the transit and eclipse (with N = 33) is shorter by 182.8 seconds (approximately three minutes), because of the relativistic effect, as compared to the one which was measured ten years ago in 2010. The values of this observable quantity given by Eq. (30) are also presented for the transits of the next seven years in Table 3; in December 2024, the difference between the eclipse and transit (with N = 49) is shorter by 271.4 seconds (4.5 minutes) as compared to the one measured in 2010.
The 2010 measurement yielded (Laughlin et al. 2009;Hébrard et al. 2011) t tr−ec (0) = (5.8491 ± 0.003) days , with good published uncertainty in the measurement, of the order of 4.5 minutes, comparable to the relativistic effect we want to detect. Therefore we conclude that for the next observations of the transit/eclipse in the coming years, we should already be on the verge of detecting the relativistic effect.

Hamiltonian integration of the motion
At the 1PN order, the Hamiltonian of the relative motion of two point masses reads (Blanchet 2014) where P = P1+P2 µ is the reduced linear momentum, which is conjugate to X, and we pose P R = P · X/R with R = |X|, and P 2 = P 2 . We shall use the Delaunay-Poincaré canonical variables. We start with the usual orbital elements {a, e, , ω}, where a and e are the semi-major axis and eccentricity, ω is the argument of the periastron, and = n(t − t P ) denotes the mean anomaly, where n = GM/a 3 by definition. 6 Since the motion takes place in a fixed orbital plane, we do not need to consider the inclination I and the longitude of the node Ω. Thus we have together with = ψ − e sin ψ, where ψ is the eccentric anomaly. Then Eq. (32) reads where we introduced X = 1 − e cos ψ for convenience. The conjugate Delaunay-Poincaré variables {λ, Λ, h, H} are then defined by and the Hamiltonian equations take the canonical form With the explicit expression of the 1PN Hamiltonian in Eq. (34) we obtain dh dt = GM n √ 1 − e 2 2c 2 e 2 aX 2 ν − 4 + (4 + ν)(1 − e 2 ) + 12 Let us now turn to the integration of the latter equations. As we see from Eq. (37) we need to obtain the indefinite integrals I n (ψ) = dψ X n and J n (ψ) = dψ sin ψ X n . Looking at Gradshteyn & Ryzhik (1980), one finds all the needed integrals: With these elementary integrals it is straightforward to obtain the complete solution of the motion. We first compute Λ(ψ) and H(ψ). To the zero-th order their values Λ 0 and H 0 are constant, so the solution is described by constant orbital elements a 0 , e 0 , and n 0 = (GM/a 3 0 ) 1/2 , and we can integrate using n 0 dt = X 0 dψ to the requested order, where X 0 = 1 − e 0 cos ψ. Using the definitions of Λ and H in Eq. (35) we obtain the solution as a = a 0 + δa and e = e 0 + δe, where In exactly the same way, dh/dt can be integrated to give ω = ω 0 + δω, with (40a) Finally, we integrate the mean anomaly = ψ − e sin ψ = λ + h. The point is that in the equation for λ, we first need to express the mean motion as n = n 0 + δn, where δn = − 3n0 2a0 δa and δa has been found in Eq. (39). We thus obtain = 0 + δ , where 0 = n 0 (t − t 0,P ) = ψ − e 0 sin ψ (with t 0,P the constant passage at periastron), and We note that we could separate, as in Sect. 3, the secular contributions given by the first terms in both δω and δ from the manifestly periodic ones. Combining all those equations, we can solve the transit condition to Newtonian order (Eq. (15)), which gives ψ 0 , and then the correction δψ coming from Eq. (22), where the modification of the impact parameter is given by perturbing Eq. (13): cos ω 0 δe 1 + e 0 cos ω 0 + sin ω 0 δω 1 + e 0 cos ω 0 (eclipsesT i ) .
Finally this gives the modification of the instants of transits due to the 1PN corrections as where now X 0 = 1 − e 0 cos ψ 0 . This method gives the same results as the post-Keplerian ones, reported in Tables 1 and 2, up to 0.5% after 33 cycles. This difference can be explained by comparing Eqs. (21) and (42): the computation of δb slightly differs depending on the method used. Moreover by replacing roughly ϕ N → ϕ N /(1 + k) in Eq. (21), the second order perturbations can be estimated to play a role at the level of ∼ 5 · 10 −6 .

Lagrangian perturbation theory
It is instructive to work out the problem by means of a Lagrangian perturbation method. In this case we shall focus on secular effects (on a timescale much longer than the orbital period) after averaging the perturbation equations. The Lagrangian of the relative motion of two point masses at the 1PN order is given by where the 1PN perturbation function R depends on the relative position x and velocity u = dx/dt of the particles (we pose r = |x| andṙ = dr/dt), and reads (see e.g. Blanchet & Iyer 2003), In the Lagrangian perturbation formalism, when looking at secular effects, we can apply the usual perturbation equations of celestial mechanics directly with the perturbation function R in the Lagrangian (see e.g. Damour & Esposito-Farèse 1994). As in Sect. 4 we choose the independent orbital elements to be {a, e, , ω}, and we pose n = (GM/a 3 ) 1/2 . The perturbation (Eq. (48) It must be remembered that in perturbation theory n is not equal to 2π/P with P the period. Indeed we have the usual definition = n(t − t P ) but the instant of passage at periastron depends on time: t P = t P (t). Instead the period will be given by averaging Eq. (46c) for . As in Sect. 4 we denote the constant orbital elements to zero-th order by a 0 , e 0 , and pose n 0 = (GM/a 3 0 ) 1/2 . To first order the equations for a and e can readily be integrated as where the perturbation function reduces in this case to (with X 0 = 1 − e 0 cos ψ) The results of Eqs. (47)-(48) coincide exactly with those of the Hamiltonian formalism, Eqs. (39). Next, the equations for and ω to first order become Finally we have to perform the orbital average of these equations, f = 1 P P 0 dt f (t). Either we can compute the average directly from Eqs. (49), or in a simpler way we can substitute R by its average R in Eqs. (49). We have hence we end up with the two main results These are in agreement with the relative modification of the mean motion or period P given by the relativistic parameter ζ in Eq. (19a), and of course with the relativistic precession parameter k in Eq. (19b). Thus the two averaged Eqs. (51) reproduce the two gauge-invariant equations of the quasi-Keplerian formalism in Sect. 3. They can also be obtained by averaging the equations of the Hamiltonian formalism in Sect. 4. The times of transit and eclipse are straightforwardly computed in perturbation theory using Eqs. (51) and the transit conditions of Eq. (8). In this way we have confirmed the results reported in Tables 1, 2, and 3, but there is a typical difference of nearly 2% after 33 cycles, due to the fact that in this perturbative approach we neglect the periodic terms. As we discussed in Sect. 3, the relativistic effect on the times of transit after a large number of orbits essentially depends on the precession k, and the modification of the orbital period ζ, with the latter effect being smaller (see Eq. (20)).

Conclusion and discussion
In this paper we investigated the relativistic effects in the orbital motion of the high eccentricity exoplanet HD 80606b orbiting around the G5 star HD 80606. We proposed a method to detect these effects, based on the accurate measurement of the elapsed time between a mid-transit instant of the planet passing in front of the parent star, and the preceding mid-eclipse instant when it passes behind the star.
We presented different computations of the relativistic effects on the transit and eclipse times. One is based on the post-Keplerian parametrisation of the orbit in Sect. 3, another one is a Hamiltonian method with Delaunay-Poincaré canonical variables in Sect. 4, and finally a Lagrangian perturbation approach for computing the secular effects in Sect. 5. A crude understanding and estimation of the main effect is also provided in Appendix A. These methods gave consistent results, which are reported in Tables 1, 2, and 3.
We found that in ten to fifteen years, corresponding to 33 to 49 orbital periods, the time difference between the eclipse and the next transit is reduced by approximately 3 to 4.5 minutes due to the relativistic effects for this planet.
We conclude that these effects should be detectable for the next observations of the full transit and eclipse of HD 80606b in coming years by comparing to previous observations done in 2010.
By comparison to short-period exoplanets, only a few long-duration eclipses and transits of HD 80606b have been observed today. Laughlin et al. (2009) (Shporer et al. 2010) and could only observe partial portions of transits with a given telescope; so their accuracies on the mid-transit times were ±310 seconds or even poorer.
Thus, the best available measurement of t tr−ec was obtained in January 2010 with Spitzer with an accuracy of ± 275 seconds. A similar accuracy could be reached in 2020 with Spitzer, but it would be probably insufficient to detect the effect, slightly below 200 seconds at that time. Thanks to its 6.5-m aperture diameter, the James Webb Space Telescope (JWST) would provide an accuracy on t tr−ec a few times better than Spitzer (0.85-m diameter). So after JWST has started its operation, hopefully in 2021, it should be feasible to detect the effect with its NIRCam or MIRI instruments, while the effect will approach 300 seconds by comparison to 2010. Spitzer and JWST are among the rare telescopes able to continuously observe the 12-hour duration transit of HD 80606b, which is mandatory to reach a good timing accuracy. Other current or future telescopes such as TESS, CHEOPS, or PLATO could also sample long-duration transits; however, they are smaller and less sensitive than Spitzer or JWST, so are not able to reach their accuracy on timing. In addition they operate in the optical domain, whereas Spitzer and JWST observe in infrared where eclipses are deeper and thus more accurately measured. JWST clearly is the best facility to attempt that detection.
Finally we quickly discuss the other perturbing effects that may affect the measurement. HD 80606b is the only planetary companion detected today in that system. More than 15 years of radial-velocity monitoring 7 did not provide a hint of any additional companion (Hébrard et al. 2011;Figueira et al. 2016). With the available data, one can put a 3-σ upper limit of 0.4 m/s/yr of a possible linear drift in addition to the radial velocity signature of HD 80606b; this allows any additional planetary companion with a sky-projected mass larger than 0.5 Jupiter mass to be excluded with an orbital period shorter than 40 years around HD 80606. On the other hand, the companion star HD 80607 is too far away (1200 AU) to produce a noticeable effect on the orbital precession of HD 80606b.
Besides the possibility that there might be some disturbing other bodies, we have to consider the precession of the orbit caused by the oblateness of the star, the tidal effects between the star and the planet, and the Lense-Thirring effect. The orbital precession rate due to the quadrupolar deformation J 2 of the star is ∆ J2 = 3πJ 2 R 2 a 2 (1 − e 2 ) 2 .
Assuming the solar value J 2 ∼ 10 −7 we obtain ∆ J2 ∼ 0.4 arcsec/century for HD 80606b, which is negligible for our purpose. As for the orbital precession induced by the tidal interaction of the planet with its parent star, it produces a Appendix A: Rough estimation of the effect It is also possible to estimate the PN effects by considering a step-by-step shift of the orbit along an approximate sequence of Keplerian ellipses. As depicted in Fig. A.1, at each step we have to consider two time delays: the one induced by the new position of the transit on the orbit, δt k , and the time to pass from the previous periastron to the current one, δt P . The first delay is due to the different orientation of the N -th Keplerian ellipse with respect to the (N − 1)-th one, shifted by the relativistic precession angle ω N − ω N −1 2πk (see Fig. A.1). It simply reads n 0 δt k N − N −1 , where n 0 is the Newtonian mean motion and N the mean anomaly solving the transit condition of the N -th orbit (for a given transit T i or eclipseT i ). Of course the mean anomaly is counted from the periastron P N of the N -th orbit, and is given in terms of the eccentric anomaly by Kepler's equation N = ψ N − e sin ψ N . The condition relating the corresponding true anomalies of the successive ellipses is thus arctan 1 + e 1 − e tan ψ N 2 − arctan 1 + e 1 − e tan ψ N −1 2 + πk 0 . (A.1) As for the second time delay, it is due to the modification of the mean motion (and thus the period) as given by n = n 0 (1 + ζ) (see Eq. (19a)). We have = n 0 (t − t P N ), which we can roughly equate to n(t − t P N −1 ), where t P N −1 and t P N are the successive instants of passage at periastron. Thus t P N − t P N −1 −ζ(t − t P N −1 ), and we roughly approximate t − t P N −1 by the period so that n 0 δt P −2πζ. The time delay between the N -th orbit and the reference one finally becomes It is clear from Fig. A.1 that the two effects are negative, and add up negatively in the total delay. Even if this method seems quite crude, we find that it gives a reasonable estimate of the time shifts. Indeed, the computed times differ from the "exact" post-Keplerian results reported in Tables 1-3 by only 9% after 33 cycles.