Tentative Evidence for Transit-timing Variations of WASP-161b

TESS takes exposures of 2 seconds duration but coadds the data into science products of three different cadences, due to data downlink limitations. The individual frames in the vicinity of sources of high interest are released with a cadence of 2 minutes, namely target pixel files (TPF). The 30 minutes cadence data for all the sources in the field of view are available as Full Frame Images (FFI; Ricker et al. 2015). 10 minutes cadence FFI is available for the extended mission data. WASP-161b has been observed in TESS Sector 7 (2019, January) with a 30 minutes cadence and Sector 34 (2021 January) with 2 minutes and 10 minutes cadence ⁵ .

We developed a photometric pipeline to obtain light curves from TESS images. The pipeline includes modules for astrometry correction, nearby source deblending, circular aperture photometry, and light curve detrending. The details of the pipelines are described in Yang et al. (2022). The light curves obtained from the pipeline are consistent with light curves generated by TESS Pre-search Data Conditioning (PDC; Smith et al. 2012) among the comparison sample of objects (Yang et al. 2022, 2021). In this work, the fitting results from our pipeline are consistent within 1σ to that derived from the PDC light curves for the 2021 data and quick look pipeline light curves (Twicken et al. 2016; Huang et al. 2020) for the 2019 data. As a result, all the parameters quoted in the paper are using our pipeline.

We detrend the TESS light curve using the published ephemeris (Barkaoui et al. 2019) as the baseline. We masked the data near the transit and fitted the data points with a linear function (Yang et al. 2022, 2021). High order polynomials (up to an order of 5) and cubic spline functions are also used to fit the trends in the data, but we find that the choice of function does not result in a significant difference to the light curve. We tested adding a trend (linear function of time) as a free parameter during transit modeling. The timing differences are within 0.2 minutes when compared to not including extra baseline parameters in light-curve modeling. Our derived TESS light curve is shown in Figure 1.

Zoom In Zoom Out Reset image size

Barkaoui et al. (2019) showed that the stellar variability does not appear to be significant enough to influence the transit modeling. The FFI data in 2019 shows the median brightness to be 9068 ± 7 e⁻ s⁻¹ while the 2021 10 minutes cadence data shows the brightness to be 9070 ± 11 e⁻ s⁻¹. We note that the brightness comparison between different epochs is performed among the same kind of data product for consistency. As a result, we do not find any significant evidence for variability in the TESS light curve either, i.e., any structure (other than transit) exceeding local median flux is less than 300 ppm. We thus conclude that the impact of stellar variability on the derivation of transit parameters is negligible.

SSO-Europa light curve is a high-quality ground-based light curve with an out-of-transit baseline both before and after transit which is fully described in Barkaoui et al. (2019). We detrend the data by masking the transit duration and fitting a linear function to the out-of-transit data points. We then fit the SSO-Europa photometry in a similar way to the process of fitting TESS light curves as described below. The fitted light curve is shown in Figure 2.

Zoom In Zoom Out Reset image size

The light curves we use are generated by our photometric pipeline, using 30 minutes cadence FFI for TESS Sector 7 and 2 minutes cadence TPF for TESS Sector 34. The detrended light curves are fitted with a transit model (Mandel & Agol 2002; Eastman et al. 2013) using a Monte Carlo Markov chain (MCMC) technique (Patil et al. 2010; Czesla et al. 2019). We choose a circular orbit for obtaining the transit midpoint, with the free parameters and priors the same as in previous work (Yang et al. 2022, 2021). The MCMC chain takes 100,000 steps with the first 50,000 steps as burn-in. We derive transit parameters by fitting the folded light curve with data from each Sector. Each transit is then separately fitted for the specific transit midpoint. The priors when fitting the individual transit are set as the derived parameters from the folded light-curve modeling. The derived transit timings are listed in Table 1. As shown in the Table, we also tried a Keplerian orbit with eccentricity as a free parameter.

