A seven-planet resonant chain in TRAPPIST-1

TL;DR

The paper reports the discovery of a seven-planet resonant chain in the TRAPPIST-1 system, with new orbital data, resonance relations, and planetary characteristics, enhancing understanding of its complex orbital dynamics.

Contribution

It provides the first measurement of TRAPPIST-1h's orbital period and confirms a resonant chain with three-body resonances linking all seven planets.

Findings

TRAPPIST-1h has an orbital period of 18.77 days.

The planetary system forms a complex resonant chain.

TRAPPIST-1's star has a rotational period of 3.3 days.

Abstract

The TRAPPIST-1 system is the first transiting planet system found orbiting an ultra-cool dwarf star. At least seven planets similar to Earth in radius and in mass were previously found to transit this host star. Subsequently, TRAPPIST-1 was observed as part of the K2 mission and, with these new data, we report the measurement of an 18.77 d orbital period for the outermost planet, TRAPPIST-1h, which was unconstrained until now. This value matches our theoretical expectations based on Laplace relations and places TRAPPIST-1h as the seventh member of a complex chain, with three-body resonances linking every member. We find that TRAPPIST-1h has a radius of 0.727 Earth radii and an equilibrium temperature of 173 K. We have also measured the rotational period of the star at 3.3 d and detected a number of flares consistent with a low-activity, middle-aged, late M dwarf.

Figures

Tables

Table 1. Table 1: The three-body resonances of TRAPPIST-1. The transit times are used to track the ϕ italic-ϕ \phi angles of each set of three adjacent planets over the dataset, assuming low eccentricities such that transits occur at a phase angle λ = 90 ∘ 𝜆 superscript 90 \lambda=90^{\circ} [ 11 ] . The ranges of three-body frequency and angle given encompass the changes — most likely librations — seen during the observations.

planets 1,2,3	$p$	$q$	$\frac{p}{P_1}-\frac{\left(p+q\right)}{P_2}+\frac{q}{P_3}$ [d-1]	$\varphi =p{\lambda }_1-\left(p+q\right){\lambda }_2+q{\lambda }_3$
b,c,d	2	3	$[-4.6,-0.3]\times {10}^{-5}$	[${176}^{\circ }$, ${178}^{\circ }$]
c,d,e	1	2	$[-5.2,4.5]\times {10}^{-5}$	[${47}^{\circ }$, ${50}^{\circ }$]
d,e,f	2	3	$[-1.9,+1.9]\times {10}^{-4}$	[$-{154}^{\circ }$, $-{142}^{\circ }$]
e,f,g	1	2	$[-1.4,+1.1]\times {10}^{-4}$	[$-{79}^{\circ }$, $-{72}^{\circ }$]
f,g,h	1	1	$[-6.0,+0.2]\times {10}^{-5}$	[${176.5}^{\circ }$, ${177.5}^{\circ }$]

Table 2. Table 2: Properties of TRAPPIST-1h, limb-darkening parameters and transit timings derived using a joint Spitzer and K2 dataset. Parameter values are the medians of the posterior distributions from the MCMC and the associated error bars are the 1-sigma credible intervals.

Parameter	Value
Transit depth ${\left(R_p/R_{\ast }\right)}^2$ [%]	0.346$\pm 0.018$
Transit duration [d]	0.0525$\pm$0.0008
Impact parameter $b$ [$R_{\ast }$]	0.45${}_{}{}^{+0.06}{}_{-0.08}{}^{}$
Mid-transit time, $T_0$ (${\mathrm{BJD}}_{\mathrm{TDB}}$)	2,457,662.55284$\pm 0.00037$
Period, P [days]	18.767${}_{}{}^{+0.004}{}_{-0.003}{}^{}$
Radius ratio, $R_p/R_{\ast }$	0.0588$\pm 0.0016$
Radius, $R_p$ [$R_{\oplus }$]	0.752${}_{}{}^{+0.032}{}_{-0.031}{}^{}$
Inclination, $i$(∘)	89.76${}_{}{}^{+0.05}{}_{-0.04}{}^{}$
Scale parameter, $a/R_{\ast }$	109$\pm$4
Equilibrium temperature [K]	173$\pm$4
Irradiation, $S_p$ [$S_{\mathrm{Earth}}$]	0.165$\pm$0.025
Limb-darkening parameters (Kepler bandpass)
$u_1$	1.00$\pm$0.02
$u_2$	-0.04$\pm$0.04
Individual transit timings from K2 (${\mathrm{BJD}}_{\mathrm{TDB}}$)
Transit 1	2,457,756.3874${}_{}{}^{+0.0013}{}_{-0.0013}{}^{}$
Transit 2	2,457,775.1539${}_{}{}^{+0.0016}{}_{-0.0016}{}^{}$
Transit 3	2,457,793.9230${}_{}{}^{+0.0024}{}_{-0.0025}{}^{}$
Transit 4	2,457,812.6987${}_{}{}^{+0.0045}{}_{-0.0042}{}^{}$

Equations

$$h\sim 250{\left(\frac{{\varphi }_{int}}{1{\text{Wm}}^{-2}}\right)}^{-1}\text{m}$$

$$k_{xy}(x_i,y_i,x_j,y_j)=A_{xy}\exp \left[-\frac{(x_i-x_j{)}^2}{L_x^2}-\frac{(y_i-y_j{)}^2}{L_y^2}\right]$$

