Temperate Earth-sized planets transiting a nearby ultracool dwarf star

Spectral type, parallax, and age of the star

TRAPPIST-1 = 2MASS J23062928-0502285 was discovered in 2000 by a search for nearby ultracool dwarfs based on photometric criteria³¹ and identified as a high proper-motion (μ_α=0.89”, μ_δ=-0.42”), moderately active (logL_Hα/L_bol=-4.61), M7.5 dwarf at ~11 parsecs. Subsequent studies converged to a spectral type of M8.0±0.5^11,32, while confirming a moderate level of activity typical of similar spectral type stars in the vicinity of the Sun^{15, 33, 34}. The spectral classification was checked by comparing a low-resolution (R~150) near-infrared spectrum of the star¹³ obtained with the SpeX spectrograph³⁵ mounted on the 3m NASA Infrared Telescope Facility to several spectral type standards, and the best-fit was obtained with the spectrum of the M8-type standard LHS 132 (Extended Data Fig. 3.a). The CTIOPI project reported the star's trigonometric parallax as π=82.6±2.6 mas¹², which translates to a distance of 12.1±0.4 parsecs. High-resolution optical spectroscopy failed to detect significant absorption at the 6708 Å Li line³⁶, inferring that the object is not a very young brown dwarf but rather a very-low-mass main-sequence star. This is in agreement with its thick disk kinematics³⁶, its relatively slow rotation (vsini=6±2 km.s^-1)¹⁵, its moderate activity, and its reported photometric stability¹⁶ that all point to an age of at least 500 Myrs¹³.

Metallicity of the star

We obtained new near-infrared (0.9-2.5 μm) spectroscopy of TRAPPIST-1 with the SpeX spectrograph on the night of 2015 November 18 (UT) during clear conditions and 0.8" seeing at K-band. We used the cross-dispersed mode and 0.3"x15" slit aligned at the parallactic angle, to acquire moderate-resolution data (λ/Δλ ≈ 2000) with dispersion 3.6 Å pixel^-1 covering the spectral range 0.9-2.5 μm in seven orders. Ten exposures of 300s each were obtained over an airmass range of 1.14-1.17, followed by observations of the A0V star 67 Aqr (V = 6.41) at an airmass of 1.19 for telluric and flux calibration, as well as internal lamp exposures. Data were reduced using the SpeXtool package version 4.04^37,38.The reduced spectrum has a median signal-to-noise of 300 in the 2.17-2.35 μm region, which is shown in Extended Data Fig. 3.b, with metallicity-sensitive NaI (2.206, 2.209 μm) and CaI (2.261, 2.263, 2.266 μm) atomic features labeled. We measured equivalent widths of these features and the H₂O-K2 index defined in Rojas-Ayala et al. (2012)³⁹, and used the mid- and late-M dwarf metallicity calibration of Mann et al. (2014)⁴⁰ to determine [Fe/H] = 0.04±0.02 (measurement) ±0.07 (systematic) for TRAPPIST-1. The quadratic sum of the two errors resulted in our final measurement [Fe/H] = 0.04±0.08.

Basic parameters of the star

A recent study¹³ derived a luminosity L_★= 0.000525±0.000036 L_⊙ for TRAPPIST-1, using as input data the trigonometric parallax and VRI magnitudes as measured by the CTIOPI project¹², 2MASS JHK magnitudes⁴¹, WISE W123 magnitudes⁴², an optical spectrum measured with KPNO/R-C Spec⁴³, and a near-IR spectrum measured by SpeX/Prism. Using this luminosity and the constraint on the age to be >500 Myr as input, the authors of the study¹³ derived from evolutionary model isochrones and the Stefan-Boltzmann law the following values for the stellar mass, radius, and effective temperature: M_★= 0.082±0.009 M_⊙, R_★= 0.116±0.004 R_⊙, and T_eff = 2557±64K, respectively. To account for the uncertainties coming from the assumptions and details of the evolutionary models, we performed a new determination of these three basic parameters with recent solar metallicity evolutionary model isochrones that consistently couple atmosphere and interior structures⁴⁴. We obtained M_★= 0.089 M_⊙, R_★= 0.112 R_⊙, and T_eff = 2615K. We then added the difference between the two determinations quadratically to the errors of Filippazzo et al.¹³, adopting finally M_★= 0.082±0.011 M_⊙, R_★= 0.116±0.006 R_⊙, and T_eff = 2555±85K. The normal distributions corresponding to these values and errors were assumed as prior probability distribution functions in the Bayesian analysis of our photometric data (see below).

