TOI-216b and TOI-216c: Two warm, large exoplanets in or slightly wide of the 2:1 orbital resonance

TL;DR

This paper confirms and characterizes a pair of warm, large exoplanets orbiting TOI-216, exploring their possible resonant configurations and implications for planetary formation theories.

Contribution

It provides the first detailed analysis of the TOI-216 system, identifying two potential orbital solutions and discussing their formation scenarios and future observational prospects.

Findings

Two families of orbital solutions identified.

Preference for the Jupiter-Neptune resonance configuration.

System's eccentricities and inclinations suggest undetected planets.

Abstract

Warm, large exoplanets with 10-100 day orbital periods pose a major challenge to our understanding of how planetary systems form and evolve. Although high eccentricity tidal migration has been invoked to explain their proximity to their host stars, a handful reside in or near orbital resonance with nearby planets, suggesting a gentler history of in situ formation or disk migration. Here we confirm and characterize a pair of warm, large exoplanets discovered by the TESS Mission orbiting K-dwarf TOI-216. Our analysis includes additional transits and transit exclusion windows observed via ground-based follow-up. We find two families of solutions, one corresponding to a sub-Saturn-mass planet accompanied by a Neptune-mass planet and the other to a Jupiter in resonance with a sub-Saturn-mass planet. We prefer the second solution based on the orbital period ratio, the planet radii, the lower free eccentricities, and libration of the 2:1 resonant argument, but cannot rule out the first. The free eccentricities and mutual inclination are compatible with stirring by other, undetected planets in the system, particularly for the second solution. We discuss prospects for better constraints on the planets' properties and orbits through follow-up, including transits observed from the ground.

Tables

Table 1. Table 1: Stellar Parameters a a As a summary statistic we report the median and 68.3% confidence interval of the posterior distribution. for TOI-216

Catalog Information
Parameters	Value	Source
R.A. (h:m:s)	04:55:55.3	Gaia DR2
Dec. (d:m:s)	$-63$:16:36.2	Gaia DR2
Epoch	2015.5	Gaia DR2
Parallax (mas)	$5.59\pm 0.03$	Gaia DR2
${\mu }_{\mathrm{ra}}$ (mas yr-1)	$-22.7\pm 0.04$	Gaia DR2
${\mu }_{\mathrm{dec}}$ (mas yr-1)	$-56.355\pm 0.05$	Gaia DR2
$g$ magnitude	12.163126
Gaia DR2 ID	4664811297844004352
TIC ID	55652896
TOI ID	216
TIC TESS magnitude	11.504
$V$ magnitudebbUsing the relationship derived by Jordi et al. (2010), we compute the $V$ magnitude from the Gaia $g$ magnitude and the Johnson-Cousins $I_C$ magnitude. We estimate the $I_C$ magnitude to be the TESS magnitude, because the two band passes have the same center (Ricker et al., 2015).	12.393
Spectroscopic properties
Parameters	Spectrum	TakedaccTakeda et al. (2007)	DartmouthddDotter et al. (2008)	+GaiaddDotter et al. (2008)
Stellar effective temperature, $T_{\mathrm{eff}}$ [K]	5045$\pm$110	50560${}_{-1120}{}^{}{}_{}{}^{+1100}$	50540${}_{-1200}{}^{}{}_{}{}^{+1030}$	50890${}_{-450}{}^{}{}_{}{}^{+430}$
Iron abundance, [Fe/H]	-0.16 $\pm$0.09	-0.16$\pm 0.08$	-0.16$\pm 0.09$	-0.15${}_{-0.09}{}^{}{}_{}{}^{+0.08}$
Surface gravity, $\log g[$cm s-2]	4.53$\pm$0.12	4.578${}_{-0.023}{}^{}{}_{}{}^{+0.02}$	4.58${}_{-0.04}{}^{}{}_{}{}^{+0.03}$	4.58600${}_{-0.0350}{}^{}{}_{}{}^{+0.003}$
Stellar mass, $M_{\star }$ [$M_{\odot }$]		0.78${}_{-0.02}{}^{}{}_{}{}^{+0.04}$	0.76${}_{-0.03}{}^{}{}_{}{}^{+0.04}$	0.77${}_{-0.03}{}^{}{}_{}{}^{+0.03}$
Stellar radius, $R_{\star }$ [$R_{\odot }$]		0.765${}_{-0.02}{}^{}{}_{}{}^{+0.023}$	0.74${}_{-0.03}{}^{}{}_{}{}^{+\pm 0.043}$	0.747${}_{-0.014}{}^{}{}_{}{}^{+0.015}$
Stellar density, ${\rho }_{\star }$ [${\rho }_{\odot }$]		1.812${}_{-0.146}{}^{}{}_{}{}^{+0.14}$	1.995${}_{-0.230}{}^{}{}_{}{}^{+0.213}$	1.84${}_{-0.15}{}^{}{}_{}{}^{+0.14}$

Table 2. Table 2: Observation Log

TOI-216	Date	Telescope${\dagger}$${\dagger}$Telescopes: LCO-CTIO-1.0: Las Cumbres Observatory - Cerro Tololo Astronomical Observatory (1.0 m) LCO-SSO-1.0: Las Cumbres Observatory - Siding Spring (1.0 m) LCO-SAAO-1.0: Las Cumbres Observatory - South African Astronomical Observatory (1.0 m) LCO-SSO-0.4: Las Cumbres Observatory - Siding Spring (0.4 m) Myers-T50: Siding Spring Observatory - T50 (0.43 m) Hazelwood-0.3: Stockdale Private Observatory - Victoria, Australia (0.32 m) El Sauce-0.36: El Sauce Private Observatory - El Sauce, Chile (0.36 m)	Filter	ExpT	Exp	Dur.	Transit	Ap. Radius	FWHM
	(UTC)			(sec)	(N)	(min)	expected coverage${\ddagger}$${\ddagger}$Observations did not detect a transit event because they were scheduled using the initial public TESS linear ephemeris. The TTV offset from the linear ephemeris is now known to be larger than the time coverage of the observations.	(arcsec)	(arcsec)
b	2018-11-22‡	LCO-SSO-0.4	i′	90	54	100	Ingress+30$\%$	8.5	7.5
	2018-12-09‡	Myers-T50	Lum	60	200	240	full	8.3	4.6
	2018-12-26‡	LCO-SSO-1.0	i′	30	85	99	Ingress+25$\%$	7.0	2.8
	2019-01-29‡	LCO-SAAO-1.0	r′	100	97	225	Full	9.3	2.4
	2019-01-29‡	LCO-SAAO-1.0	i′	25	181	198	Full	5.8	2.2
	2019-02-15	LCO-SSO-1.0	Zs	60	160	236	Out-of-Transit	4.7	2.0