$$k_{xy}(t_i,t_j)=A_t\exp \left[-\frac{(t_i-t_j{)}^2}{L_t^2}\right]$$

$$K_{ij}=k_{xy}(x_i,y_i,x_j,y_j)+k_t(t_i,t_j)+{\sigma }^2{\delta }_{ij}$$

$$\log \mathcal{L}=-\frac{1}{2}(\vec{y}-{\vec{m}}_{t_0}{)}^{\top }{\mathbf{\Sigma }}^{-1}(\vec{y}-{\vec{m}}_{t_0})+C$$

$$\mathbf{\Sigma }=\mathbf{X\Lambda }{\mathbf{X}}^{\top }+\mathbf{K}$$

$$\Delta {\chi }_{\mathrm{cond}}^2=\Delta {\chi }^2-{\left(\frac{d-d_0}{{\sigma }_d}\right)}^2$$

$$d={\sigma }_d^2{\vec{m}}_{t_0}^{\top }{\mathbf{\Sigma }}^{-1}\vec{y}$$

$${\sigma }_d^2=({\vec{m}}_{t_0}^{\top }{\mathbf{\Sigma }}^{-1}{\vec{m}}_{t_0}{)}^{-1}$$

$$\lambda =\frac{360^{\circ }}{P}t,$$

$$\lambda =360^{\circ }\left(\frac{1}{4}+\frac{t-T_n}{P}\right)$$

$$\varphi =360^{\circ }\left[-p\frac{T_1}{P_1}+(p+q)\frac{T_2}{P_2}-q\frac{T_3}{P_3}+p\frac{t}{P_1}-(p+q)\frac{t}{P_2}+q\frac{t}{P_3}\right].$$

Full text

{affiliations}

Astronomy Department, University of Washington, Seattle, WA, 98195, USA

NASA Astrobiology Institute’s Virtual Planetary Laboratory, Seattle, WA, 98195, USA

University of Bern, Center for Space and Habitability, Sidlerstrasse 5, CH-3012, Bern, Switzerland

Laboratoire AIM Paris-Saclay, CEA/DRF - CNRS - Université Paris Diderot, IRFU/SAp Centre de Saclay, 91191, Gif-sur-Yvette

Department of Astronomy and Astrophysics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA

Space Sciences, Technologies and Astrophysics Research (STAR) Institute, Université de Liège, Allée du 6 Août 19C, B-4000 Liège, Belgium

Center for Astrophysics and Space Science, University of California San Diego, La Jolla, CA, 92093, USA

IPAC, Mail Code 314-6, California Institute of Technology, 1200 E California Boulevard, Pasadena, CA 91125, USA

Laboratoire d’Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France

Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK

NASA Ames Research Center, Moffett Field, CA 94043, USA

Bay Area Environmental Research Institute, 25 2nd St. Ste 209 Petaluma, CA 94952, USA

NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

Cavendish Laboratory, J J Thomson Avenue, Cambridge, CB3 0HE, UK

Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA

Jeremiah Horrocks Institute, University of Central Lancashire, Preston PR1 2HE, UK

NASA Johnson Space Center, 2101 NASA Parkway, Houston, Texas, 77058, USA

Laboratoire de Météorologie Dynamique, Sorbonne Universités, UPMC Univ Paris 06, CNRS, 4 place Jussieu, 75005 Paris, France

Space and Astronomy Department, Faculty of Science, King Abdulaziz University, 21589 Jeddah, Saudi Arabia

King Abdullah Centre for Crescent Observations and Astronomy (KACCOA), Makkah Clock, Saudia Arabia

LPHEA laboratory, Oukaïmeden Observatory, Cadi Ayyad University FSSM, BP 2390, Marrakesh, Morocco

Observatoire de Genève, Université de Genève, 51 chemin des Maillettes, CH-1290 Sauverny, Switzerland

Sagan Fellow

The star TRAPPIST-1 (EPIC 246199087) was observed for 79 days by NASA’s Kepler Space Telescope in its two-reaction wheel mission[4] (K2) as part of Campaign 12, starting on 2016 Dec 15 and ending on 2017 Mar 04. The spacecraft was in safe mode between 2017 Feb 1 and 2017 Feb 6, resulting in a 5-day data loss. Typically upon downlink from the spacecraft, the raw cadence data are calibrated with the Kepler pipeline[5], a lengthy procedure that includes background subtraction, smear removal, and undershoot and nonlinearity corrections. However, given the unique science drivers in this dataset, the raw, uncalibrated data for Campaign 12 were made publicly available on 2017 Mar 8 shortly after downlink. We download and calibrate the long cadence ($t_{\mathrm{exp}}=30\,\mathrm{min}$) and short cadence ($t_{\mathrm{exp}}=1\,\mathrm{min}$) light curves using a simple column-by-column background subtraction, which also removes smear and dark noise (see Methods). Because of its two failed reaction wheels, rolling motion of the Kepler spacecraft due to torque imbalances introduces strong instrumental signals, leading to an increase in photometric noise by a factor of ${\sim}3-5$ compared to the original mission. Since TRAPPIST-1 is a faint M8 dwarf with Kepler magnitude $Kp\sim 16-17$ (see Methods), these instrumental signals must be carefully removed to reach the $\sim$0.1% relative photometric precision required to detect Earth-size transits[6]. To this end, we detrend the long cadence light curve for TRAPPIST-1 using both EVEREST[7, 8] and a Gaussian process-based pipeline, achieving an average 6-hr photometric precision of 281.3 ppm, a factor of 3 improvement over the raw light curve. After analysis of the long cadence light curve, we detrend the short cadence light curve in the vicinity of features of interest, achieving comparable or higher 6-hr precision (see Methods).