The Keplerian orbit model has as free parameters, the radius ratio of the planet to the host star (R_p /R_*), the orbital inclination (i), the semimajor axis (in the unit of stellar radii; a/R_*), the time of periapse passage, the longitude of the ascending node, the argument of periapsis (ω), orbital eccentricity (e), and the limb darkening coefficients (LD; u1 for linear coefficient, u2 for quadratic coefficient). The parameterizations of e and ω are $\sqrt{e}$ sinω and $\sqrt{e}$ cosω, as the same of Eastman et al. (2013). Except for the LD priors, the free parameters take uniform priors. The LDs assume Gaussian priors, centered at the stellar model prediction with a σ of 0.05 (Claret 2018; Yang et al. 2021). Applying the stellar parameters from Barkaoui et al. (2019), the priors of u1 and u2 are 0.31 ± 0.05 and 0.22 ± 0.05. The same steps are executed for circular orbit fitting. The 2 minutes cadence TESS data provides the strongest constraints on LD parameters since it samples the transit ingress and egress the best. It yields a fitted e of 0.28 ± 0.14 (as shown in Table 1).

The fitting results and uncertainty are checked by simulating synthetic light curves from the data where the data points are randomly scattered by their Gaussian uncertainty (the simulation details are similar to Yang et al. 2022, 2021). We generate 1000 such simulated light curves and assess the distribution of fit parameters with them. The simulation indicates that most of the transit parameters have consistent parameter estimates between MCMC obtained values and the median results from the simulations. However, uncertainties of e and a/R_* are underestimated by MCMC fitting by a factor of ∼2. MCMC fitting inevitably has difficulties in global parameter estimates when there are too many free parameters and some of which are not tightly constrained by the data (Mackay 2003; Hogg & Foreman-Mackey 2018; Yang et al. 2021). We cite the uncertainties of e and a/R_* as the standard deviation of our simulation fit results, rather than the MCMC derived values, in this work.

We adopt double the timing uncertainty derived from the folded light curve, for the timing uncertainty of each epoch. We also provide timing uncertainties obtained from light curves of single visits (Table 1). For single visit light curve, we find that timing uncertainty inferred directly from MCMC fitting can be overestimated. The time series of the parameters in the chain present clumps when fitting single transits, indicating imperfect MCMC fits. Timing uncertainties derived from single-transit light curves are ∼four times (for 30 minutes data) and three times larger (for 2 minutes data) than the result from folded light curves among one TESS Sector. According to error propagation, the uncertainty from the folded light curve should be 1/2 of the single visit given that every TESS sector contains four successful transits. Moreover, citing timing uncertainties from MCMC fitting to a single visit would lead to a significantly underestimated reduced χ² of 0.13 when modeling timing evolution.

We note that the correlated noise is not included in our light-curve modeling though a proper correlated noise model can potentially reduce the systematic errors (Maxted et al. 2021; Patel & Espinoza 2022). The Gaussian process is commonly utilized when modeling the correlated noise (Foreman-Mackey et al. 2017). TESS data is stable enough that the noise dependency toward time (within 1 day) can be ignored in our precision of transit fitting (Guerrero et al. 2021; Ivshina & Winn 2022). TESS PDC light curves are modeled with correlated noise. We have compared our generated light curves and PDC light curves in previous work and find the difference is negligible in transit analyzing (Yang & Wei 2021; Yang et al. 2021, 2022). The timings of WASP-161b obtained by our generated light curve and PDC light curve are consistent within 0.1 minutes and the uncertainty from our generated light curve is 1.4 times larger.

Some orbital parameters that are free parameters in light-curve fitting are better constrained by the RV data points. It, therefore, makes sense, if additional constraints on parameters can be obtained by undertaking joint fits to the RV and light-curve data. We constrain the value of e by a combined fit to TESS 2 minutes cadence light curve and the RV curve from Barkaoui et al. (2019). The combined fit connects the constraints of both the LC and RV curves by adding the log-like likelihood of RV and LC models together (Czesla et al. 2019). The likelihood weights between the RV model and LC model are set equal. The combined model takes the same MCMC technique. The free parameter of the combined model includes all the free parameters of LC fitting and additional parameters for the RV curve, i.e., offset in radial velocity (γ), and the RV semiamplitude (K). The fit gives an e as 0.34 ± 0.03. The fitting results of key parameters are as shown in Table 1.

Moreover, the MCMC fit to the RV curve alone using the Keplerian orbit model yields e = 0.34 ± 0.03 indicating that the RV data dominates the value of e in the combined fit. Our derived e is consistent with the 3σ upper limit of 0.43 reported by Barkaoui et al. (2019).

We also attempt to fit a circular orbit model to the RV curve and find the fit is not as good as the eccentric Keplerian orbit model. The reduced χ² of the circular model is ∼3, compared to the reduced χ² of the Keplerian model ∼1. Applying the Bayesian information criterion (BIC; details in Kass & Raftery 1995), the ΔBIC between the circular orbit and Keplerian orbit models is 27, preferring a Keplerian orbit model. The fitted light and radial velocity curves are present in Figure 1.