c	2018-12-16	LCO-SAAO-1.0	i′	90	75	180	Egress+60$\%$	5.8	2.5
	2018-12-16	LCO-SAAO-1.0	i′	39	331	450	Full	5.8	2.1
	2019-01-20	Hazelwood-0.3	g′	240	101	449	Egress+70$\%$	5.5	3.2
	2019-02-23	LCO-SAAO-1.0	Zs	60	148	212	Out-of-Transit	6.2	2.5
	2019-02-24	LCO-CTIO-1.0	Zs	60	150	213	In-Transit	6.2	2.5
	2019-02-24	El Sauce-0.36	Rc	30	514	303	Egress+90$\%$	5.9	3.7
	2019-02-24	LCO-SSO-1.0	Zs	60	81	117	Out-of-Transit	6.2	2.5

Table 3. Table 3: Planet Parameters for TOI-216b and TOI-216c Derived from the Light Curves

Parameter	ValueaaAs a summary statistic we report the median and 68.3% confidence interval of the posterior distribution.
TOI-216b
Planet-to-star radius ratio, $R_p/R_{\star }$	0.11	${}_{-0.02}{}^{}{}_{}{}^{+0.04}$
Planet radius, $R_p$ [$R_{\oplus }$]	8.6	${}_{-1.9}{}^{}{}_{}{}^{+2.9}$
Light curve stellar density, ${\rho }_{\mathrm{circ}}$ [${\rho }_{\odot }$]	1.13	${}_{-0.19}{}^{}{}_{}{}^{+0.29}$
$a/R_{\star }$bbIf the planet’s orbit is not circular, this corresponds to the average planet-star-separation during transit divided by the stellar radius.	29.1	${}_{-1.8}{}^{}{}_{}{}^{+2.3}$
Impact parameter, $\left|b\right|$	0.99	${}_{-0.04}{}^{}{}_{}{}^{+0.05}$
Sky-plane inclination, $i_{\mathrm{sky}}$ ${[}^{\circ }]$	88.0	${}_{-0.2}{}^{}{}_{}{}^{+0.2}$
Mid-transit times	1325.328	${}_{-0.004}{}^{}{}_{}{}^{+0.003}$
	1342.431	${}_{-0.003}{}^{}{}_{}{}^{+0.003}$
	1359.539	${}_{-0.003}{}^{}{}_{}{}^{+0.003}$
	1376.631	${}_{-0.003}{}^{}{}_{}{}^{+0.003}$
	1393.723	${}_{-0.003}{}^{}{}_{}{}^{+0.003}$
	1427.879	${}_{-0.003}{}^{}{}_{}{}^{+0.003}$
	1444.958	${}_{-0.003}{}^{}{}_{}{}^{+0.003}$
	1462.031	${}_{-0.003}{}^{}{}_{}{}^{+0.003}$
	1479.094	${}_{-0.003}{}^{}{}_{}{}^{+0.003}$
	1496.155	${}_{-0.003}{}^{}{}_{}{}^{+0.003}$
	1513.225	${}_{-0.003}{}^{}{}_{}{}^{+0.003}$
TOI-216c
Planet-to-star radius ratio, $R_p/R_{\star }$	0.1236	${}_{-0.0008}{}^{}{}_{}{}^{+0.0008}$
Planet radius, $R_p$ [$R_{\oplus }$]	10.2	${}_{-0.2}{}^{}{}_{}{}^{+0.2}$
Light curves stellar density, ${\rho }_{\mathrm{circ}}$ [${\rho }_{\odot }$]	1.75	${}_{-0.06}{}^{}{}_{}{}^{+0.04}$
$a/R_{\star }$bbIf the planet’s orbit is not circular, this corresponds to the average planet-star-separation during transit divided by the stellar radius.	53.8	${}_{-0.6}{}^{}{}_{}{}^{+0.4}$
Impact parameter, $\left|b\right|$	0.11	${}_{-0.00}{}^{}{}_{}{}^{+0.09}$
Sky-plane inclination, $i_{\mathrm{sky}}$ ${[}^{\circ }]$	89.89	${}_{-0.10}{}^{}{}_{}{}^{+0.08}$
Mid-transit times	1331.2851	${}_{-0.0007}{}^{}{}_{}{}^{+0.0007}$
	1365.8245	${}_{-0.0007}{}^{}{}_{}{}^{+0.0007}$
	1400.3686	${}_{-0.0007}{}^{}{}_{}{}^{+0.0007}$
	1434.9227	${}_{-0.0007}{}^{}{}_{}{}^{+0.0007}$
	1469.4773	${}_{-0.0007}{}^{}{}_{}{}^{+0.0007}$
LCO	1469.4781	${}_{-0.0004}{}^{}{}_{}{}^{+0.0004}$
Hazelwood	1504.037	${}_{-0.002}{}^{}{}_{}{}^{+0.002}$
El Sauce	1538.5939	${}_{-0.0015}{}^{}{}_{}{}^{+0.0015}$
Minimum mutual inclination ${[}^{\circ }]$	$1.8$	${}_{-0.2}{}^{}{}_{}{}^{+0.2}$

Table 4. Table 4: Light Curve Parameters for the TOI-216 system

ParameteraaAs a summary statistic we report the mode and 68.3% confidence interval of the posterior distribution.	TESS		El Sauce		LCO		Hazelwood
Limb darkening coefficient, $q_1$	$0.33$	${}_{-0.09}{}^{}{}_{}{}^{+0.12}$	$0.5$	${}_{-0.2}{}^{}{}_{}{}^{+0.2}$	$0.52$	${}_{-0.12}{}^{}{}_{}{}^{+0.15}$	$0.50$	${}_{-0.15}{}^{}{}_{}{}^{+0.23}$
Limb darkening coefficient, $q_2$	$0.32$	${}_{-0.11}{}^{}{}_{}{}^{+0.14}$	$0.30$	${}_{-0.16}{}^{}{}_{}{}^{+0.26}$	$0.21$	${}_{-0.08}{}^{}{}_{}{}^{+0.08}$	$0.7$	${}_{-0.2}{}^{}{}_{}{}^{+0.2}$
Red noise, ${\sigma }_r$ [ppm]	3000	${}_{-900}{}^{}{}_{}{}^{+800}$	10000	${}_{-4000}{}^{}{}_{}{}^{+4000}$	1500	${}_{-1000}{}^{}{}_{}{}^{+1600}$	4000	${}_{-3000}{}^{}{}_{}{}^{+3000}$
White noise, ${\sigma }_w$ [ppm]	2367	${}_{-17}{}^{}{}_{}{}^{+17}$	3140	${}_{-80}{}^{}{}_{}{}^{+80}$	1060	${}_{-40}{}^{}{}_{}{}^{+40}$	2450	${}_{-190}{}^{}{}_{}{}^{+190}$

