A search for multi-planet systems with TESS using a Bayesian N-body retrieval and machine learning

TL;DR

This paper introduces a Bayesian N-body and machine learning framework to detect multi-planet systems with TESS data by analyzing transit timing variations and vetting signals efficiently.

Contribution

It presents a novel combination of Bayesian N-body simulations and machine learning to identify multi-planet candidates and characterize perturbing bodies in TESS data.

Findings

Evidence for 3 new multi-planet candidates with non-transiting companions.

A machine learning model predicts perturbing planet parameters from TTVs.

The method expedites detection of multi-planet systems in large TESS datasets.

Abstract

Transiting exoplanets in multi-planet systems exhibit non-Keplerian orbits as a result of the gravitational influence from companions which can cause the times and durations of transits to vary. The amplitude and periodicity of the transit time variations (TTV) are characteristic of the perturbing planet's mass and orbit. The objects of interest (TOI) from the Transiting Exoplanet Survey Satellite (TESS) are analyzed in a uniform way to search for TTVs with sectors 1-3 of data. Due to the volume of targets in the TESS candidate list, artificial intelligence is used to expedite the search for planets by vetting non-transit signals prior to characterizing the light curve time series. The residuals of fitting a linear orbit ephemeris are used to search for transit timing variations. The significance of a perturbing planet is assessed by comparing the Bayesian evidence between a linear and non-linear ephemeris, which is based on an N-body simulation. Nested sampling is used to derive posterior distributions for the N-body ephemeris and in order to expedite convergence custom priors are designed using machine learning. A dual input, multi-output convolutional neural network is designed to predict the parameters of a perturbing body given the known parameters and measured perturbation (O-C). There is evidence for 3 new multi-planet candidates (WASP-18, WASP-126, TOI-193) with non-transiting companions using the 2 minute cadence observations from TESS. This approach can be used to identify multi-planet systems and stars in need of longer radial velocity and photometric follow-up than those already performed.

Tables

Table 1. Table 1: TTV Estimates for the Solar System

Planet	Mplanet	Period	Eccentricity	TTVmax
	(Earth)	(day)		(min)
Mercury	0.055	87.97	0.205	1.3
Venus	0.815	224.70	0.007	4.5
Earth	1	365.24	0.017	6.5
Mars	0.107	686.98	0.094	42
Jupiter	317.83	4332.59	0.049	1467
Saturn	95.16	10759.22	0.056	5652
Uranus	14.54	30588.74	0.046	11055
Neptune	17.15	59799.90	0.0113	6538

Table 2. Table 2: Stellar Parameters

TIC ID	T mag	Rs [Sun]	Teff [K]	log(g)	[Fe/H]	u1	u2
25155310	10.56	1.27	5800	4.28	-0.06	0.43	0.14
100100827	8.83	1.23	6400	4.37	0.1	0.39	0.15
183985250	9.10	0.99	5422	4.42	0.0	0.47	0.12
…			…				…

Table 3. Table 3: Results of AI Vetting Algorithm

TIC ID	Transit Probability
120916706	1
122612091	1
159951311	1
164767175	0.0068
166836920	1
219253008	1
220479565	0.0016
…	…

Table 4. Table 4: Phase folded light curve parameters

	WASP-18 b	TOI 193	WASP-126 b
TIC ID	100100827	183985250	25155310
Period [day]	0.9415 $\pm$ 1.34e-06	0.7921 $\pm$ 1.80e-05	3.2888 $\pm$ 1.94e-05
Tmid	1354.4579 $\pm$ 4.15e-05	1354.2151 $\pm$ 3.56e-04	1327.5207 $\pm$ 2.57e-04
Rp/Rs	0.0964 $\pm$ 1.08e-04	0.0451 $\pm$ 5.67e-04	0.0780 $\pm$ 2.17e-04
a/Rs	3.570 $\pm$ 1.55e-02	3.671 $\pm$ 1.27e-01	7.887 $\pm$ 3.98e-02
Inclination	85.68 $\pm$ 2.79e-01	75.40 $\pm$ 6.12e-01	89.51 $\pm$ 4.35e-01
${\sigma }_{res}$ [ppm]	565	611	1302
SNR	19.16	2.65	5.59
TTVmax (min)	0.4	3.0	1.5
Transits	47	32	21

Table 5. Table 5: Ephemeris Parameters

	WASP-18 b	WASP-126 b	TOI 193	TOI 193
TIC ID	100100827	25155310	183985250	183985250
M∗ [M.]	1.22	1.12	1	1
M1 [MEarth]	3314.82	90.3	22.0	22.0
Linear Evidence	-45.0	-18.9	-24.4	-24.4
L. Period [day]	0.941452 $\pm$ 1.34e-06	3.288782 $\pm$ 1.94e-05	0.792043 $\pm$ 1.82e-05	0.792043 $\pm$ 1.82e-05
L. Tmid	1354.457943 $\pm$ 4.15e-05	1327.520692 $\pm$ 2.57e-04	1354.215278 $\pm$ 3.54e-04	1354.215278 $\pm$ 3.54e-04
Non-linear Evidence	-34.1	-15.9	-21.1	-21.5
NL. P1 [day]	0.941240 $\pm$ 2.34e-06	3.28883 $\pm$ 3.71e-05	0.792074 $\pm$ 2.39e-05	0.792169 $\pm$ 3.44e-05
NL. Tmid	1354.458393 $\pm$ 5.45e-05	1327.5203 $\pm$ 3.85e04	1354.216555$\pm$ 6.60e-04	1354.215943$\pm$ 8.09e-04
NL. Mass 2 [MEarth]	55.2$\pm$12.3	64.2$\pm$24.6	39.4$\pm$9.5	100
NL. Period 2 [day]	2.1558$\pm$0.0064	7.63$\pm$0.17	1.516$\pm$0.021	1.65$\pm$0.04
NL. Eccentricity 2	0.015	0.05	0.02	0.1
NL. Omega 2 [rad]	4.74$\pm$0.24	3.197$\pm$0.59	3.19$\pm$0.51	1.92$\pm$0.78