We conduct three separate transit searches on the long cadence light curve, aiming to constrain the period of TRAPPIST-1h, which had only been observed to transit once[2], as well as to detect additional planets in the system. A dynamical analysis made by our team prior to the release of the K2 data suggested certain values of the period of TRAPPIST-1h based on the presence of three-body resonances among the planets. Three-body resonances satisfy $pP_{1}^{-1}-(p+q)P_{2}^{-1}+qP_{3}^{-1}\approx 0$ and $p\lambda_{1}-(p+q)\lambda_{2}+q\lambda_{3}=\phi$ for integers $p,q$ and where $P_{i}$ and $\lambda_{i}$ are the period and mean longitude of the $i^{\rm th}$ planet[9, 10] and $\phi$ is the 3-body angle which librates about a fixed value. Such resonances occur both in our Solar System — the archetypical case being the Laplace resonance among Jupiter’s satellites, satisfying $(p,q)=(1,2)$ — and in exoplanet systems, two of which were recently observed to have resonant chains among four planets: Kepler-223[11] with $(p,q)=(1,1)$ and Kepler-80[12] with $(p,q)=(2,3)$. Among the inner six planets in TRAPPIST-1, there are four adjacent sets of three planets that satisfy this relation for $1\leq p\leq 2$ and $1\leq q\leq 3$ (Table 1). This suggested that the period of planet TRAPPIST-1h may also satisfy a three-body resonance with TRAPPIST-1f and g. The six potential periods of TRAPPIST-1h that satisfy three-body relations with $1\leq p,q\leq 3$ are 18.766 d ($p=q=1$), 14.899 d ($p=1,q=2$), 39.026 d ($p=2,q=1$), 15.998 d ($p=2,q=3$), 13.941 d ($p=1,q=3$), and 25.345 d ($p=3,q=2$). We examined $\sim 1000$ hours of ground-based data taken prior to the Spitzer dataset[2] and found a lack of obvious additional transits at the expected times for all of these periods save 18.766 d. The period of 18.766 d corresponds to prior transit times in windows that were missed by the previous ground-based campaigns; hence, this was the only period that could not be ruled out. Furthermore, as this value is consistent with the period estimate of $20^{+15}_{-6}$ d based on the duration of the Spitzer transit, we had reason to believe it was the correct period for TRAPPIST-1h. To test this hypothesis, in our first transit search we simply fold the long cadence light curve at the four expected times of transit given this period and the single Spitzer transit time, finding evidence for a transiting planet at that period. Follow-up with detrended short cadence data confirms the transit-like shape of each of the four events and a depth consistent with that of TRAPPIST-1h (see Methods).

To prove the uniqueness of this detection, in a second analysis, we search the detrended K2 light curve after subtracting a transit model including all known transits of planets b–g, based on published ephemerides and planet parameters[2]. We use the photometric residuals as input to a Box-fitting Least Squares (BLS) algorithm (see Methods) to search for additional transit signals. In this search, we do not impose prior information on TRAPPIST-1h. We find a periodic signal at $\sim 18.77$ d with a transit center at $\mathrm{BJD}=2,457,756.39$, which matches the single transit observed by Spitzer[2].

Independently, we perform a joint instrumental/transit model fit to the data after subtracting a model for planets TRAPPIST-1b to g based on the Spitzer parameters (see Methods). We compute the relative likelihood (delta-chi squared) of a transit model with the best-fit Spitzer parameters of TRAPPIST-1h centered at every long cadence and sum the delta-chi squared values at the transit times corresponding to different periods in the range $1-50$ d. Strong peaks emerge at $18.766$ d and its aliases, corresponding to four transit-like events consistent with the parameters of TRAPPIST-1h at the times recovered in the previous searches.

We use the orbital period of TRAPPIST-1h determined in the previous step along with the parameters[2] of planets TRAPPIST-1b to g to determine whether a model including TRAPPIST-1h is favoured. This is achieved through Markov Chain Monte Carlo model fits with and without TRAPPIST-1h. We find a Bayes Factor of 90 in favour of a model that includes TRAPPIST-1h (see Methods), supporting the photometric detection of this seventh planet in the K2 dataset. The detection of TRAPPIST-1h is thus supported by 1) the three transit search analyses that recovered both the orbital phase from the Spitzer ephemeris[2] and the period of 18.766 d, 2) the Bayes Factor in favour of the 7-planet model and 3) the orbital period that is the exact value predicted by Laplace relations. Figure 1 shows the full light curve, the newfound transits of TRAPPIST-1h, as well as an update to the geometry of the orbits given the new orbital period. In Table 2 we report the properties of the planet derived in this study.