Table 5. Table 5: Planet Parameters for TOI-216b and TOI-216c Derived from TTVs

Parameter	Soln 1a,ba,bfootnotemark:		Soln 2a,ca,cfootnotemark:
$M_b$ ($M_{\mathrm{Jup}}$)	0.05	${}_{-0.03}{}^{}{}_{}{}^{+0.023}$	0.10	${}_{-0.02}{}^{}{}_{}{}^{+0.03}$
$M_b/M_c$	0.149	${}_{-0.012}{}^{}{}_{}{}^{+0.011}$	0.133	${}_{-0.010}{}^{}{}_{}{}^{+0.010}$
$e_b$	0.214	${}_{-0.048}{}^{}{}_{}{}^{+0.154}$	0.15	${}_{-0.03}{}^{}{}_{}{}^{+0.04}$
${\varpi }_b$ (deg.)	240	${}_{-30}{}^{}{}_{}{}^{+40}$	293	${}_{-10}{}^{}{}_{}{}^{+7}$
$M_c$ ($M_{\mathrm{Jup}}$)	0.26	${}_{-0.17}{}^{}{}_{}{}^{+0.14}$	0.57	${}_{-0.16}{}^{}{}_{}{}^{+0.21}$
$e_c$	0.06	${}_{-0.03}{}^{}{}_{}{}^{+0.11}$
${\varpi }_c$ (deg.)	-30	${}_{-60}{}^{}{}_{}{}^{+30}$
$\Delta \varpi$ (deg.)	-80	${}_{-30}{}^{}{}_{}{}^{+30}$
$2{\lambda }_c-{\lambda }_b-{\varpi }_c$ (deg).	-20	${}_{-30}{}^{}{}_{}{}^{+40}$
$2{\lambda }_c-{\lambda }_b-{\varpi }_b$ (deg).	60	${}_{-14}{}^{}{}_{}{}^{+11}$	41	${}_{-6}{}^{}{}_{}{}^{+7}$

Table 6. Table 6: Missed and future transit times

Solution 1a,ba,bfootnotemark:		Solution 2a,ca,cfootnotemark:
TOI-216b
1530.286	${}_{-0.004}{}^{}{}_{}{}^{+0.006}$	1530.295	${}_{-0.007}{}^{}{}_{}{}^{+0.011}$
1547.351	${}_{-0.007}{}^{}{}_{}{}^{+0.009}$	1547.363	${}_{-0.010}{}^{}{}_{}{}^{+0.013}$
1564.413	${}_{-0.010}{}^{}{}_{}{}^{+0.013}$	1564.430	${}_{-0.015}{}^{}{}_{}{}^{+0.020}$
1581.479	${}_{-0.014}{}^{}{}_{}{}^{+0.019}$	1581.50	${}_{-0.02}{}^{}{}_{}{}^{+0.02}$
1598.54	${}_{-0.02}{}^{}{}_{}{}^{+0.03}$	1598.58	${}_{-0.03}{}^{}{}_{}{}^{+0.03}$
1615.61	${}_{-0.02}{}^{}{}_{}{}^{+0.04}$	1615.65	${}_{-0.04}{}^{}{}_{}{}^{+0.04}$
TOI-216c
1573.09	${}_{-0.03}{}^{}{}_{}{}^{+0.03}$	1573.16	${}_{-0.03}{}^{}{}_{}{}^{+0.04}$
1607.63	${}_{-0.04}{}^{}{}_{}{}^{+0.04}$	1607.71	${}_{-0.04}{}^{}{}_{}{}^{+0.05}$
1642.18	${}_{-0.05}{}^{}{}_{}{}^{+0.04}$	1642.26	${}_{-0.04}{}^{}{}_{}{}^{+0.05}$
1676.72	${}_{-0.05}{}^{}{}_{}{}^{+0.05}$	1676.82	${}_{-0.05}{}^{}{}_{}{}^{+0.06}$
1711.26	${}_{-0.06}{}^{}{}_{}{}^{+0.05}$	1711.37	${}_{-0.05}{}^{}{}_{}{}^{+0.07}$
1745.81	${}_{-0.06}{}^{}{}_{}{}^{+0.06}$	1745.92	${}_{-0.06}{}^{}{}_{}{}^{+0.07}$

Full text

1 Introduction

Warm large exoplanets, giant planets with 10–100 day orbital periods, pose a major challenge to our understanding of how planets form and evolve. Origins hypotheses developed and fine-tuned to account for the more readily discovered hot Jupiters (orbital periods $<10$ days) and the far more abundant warm sub-Neptunes find it challenging to account for warm, large exoplanets’ occurrence rates, eccentricities, masses, and companion properties (e.g., Wu & Lithwick 2011; Beaugé & Nesvorný 2012; Petrovich 2015; Dawson et al. 2015a; Huang et al. 2016; see Section 4.3 of Dawson & Johnson 2018 for a review). Although rarer than smaller planets and more distant giants, warm, large exoplanets are an outcome of physical processes that likely sculpt many planetary systems.

Recently some have argued for two origins channels for warm, large exoplanets (e.g., Dawson & Murray-Clay 2013; Dong et al. 2014; Dawson et al. 2015b; Petrovich & Tremaine 2016): high eccentricity tidal migration, and a second channel that may involve disk migration and/or in situ formation. Under the hypothesis of high eccentricity tidal migration, warm, large exoplanets are planets caught in the act of migration: they began further from the star, were disturbed onto highly elliptical orbits, and are tidal circularizing to short orbital periods. However, a key piece of evidence supporting the second channel is the handful of warm, large exoplanets with nearby planets, which are incompatible with high eccentricity migration and are not on route to becoming hot Jupiters. Fig. 1 shows all confirmed systems with a warm, large exoplanet (mass greater than 0.25 $M_{\rm Jup}$ or radius greater than 8 Earth radii; period less than 100 days) and a companion with a $<100$ day orbital period. It is striking that most of these systems are in or near an orbital resonance, and almost all contain a known small planet on a $<10$ day orbital period, despite the low occurrence rate of such short period planets in general (e.g., Mulders et al. 2015). They also happen to be some of the most iconic, well-studied exoplanet systems, probably because large and/or massive planets with short orbital periods are most amenable to transit and radial velocity characterization. Discovering and characterizing more warm, large exoplanets with nearby planets could help shed light on the nature of this second channel.

