TOI-216b and TOI-216 c: Two Warm, Large Exoplanets in or Slightly Wide of the 2:1 Orbital Resonance

We use the publicly available 2 minute cadence data from the TESS Alerts, which is processed with the Science Processing Operations Center pipeline. We downloaded the publicly available data from the Mikulski Archive for Space Telescopes (MAST). 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 events at periods 34.54 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).

We used the resources of the TESS Follow-up Observing Program (TFOP) Working Group Sub Group 1 (SG1)²⁴ 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). All photometric time series are publicly available on the Exoplanet Follow-up Observing Program for TESS (ExoFOP-TESS) website.

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 the Las Cumbres Observatory (LCO) telescope network (Brown et al. 2013), the Hazelwood Observatory, the Myers-T50 Telescope, and the El Sauce Observatory. All LCO 1 m telescopes are equipped with the Sinistro camera, with a 4k × 4k pixel Fairchild back illuminated CCD and a 26.5 × 26.5 arcmin FOV. The LCO 0.4 m telescopes are mounted with an SBIG STX6303 2048 × 3072 pixels CCD with a 19 × 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 × 1.5K CCD, giving a 20' × 13' 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 155 × 155 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 × 1.0K CCD, giving a 185 × 123 field of view.

Table 2.  Observation Log

TOI-216	Date	Telescopea	Filter	ExpT	Exp	Dur.	Transit	Ap. Radius	FWHM
	(UTC)			(s)	(N)	(minutes)	Expected Coverageb	(arcsec)	(arcsec)
b	2018 Nov 22b	LCO-SSO-0.4	i'	90	54	100	Ingress+30%	8.5	7.5
	2018 Dec 09b	Myers-T50	Lum	60	200	240	full	8.3	4.6
	2018 Dec 26b	LCO-SSO-1.0	i'	30	85	99	Ingress+25%	7.0	2.8
	2019 Jan 29b	LCO-SAAO-1.0	${{\rm{r}}}^{{\prime} }$	100	97	225	Full	9.3	2.4
	2019 Jan 29b	LCO-SAAO-1.0	i'	25	181	198	Full	5.8	2.2
	2019 Feb 15	LCO-SSO-1.0	Zs	60	160	236	Out-of-transit	4.7	2.0
c	2018 Dec 16	LCO-SAAO-1.0	i'	90	75	180	Egress+60%	5.8	2.5
	2018 Dec 16	LCO-SAAO-1.0	i'	39	331	450	Full	5.8	2.1
	2019 Jan 20	Hazelwood-0.3	g'	240	101	449	Egress+70%	5.5	3.2
	2019 Feb 23	LCO-SAAO-1.0	Zs	60	148	212	Out-of-transit	6.2	2.5
	2019 Feb 24	LCO-CTIO-1.0	Zs	60	150	213	In-transit	6.2	2.5
	2019 Feb 24	El Sauce-0.36	Rc	30	514	303	Egress+90%	5.9	3.7
	2019 Feb 24	LCO-SSO-1.0	Zs	60	81	117	Out-of-transit	6.2	2.5

Notes.

^aTelescopes: 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). ^bObservations 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.

Download table as:  ASCIITypeset image

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 ∼6''. We conducted five TOI-216b observations at four transit epochs and ruled out the ∼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).