In addition, we jointly fit the light curves and radial velocities under the assumption of a circular orbit. We find that the derived timings from the joint fit are consistent (within 0.4σ) with timings from fitting the light curve only with a circular orbit. We note that the reduced χ² of the RV curve is ∼3, and is indifferent to whether we do a joint fit or fit separately.

We conclude that the combined LC and RV fit with free eccentricity is the optimal model for the orbital parameters. A circular orbit is rejected by the RV curve. Fitting the LC and RV separately leads to inconsistency in the derived eccentricity (as shown in Table 1).

The initial ephemeris is taken from Barkaoui et al. (2019) which gives a transit midtime (T_c ) at 2457416.5289 ± 0.0011 days (Heliocentric Julian Date; HJD) and a period of 5.4060425 ± 0.0000048 days. Table 1 shows the transit midpoint in HJD, derived from the fits to the TESS and SSO-Europa light curves described earlier. The difference in minutes of the transit midpoint, with respect to the original ephemeris, is also shown in Table 1. The transit midpoint in SSO-Europa light curve is consistent at the 1.6σ level with that in Barkaoui et al. (2019). The uncertainty is estimated as the quadrature sum of the ephemeris uncertainty from Barkaoui et al. (2019) and timing uncertainty from SSO-Europa light curve. We note the previously reported ephemeris is indeed the best result given the uncertainties in the archival observations. However, the TESS data taken at later epochs, show that the difference in the transit midpoint becomes increasingly larger. We note that the timings from Barkaoui et al. (2019) are in HJD and TESS timings are in Barycentric Julian Date. The difference is within ±4 s which is negligible with the timing precision in our discussion.

We assess the accuracy of our TESS timing by applying our light-curve generation and fitting analysis to WASP-58b and HAT-P-31b. The TESS timings are consistent within 1σ with respect to the published ephemeris (Mallonn et al. 2019). The targets for TESS timing check quoted here are selected randomly. In a separate work we discuss other comparisons of TESS timing precision (Shan et al. 2021). We also note that TESS timing has been used for reliable long-term monitoring of WASP-12b (Turner et al. 2021) and to detect TTVs in WASP-4b and XO-3b (Bouma et al. 2019; Yang & Wei 2021).

We find that if the TESS timings are accurate, the transit midpoint of SSO-Europa light curve is earlier by 6.97 minutes relative to the constant period assumption (as shown in Figures 2 and 3). This is at a significance level of 4.63σ calculated as the sum in quadrature of the timing uncertainty from SSO-Europa light curve and TESS ephemeris uncertainty (as shown in Figure 3). The TESS ephemeris uncertainty is evaluated as the 1σ region of the MCMC chain from a linear fit to the TESS timings.

Zoom In Zoom Out Reset image size

Conversely, the timing from Barkaoui et al. (2019) yields a 132.30 minutes deviation at 53.03σ when compared to TESS timings. TESS observations indicate transit midpoints which are earlier by ∼67.24 ± 2.22 minutes at 2019 January, and ∼202.71 ± 0.91 minutes at 2021 January, compared to the ephemeris from Barkaoui et al. (2019). The values are the mean values among the four transits in each period of TESS observation, each of which are shown in Table 1. The uncertainty is the mean of the quadrature sum of the timing uncertainties.

We can not distinguish whether the timing from SSO-Europa light curve or the timing from Barkaoui et al. (2019) is more accurate because of the constant period assumption in their work and the absence of the timings from the individual epochs in Barkaoui et al. (2019). We note that their timings are consistent with our analysis of the timing from SSO-Europa light curve, under the assumption of a constant period, which, however, seems unlikely.

Furthermore, we can conclude that the transit midpoint in the SSO-Europa data is statistically different compared to the transit midpoint in the TESS data (Figure 3) and similarly, that the TESS data show transit midpoints which are different from the published ephemeris of WASP-161b by large amounts.

We fit the evolution in transit midpoints with two scenarios: a period evolution $\dot{P}$, which is zero and a $\dot{P}$ which is constant. We refer to these as the linear model and the quadratic model, respectively. We find that a linear model is not as good as a quadratic model when explaining the transit midpoints from both the SSO-Europa and TESS observations (as shown in Figure 4). The BIC of the quadratic model is 10.33 smaller than the linear model based on our analysis of the archival data. Furthermore, if we use the ephemeris from the Barkaoui et al. (2019), we find that the BIC of the quadratic model is 6,500 smaller than any linear model. The quality of fits are shown in Table 2.

