Transit observations and analysis
We analysed a heterogeneous dataset of light curves from space- and ground-based telescopes (Supplementary Table 1) to measure transit times for the ultimate purpose of modelling TTVs. We used PyMC344, exoplanet (https://docs.exoplanet.codes/en/stable/)45 and starry46 to fit the light curve, incorporating tailored models for correlated noise and instrumental systematics appropriate for each dataset.
Our analysis of the K2 and Transiting Exoplanet Survey Satellite (TESS) light curves involved two distinct approaches with different noise models. For the joint analysis of all transits in both light curves (described below), we modelled stellar variability as a Gaussian process47. By contrast, for measuring individual transit times (see below), a third-order basis spline was sufficient to model the local correlated noise.
To account for systematics in the Spitzer data, we used pixel-level decorrelation (PLD)48, which uses a linear model with a design matrix formed by the PLD basis vectors (see below for more details). For the ground-based datasets, we included a linear model with a design matrix formed by airmass, pixel centroids, and the pixel response function peak and width covariates, when available.
The limb-darkening coefficients were calculated using stellar parameters from ref. 2 by interpolation of the parameters tabulated by refs. 49,50. These were fixed for individual transit fits but sampled with uninformative priors in the joint K2 and TESS analysis described below.
We used Broyden–Fletcher–Goldfarb–Shanno optimization51 as implemented in scipy.optimize for initial parameter estimates, followed by posterior sampling with the No-U-Turn Sampler52, an efficient gradient-based Hamiltonian Monte Carlo sampler implemented in PyMC3. The chains were well mixed (Gelman–Rubin statistic less than about 1.01) with negligible sampling error.
We first performed a joint fit of K2 and TESS data assuming a linear ephemeris (see below). We then measured all individual transit times uniformly using Gaussian priors from the joint fit for R p /R ⋆ , b (the transit impact parameter), and T 14 (the total transit duration), and uniform priors for T c (the transit centre time) centred on predicted times. We verified that T c posteriors were Gaussian and isolated well from prior edges.
In all individual transit fits, we assumed Gaussian independent and identically distributed noise and included a jitter parameter σ jit to account for underestimated photometric uncertainties. The log-likelihood was, thus,
$$\text{ln}{\mathcal{L}}=-\frac{1}{2}\text{ln}| \varSigma | -\frac{1}{2}{{\bf{r}}}^{{\rm{T}}}{\varSigma }^{-1}{\bf{r}}+{\rm{const.}},$$
where Σ is the diagonal covariance matrix with entries equal to the total variance (that is, the ith entry is \({\sigma }_{{\rm{tot,}}i}^{2}={\sigma }_{{\rm{obs,}}i}^{2}+{\sigma }_{{\rm{jit}}}^{2}\), where σ obs,i is the observational uncertainty of the ith data point), and r is the residual vector (\({\bf{r}}=[{\widehat{y}}_{1}-{y}_{1},{\widehat{y}}_{2}-{y}_{2},\ldots ,{\widehat{y}}_{n}-{y}_{n}]\), where \(\widehat{y}\) is the model and y is the data consisting of n measurements). When a given transit event was observed by several telescopes (for example, at Las Cumbres Observatory (LCO))53, or several band-passes from the same instrument (for example, from MuSCAT3)54, we jointly fitted all light curves covering the same event.
... continue reading