The TESS pipeline (Jenkins et al., 2016; Twicken et al., 2018; Li et al., 2018) recently discovered a pair of warm, large planet candidates orbiting TOI-216. Like the other systems in Fig. 1, the putative planets are in or near an orbital resonance. Their proximity to resonance leads to detectable transit timing variations (TTVs). Based on expected TESS planet yields, Hadden et al. (2018) predicted that significant mass constraints from TTVs would be possible for of order five planets. Here we seek to validate and characterize the TOI-216 planet candidates and assess what additional follow-up is necessary to test theories for their origin. We characterize the host star in Section 2. In Section 3, we describe our analysis of the TESS data and extraction of planet parameters. We rule out most astrophysical false positive scenarios in Section 4. We constrain the system’s orbital architecture in Section 5 – including mutual inclination, TTVs, eccentricities, and additional transit signals – and the planets’ masses sufficiently to confirm the planets. We present our conclusions in Section 6.

2 Stellar Characteristics

TOI-216 is an 11.5 TESS apparent magnitude, main sequence K-dwarf. To better refine its parameters – particularly the metallicity – we obtained seven spectra of TOI-216 with the ANU 2.3m Echelle spectrograph over a period of 11 days in 2018 Nov. These observations were also made to broadly constrain the mass of the planets and to check for obvious astrophysical false positive scenarios, such as line blending due to background stars. The ANU 2.3m/Echelle is located at Siding Spring Observatory, Australia. The spectrograph has a spectral resolution of $\lambda/\Delta\lambda\equiv R=23000$, covering the wavelength region of 3900–6700 Å. Observations are bracketed by ThAr arc lamp exposures for wavelength calibration. Instrument stability issues limit the radial velocities to a typical precision of only $\sim 500\,\mathrm{m\,s}^{-1}$ for this facility. Stellar parameters for TOI-216 were derived using SpecMatch (Yee et al., 2017) on the ANU 2.3m/Echelle spectra, yielding atmospheric parameters of $T_{\mathrm{eff}}=5045\pm 110$ K, $\log g=4.53\pm 0.12$ dex, and $\mathrm{[Fe/H]}=-0.16\pm 0.09$ dex.

We use the approach described by Dawson et al. (2015b) to fit the observed stellar properties using the Takeda et al. (2007) and Dartmouth (Dotter et al., 2008) stellar evolution models. We perform an additional fit using the Dartmouth models to both the spectrum properties and the Gaia DR2 parallax and apparent $g$ magnitude (Gaia Collaboration et al., 2016, 2018). We find that the measured atmospheric parameters are consistent with a main-sequence K-dwarf and list the derived stellar mass, radius, and density in column 2 of Table 1. The resulting values are in agreement with the TESS Input Catalog (TIC; Stassun et al. 2018) but more precise. We choose to use the Dartmouth values hereafter because the posteriors extend to a lower mass ($M_{\star}<0.7M_{\odot}$) than covered by the Takeda et al. (2007) models and because they allow us to fit the Gaia DR2 parameters.

3 Light Curve Analysis

TOI-216 is located near the southern ecliptic pole, and is scheduled to be observed for 12 sectors of the first year of the TESS Primary Mission. This paper is based on data from Sectors 1, 2, 3, 4, 5, and 6 (2018 July 25 – 2019 January 7), during which TOI-216 was observed with CCD 1 on Camera 4, and from ground-based observatories.

3.1 Data from TESS Mission

We use the publicly available 2-min cadence data from the TESS Alerts, which is processed with the Science Processing Operations Center pipeline. The pipeline, a descendant of the Kepler mission pipeline based at the NASA Ames Research Center (Jenkins et al., 2002, 2010, 2016), analyzes target pixel postage stamps that are obtained for pre-selected target stars. For TOI-216, the short cadence pipeline detected two threshold crossing event at periods 34.54 days and 17.1 days with high signal-to-noise. The candidates were also detected by the long cadence MIT Quick Look Pipeline (Sha et al., 2019).

3.2 Ground-based Photometric Follow-up

We used the resources of the TESS Follow-up Observing Program (TFOP) Working Group (WG) Sub Group 1 (SG1)111https://tess.mit.edu/followup/ to collect seeing-limited time-series photometric follow-up of TOI-216. The transit depths of both TOI-216 planet candidates, as predicted by the TESS light curves, are deep enough to detect from the ground at high significance. Therefore our primary goal was to attempt to detect the transits using our higher spatial resolution ground-based imaging and a photometric aperture that is small enough to exclude the flux from known nearby stars that are bright enough to cause the TESS detected events. The secondary goal was to identify or rule out potential nearby eclipsing binaries (Section 4).

We used the TESS Transit Finder, which is a customized version of the Tapir software package (Jensen, 2013), to schedule photometric time-series follow-up observations. We initially scheduled observations for both planet candidates according to the public linear ephemerides derived from Sectors 1 and 2 TESS data. Our eight time-series follow-up observations are listed in Table 2. We used the AstroImageJ software package (Collins et al., 2017) for data reduction and aperture photometry for all of our follow-up photometric observations. The facilities used to collect the TOI-216 observations are: Las Cumbres Observatory (LCO) telescope network (Brown et al., 2013); Hazelwood Observatory; the Myers-T50 Telescope; and El Sauce Observatory. All LCO 1 m telescopes are equipped with the Sinistro camera, with a 4k x 4k pixel Fairchild back illuminated CCD and a 26.5 x 26.5 arcmin FOV. The LCO 0.4 m telescopes are mounted with an SBIG STX6303 2048 x 3072 pixels CCD with a 19 x 29 arcmin FOV. Hazelwood is a private observatory with an f/8 Planewave Instruments CDK12 0.32 m telescope and an SBIG STT3200 2.2K$\times$1.5K CCD, giving a $20\arcmin\times 13\arcmin$ field of view. The Myers-T50 is an f/6.8 PlaneWave Instruments CDK17 0.43 m Corrected Dall-Kirkham Astrograph telescope located at Siding Spring, Australia. The camera is a Finger Lakes Instruments (FLI) ProLine Series PL4710 - E2V, giving a $15\farcm 5\times 15\farcm 5$ field of view. El Sauce is a private observatory with a Planewave CDK14 0.36 m telescope on a MI500/750F fork mount. The camera is an SBIG STT1603-3 1.5K$\times$1.0K CCD, giving a $18\farcm 5\times 12\farcm 3$ field of view.

We observed five transits of TOI-216c at three epochs and confirmed that the transit events occur on target using follow-up apertures with radius $\sim 6\arcsec$. We conducted five TOI-216b observations at four transit epochs and ruled out the $\sim 4$ parts per thousand transit events at the public linear ephemeris. However, with the later addition of data from TESS sectors 3 and 4 to the TTV analysis, we determined that the large TTV signal caused the transit events to egress before our follow-up observations started. We then observed an out-of-transit sequence that occurred just prior to the newly determined transit ingress time to help constrain the TTV model (since the time of transit was not observable from our available facilities).

3.3 Light curve fits