Possible binarity of the star

High-resolution imaging from the ground^45,46,47 and from space with HST⁴⁸ discarded the existence of a companion down to an angular distance of 0.1", corresponding to a projected physical distance of 1.2 au at 12 parsecs, in good agreement with the reported stability of the radial velocity of the star at the ~10 ms^-1 level over a week⁴⁹ and at the ~150 ms^-1 level over ~10 weeks⁵⁰. We performed spectral binary template fitting⁵¹ to the IRTF/SpeX spectroscopy, and statistically reject the presence of an L- or T-type brown dwarf companion that would be visible in a blended-light spectrum. TRAPPIST-1 can thus be considered in all probability as an isolated star.

Upper magnitude limits on a background eclipsing binary (BEB)

We measured the J2000 equatorial coordinates of TRAPPIST-1 in the 2015 TRAPPIST images, using 29 stars from the UCAC2 catalog⁵² and the Pulkovo Observatory Izmccd astrometric software⁵³. We obtained coordinates of RA=23h06m30.34s, DEC=-05d02m36.44s. Due to the high proper motion of TRAPPIST-1 of ~1"/yr, the presence of a possible background object could be assessed by examining this exact position in several previous images taken from the POSS⁵⁴ (1953) and 2MASS⁴¹ (1998) image catalogs. No possible additional source was detected at its position in any of these images. The faintest stars detected at other positions in the 2MASS images have J-band magnitudes of ~17. We adopt this value as an absolute lower threshold for the J-band magnitude of a background source blended with TRAPPIST-1 in our TRAPPIST 2015 images. TRAPPIST-1 has a J-band magnitude of 11.35⁴², and the achromaticity of the transits of TRAPPIST-1b as observed from 0.85 to 2.09 μm imposes that, if they originated from a BEB, the latter would have to be a very red object with a spectral type similar to TRAPPIST-1. Combining these two facts, the BEB scenario would require an unphysical eclipse depth >100% in the photometric bands probed by our observations to match the ~0.8% depths measured after dilution by the light of TRAPPIST-1. The BEB scenario is thus firmly discarded.

Photometric observations and analysis