Equations

$$F_{obs}=\left(a_0+a_1\ast t+a_2\ast t^2\right)F_{transit}.$$

$$t_{next}=n_{orbit}\ast P+T_{mid}$$

$$t_{1,next}=f_{Nbody}(n_{orbit},M_{\ast },M_1,P_1,e_1,{\Omega }_1,P_2,M_2,e_2,{\Omega }_2)+T_{1,mid}$$

Full text

methods: data analysis – planets and satellites: detection – techniques: photometric

1 Introduction

The Transiting Exoplanet Survey Satellite (TESS) is conducting an all-sky photometric survey to discover hundreds of transiting planets around bright stars that are most suitable for mass measurements through radial velocity (RV) observations (Ricker et al. 2015; Ricker et al. 2016). TESS acquires observations on a 30-minute cadence of all objects in the field of view but increases the cadence to 2-minutes for select bright stars with the goal of detecting small transiting planets (Stassun et al., 2018). Each sector of TESS data is observed for 27.4 days (two spacecraft orbits), allowing the survey to cover most of the sky in 2 years. Simulations have shown that TESS is capable of detecting hundreds of small planets (Sullivan et al. 2015; Bouma et al. 2017; Barclay et al. 2018;). To date, TESS has discovered only a few Earth-sized planets (Huang et al. 2018; Gandolfi et al. 2018; Dragomir et al. 2018; Vanderspek et al. 2019) and a some larger ones too (Nielsen et al. 2018; Trifonov et al. 2019; Wang et al. 2019; Huber et al. 2019; Rodriguez et al. 2019). Other discoveries have found new multiplanet systems (Leleu et al. 2019; Quinn et al. 2019) with only two discoveries yielding measurable transit timing variations so far (Kipping et al. 2019; Bouma et al. 2019). Hadden et al. 2018 suggests within the first 2 years of the TESS mission there will be of order 1-10 planet masses precisely measurable with transit timing variations. This dreary estimate contradicts the $\sim$5$\%$ estimate in Ballard (2019) and could be a result of the large constraint on signal-to-noise (7.3) or imperfect knowledge of the underlying distribution of multi-planet system architectures in Hadden et al. 2018.

Transiting exoplanets in multi-planet systems exhibit non-Keplerian orbits as a result of the gravitational influence from companions which can cause the times and durations of transits to vary. Transit timing variations (TTVs) are characteristic of the perturbing planet’s mass, eccentricity and period and can be used to constrain the planetary parameters of the perturbing body (Agol et al. 2005; Holman & Murray 2005; Nesvorný & Morbidelli 2008; Agol & Fabrycky 2018). Understanding planet-planet perturbations aided in the historical discovery of Neptune by examining the acceleration, then deceleration of Uranus as it passed Neptune (Adams 1847; Le Verrier 1877). Perhaps it is a coincidence but Uranus also exhibits the largest TTV signal of all the Solar System planets (see Table 1).

Past detections of transit timing variations have come mainly from the Kepler mission due to its long term photometric monitoring of the same field of stars for $\sim$4 years (Borucki et al., 2010). Timing variations enable mass and eccentricity measurements of small planets that would otherwise be difficult to measure with radial velocity observations (e.g., Jontof-Hutter et al. 2016; Hadden & Lithwick 2016). In some cases, the timing variations can be used to assess the presence of additional bodies that may not transit yet still allow for constraints on their mass. Since one of the primary objectives for TESS is to measure masses for 50 planets smaller than 4 Rearth timing variations seem like an intriguing pursuit. However, one large constraint on the search for timing variations is TESS’s significantly shorter baseline compared to Kepler. Only select areas of the sky111https://heasarc.gsfc.nasa.gov/docs/tess/primary-science.html will have data for more than 27.4 days and based on prior studies $\sim$20 transits are generally needed to ensure a unique solution (Nesvorný & Morbidelli, 2008). This limits the search to planets in compact orbital resonances if they have only 1 sector of data.

A uniform analysis is conducted on the TESS objects of interest catalog in search of transit timing variations. The two minute cadence data from TESS will help determine prime candidates for RV follow-up that would benefit from mass and orbit constraints. In this paper the observations are described in section 2, followed by the data analysis and transit vetting algorithm in section 3. Section 4 explains the analysis of the light curves and derivation of the linear and non-linear ephemeris models. The results are presented in section 5 and the paper concludes in section 6.

2 Observations

The results presented in this paper use data from Sectors 1, 2 and 3 of the TESS spacecraft. Light curves of 74 targets from the TOI catalog (Stassun et al., 2018) have their time series photometric measurements analyzed. The targets are observed with a 2-minute cadence using a 11$\times$11 pixel subarray centered on the target. The photometric data were processed using the Science Processing Operations Center (SPOC) pipeline (Jenkins et al., 2016), which is based on the predecessor Kepler mission pipeline (Jenkins et al. 2010; Jenkins 2017). Each target is analyzed in the same manner unless otherwise stated. The flux is extracted from the target pixel using apertures centered on the star extending 2–3 pixels in radius. A background subtraction is performed by computing the average background flux from the lower 50th percentile of pixels within an annulus around the target (see Figure 1). The target pixel files include quality flags that indicate when photometric measurements may have been compromised by non-optimal operating conditions on the spacecraft. All flagged data points are removed from the light curves. Additionally, 3 $\sigma$ outliers from a rolling median filter with a bin size of 32 minutes are removed.