We fit the transit light curves (Fig. 2) using the TAP software (Gazak et al., 2012), which implements Markov Chain Monte Carlo using the Mandel & Agol (2002) transit model and the Carter & Winn (2009) wavelet likelihood function, with the modifications described in Dawson et al. (2014). The results are summarized in Table 3. We use the presearch data conditioned (PDC) flux, which is corrected for systematic (e.g., instrumental) trends using cotrending basis vectors (Smith et al., 2012; Stumpe et al., 2014); the Carter & Winn (2009) wavelet likelihood function (which assumes frequency*-1* noise) with free parameters for the amplitude of the red and white noise; and a linear trend fit simultaneously to each transit light curve segment with other transit parameters. We assign each instrument (TESS, Hazelwood, LCO, El Sauce) its own set of limb darkening parameters because of the different wavebands. We use different noise parameters for TESS, Hazelwood, LCO, and El Sauce. We adopt uniform priors on the planet-to-star radius ratio ($R_{p}/R_{\star}$), the log of the light curve stellar density $\rho_{\rm circ}$ (i.e., equivalent the light curve parameter $d/R_{\star}$, where $d$ is the planet-star separation, converted to stellar density using the planet’s orbital period and assuming a circular orbit), the impact parameter $b$ (which can be either negative or positive; we report $|b|$), the mid transit time, the limb darkening coefficients $q_{1}$ and $q_{2}$ (Kipping, 2013), and the slope and intercept of each transit segment’s linear trend. For the Hazelwood, LCO, and El Sauce observations, we fit a linear trend to airmass instead of time.

The inner planet candidate’s transits are grazing, so the planet-to-star radius ratio $R_{p}/R_{\star}$ is not well-constrained. We impose a uniform prior from 0 to 0.17, with the upper limit corresponding to a radius of 0.13 solar radii. Fig. 3 shows the covariance between $R_{p}/R_{\star}$ and the light curve stellar density $\rho_{\rm circ}$ and impact parameter $b$. The larger the planet, the larger the impact parameter required to match the transit depth. The larger the impact parameter, the shorter the transit chord and the lower the light curve stellar density (which correlates with the transit speed) required to match the transit duration. Through its affect on $|b|$ and $\rho_{\rm circ}$, the upper limit on $R_{p}/R_{\star}$ affects our inference of the inner planet’s eccentricity and the mutual inclination between the planets; in Section 5, we will assess the sensitivity to this upper limit.

3.4 Search for additional transit signals

We ran the box car least squares algorithm on the residuals of the light curve after removing the transit signal of TOI-216b and TOI-216c. We used a duration of 2.5 hours, which corresponds to an impact parameter equal to planet c’s at an orbital period of 3 days. We did not find any signal with signal-to-noise larger than 7.3. Using per-point rms precision of 0.00233, this limit rules out any planets interior to TOI-216b with a radius larger than 2.18 $R_{\oplus}$ or planets with periods less than 3 days and radii larger than 1.17 $R_{\oplus}$. With future TESS data from 12 sectors in total, the detection threshold for all planets interior to TOI-216b will be lowered to 1.13 $R_{\oplus}$ planets.

4 Validation

Here we seek to validate the planet candidates by ruling out false positive scenarios using follow up observations and dynamical arguments. In Section 4.1, we consider and rule out unblended astrophysical false positive scenarios using radial velocity (RV) measurements. In 4.2, we consider and rule out most blended false positive scenarios using photometry. We summarize the results in Section 4.3.

4.1 Low precision radial-velocity follow up to rule out stellar companions to TOI-216

One or both transiting objects could be brown dwarf or stellar companions to TOI-216. The following astrophysical false positive scenarios can be tested through radial velocity follow up: TOI-216b and/or TOI-216c is a brown dwarf; TOI-216b (which has a poorly constrained transit depth) is an unblended stellar companion; or TOI-216b and/or TOI-216c is a blended stellar companion to TOI-216 with a background or bound star diluting the transit depth.

If both objects are transiting TOI-216, but one or both is of brown dwarf or stellar mass, the system would be unstable if the objects are not in resonance or, if in resonance, the mass of the secondary would cause large TTVs incompatible with those observed (Section 5.8). Furthermore, the brown dwarf scenario is less likely a priori. Grieves et al. (2017) find that the occurrence rate of brown dwarfs with orbital periods less than 300 days is about 0.56%, compared to 4.0% for planets $>0.3M_{\rm Jup}$ (Cumming et al., 2008).

We use radial velocity (RV) measurements to put mass limits on any companion to TOI-216. The spectra described in Section 2 shows no large radial velocity variations, with the measurements exhibiting a scatter of $470\,\mathrm{m\,s}^{-1}$. From these velocities, we derive the $3\sigma$ upper limit on the masses of the inner planet to be $\sim 18\,M_{J}$ and the outer planet to be $\sim 25\,M_{J}$. The upper limits rule out any scenario involving a stellar companion to TOI-216. The constraints also support our limit on $R_{p}/R_{\star}$ for the light curve fits for TOI-216b (Section 3) corresponding to 1.3 $R_{J}$ because radii only start to increase above $\sim 1R_{J}$ at around $60M_{J}$ (e.g., Hatzes & Rauer 2015, Fig. 2). The scenario in which one or both objects are brown dwarf companions is not ruled out by the RVs but will be ruled out by the TTVs in Section 5.8.

4.2 Photometry rules out most blended false positive scenarios

Analysis of systems with multiple transiting planet candidates from Kepler has shown that the transit-like events have a higher probability of being caused by bona fide planets (e.g., Lissauer et al., 2012) compared to single-planet candidate systems, lending credibility to the planetary nature of the transit-like events associated with TOI-216. However, the pixel scale of TESS is larger than Kepler’s (${21\arcsec}$ for TESS vs. ${4\arcsec}$ for Kepler) and the point spread function of TESS could be as large as ${1\arcmin}$, both of which increase the probability of contamination of the TESS aperture by an nearby eclipsing binary. For example, a deep eclipse in a nearby faint eclipsing binary might cause a shallow transit-like detection by TESS on the target star due to the dilutive effect of blending in the TESS aperture.

A scenario in which both TOI-216b and TOI-216c are orbiting the same background binary is ruled out by the TTVs (Section 5.8). One object could be a planet-mass companion to TOI-216 and the other a background binary. Alternatively, both objects could be background binaries.

From a single sector of TESS data, the one standard deviation centroid measurement uncertainty is $2\farcs 58$ for TOI-216b and $3\farcs 3$ for TOI-216c. TOI-216c would need to fully eclipse a star with Tmag 15.85 to cause the blend, and TOI-216b would need to fully eclipse a star with Tmag 17.5 to cause the blend. The brightest Gaia D2 object within $40\arcsec$ has Gaia $rp$ magnitude of 16.8 and is $3\farcs 768$ away and therefore is marginally compatible with a blend scenario for TOI-216b. The second brightest Gaia object within $40\arcsec$ has Gaia $rp$ magnitude of 17.94, which cannot cause either of the transit signals we see.