To characterize the three-body resonance, we use the transit timing data to identify $\phi$ (see Methods) for each set of three planets. A full transit-timing cycle has not elapsed within the data, so we cannot estimate the libration center for each $\phi$. However, we report in Table 1 the values of $\phi$ represented in the dataset. In the case of Jupiter’s satellites, $\phi=180^{\circ}$, but due to the complexity of the multi-planetary system, we make no predictions for TRAPPIST-1 at this time. Migration and damping simulations applied to Kepler-80[12] naturally predicted the values of the libration centers in that system; the measured values for TRAPPIST-1b–h call for future theoretical work interpreting them.

The resonant structure of the system suggests that orbital migration may have played a role in its formation. Embedded in gaseous planet-forming disks, planets growing above $\sim 1$ MMars create density perturbations that torque the planets’ orbits and trigger radial migration[13]. One model for the origin of low-mass planets found very close to their stars proposes that Mars- to Earth-sized planetary embryos form far from their stars and migrate inward[14]. The inner edge of the disk provides a migration barrier[15] such that planets pile up into chains of mean motion resonances[16, 17, 18]. This model matches the observed period ratio distribution of adjacent super-Earths[19] if the vast majority ($\sim 90\%$) of resonant chains become unstable and undergo a phase of giant impacts[20]. Some resonant chains do survive, and a handful of multiple-resonant super-Earth systems have indeed been characterized[21, 11]. The TRAPPIST-1 system may thus represent a pristine surviving chain of mean motion resonances. Given that TRAPPIST-1’s planet-forming disk was likely low in mass[22] and the planets themselves are low-mass, their migration was likely relatively slow. This may explain why TRAPPIST-1’s resonant chain is modestly less compact than chains in systems with more massive planets[21, 11, 20], which may have protected it from instability[23].

Given how close the planets orbit TRAPPIST-1, tidal interactions are likely to be important in the planets’ orbital evolution[24]. Tidal simulations of the system (see Methods) show that within a few Myr the eccentricity of each planet is damped to less than 0.01. Nonetheless, tidal heating is significant: all planets except TRAPPIST-1f and h have a tidal heat flux higher than Earth’s total heat flux.

The incident stellar flux on planet TRAPPIST-1h, 200 W m*-2*, is below the 300 W m*-2* required to sustain surface liquid water under a N2-CO2-H2O atmosphere (see Methods). To obtain the missing 100 W m*-2* from tidal heating would require a high eccentricity strictly incompatible with the orbits of the other planets. Our simulations show that the stellar input is also too low to sustain a thick CO2 atmosphere due to CO2 condensation. In particular, CO2 levels cannot exceed 100 ppm within a 1 bar N2 atmosphere. Alternatively, a liquid water ocean is possible under a layer of ice. The minimum thickness $h$ of this layer depends on the internal heat flux $\Phi_{int}$:

[TABLE]

Assuming the Earth’s current geothermal flux, a layer of 2.8 km (the mean depth of Earth’s oceans) would be necessary.

While the long spin-down times of ultra-cool dwarfs prevent derivation of a robust gyrochronology relation[25], the rotational period of TRAPPIST-1 can be used to derive a provisional age estimate for the system. Fourier analysis of the detrended K2 data (Fig. 2), which is visibly modulated by star spots, leads to the determination of a rotational period of ${\sim}3.3$ d for the host star (see Methods). This rotation corresponds to an angular momentum of about 1% that of the Sun. It is roughly in the middle of the period distribution of nearby late M dwarfs[26], suggesting an age in the range $3-8$ Gyr based on a star formation history that declines slightly with time[27]. Such an age is consistent with the star’s solar metallicity[1] and borderline thin disk/thick disk kinematics[28]. The amplitude of the modulation due to star spots and infrequent weak optical flares (0.26 d*-1* for peak fluxes above 1% of the continuum, 40 times less frequent than active M6-M9 dwarfs[29]) are consistent with a low-activity M8 star, also arguing in favor of a relatively old system. Near the end of the K2 campaign, a very energetic flare erupted, and was observed by Kepler. Full modeling of flares will be presented in a forthcoming paper.

The K2 observations of the TRAPPIST-1 star have enabled the detection of the orbital period of TRAPPIST-1h, continuing the pattern of Laplace resonance amongst adjacent triplets of planets. We search for but do not detect additional planets in the system. Compared with previous ground-based and Spitzer observations, the continuous coverage, high precision, and shorter wavelength of the K2 observations enable a robust estimate of the rotational period and flare activity of the star, motivating further study of the atmospheres and dynamical evolution of the planetary system.

0.1 References

0.2 Light curve preparation.

We use the package kadenza to generate a target pixel file (TPF) from the Campaign 12 raw data for TRAPPIST-1 (see Code Availability). In addition to EPIC ID 246199087, which corresponds to a standard-size postage stamp centered on TRAPPIST-1, the Kepler GO Office also made available a larger, $11\times 11$ pixel custom mask with a different ID (200164267), which we use for the purposes of this study. We manually select a rectangular $6\times 6$ pixel (24” $\times$ 24”) aperture centered on TRAPPIST-1 and use the median values of the remaining pixels to perform a column-by-column background subtraction. This process removes dark current and background sky signals while mitigating smear due to bright stars in the same CCD column as the target. Based on simple error function fits to the stellar image, we find that our aperture encloses $>0.9996$ of the flux of TRAPPIST-1h throughout the entire timeseries. We then detrend and perform photometry on the pixel-level data using two independent pipelines: EVEREST and a Gaussian process-based detrender. For computational expediency, we perform preliminary analyses on the long cadence data, following up on features of interest in the short cadence data (see below).