2.1 Aperture Photometry

Two different apertures are used to extract the stellar flux and are based on the aperture from the SPOC pipeline. The SPOC aperture is given in the target pixel files and another aperture is a dilated version of the original SPOC aperture. The larger aperture is dilated by one pixel which is used to incorporate more flux from the star and has been found to reduce pointing induced systematics for some stars (see Figure 2). The aperture which produces the lowest scatter in the residuals of the phase folded light curve fit is used to model the transit parameters. Figure 2 shows an example of the aperture selection for the target Pi Mensae b because it has been characterized by two different groups (Huang et al. 2018; Gandolfi et al. 2018;) The original SPOC aperture for Pi Mensae b is a little small and because of that it is sensitive to sub-pixel shifts of the star on the detector. The flux for Pi Men b can vary up to a few hundred ppm with just a 0.1 pixel shift in the centroid. Increasing the aperture size decreased the scatter in the residuals by $\sim$30$\%$ (see Figure 2).

3 Data Analysis

A custom pipeline is designed to analyze the time series measurements with artificial intelligence being used to vet data that does not exhibit a transit-like shape. A vetting algorithm ensures each target is processed in a homogeneous manner which reduces the removal of arbitrary choices made by human vetting. Vetting targets saves computational time because full parameter posteriors are derived for each light curve fit using nested sampled ($\sim$10-30 seconds per light curve fit to achieve $\sim$10000 evaluations of the parameter space). Since TESS will observe a plethora of transit events from 1000s of TOIs, some of which will have low signal, it’s crucial to dedicate the time to model each light curve with a global optimization scheme (e.g. nested sampling and not Levenberg-Marquardt) in order to produce robust uncertainties from posteriors. This study focuses only on planets from the TOI catalog since they are prioritized for follow-up observations. Additionally, all existing multi-planet systems are ignored and will be analyzed in the future and contributed to community driven follow-up.

3.1 Transit Vetting with Artificial Intelligence

The ideal algorithm for detecting planets should be fast, robust to noise and capable of learning and abstracting highly non-linear systems. A neural network (or deep net) trained to recognize planets with photometric data provides the ideal platform. Deep nets are composed of layers of “neurons”, each of which are associated with different weights to indicate the importance of one input parameter compared to another. A neural network is designed to make decisions, such as whether or not an observation detects a planet, based on a set of input parameters that treat, e.g. the shape and depths of a light curve, the noise and systematic error, such as star spots. The discriminative nature of this deep net can only make a qualitative assessment of the candidate signal by indicating the likelihood of finding a transit within a subset of the time series. The advantage of a deep net is that it can be trained to identify very subtle features in large data sets. This learning capability is accomplished by algorithms that optimize the weights in such a way as to minimize the difference between the output of the deep net and the expected value from the training data. The network does not rely on hand designed metrics to search for planets, instead it will learn the optimal features necessary to detect a transit signal from training data. Neural networks have been used in a few planetary science applications including those for transit detection and atmospheric characterization (Waldmann 2016; Kipping & Lam 2017; Pearson et al. 2018). A notable example used a neural network to validate the detection of a few Kepler planets by assessing the transit-like shape of the data (Shallue & Vanderburg, 2018). This study adopts a different approach and uses deep nets to vet transit signals prior to fitting the light curves. The extensive TESS target list requires an efficient analysis in order to prioritize follow-up projects.

A convolutional neural network (CNN) is used to analyze TESS time series data and and it is based on the CNN design in Pearson et al. (2018). A convolutional neural network (CNN) is well suited for processing time series measurements since the input data is correlated to one another through time. The CNN computes local features in the data by convolving it with multiple filters. Afterwards, the data is down sampled using an average pooling layer which mimics the effect of binning in time to reduce photometric scatter and decreases the number of weights in the next layer. After passing through two convolutional layers the data is processed using 3 layers of neurons which are connected fully. The final CNN network has $\sim$23,000 free parameters (i.e. the weights within all the neurons) that are optimized from training data. Increasing the number of weights can cause a neural network to over-fit where it learns to classify examples based on features not associated with the signal (e.g. the noise). Therefore, a set of new set of light curves are created only to validate and test the neural network. Optimizing the weights of the CNN is done by training the network on real TESS data with injected transits. Afterwards, a recovery test is performed to examine the sensitivity of detecting transits at low SNR. The CNN was trained using 20000 random light curves with varying SNR where the injected noise was based on the residuals to TESS light curves.