The TRAPPIST^8,55 observations in which the transits were detected consisted of 12,295 exposures of 55-s gathered with a thermoelectrically-cooled 2kx2k CCD camera (field of view of 22' x 22', pixel scale of 0.65"). Most of the observations were obtained through an 'I+z' filter having a transmittance >90% from 750 nm to beyond 1100 nm, the red end of the effective bandpass as defined by the spectral response of the CCD. Based on the spectral efficiency model for TRAPPIST and an optical spectrum of a spectroscopic standard M8V star (VB10), we compute an effective wavelength of 885±5 nm for these observations. For the nights of 20 November and 19 December 2015, the target was close to the full Moon and the observations were performed in the Sloan z' filter to minimize the background. After a standard pre-reduction (bias, dark, flat-field correction), the TRAPPIST automatic pipeline extracted the stellar fluxes from the images using the DAOPHOT aperture photometry software⁵⁶ for eight different apertures. A careful selection of both the photometric aperture size and of stable comparison stars was then manually performed to obtain the most accurate differential light curves of the target.

Photometric follow-up observations were performed with the HAWK-I near-IR imager⁵⁷ on the ESO 8m Very Large Telescope (Chile), with the HFOSC optical spectro-imager⁵⁸ on the 2m Himalayan Chandra Telescope (India), and with the WFCAM Wide-Field infrared CAMera⁵⁹ located at the prime focus of the 3.8m UKIRT telescope (Hawaii).

The VLT/HAWK-I observations of a transit of planet b were performed during the night of 8 November 2015. HAWK-I is composed of four Hawaii 2RG 2048 x 2048 pixels detectors (pixel scale = 0.106"). Its total field of view on the sky is 7.5'x7.5'. The transit was observed through the narrowband filter NB2090 (λ=2.095μm, width=0.020μm). 185 exposures composed of 17 integrations of 1.7s each were acquired during the run in stare mode, i.e. without applying a jitter pattern. After standard calibration of the images, stellar flux measurement was performed by aperture photometry¹⁴.

The HCT/HFOSC observations of a transit of TRAPPIST-1b were performed on 18 November 2015. The imager in the HFOSC CCD detector is an array of 2048x2048 pixels corresponding to a field of view of 10'x10' on-sky (pixel scale=0.3"). The observations consisted of 104 exposures, each 20-s taken in stare mode and in the I filter, centered on the expected transit time. After a standard calibration of these images and their photometric reduction with DAOPHOT, differential photometry was performed. We estimate the effective wavelength of these observations to be 840±20 nm, based on the spectral response of HFOSC and an optical spectrum of the M8V standard star VB10.

The UKIRT/WFCAM observations of 2 transits of planet b and 1 transit of planet c consisted of 3 runs of 4 hrs each, performed on 5, 6 and 8 December 2015 in J-band. WFCAM is composed of four HgCdTe detectors of 2048x2048 pixels each, with a pixel scale of 0.4" resulting in a field of view of 13.65'x13.65' for each detector. On 5 December 2015, 1365 exposures composed of 3 integrations of 2s each were performed in stare mode. For the runs on 6 and 8 December 2015, respectively, 1181 and 1142 exposures composed of five 1s-exposures were performed, again in stare mode and using the same pointing as 5 Dec 2015. Differential aperture photometry was performed with DAOPHOT on all calibrated images.

Global analysis of the photometry

We inferred the parameters of the three detected planets transiting TRAPPIST-1 from the analysis of their transit light curves (Extended Data Fig. 1 and Extended Data Table 1) with an adaptive Markov-Chain Monte Carlo (MCMC) code¹⁴. We converted each UT time of mid-exposure to the BJD_TDB time system⁶⁰. The model assumed for each light curve was composed of the eclipse model of Mandel & Agol (2002)⁶¹ multiplied by a baseline model aiming to represent the other astrophysical and instrumental mechanisms able to produce photometric variations. Assuming the same baseline model for all light curves, and minimizing the Bayesian Information Criterion (BIC)⁶², we selected a second-order time polynomial as a baseline model to represent the curvature of the light curves due to the differential extinction and the low-frequency variability of the star, and added an instrumental model composed of a second-order polynomial function of the positions and widths of the stellar images.

Stellar metallicity, effective temperature, mass, and radius were four free parameters in the MCMC for which prior probability distribution functions (PDFs) were selected as input. Here, the normal distributions N(0.04, 0.08²) dex, N(2555, 85²) K, N(0.082, 0.011²) M_⊙, and N(0.114, 0.006²) R_⊙ were assumed based on a priori knowledge of the stellar properties (see Basic parameters of the star). Circular orbits were assumed for all transiting objects. For each of them, the additional free parameters in the MCMC included (1) the transit depth dF defined as (R_p/R_★)², with R_p and R_★ the planetary and stellar radii, respectively (2) the transit impact parameter b =a cosi/R_★, with a and i the planet's semi-major axis and orbital inclination, respectively (3) the orbital period P, (4) the transit width W defined as (P R_★/a) [(1+R_p/R_★)² - b²]^1/2/π, and (5) the mid-transit time (time of inferior conjunction) T₀. Uniform prior distributions were assumed for each of these free parameters. At each step of the MCMC, values for R_p, a, i, were computed from the values for the transit and stellar parameters, as well as values for the irradiation of the planet in Earth units and for its equilibrium temperatures assuming a Bond albedo of 0 and 0.75, respectively. A quadratic limb-darkening law⁶² was assumed for the star. For each bandpass, values and errors for the limb-darkening coefficients u₁ and u₂ were derived from the Claret & Bloemen tables (2011)⁶³ (see Extended Data Table 2), and the corresponding normal distributions were used as prior PDFs in the MCMC. u₁ and u₂ were free parameters under the control of these PDFs in the MCMC.

We divided our analysis into three phases. The first phase focused on the two inner planets, for which the period is firmly determined. A circular orbit was assumed for both planets. All transit light curves of the two planets were used as input data for this first phase, except the TRAPPIST light curve of 11 December 2015 for which the transit of planet c is blended with a transit of planet d. A preliminary MCMC analysis composed of 1 chain of 50,000 steps was first performed to estimate the need to rescale the photometric errors¹⁴. Then a longer MCMC analysis was performed, composed of 5 chains of 100,000 steps whose convergence was checked using the statistical test of Gelman & Rubin⁶⁴. The parameters derived from this analysis for the star and its two inner planets are shown in Table 1. We performed a similar analysis assuming a uniform prior PDF for the stellar radius to derive the value of the stellar density constrained only by the transit photometry⁶⁵. It resulted in a stellar density of ${49.3}_{-8.3}^{+4.1}{\rho }_{\odot },$ in excellent agreement with the density of 55.3±12.1 ρ_⊙ derived from the a priori knowledge of the star, bringing thus a further validation of the planetary origin of the transit signals.

In the second phase of our analysis, we performed 11 global MCMC analyses of all transit light curves, each of them consisting of one chain of 50,000 steps and corresponding to one of the possible values of the period of TRAPPIST-1d (see Table 1) for which a circular orbit was assumed. We then repeated the 11 analyses under the assumption of an eccentric orbit for planet d. We used the medians of the BIC posterior distributions to compare the relative posterior probability of each orbital model through the formula $P1/P2=e^{(BI{C}_2-BI{C}_1)/2}$. The resulting relative probabilities are given in Extended Data Table 3. It shows that our data favour (relative probability > 10%) a circular orbit and an orbital period between 10.4 and 36.4 days, the most likely period being 18.4 days.

In the final phase, we performed individual analyses of the light curves to measure the mid-eclipse time of each transit to support future TTV studies of the system^{22, 66}. The resulting timings are shown in Extended Data Table 4. They do not reveal any significant TTV signal, which is not surprising given the amplitude of the expected periodicity departures (see below) combined with the limited timing precision of the TRAPPIST photometry.

Extended Data Fig. 1 and 2 show the raw and detrended light curves, respectively, and for each of these the best-fit eclipse + baseline model is overplotted in the figures. The phased-time detrended light curves are shown for each planet and bandpass in Fig. 1.

Photometric variability of the star

The TRAPPIST dataset was used to assess the photometric variability of the star at ~900 nm. On the timescale of a few hours, corresponding to the typical duration of our observing runs, the star appears to be relatively stable, except for the transits and for 4 sharp, low-amplitude (1 to a few %) increases of brightness followed by exponential-type decreases to the normal levels within 10-15 minutes (Extended Data Fig. 4) that we attribute to flares⁶⁷. The low amplitude and inferred low frequency (1/60 hr^-1) of these flares is consistent with the reported low level of activity of the star^15,33,34, strengthening the inference that the system is not young.

To assess the lower-frequency variability of TRAPPIST-1, we built its global differential light curve in the I+z filter, using four stable stars of similar brightness in the TRAPPIST images as comparison stars. We filtered out the flares, transits, and measurements taken in cloudy conditions to create the resulting light curve consisting of 12,081 photometric measurements. It is compared to the light curve for the comparison star 2MASSJ23063445-0507511 in Extended Data Fig. 5.a. It clearly shows some variability at the few % level, which is consistent with previous photometric results obtained in the I-band¹⁶. A Lomb-Scargle (LS) periodogram⁶⁸ analysis of the light curve, filtered out of low-frequency variations and differential extinction by division of the best-fit 4th-order polynomial in time and airmass, reveals a power excess with a period of 1.4 days (see Extended Data Fig. 5.b). Cutting the light curves in two, and in four in a second test, and performing a LS analysis of each fraction, revealed a power excess at ~1.4 days for all of them, supporting a genuine periodic signal of astrophysical origin. Associating it with the stellar rotation period, the resulting equatorial rotation speed of 4.1 km.s^-1 (assuming R_★= 0.114 R_⊙) is consistent with the literature measurement¹⁵ vsini=6±2 km.s^-1, making this association physically meaningful. Based on the scatter of the peak values obtained in the LS analyses of the light curve fractions, we estimate the error bar on the rotation period of 1.40 days to be 0.05 day. In summary, the photometric variability of the star appears thus to be dominated by the rotation and evolution of photospheric inhomogeneities (spots) combined with rare flares.

Dynamics of the system

We computed the tidal circularization timescales⁶⁹ $t_{circ}=\frac{2PQ}{63\pi }\times \frac{M_p}{M_{\ast }}\times {\left(\frac{a}{R_p}\right)}^5$ of the three planets, assuming planetary masses M_p ranging from 0.45 Earth masses (pure ice composition) to 3 Earth masses (pure iron composition)²² and a tidal quality factor⁷⁰ Q of 100, corresponding to the maximum value derived for terrestrial planets and satellites of the solar system⁷⁰. For planets b and c, the computed values range from 22 Myr to 145 Myr and from 177 Myr to 1.1 Gyr, respectively. Taking into account that the system is apparently not very young and that the orbits have weak mutual perturbations as they are not close to any mean-motion resonance, our assumption of circular orbits for the two inner planets is thus reasonable. On the other hand, the same computations result in values ranging from a few to tens of Gyr for planet d, making a significant orbital eccentricity possible from a tidal theory perspective. Still, a nearly-circular orbit for this outer planet is still a reasonable hypothesis when considering the strong anticorrelation of orbital eccentricity and multiplicity of planets detected by radial velocities⁷¹, and is favoured by our global analysis of the transit photometry (see above).

We used the Mercury software package⁷² to assess the dynamical stability of the system over 10,000 years for all possible periods of planet d. Instabilities appeared in our simulations only for the unlikely scenarios of planet d on a significantly eccentric (e≥0.4) 4.5 or 5.2-day orbit.

To assess the potential of the TTV method^{24, 66} to measure the masses of the planets, we integrated the dynamical evolution of the system at high sampling over 2 years, assuming Earth-masses for the three planets and an 18.4 day circular orbit for TRAPPIST-1d. These simulations resulted in TTV amplitudes of several tens of seconds and led us to conclude that with an intensive transit monitoring campaign with instruments able to reach timing precisions of a few tens of seconds (e.g. VLT/HAWK-I, UKIRT/WFCAM, see Extended Table 4), it should be possible to constrain the planetary masses.

Estimation of the suitability of the planets for detailed atmospheric characterization

We estimated the typical signal amplitude in transit transmission spectroscopy for all the transiting exoplanets with a size equal to or smaller than the mini-Neptune GJ1214b⁷³. We computed it as $2{\mathrm{R}}_{\text{p}}{\mathrm{h}}_{\text{eff}}/{\text{R}}_{\ast }^2$ where R_p is the planetary radius, h_eff the effective atmospheric height (i.e., the extent of the atmospheric annulus), and R_∗ is the stellar radius. The effective atmospheric height is directly proportional to the atmospheric scale height, H=kT/(μg) where k is Boltzmann’s constant, T the atmospheric temperature, μ the atmospheric mean molecular mass, and g the surface gravity. The ratio h_eff/H for a clear atmosphere^24,74 is typically between 6 and 10 and is thus strongly dependent on the presence of clouds and the spectral resolution and range covered. Our estimates, presented in Fig. 2, are based on an h_eff/H ratio of 7 and the conservative assumption of a volatile-dominated atmosphere (μ =20) with a Bond albedo of 0.3. All other parameters for the planets were derived from exoplanets.org⁷⁵. As illustration, the maximum transit depth variations projected under those assumptions for GJ1214b are of ~250 ppm, in agreement with independent simulations⁷⁶.

We also derived for the same planets sample typical SNRs in transit transmission spectroscopy from the ratio of our computed signal amplitudes over the square root of the flux (determined from the J-band magnitudes of the host stars). The SNRs of TRAPPIST-1’s planets in transmission are expected to range between 0.22 and 0.55 times GJ 1214b’s, implying that these planets are well-suited for atmospheric studies with HST/WFC3 similar to those previously targeting GJ1214b^76,77.

Based on literature simulations for terrestrial planets²⁴, we estimated that the characterization of planets b, c, and d should require up to 70, 90, and 270 hrs of in-transit observations with JWST, respectively, and yield the atmospheric temperatures with relative uncertainties below 15% and the abundances within a factor of 4. Assuming the atmospheres of TRAPPIST-1’s planets are not depleted and do not harbor a high-altitude cloud deck, JWST should notably yield constraints on the abundances of molecules with large absorption bands such as H₂O, CO₂, CH₄, CO, and O₃ if their abundances are ≥10-ppm level.

We also assessed the potential of the cross-correlation technique⁷⁸ to constrain the atmospheric properties of the TRAPPIST-1 planets, following a literature formalism⁷⁹. We inferred that detecting O₂ in TRAPPIST-1’s planets would require up to 80 transit observations with one of the next generation giant ground-based telescopes. Taking in account the limited fraction of transits visible at low airmass, such an endeavor could be reached in 5 to 15 years.

Code availability

Equivalent widths and H₂O-K2 index measurements in the SpeX spectra were made using the IDL program created by A. Mann and distributed at http://github.com/awmann/metal. The conversion of the UT times of the photometric measurements to the BJD_TDB system was performed using the online program created by J. Eastman and distributed at http://astroutils.astronomy.ohio-state.edu/time/utc2bjd.html. IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. The MCMC software used to analyse the photometric data is a custom Fortran 90 code that can be obtained upon request to the first author.