0.3 Light curve detrending – EVEREST.

The EVEREST K2 pipeline[7] uses a variant of pixel level decorrelation\citeMethodsDeming2015 (PLD) to remove instrumental systematics from stellar light curves. Given a stellar image spread out over a set of pixels $\{\vec{p_{i}}\}$, EVEREST regresses on polynomial functions of the fractional pixel fluxes, $\vec{p_{i}}/\sum_{j}\vec{p_{j}}$, identifying the linear combination of these that best fits instrumental signals present in the light curve. Because astrophysical signals (such as transits) are equally present in each of the pixels, whereas instrumental signals are spatially variable, PLD excels at removing instrumental noise while preserving astrophysical information. EVEREST uses a Gaussian process (GP) to model correlated astrophysical noise and employs an L2-regularized regression scheme to minimize overfitting.

Since PLD may overfit in the presence of bright contaminant sources in the target aperture[7], we manually inspect the TRAPPIST-1 K2 postage stamp and high resolution images of TRAPPIST-1 taken with the APO ARC 3.5m telescope in the SDSS z band to verify that there are no other targets brighter than the background level in our adopted aperture. Given the faint magnitude of TRAPPIST-1 in the Kepler band, we use the PLD vectors of 14 nearby bright stars (EPIC IDs 246177238, 246165150, 246211745, 246171759, 246127507, 246228828, 206392586, 246121678, 246229336, 246196866, 246217553, 246239441, and 246144695) generated using the same method as above to improve the signal-to-noise ratio of the instrumental model[8]. We further mask the data in the vicinity of all transits of planets TRAPPIST-1b–g and the potential transits of TRAPPIST-1h when computing the model to prevent transit overfitting. Finally, we divide the light curve into three roughly equal segments and detrend each separately to improve the predictive power of the model. Following these steps, we obtain a detrended light curve for TRAPPIST-1 with a 6-hr photometric precision\citeMethodsChristiansen2012,[7] of 281.3 ppm, a factor of 3 improvement on that of the raw light curve (884.4 ppm); see Fig. 1. Before the K2 observation, we estimated the $Kp$ magnitude of TRAPPIST-1 to be $17.2\pm 0.3$ based on a fit to a corrected blackbody spectrum. However, the photometric precision we achieve with EVEREST is inconsistent with a target dimmer than $Kp\sim 17$. Our detrending therefore suggests that the magnitude of TRAPPIST-1 in the Kepler band is $16<Kp<17$.

0.4 Light curve detrending – Gaussian process (GP) model.

Independently, we also detrend the data with a GP-based pipeline. To perform aperture photometry, we locate the star using a centroid fit and apply a circular top-hat aperture following the star’s centroid coordinates. We then use a GP model to remove the pointing drift systematics using an additive kernel with separate spatial, time and white noise components\citeMethodsaigrain2015, aigrain2016:

[TABLE]

[TABLE]

[TABLE]

where $x$ and $y$ are the pixel coordinates of the centroid, $t$ is the time of the observation, and the other variables ($A_{xy}$, $L_{x}$, $L_{y}$, $A_{t}$, $L_{t}$, $\sigma$) are hyperparameters of the GP model. We use the GEORGE package\citeMethodsambikasaran2014 in PYTHON to implement the GP model. To find the maximum likelihood hyperparameters we use a differential evolution algorithm\citeMethodsStorn1997, followed by a local optimization. This method was tested on magnitude $16-18$ stars observed in Campaign 10 of K2, and we use the results of those tests, and of previous GP applications to K2 data\citeMethodsaigrain2016, to inform our priors on the hyperparameters.

For the TRAPPIST-1 data, we use an iterative sigma-clipping method to remove outliers and prevent the time component from overfitting. This method has been previously used in the k2sc\citeMethodsaigrain2016 pipeline. First, using fiducial hyperparameter values based on analysis of a Campaign 10 target, we remove all measurements with residuals greater than $3\sigma$ from the mean GP prediction. With the remaining measurements, we update the hyperparameters by maximizing the GP likelihood. Using these parameters, we once again clip all $3\sigma$ outliers and maximize the GP likelihood using only the remaining measurements. The final detrending is calculated for all points, including outliers.

0.5 Photometric analysis I.

We fold the long cadence data on the dynamically predicted orbital period of 18.765 d for TRAPPIST-1h with a time of first transit constrained by the Spitzer observation, revealing a feature consistent with a transit in both the EVEREST and GP-based light curves. In order to confirm the planetary nature of this signal, we analyze the short cadence data. To this end, we use kadenza to generate a short cadence TPF of TRAPPIST-1 and detrend it in windows of 1.5–2 d centered on each of the four features using EVEREST. We use the PLD vectors of five bright stars observed in short cadence mode (EPIC IDs 245919787, 246011640, 246329409, 246331757, and 246375295) to aid in the detrending. When generating these light curves, we explicitly mask large flares so that these do not inform the fit. Following this procedure, we obtain binned 6-hr photometric precision of 266.6, 176.1, 243.4, and 219.3 ppm in each of the four windows. The short and long cadence data in these windows is shown in Supplementary Fig. 1. For transits 3 and 4, additional correction of the light curve is necessary, since the transit of TRAPPIST-1h coincides with a transit of TRAPPIST-1b (panel 3a) and a small flare (panel 4b). The transit of TRAPPIST-1b is subtracted out using a transit model\citeMethodsMandelAgol2002 with the Spitzer parameters and a mid-transit time determined from the data, yielding the light curve in panel 3b. The flare is fit using a 3-parameter flare model\citeMethodsDavenport2014 for stars observed with Kepler, yielding the light curve in panel 4b (see also Supplementary Fig. 2). In both cases, the transit of TRAPPIST-1h is visible in the residuals. In Supplementary Fig. 3 we show the folded short cadence data after accounting for TTVs, with a transit model based solely on the Spitzer parameters.