The training data is based on TESS residuals which have the transit signal removed but variations due to instrumental and stellar variability still included. Residuals were computed from TOI candidates with SNRs greater than 3 because removing the transit signal is more reliable and less prone to biases from the noise. Light curves with non-detections were avoided because extra processing is required to evaluate if the data includes e.g. a highly variable star, a neighboring and contaminating light source or an instrumental effect. All of which can perturb the data in a manner similar to a transit and thus bias the algorithm during training. The neural network was trained negatively with data of SNR varying between -1 and 0.25. The negative training refers to the training label being 0, indicating no transit is present in the data. The negative SNR indicates an inverted transit signal meant to simulate a systematic. The network is trained negatively on transits with an SNR between 0 and 0.25 because it decreases the false positive rate due to small systematic effects that decrease flux in a manner similar to a transit (e.g. pointing induced variations, stellar effects). Consequently the positive training samples come from transits with an SNR between 1 and 2, typically these transits are easily visible to the human eye. The neural network is trained on batches of 32 light curves from a sample of 10,000 positive and 10,000 negative samples. Afterwards, the network is tested on data it has not seen before in order to assess the sensitivity of detecting low SNR transits. Figure 3 shows the results of the transit recovery test. The neural network was cautious in selecting transits and has a 15.7$\%$ false negative rate because a high tolerance (0.8) is used for classifying the transit probability. Additionally, the network has a 1.35$\%$ false positive rate which considers inverted transit shapes and other sources of systematics perturbing the light curve. The network exhibits a smaller false positive rate if the data is detrended prior to prediction. The CNN algorithm has an average detection accuracy of 98$\%$ or more for transits that are at least 1.25 $\times$ greater than the noise. For more information on the training and design of the deep net here please see the online material222https://github.com/pearsonkyle/Exoplanet-Artificial-Intelligence.

Integrating the convolutional neural network into the analysis is essential to expedite the processing of optimal looking transits that would allow for robust constraints on the system parameters. This vetting algorithm is applied to TESS data by computing a binned phase folded light curve and using the CNN to evaluate the probability a transit is present. If the transit probability is greater than 80$\%$ the light curve will be characterized. The phase folded data is truncated to surround the primary transit with $\pm 1$ transit duration before and after ingress, respectively. The data is then interpolated onto a grid of 180 data points because the neural network has a fixed input size of 180. Pearson et al. (2018) found that interpolation of a longer time series did not significantly affect the performance of a CNN for simulated data and Kepler light curves, the same is assumed to be true for TESS. A portion of the TOI candidates that were removed from the analysis exhibit either low SNR transits ($<\sim$1) or had systematic effects in the light curves that would not allow for a precise transit measurement.

3.2 Light Curve Model and Nested Sampling

After the light curves are vetted, the phase folded light curve is characterized and the transit parameters are derived. If the SNR of the phase folded light curve is above 2, the individual transits in the time series are fit in order to derive a new ephemeris. The SNR is defined as the amplitude of the transit depth normalized by the standard deviation of the residuals to the phase folded light curve. In most cases, the SNR is approximately $(rp/rs)^{2}/\sigma_{res}$ unless the transit is highly inclined in which case the amplitude of the transit depth comes from the transit model itself. The transit parameters are derived by optimizing a model’s parameters using observations and nested sampling. The following function below is used to maximize the likelihood of a transit model and non-linear correction simultaneously:

[TABLE]

Here $F_{obs}$ is the flux recorded on the detector, $F_{transit}$ is the actual astrophysical signal (i.e. the transit light curve, given by a Mandel & Agol (2002) model light curve), $a_{i}$ are quadratic correction coefficients. The quadratic term removes non-linear variations in the light curve that would otherwise bias the mid-transit estimate. As the planet transits in front of the host star, brightness contrasts between the stellar limb to modulate the shape of the transit. A quadratic limb darkening model is used and computed from the tool LDTk and the stellar parameters in Table 2 (Parviainen & Aigrain, 2015).

The transit parameters are optimized using the multimodal nested sampling algorithm called MultiNest (Skilling 2006; Feroz & Hobson 2008; Feroz et al. 2009). MultiNest is a Bayesian inference tool that uses the Monte Carlo strategy of nested sampling to calculate the Bayesian evidence allowing simultaneous parameter estimation and model selection. A nested sampling algorithm is efficient at probing parameter spaces which could potentially contain multiple modes and pronounced degeneracies in high dimensions; a regime in which the convergence for traditional Markov Chain Monte Carlo (MCMC) techniques becomes incredibly slow (Skilling 2004; Feroz & Hobson 2008). Optimization of the hyperparameters for nested sampling enable an efficient search for the global solution while producing a numerical uncertainty from sampling the posterior distribution. 500 live points are used with an evidence tolerance of 0.1 and a sampling efficiency of 50% to ensure enough points in the prior space are sampled for convergence and to reveal any multimodal posterior distributions. When fitting the phase folded light curve the following parameters are allowed to vary: $R_{p}/R_{*}$, $a/R_{*}$, $inc$, $T_{mid}$, $a_{0}$, $a_{1}$, $a_{2}$. After the phase folded light curve is characterized, individual transits in the time series are characterized for only $a/R_{*}$, $T_{mid}$, $a_{0}$, $a_{1}$, $a_{2}$ while leaving the other values fixed and adopting them from the phase folded light curve fit (See Table 4).

3.3 Linear Ephemeris Model

The orbits of bodies in multi-planet systems are subject to additional gravitational forces than just the star which can cause time-dependent deviations in the orbit parameters. These deviations yield measurable differences in orbital periods and transit durations. For transiting exoplanets, measuring the orbital period is done by fitting a linear function to mid-transit values. The next mid-transit or mid-eclipse time can be calculated using:

[TABLE]

where $t_{next}$ is the time of the next transit, $P$ is the orbital period of the planet, $T_{mid}$ is the planet’s epoch of mid transit and $n_{orbit}$ is the number of orbits that have occurred between $t_{next}$ and $T_{mid}$. Equation 2 is used to derive the period and epoch by fitting for those parameters using nested sampling. Discussed more in the previous section, MultiNest is a Bayesian inference tool that uses nested sampling to calculate the Bayesian evidence alongside enabling posterior inference. The Bayesian evidence is compared between a linear and non-linear ephemeris, which accounts for perturbations induced by planetary companions.

