The apparent tidal decay of WASP-4 b can be explained by the Rømer effect

: Tidal orbital decay plays a vital role in the evolution of hot Jupiter systems. As of now, this was only observationally confirmed for the WASP-12 system. There are a few other candidates, including WASP-4 b, but no conclusive result could be obtained for these systems as of yet. In this study, we present an analysis of new TESS data of WASP-4 b together with archival data, taking the light-time effect (LTE), induced by the second planetary companion, into account as well. We make use of three different Markov-Chain-Monte-Carlo models; a circular orbit with a constant orbital period, a circular orbit with a decaying orbit, and an elliptical orbit with apsidal precession. This analysis is repeated for four cases. The first case features no LTE correction, with the remaining three cases featuring three different timing correction approaches. Comparison of these models yields no conclusive answer to the cause of WASP-4b’s apparent transit timing variations. A broad range of values of the orbital decay and apsidal precession parameters are possible, depending on the LTE correction. This work highlights the importance of continued photometric and spectroscopic monitoring of hot Jupiters.


Introduction
The study of exoplanets has unveiled a great diversity in terms of their physical properties, orbital characteristics, and formation mechanisms.Among the most intriguing exoplanet classes are the hot Jupiters.These relatively rare companions (e.g.[1]) stand out as a group of gas giants that orbit their host stars at exceptionally close distances and have challenged our understanding of planetary formation and evolution.One of the key phenomena that has piqued the interest of astronomers is the tidal orbital decay of Hot Jupiters, where as of now only WASP-12 b could be observationally confirmed to experience this effect [2][3][4].However, there seems to be evidence for the occurrence rates of hot Jupiters in orbit around Sun-like stars to decrease with stellar age, as opposed to the occurrence rates of cold Jupiters, as found by Miyazaki and Masuda [5].This is supported by the likely observation of planetary engulfment by a Sun-like star made by the Zwicky Transient Facility [6].
Due to the close proximity of hot Jupiters to their stellar hosts, there are strong tidal interactions between the two bodies.These gravitational interactions can manifest themselves in the planet raising a bulge on the surface of the star.Depending on the rotational period of the star and the planetary orbital period, the viscosity of the stellar plasma can lead to a lag between the position of the tidal bulge and the virtual line connecting the stellar and planetary centers.If the star is rotating slower than the planet is orbiting around it, orbital angular momentum of the planet will be transferred to the star (equilibrium tide).The dynamical tide, arising from stellar oscillations, is also contributing to this [7].This means that the star will spin up and the planetary orbit shrink gradually [8,9].This provides us with insights into the long-term stability of these system.
WASP-4 b is a hot Jupiter discovered by Wilson et al. [10].It shows TTVs that have been examined before, and other effects, mimicking the orbital decay signature, were ruled out previously, like the Applegate mechanism [11,12].However, apsidal precession due to a small eccentricity of the planetary orbit, or TTVs arising from the Rømer effect (or light-time effect, LTE hereafter, Irwin 13), due to a companion in the system, have not yet been ruled out as the cause of the TTVs, even after the discovery of the companion candidate.Previous studies examining the decay rate in the case of tidal orbital decay, found values in the range from −2.4 ms yr −1 to −12.6 ms yr −1 [12,[14][15][16][17][18], with the most recent estimate being Ṗ = (−6.2± 1.2) ms yr −1 [19].
In this work, we re-examine the TTVs in this system by making use of recently acquired TESS photometry and combine them with archival data from previous works.In particular, we account for the time shift due to the LTE induced by the additional planet candidate, as discovered by Turner et al. [18].The observations are described in Sect.2, with our modelling and results thereof being described in Sect. 3 and 4, respectively.The latter are discussed in Sect. 5 and our final conclusions can be found in Sect.6.

Observations
For our analysis, we made use of the previously described data set in Harre et al. [19] for WASP-4 b.In short, this data set is mainly based on the homogeneous re-analysis of previously published light curves from Baluev et al. [16], with the addition of TESS Sectors 28 and 29, as well as a re-analysis of the TESS Sector 2 data.Furthermore, in [19], we re-fitted the publicly available light curves from the ExoClock project [20] and from WASP [10], and added eight transit observations taken with the CHEOPS space telescope [21].Included in this data set are also four occultation timings from the literature [22][23][24].All these timings, including three different corrections, can be found in Appendix B, Table A1 for the transit timings and Table A2 for the occultation timings.
To this data set, we add new TESS observations from Sector 69 at a cadence of 120 s.We make use of the Presearch Data Conditioning Simple Aperture Photometry (PDCSAP) flux for the analysis of the light curves, which consists of data produced by the TESS Science Processing Operations Center (SPOC) at NASA Ames Research Center [25].This data is publicly available at MAST 1 .

Light curve modelling
For the analysis of the TESS transits, we made use of the Transit and Light Curve Modeller (TLCM, Csizmadia 26) as described in Harre et al. [19].In a first run, we fitted all transits together to get the combined shape of all transits to reduce the impact of stellar activity on the transit timings.The respective priors are shown in Table 1.During a second run, we fixed the shapes of the transits to those of the combined model and fitted all transits individually, with only the transit ephemerides being free.In both cases, we used the median solution for the final results.

Transit timing variation analysis
For the analysis of the mid-transit times, obtained from our light curve modelling, we employ the same three models as Harre et al. [19].The first is a model assuming a circular Keplerian orbit, describing a linear ephemeris with a constant orbital period: where t tra (N) and t occ (N) are the calculated mid-transit and mid-occultation timings at the epoch N, T 0 denotes the reference mid-transit time, and P the planetary orbital period.This gives us two free parameters in the fit.
The second model describes the case of a decaying orbit due to the transfer of angular momentum (see e.g.Counselman 8,Rasio et al. 9).These are quadratic functions with a constant change in the orbital period of the planet: This constant change in the orbital period is denoted by the decay rate dP dN , which can be converted to the period derivative Ṗ = dP dt = 1 P dP dN .There are three free parameters (T 0 , P, and dP dN ) when fitting this model.The period derivative is linked to the stellar tidal modified quality factor Q ′ ⋆ via the constant-phase lag model, as defined in Goldreich and Soter [27]: where f denotes the tidal factor, M p and M ⋆ the planetary and stellar masses, R ⋆ the stellar radius, and a the semi-major axis of the planetary orbit.Depending on the ratio of the planetary orbital period to the stellar rotational period and the orbital (mis-)alignment, f takes different values.If the planetary orbital period is shorter than the rotation period of the star, as is the case for WASP-4 b [28], we get f = − 27 2 .Refer to Sect. 4 of Harre et al. [19] for a more detailed description, also including the case of inclined orbits and different planet-to-star period ratios.
The third model we are using is an apsidal precession orbit, where a small eccentricity e leads to the planetary orbit to precess around the star.This can, on relatively short timescales, induce the same TTV signature as orbital decay, which makes it hard to differentiate the two models if the orbital eccentricity is only loosely constrained.This sinusoidal model follows the descriptions of Giménez and Bastero [29]: t occ (N) = t 0 + P a 2 + N P s + e P a π cos ω(N), (7) where P s is the sidereal period, P a the anomalistic period, and ω the argument of pericenter.The sidereal and anomalistic periods are related via: with a constant change in the argument of pericenter dω dN .The relationship between ω and N can be described as follows: with w 0 being the argument of pericenter at the reference time T 0 .In total, this model gives us five free parameters in the fitting process.These three models are then fitted to the data via MCMC optimization using the emcee Python package [30].Per model, we use 100 walkers and use a burn-in period of 10.000 steps with a total chain length of 75.000 steps.Convergence is ensured by checking that the chains are longer than 50 times the autocorrelation time of the parameters.The priors for our models can be found in Table 2. U (0.5 ± 1.0) 0.22 ± 0.14

Light-time effect
The presence of a candidate planetary companion in a 7000 d orbit, as established by Turner et al. [18], would induce an orbital motion onto WASP-4.Due to the high mass of this candidate ("planet c" hereafter), the system's center of mass is shifted by about 8.9 times the radius of WASP-4, see Fig. 1.This will have a significant impact on the observed mid-transit times because of the time difference that it takes the light to travel from the far side to the near side of the orbit from our point of view.In some cases, this could induce TTVs akin to the imprint of tidal decay on short time scales, depending on the orbital period of planet c.To include a correction of the LTE, we apply the formula of Schneider [31] to find the maximum time-shift from this effect for circular orbits: where ∆T max is the maximum resulting time-shift due to the LTE, i is the inclination of the planetary orbit, and c is the speed of light.For WASP-4, using the planetary parameters of planet c, and assuming a circular orbit with i = 90 • , we obtain a maximum time-shift of ∆T max = 37.7 s.To validate this result, we simulate the system using the N-body code REBOUND [32], and measure the orbit of the star around the system's center of mass, as shown in Fig. 1.From this, we obtain ∆T max = 37.6 s.Using the non-circular solution from Turner et al. [18] for planet c, we get ∆T max, ell = 41.1 s.However, since the non-circular solution is not preferred in their paper, we adopt the circular solution and correct our transit times according to: where ∆T(t) is the time-shift due to the LTE at the time t, T 0,c is the time of inferior conjunction of planet c, and P c = 7001.0d is the orbital period of planet c.This formula arises because the star is at its furthest point from us at the time of inferior conjunction of planet c.Due to the uncertainty of the time of conjunction of planet c (T 0,c = 2455059 +2300 −2100 d [BJD TDB ]), we examine three cases of the LTE correction.Firstly, assuming the best-fit value of the time of conjunction, secondly, the minimum value of the 1σ interval, and lastly, the maximum value of the 1σ interval.

Transit fitting with TLCM
The priors for and parameters resulting from the combined TLCM fit of all TESS transits from Sector 69 to constrain the transit shape, are listed in Table 1, with the fit of the median model to the data shown in Fig. 2. The resulting mid-transit times can be found in Appendix B in Table A1.Table 2. Priors and results from our MCMC modelling of the transit timings without any LTE correction.T 0 is given in BJD TDB − 2450000.Units are given in square brackets if applicable.Priors were chosen according to the results of Harre et al. [19] with enough flexibility for the walkers to sufficiently explore the parameter space.

TTV fits without LTE correction
The priors for the MCMC modelling of the circular orbit, orbital decay and apsidal precession models and their resulting parameters are listed in Table 2.
We find an orbital decay rate of Ṗ = dP dt = (−5.75± 0.52) ms yr −1 , leading to a stellar modified quality factor for WASP-4 of Q ′ ⋆ = (6.10 ± 0.55) × 10 4 and a decay timescale of τ = (20.2± 1.8) Myr from the median solution.The BIC values of our solutions are obtained via: with k the number of free parameters and N the number of data points.For the orbital decay fit, we find ∆BIC decay = 119.4 between the linear and tidal decay models, in favour of the latter.In the case of apsidal precession, we find a precession rate of ω = dw dt = 1 P dw dN = (0.033 ± 0.007) • d −1 and a small eccentricity of e = 0.0013 ± 0.0005.The fit yields ∆BIC precession = 110.7 between the linear and apsidal precession models.This leads to a difference of ∆BIC dec.,prec.= 8.7 in favour of the orbital decay model.
The result of the MCMC fit to the transit and occultation timing data including the final median models is shown in Fig. 3. Especially the latest data points highlight the deviation from a linear ephemeris in this system.

TTV fits with LTE correction
In case the candidate planet c does exist, we corrected for the induced LTE with a maximum time-shift of about 38 s in three different cases.These cases are (1) planet c has the best-fit time of inferior conjunction from model #3 of Turner et al. [18], (2) the time of inferior conjunction is at the lower boundary of the respective 1σ confidence interval, and (3) the time of inferior conjunction is at the upper boundary of the 1σ confidence interval.The corrections, according to Eq. 11, are subtracted for each transit and occultation timing.
The results in Table 3 show that the orbital decay model provides the best fit to the data in case (1), with the nominal LTE correction, although this model is only slightly preferred over the apsidal precession model.In this case, the decay rate would be Ṗ = (−7.04± 0.52) ms yr −1 , leading to Q ′ ⋆ = (4.98 ± 0.37) × 10 4 and a decay time-scale of τ = (16.5 ± 1.2) Myr.The apsidal precession model yields ω = (0.024 ± 0.007) • d −1 with an eccentricity of e = 0.0029 ± 0.0019.This case is displayed in Fig. 4.   According to Turner et al. [18], the TTVs induced by planet c onto planet b should not exceed 2 s, which is why they are neglected here.

TTVs of WASP-4 b
As is evident from the results of our TTV modelling in Tables 2 and 3, the consideration of the LTE is very important in this system.The results of the orbital decay and apsidal precession modelling are only comparable to each other to a certain degree for the cases that were examined here.
If the detection of WASP-4 c should turn out spurious, then the measured decay rate value is comparable, albeit slightly smaller than those from previous studies [12,[14][15][16]18,19,33]. In this case, the apsidal precession model would be disfavoured by ∆BIC = 8.5, which is a bit higher than in the latest study [19].The decay rates are consistent within 1σ, with the eccentricity in the precession model being doubled, compared to the previous result, which is still in agreement with our value and uncertainty.The linear model is disfavoured by ∆BIC ≈ 110 − 120.
The three different LTE corrections lead to three vastly different results.The nominal LTE correction agrees within 3σ with the results without the LTE correction for most of the parameters.The other solutions do not.Depending on the correction, we get decay rates in the range from (−9.35 ± 0.52) × 10 4 ms yr −1 to (−1.45 ± 0.50) × 10 4 ms yr −1 .This means that the LTE induced by a planetary companion can solely explain the observed TTVs.The quite broad range of results can only be further constrained with more radial velocity measurements to reduce the size of the error bars on the parameters of WASP-4 c.However, a value near the middle of this range is more probable.From these results, we get modified stellar quality factors in the range from 3.8 × 10 4 to 2.4 × 10 5 , which includes the value of Q ′ ⋆ = 1.8 × 10 5 obtained for WASP-12 b [2].Moreover, this leads to decay timescales ranging from 12 Myr to 80 Myr.The apsidal precession models show a similar range of possible values like the orbital decay models.However, with more observations, either transits or especially occultations, the models should be relatively easy to distinguish.Occultation observations with e.g.JWST could perform a double duty in this case.They can help to refine atmospheric properties and provide very accurate occultation timing measurements, which might be able to rule out either of the TTV models.
These results highlight the importance of continued radial velocity observations, even, or especially, in hot Jupiter systems.One theory of how hot Jupiters get into their tight orbits is high-eccentricity migration, which necessitates the presence of a massive second companion for the excitation of the high eccentricities that are required for this migration pathway [34,35].Finding such a body requires long observational baselines due to the distance of the perturber to the host star.

Optimal observing strategy
For hot Jupiter systems in general, the most favourable tidal orbital decay systems have, according to Eq. 5 and Table 2 in Harre et al. [19], first and foremost with the highest impact, large stellar radii and short orbital separations.Small stellar tidal quality factors, high planetary masses and low stellar masses are also beneficial.However, stellar tidal quality factors can only be constrained after the observations.Furthermore, orbital periods of the companion bodies smaller than the stellar rotation periods are essential, but this only applies to equatorial orbits.For polar orbits, orbital decay should always be present, even for fast-rotating stars.Still, slow stellar rotation rates are beneficial here as well.In summary, as is known, hot Jupiters are prime targets to examine the effect of tidal orbital decay.Moreover, brown dwarfs on close-in orbits should even experience a more pronounced effect due to their higher masses.In addition, companions orbiting evolved or (sub-) giant stars should present the best laboratories to explore orbital decay.Yet, the expected lifetimes of these close-by companions will be relatively short, and observing them is only possible in just the right time window, providing us with only a very small sample of targets (see e.g.Grunblatt et al. 36).
In terms of radial velocity monitoring, after sufficient data has been accumulated for the characterisation of hot Jupiters, it would be optimal for the target to be re-observed every few months so that even farther out companions could be detected in principle.This would not only shed light onto the architecture of these systems, but also on their formation and migration pathways.This case study is an excellent example of this, where an outer companion candidate has been detected [18], but the whole detection hinges on the latest set of observations to not be spurious or have an incorrect offset in relation to the previously accumulated RV data.Due to the possibility for the outer planet to induce a small eccentricity onto the orbit of the inner planet, a non-detection constraining possible planet masses and distances from the host star could rule out the apsidal precession models in some cases.This could be achieved by comparing the eccentricity damping timescale to the orbital decay timescale.
For photometric observations, it would be ideal to obtain high-precision transit observations every few months.Long-term monitoring is essential to detect orbital decay signatures with a baseline of more than 15 years having been necessary for WASP-12 b [2].Occultation observations would be even more helpful to differentiate the orbital decay and apsidal precession models, not only from the point-of-view of TTV fitting, but also because they can indicate eccentricities directly.In addition to the timing, secondary eclipse observations can give clues about e.g. the atmopsheric composition or the thermal structure of the atmosphere and the presence of clouds (see e.g.The PLATO mission [42] has the potential to revolutionise the field with its long, uninterrupted observations.Due to the expected excellent precision of the observations, occultations of hot Jupiters are likely to be detected as well.This should allow discrimination between the orbital decay and apsidal precession models.Besides that, even relatively far-out companions could be detected with radial velocity follow-up from PLATOs ground segment, should they not transit the host star.As well as PLATO, ESA's Gaia mission [43], starting with its fourth data release, has the potential to improve our understanding of already known hot Jupiter systems.Its precise astrometric measurements will allow the detection of giant companions in these systems and provide detailed characterisation of their orbits.This will provide information on the formation and migration history of each system.

Conclusions
By analysing new TESS observations of WASP-4 in combination with archival timing data, and taking into account the additional planet candidate in the system, we get a broad range of results for the cause of the TTVs of WASP-4 b.We examined a total of four cases.The first case does not take into account the presence of planet c, and yields comparable results to the hitherto published results for this system in terms of the orbital decay parameters.In addition, the remaining cases consider, for the first time, the light-time effect from the host star's orbital motion around the system's center of mass, induced by planet c.Depending on the application of the timing correction, we find results that range from a slight preference (∆BIC ≈ 4) of the orbital decay model in the nominal case (using the nominal value for the time of inferior conjuction of planet c from Turner et al. [18]), over a strong preference (∆BIC ≈ 20) of the apsidal precession model using the lower boundary of the respective 1σ confidence interval, to a nearly indistinguishable result between the linear ephemeris and orbital decay models (∆BIC ≈ 2) using the upper boundary of the respective interval for the time of inferior conjunction of planet c.These results leave us with no conclusive answer to the question of what the true origin of the TTVs of WASP-4 b is.We need more radial velocity observations to better constrain the phase of planet c.Only then will we be able to determine whether the LTE solely explains the observed TTVs, or whether other mechanisms, like tidal decay or apsidal precession are present.It is, however, likely to be a mix of all of the effects mentioned here.
This case study highlights the importance of continued monitoring of hot Jupiter systems, in terms of photometric and radial velocity measurements.Only in this way, small effects, like orbital decay or apsidal precession, can be measured and differentiated.An additional benefit is given by the chance of discovering and characterizing companions to hot Jupiters, providing hints of the formation and migration scenarios that lead to these special systems.This is a necessary step to inform our theoretical models of planet formation and migration, and towards the final answer to the question of the origin of these planets.
Table A1.Transit timings of WASP-4 b including three different corrections for the LTE.Time is given in BJD TDB − 2450000."Time" denotes the mid-transit times without any light-time correction, "Time LTE" denotes LTE corrected timings using the nominal value for the time of inferior conjunction of planet c from Turner et al. [18], "Time LTE lower " denotes the timings corrected using the lower boundary of the respective 1σ interval, and "LTE upper " denotes the timings corrected using upper boundary of the interval.The "Source" column denotes the source of the timings, with 0-4 defined as the homogeneous reanalysis from Baluev et al. [16] (0), TESS timings (1) from [19] and this work (last sector), the ExoClock project (2) [20] as fitted by Harre et al. [19], as well as their CHEOPS timings (3), and WASP timings (4), respectively.

Figure 1 .
Figure1.The orbits of WASP-4 (star symbol, dashed line) and WASP-4 b (black dot, solid line) with reference to the system's center of mass (black plus symbol) from REBOUND.The x-and y-axes lie within the orbital plane of the planets, which are assumed to be coplanar.The orbit of planet c is not visible in this view, since it's semi-major axis is assumed to be 6.82 AU.

Figure 2 .
Figure 2. Phase-folded TESS light curve of WASP-4 b from Sector 69.The data are shown as black dots after modelling and correction with TLCM, with the median solution model shown as the red line.

Figure 3 .
Figure 3. O-C plot showing all transit timing (top) and occultation data (bottom) together with the orbital decay and apsidal precession models.The transit number from the reference epoch is shown on the x-axis, with the deviation from the median linear ephemeris is shown on the y-axis.Colors according to the legend, with the pink shaded area showing the 1σ interval around the median solution orbital decay fit.For the available TESS data, the weighted mean timings with their respective error bars for each sector have been added.

Figure 4 .
Figure 4. Same as Fig. 3, but with the nominal LTE correction applied to the timing data.

Table 1 .
The resulting parameters from the TLCM fit to the TESS S69 data.The impact parameters is given by b and u a and u b describe the quadratic limb-darkening parameters.T 0 is given in [BJD TDB − 2450000] and U denotes a uniform prior.

Table 3 .
Results from our MCMC modelling for the three cases of the LTE correction.The priors are the same as those given in Table2.T 0 is given in BJD TDB − 2450000.Units are given in square brackets if applicable.LTE refers to case (1), LTE low and LTE up refer to cases (2) and (3), respectively.∆BIC is defined relative to the circular Keplerian orbit case for each LTE correction.Note that the ∆BIC values should only be compared within a column, not between columns.