0.6 Photometric analysis II.

We use the long-cadence detrended light curve to perform a transit search with a box least-squares fitting algorithm (BLS)\citeMethodsKovacs:2002. We set the BLS to orbital periods ranging from 10 d to 50 d. The ratio of the transit duration over the planet orbital period is set between 0.0007 and 0.06 to include a wide range of orbital periods, eccentricities, and impact parameters for additional planets in the system. Due to Kepler’s 30-minute cadence, transits of planets orbiting TRAPPIST-1 appear significantly smeared out, which we take into account in our analysis. The highest peak in the periodogram corresponds to a signal with a 15.44-day period. Its origin stems from residuals in blended transits with TRAPPIST-1c and d and from one outlier in the data. No signal is seen at the two other epochs where a transit should have appeared, which confirms the signal being spurious. The next highest peak in the periodogram corresponds to a $\sim$18.77-day period and a transit center at $\mathrm{BJD}=2,457,756.39$ that is consistent with the single transit seen with Spitzer.

We then use a Markov Chain Monte Carlo (MCMC) algorithm previously described in the literature\citeMethodsGillon:2012 to derive the transit parameters of TRAPPIST-1h from the detrended light curve. Each photometric data point is attached to a conservative error bar that accounts for the uncertainties in the detrending process presented in the previous section. We impose normal priors in the MCMC fit on the orbital period, transit mid-time center and impact parameter for planets TRAPPIST-1b to g to the values recently published[2]. We further assume circular orbits for all planets[1, 2]. We also include normal priors for the stellar properties, which are $\mathcal{N}$(0.080, 0.0072) M*⊙* for the mass, $\mathcal{N}$(0.117, 0.0042) R*⊙* for the radius, $\mathcal{N}$(2,555, 852) K for the effective temperature and $\mathcal{N}$(0.04, 0.082) dex for the metallicity[2]. We use these stellar parameters to compute the quadratic limb-darkening coefficients $u_{1}$ and $u_{2}$ in the Kepler bandpass from theoretical tables\citeMethodsClaret:2011. In a first MCMC fit, we use a 7-planet model that includes all seven planets, with no prior information on the orbital period or $t_{0}$ of TRAPPIST-1h. In the second fit we employ a 6-planet model that excludes TRAPPIST-1h. We use the results from both MCMC fits to compute the Bayesian and Akaike Information Criteria (BIC and AIC, respectively) to determine which model is favoured. We find BIC values of 2888 and 2897 for the 7 and 6-planet models, respectively. This corresponds to a Bayes Factor $e^{(BIC_{1}-BIC_{2})/2}=90$ in favour of the 7-planet model. Similarly, we find AIC values of 2691 and 2725 for the 7 and 6-planet models, respectively. We perform a third MCMC fit to refine the transit parameters of TRAPPIST-1h. For this fit, we use as input data the K2 short-cadence data centered on the 4 transits of TRAPPIST-1h (Fig. 1) and the single transit light-curve previously obtained with Spitzer. This fit includes a model for TRAPPIST-1b and a flare that both affect the transit shape of TRAPPIST-1h in the K2 short-cadence data. This fit also allows for TTVs for the individual transit timings. We find photometric precisions of 365 ppm and $\sim$1100 ppm per 10 min for Spitzer and K2 respectively. We report the median and 1-sigma credible intervals of the posterior distribution functions for the transit parameters of TRAPPIST-1h in Table 2, along with the individual transit times.

0.7 Photometric analysis III.

In order to prevent the overfitting of transit features, we mask all transits of $b-h$ when detrending with EVEREST. However, this inevitably results in a lower detrending power during transits. A powerful alternative to the detrend-then-search method employed above is to simultaneously fit the instrumental and transit signals, without masking those features\citeMethodsForeman-Mackey2015. We therefore conduct a second, separate blind search on the EVEREST light curve specifically for TRAPPIST-1h. Given a raw light curve $\vec{y}$, a data covariance matrix $\mathbf{\Sigma}$, and a single-transit model $\vec{m}_{t_{0}}$ centered at $t=t_{0}$, the log likelihood of the transit fit is

[TABLE]

where $C$ is a constant. The data covariance matrix, $\mathbf{\Sigma}$, is the sum of the astrophysical covariance and the L2-regularized PLD covariance and is given by

[TABLE]