3.4 Non-Linear Ephemeris using N-body simulations

Planetary companions will perturb the orbits of a transiting planet in a manner indicative to the period and mass of the perturber. To predict the transit times in a system with multiple planets a non-linear ephemeris is computed by:

[TABLE]

where $f_{Nbody}$ is the transit time returned from an N-body simulation, $P_{1}$ is the period of the transiting planet, $P_{2}$ is the period of the companion planet, $M_{2}$ is the mass of the companion planet and $\Omega_{2}$ is the argument of periapsis and $e_{2}$ is the eccentricity of the companion planet and $T_{1,mid}$ is the epoch of mid-transit for the transiting exoplanet.

Computing orbital perturbations on the order of a minute or less requires precise calculations therefore a full N-body simulation is used. REBOUND is an open-source N-body code that features the IAS15 integrator, a 15th-order integrator well suited to simulate gravitational dynamics (Rein & Liu, 2012). The IAS15 integrator is designed to handle close encounters and high-eccentricity orbits by preserving the symplecticity of Hamiltonian systems using an adaptive time step and a 15th-order modified Runge-Kutta integrator(Rein & Spiegel, 2015). 30 minutes is chosen as a default time step for the N-body simulation. A smaller time step of 1 minute was compared against 30 minutes and found to have perturbations within a few seconds of one another. When computing the transit timing variations the position information of the planet is interpolated linearly between time steps to achieve a precision smaller than 30 minutes. Past literature suggests using a time step at least 1/20 the period therefore 30 minutes is chosen to optimize for speed however it limits the analysis to transits with orbital periods longer than 10 hours (Nesvorný & Morbidelli, 2008). An example N-body simulation is shown in Figure 5.

N-body simulations contain numerous parameters per planet and the optimization of them can yield millions of combinations to search through. Not every parameter influences the perturbations in a measurable way given the observational precision and therefore certain assumptions can be made about the system to optimize a retrieval. Figure 6 shows how each parameter in an N-body simulation perturbs the periodicity of a transiting planet. TTV signals within the phase space tested do not exhibit a strong sensitivity to inclination of either planet. Therefore, the inclination is left as a fixed parameter during an N-body retrieval. The inclination value of the transiting planet is adopted from the phase folded light curve fit while the inclination of the perturbing planet is fixed to 90 degrees. The eccentricity of the perturbing planet is however left as a free parameter. The eccentricity of the inner planet is fixed to the literature value when known or assumed to be 0, yielding an upper limit on the perturbing planet’s mass. The eccentricity of the perturbing planet has an interesting affect on the structure of the TTV signal by making the peaks wider and troughs narrower or vice versa (see Figure 6). There is a degeneracy between the eccentricity of the transiting planet and the mass of a perturbing planet. Essentially, the TTV signal could be from a more massive perturbing planet or an eccentric transiting planet. The easiest way to constrain this degeneracy is by constraining the eccentricity of the transiting planet via secondary eclipse or radial velocity measurement. The N-body retrieval is computed using the Bayesian inference tool Multinest. Typically, the retrieval computes between 5000 and 10000 simulations before converging. Unstable orbits cause the period of the perturber to vary significantly which affects the measured TTV signal up to hours or more. If the solution is unstable, it is returned to the retrieval with a chi-squared value twice as large, this helps constrain the search space for new solutions.

3.5 Estimation of N-body Priors with Machine Learning

Matching data with a perturbation model requires testing thousands of N-body simulations in order to find the best-fit parameters and their uncertainties for a perturber. Machine learning can estimate where to search in parameter space before running a more extensive analysis (e.g. using nested sampling). The N-body model depends on at least a few basic parameters; Mass of the star (M∗), mass of the inner planet (M1), period of the inner planet (P1), mass of the outer planet (M2), period of the outer planet (P2), eccentricity of the outer planet (e2), and argument of periastron for the outer planet (w2). Typically, parameters pertaining to the star and inner planet are known ahead of time. While figure 6 shows that TTVs of the transiting planet are not influenced by its mass, that may not hold for non-zero eccentricities or variation of the mass of one order of magnitude (e.g. from super-Earths to Jupiter-like planets), and M2 will be influenced by the choice of M1. When the mass of M1 is not known aprior it is estimated using an empirical mass radius relation. The empirical mass radius relation uses all of the transiting exoplanets with mass measurements from the NASA Exoplanet Archive. A Gaussian Process is used to estimate a general trend between mass and radius (see Figure 7).

A machine learning algorithm is designed to predict the parameters of the perturbing body (M2, P2, e2, w2) given the known parameters (M∗, M1, P1, e1, w1) and measured perturbation (i.e. measured mid-transit values subtracted from a linear fit). To do this a dual-input, multi-output regression model is created. The independent features (M∗, M1, P1) are analyzed using a fully connected neural network while the time-dependent features are analyzed using a convolutional neural network. The output of each branch is piped into a fully connected neural network and finally, 4 parameters are predicted representing the planetary and orbital elements of the perturbing body.

Training data is simulated in order to optimize the neural network for analyzing high cadence data from TESS observations. The TESS mission is conducting an all sky survey and will observe portions of the sky for as many days at 180 but also as little as 27 depending on the RA and Dec 333https://heasarc.gsfc.nasa.gov/docs/tess/operations.html. Therefore random simulations of short-period planets around solar type stars are created to train the neural network. Stellar masses range between 0.85 and 1.15 solar since the targets of interest in this study are around solar type stars. Inner planets are generated from a uniform distribution between 0.66-150 MEarth and a period between 0.75 and 10 days. The outer planet is generated from a conditional probability distribution in order to avoid randomly generating unstable orbits. The outer planet is randomly generated with a mass ratio between 0.25 and 3 with respect to the inner planet. The period of the outer planet is generated with a period ratio between 1.25 to 2.5 but will be larger if the random period is within the hill sphere of the first planet. One last constraint is in place to prevent generating unstable orbits; only simulations producing TTVs less than 10 minutes are used.