We use higher spatial resolution ground-based time-series imaging to attempt to detect the transit-like events on target and/or to identify or rule out potential nearby ecliping binaries out to ${2.5\arcmin}$ from TOI-216. The higher spatial resolution and smaller point spread function of the ground-based observations facilitates the use of much smaller photometric apertures compared to the TESS aperture to isolate a possible transit or eclipse signal to within a few arcseconds of the center of the follow-up aperture. From the ground, follow-up apertures exclude the flux of all known neighboring stars, except the two $\sim 4\arcsec$ Gaia DR2 neighbors. We collected observations of TOI-216c in both g′ and i′ filters (Section 3) and found no obvious filter-dependent transit depth, which strengthens the case for a planetary system.

4.3 Validation summary

In summary, we can rule out all astrophysical false positive scenarios with a couple exceptions. First, TOI-216b could be a blended binary orbiting the 16.8 $rp$ magnitude Gaia DR2 object, in which case it would need a 53% transit depth. Second, TOI-216b and/or TOI-216c could be a binary orbiting a star located at the same sky position as TOI-216, creating a blend not resolved by Gaia. However, we will show in Section 5.8 that the two transiting objects are fully compatible with causing each other’s TTVs and that the TTVs have concavity in opposite direction (i.e., one planet loses orbital energy as the other gains). This false positive scenario would require the extremely unlikely configuration in which both objects happen to have an orbital period ratio near 2:1, happen to have non-transiting companions in or near orbital resonance causing their TTVs, and the TTVs happen to have opposite sign. Therefore we consider the system to be validated.

5 Orbital Architecture

Here we explore the orbital architecture of the TOI-216 system through analysis of the transiting timing variations (TTVs), transit shape and duration, and limits on additional transiting planets.

5.1 TTV overview

Both candidates exhibit significant deviations from a linear transit time ephemeris (Fig. 4), evidence for their mutual gravitational perturbations. These transit timing variations (TTVs) occur on two timescales. The first is the synodic timescale, $\tau_{\rm syn}=P_{c}/(P_{c}/P_{b}-1)$, which is the interval of time between successive planetary conjunctions. The second – for planets near the 2:1 resonance – is the super-period222The super-period may be longer or shorter for planets in orbital resonance experiencing fast precession., $\tau_{\rm s-p}\approx|P_{c}/(2-P_{c}/P_{b})|$, the timescale over which the planets have their conjunctions at the same longitude; $\tau_{\rm s-p}$ depends on the proximity of the ratio of the orbital periods to 2.

The synodic TTV signal, known as the chopping effect because it produces a saw-tooth like pattern (see Deck & Agol 2015 and references therein), depends on the perturbing planet’s mass, which determines the strength of the kick at conjunction. To first order, the chopping effect does not depend on eccentricity.

The super-period TTV signal, known as the near-resonant effect (e.g., Lithwick et al. 2012), has a sinusoidal shape. The near-resonant effect generates a forced eccentricity for each planet, and the free eccentricity is an extra component that contributes to the total eccentricity. The near-resonant TTV amplitude depends on the perturbing planet’s mass and the free eccentricity of the transiting and perturbing planets. To first order, the ratio of near-resonant signal amplitudes depends only on the planets’ mass ratio (e.g., Lithwick et al. 2012’s Eqn. 14–15). Therefore TTVs covering a significant fraction of the super-period can provide a good estimate of the mass ratio.

For planets near resonance, the amplitude of the near-resonant effect is typically much larger than the amplitude of the chopping effect. Measuring the chopping and near-resonant signals for both transiting planets – assuming there are no additional planets in the system contributing significantly to the TTVs – would allow us to uniquely constrain their masses and eccentricities.

5.2 Evidence for free eccentricity

The phasing of the TTVs allows to diagnose at least one planet likely has significant free eccentricity. In Fig. 5 we plot the TTVs as a function of phase. The top panel shows the TTVs of the inner planet phased with $2(\lambda_{b}-\lambda_{c}$), where $\lambda$ is the mean longitude (Section 5.3). If the free eccentricities are zero, the TTVs should follow a sinusoid with no phase shift (Deck & Agol, 2015). The non-phase shifted sinusoid is inconsistent with the observed TTVs of TOI-216b, so we infer at least one planet has free eccentricity. [For the outer planet, no phase shift in $\lambda_{b}-\lambda_{c}$ (Fig. 5, row 2) is necessary.] We also follow Lithwick et al. (2012) and plot the TTVs phased to $2\lambda_{c}-\lambda_{b}$ (Fig. 5, row 3; equivalent to rows 1–2 because transit times are sampled at the planets’ orbital period) and find that again a phase shift is necessary to match the inner planet’s observed TTVs, indicating free eccentricity for one or both planets.

The chopping signal would appear as additional harmonics, i.e., $\lambda_{b}-\lambda_{c}$, $3(\lambda_{b}-\lambda_{c})$, etc. for TOI-216b and $2(\lambda_{b}-\lambda_{c})$ and $3(\lambda_{b}-\lambda_{c})$, etc. for TOI-216c. The fact that a sinusoid goes through the data points in Fig. 4 without these additional harmonics gives us a sense that the chopping signal will not be easily measured in this dataset. There will be a degeneracy between planet masses and free eccentricity.

5.3 A large range of best-fit planet masses

We fit the transit times using our N-body TTV integrator model (Dawson et al., 2014). Our model contains five parameters for each planet: the mass $M$, orbital period $P$, mean longitude at epoch $\lambda$, eccentricity $e$, and argument of periapse $\omega$. For each planet, we fix the sky plane inclination $i_{\rm sky}$ to the value in Table 3 and set the longitude of ascending node on the sky to $\Omega_{\rm sky}=0$. We use the conventional coordinate system where the $X-Y$ plane is the sky plane and the $Z$ axis points toward the observer. See Murray & Correia (2010) for a helpful pedagogical description of the orbital elements.