where $\mathbf{K}$ is the astrophysical covariance given by the EVEREST GP model, $\mathbf{X}$ is the matrix of PLD regressors (the design matrix), and $\mathbf{\Lambda}$ is the prior covariance of the PLD weights (the regularization matrix) which we obtain by cross-validation[8]. Since the transit shape and duration of TRAPPIST-1h are known[2], the only free parameter in the search was $t_{0}$, the time of transit. We therefore evaluate $\vec{m}_{t_{0}}$ multiple times, centering the transit model at each long cadence and computing the likelihood of the transit model fit as a function of cadence number. We then subtract these values from the log likelihood of the data with no transit model ($\vec{m}_{t_{0}}=0$) and multiply by 2 to get the delta-chi squared ($\Delta\chi^{2}$) metric, which measures the decrease in the $\chi^{2}$ value of the fitted light curve for a transit of TRAPPIST-1h centered at each cadence. Finally, we also compute $\Delta\chi^{2}$ conditioned on the known “true” transit depth of TRAPPIST-1h, $d_{0}=0.00352\pm 0.000326$:

[TABLE]

where

[TABLE]

is the maximum likelihood depth of the transit model and

[TABLE]

is the variance of the depth estimate. Positive peaks in $\Delta\chi^{2}$ indicate features that are well described by the transit model, while positive peaks in $\Delta\chi^{2}_{\mathrm{cond}}$ reveal features that are well described by the transit model with depth $d=d_{0}$. In Supplementary Fig. 4 we show the two $\Delta\chi^{2}$ metrics across the full TRAPPIST-1 light curve after subtracting a transit model for planets $b-g$ based on their Spitzer parameters. The strongest features in the $\Delta\chi^{2}$ plot (top) are flares, as these can be fitted out with an inverted transit model. When conditioning on the true depth of TRAPPIST-1h (bottom), the significance of most of the flare features decreases, revealing the four peaks of TRAPPIST-1h (red arrows).

In order to assess the robustness of our detection, we compute the total $\Delta\chi^{2}_{\mathrm{cond}}$ as a function of orbital period. Starting from the time of transit of TRAPPIST-1h in the Spitzer dataset, we compute the times of transit in the K2 light curve for 500,000 values of the orbital period evenly spaced between 1 and 50 d and sum the values of $\Delta\chi^{2}_{\mathrm{cond}}$ at each transit time to produce the total $\Delta\chi^{2}_{\mathrm{P}}$. We sum these in two different ways. First, we linearly interpolate the grid of $\Delta\chi^{2}_{\mathrm{cond}}$ to each transit time to obtain $\Delta\chi^{2}_{\mathrm{P}}$ for perfectly periodic transits. Next, to allow for TTVs of up to one hour, we take the largest value of $\Delta\chi^{2}_{\mathrm{cond}}$ in the vicinity of each transit time and sum them for each period. We adopted a tolerance of 2 cadences, corresponding to maximum TTVs of 1 hr. Our results are shown in Supplementary Fig. 5, where we plot the periodic $\Delta\chi^{2}_{\mathrm{P}}$ (top) and the $\Delta\chi^{2}_{\mathrm{P}}$ allowing for TTVs (bottom). In both cases, the period of TRAPPIST-1h (18.766 d) and its $1/2$ and $1/3$ period aliases emerge as the three strongest peaks. The peak at 18.766 d is the strongest signal in the period range constrained by the Spitzer transit[2] and confirms our detection of TRAPPIST-1h.

0.8 Three-body angles.

The mean longitude of a planet with orbital period $P$ is an angular variable that progresses at a constant rate with respect to time $t$, which is measured from the time the planet passes a given reference direction:

[TABLE]

with $\lambda$ measured in degrees. For transiting planets, the reference direction is taken as the plane perpendicular to the observer’s line of sight, as the planet is progressing towards the transiting configuration. We assume orbits with negligible eccentricities[11], for which $\lambda=90^{\circ}$ at transit mid-time, so that we may write

[TABLE]

for each planet, where $T_{n}$ is the time of transit of the $n^{\mathrm{th}}$ planet. For a 3-body resonance $(p,q)$, we may therefore express the three-body angle as

[TABLE]

The state of a $\phi$ value is assessed when individual transit times of three planets are taken near each other. For instance, for the transit times[2] $T_{f}=7662.18747$, $T_{g}=7665.35151$, $T_{h}=7662.55463$, and $(p,q)=(1,1)$, we compute $\phi=177.4^{\circ}$ at $t=7664$ ($\mathrm{BJD}-2,450,000$).

0.9 Tidal simulations.

Tidal interactions with the star are important for all 7 planets. We perform N-body simulations of the system including an equilibrium tidal dissipation formalism\citeMethodseggleton98,leconte10 using the Mercury-T code[Bolmont2015]. We use orbital parameters from the discovery paper[2] and a period of TRAPPIST-1h of 18.765 d (a near 2:3 resonance configuration with planet $g$). We consider the planets’ spins to be tidally synchronized with small obliquities. This is justified because even if the age of the system is 400 Myr (a lower estimate for the age of TRAPPIST-1), planetary tides would have had time to synchronize the spins\citeMethodsleconte2015 (even when atmospheric tides are accounted for). We test different initial eccentricities and different values for the planets dissipation factors (from 0.01 to 10 times the Earth’s value\citeMethodsDeSurgyLaskar1997).