Zoom In Zoom Out Reset image size

Table 2. Models for Timing Deviation

Linear Model: td = al × t + bl
Period Decaying Model: td = ad × t2 + bd × t + cd
Precession Model: td = ap 1 × cos(ωt0 + ωt × t) + bp × t + cp
Parameters	Only TESS Timing	SSO-Europa + TESS Timing	Barkaoui et al. (2019) reported + TESS Timing
Linear Model
al	−0.1836 ± 1.5 × 10−6	−0.1796 ± 3.1 × 10−6	−0.1260 ± 0.0003
bl	132.3 ± 3.5	125.7 ± 5.9	43.1 ± 550.2
reduced χ2	⋯	2.03	940.93
BIC	⋯	18.67	6590.90
Period Decaying Model
ad	⋯	−1.54 × 10−5 ± 3.0 × 10−6	−6.67 × 10−5 ± 1.3 × 10−6
bd	⋯	−0.1387 ± 0.0080	0.0106 ± 0.0027
cd	⋯	101.7 ± 4.8	0.0 ± 1.4
$\dot{P}$ (days per day)	⋯	−1.16 × 10−7 ± 2.2 × 10−8	−5.01 × 10−7 ± 9 × 10−9
reduced χ2	⋯	0.29	0.33
BIC	⋯	8.34	8.58
Precession Model
ap	⋯	−842.49 ± 0.57	−842.50 ± 0.57
ωt0	⋯	1.689 ± 0.008	2.114 ± 0.043
ωt	⋯	3.22 × 10−4 ± 1.6 × 10−5	4.42 × 10−4 ± 1.0 × 10−5
bp	⋯	−0.409 ± 0.008	−0.323 ± 0.014
cp	⋯	0.12 ± 5.78	−435.50 ± 31.79
$\tfrac{d\omega }{{dN}}$ (rad)	⋯	1.74 × 10−3 ± 0.08 × 10−3	2.39 × 10−3 ± 0.06 × 10−3
reduced χ2	⋯	0.36	0.46
BIC	⋯	10.61	11.10

Note. (1) a_p has a uniform prior of (−843.5, −841.5), according to Equation (2). t_d is the timing difference relative to the zero line present in Figures 3 and 4. t is the transit time relative to the transit zero-point. a, b, c, and $\omega$ with subscripts are all coefficients.

Download table as:  ASCIITypeset image

The origin of orbital period decay in planet-satellite systems and hot Jupiters have been discussed in the literature (Goldreich & Soter 1966; Murray & Dermott 1999; Ragozzine & Wolf 2009; Correia & Laskar 2010; Hansen 2010, 2012; Yee et al. 2020, and references therein). Gravitational tides circularize elliptical orbits on timescales of ∼1 Myr unless the eccentricity is excited by the gravitational influence of an outer planetary companion (Mardling 2007). For tides to cause a decay in the orbital period of the planet, the star must be rotating slower than the orbital velocity of the planet (Hut 1981; Murray & Dermott 1999; Leconte et al. 2010; Millholland & Laughlin 2018). In the case of WASP-161b, Barkaoui et al. (2019) has measured the stellar rotation to be faster than the orbital velocity, with a rotation period of ∼4.8 days compared to an orbital period of ∼5.4 days. Thus, it seems unlikely that tidal dissipation can explain the TTVs unless the planet is in a retrograde orbit.

Alternately, the rotational and tidal bulges on the planet can introduce apsidal precession; this is, however, thought to affect the light curve shape to a greater degree than transit timings (Pál & Kocsis 2008; Ragozzine & Wolf 2009). We do not detect any clear evidence for variation in the transit duration.

We plot our transit midpoint variations and compare them to those expected from the period decay and precession models in Figure 4.

The period decay model yield transit times (t_tra)for the Nth transit relative to the transit zero-point t₀ of Patra et al. (2017); Yee et al. (2020):

$$\begin{eqnarray}&&{t}_{\mathrm{tra}}(N)={t}_{0}+{NP}+\displaystyle \frac{1}{2}\displaystyle \frac{{dP}}{{dN}}{N}^{2}\end{eqnarray} \tag{ 1 }$$