To explore the degeneracy between mass and eccentricity, we use the Levenberg-Marquardt alogrithm implemented in IDL mpfit (Markwardt, 2009) to minimize the $\chi^{2}$ on a grid of $(M_{c},e_{b})$. We report the total $\chi^{2}$ for eighteen transit times and ten free parameters, i.e., eight degrees of freedom. The resulting contour plot is shown in Fig. 6. The lowest $\chi^{2}$ fits, i.e., those with $13<\chi^{2}<18$, are possible for a range of outer planet masses ($M_{c}<3.0\,M_{\rm Jup}$). However, for small outer planet masses, a large range of inner planet eccentricities allow for a good fit, whereas a particular value of the eccentricity ($e_{b}\sim 0.13$ is necessary for larger planet masses. (See also discussions by Hadden & Lithwick 2017 and Migaszewski & Goździewski 2018.) Because there is so much more “real estate” in parameter space at low outer planet masses, a Markov Chain Monte Carlo (MCMC) will identify this type of solution as most probable. However, if we have a priori reason to suspect the outer planet is massive – like a large transit depth – and/or that free eccentricities are low, we could be misled.

Fig. 7 shows how other parameters correlate with the outer planet’s mass $M_{c}$. The mass ratio, $M_{b}/M_{c}$, of the planets is about 0.17 for $M_{c}<0.5M_{\rm Jup}$ and decreases for larger $M_{c}$. Solutions with $M_{c}<0.5M_{\rm Jup}$ have larger values for the eccentricity of planet c. (Note that the eccentricity plotted in Fig. 6 and 7 is the osculating eccentricity; we will explore how these solutions translate to free and forced eccentricities in Section 5.4.) The arguments of periapse $\omega_{b}$ and $\omega_{c}$ for planets b and c also correlate with planet c’s mass.

5.4 Longterm behavior of best-fit solutions

We integrate the $\chi^{2}<18$ solutions for $10^{6}$ days using mercury6 (Chambers et al., 1996) to assess the longer term behavior (Fig. 8). We find that resonant argument $2\lambda_{c}-\lambda_{b}-\varpi_{b}$ librates for the high $M_{c}$ ($M_{c}\gtrapprox 0.3\,M_{\rm Jup}$) solutions but not for the lower $M_{c}$ solutions. Larger $M_{c}$ solutions have lower free and forced eccentricities for both planets (Fig. 8, rows 2–3). Period ratios $P_{c}/P_{b}$ are wider of the 2:1 for the higher $M_{c}$ solutions. We extend the simulations to 10 Myr and find that all configurations remain stable over that interval.

5.5 Transit exclusion intervals

We use ground-based observations in which an ingress or egress for TOI-216b is excluded (Table 2) to check solutions. Before the TESS Sector 6 data were available, the exclusion interval on the December 26 observation ruled out some solutions. Almost all solutions based on Sector 1–6 are consistent with no ingress or egress during the intervals in Table 2.

5.6 Ruling out the lowest-mass solutions with the “photoeccentric” effect

The light curve stellar densities (Table 3) are similar to the true stellar density (Table 1), consistent with the planets being on nearly circular orbits. We follow Dawson & Johnson (2012) to estimate the candidates’ eccentricities from the light curve using the “photoeccentric effect,” but instead of applying the approximations appropriate for a grazing transit, we use the full Eqn. 15 from Kipping (2010). We find eccentricities that could be low for both candidates; their modes and 68.3% confidence intervals a are: $e_{b}=0.20_{-0.06}^{+0.48},e_{c}=0.025_{-0.004}^{+0.490}$ (Fig. 9). The medians and their 68.3% confidence intervals are $e_{b}=0.30_{-0.16}^{+0.38},e_{c}=0.10_{-0.08}^{+0.41}$. High eccentricities are not ruled out, e.g., the posterior probability of $e>0.5$ is 28% for TOI-216b and 16% TOI-216 c. The posterior probability of an eccentricity less than 0.01 is 0.7% for TOI-216b and 8% for TOI-216c.

The constraints on the eccentricity from the light curve allow us to rule out the lowest-mass solutions (Fig. 10). These solutions – which correspond to an eccentric TOI-216c with its apoapse near our line of sight – would produce a transit duration that is too long. Some higher-mass solutions that correspond to an eccentric TOI-216c with its periapse near our line of sight are also ruled out.

5.7 MCMC fits

Following Dawson et al. (2014), we derive posteriors for the parameters using Markov Chain Monte Carlo with the Metropolis-Hastings algorithm. We incorporate the transit exclusion intervals and light curve stellar density (i.e., combining the $\rho_{\rm circ}$ posterior from the light curve and $\rho_{\star}$ posterior from the Dartmouth models) into the MCMC. Instead of including the orbital period and mean longitude at epoch as parameters in the MCMC, we optimize them at each jump using the Levenberg-Mardquardt algorithm. We visually inspect each parameter for convergence.

We perform two fits with different priors to explore both ends of the parameter degeneracy evident in the grid of outer mass vs. inner eccentricity (Fig. 6). The first solution (Table 5, column 1) imposes uniform priors on eccentricities and log uniform priors on mass (i.e., priors that are uniform in log space); the second (Table 5, column 2) imposes uniform priors on mass and sets $e_{c}=0$ (which we found to yield indistinguishable results from an eccentricity prior that is uniform in log space). All other fitted parameters (orbital period, mean longitude, argument of periapse) have uniform priors. The uniform priors on mass favors the higher-mass solutions seen in Fig. 6–8, whereas the log uniform prior on mass favor the lower-mass solutions.

Because the results are so prior-dependent (every parameter in Table 5 differs significantly between the two solutions except TOI-216b’s eccentricity of $\sim 0.2$), we do not recommend currently adopting either solution. Instead, the MCMC approach is a way to formally separate the two types of solutions seen in the grid search and to incorporate the light curve stellar densitities and transit exclusion windows into the likelihood function.

5.8 Mass-radius

We plot the two solutions on a mass-radius plot in Fig. 11. TOI-216c’s radius is comparable to other known exoplanets for both mass solutions. The same is true for TOI-216b if its radius is close to the lower limit derived from its grazing transits. However, if its radius is somewhat larger than the lower limit, the lower-mass solution would correspond to a very low density.

5.9 Predictions for future transits

In Table 6, we tabulate the predicted times for missed and future transits. For the inner planet, the predictions of the two solutions overlap within one standard deviation for each transit. However, the outer planet’s transits differ between the solutions, so the next few sectors of TESS data may help distinguish between them.

5.10 Mutual inclination

A larger impact parameter for an inner planet than an outer planet points to at least a small mutual inclination between their orbits. The difference in the TOI-216b and TOI-216c’s sky plane inclination (Table 3) corresponds to a mutual inclination of at least $1^{\circ}.90^{+0.15}_{-0.34}$ (mode; the median is $1^{\circ}.8^{+0.2}_{-0.3}$). This value is a minimum because we do not know the component of the mutual inclination parallel to the sky plane. Future observations of transit duration variations – and depth changes for the grazing transit – may allow for constraints on the full 3D orbital architecture.

5.11 Comparison to other work

While this manuscript was in preparation, we learned of a submitted paper by Kipping et al. (2019) on this system. We conducted the work here independently. After submitting this manuscript and revising in response to the referee report, we read Kipping et al. (2019) study in order to compare our results. Our solutions are generally consistent. We infer a larger range of possible masses and eccentricities. We find a smaller radius for the outer planet due to our different stellar parameters derived from ground based spectroscopy and a larger range of possible radii and impact parameters for the inner planet. Ground-based transits aided our work by extending the TTV baseline and filling in transit times that were missed by TESS.

5.12 Summary

From the TTVs alone, we end up with solutions that occupy two qualitatively different parts of parameter space. The first corresponds to a sub-Saturn-mass planet and Neptune-mass planet with larger free eccentricities, period ratio near 2.00, and near but not in orbital resonance. The second corresponds to a Jupiter accompanied by a sub-Saturn with smaller free eccentricities, period ratios near 2.02, and librating in orbital resonance. Although the masses are not precisely constrained due to the degeneracy with eccentricity, we narrow the range of possible masses sufficiently to consider these candidates now confirmed as planets.

Although we cannot yet rule out the former solution, the latter solution has several appealing features. The period ratio falls outside the observed gap among Kepler multis (Fabrycky et al., 2014). The lower free eccentricities and libration of the resonant argument are suggestive of a dissipative process, such as disk migration, capturing the planets into resonance so that we observe them near a 2:1 period ratio. The masses are more typical of the observed radii (Fig. 11).

6 Discussion

TOI-216 is a system of two known transiting candidates in or near a 2:1 orbital resonance with accuracy-to-minutes constraints on their mid-transit times. Unlike most333See Dawson et al. 2014 for an example of a Kepler warm Jupiter with ground-based mid-transit times. Kepler systems, the 12.393 $V$ magnitude star is sufficiently bright for ground-based follow up to play an important role in supplying additional transits and transit exclusion intervals. From the phases of the TTVs, we identified that the pair contains significant free eccentricity that leads to degeneracy between eccentricities and masses. We ruled out the lowest-mass solutions using the “photoeccentric” effect and the highest-mass solutions using transit exclusion intervals from missed ground-based transits. Their mutual inclination may be modest (minimum $1^{\circ}.90^{+0.15}_{-0.34}$) but the component parallel to the sky plane is unknown. We identified two families of solutions. One solution family corresponds to lower masses (a sub-Saturn-mass planet and Neptune-mass planet), larger eccentricities, period ratio near 2, planets near but not in resonance, and puffy radii. The other corresponds to larger masses (Jupiter-mass planet and sub-Saturn-mass planet), lower eccentricities, a period ratio of 2.02, masses typical of the planets’ sizes, and orbital mean motion resonant libration. We prefer the second family of solutions but cannot yet rule out the first.

6.1 Formation and evolution

TOI-216 joins the population of systems featuring warm, large exoplanets that could not have achieved their close-in orbits through high eccentricity tidal migration (Fig. 1). They may have formed at or near their current locations (e.g., Huang et al. 2016), or formed at wider separations and migrated in (e.g., Lee & Peale 2002). Both scenarios could lead to planets in or near resonance (e.g., Dong & Dawson 2016; MacDonald & Dawson 2018). The in situ scenario would require the planets to coincidentally form with a period ratio close to 2, but in situ formation sculpted by stability can produce ratios near this value (e.g., Dawson et al. 2016). For the lowest-mass solutions, formation beyond the snow line may be necessary to account for the large radii (Lee & Chiang, 2016).

The planets have at least small and possibly moderate free eccentricities and mutual inclination. The free eccentricities and inclinations might result from dynamical interactions with other undetected planets in the system. For the higher-mass/low eccentricity solution, the eccentricities/inclinations are small enough to be consistent with self-stirring (e.g., Petrovich et al. 2014) by Neptune-mass or larger planets. The free eccentricities could even be generated by the gas disk (e.g., Duffell & Chiang 2015). However, the free eccentricities in the lower-mass solution would require nearby, undetected giant planets to accompany the observed sub-Saturn-mass planet and Neptune-mass planet pair.

Among the eleven systems featuring a warm, large exoplanet with companions with $<100$ day orbital period (Fig. 1), only TOI-216 and Kepler-30 lack a detected small, short period planet (Section 3.4). Whatever formation and migration scenario led to the short period planets in the other systems may not have operated here, or the planet may have been lost through stellar collision or tidal disruption. If present but non-transiting, such a planet would need to be mutually inclined to the rest of the system (for example, a non-transiting 3 day TOI-216 d would need to be inclined by $5^{\circ}$ with respect to TOI-216 c). The same stirring environment that led to free eccentricities could also have generated a mutual inclination for this interior planet. (Of course, it may be that no planet formed or migrated interior to TOI-216 b.) More generally, the mutual inclination between b and c makes it plausible that there are non-transiting planets in the system.

6.2 Future observations

Future TESS sectors will allow for additional transit timing measurements. As shown in Section 5.9, distinguishing between the two families of solutions may be possible with additional transits of the outer planets. Moreover, we can likely distinguish between the two families of solutions by measuring the masses through RV follow up: the radial velocity amplitudes are $\sim 53$ m/s and $\sim$ 2015 m/s for planet b and c respectively in Solution 1 (Table 5) and $\sim$ 10 m/s and $\sim$ 670 m/s for planet b and c respectively in Solution 2 (Table 5). We caution that because of the planets’ period ratio and mass ordering, the RV signal alone is subject to significant degeneracy between the inner planet’s mass and the outer planet’s eccentricity (Anglada-Escudé et al., 2010). Combining TTVs and RVs can break this degeneracy.

Unfortunately TOI-216 does not fall within the observable part of the sky for CHEOPS. Other space-based follow up possibilities, particularly to detect a change in transit depth/impact parameter for the inner planet due to its precession, include Spitzer.

We expect ground-based observations to play an essential role in follow up of TOI-216. As demonstrated here, ground-based observations can provide accurate and precise transit times for this bright star with two large transiting planets. For the larger planet in particular, ground-based transits can yield transit times that are more precise than from TESS data (e.g., the transit observed by LCO in Table 3). We can identify in advance which transit epoch(s) would be most valuable for distinguishing among models (Goldberg et al., 2018). Ground-based transits will allow for a long baseline of observations for better constraining the planets’ masses and eccentricities and possibly even detect precession of the planets’ orbits.

We thank the TESS Mission team and follow up working group for the valuable dataset. We acknowledge the use of public TESS Alert data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center. This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST). This research has made use of the Exoplanet Follow-up Observation Program website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This work has made use of observations from the Las Cumbres Observatory network. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products. We gratefully acknowledge support by NASA XRP NNX16AB50G and NASA TESS GO 80NSSC18K1695. The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. T.D. acknowledges support from MIT’s Kavli Institute as a Kavli postdoctoral fellow. K.H. acknowledges support from STFC grant ST/R000824/1. MŽ acknowledges funding from the Australian Research Council (grant DP170102233). We thank Samuel Hadden and Sarah Morrison for helpful discussions. We thank the referee for a helpful report that improved the clarity of the paper.