Precise transit and radial-velocity characterization of a resonant pair: a warm Jupiter TOI-216c and eccentric warm Neptune TOI-216b
Rebekah I. Dawson, Chelsea X. Huang, Rafael Brahm, Karen A. Collins, Melissa J. Hobson, Andrés Jordán, Jiayin Dong, Judith Korth, Trifon Trifonov, Lyu Abe, Abdelkrim Agabi, Ivan Bruni, R. Paul Butler, Mauro Barbieri, Kevin I. Collins, Dennis M. Conti, Jeffrey D. Crane, Nicolas Crouzet, Georgina Dransfield, Phil Evans, Néstor Espinoza, Tianjun Gan, Tristan Guillot, Thomas Henning, Jack J. Lissauer, Eric L. N. Jensen, Wenceslas Marie Sainte, Djamel Mékarnia, Gordon Myers, Sangeetha Nandakumar, Howard M. Relles, Paula Sarkis, Pascal Torres, Stephen Shectman, François-Xavier Schmider, Avi Shporer, Chris Stockdale, Johanna Teske, Amaury H.M.J. Triaud, Sharon Xuesong Wang, Carl Ziegler, G. Ricker, R. Vanderspek, David W. Latham, S. Seager, J. Winn, Jon M. Jenkins, L. G. Bouma, Jennifer A. Burt, David Charbonneau, Alan M. Levine, Scott McDermott, Brian McLean, Mark E. Rose, Andrew Vanderburg, Bill Wohler
aa r X i v : . [ a s t r o - ph . E P ] F e b Draft version February 16, 2021
Typeset using L A TEX twocolumn style in AASTeX63
Precise transit and radial-velocity characterization of a resonant pair: a warm Jupiter TOI-216c andeccentric warm Neptune TOI-216b
Rebekah I. Dawson, Chelsea X. Huang,
2, 3
Rafael Brahm,
4, 5
Karen A. Collins, Melissa J. Hobson,
5, 7
Andr´es Jord´an,
4, 5
Jiayin Dong,
1, 8
Judith Korth, Trifon Trifonov, Lyu Abe, Abdelkrim Agabi, Ivan Bruni, R. Paul Butler, Mauro Barbieri, Kevin I. Collins, Dennis M. Conti, Jeffrey D. Crane, Nicolas Crouzet, Georgina Dransfield, Phil Evans, N´estor Espinoza, Tianjun Gan, Tristan Guillot, Thomas Henning, Jack J. Lissauer, Eric L. N. Jensen, Wenceslas Marie Sainte, Djamel M´ekarnia, Gordon Myers, Sangeetha Nandakumar, Howard M. Relles, Paula Sarkis, Pascal Torres,
Stephen Shectman, Franc¸ois-Xavier Schmider, Avi Shporer, Chris Stockdale, Johanna Teske,
Amaury H.M.J. Triaud, Sharon Xuesong Wang, Carl Ziegler, G. Ricker, R. Vanderspek, David W. Latham, S. Seager,
2, 30,31
J. Winn, Jon M. Jenkins, L. G. Bouma, Jennifer A. Burt, David Charbonneau, Alan M. Levine, Scott McDermott, Brian McLean, Mark E. Rose, Andrew Vanderburg, and Bill Wohler Department of Astronomy & Astrophysics, Center for Exoplanets and Habitable Worlds, The Pennsylvania State University, UniversityPark, PA 16802, USA Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA02139, USA Juan Carlos Torres Fellow Facultad de Ingenier´ıa y Ciencias, Universidad Adolfo Ib´a˜nez, Av. Diagonal las Torres 2640, Pe˜nalol´en, Santiago, Chile Millennium Institute for Astrophysics, Chile Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138, USA Instituto de Astrof´ısica, Facultad de F´ısica, Pontificia Universidad Cat´olica de Chile Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA Rheinisches Institut f¨ur Umweltforschung, Abteilung Planetenforschung an der Universit¨at zu K¨oln, Universit¨at zu K¨oln, Aachenerstraße209, 50931 K¨oln, Germany Max Planck Institute for Astronomy, Koenigstuhl 17, D-69117 Heidelberg, Germany Universit´e Cˆote d’Azur, Observatoire de la Cˆote d’Azur, CNRS, Laboratoire Lagrange Bd de l’Observatoire, CS 34229, 06304 Nicecedex 4, France INAF OAS, Osservatorio di Astrofisica e Scienza dello Spazio di Bologna via Piero Gobetti 93/3, Bologna Carnegie Institution for Science, Earth & Planets Laboratory, 5241 Broad Branch Road NW, Washington DC 20015, USA INCT, Universidad de Atacama, calle Copayapu 485, Copiap´o, Atacama, Chile George Mason University, 4400 University Drive, Fairfax, VA, 22030 USA American Association of Variable Star Observers, 49 Bay State Road, Cambridge, MA 02138, USA Observatories of the Carnegie Institution for Science, 813 Santa Barbara Street, Pasadena, CA 91101 European Space Agency (ESA), European Space Research and Technology Centre (ESTEC), Keplerlaan 1, 2201 AZ Noordwijk, TheNetherlands University of Birmingham, School of Physics & Astronomy, Edgbaston, Birmingham B15 2TT, United Kingdom El Sauce Observatory, Coquimbo Province, Chile Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA Department of Astronomy, Tsinghua University, Beijing 100084, China NASA Ames Research Center, Moffett Field, CA 94035 Dept. of Physics & Astronomy, Swarthmore College, Swarthmore PA 19081, USA Institut Paul Emile Victor, Concordia Station, Antarctica AAVSO, 5 Inverness Way, Hillsborough, CA 94010, USA Hazelwood Observatory, Australia NASA Hubble Fellow Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, Ontario M5S 3H4, Canada Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Department of Aeronautics and Astronautics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Corresponding author: Rebekah I. [email protected]
Dawson et al. Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08540, USA Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove drive, Pasadena CA 91109, USA Proto-Logic LLC, 1718 Euclid Street NW, Washington, DC 20009, USA Mikulski Archive for Space Telescopes, Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA Department of Astronomy, The University of Texas at Austin, Austin, TX 78712, USA NASA Sagan Fellow SETI Institute, Mountain View, CA 94043, USA (Received July 31, 2020; Revised December 12, 2020; Accepted January 4, 2021)
Submitted to AJABSTRACTTOI-216 hosts a pair of warm, large exoplanets discovered by the
TESS
Mission. These planets werefound to be in or near the 2:1 resonance, and both of them exhibit transit timing variations (TTV).Precise characterization of the planets’ masses and radii, orbital properties, and resonant behaviorcan test theories for the origins of planets orbiting close to their stars. Previous characterization ofthe system using the first six sectors of
TESS data suffered from a degeneracy between planet massand orbital eccentricity. Radial velocity measurements using HARPS, FEROS, and PFS break thatdegeneracy, and an expanded TTV baseline from
TESS and an ongoing ground-based transit observingcampaign increase the precision of the mass and eccentricity measurements. We determine that TOI-216c is a warm Jupiter, TOI-216b is an eccentric warm Neptune, and that they librate in the 2:1resonance with a moderate libration amplitude of 60 +2 − degrees; small but significant free eccentricityof 0 . +0 . − . for TOI-216b; and small but significant mutual inclination of 1.2–3.9 ◦ (95% confidenceinterval). The libration amplitude, free eccentricity, and mutual inclination imply a disturbance ofTOI-216b before or after resonance capture, perhaps by an undetected third planet. INTRODUCTIONWarm Jupiters – defined here as giant planets with10–100 day orbital periods – in systems with othernearby planets are an important population for theinvestigation of giant planets close to their stars.High eccentricity tidal migration – a strong contenderfor the origins of many close-in giant planets (seeDawson & Johnson 2018 for a review) is unlikely tohave been at work in these systems (e.g., Mustill et al.2015). Disk migration has been a persistently pro-posed explanation (e.g., Lee & Peale 2002) – particu-larly for systems in or near mean motion resonance –but recent studies have argued that in situ formationmay be more consistent with these planets’ observedproperties (e.g., Huang et al. 2016; Frelikh et al. 2019;Anderson et al. 2020) and could possibly be consis-tent with resonant configurations (e.g., Dong & Dawson2016; MacDonald & Dawson 2018; Choksi & Chiang2020). Precise characterization of the orbital proper-ties and resonant behavior of individual systems can testorigins scenarios, complementary to population studies.One warm Jupiter system potentially amenableto such detailed characterization is TOI-216, whichhosts a pair of warm, large exoplanets discovered bythe
TESS
Mission (Dawson et al. 2019; Kipping et al.2019). Their masses and orbits were characterized us-ing transit timing variations (TTVs). However, previ- ously the
TESS
TTV dataset was not sufficiently preciseto break the degeneracy between mass and eccentric-ity that arises when we can measure the near-resonantTTV signal but not the chopping TTV signal (e.g.,Lithwick & Naoz 2011; Deck & Agol 2015). Differentpriors on mass and eccentricity led to two qualitativelydifferent solutions for the system (Dawson et al. 2019):a Jupiter-mass planet accompanied by a Saturn-massplanet librating in orbital resonance, and a puffy sub-Saturn-mass planet and puffy Neptune-mass planet nearbut not in orbital resonance. The solutions also differedin the planets’ free eccentricity, the eccentricity not as-sociated with the proximity to resonance. Since anyorigins scenario under consideration needs to be ableto account for the planets’ masses, free eccentricities,and resonant behavior, detailed characterization of theseproperties with a more extended dataset is warranted.Since our previous study (Dawson et al. 2019), sevenmore
TESS transits of planet b, three more
TESS tran-sits of planet c, a ground-based transit of planet b, andfive more ground-based transits of planet c have beenobserved, including recent ground-based transit observa-tions that significantly extend the baseline for measuringtransit timing variations beyond the observations con-ducted during the first year of the
TESS primary mis-sion. Furthermore, we have been conducting a ground-based radial-velocity campaign using HARPS, FEROS,
OI-216
TESS light curves, ground-basedlight curves, and ground-based radial velocities to pre-cisely characterize TOI-216b and c. In Section 2, wedescribe our analysis of the
TESS data and the obser-vation and analysis of ground-based light curves. Weidentify a weak stellar activity periodicity in the
TESS data that also shows up in the radial velocities. In Sec-tion 3, we present radial velocity measurements fromHARPS, FEROS, and PFS and show that they imme-diately confirm the higher mass, lower eccentricity so-lution. We investigate additional weak periodicities inthe radial-velocities and argue that they are caused bystellar activity. We jointly fit the transit timing vari-ations and radial velocities in Section 4 and determinethat the planet pair is librating in resonance, that theinner planet has a significant free eccentricity, and thatplanets have a small but significant mutual inclination.We present our conclusions – including implications forthe system’s origins – in Section 5. LIGHT CURVE ANALYSISThis paper is based on data from TESS Sectors 1–13 (2018 July 25 – 2019 July 17) and Sectors 27–30 (2020 July 4 – 2020 October 21), during whichTOI-216 was observed with CCD 1 on Camera 4, andfrom ground-based observatories. We use the publiclyavailable 2-min cadence
TESS data, which is processedwith the Science Processing Operations Center pipeline(Jenkins et al. 2016). We download the publicly avail-able data from the Mikulski Archive for Space Tele-scopes (MAST). The pipeline, a descendant of the
Ke-pler mission pipeline based at the NASA Ames ResearchCenter (Jenkins et al. 2002, 2010), analyzes target pixelpostage stamps that are obtained for pre-selected targetstars.We use the resources of the
TESS
Follow-up Ob-serving Program (TFOP) Working Group (WG) SubGroup 1 (SG1) to collect seeing-limited time-series pho-tometric follow-up of TOI-216. All photometric timeseries are publicly available on the Exoplanet Follow-up Observing Program for TESS (ExoFOP-TESS) web-site . Light curves observed on or before February24, 2019 are described in Dawson et al. (2019). Ournew time-series follow-up observations are listed in Ta-ble 1. We used the TESS Transit Finder , which isa customized version of the
Tapir software package(Jensen 2013), to schedule our transit observations. Un- https://tess.mit.edu/followup/ https://exofop.ipac.caltech.edu/tess/ less otherwise noted, the photometric data were ex-tracted using the AstroImageJ ( AIJ ) software pack-age (Collins et al. 2017). The facilities used to col-lect the new TOI-216 observations published here areLas Cumbres Observatory Global Telescope (LCOGT)network (Brown et al. 2013), Hazelwood Observatory(Churchill, Vic, Australia), El Sauce Observatory (Co-quimbo Province, Chile), and the Antarctic Searchfor Transiting ExoPlanets (ASTEP) observatory (Con-cordia Station, Antarctica). All LCOGT 1 m tele-scopes are equipped with the Sinistro camera, witha 4k ×
4k pixel Fairchild back illuminated CCD and a26 ′ × ′ FOV. The LCOGT images were calibrated usingthe standard LCOGT BANZAI pipeline (McCully et al.2018). Hazelwood is a private observatory with an f/8Planewave Instruments CDK12 0.32 m telescope and anSBIG STT3200 2.2k × ′ × ′ fieldof view. El Sauce is a private observatory which hostsa number of telescopes; the observations reported herewere carried out with a Planewave Instruments CDK140.36 m telescope and an SBIG STT-1603-3 1536k × ′ × ′ field of view. ASTEP is a 0.4 mtelescope with an FLI Proline 16801E 4k ×
4k CCD, giv-ing a 63 ′ × ′ field of view. The ASTEP photometricdata were extracted using a custom IDL-based pipeline.We fit the transit light curves (Fig. 1) using the TAPsoftware package (Gazak et al. 2012), which implementsMarkov Chain Monte Carlo using the Mandel & Agol(2002) transit model and the Carter & Winn (2009)wavelet likelihood function, with the modifications de-scribed in Dawson et al. (2014). The results are summa-rized in Table 2. For TESS light curves, we use the pre-search data conditioned (PDC) flux, which is correctedfor systematic (e.g., instrumental) trends using cotrend-ing basis vectors (Smith et al. 2012; Stumpe et al.2014). For all light curves, we use Carter & Winn (2009)wavelet likelihood function (which, for the red noisecomponent, assumes noise with an amplitude that scalesas frequency − ) with free parameters for the amplitudeof the red σ r and white noise σ r and a linear trend fitsimultaneously to each transit light curve segment withother transit parameters. For the ground-based obser-vations, we fit a linear trend to airmass instead of time.We assign each instrument and filter ( TESS , Hazelwoodg’ and Rc, LCOGT i’ and I, El Sauce Rc, and ASTEPRc) its own set of limb darkening parameters because ofthe different wavebands. We use one set of noise param-eters for all the
TESS light curves and an additionalset for each ground-based light curve. We adopt uni-form priors on the planet-to-star radius ratio ( R p /R ⋆ ),the impact parameter b (which can be either negative orpositive; we report | b | ), the mid transit time, the limb Dawson et al.
Table 1.
Observation Log: TOI-216/TIC 55652896
TOI-216 Date Telescope † Filter ExpT Exp Dur. Transit Ap. Radius FWHM(UTC) (sec) (N) (min) Coverage (arcsec) (arcsec)b 2019-10-30 LCOGT-CTIO-1.0 I 40 220 241 Full 5.8 2.42020-05-23 ASTEP-0.4 ∼ Rc 120 376 924 Full 10.2 5.3c 2018-12-16 LCOGT-SAAO-1.0 i ′
90 331 450 Full 5.8 2.12019-01-20 Hazelwood-0.3 g ′
240 101 449 Egress+70% 5.5 3.22019-02-24 El Sauce-0.36 Rc 30 514 303 Egress+90% 5.9 3.42019-06-07 ASTEP-0.4 ∼ Rc 120 549 1440 Full 12.0 5.02019-11-27 Hazelwood-0.3 Rc 60 114 255 Egress+60% 5.5 3.12019-12-31 LCOGT-SAAO-1.0 Ic 40 316 345 Ingress+90% 4.7 2.02020-02-24 Hazelwood-0.3 Rc 120 161 380 Egress+60% 5.5 3.12020-03-09 LCOGT-CTIO-1.0 I 40 146 176 Egress+30% 8.2 2.62020-06-21 ASTEP-0.4 ∼ Rc 120 272 656 Full 11.0 5.0 † Telescopes:LCOGT-CTIO-1.0: Las Cumbres Observatory - Cerro Tololo Inter-American Observatory (1.0 m)LCOGT-SAAO-1.0: Las Cumbres Observatory - South African Astronomical Observatory (1.0 m)LCOGT-SSO-1.0: Las Cumbres Observatory - Siding Spring Observatory (1.0 m)Hazelwood-0.3: Stockdale Private Observatory - Victoria, Australia (0.32 m)El Sauce-0.36: El Sauce Observatory - Coquimbo Province, Chile (0.36 m)ASTEP-0.4: Antarctic Search for Transiting ExoPlanets - Concordia Station, Antarctica (0.4 m) darkening coefficients q and q (Kipping 2013a), andthe slope and intercept of each transit segment’s lineartrend. For the grazing transits of the inner planet, weimpose a uniform prior on R p /R ⋆ from 0 to 0.17, withthe upper limit corresponding to a planet radius of 0.13solar radii (see Dawson et al. 2019 for details and justifi-cation). Despite a well-constrained transit depth (Table2), the inner planet’s radius ratio is highly uncertaindue to degeneracy between the radius ratio and impactparameter (see Fig. 3 of Dawson et al. 2019). We alsoimpose a uniform prior on the log of the light curvestellar density ρ circ . We use ρ circ , the stellar density de-rived from the light curve assuming a circular orbit, tocompute the Mandel & Agol (2002) model normalizedseparation of centers z = d/R ⋆ , assuming M p << M ⋆ and a circular orbit: z = ( ρ circ /ρ ⊙ ) / ( P/P ⊕ ) / ( a ⊕ /R ⊙ ) (1)where P is the orbital period, the subscript ⊕ denotesthe Earth, and the subscript ⊙ denotes the Sun. Wewill later combine the posteriors ρ circ from the lightcurve and ρ ⋆ (from Dawson et al. 2019’s isochrone fit)to constrain the orbital eccentricity (Section 4.2). Weperform an additional set of fits where we allow for adilution factor for the TESS light curves and find theresults to be indistinguishable. We measure a dilutionfactor of 1 . +0 . − . .We also perform two additional fits to look for tran-sit duration variations, following Dawson (2020). In thefirst fit, we allow the impact parameter to vary with the prior recommended by Dawson (2020). We find ten-tative evidence for the variation in the transit impactparameter for the inner planet, with an impact param-eter change scale of 0.007 +0 . − . . In the second fit, weallow ρ circ (Eqn. 1) and b to vary for each transit with aprior that corresponds to a uniform prior in transit du-ration (Dawson 2020). Again the inner planet exhibitstentative evidence for variation in its transit durations.We fit a line to the transit durations as a function ofmid-transit time and determine a three-sigma confi-dence interval on the slope of -0.4 – 4 seconds/day forthe inner planet and -1 – 0.4 seconds/day for the outerplanet.For comparison and to ensure that the results arenot sensitive to the correlated noise treatment, wealso fit the light curves using the exoplanet pack-age (Foreman-Mackey et al. 2019), which uses Gaussianprocess regression. Each set of TESS light curves alongwith eight ground-based light curves is modeled with alight curve transit model built from starry (Luger et al.2019) plus a Matern-3/2 Gaussian process kernel witha white noise term. Seven sets of limb-darkening pa-rameters (Kipping 2013b) are used for observations con-ducted in seven different filters. We use a log-uniformprior on stellar densities ( ρ circ ), log-normal prior ontransit depths, uniform prior on the impact parameters,and uniform prior on mid-transit times. We infer pos-teriors for each parameter using this approach that areconsistent within 1 sigma to our nominal fit above. OI-216 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5O-C (days)0.50.60.70.80.91.0 no r m a li z ed de t r ended f l u x TOI-216b
LCOGTASTEP -0.5 0.0 0.5O-C (days)0.40.50.60.70.80.91.0 no r m a li z ed de t r ended f l u x TOI-216c
LCOGTHazelwoodEl SauceASTEPHazelwoodLCOGTHazelwoodLCOGTASTEP
Figure 1.
Detrended light curves, spaced with arbitrary vertical offsets, and with a model light curve overplotted. The lightcurves are phased based on a constant orbital period linear ephemeris to show the TTVs. TESS data are publicly available fromMAST and ground-based data from the ExoFOP-TESS website.
Dawson et al.
Table 2 . Planet Parameters for TOI-216b and TOI-216c De-rived from the Light Curves
Parameter Value a TOI-216bPlanet-to-star radius ratio, R p /R ⋆ +0 . − . Transit depth [ppm] 4750 +140 − Planet radius, R p [ R ⊕ ] 8 +3 − Light curve stellar density b , ρ circ [ ρ ⊙ ] 1.1 +0 . − . Impact parameter, | b | +0 . − . Mid-transit times (days c,d ) 1325.328 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . LCOGT 1786.698 +0 . − . ASTEP 1993.486 +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . +0 . − . Best fit linear ephemeris:Period (days) 17.16073Epoch (days) 1324.5911TOI-216cPlanet-to-star radius ratio, R p /R ⋆ +0 . − . Transit depth [ppm] 18310 +120 − Planet radius, R p [ R ⊕ ] 10.1 +0 . − . Light curves stellar density, ρ circ [ ρ ⊙ ] 1.73 +0 . − . Impact parameter, | b | +0 . − . Mid-transit times (days c,d ) 1331.2850 +0 . − . +0 . − . +0 . − . Table 2 continued
Table 2 (continued)
Parameter Value a +0 . − . +0 . − . LCOGT 1469.4781 +0 . − . Hazelwood 1504.036 +0 . − . El Sauce 1538.5938 +0 . − . +0 . − . +0 . − . +0 . − . ASTEP 1642.2595 +0 . − . +0 . − . Hazelwood 1814.939 +0 . − . LCOGT 1849.4526 +0 . − . Hazelwood 1883.955 +0 . − . LCOGT 1918.4506 +0 . − . ASTEP 2021.8887 +0 . − . +0 . − . +0 . − . +0 . − . Best fit linear ephemeris:Period (days) 34.525528Epoch (days) 1331.4110 a As a summary statistic we report the median and 68.3% confidenceinterval of the posterior distribution. b Eqn. 1 c BJD - 2457000.0 days d Mid-transit times not otherwise noted are from
TESS light curves.
OI-216 eleanor (Feinstein et al. 2019) with a 15 x 15 tar-get pixel file, background size of 100, and custom squareaperture of 3x3 pixels centered on TOI-216. We follow eleanor ’s recommendation for background subtraction:the 1D postcard background for sectors 1, 2, 4, 6, 7, 8,11 and 13; the 1D target pixel file background for sec-tors 3 and 5; and the 2D target pixel file backgroundfor sectors 9 and 12. Then we mask out the transits andcompute a discrete autocorrelation function (DCF; Eqn.2 of Edelson & Krolik 1988), plotted in Fig. 2.The peak at 6.5 day and valley at approximately halfthat value are consistent with a 6.5 day periodicity.This periodicity could represent the rotation period ora shorter harmonic. A rotation period of ∼ − M ⊙ star on themain sequence (e.g., McQuillan et al. 2014), so the pe-riodicity could plausibly be an integer fraction of therotation period (e.g., one half, one third). The stel-lar radius of 0.747 +0 . − . R ⊙ (Dawson et al. 2019) and v sin i = 0 . ± .
70 km/s from the FEROS spectra cor-respond to a rotation period of 45 +225 − days assumingan edge-on orientation. A rotation period of 26 dayswould be compatible with the v sin i measurement, ora shorter rotation period would indicate a spin-orbitmisalignment. A periodicity of 13 days or 26 days ischallenging to detect in the TESS light curve (e.g.,Canto Martins et al. 2020). 13 days is close to
TESS ’s13.7 day orbital period and thus susceptible to correc-tions for systematics. Both are a significant fractionof the sector (27 days) and subject to imprecision institching together different segments. Given the manybumps and wiggles in the DCF, we do not consider thislightcurve detection of a stellar rotation harmonic defini-tive. RADIAL VELOCITY ANALYSISTOI-216 was monitored with three different high res-olution echelle spectrographs over a time span of 16months with the goal of obtaining precision radial veloc-ities to further constrain the masses and orbital parame-ters of the giant planets present in the TOI-216 system.These observations were performed in the context of theWarm gIaNts with tEss (WINE) collaboration, whichfocuses on the systematic characterization of TESS tran-siting giant planets with orbital periods longer than ≈ DC F Figure 2.
Discrete autocorrelation function (Eqn. 2 ofEdelson & Krolik 1988) of eleanor long cadence TOI-216light curve.
Table 3.
Light Curve Parameters a for the TOI-216 system q q σ r [ppm] σ w [ppm]TESS 0 . +0 . − . . +0 . − . +600 − +14 −
20 s 15000 +2000 − +26 − LCOGT i’ 0 . +0 . − . . +0 . − . +1970 − +40 − Hazelwood g’ 0 . +0 . − . . +0 . − . +3000 − +190 − El Sauce 0 . +0 . − . . +0 . − . − +80 − Hazelwood Rc: 0 . +0 . − . . +0 . − . Nov. 2019 9000 +3000 − +200 − Feb. 2020 10 − +160 − LCOGT I: 0 . +0 . − . . +0 . − . Oct. 2019 1400 +1500 − +110 − Dec. 2019 1200 +1100 − +40 − Mar. 2020 1600 +1400 − +50 − ASTEP 0 . +0 . − . . +0 . − . June 2019 13000 +2000 − +100 − May 2020 3800 +1300 − +50 − June 2020 11900 +1500 − +70 − a As a summary statistic we report the mode and 68.3% confidence interval of theposterior distribution.
We obtained 27 spectra with the Fibre-fed, ExtendedRange, ´Echelle Spectrograph (FEROS, Kaufer et al.1999) between November of 2018 and March of 2019.FEROS is mounted on the MPG 2.2 m telescope at theESO La Silla Observatory, and has a resolving power ofR ≈
48 000. Observations were performed with the si-multaneous calibration mode for tracing radial velocityvariations produced by environmental changes in the in-strument enclosure. The adopted exposure time of 1200seconds yielded spectra with signal-to-noise ratios in therange from 40 to 75. FEROS spectra were processed
Dawson et al. from raw data with the ceres pipeline (Brahm et al.2017), which delivers precision radial velocities and linebisector span measurements via cross-correlation witha binary mask resembling the spectral properties of aG2-type star. The radial velocity uncertainties for theFEROS observations of TOI-216 range between 7 and15 m s − . We remove two outliers from the FEROSradial-velocity time series at 1503.75 and 1521.57 days.All subsequent results do not include these outliers. Wehave checked that no results except the inferred jitterfor the FEROS dataset are sensitive to whether or notthe outliers are included.We observed TOI-216 on 15 different epochs be-tween December of 2018 and October of 2019 with theHigh Accuracy Radial velocity Planet Searcher (HARPSMayor et al. 2003) mounted on the ESO 3.6 m telescopeat the ESO La Silla Observatory, in Chile. We adoptedan exposure time of 1800 seconds for these observationsusing the high radial velocity accuracy mode (HAM,R ≈
115 000), which produced spectra with signal-to-noise ratios of ≈
40 per resolution element. As in thecase of FEROS, HARPS data for TOI-216 was processedwith the ceres pipeline, delivering radial velocity mea-surements with typical errors of ≈ − .TOI-216 was also monitored with the Planet FinderSpectrograph (PFS; Crane et al. 2006, 2008, 2010)mounted on the 6.5 m Magellan II Clay Telescope at LasCampanas Observatory (LCO), in Chile. These spectrawere obtained on 18 different nights, between Decem-ber of 2018 and March of 2020, using the 0.3 ′′ × ′′ slit,which delivers a resolving power of R = 130 000. Due toits moderate faintness, TOI-216 was observed with the3 × ≈ − . An iodine cell wasused in these observations as a reference for the wave-length calibration. The PFS data were processed witha custom IDL pipeline (Butler et al. 1996). Three con-secutive 1200 second iodine-free exposures of TOI-216were obtained to construct a stellar spectral templatefor disentangling the iodine spectra from the stellar onefor computing the radial velocities.It is immediately evident that the radial velocitiesshow good agreement with Dawson et al. (2019)’s highermass solution, which was fit to the earlier, transit timeonly dataset (Fig. 3), and with Kipping et al. (2019)’ssolution based on transit times from the first six sec-tors. The bottom panel of Fig. 3 shows the RVs phasedto the outer planet’s orbital period. We compute ageneralised Lomb-Scargle periodogram (Cumming et al.1999; Zechmeister & K¨urster 2009) in Fig. 4 and 5. Thex-axis is frequency f in cycles per day. The y-axis is the square root of the power, where we define power asPower f = Σ i ( v i , f − v i, ) σ i − Σ i v i, σ i i σ i , (2)where σ i is the reported uncertainty (Table 5) and v i, is the radial-velocity (Table 5) with the error weightedmean Σ i v i /σ i σ i for each of the three datasets (PFS,HARPS, and FEROS) subtracted. The sinusoidal func-tion v i , f is v i , f = A cos[2 πf ( t i − t )]+ BA sin[2 πf ( t i − t )]+ C k (3)where A and B are computed followingZechmeister & K¨urster (2009), k is each of the threedatasets, t is the time of the first observed radial ve-locity, and C k = − Σ i ( A cos[2 πf ( t i − t )]+ BA sin[2 πf ( t i − t )]) σ i Σ i σ i . (4)We see a peak in the periodogram at planet c’s or-bital period (Fig. 4, left panel). The noise-free c model(overplotted, dashed; a Keplerian signal computed us-ing the parameters of planet c from Table 4 sampled atthe observed times in Table 5) shows that many of theother peaks seen in the three periodograms are aliasesof planet c’s orbital period ( i.e., they caused by the ob-servational time sampling of the planet’s signal). Planetb’s signal is below the noise level (Fig. 3, bottom rightpanel; Figure 4, right panels). In the PFS dataset, in-cluding planet b improves the chi-squared from 403 to328 (for 18 datapoints and five additional parameters);for the HARPS and FEROS datasets, there is no im-provement in chi-squared. Given that planet b is barelydetected, we cannot rule out other planets in the systemwith smaller orbital radial velocity amplitudes.Systems containing an outer planet in or near a2:1 mean motion resonance with a less massive innerplanet can be mistaken for a single eccentric planet,because the first eccentricity harmonic appears at halfthe planet’s orbital period (e.g., Anglada-Escud´e et al.2010; K¨urster et al. 2015). In the case of the systemTOI-216, the lack of detection of TOI-216b is not dueto this phenomenon because solution we subtract off forplanet c has near 0 eccentricity. However, without priorknowledge that the system contains a resonant pair, ifwe only had the radial-velocity datasets and the datasetswere less noisy (or had more data points), we might besensitive to planet b’s signal but mistakenly interpret itas planet c’s eccentricity. OI-216 R V ( m / s ) HARPSFEROSPFSmodel (D+20)model 2 (D+19)model 1 (D+19) R V ( m / s ) Phased to c R V ( m / s ) Phased to b
Figure 3.
Top: Radial velocity measurements of TOI-216 from HARPS (gray), FEROS (red; outliers in lighter red), andPFS (blue). Best-fit models from this work (Table 4) and from Dawson et al. (2019)’s earlier analysis of transit times only areoverplotted. Row 2, left: RVs phased to planet c’s orbital period with same models as row 1 for planet c only (i.e., the RVvariation only due to planet c) overplotted. Row 2, right: Residuals of Table 4 model for planet c only, with b component ofthe model overplotted, phased to planet b’s orbital period.
The residuals of the two-planet fit (and one-planet fit)show evidence of a signal that we attribute to stellar ac-tivity. We examined the periodograms of the residualsof each of three datasets; the bisectors for the HARPSand FEROS datasets; and H alpha for the FEROSdataset, computed with the ceres pipeline, followingBoisse et al. (2009). Periodograms of the PFS residu-als and FEROS H alpha are shown in Figure 5. Thereare no strong peaks in any of the residual or activity in-dicator periodograms. Some residuals datasets exhibit(weak) peaks near the 6.5 day periodicity (HARPS resid-uals and bisectors; FEROS residuals) identified in the light curve (Section 2), or a multiple of 6.5 days (PFS;FEROS bisectors and H alpha). In Fig. 6, we plot theresiduals of the two-planet fit with best-fit sinusoids for6.5, 13, and 26 day periodicities for each dataset. Com-paring the three datasets, best-fit sinusoids are out ofphase with each other. The similarity of the periodici-ties to those seen in the light curve, the appearance ofthe 6.5 day periodicity in the HARPS bisectors, the 6.5and 26 day periodicities in the FEROS H alpha, and the13 day and 6.5 day periodicities in the FEROS bisectors,and the difference in phase among datasets suggest the0
Dawson et al. P o w e r / c Combined dataModel (noise free) c P o w e r / b Combined dataModel (noise free) b
Figure 4.
Left: Periodograms of the combined RV dataset (black solid), with periodogram of noise-free planet c only model (i.e., a Keplerian signal computed using the parameters of planet c from Table 4 sampled at the observed times in Table 5; graydashed) overplotted. Right: Same for residuals to planet c only model, with noise-free planet b only model overplotted (i.e., aKeplerian signal computed using the parameters of planet b from Table 4 sampled at the observed times in Table 5). P o w e r / PFS P o w e r / FEROS H alpha
Figure 5.
Periodograms of residuals of the best two planetfit (Table 4) for PSF and H alpha for FEROS. The 6.5, 13,and 26 day periodicities (associated with the periodicity inthe light curve, Fig. 2) are overplotted as dotted lines. weak periodicities do not arise from additional planets.They are possibly caused by stellar variability. JOINT FIT AND ORBITAL ARCHITECTUREHere we jointly fit the transit light curves (Section2) and radial-velocity measurements (Section 3) to pre-cisely measure the orbital parameters and masses ofboth planets. In Section 4.1, we describe our analysisof the transit timing measurements. In Section 4.2, wejointly fit the transit and radial-velocity measurements.4.1.
Transit timing variation analysis
TOI-216b and c exhibit transit timing variations –plotted in Fig. 7 – due to the near-resonant effect (e.g.,Agol et al. 2005; Lithwick et al. 2012). We have not yetobserved a full super-period (which depends on the plan- ets’ proximity to the 2:1 resonance). The amplitude de-pends on the perturbing planet’s mass and the free ec-centricity of the transiting and perturbing planets. Wepreviously found a significant free eccentricity, with theexact value and partition between planets degeneratewith planet mass (Dawson et al. 2019).We fit the transit times using our N-body integratorTTV model (Dawson et al. 2014). Our model containsfive parameters for each planet: the mass M , orbitalperiod P , mean longitude λ , eccentricity e , and argu-ment of periapse ω . All orbital elements are osculatingorbital elements at the epoch 1325.3279 days. For eachplanet, we fix the impact parameter b to the value inTable 2 and set the longitude of ascending node on thesky to Ω sky = 0. We use the conventional coordinatesystem where the X − Y plane is the sky plane and the Z axis points toward the observer. We will explore otherpossibilities for b and Ω sky in Section 4.2.We begin by fitting the transit times only. To ex-plore the degeneracy between mass and eccentricity, weuse the Levenberg-Marquardt alogrithm implemented inIDL mpfit (Markwardt 2009) to minimize the χ on agrid of ( M c , e b ). We first use a dataset consisting ofall the ground-based data and the TESS Year 1 data.We contour the total χ for thirty-four transit timesand ten free parameters, i.e., twenty-four degrees of free-dom in Fig. 8. Even with the additional transits sinceDawson et al. (2019), the fits still suffer from a degener-acy between mass and eccentricity, though the relation-ship is much tighter. We integrate a random sample of100 solutions with χ <
50 and find that they all libratein the 2:1 resonance with the resonant angle involvingthe longitude of periapse of the inner planet. Next weadd TESS Year 3 transits from the TESS Extended Mis-sion, which became available while this manuscript wasunder review. This very extended baseline for the TTVs
OI-216 R V ( m / s ) HARPSFEROSPFS R V ( m / s ) R V ( m / s ) Figure 6.
Residuals to the best two planet fit (Table 4)phased to periodicities of 26 days (top), 13 days (middle),and 6.5 days (bottom). Best-fit sinusoids for these periodic-ities are overplotted as dashed lines. Comparing the threedatasets, the best-fit sinusoids are out of phase with eachother, indicating that these signals are unlikely to arise fromadditional planets. -2000-1000010002000 O - C ( m i nu t e s ) TESSLCOGTHazelEl SauceASTEPmodel planet b1400 1600 1800 2000time (days)-400-2000200 O - C ( m i nu t e s ) planet c Figure 7.
Observed mid-transit times (diamonds) withsubtracted linear ephemeris for TOI-216b (top) and TOI-216c (bottom), with the best fit model overplotted (dia-monds, dotted line). substantially reduces the degeneracy between mass andeccentricity. We also contour the total χ for these forty-three transit times and ten free parameters, i.e., thirty-three degrees of freedom in Fig. 8.A resonant or near-resonant planet’s total eccentricityis the vector sum of its free and forced eccentricity. Theforced eccentricity vector is dictated by the resonant dy-namics, and the free eccentricity vector oscillates aboutthe tip of the forced vector. Dissipation (e.g., eccen-tricity damping from the disk) tends to damp the freeeccentricity, and other perturbations (e.g., from a thirdplanet) can excite it. We estimate the eccentricity com-ponents from simulations as e free = ( e max − e min ) /
2) and e forced = ( e max + e min ) / Joint transit and radial velocity fit
We perform a joint N-body fit to the transit times andradial velocities, imposing two additional constraints.One constraint is the transit exclusion intervals re-ported in Dawson et al. (2019). The other is the lightcurve joint posterior of stellar density ρ circ vs. im-pact parameter b (Table 2; Eqn. 1) combined withthe ρ ⋆ posterior from Dawson et al. (2019). Following2 Dawson et al. c (M Jup )0.000.050.100.150.200.250.30 e b TTVs only, Year 1/3TTVs only, Year 1Same model, RVs only
Figure 8.
Contours of χ for the fit to transit times only.With ground-based transits and Year 1 of TESS data (blue-dashed), fits to the transit times only result in degeneracybetween the inner planet’s (osculating) eccentricity and theouter planet’s mass. The levels are χ = 42, 47, 52, 62,and 77 and the best-fit solutions occupy the innermost con-tour. The black contours show χ for these same solutionscompared to the RV measurements, with an RV offset foreach instrument as the only free parameter. The levels are χ =540, 560, 600, and 650 and the best-fit solutions occupythe innermost contour. The RV measurements break thedegeneracy between mass and eccentricity. The addition ofYear 3 TESS Extended Mission data (through sector 30; redsolid) reduces the degeneracy between mass and eccentricityand shows good agreement with the RV data. The levels are χ = 65, 70, 75, 85, and 100 and the best-fit solutions occupythe innermost contour. Dawson & Johnson (2012), we convert the ρ circ vs. b toa g vs. b posterior according to g = ( ρ circ /ρ ⋆ ) / . Weuse the g vs. b posterior to add a term to the likelihoodbased on e and ω using g = (1 + e sin ω ) / (1 − e ) / .We compute the sky-plane inclination for the dynam-ical model from b according to sin i sky a/R ⋆ = (1 + e sin ω ) / (1 − e ).Following Dawson et al. (2014), we derive posteriorsfor the parameters using Markov Chain Monte Carlowith the Metropolis-Hastings algorithm. Instead of in-cluding the orbital period and mean longitude at epoch1325.3279 days as parameters in the MCMC, we opti-mize them at each jump (i.e., each MCMC step) usingthe Levenberg-Marquardt algorithm (i.e., optimizing thedynamical model). We also optimize the radial velocity e f r ee b RV + TTVsc RV + TTVsb TTVs 1/3c TTVs 1/3b TTVs 1c TTVs 1 c (M Jup )0.0010.0100.100 e f o r c ed Figure 9.
Long-term (10 days) behavior of solutions with χ <
60 (TTV-only fits with ground-based transits and Year1 of TESS data, gray and and black, which have discrete val-ues because they are performed on a grid); χ <
80 (TTV-only fits with addition of Year 3 TESS Extended Missiondata, orange and light blue); and from the MCMC poste-rior for the joint RV-TTV fit (red and blue, Table 4). Row1: e free (approximated as ( e max − e min ) / e forced (approximated as ( e max + e min ) / offset for each instrument at each jump. We fit for aradial velocity jitter term for each of the three instru-ments; even though we identified possible periodicitiesin Fig. 5, none are strong enough to justify to explic-itly modeling (e.g., with false alarm probabilities of 0.2,0.002, and 0.02 for the highest peak in the periodogramfor PFS, FEROS, and HARPS respectively; Cumming2004). We visually inspect each parameter for conver-gence. We know there is a mutual inclination perpendic-ular to the sky plane because the inner planet exhibitsgrazing transits and the outer planet does not, To alsoallow for a mutual inclination parallel to the sky plane,we fit for ∆Ω sky , the difference in longitude of ascendingnode.Following Dawson et al. (2019), we perform two fitswith different priors to assess the sensitivity to these pri-ors. The first solution (Table 4) imposes uniform priorson eccentricities and log uniform priors on mass (i.e., pri-ors that are uniform in log space). The second imposesuniform priors on mass and log-uniform priors on eccen-tricity. All other fitted parameters (orbital period, meanlongitude, argument of periapse, radial velocity jitters,difference in longitude of ascending node) have uniformpriors. With the dataset in Dawson et al. (2019), theresults were very sensitive to the priors; with the new OI-216 ◦ (95% confidence interval).To ensure that our results are robust and that theparameter space has been thoroughly explored by thefitting algorithm (particularly the degeneracy betweeneccentricity and mass), we carry out a fit using a dif-ferent N-body code and fitting algorithm. We usethe Python Tool for Transit Variations ( PyTTV ; Korth(2020)) to fit the transit times (Table 2), RVs (Ta-ble 5), and the stellar parameters reported in Dawson(2020). The parameter estimation is carried out bya joint-N-body fit using
Rebound with the IAS15 in-tegrator (Rein & Liu 2012; Rein & Spiegel 2015) and
Reboundx (Tamayo et al. 2020), to model all the observ-ables without approximations. Within the simulation,carried out in barycentric coordinates, a common coor-dinate system was chosen where the x − y plane is theplane of sky. The x -axis points to the East, the y -axispoints to the North, and the z -axis points to the ob-server. Ω is measured from East to North, and ω ismeasured from the plane of sky. For the parameter es-timation with PyTTV , the gravitational forces betweenthe planets and the influence of general relativity (GR)were considered, in case the influence of GR is signifi-cant enough to be visible in the TTVs; for this system,general relativistic precession does not significantly af-fect the TTVs. The parameter estimation is initializedusing Rayleigh priors on eccentricities (Van Eylen et al.2019) and uniform priors on log values for the plane-tary masses. The estimation of physical quantities fromthe TTV signal is done in two steps. First, the posteriormode is found using the Differential Evolution algorithmimplemented in
PyTransit (Parviainen 2015). The opti-mization is carried out varying the planetary masses, or-bital periods, inclinations, eccentricities and argumentsof periastron. The eccentricity and argument of peri-astron are mapped from sampling-parameters √ e cos ω and √ e sin ω . The longitudes of the ascending nodes arefixed for both planets. After the posterior mode is found,a sample from the posterior is obtained using the affine-invariant Monte Carlo Markov Chain ensemble sampler emcee (Foreman-Mackey et al. 2013). The parametersand uncertainties are consistent with those reported inTable 4. 4.3. Dynamics and origin
We randomly draw 1000 solutions from the poste-rior for longer integrations of 10 days using mercury6 Table 4.
Planet Parameters (OsculatingOrbital Elements at Epoch 1325.3279 days)for TOI-216b and TOI-216c Derived fromjoint TTV/RV fit
Parameter Value a M ⋆ ( M ⊙ ) b R ⋆ ( R ⊙ ) b M b ( M Jup ) 0.059 +0 . − . P b +0 . − . e b +0 . − . ̟ b (deg.) 291.8 +0 . − . λ b (deg) 82.5 +0 . − . Ω b, sky (deg) 0 i b, sky (deg) 88.60 +0 . − . M c ( M Jup ) 0.56 +0 . − . P c +0 . − . e c +0 . − . ̟ c (deg.) 190 +30 − λ c (deg) 27.8 +1 . − . ∆Ω sky = Ω c, sky (deg) -1 +2 − i b, sky (deg) 89.84 +0 . − . λ c − λ b − ̟ b (deg). 42.5 +0 . − . i mut (deg) c +1 . − . Jitter (m/s):HARPS 8 +3 − FEROS 22 +5 − PFS 7.7 +1 . − . a As a summary statistic we report the medianand 68.3% confidence interval of the poste-rior distribution. An example fit with highprecision suitable for numerical integration isgiven in Table 6. b Stellar parameters fixed to the val-ues reported by Dawson et al. (2019):0.77 ± . M ⊙ and 0.748 ± . R ⊙ . Uncer-tainties in estimated planetary masses onlyaccount for the dynamical fitting, i.e., they do not include uncertainties in the star’smass. c
95% confidence interval is 1.2–3.9 ◦ . The99.7% confidence interval is 1.1–4.3 ◦ . (Chambers et al. 1996). We compute the libration am-plitudes for the 2:1 resonant angle 2 λ c − λ b − ̟ b , wherethe longitude of periapse ̟ b = ω b + Ω b . We can nowdefinitively determine that the system is librating inresonance. The libration amplitude is well-constrained4 Dawson et al. to 60 +2 − degrees. The very small uncertainty in thelibration amplitude leads us to believe that its mod-erate value is real, not an artifact of imperfect char-acterization (e.g., Millholland et al. 2018). Fig. 10shows a trajectory demonstrating libration about a fixedpoint. In our solutions, the other 2:1 resonant an-gle 2 λ c − λ b − ̟ c always circulates. However, theouter planet’s eccentricity and longitude of periapseare poorly constrained, so it is possible that solutionswhere this angle librates are also consistent with thedata. The inner planet has a forced eccentricity (Fig.9) of 0 . +0 . − . and a modest but significant free ec-centricity of 0 . +0 . − . . The outer planet has verysmall forced and free eccentricities of 0 . +0 . − . and0 . +0 . − . respectively. In our integrated solutions,we find that the timescale for the biggest variations inboth planets’ semi-major axes and the inner planet’seccentricity is the resonant libration timescale (approx-imately 5 years), and the outer planet’s eccentricityvaries on apsidal alignment timescale (approximately 25years). This system falls into the regime of resonantTTVs (e.g., Nesvorn´y & Vokrouhlick´y 2016) – ratherthan the more commonly characterized near-resonantTTVs (e.g., Lithwick et al. 2012)– and we expect theTTVs to oscillate on the libration timescale, which hasnot yet been fully covered by the observational baseline(Fig. 7). We integrate a random subset of 200 solutionsfor 1 Myr and find that all are stable.One possible scenario for establishing the resonance isconvergent migration of planet c preceded or followed byperturbations of planet b by another body (Fig. 11). Weshow example simulations that apply a migration forceusing the user-defined force feature of mercury6 (as de-scribed in Wolff et al. 2012); the simulations are done in3D with initial inclinations set to the present day values.Planet c migrates a short distance (0.8% of its initialsemi-major axis) toward planet b and captures b intoresonance. In the first example (top), planet b startswith a modestly eccentric (bottom) orbit ( e = 0 . e = 0 .
02) and is captured into resonancewith a tight libration amplitude. At about 80,000 years,planet b’s orbit is disturbed, which we approximate asan instantaneous change in the magnitude and direc-tion of its eccentricity vector. The first history resultsin large amplitude libration of the resonant angle involv-ing ̟ c and the second in circulation; because we are notconfident that only circulating resonant angles are com-patible with the data, we cannot use this distinction tofavor one history over the other.These example dynamical histories are compatiblewith in situ formation, but are potentially compatible -0.2 -0.1 0.0 0.1 0.2e b cos(2 λ c - λ b - ϖ b )-0.2-0.10.00.10.2 e b s i n ( λ c - λ b - ϖ b ) Figure 10.
Example trajectory for TOI-216 solution. Thetrajectory does not pass through the origin, indicating libra-tion about a fixed point. The offset from the origin is theforced eccentricity. with long distance migration as well. In the in situ for-mation scenario, three or more planets form in situ, andplanet c migrates a tiny distance toward b. A thirdplanet jostles b – exciting its eccentricity and mutualinclination – before or after c’s migration. In the longdistance migration scenario, the disturbing third planetcould have migrated earlier. As future work, orbital dy-namics simulations could place limits on the propertiesof this possible third planet. Although we have invokeda third planet in the example scenarios, we have notruled out the possibility that planet b itself disturbedplanet c, in a process separate from migration.Another hypothesis for the mutual inclinationis that it was excited when the planets passedthrough the 4:2 inclination resonance during migra-tion (Thommes & Lissauer 2003). Precession from theproto-planetary disk can separate the 4:2 inclination res-onances from the 2:1 eccentricity resonance. Figure 12shows a proof of concept simulation where the poten-tial of the proto-planetary disk is approximated as a J oblateness term for the stellar potential. A small butsignificant mutual inclination is excited.We place planet b and c on a mass-radius plot of warmexoplanets in Fig. 13. TOI-216c has a typical radius forits mass; its bulk density is 0.885 +0 . − . gcm − . TOI-216b likely also has a typical radius for its mass, butbecause of its grazing transit (Section 2) and the result- OI-216 -180-90090180 λ c - λ b - ϖ b ( deg ) λ c - λ b - ϖ c ( deg ) a ( A U ) e -180-90090180 λ c - λ b - ϖ b ( deg ) λ c - λ b - ϖ c ( deg ) a ( A U ) e Figure 11.
Examples of resonance capture through short distance convergent migration of planet c. Rows 1–2: resonantangle; row 3: semi-major axis of planet c (top black line) and planet b (bottom black line); row 4: eccentricity of planet b(top black line) and planet b (bottom black line). Values from simulations with the present-day observed orbits (i.e., range ofoscillation) are marked with a dotted red line. ing degeneracy with impact parameter, we cannot ruleout a large radius that would result in a low density forits mass compared to similar mass planets. Its poorly-constrained bulk density is 0.17 +0 . − . gcm − . CONCLUSIONTOI-216b and c are now a very precisely character-ized (with the exception of planet b’s radius) pair ofwarm, large exoplanets. Radial velocity measurementsusing HARPS, FEROS, and PFS broke a degeneracy be-tween mass and eccentricity in the TTV-only fits thatwas particularly severe before the TESS Extended Mis-sion observations, and an expanded TTV baseline from
TESS and an ongoing ground-based transit observingcampaign increased the precision of the fits. We cannow better assess the consistency of its properties withdifferent theories for formation and evolution of giantplanets orbiting close to their stars.We now know that TOI-216c is a warm Jupiter (0 . ± .
02 Jupiter masses) with a mass and radius typical of other 10–200 day giant planets (Fig. 13). TOI-216b hasa mass similar to Neptune’s (18.4 ± . M ⊕ ). Its radius isnot well constrained due to its grazing transits; its radiusmay very well be typical for its mass (Fig. 13). Giventhe large uncertainty in planet b’s radius, no inflationmechanisms or scenarios requiring formation beyond theice line (e.g., Lee & Chiang 2016) are required to explaineither planet’s radius.Furthermore, we know now that TOI-216b and care not just near but librating in the 2:1 resonance.The argument involving the longitude of periapse ofTOI-216b librates with an amplitude of 63 +3 − degrees.TOI-216b has a small but significant free eccentricity0 . +0 . − . . The mutual inclination with respect toTOI-216c is between 1.2–3.9 ◦ (95% confidence inter-val). The libration amplitude, free eccentricity, and mu-tual inclination imply a disturbance of TOI-216b beforeor after resonance capture, perhaps by an undetectedthird planet. The orbital properties can be consistenteither with in situ formation, with resonance capture6 Dawson et al. -180-90090180 λ c - λ b - ϖ b ( deg ) i ( deg ) a ( A U ) e Figure 12.
Examples of excitation of a mutual inclinationduring short distance convergent migration of planet c. Row1: eccentricity resonant angle; row 2: mutual inclination;row 3: semi-major axis of planet c (top black line) and planetb (bottom black line); row 4: eccentricity of planet b (topblack line) and planet b (bottom black line). Values fromsimulations with the present-day observed orbits (i.e., rangeof oscillation) are marked with a dotted red line. through a very short distance migration, or long dis-tance migration. Future origins scenarios must matchthese precisely constrained properties. Future atmo-spheric characterization by the James Webb Space Tele-scope (JWST) may help distinguish between these ori-gin scenarios. If the planet formed outside the watersnow line, we expect its C/O ratio to be significantlysmaller than that of its host star; in-situ formation, onthe other hand, would imply C/O ratios closer to thestar (Espinoza et al. 2017). Simulations performed withPandExo (Batalha et al. 2017) show that water featuresin the spectra that constrain the C/O ratio may be de-tectable for TOI-216c with the JWST Near Infrared Im-ager and Slitless Spectrograph (NIRISS), even beneatha moderate cloud layer. p (M Earth ) TOI-216b TOI-216c R p ( R E a r t h ) Figure 13.
Warm (10–200 day orbital period) planetswith both mass and radius measurements (exoplanet.eu;Schneider et al. 2011), including TOI-216 (orange).
This system will benefit from continued long-termradial-velocity and transit monitoring. Long-term ra-dial velocity monitoring could reveal the presence of ad-ditional planets in the system, though stellar activityposes a challenge for detecting low amplitude signals.Observations by the
TESS
Extended Mission are con-tinuing and will further increase the baseline of obser-vations. The new
TESS observations may allow us tobetter constrain the change in impact parameter tenta-tively detected here for TOI-216b, allowing for tighterconstraints on the mutual inclination.ACKNOWLEDGMENTSWe thank the
TESS
Mission team and follow up work-ing group for the valuable dataset. We acknowledge theuse of public
TESS
Alert data from pipelines at the
TESS
Science Office and at the
TESS
Science Process-ing Operations Center. This paper includes data col-lected by the
TESS mission, which are publicly avail-able from the Mikulski Archive for Space Telescopes(MAST). Resources supporting this work were providedby the NASA High-End Computing (HEC) Programthrough the NASA Advanced Supercomputing (NAS)Division at Ames Research Center for the production ofthe SPOC data products.This research has made use of the ExoplanetFollow-up Observation Program website, which isoperated by the California Institute of Technology,under contract with the National Aeronautics and
OI-216
Gaia
Gaia
Gaia
Multilateral Agreement. We acknowledgesupport for ASTEP from the French and Italian PolarAgencies, IPEV and PNRA, and from Universit´e Cˆoted’Azur under Idex UCAJEDI (ANR-15-IDEX-01). Thispaper includes data gathered with the 6.5 meter Mag-ellan Telescopes located at Las Campanas Observatory,Chile.RID and JD gratefully acknowledge support by NASAXRP NNX16AB50G, NASA XRP 80NSSC18K0355,NASA
TESS
GO 80NSSC18K1695, and the Alfred P.Sloan Foundation’s Sloan Research Fellowship. TheCenter for Exoplanets and Habitable Worlds is sup-ported by the Pennsylvania State University, the EberlyCollege of Science, and the Pennsylvania Space GrantConsortium. JD gratefully acknowledges support andhospitality from the pre-doctoral program at the Cen-ter for Computational Astrophysics, Flatiron Institute.Research at the Flatiron Institute is supported by the Si-mons Foundation. This research made use of computingfacilities from Penn State’s Institute for CyberScienceAdvanced CyberInfrastructure.We thank Caleb Ca˜nas and David Nesvorn´y for helpfulcomments and discussions. We thank the referee for ahelpful report.R.B. acknowledges support from FONDECYTProject 11200751 and from CORFO project N ◦ ◦ ◦ ST/S00193X/1). J. Korth acknowledges support byDFG grants PA525/19-1 within the DFG SchwerpunktSPP 1992, Exploring the Diversity of Extrasolar Plan-ets. Part of this research was carried out at the JetPropulsion Laboratory, California Institute of Technol- ogy, under a contract with the National Aeronauticsand Space Administration (NASA).
Software: mpfit (Markwardt 2009),TAP(Gazak et al. 2012), Tapir (Jensen 2013),
TESS pipeline (Jenkins et al. 2016; Twicken et al. 2018;Li et al. 2019), AstroImageJ (Collins et al. 2017), exoplanet (Foreman-Mackey et al. 2019), astropy (Astropy Collaboration et al. 2013, 2018), celerite (Foreman-Mackey et al. 2017; Foreman-Mackey 2018), starry (Luger et al. 2019), pymc3 (Salvatier et al. 2016), theano (Theano Development Team 2016), ceres (Brahm et al. 2017),
PyTTV , PyTransit (Parviainen2015), emcee (Foreman-Mackey et al. 2013), eleanor (Feinstein et al. 2019)8
Dawson et al.
REFERENCES
Agol, E., Steffen, J., Sari, R., & Clarkson, W. 2005,MNRAS, 359, 567, doi: 10.1111/j.1365-2966.2005.08922.xAnderson, K. R., Lai, D., & Pu, B. 2020, MNRAS, 491,1369, doi: 10.1093/mnras/stz3119Anglada-Escud´e, G., L´opez-Morales, M., & Chambers, J. E.2010, ApJ, 709, 168, doi: 10.1088/0004-637X/709/1/168Astropy Collaboration, Robitaille, T. P., Tollerud, E. J.,et al. 2013, A&A, 558, A33,doi: 10.1051/0004-6361/201322068Astropy Collaboration, Price-Whelan, A. M., Sip˝ocz, B. M.,et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4fBatalha, N. E., Mandell, A., Pontoppidan, K., et al. 2017,PASP, 129, 064501, doi: 10.1088/1538-3873/aa65b0Boisse, I., Moutou, C., Vidal-Madjar, A., et al. 2009, A&A,495, 959, doi: 10.1051/0004-6361:200810648Brahm, R., Jord´an, A., & Espinoza, N. 2017, PASP, 129,034002, doi: 10.1088/1538-3873/aa5455Brahm, R., Espinoza, N., Jord´an, A., et al. 2019, AJ, 158,45, doi: 10.3847/1538-3881/ab279aBrown, T. M., Baliber, N., Bianco, F. B., et al. 2013,Publications of the Astronomical Society of the Pacific,125, 1031, doi: 10.1086/673168Butler, R. P., Marcy, G. W., Williams, E., et al. 1996,PASP, 108, 500, doi: 10.1086/133755Canto Martins, B. L., Gomes, R. L., Messias, Y. S., et al.2020, arXiv e-prints, arXiv:2007.03079.https://arxiv.org/abs/2007.03079Carter, J. A., & Winn, J. N. 2009, ApJ, 704, 51,doi: 10.1088/0004-637X/704/1/51Chambers, J. E., Wetherill, G. W., & Boss, A. P. 1996,Icarus, 119, 261, doi: 10.1006/icar.1996.0019Choksi, N., & Chiang, E. 2020, MNRAS,doi: 10.1093/mnras/staa1421Collins, K. A., Kielkopf, J. F., Stassun, K. G., & Hessman,F. V. 2017, AJ, 153, 77, doi: 10.3847/1538-3881/153/2/77Crane, J. D., Shectman, S. A., & Butler, R. P. 2006, Societyof Photo-Optical Instrumentation Engineers (SPIE)Conference Series, Vol. 6269, The Carnegie Planet FinderSpectrograph, 626931, doi: 10.1117/12.672339Crane, J. D., Shectman, S. A., Butler, R. P., et al. 2010,Society of Photo-Optical Instrumentation Engineers(SPIE) Conference Series, Vol. 7735, The CarnegiePlanet Finder Spectrograph: integration andcommissioning, 773553, doi: 10.1117/12.857792Crane, J. D., Shectman, S. A., Butler, R. P., Thompson,I. B., & Burley, G. S. 2008, Society of Photo-OpticalInstrumentation Engineers (SPIE) Conference Series,Vol. 7014, The Carnegie Planet Finder Spectrograph: astatus report, 701479, doi: 10.1117/12.789637 Cumming, A. 2004, MNRAS, 354, 1165,doi: 10.1111/j.1365-2966.2004.08275.xCumming, A., Marcy, G. W., & Butler, R. P. 1999, ApJ,526, 890, doi: 10.1086/308020Dawson, R. I. 2020, AJ, 159, 223,doi: 10.3847/1538-3881/ab7fa5Dawson, R. I., & Johnson, J. A. 2012, ApJ, 756, 122,doi: 10.1088/0004-637X/756/2/122—. 2018, ARA&A, 56, 175,doi: 10.1146/annurev-astro-081817-051853Dawson, R. I., Johnson, J. A., Fabrycky, D. C., et al. 2014,ApJ, 791, 89, doi: 10.1088/0004-637X/791/2/89Dawson, R. I., Huang, C. X., Lissauer, J. J., et al. 2019,AJ, 158, 65, doi: 10.3847/1538-3881/ab24baDeck, K. M., & Agol, E. 2015, ApJ, 802, 116,doi: 10.1088/0004-637X/802/2/116Dong, R., & Dawson, R. 2016, ApJ, 825, 77,doi: 10.3847/0004-637X/825/1/77Edelson, R. A., & Krolik, J. H. 1988, ApJ, 333, 646,doi: 10.1086/166773Espinoza, N., Fortney, J. J., Miguel, Y., Thorngren, D., &Murray-Clay, R. 2017, ApJ, 838, L9,doi: 10.3847/2041-8213/aa65caFeinstein, A. D., Montet, B. T., Foreman-Mackey, D., et al.2019, PASP, 131, 094502, doi: 10.1088/1538-3873/ab291cForeman-Mackey, D. 2018, Research Notes of the AmericanAstronomical Society, 2, 31,doi: 10.3847/2515-5172/aaaf6cForeman-Mackey, D., Agol, E., Ambikasaran, S., & Angus,R. 2017, AJ, 154, 220, doi: 10.3847/1538-3881/aa9332Foreman-Mackey, D., Czekala, I., Luger, R., et al. 2019,dfm/exoplanet: exoplanet v0.2.1,doi: 10.5281/zenodo.3462740Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman,J. 2013, PASP, 125, 306, doi: 10.1086/670067Frelikh, R., Jang, H., Murray-Clay, R. A., & Petrovich, C.2019, ApJ, 884, L47, doi: 10.3847/2041-8213/ab4a7bGazak, J. Z., Johnson, J. A., Tonry, J., et al. 2012,Advances in Astronomy, 2012, 697967,doi: 10.1155/2012/697967Huang, C., Wu, Y., & Triaud, A. H. M. J. 2016, ApJ, 825,98, doi: 10.3847/0004-637X/825/2/98Jenkins, J. M., Caldwell, D. A., & Borucki, W. J. 2002,ApJ, 564, 495, doi: 10.1086/324143Jenkins, J. M., Caldwell, D. A., Chandrasekaran, H., et al.2010, ApJ, 713, L87, doi: 10.1088/2041-8205/713/2/L87Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, inProc. SPIE, Vol. 9913, Software and Cyberinfrastructurefor Astronomy IV, 99133E, doi: 10.1117/12.2233418
OI-216 Dawson et al.
Table 5 . HARPS,FEROS, and PFS radial velocity measurements ofTOI-216.BJD RV σ RV BIS σ BIS
Inst. Note-2450000 (m s − ) (m s − ) (m s − ) (m s − )8449.70722 36690.0 7.8 -16.0 12.0 FEROS8450.62181 36709.1 8.0 -12.0 12.0 FEROS8450.75019 36686.8 8.3 -35.0 12.0 FEROS8451.72708 36706.3 7.7 -21.0 12.0 FEROS8451.74212 36696.3 7.2 -44.0 11.0 FEROS8451.75750 36694.2 7.2 -25.0 11.0 FEROS8452.77587 36713.0 7.6 -48.0 11.0 FEROS8464.76381 36754.4 8.7 11.0 11.0 HARPS8464.81000 36744.7 10.3 19.0 13.0 HARPS8465.73524 36759.3 4.1 9.0 5.0 HARPS8465.75656 36753.4 4.6 5.0 6.0 HARPS8466.72641 36733.7 6.4 -5.0 8.0 HARPS8466.74802 36730.1 7.4 -11.0 10.0 HARPS8467.74159 36714.5 7.9 -35.0 12.0 FEROS8468.67914 10.49 1.11 . . . . . . . . . . . . . . PFS8468.74406 36727.7 7.7 0.0 12.0 FEROS8472.76399 1.38 1.23 . . . . . . . . . . . . . . PFS8476.74101 -23.94 1.18 . . . . . . . . . . . . . . PFS8479.63075 -33.40 1.15 . . . . . . . . . . . . . . PFS8480.71817 -33.50 1.33 . . . . . . . . . . . . . . PFS8481.63918 36681.1 5.6 3.0 7.0 HARPS8481.66088 36670.1 4.9 -7.0 6.0 HARPS8482.63750 36690.5 4.1 2.0 5.0 HARPS8482.65892 36684.0 4.3 -1.0 6.0 HARPS8483.63626 36702.2 4.3 -8.0 6.0 HARPS8483.65821 36704.1 4.1 11.0 5.0 HARPS8483.76291 36668.1 8.0 -1.0 12.0 FEROS8485.70826 36650.9 13.9 47.0 19.0 FEROS8493.72558 36735.9 9.0 -11.0 13.0 FEROS8497.62083 36696.5 7.8 -40.0 12.0 FEROS8501.68699 34.11 1.72 . . . . . . . . . . . . . . PFS8503.65356 36813.1 14.3 64.0 19.0 FEROS Outlier8504.65888 36684.4 9.0 -61.0 13.0 FEROS8507.66569 -11.68 2.42 . . . . . . . . . . . . . . PFS8510.68697 -27.18 1.88 . . . . . . . . . . . . . . PFS8521.56775 36625.2 7.6 5.0 11.0 FEROS Outlier8528.60948 46.62 1.76 . . . . . . . . . . . . . . PFS8529.58589 45.56 1.65 . . . . . . . . . . . . . . PFS8542.55798 36659.1 8.7 -44.0 13.0 FEROS8543.54347 36640.3 7.4 -6.0 11.0 FEROS8544.61674 36663.3 7.8 -26.0 12.0 FEROS8545.57452 36676.3 7.6 -3.0 11.0 FEROS8546.58251 36714.8 9.7 38.0 14.0 FEROS8547.53676 36635.1 7.9 -17.0 12.0 FEROS8548.58309 36652.2 7.7 23.0 12.0 FEROS OI-216 Table 5 . Continued.BJD RV σ RV BIS σ BIS
Inst. Note-2450000 (m s − ) (m s − ) (m s − ) (m s − )8550.61335 36613.1 8.2 -62.0 12.0 FEROS8551.62488 36637.6 11.4 -26.0 16.0 FEROS8554.56409 36707.0 8.9 -38.0 13.0 FEROS8557.56708 36713.0 9.2 -4.0 13.0 FEROS8708.90979 28.50 2.14 . . . . . . . . . . . . . . PFS8711.89325 12.38 2.62 . . . . . . . . . . . . . . PFS8714.92414 -16.98 2.14 . . . . . . . . . . . . . . PFS8742.88505 19.06 2.19 . . . . . . . . . . . . . . PFS8764.84159 36736.7 5.2 -1.0 7.0 HARPS8766.79176 36750.6 6.4 3.0 8.0 HARPS8767.86802 29.71 1.98 . . . . . . . . . . . . . . PFS8777.79725 36695.8 23.5 23.0 31.0 HARPS8914.55281 21.41 2.73 . . . . . . . . . . . . . . PFS8920.52884 -5.82 1.77 . . . . . . . . . . . . . . PFS8923.53817 -33.84 1.94 . . . . . . . . . . . . . . PFS Dawson et al.
Table 6.
Example Fit with HighPrecision Parameters (Osculat-ing Orbital Elements at Epoch1325.3279 days) for TOI-216band TOI-216c Derived from jointTTV/RV fit
Parameter Value M ⋆ ( M ⊙ ) 0.77 M b ( M Jup ) 0.061381 P b e b ̟ b (deg.) 292.55926 λ b (deg) 82.329609Ω b, sky (deg) 0 i b, sky (deg) 88.594066 M c ( M Jup ) 0.55454 P c e c ̟ c (deg.) 200.93260 λ c (deg) 27.931018∆Ω sky = Ω c, sky (deg) -0.20544800 i b, skysky