We fit the transit light curves (Figure 2) using the TAP software (Gazak et al. 2012), which implements Markov Chain Monte Carlo (MCMC) using the Mandel & Agol (2002) transit model and the Carter & Winn (2009) wavelet likelihood function, with the modifications described in Dawson et al. (2014) and in Table 4. The results are summarized in Table 3. We use the presearch data conditioned 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⁻¹ 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_⋆), the log of the light-curve stellar density ${\rho }_{\mathrm{circ}}$ (i.e., equivalent the light-curve parameter d/R_⋆, 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 midtransit time, the limb darkening coefficients q₁ and q₂ (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.

Zoom In Zoom Out Reset image size

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

Parameter	Valuea
TOI-216b
Planet-to-star radius ratio, Rp/R⋆	0.11	${}_{-0.02}^{+0.04}$
Planet radius, Rp [R⊕]	8.6	${}_{-1.9}^{+2.9}$
Light-curve stellar density, ${\rho }_{\mathrm{circ}}$ [ρ⊙]	1.13	${}_{-0.19}^{+0.29}$
a/R⋆b	29.1	${}_{-1.8}^{+2.3}$
Impact parameter, $| b|$	0.99	${}_{-0.04}^{+0.05}$
Sky-plane inclination,  isky [°]	88.0	${}_{-0.2}^{+0.2}$
Midtransit 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, Rp/R⋆	0.1236	${}_{-0.0008}^{+0.0008}$
Planet radius, Rp [R⊕]	10.2	${}_{-0.2}^{+0.2}$
Light curves stellar density, ${\rho }_{\mathrm{circ}}$ [ρ⊙]	1.75	${}_{-0.06}^{+0.04}$
a/R⋆b	53.8	${}_{-0.6}^{+0.4}$
Impact parameter, $| b|$	0.11	${}_{-0.00}^{+0.09}$
Sky-plane inclination, isky [°]	89.89	${}_{-0.10}^{+0.08}$
Midtransit 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 [°]	1.8	${}_{-0.2}^{+0.2}$

Notes.

^aAs a summary statistic we report the median and 68.3% confidence interval of the posterior distribution. ^bIf the planet's orbit is not circular, this corresponds to the average planet–star-separation during transit divided by the stellar radius.

Download table as:  ASCIITypeset image

The inner planet candidate's transits are grazing, so the planet-to-star radius ratio R_p/R_⋆ 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. Figure 3 shows the covariance between R_p/R_⋆ and the light-curve stellar density ${\rho }_{\mathrm{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 }_{\mathrm{circ}}$, the upper limit on R_p/R_⋆ 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.

Zoom In Zoom Out Reset image size

Table 4.  Light-curve Parameters for the TOI-216 System

Parametera	TESS		El Sauce		LCO		Hazelwood
Limb darkening coefficient, q1	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, q2	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, σr [ppm]	3000	${}_{-900}^{+800}$	10000	${}_{-4000}^{+4000}$	1500	${}_{-1000}^{+1600}$	4000	${}_{-3000}^{+3000}$
White noise, σw [ppm]	2367	${}_{-17}^{+17}$	3140	${}_{-80}^{+80}$	1060	${}_{-40}^{+40}$	2450	${}_{-190}^{+190}$

Note.

^aAs a summary statistic we report the mode and 68.3% confidence interval of the posterior distribution.

Download table as:  ASCIITypeset image

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 hr, 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_⊕ or planets with periods less than 3 days and radii larger than 1.17 R_⊕. 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_⊕ planets.

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 RV 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_Jup (Cumming et al. 2008).

We use RV measurements to put mass limits on any companion to TOI-216. The spectra described in Section 2 show no large RV variations, with the measurements exhibiting a scatter of $470\,{\rm{m}}\,{{\rm{s}}}^{-1}$. From these velocities, we derive the 3σ upper limit on the masses of the inner planet to be ∼18 M_J and the outer planet to be ∼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_⋆ for the light-curve fits for TOI-216b (Section 3) corresponding to 1.3 R_J because radii only start to increase above ∼1R_J at around 60M_J (e.g., Hatzes & Rauer 2015, Figure 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.

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'' for TESS versus 4'' for Kepler) and the point-spread function of TESS could be as large as 1', both of which increase the probability of contamination of the TESS aperture by a 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 258 for TOI-216b and 33 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'' has Gaia rp magnitude of 16.8 and is 3768 away and therefore is marginally compatible with a blend scenario for TOI-216b. The second brightest Gaia object within 40'' has a 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 eclipsing binaries out to 25 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 ∼4'' 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.

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 nontransiting 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.

Both candidates exhibit significant deviations from a linear transit time ephemeris (Figure 4), evidence for their mutual gravitational perturbations. These transit timing variations (TTVs) occur on two timescales. The first is the synodic timescale, ${\tau }_{\mathrm{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-period,²⁵ ${\tau }_{{\rm{s}}-{\rm{p}}}\approx | {P}_{c}/(2-{P}_{c}/{P}_{b})|$, the timescale over which the planets have their conjunctions at the same longitude; τ_s–p depends on the proximity of the ratio of the orbital periods to 2.

Zoom In Zoom Out Reset image size

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 Equations (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.

The phasing of the TTVs allows us to diagnose that at least one planet likely has significant free eccentricity. In Figure 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 λ 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 nonphase shifted sinusoid is inconsistent with the observed TTVs of TOI-216b, so we infer that at least one planet has free eccentricity. (For the outer planet, no phase shift in ${\lambda }_{b}-{\lambda }_{c}$ (Figure 5, row 2) is necessary.) We also follow Lithwick et al. (2012) and plot the TTVs phased to $2{\lambda }_{c}-{\lambda }_{b}$ (Figure 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.

Zoom In Zoom Out Reset image size

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})$, $3({\lambda }_{b}-{\lambda }_{c})$, etc. for TOI-216c. The fact that a sinusoid goes through the data points in Figure 4 without these additional harmonics gives us a sense that the chopping signal will not be easily measured in this data set. There will be a degeneracy between planet masses and free eccentricity.

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 λ, eccentricity e, and argument of periapse ω. For each planet, we fix the sky-plane inclination i_sky to the value in Table 3 and set the longitude of ascending node on the sky to Ω_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 χ² on a grid of $({M}_{c},{e}_{b})$. We report the total χ² for 18 transit times and 10 free parameters, i.e., eight degrees of freedom. The resulting contour plot is shown in Figure 6. The lowest χ² fits, i.e., those with 13 < χ² < 18, are possible for a range of outer planet masses (${M}_{c}\lt 3.0$ M_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 the 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, an 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.

Zoom In Zoom Out Reset image size

Figure 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}\lt 0.5{M}_{\mathrm{Jup}}$ and decreases for larger ${M}_{c}$. Solutions with ${M}_{c}\lt 0.5{M}_{\mathrm{Jup}}$ have larger values for the eccentricity of planet c. (Note that the eccentricity plotted in Figures 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 ω_b and ω_c for planets b and c also correlate with planet c's mass.

Zoom In Zoom Out Reset image size

We integrate the χ² < 18 solutions for 10⁶ days using mercury6 (Chambers et al. 1996) to assess the longer term behavior (Figure 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}_{\mathrm{Jup}}$) solutions but not for the lower ${M}_{c}$ solutions. Larger ${M}_{c}$ solutions have lower free and forced eccentricities for both planets (Figure 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.

Zoom In Zoom Out Reset image size

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.

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 non-grazing transit, we use the full Equation (30) from Kipping (2010). We find eccentricities that could be low for both candidates; their modes and 68.3% confidence intervals are ${e}_{b}={0.20}_{-0.06}^{+0.48},{e}_{c}={0.025}_{-0.004}^{+0.490}$ (Figure 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.

Zoom In Zoom Out Reset image size

The constraints on the eccentricity from the light curve allow us to rule out the lowest-mass solutions (Figure 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.

Zoom In Zoom Out Reset image size

Following Dawson et al. (2014), we derive posteriors for the parameters using MCMC with the Metropolis–Hastings algorithm. We incorporate the transit exclusion intervals and light-curve stellar density (i.e., combining the ${\rho }_{\mathrm{circ}}$ posterior from the light curve and ρ_⋆ 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 versus inner eccentricity (Figure 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 favor the higher-mass solutions seen in Figures 6–8, whereas the log uniform prior on mass favors the lower-mass solutions.

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

Parameter	Soln 1a,b		Soln 2a,c
${M}_{b}$ (MJup)	0.05	${}_{-0.03}^{+0.023}$	0.10	${}_{-0.02}^{+0.03}$
Mb/Mc	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}$
Mc (MJup)	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}$	⋯	⋯
Δϖ (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}$

Notes.

^aAs a summary statistic we report the median and 68.3% confidence interval of the posterior distribution. ^bUniform prior on eccentricity and log uniform prior on mass. ^ce_c = 0 and uniform prior on mass.

Download table as:  ASCIITypeset image

Because the results are so prior-dependent (every parameter in Table 5 differs significantly between the two solutions except TOI-216b's eccentricity of ∼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 densities and transit exclusion windows into the likelihood function.

We plot the two solutions on a mass–radius plot in Figure 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.

Zoom In Zoom Out Reset image size

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.

Table 6.  Missed and Future Transit Times

Solution 1a,b		Solution 2a,c
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}$

Notes.

^aAs a summary statistic we report the median and 68.3% confidence interval of the posterior distribution. ^bUniform prior on eccentricity and log uniform prior on mass. ^cLog uniform prior on eccentricity and uniform prior on mass.

Download table as:  ASCIITypeset image

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\buildrel{\circ}\over{.} {90}_{-0.34}^{+0.15}$ (mode; the median is $1\buildrel{\circ}\over{.} {8}_{-0.3}^{+0.2}$). 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.

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.'s (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.

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 and period ratios near 2.00 that are near but not in orbital resonance. The second corresponds to a Jupiter accompanied by a sub-Saturn with smaller free eccentricities and period ratios near 2.02 that are 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 (Figure 11).

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 (Figure 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 11 systems featuring a warm, large exoplanet with companions with <100 day orbital periods (Figure 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 nontransiting, such a planet would need to be mutually inclined to the rest of the system (for example, a nontransiting 3 day TOI-216 d would need to be inclined by 5° 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 nontransiting planets in the system.

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 RV amplitudes are ∼53 m s⁻¹ and ∼2015 m s⁻¹ for planets b and c, respectively, in Solution 1 (Table 5) and ∼10 m s⁻¹ and ∼670 m s⁻¹ for planets 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. 2019). 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.