where t_tra(N) is transit timing at Nth transit relative to the timing zero-point t₀, P is the period at t₀. t_d is the difference in transit time relative to that expected from a constant orbital period i.e., t_d = t_tra(N) − t₀ − NP. The timing variance is thus a quadratic function (as shown in Table 2), and leads to a $\dot{P}$ of −1.16 × 10⁻⁷ ± 2.25 × 10⁻⁸ days per day (or −3.65 s yr⁻¹). The P/$\dot{P}$ is 0.128 Myr. Given the derived eccentricity of WASP-161b, and the expected circularization of the orbit on relatively short timescales, a companion in an outer orbit is needed to excite and maintain the high eccentricity (Naoz et al. 2011; Dawson & Johnson 2018). Yu et al. (2021) present a nonlinear interaction, very effective in tidal dissipation, enabling a dissipation timescale as short as 10⁴ yr.

For apsidal precession, the transit timing is predicted as (Patra et al. 2017; Yee et al. 2020):

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{tra}}(N) & = & {t}_{0}+{{NP}}_{s}-\displaystyle \frac{{{eP}}_{a}}{\pi }\cos \omega (N)\\ \omega (N) & = & {\omega }_{0}+\displaystyle \frac{d\omega }{{dN}}N\\ {P}_{s} & = & {P}_{a}\left(1-\displaystyle \frac{1}{2\pi }\displaystyle \frac{d\omega }{{dN}}\right),\end{array}\end{eqnarray} \tag{ 2 }$$

where P_s is the sidereal period, P_a is the anomalistic period and ω(N) is the precession amplitude of the Nth orbit. As shown by Ragozzine & Wolf (2009), the precession amplitude is proportional to the Love numbers of the planet and the star.

The observational precession rate is derived from TTVs fitting. The TTVs are fitted with a combined-formula including a cosinusoidal and a linear formula (as shown in Table 2). We use the best-fitting TTVs from our generated SSO-Europa and TESS light curves. However, we find that the Barkaoui et al. (2019) reported values and TESS timing give a similar result (within 30%). We obtain a precession rate of 3.22 × 10⁻⁴ ± 1.6 × 10⁻⁵ rad per day (or ∼6° yr⁻¹), corresponding to a cycle of ${53}_{-2.5}^{+2.7}$ yr.

Both models provide reasonable fits to the observational TTVs, though the precession model has more free parameters (as shown in Table 2) and results in abnormal Love Numbers of 141 and 123 for the planet and star, respectively. The tidal model clearly also has issues given the rotation rate of the star. Further data is clearly required to understand the unusual orbital parameters of WASP-161b.

We have presented the evidence for TTVs of WASP-161b and discussed plausible interpretations of the measurements. It is interesting to note that of the four observed planets (WASP-161b, WASP-12b, WASP-4b, XO-3b) with TTVs, all have earlier timing than the constant period model. Apsidal precession should result in an unsigned timing variation which is not the case. Thus, the declining period model with the presence of an undetected outer planet to explain the significant eccentricity, is the preferred explanation for the TTVs in WASP-161b. We note that the magnitude of the transit time difference in SSO-Europa light curve is well matched to the difference between HJD and JD on the date of the SSO-Europa observations; although Barkaoui et al. (2019) states SSO-Europa light curve is in HJD in their data release notes. ^{6 7}

Further observations with more precise timing are needed to distinguish between the different scenarios presented. We have been approved for observations with the space telescope CHEOPS (Benz et al. 2021) for two transits around 2022 January.

The transit midpoints predicted from the different models can be distinguished with the CHEOPS observation (Benz et al.2021), given the quoted precision of CHEOPS photometry and sampling rate (∼700 ppm precision with an exposure time of 30 s). The CHEOPS observation would present a TTV larger than 200 minutes when compared to the ephemeris from Barkaoui et al. (2019) if the period decay model is as derived from the Barkaoui et al. (2019) and TESS timings (Figure 4). In the favored scenario, the period decay model obtained from SSO-Europa and TESS timings would yield a TTV of ∼5 minutes when compared to the constant period timing predicted from TESS observations. The difference between a period decay model fit to the Barkaoui et al. (2019) ephemeris and TESS data would be 22 minutes relative to the model fit to the SSO-Europa and TESS data. Thus, we conclude that forthcoming observations will be able to enhance our understanding of the orbital evolution of WASP-161b.