The constraints on the training data are tailored for this specific study but the code can be easily modified to account for larger TTVs or smaller mass stars. Two variations are made to the training data in order to improve the network’s accuracy during training. For every parameter combination (M2, P2, e2) the argument of periastron is uniformly distributed between -$\pi$ and +$\pi$ with 16 points. Additionally, for each parameter combination the O-C data is varied according to a Gaussian distribution with a standard deviation 10$\%$ of the TTV amplitude. Varying the O-C data simulates acquiring different observations. The neural network is trained with 100,000 samples and 100 epochs using the mean squared error as the loss function. After training the neural network, the error and range of reliability are characterized as a function of orbital period ratio and companion mass (See figure 8). The neural network errors are used to estimate the priors in a MultiNest Fit (described in section 3.4) by extending the errors $\pm$2 sigma (see Figure 9).

4 Results and Discussion

Searching through the TESS objects of interest (TOI) catalog yielded information about the planetary ephemerides and transit parameters of 74 targets. A subset of the catalog is studied for single planet systems with high cadence data, more than 3 orbits, transits longer than 30 minutes, transit probability greater than 80$\%$ from a classifying neural network (see Table 3) and a transit SNR greater than at least 2. These constraints enable high precision measurements of multiple transit events which are used to search for planets by measuring transit timing variations. 148 different targets entered the pipeline, 30 failed the vetting algorithm and of the 118 left only 74 had an SNR greater than 2 and at least 3 transit measurements. Figure 10 shows the photometric variability of TESS stars in this sample as a function of their magnitude. A time series analysis is conducted to derive a new ephemeris for each target. Ephemeris residuals that exhibit periodic structure undergo an in depth analysis to explore perturbing bodies. A neural network is used to estimate the priors for a retrieval and a nested sampling algorithm derives posteriors for an N-body model.

4.1 WASP-126

The first target of interest is WASP-126 because it shows subtle periodic structure in the residuals of the linear ephemeris. WASP-126 b is a 0.28 MJup mass planet at an orbital period of 3.28 day around a G2 5800K star (Maxted et al., 2016). Not much is known beyond the initial discovery of WASP-126 b however it does have a large scale height ($\sim$650 km) and would be an interesting candidate for atmospheric characterization. The best fit parameters for the phase folded light curve are reported in Table 4. The time series analysis of WASP-126 b is shown in figure 4.

A search for contaminating light sources affecting the transit depth is performed by cross-referencing the TIC catalog for the position and magnitude of neighboring stars. The plate scale of TESS is 21”x21” which corresponds to 50$\%$ and 90$\%$ of a star’s flux being contained within a 1×1 and 4×4 pixel region around the centroid, respectively (Ricker et al., 2015). The large pixel scale helps TESS in their all sky survey of transiting planets but limits the ability to spatially resolve stars therefore pixels may contain one or more light sources. WASP-126 (T = 10.55) has two stars within an arc minute, TIC 25155316 (T = 14.14) and TIC 25155311 (T = 16.06) at a distance of 26.5” and 7.9” respectively. Compared to WASP-126 the nearby sources are 27 and 160 times fainter, respectively. Given the proximity to WASP-126 these objects contribute $\sim$1$\%$ of the star’s flux within the aperture used which is comparable to the transit signal of $\sim$0.6$\%$. It is unlikely the transit signal is a systematic effect induced by contamination because the stars would have to completely disappear or move out of the aperture in order to simulate the transit. The stars are at most 1.26 pixels away from the target and the aperture size is roughly 3 pixels in radius. Since the median displacement of the centroid in time is 0.02 pixels with the max being 0.05 pixels, it unlikely the transit signal is being modulated by pointing induced systematics from the contaminating stars.

The transit times for WASP-126 b show periodic structure (see O-C data in figure 4) and to test the significance of a perturbing planet the data are analyzed with an N-body retrieval described in Section 3.4. The posteriors from the N-body analysis are shown in Figure 11. The Bayesian evidence for a linear ephemeris is -18.9 and the Bayesian evidence for a non-linear ephemeris is -15.9, suggesting a non-linear ephemeris is a better fit. The greatest evidence is for a perturbing planet at an orbital period of 7.65 $\pm$ 0.27 days with a mass of 64.0$\pm$26.9 MEarth. The data is weakly sensitive to the eccentricity of the perturbing body, which is estimated to be $<$ 0.05 (see Table 5). The upper limit on the prior was determined by the machine learning algorithm discussed in Section 3.5. The perturbing body should be detectable with future RV observations and it corresponds to an RV semi-amplitude of $\sim$11–27 m/s.

There is no evidence for another transiting companion in the WASP-126 system (see Figure 4). A transit periodogram analysis is conducted on the time series after removing the primary transit signal from every epoch. The large spike in the periodogram is a result of motions from the space craft and seen in every target analyzed, sadly it’s not a planet signal. In order for the planet candidate, WASP-126 c, to not transit it must have an inclination less than $\sim$85.5. The inclination constraint places a lower limit on the mutual inclination of the system to be 4 degrees.

4.2 WASP-18