Our simulations show that the planets’ orbital eccentricities are likely to be low. In just a few Myr all eccentricities decrease to below 0.01 for dissipation factors $\geq 0.1$ times the Earth’s value. Due to planet-planet interactions, the eccentricities do not decrease to zero but instead reach an equilibrium value determined by the competition between tidal damping and planet-planet eccentricity excitation\citeMethodsBolmont2013. All planets stay in resonance during the evolution towards tidal equilibrium. The small equilibrium eccentricities are sufficient to generate significant tidal heating. Assuming the TRAPPIST-1 planets have a tidal dissipation equal to the Earth’s, we find that TRAPPIST-1b might have tidal heat flux similar to Io’s\citeMethodsSpencer2000 ($\sim 3$ W m*-2*), with peaks at more than 10 W m*-2* (corresponding to $\sim 10^{4}$ TW) when the eccentricity oscillation is at its maximum. Planets $c$ through $e$ have tidal heat fluxes higher than Earth’s internal (primarily radioactive) heat flux\citeMethodsDavies2010 ($\sim$0.08 W m*-2*) but lower than Io’s. TRAPPIST-1f, g and h have a tidal heat flux inferior to Earth’s. Supplementary Fig. 6 shows a possible snapshot of the system’s evolution over the course of 40 yr. This very high flux could plausibly generate intense volcanism on the surfaces of the inner planets, with potential consequences for their internal structures.

0.10 Planet habitability.

We calculate the minimum stellar flux required for liquid water with the LMD 1D/3D Global Climate Model\citeMethodsTurbet2016 using a synthetic spectrum of TRAPPIST-1 based on its reported $T_{\mathrm{eff}}$, $\log g$, metallicity, and bolometric luminosity[2], obtaining a value of 300 $\mathrm{W/m^{2}}$, which is 100 $\mathrm{W/m^{2}}$ higher than the planet’s present-day instellation. Our results are in agreement with habitable zone boundaries computed for a 3000 K star\citeMethodsKopparapu2013. Assuming zero albedo, we find that the equilibrium temperature of TRAPPIST-1h is 169$\pm 4$ K.

Whether or not TRAPPIST-1h presently hosts an atmosphere is unclear. Given its radius measurement and a range of possible compositions (from pure water ice to pure iron), the mass of TRAPPIST-1h is likely in the range 0.067–0.863 $M_{\oplus}$; if TRAPPIST-1h has an Earth-like composition, this value is 0.33 $M_{\oplus}$\citeMethodsSotin2007. Assuming TRAPPIST-1h migrated to its current location quickly, the planet’s low surface gravity could have led to vigorous hydrodynamic escape of a primordial atmosphere in the first few 100 Myr after its formation, since at that time TRAPPIST-1 was significantly brighter and TRAPPIST-1h would have been interior to the habitable zone\citeMethodsLuger2015a,Luger2015b. The presence of a subsequently outgassed atmosphere, however, cannot be ruled out.

In theory, the surface of TRAPPIST-1h could harbor liquid water under such an outgassed atmosphere if it is H2-rich. Atmospheres made of H2\citeMethodsStevenson1999,Pierrehumbert_Gaidos2011, N2-H2\citeMethodsWordsworth_Pierrehumbert2013 or CO2-H2\citeMethodsRamirez_Kaltenegger2017,Wordsworth_2016 have been shown to provide a sufficient greenhouse effect and internal heat blanketing for even lower instellation levels. Unless prevented by high altitude clouds, transit spectroscopy with the Hubble Space Telescope may be able to reveal or rule out such H2-rich atmospheres.

0.11 Stellar variability.

To establish the rotation period of TRAPPIST-1, we clean the long cadence EVEREST light curve (Fig. 2) of remaining outliers, transits, and flares. We iteratively fit and remove low frequencies that remain in the light curve from the detrending process. To extract the rotation period, we calculate the Discrete Fourier Transform (DFT) following a previously established method,\citeMethodsKurtz85 and proceed to fit the detected frequency with non-linear least-squares. The rotation period is determined to be $3.30\pm 0.14$ d.

To determine the occurrence rate of stellar flares, we again take the EVEREST detrended long cadence light curve and remove transits and outliers. We require flares to have peak emission 1% above the normalized continuum flux and two consecutive signals above the continuum. We detect a total of 19 flare events, corresponding to an average rate of 0.26 d*-1*. A more refined determination of flaring rate and energies requires analysis of the short cadence data, which will be presented in a forthcoming paper.

0.12 Code availability.

The kadenza code we use to generate pseudo-target pixel files for all K2 targets was downloaded from https://doi.org/10.5281/zenodo.344973. The code we use to generate and analyze the EVEREST light curves for TRAPPIST-1 is openly available at https://github.com/rodluger/trappist1. A static version of the repository has been archived at https://doi.org/10.5281/zenodo.437548.

0.13 Data availability.

The $K2$ raw cadence data used in this study is available for download at https://archive.stsci.edu/missions/k2/c12_raw_cadence_data/. The pseudo-target pixel files for TRAPPIST-1 and its neighboring stars generated with kadenza are archived at https://doi.org/10.5281/zenodo.437876. The detrended TRAPPIST-1 long cadence light curve and segments of the short cadence light curve in the vicinity of the transits of TRAPPIST-1h are available at https://doi.org/10.5281/zenodo.437548. All other data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.

\bibliographystyleMethods

naturemag \bibliographyMethodsk2_trappist1_methods

\supplement