WASP-18 is massive (10.4 MJup) Jupiter-sized planet orbiting a 6400 K F type star at an orbital period of 0.94 day (Hellier et al. (2009); Southworth et al. (2009)). WASP-18 b is an interesting candidate for studying tidal interactions and star-planet interactions because of its short orbital period and extreme mass (Csizmadia et al. (2019); Barker & Ogilvie (2010); Arras et al. (2012); Miller et al. (2012); Albrecht et al. (2012); Wilkins et al. (2017)). WASP-18 b was observed with Spitzer and found to have a brightness temperature of $\sim$3100 K which yields a secondary eclipse depth of 0.38$\%$ (Nymeyer et al. (2011); Maxted et al. (2013)). Additionally, the transit times from Maxted et al. (2013) show deviations on the order of $\sim$60 seconds, the same level as those seen here with TESS. Recently, a full-orbit phase curve in the optical was measured with TESS (Shporer et al., 2018). The optical phase curve yielded a low albedo, lack of atmospheric phase shift and inefficient heat distribution from day to night hemispheres. The phase curve variations affect the measurement of the transit parameters since it adds a non-linear component into the baseline (e.g. see the phased-folded data in Figure 13). As opposed to modeling the phase curve, which is beyond the scope of this project, the non-linearity is accounted for by simultaneously fitting a quadratic function along with the light curve model (see Equation 1). This approach yields a transit depth of 9293 $\pm$ 21 ppm which is inconsistent to 4.26 $\sigma$ with the 9439${}^{27}_{26}$ ppm derived in Shporer et al. (2018). The inconsistency is due to differences in limb darkening parameters considering the other parameters like inclination and a/Rs are consistent. A separate analysis with the same limb darkening parameters as Shporer et al. (2018) yielded results consistent within 0.5 $\sigma$.

The effects of contaminating sources can alter the transit depth of a planet by making it appear smaller than it actually is. Shporer et al. (2018) indicated two moderately bright stars in the TESS catalog that are within the vicinity of WASP-18, one of targets in this study. The stars TIC 100100823 and TIC 100100829 are 73“ and 83“ away from WASP-18 and have TESS magnitudes of T = 12.65 mag and T = 12.50 mag, respectively. Compared to WASP-18 (T = 8.83 mag) the nearby sources are 34 and 29 times fainter, respectively. Shporer et al. (2018) found the stars would contribute minimally (0.01$\%$) to the flux of WASP-18 so the same aperture size (that of SPOC) is adopted here as the former study on WASP-18 b’s phase curve.

The residuals of a linear ephemeris display transit timing variations on the order of 30 seconds with measurements being correlated in time. Shporer et al. (2018) also searched for TTVs in the WASP-18 system and concluded there was more evidence for no TTV. The findings here contradict their result because different models are used to assess the signal in the data. Shporer et al. (2018) use a parameterized sin wave to fit the residuals of their linear ephemeris. Here, the transit times themselves are fit with a more complicated model rather than to the residuals in a disjointed approach. The Bayesian evidence is then compared between a linear and non-linear ephemeris to assess the significance of another planet being present in the system. For consistency reasons with respect to the data acquisition and time keeping, past measurements of WASP-18 b are ignored and only the TESS measurements are used to assess the TTV signal. The N-body analysis uses an eccentricity value of 0 for WASP-18 b which is consistent with the value measured from Spitzer of 0.0091 $\pm$ 0.0012 (Nymeyer et al., 2011). The best fit companion from the retrieval has a mass of 55.2 $\pm$ 12.3 $M_{Earth}$ and a period of 2.155 $\pm$ 0.006 day. The global log evidence from MultiNest is reported in the legend of Figure 15 and is greater by 10.9 for the non-linear ephemeris. Despite the close orbital distances and extreme mass of WASP-18b, the hill spheres of WASP-18 and the perturber do not interact. The hill sphere for WASP-18 b is 0.003 AU and it has an orbital distance of 0.02 AU, where as the companion has a hill radius of 0.0015 AU and an orbital distance of 0.035 AU.

There is no evidence for additional transiting companions in the system based on a periodogram analysis. The periodogram shows a spike around 13.8 days consistent with the orbits of the spacecraft and some minor signals around 5 days, which are suspect to come from the star. There is a rotation period for the star measured around 5.6 days (Hellier et al. 2009; Maxted et al. 2013). The inclination limit for a transiting planet at 2.12 days is 80.6 degrees and is an upper limit for the true value. The companion would cause a RV signal of $\sim$18–29 m/s which is small compared to the $\sim$ 1800 m/s value for WASP-18 b. The companion’s RV signal is interestingly close to the predicted amplitude of RV jitter due to tidal distortions in the host star from WAPS-18 b (Arras et al., 2012). It seems possible to detect a companion in this system considering the 9 RV measurements from Hellier et al. (2009) have an error of $\sim$ 10 m/s. The most recent analysis from Bonomo et al. (2017) reports a semi-amplitude of WASP-18 b to be 1817.2 $\pm$ 2.3 m/s with an upper limit on the RV jitter of the star to be 3.5 m/s. Additionally, the measurements in Triaud et al. (2010) exhibit scatter on the order of $\sim$10-25 m/s. However, it’s difficult to say if the past measurements resolve the RV peak from the companion because of aliasing and whether they’re consistent with the results here. In which case, it would be worth performing a joint simultaneous fit between the TTV signal and RV measurements in order to further constrain the companion’s mass and orbit. The orbital constraints placed on the outer candidate here can be used to optimize follow-up RV measurements and prevent errors from aliasing the signal.

4.3 TOI 193

The last object of interest is a candidate for a multi-planet system, TOI 193 or TIC 183985250. TOI 193 is a G type star with a radius of 0.99 RSun. Since TOI 193 is an object of interest to the TESS community and active follow-up is being conducted to constrain the system properties a few assumptions are made. For the N-body simulation the star is assumed to have a mass of 1 MSun and the transiting planet is assumed to have an eccentricity of 0. The phase folded light curve yields a radius of 4.9 REarth which corresponds to a mass of $\sim$22 MEarth based on the Gaussian process in Figure 7. The detection of the transit was a little more tenuous then the others in this study because it has an SNR of 2.65. The time series analysis is shown in Figure 16. A search for contaminating light sources around TOI 193 (TIC 183985250) yielded one star within 60“ (about the photometric aperture size). The neighboring star (TIC 183985251) is 8.2“ away with a magnitude of T = 13.957. The star is 88 times dimmer but given the distance to the target it corresponds to $\sim$1.1$\%$ of the flux in the aperture. The average displacement in x and y of the target TOI 193 is 0.025 pixels with the max being 0.11 pixels. The contaminating star likely has a similar point error and is going to be in the same pixel as the target most of the time. Therefore shifts on the CCD are unlikely to induce the transit signal seen because the contaminating star would have to drift out of the aperture and back in within a transit duration.

There is more evidence for the non-linear ephemeris compared to a linear ephemeris however the non-linear ephemeris is degenerate. The final parameters for each ephemeris are shown in Table 5. Two degenerate TTV modes are plotted in the posteriors (Figure 17) because each solution has a Bayesian evidence greater than the linear ephemeris. Each perturbation solution is close to the orbital resonance of 2:1. The two solutions look similar (see Figure 18) but yield noticeable different RV signals. Future RV observations will be able distinguish the correct mode because the companion can cause as much variation in the RV semi-amplitude as 54 m/s but as little as $\sim$15, depending on the mass. There is no evidence for a transiting companion at an SNR greater than 1 in the residuals of the time series. The periodic structure in the transit periodogram is an aliasing effect of the short period transit being phase folded into the data.

5 Conclusion

The Transiting Exoplanet Survey Satellite (TESS) objects of interest (TOI) are analyzed in search of transit time variations. Data from sectors 1-3 are used to analyze 74 planet candidates from the TOI catalog. The average photometric precision of candidate hosting stars are measured. Stars with TESS magnitudes between $\sim$9–11 exhibit noise on the order of $\sim$1000 ppm. The extensive TESS candidate list of over 3 million targets requires an efficient analysis in order to produce results in a timely manner which can be used to refine follow-up projects (e.g. Exoplanet Transit Survey; Zellem et al. (2019)). Artificial intelligence is used to vet a portion of non-transit signals in order to expedite the characterization of numerous transits and ephemerides. A convolutional neural network is trained on time series measurements from TESS and has a transit detection accuracy of 99$\%$ with transits over an SNR of 2. Preliminary excavations using the AI vetting algorithm on the TESS candidate target list (CTL) yielded a large amount of eclipsing binary signals. Means to distinguish the eclipsing binaries from transit signals will be implemented in future versions of the AI. After characterizing individual transit measurements, residuals of a linear ephemeris are used to search for transit timing variations. To assess the significance of a perturbing planet the Bayesian evidence is compared between a linear and a non-linear ephemeris, which is based on an N-body simulation. There is evidence for 3 new multi-planet systems with non-transiting companions using the 2 minute cadence observations from TESS. Two previously known exoplanets, WASP-18 b and WASP-126 b are studied and there is evidence for WASP-18 c at a period of 2.155 $\pm$ 0.006 day with a mass of 55.2 $\pm$ 12.3 MEarth and WASP-126 c a 64.2 $\pm$24.6 MEarth planet at an orbital period of 7.63$\pm$0.17 day. The third system with a significant TTV is TOI 193.01 with a candidate planet, TOI 193.02, exhibiting a degenerate solution near the 2:1 resonance with a mass of either 39.4$\pm$9.5 or 73.4$\pm$28.1. Radial velocity observations in the future will be able to distinguish between the degenerate solutions because the RV semi-amplitude will vary between at least $\sim$15–54 m/s depending on the perturbing planet’s mass. Results for these targets will be uploaded to the exoplanet TESS follow up program to encourage further analysis and collaboration with the community. TESS provides us with an exciting opportunity to perform high cadence observations of many new and existing exoplanet systems during the all sky survey.

6 Data Availability

Photometric data from the Transiting Exoplanet Survey Satellite is available online: https://archive.stsci.edu/tess/bulk_downloads.html.

7 Code Availability

The machine learning algorithm used to vet transit candidates is available here: https://github.com/pearsonkyle/Exoplanet-Artificial-Intelligence. The algorithm used to model exoplanet transits with nested sampling is available online: https://github.com/pearsonkyle/Exoplanet-Light-Curve-Analysis. The N-body retrieval algorithm is available online: https://github.com/pearsonkyle/Nbody-AI.

8 Acknowledgements

I would like to thank my advisor, Caitlin Griffith, for supplying me with my first research assistantship (RA) so I could focus on writing my thesis proposal even though most of my time was spent on this study. I thank Eli Shneider and the whole Bently’s crew for fueling parts of this research and for providing me with a location to work when campus was too distracting. Daynal Tansu and Avi Shporer for their discussion about WASP-18 b. Adam Sutherland for putting up with my questions about orbital dynamics while on the ski lift. Last but not least, the anonymous referee for providing helpful comments which strengthened the paper.