EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 613 (May 2018) A&A, 613 (2018) A68 Full HTML
Issue
A&A
Volume 613, May 2018
Article Number A68
Number of page(s) 21
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201732233
Published online 04 June 2018

© ESO 2018

1 Introduction

The TRAPPIST-1 system, which harbours seven Earth-size exoplanets orbiting an ultra-cool star (Gillon et al. 2017), represents a fascinating setting to study the formation and evolution of tightly packed small planet systems. While the TRAPPIST-1 planet sizes are all known to better than 5%, their densities suffer from significant uncertainty (between 28% and 95%) because of loose constraints on planet masses. This critically impacts in turn our knowledge of the planetary interiors, formation pathway (Ormel et al. 2017; Unterborn et al. 2018) and long-term stability of the system. So far, most exoplanet masses have been measured using the radial-velocity technique. But because of the TRAPPIST-1 faintness (V = 19), precise constraints on Earth-mass planets are beyond the reach of existing spectrographs.

Thankfully, the resonant chain formed by the planetary septet (Luger et al. 2017b) dramatically increases the exchange of torque at each planet conjunction, resulting in transit timing variations (TTVs; Holman 2005; Agol et al. 2005) that are well above our demonstrated noise limit for this system. Presently, the TTV approach thus represents the only avenue to characterise the physical properties of the system. The TRAPPIST-1 system shows dynamical similarities to Kepler-90 system (Cabrera et al. 2014) which contains also seven planets and resonances conditions between pairs of them.

Planetary masses published in the TRAPPIST-1 discovery paper (Gillon et al. 2017) were bearing conservative uncertainties because the different techniques used by the authors suggested a non-monotonous parameter space with the absence of a single global minimum. Subsequent studies have adequately invoked the requirement for long-term stability to refine these masses further (Quarles et al. 2017; Tamayo et al. 2017), but the parameter space allowed by this additional constraint may still be too large to precisely identify the planet physical properties. The recent K2 observations of TRAPPIST-1 (Luger et al. 2017b) enabled another team to compute updated masses for the system using the K2 data combined to archival data (Wang et al. 2017). Their approach relies on the TTVFast algorithm (Deck et al. 2014), which uses low-order symplectic co-ordinates and an approximate scheme for finding transit times to increase efficiency. It is however unclear from that paper how the correlations between parameters are taken into account and how comprehensive the search of the parameter space is. Only a full benchmarking of this approach with more accurate integrators for this specific system could validate their results.

In the present paper, we have used a novel approach that combines an efficient exploration of the parameter space based on a genetic algorithm with an accurate N-body integration scheme. The associated complexity being compensated by more computing resources. The philosophy of this approach could be considered “brute force” but still represents a useful avenue to appreciate the degeneracy of the problem without doubting the accuracy of the numerical integration scheme.

2 Observations

2.1 Published data

This study is based on 284 transit timings obtained between September 17, 2015 and March 27, 2017 through the TRAPPIST and SPECULOOS collaboration. The input data for our transit-timing analysis includes 107 transits of planet b, 72 of c, 35 of d, 28 of e, 19 of f, 16 of g, and 7 of h. In addition to the TRAPPIST-1 transit timings already presented in the literature (Gillon et al. 2016, 2017; de Wit et al. 2016; Luger et al. 2017b), we have included new data from the Spitzer Space Telescope (PID 12126, 12130, and 13067) and from Kepler and K2 (PID 12046). Transit timing uncertainties range from 8 s to 6.5 min with a median precision across our dataset of 55 s. The analysis of the K2 data is presented below while the new Spitzer data obtained between February and March 2017 are presented in a separate publication (Delrez et al. 2018). We include the full list of transit timings used in this work in Tables A.1 through G.1.

2.2 K2 short-cadence photometry

For the purpose of this analysis we have included the transit timings derived from the K2 photometry (Howell et al. 2014), which observed TRAPPIST-1 during Campaign 12 (Luger et al. 2017b). We detail in the following the data reduction of this dataset. We used the K2’s pipeline-calibrated short-cadence target pixel files (TPF) that includes the correct timestamps. The K2 TPF TRAPPIST-1 (EPIC ID 246199087) aperture is a 9 × 10 postage stamp centred on the target star, with 1-minute cadence intervals. We performed the photometric reduction by applying a centroiding algorithm to find the (x, y) position of the PSF centre in each frame. We used a circular top-hat aperture, centred on the fitted PSF centre of each frame, to sum up the flux. We find this method to produce a better photometric precision compared to apertures with fixed positions.

The raw lightcurve contains significant correlated noise, primarily from instrumental systematics due to K2’s periodic roll angle drift and stellar variability. We have accounted for these systematic sources using a Gaussian-processes (GP) method, relying on the fact that the instrumental noise is correlated with the satellite’s roll angle drift (and thus also the (x, y) position of the target) and that the stellar variability has a much longer timescale than the transits. We used an additive kernel with separate spatial, time and white noise components (Aigrain et al. 2015, 2016; Luger et al. 2017b):

where x and y are the pixel co-ordinates of the centroid, t is the time of the observation, and the other variables (Axy, Lx, Ly, At, Lt, σ) are hyperparameters in the GP model (Aigrain et al. 2016). We used the GEORGE package (Ambikasaran et al. 2016) to implement the GP model. We used a differential evolution algorithm (Storn & Price 1997) followed by a local optimisation to find the maximum-likelihood hyperparameters.

We optimised the hyperparameters and detrended the lightcurve in three stages. In each stage, using the hyperparameter values from the previous stage (starting with manually chosen values), we fitted the GP regression and flag all points further than 3σ from the mean as outliers. The next GP regression and hyperparameter optimisation is performed by excluding the outlier values. Points in and around transits are not included in the fit. Due to the large numbers of points in the short cadence lightcurve, only a random subset of the points is used to perform each detrending and optimisation step to render the computation less intensive. We achieved a final RMS of 349 ppm per 6 h (excluding in-transit points and flares). Out of the entire dataset, we discarded only two transits because of a low signal-to-noise ratio at the following BJDTDB: 7795.706 and 7799.721. Both eclipses correspond to transits of planet d.

We then used a Markov chain Monte Carlo (MCMC) algorithm previously described in the literature (Gillon et al. 2012) to derive the individual transit timings of TRAPPIST-1b, c, d, e, f, g, and h from the detrended K2 lightcurve. Each photometric data point is attached to a conservative error bar that accounts for the uncertainties in the detrending process presented above. We have imposed normal priors in the MCMC fit on the orbital period, transit mid-time centre and impact parameter for all planets to the published values (Gillon et al. 2017). We computed the quadratic limb-darkening coefficients u1 and u2 in the Kepler bandpass from theoretical tables (Claret & Bloemen 2011) and employ the transit model of Mandel & Agol (2002) for our fits. We derived the TTVs directly from our MCMC fit for all TRAPPIST-1 planets. We report the median and 1-sigma credible intervals of the posterior distribution functions for the 124 K2 transit timings in Tables A.1–G.1. The resulting K2 short-cadence stacked lightcurves are shown in Fig. 1.

thumbnail Fig. 1

Folded short-cadence lightcurves extracted from K2 data, corrected for their corrected for their TTVs. Relative fluxes are shown for each planet and are shifted vertically for clarity. The short-cadence data are binned to produce the coloured points, with five cadences per point. The white points are further binned, with ten points taken from the folded curve per bin.

3 Methodology

3.1 Dynamical modelling

Dynamical studies of the TRAPPIST-1 system are challenging due to the 7-planet, Laplace resonance-chain architecture with tight pair period ratios of 8:5, 5:3, 3:2, 3:2, 4:3, and 3:2. This configuration requires computationally expensive orbital integrations over a large parameter space. The observed transit times are distributed over more than 550 days, corresponding to over 370 orbits of the innermost planet b. In order to accommodate the order of one billion time steps needed in an MCMC method to match the timing precisions with the timespan of observations, we used graphics processing units (GPUs) and the GPU N-body code GENGA (Grimm & Stadel 2014). We calculated the orbits of all seven planets and determine the TTVs through MCMC techniques. GENGA uses a hybrid symplectic integration scheme (Chambers 1999) to run many instances of planetary systems in parallel on the same GPU. We have extended the GENGA code by implementing a GPU version of the parallel Differential Evolution MCMC (DEMCMC) technique (Braak 2006; Vrugt et al. 2009). DEMCMC deploys multiple Markov chains simultaneously that efficiently sample the highly correlated multi-dimensional parameter space that is typical of TTV inversion problems (Mills et al. 2016). In addition we have modified the DEMCMC sampling method, such that it works more efficiently with the large correlations impacting the masses, semi-major axes and mean anomalies of the different planets.

3.2 Transit timing calculations

Following (Fabrycky 2010), we defined the X–Y plane as the plane of the sky, while planets transit in front of the star at positive values of the Z co-ordinate. The mid-transit times are found by minimising the value of the function (4)

which can be solved by setting the next time step of our numerical integrator to (5)

with (6)

The quantities xi and yi are the astrocentric co-ordinates of the planet i.

We useda pre-checker in the integrator to determine if a transit is expected to happen during the next time step. This will happen if the value of gi moves from a negative value to a positive value during the next time step, and if zi > 0. In addition, we used the conditions |g∕ġ| < 3.5 ⋅ dt and rsky < (Rstar + Ri) + |vi|⋅dt to refine the pre-checker, where |rsky| = |(xi, yi)| is the radial co-ordinate on the sky plane, |vi| is the norm of the velocity, and Ri the radius ofplanet i; and R⋆ the radius of the star.

Since the integration is performed with a symplectic integrator, the co-ordinates of the position and velocity of a planet are not simultaneously known, which leads to a small error in the calculation of g. If a transit occurs very close to a time step, then it can happen that the transit is reported in both successive time steps with a slightly different mid-transit time. But when the time step is small enough, this error can be safely neglected. Also, in highly eccentric orbits, the described pre-checker may not work properly, because g changes too quickly between each time step. We thus restricted ourselves in this work to eccentricities smaller that 0.2 and used the fourth-order integrator scheme with a time step of 0.08 days. When the pre-checker has detected a transit candidate, then all planets are integrated with a Bulirsh–Stoer direct N-body method for a time step and the Eq. (5) is iterated until δt is smaller than a tolerance value. A transit is reported if rsky < Rstar + Ri.

3.3 Orbital parameter search

To determine the best orbital parameters, we used a DEMCMC method (Braak 2006; Vrugt et al. 2009). We used N parallel Markov chains, where each chain consists of a d dimensional parameter vector xi. To update the population of N chains, each x is updated by generating a proposal (7)

with i≠j, i≠ k, j≠ k, and a small perturbation e. The proposal is accepted with a probability p = min(1, π(xp)∕π(xi)). When the proposal is accepted, then xi is replaced by xp in the next generation; otherwise the state remains unchanged. In each 30th generation, we set γ = 0.98 to allow jumps between multimodal solutions (Braak 2006). In addition, we set γ = 0.01 and x i = xl≠i during the burn-in phase to eliminate outliers. Alternatively, we also tested the affine invariant ensemble walker MCMC method (Goodman & Weare 2010; Foreman-Mackey et al. 2013), which yields a comparable performance.

For the probability density function, π(xi) we used (8)

where the running index t refers to the transit epoch, Tcalc are the calculated mid transit times, Tobs are the observed mid-transit times and σ are the observation uncertainties. Using Eq. (8), we rewrite the acceptance probability as (9)

where τ is the MCMC sampling “temperature”. In this work, we used values for τ between 1000 and 1. Using a large value for τ increases the acceptance probability of the next DEMCMC step and allows the walkers to explore more easily a large parameter space, but the obtained probability distribution does not correspond to the likelihood. Resampling the so-obtained likelihood region with a smaller value of τ, and startingfrom the median values of the previous runs, refines the sampling more accurately. Using different values of τ in an iterative order allows us to sample a large parameter space, with good accuracy in the most likely region. The median and the standard deviation are calculated with a value of τ = 1. Figure 4 shows the obtained posterior probability distribution of the masses and eccentricities. We note that using a value for τ < 1 in Eq. (9) would lead to smaller standard deviations.

According to (Goździewski et al. 2016), we used the following fitting parameters for the orbital elements:

with the period Pi, the time of the first transit Ti, the start time of the simulation t, the mean anomaly at the first transit , the mean anomaly Mi, the eccentricity ei, the argument of perihelion ωi and the Jacobi mass for each planet i. The Jacobi mass of planet i includes also the masses of all objects with a smaller semi-major axis. We used the square root of the eccentricity in the parameters xi and yi to favour low eccentricity solutions. We set the longitude of the ascending node Ω i to zero and the inclination of all planets to π∕2, which allows us to calculate through the true anomaly at the transit 1.

Assuming coplanarity is motivated by the fact that the standard deviation of the derived inclinations for all seven planets with respect to the sky plane is 0.08° only (Gillon et al. 2017). If the longitudes of nodes were distributed randomly on the sky, then the probability that all planets transit would be very small (most observers would see only one planet transit if this were the case). Thus, the three-dimensional mutual inclination can be constrained by simulating the angular momentum vectors of the planets drawn from an 3D Gaussian inclination distribution of width σθ, allowing the density of the star to vary, ρ*, and determining which set of parameters matches the transit durations most precisely from observers drawn from random locations on the sphere. This yields a constraint on the three dimensional inclination distributions of σθ < 0.3° at the 90% confidence level (Luger et al. 2017a). TTVs depend very weakly on the mutual inclinations of the planets (Nesvorný & Vokrouhlický 2014), and since these planets are constrained to be coplanar to a high degree based upon the argument in Luger et al. (2017b), our model is justified in neglecting mutual inclinations of the planets for our dynamical analysis.

As initial conditions of the DEMCMC parameter search, we have randomly distributed the parameters in the range mi ∈[0, 6 × 10−6M⊙], ei ∈ [0, 0.05]2, ωi ∈ [0, 2π] and Mi ∈ [0, 2π]. We assumed uniform priors on the parameters, but restrict the eccentricity to e < 0.2.

A difficulty in sampling the orbital parameters of TTVs is, that there can exist strong correlations between some of the parameters, especially between m and a or between e and ω. The correlation between m and a can be explained, because the period Pi in Eq. (10) depends on all the masses of the more inner planets. The correlation between e and ω is caused by the resonance configuration, and a change in one of these parameters must be compensated by the others to get a similar time of closest approach between the planets. When the different walkers of the DEMCMC method are spread out over a large region in the parameter space, then the acceptance ratio of the DEMCMC steps will get very low, due to inaccurate guesses of the semi-major axes and mean anomalies.

3.3.1 Sub-step optimisation

The DEMCMC algorithm generates proposal values of the parameters, by linearly interpolating between two accepted values, but often the parameters in the optimisation problem show a non-linear dependency. This means that the DEMCMC approach always deviates from an optimal choice of the new proposal steps. Since the value of χ2 is very sensitive to small perturbations of the semi-major axis and the mean anomaly, inaccurate guesses on these parameters will lead to dramatic high values of χ2, and the proposal step is very unlikely to be accepted. The consequence is an acceptance ratio going towards zero as the parameters begin to populate a broader region in the parameter space. To improve this issue, we introduced a sub-step optimisation scheme to find the optimal values for the semi-major axis a and the mean anomaly M. The sub-step scheme is applied after each DEMCMC step, and has to be performed for each planet in a serial way. When the value of a or M is changed for only a single planet, then the χ2 shows a parabolic behaviour, which enables the use of a quadratic estimator to find its optimal values x, based on three guesses x1, x2, and x3 as follows:

where x means either a or M, and yj means the values of χ2 at locations xj. The cost of the described sub-step optimisation scheme is, that three times as many walkers are needed to generate the values y1, y2, and y3, and each DEMCMC step has to be followed by 14 sub-steps of computing the TTVs to adjust a and M for each planet. But even if this scheme is expensive to compute, it allows us to achieve an acceptance rate that remains >20–30% for a much larger number of DEMCMC steps. The best χ2 obtained is 342, details are listed in Table 1. The evolution of the masses and the autocorrelation functions of 5000 DEMCMC steps with the described sub-step sampling for 100 chains are shown in Fig. 5, which shows efficient convergence of the chains.

Table 1

Number of observations, number of degrees of freedom Ndof = Nobservations − Nparameters, χ2 from Eq. (8) and reduced χ2 = χ2∕Ndof, for all planets separately and for all planets together.

3.3.2 Independent analysis

We also carried out an independent transit timing analysis with a new version of TTVFast (Deck et al. 2014) which utilises a novel symplectic integrator based upon a fast Kepler solver (Wisdom & Hernandez 2015; Hernandez & Bertschinger 2015). The integrator uses a time step of 0.05 days, assumes a plane-parallel geometry, and alternates between drifts and universal Kepler steps between pairs of planets. A drift in Cartesian co-ordinatesis defined as an update of some or all positions assuming constant velocities. The initial conditions are constructed with Jacobi co-ordinates, and the integration uses Cartesian co-ordinates in a center-of-mass frame. The transit times are found by tracking the projected sky position and velocity, and finding when the dot product changes sign, using Eqs. 4-6. The transit centre is found by bracketing and interpolating the time steps (Deck et al. 2014). This yields timing precisions of better than a few seconds, which is sufficient to model the data given the observational uncertainties. We modelled the transit times with this code, and obtained identical masses for the maximum likelihood (within the uncertainties), as well as broadly consistent eccentricities. Since this analysis was carried out with a different code in a different language (Julia) with a different integration technique, the fact that a similar maximum a posteriori likelihood gives confidence that we have found a unique solution for the mass ratios of this planet system.

4 Results

The TTV values of 1000 MCMC posterior samples with τ = 1 are shown in Fig. 2, in comparison to the observed transit times. The TTV residuals for each transit are shown in Fig. 3. We compute the planet densities ρp independently from the stellar mass, by using the planetary radii- and mass-ratio posteriors and , along with the well constrained stellar density ρ⋆ determined photometrically from the Spitzer dataset (Seager & Mallén-Ornelas 2003). The planet density then is (Jontof-Hutter et al. 2014). Using the stellar density from the photometry is valid in our case because the planets eccentricities are found to be small. Determining planetary densities using this approach effectively removes any inaccuracy from the stellar models and improves our constraints on the planetary interiors. To transform our results into physical masses and radii, we use the most recent stellar mass estimate of 0.089 ± 0.007 M⊙ (Van Grootel et al. 2018).

Our resulting posterior distribution functions for the masses and eccentricities for all seven planets are shown in Fig. 4. To perform the search over a large parameter space, different sampling “temperatures” τ are used in an iterated order. Through our extensive exploration of the parameter space we find that the masses and eccentricities for all planets are reasonably constrained (3–9% for mass, and 6–25% for eccentricity). Table 2 summarises the planetary physical parameters (mass and radius ratios) while Table 3 displays the planets’ orbital parameters. A full posterior distribution between all mutual pairs of parameters is shown in Fig. 6.

5 Comparison with other studies

At the time of writing, we are aware of a preprint by Wang et al. (2017) that re-analysed our initial Spitzer and ground-based datasets, and added in K2 transit timing measurements. Wang et al. (2017) utilises TTVFast (Deck et al. 2014), which also employs a symplectic integrator. Our analysis improves upon their initial analysis in several ways:

  • 1.

    We account for the correlations in mass and radius when computing the constraints on the planets in the mass–radius plane. The density is better constrained than either the mass or radius, leading to a strong correlation between planetary mass and radius (Fig. 10), which parallels the isocomposition contours, and so our approach should improve the constraint upon composition.

  • 2.

    We utilise a new set of Spitzer transit times, which, compared with the K2 times, are superior in timing precision by about a factor of two, cover a longer time duration (>100 days), and are less affected by stellar variability.

  • 3.

    Our fast, parallel GPU integration scheme coupled with a parallel MCMC algorithm allows a thorough exploration of parameter space. The reduced χ2 of our best fits are near unity, while Wang et al. (2017) utilise models with high reduced χ2 in their analysis.

These improvements over the Wang et al. (2017) analysis should lead to more robust and accurate constraints upon the compositions, masses and radii of the planets.

thumbnail Fig. 2

Calculated TTVs for all planets and for 1000 different MCMC samples (grey lines). Measured transit times with the corresponding uncertainties are indicated by coloured symbols, according to the used telescope. A detailed list of all transits is given in Appendix. The differences between the solutions reflect the distribution shown in Fig. 4 for τ = 1. The two panels at the bottom show a zoomed region for planet b and c.
Table 2

Updated masses, radii (Delrez et al. 2018), correlation coefficients cm,r of masses and radii as well as the densities and the surface gravity of all seven planets.

Table 3

Median and standard deviation σ of the semi-major axis a, eccentricity e, argument of periapsis ω, and mean anomaly M of the seven planets.

thumbnail Fig. 3

Residuals of the TTVs shown in Fig. 2. The transit data index corresponds to the column index of Tables A.1–G.1.

6 Dynamical properties of the TRAPPIST-1 system

In this section, we use the results of our dynamical modelling of the system to investigate the degeneracies in the planetary parameters and the stability of the system over long time-scales.

6.1 Tackling degeneracies

Degeneracies commonly plague TTV inversion problems (Agol & Fabrycky 2017), in particular, between mass and eccentricity (Lithwick et al. 2012). However, Fig. 4 shows that the eccentricities and masses are well constrained for all seven TRAPPIST-1 planets. A combination of the 1) high-precision Spitzer photometry, 2) K2’s 80-day long quasi-continuous coverage and the 3) resolved high-frequency component of the TTV pattern known as “chopping” (Holman et al. 2010, Deck & Agol 2015 all contribute to mitigate the mass-eccentricity degeneracy. The chopping patterns for all planets except for h are detected in the data (Fig. 2), while their periodicities encode the timespan between successive conjunctions of pairs of successive planets whose amplitudes yield the masses of adjacent perturbing planets.

6.2 Long-term stability

We study in the following the temporal evolution of the eccentricities and arguments of periapsis, respectively, for all seven planets over 10 Myr. The discussion on the evolution of the argument of periapsis and eccentricities for all planets are based on Figs. 7 and 8. We show in Fig. 9 the temporal evolution of the Laplace three body resonant angles ϕ. Eccentricities remain below 0.025 for all planets and show a very regular evolution for 2 Myr (i.e. 487 millions orbits of planet b and close to 30 millions of planet h). After that, the eccentricities evolve more irregularly, but are still bound to low values. The arguments of periapsis, however, exhibit different behaviours. All planets show an average precession of ≈2π/300 yr. Planets d, e, f, and g show a more sporadic evolution, with ω undergoing periodic phases of fast precession during which yr. This behaviour arises from the small eccentricities and strong mutual gravitational perturbations between theplanets. We also study the evolution with time of the three-body Laplace angles ϕ. The initial values of ϕ agree well with reported values(Luger et al. 2017b) but after 2 Myr, the resonance chain is perturbed. This behaviour reflects well the evolution of the eccentricities and is an indication that the initial conditions we are using are not known sufficiently accurately to assess resonances over very-long timescales. The exact solution should survive for several Gyrs in the resonance chain as the TRAPPIST system is comprised of a suite of resonance chains (Luger et al. 2017b).

Tides could in addition be particularly important in this closely packed system (Luger et al. 2017b). Tides damp eccentricity and stabilisethe dynamical perturbations. Preliminary results using the Mercury-T code (Bolmont et al. 2015) show that a small tidal dissipation factor at the level 1% that of the Earth (Neron de Surgy & Laskar 1997) would generate a tidal heat flux of 10–20 W m−2, which is significantly higher than Io’s tidal heat flux (Spencer et al. 2000) of 3 W m−2. Given the estimated age of the system, most of the eccentricity of planet b at least should have been damped.

thumbnail Fig. 4

Posterior probability distribution of the mass and eccentricities of all seven planets, for different MCMC sampling “temperatures” τ, assuming a stellar mass of M⋆ = 0.09 M⊙ (Van Grootel et al. 2018). The contours correspond to significance levels of 68%, 95% and 99% for τ = 1. Sampling “temperatures” greater than one are used to cover a larger parameter space.

7 The nature of the TRAPPIST-1 planets

Our improved masses and densities show that TRAPPIST-1 c and e likely have largely rocky interiors, while planets b, d, f, g, and h require envelopes of volatiles in the form of thick atmospheres, oceans, or ice, in most cases with water mass fractions ≲5% (Fig. 10). These values are close to the ones inferred by another recent study based on mass–radius modelling (Unterborn et al. 2018). It is also consistent with accretion from embryos that grew past the snow line and migrated inwards through a region of rocky planet growth.

7.1 Planetary interiors

We calculate the theoretical mass—radius curves shown in Fig. 10 for interior layers of rock, ice, and water oceanfollowing the thermodynamic model of Dorn et al. (2017). This model uses the equation of state (EoS) model for iron by Bouchet et al. (2013). The silicate-mantle model of Connolly (2009) is used to compute equilibrium mineralogy and density profiles for a given bulk mantle composition. For the water layers, we follow Vazan et al. (2013), using a quotidian equation of state (QEOS) for low pressure conditions and the tabulated EoS from Seager et al. (2007) for pressures above 44.3 GPa. We assume an adiabatic temperature profile within core, mantle and water layers.

The mass–radius diagram in Fig. 10 compares the TRAPPIST-1 planets with theoretical mass–radiusrelations calculated with published models (Dorn et al. 2017). We estimate the probability pvolatiles of each planet to be volatile-rich by comparing masses and radii to the idealised composition of Fe/Mg = 0 and Mg/Si = 1.02 (Fig. 10), which represents a lower bound on bulk density for a purely rocky planet. Planets b, d, f, g, and h very likely contain volatile-rich layers (with pvolatiles of at least 0.96, 0.99, 0.66, 1, and 0.71, respectively). Volatile-rich layers could comprise atmosphere, oceans, and/or ice layers. In contrast, planets c and e may be rocky (with pvolatiles of at least 0.24 and <0.01, respectively).

The comparison of masses and radii with the idealised interior end-member of mwater∕M = 0.05 (fixed surface temperature of 200 K, Fig. 10, blue solid line) suggest that all planets, except planet b and d, do very likely not contain more than 5% mass fraction in condensed water. For planets b and d, we thereby estimate a probability of 0.7 and 0.5 to contain less than 5% water mass fraction. However, because the runaway greenhouse limit for tidally locked planets lies near the orbit of planet d (Kopparapu et al. 2016; Turbet et al. 2018), the large amount of volatiles needed to explain the radius of the most irradiated planet b is likely to reside in the atmosphere (possibly as a supercritical fluid), reducing the mass fraction needed in a condensed phase.

thumbnail Fig. 5

Left panels: evolution of the masses of all planets of the entire ensemble of 100 walker chains in blue, and the evolution of a single chain in green. Right panels: auto-correlation function of the masses for the ensemble and a single chain.

7.2 Limits for the possible atmosphere-scenarios

We use the LMD-G one-dimensional cloud-free numerical climate model (Wordsworth et al. 2010)to simulate the vertical temperature profiles of the seven TRAPPIST-1 planets. Calculations are performed using a synthetic spectrum of TRAPPIST-1 and molecular spectroscopic properties from Turbet et al. (2018). We assume atmospheric compositions that range from pure H2, H2 –CH4, H2 –H2O to pure CO2. We further assume a volume molecular mixing ratio of 5 × 10−4 for methane and 1 × 10−3 for water, tomatch solar abundances in C and O (Asplund et al. 2009). We finally assume a core composition of Fe/Mg = 0.75 and Mg/Si = 1.02 (Unterborn et al. 2018).

For each planet, atmospheric composition and a wide range of surface pressures (from 10 mbar to 103 bar, see Table 4), we decompose the thermal structure of the atmosphere into 500 log-spaced layers in altitude co-ordinates. We estimate the transit radii of the planets, as measured by Spitzer in the 4.5 μm IRAC band, solving radiative transfer equations including molecular absorption, Rayleigh scattering, and various other sources of continuua, that are Collision Induced Absorptions (CIA) and/or far line wing broadening, when needed and available. Radiative transfer equations are solved through the 500 layers in spherical geometry (Waldmann et al. 2015) to determine the effective transit radius of a given configuration in the Spitzer band.

A fit was found when the surface pressure resides at the nominal transit radius of the Spitzer observations and thus corresponds to the maximum surface pressure. It is maximal in the sense that any reservoir of volatiles at the surface would yield a higher core radius and reduce the mass of the atmosphere needed to match the observed radius.

For atmospheres with a higher mean molecular weight, the inferred pressures are too large for a perfect gas approximation. Assuming that the mass of the atmosphere is much lower than the core, the surface pressure Psurf can be estimated by integrating the hydrostatic equation, which yields: (16)

where Ptransit is the pressure at the transit radius Rtransit, Rcore the radius at the solid surface, the mean atmospheric scale height, k is Boltzmann’s constant, T the surface temperature, μ the mean molecular weight, and mH the mass of the hydrogen atom. This relation demonstrates that the surface pressure increases exponentially with ; and is why, at the low temperatures expected for the planets beyond d, it is difficult to explain the observed radii with an enriched atmosphere above a bare core (with a standard Earth-like composition) without unrealistically large quantities of gas.

For the colder, low density planets (f, g, and h), explaining the radius with only a CO2 atmosphere is difficult due to the small pressure scale height and the fact that CO2 should inevitably collapse on the surface beyond the orbit of planet g (Turbet et al. 2018). We acknowledge that these various results could be challenged with a significantly different rock composition or thermal state of the planets.

Planet b, however, is located beyond the runaway greenhouse limit for tidally locked planets (Kopparapu et al. 2016; Turbet et al. 2018)and could potentially reach – with a thick water vapour atmosphere – a surface temperature up to 2000 K (Kopparapu et al. 2013). Assuming more realistic mean temperatures of 750–1500 K, the above estimate yields pressures of water vapour of the order of 101 –104 bar, which could explain its relatively low density (assuming Ptransit = 20 mbar). As such, TRAPPIST-1 b is the only planet above the runaway greenhouse limit which seems to require volatiles.

Given the density constraints and assuming a standard rock composition (Unterborn et al. 2018), planets b to g cannot accommodate H2 -dominated atmospheresthicker than a few bars. Within these assumptions and considering the expected intense atmospheric escape around TRAPPIST-1 (Bolmont et al. 2017; Wheatley et al. 2017; Bourrier et al. 2017a,b), the lifetime of such atmospheres would be very limited, making this scenario rather unlikely. For heavier molecules, the surface pressure needed to match a given radius varies roughly exponentially with the mean molecular weight, which imposes enormous surface pressures for which a more detailed equation of state would be needed.

thumbnail Fig. 6

Posterior probability distribution between all pairs of parameters, the masses m, the semi-major axes a, the eccentricities e, the arguments of perihelion ω, and the mean longitudes l = Ω + ω + M of all seven planets. It is showing strong correlations between many pairs of parameters, especially between m and a and between e and ω.

7.3 Mass-ratios, densities and irradiation

The mass ofa celestial object is its most fundamental property. We now compare the masses of the TRAPPIST-1 planets and place them into wider context. Exoplanet discoveries are heavily biased towards single Sun-like stars. TRAPPIST-1 provides a glimpse of what results around stars, and within disks, that are one order of magnitude lighter than the norm.

Figure 11 displays a mass-ratio versus period diagram comparing the TRAPPIST-1 system to other exoplanets (where the orbital period is used to separate various sub-populations). We also push the comparison to include the planets of the Solar system, and the principal moons of Jupiter. TRAPPIST-1’s planets cover mass-ratio a range 10−4–10−5, which is shared with the sub-Neptune and super-Earth exoplanets population that orbits Sun-like stars. The similarity may suggest a similar formation mechanism, or at minimum a comparable scaling in protoplanetary disk mass. Notably, sub-Neptunes and super-Earths are the most abundant planet types for Sun-like stars (Mayor et al. 2011; Howard et al. 2012). This region also encompasses the Galilean satellite system, and thus reflects a host or satellite configuration that spans over three orders of magnitude in mass, like has been noticed in the literature (Canup & Ward 2006).

The irradiation of the TRAPPIST-1 planets plays an important role in their evolution. It is thus also insightful to compare the TRAPPIST-1 masses and incident flux in the context of the currently known exoplanet population. Figure 12 shows a focus on planets receiving irradiation that spans Mars to Venus and 0.1–4.5 M⊕. The upper mass limit is set to 4.5 M⊕ because this corresponds to 1.5 R⊕, an indicative limit for rocky worlds (Rogers 2015; Fulton et al. 2017). The lower mass limit was arbitrarily set to 0.1 M⊕, corresponding to Mars, to represent objects that have difficulties retaining an atmosphere. It is interesting to note that TRAPPIST-1d & e are the only transiting exoplanets in this region. The closest other known transiting planets are LHS 1140 b (0.4, 6.65) (Dittmann et al. 2017) & Kepler-138b (0.64, 2.33) (Jontof-Hutter et al. 2015), and both of these have much larger uncertainties on their densities (27% and 92% respectively).

We complete this section by comparing the density of the TRAPPIST-1 planets with their irradiation level in Fig. 13. We note that the density trend vs. irradiation increases until TRAPPIST-1e, where it peaks and then decreases for the outer planets. One object stands out of that pattern, TRAPPIST-1d, which interestingly, has a similar bulk density as the Moon. These aspects will be particularly relevant to better understand the formation pathway of the TRAPPIST-1 system.

thumbnail Fig. 7

Temporal evolution of the eccentricities of all planets. The eccentricities are well constrained to values below 0.015 and show a regular behaviour for 2 Myr. After that, the systems show irregularities, caused by small uncertainties in the initial conditions. We also show the relative energy error and the relative angular momentum error of the integrator.
thumbnail Fig. 8

Temporal evolution of ω, showing a fast precession and strong mutual orbital perturbations between the planets.
thumbnail Fig. 9

Temporal evolution of the three body resonant angle ϕi = pλi−1 − (p + q)λi + qλi+1, where the values of p and q for consecutive triples of planets are (2,3), (1,2), 2,3), (1,2) and (1,1) (Luger et al. 2017b). After 2 Myr, the systems show irregularities, caused by small uncertainties in the initial conditions.

7.4 Migration and composition

The combination of a planetary system’s orbital architecture and the planets’ densities can constrain where the planets grew (or at least their feeding zones) as well as their orbital histories (Raymond et al. 2008).

Based on cosmochemical abundances (Lodders 2003), one would expect objects that condensed past the snow line to contain a significant fraction (up to 50%) of ice. However, other effects can reduce planets’ water contents, such as desiccation of constituent planetesimals by short-lived radionuclides (Grimm & McSween 1993), water loss during giant impacts between embryos (Genda & Abe 2005), and heating during the very rapid growth expected around low-mass stars (Lissauer 2007; Raymond et al. 2007). In addition, embryos do not migrate inwards through an empty expanse towards the inner parts of the disk. Rather, they likely migrate through a region in which rocky material is already growing (Izidoro et al. 2014). The water contents of migrating icy embryos are likely to be diluted by impacts with rocky embryos. In some situations, migrating icy embryos can stimulate the growth of large rocky embryos, creating a large density contrast between neighbouring planets (Raymond et al. 2006).

Gravitational interactions with the gaseous disk cause ~Earth-mass planetary embryos to migrate, usually inwards (Ward 1997; Baruteau et al. 2014). Current models invoke two steps in the growth of these embryos. First, as dust coagulates and drifts through the disk (Güttler et al. 2010; Birnstiel et al. 2012), 10–100 km scale planetesimals form by a hydrodynamical instability (such as the streaming instability; Youdin & Goodman 2005; Johansen et al. 2014) in regions where the local solids-to-gas ratio is sufficiently high (Dra̧żkowska & Dullemond 2014; Carrera et al. 2015). These conditions are expected to be met first beyond the snow line (Armitage et al. 2016; Schoonenberg & Ormel 2017). Next, the largest planetesimals grow rapidly by accreting inwards-drifting pebbles (Ormel & Klahr 2010; Lambrechts & Johansen 2012). Models of pebble accretion find that large embryos preferentially grow beyond the snow line (Morbidelli et al. 2015; Ormel et al. 2017). However, embryos’ growth is self-limiting. When an embryo reaches a critical mass, it creates a pressure bump in the gas disk exterior to its orbit, which traps drifting pebbles and shuts off pebble accretion (Lambrechts et al. 2014). This critical mass depends on the disk structure but was calculated by Ormel et al. (2017) to be ~0.7 M⊕ for the case of TRAPPIST-1 (for a specific disk model), close to the actual planet masses. In their model, embryos form sequentially, reaching this critical mass before migrating inwards.

The resonant structure of the TRAPPIST-1 system (Luger et al. 2017b) is a telltale sign of orbital migration (Terquem & Papaloizou 2007; Ogihara & Ida 2009). The fact that all seven planets form a single resonant chain indicates that the entire system migrated in concert (Cossou et al. 2014; Izidoro et al. 2017). Indeed, orbital solutions generated by disk-driven migration have been shown to be more stable than other solutions (Tamayo et al. 2017). Whereas most resonant systems are likely to be unstable (Izidoro et al. 2017; Matsumoto et al. 2012), the TRAPPIST-1 can be interpreted as a system that underwent a relatively slow migration creating a long-lived resonant system.

thumbnail Fig. 10

Mass–radius diagram for the TRAPPIST-1 planets, Earth, and Venus. Curves trace idealised compositions of rocky and water-rich interiors (surface temperature fixed to 200 K). Median values are highlighted by a black dot. Coloured contours correspond to significance levels of 68% and 95% for each planet. The interiors are calculated with model II of Dorn et al. (2017). Rocky interiors are composed of Fe, Si, Mg, and O, assuming different bulk ratios of Fe/Mg and Mg/Si. U17 refers to Unterborn et al. (2018). Equilibrium temperatures for each planet are indicated by the coloured contours.
Table 4

Planetary characteristics of the TRAPPIST-1 planets derived from our 1D numerical climate simulations.

thumbnail Fig. 11

Ratio between the mass of a planet or moon and that of its host as a function of orbital period. We represent the known population of exoplanets (from exoplanet.eu, Schneider et al. 2011) from three different detection methods. The Solar system is also highlighted. The TRAPPIST-1 system, like the Galilean moons of Jupiter, share a similar parameter space with the sub-Neptune and super-Earth population. Orbital periods are used to reveal the various sub-populations.
thumbnail Fig. 12

Masses and the incident fluxes received by the TRAPPIST-1 planets, compared to other exoplanets found by the transit method (yellow), via the radial-velocity technique (blue), and to the terrestrial worlds of the Solar system (grey). We highlight two ranges of interest in mass and incident flux. All objects contained in the exoplanet.eu catalogue (Schneider et al. 2011), found by the RV, or the transit method, are included in this figure.
thumbnail Fig. 13

Densities, and the incident fluxes received by the TRAPPIST-1 planets (red), compared to the Solar System’s telluric planets and the Moon (grey). Other exoplanets with reported uncertainties on mass and radius smaller than 100% are shown in yellow and originate from the TEPCAT catalogue (Southworth 2011).

8 Conclusions

In this paper, we have used the most recent set of transit timings of the TRAPPIST-1 system to constrain the masses, densities of the seven Earth-size planets found earlier this year. Our purpose-built TTV code enables an extensive exploration of the parameter space combined to a full n-body integration scheme. Our results yield a significant improvement in our knowledge of the planetary bulk density, with corresponding uncertainties ranging between 5 and 12%. This level of precision is unprecedented for exoplanets receiving modest irradiation in this mass range. Our conclusions regarding the nature of the TRAPPIST-1 planets are the following:

  • The TRAPPIST-1 planets display densities ranging from 0.6 to 1.0 ρ⊕.

  • TRAPPIST-1 c and e likely have largely rocky interiors.

  • TRAPPIST-1 b, d, f, g and h require envelopes of volatiles in the form of thick atmospheres, oceans, or ice, in most cases with water mass fractions ≲5%. For comparison, the Earth’s water content is <0.1%.

  • TRAPPIST-1 d, f, g and h are unlikely to have an enriched atmosphere (e.g. CO2) above a bare core (assuming a standard Earth-like composition) without invoking unrealistically large quantities of gas.

  • TRAPPIST-1 b is the only planet above the runaway greenhouse limit that seems to require volatiles, with pressures of water vapour of the order of 101–104 bar.

These updated mass and density measurements represent key information for upcoming studies straddling astrophysics, planetary sciences, and geophysics aimed at an improved understanding of the interiors of temperate, Earth-sized planets.

Acknowledgements

We are grateful to the referee for an helpful review that improved the manuscript. We thank Robert Hurt for suggesting the inclusion of Fig. 13 as well as Yann Alibert, Gavin Coleman, Apurva Ozaand Christoph Mordasini for insightful discussions on the TRAPPIST-1 system. B.-O.D. acknowledges support from the Swiss National Science Foundation (PP00P2-163967). This work has been carried out within the frame ofthe National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 679030/WHIPLASH), and FP/2007–2013 grant agreement No. 336480, as well as from the ARC grant for Concerted Research Actions, financed by the Wallonia-Brussels Federation. M. Gillon, and V. Van Grootel are Belgian F.R.S.-FNRS Research Associates, E. Jehin is F.R.S.-FNRS Senior Research Associate. S.N.R. thanks the Agence Nationale pour la Recherche for support via grant ANR-13-BS05-0003-002 (grant MOJO). This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Support for this work was provided by NASA through an award issued by JPL/Caltech. E.A. acknowledges NSF grant AST-1615315, NASA grant NNX14AK26G and from the NASA Astrobiology Institute’s Virtual Planetary Laboratory Lead Team, funded through the NASA Astrobiology Institute under solicitation NNH12ZDA002C and CooperativeAgreement Number NNA13AA93A This paper includes data collected by the K2 mission. Funding for the K2 mission is provided by the NASA Science Mission directorate. Calculations were performed on UBELIX (http://www.id.unibe.ch/hpc), the HPC cluster at the University of Bern.

Appendix A Timings – planet b

Table A.1

Planet b.

Table A.1

continued.

Appendix B Timings – planet c

Table B.1

Planet c.

Table B.1

continued.

Appendix C: Timings – planet d

Table C.1

Planet d.

Appendix D Timings – planet e

Table D.1

Planet e

Appendix E Timings – planet f

Table E.1

Planet f.

Appendix F Timings – planet g

Table F.1

Planet g.

Appendix G: Timings – planet h

Table G.1

Planet h.

References

  1. Agol, E., & Fabrycky, D. 2017, Handbook of Exoplanets, eds. H. J., Deeg & J. A., Belmonte (Basel: Springer), 7 [Google Scholar]
  2. Agol, E., Steffen, J., Sari, R., & Clarkson, W. 2005, MNRAS, 359, 567 [NASA ADS] [CrossRef] [Google Scholar]
  3. Aigrain, S., Hodgkin, S. T., Irwin, M. J., Lewis, J. R., & Roberts, S. J. 2015, MNRAS, 447, 2880 [NASA ADS] [CrossRef] [Google Scholar]
  4. Aigrain, S., Parviainen, H., & Pope, B. J. S. 2016, MNRAS, 459, 2408 [NASA ADS] [Google Scholar]
  5. Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2016, IEEE Trans. Pattern. Anal., 38, 252 [CrossRef] [Google Scholar]
  6. Armitage, P. J., Eisner, J. A., & Simon, J. B. 2016, ApJ, 828, L2 [NASA ADS] [CrossRef] [Google Scholar]
  7. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  8. Baruteau, C., Crida, A., Paardekooper, S.-J., et al. 2014, Protostars and Planets VI, eds. H, Beuther, R S., Klessen, C. P., Dullemond, & T., Henning (Tucson: University of Arizona), 667 [Google Scholar]
  9. Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Bolmont, E., Raymond, S. N., Leconte, J., Hersant, F., & Correia, A. C. M. 2015, A&A, 583, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Bolmont, E., Selsis, F., Owen, J. E., et al. 2017, MNRAS, 464, 3728 [Google Scholar]
  12. Bouchet, J., Mazevet, S., Morard, G., Guyot, F., & Musella, R. 2013, Phys. Rev. B, 87, 094102 [NASA ADS] [CrossRef] [Google Scholar]
  13. Bourrier, V., Ehrenreich, D., Wheatley, P. J., et al. 2017a, A&A, 599, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Bourrier, V., de Wit, J., Bolmont, E., et al. 2017b, AJ, 154, 121 [NASA ADS] [CrossRef] [Google Scholar]
  15. Braak, C. J. F. T. 2006, Stat. Comput., 16, 239 [Google Scholar]
  16. Cabrera, J., Csizmadia, S., Lehmann, H., et al. 2014, ApJ, 781, 18 [NASA ADS] [CrossRef] [Google Scholar]
  17. Canup, R. M., & Ward, W. R. 2006, Nature, 441, 834 [NASA ADS] [CrossRef] [Google Scholar]
  18. Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Chambers, J. E. 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  20. Claret, A., & Bloemen, S. 2011, A&A, 529, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Connolly, J. 2009, Geochem. Geophys. Geosyst., 10 [Google Scholar]
  22. Cossou, C., Raymond, S. N., Hersant, F., & Pierens, A. 2014, A&A, 569, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Deck, K. M., & Agol, E. 2015, ApJ, 802, 116 [NASA ADS] [CrossRef] [Google Scholar]
  24. Deck, K. M., Agol, E., Holman, M. J., & Nesvorný, D. 2014, ApJ, 787, 132 [CrossRef] [Google Scholar]
  25. Delrez, L., Gillon, M., Triaud, A. H. M. J., et al. 2018, MNRAS, 475, 3577 [NASA ADS] [CrossRef] [Google Scholar]
  26. de Wit, J., Wakeford, H. R., Gillon, M., et al. 2016, Nature, 537, 69 [NASA ADS] [CrossRef] [Google Scholar]
  27. Dittmann, J. A., Irwin, J. M., Charbonneau, D., et al. 2017, Nature, 544, 333 [NASA ADS] [CrossRef] [Google Scholar]
  28. Dorn, C., Venturini, J., Khan, A., et al. 2017, A&A, 597, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Dra̧żkowska, J., & Dullemond, C. P. 2014, A&A, 572, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Fabrycky, D. C. 2010, Non-Keplerian Dynamics of Exoplanets, ed. S., Seager (Tuscon: University of Arizona Press), 217 [Google Scholar]
  31. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [CrossRef] [Google Scholar]
  32. Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 [NASA ADS] [CrossRef] [Google Scholar]
  33. Genda, H., & Abe, Y. 2005, Nature, 433, 842 [NASA ADS] [CrossRef] [Google Scholar]
  34. Gillon, M., Demory, B.-O., Benneke, B., et al. 2012, A&A, 539, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Gillon, M., Jehin, E., Lederer, S. M., et al. 2016, Nature, 533, 221 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  36. Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  37. Gimenez, A., & Garcia-Pelayo, J. M. 1983, Ap&SS, 92, 203 [NASA ADS] [CrossRef] [Google Scholar]
  38. Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65 [Google Scholar]
  39. Goździewski,K., Migaszewski, C., Panichi, F., & Szuszkiewicz, E. 2016, MNRAS, 455, L104 [Google Scholar]
  40. Grimm, R. E., & McSween, H. Y. 1993, Science, 259, 653 [NASA ADS] [Google Scholar]
  41. Grimm, S. L.,& Stadel, J. G. 2014, ApJ, 796, 23 [NASA ADS] [CrossRef] [Google Scholar]
  42. Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Hernandez, D. M., & Bertschinger, E. 2015, MNRAS, 452, 1934 [NASA ADS] [CrossRef] [Google Scholar]
  44. Holman, M. J. 2005, Science, 307, 1288 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  45. Holman, M. J., Fabrycky, D. C., Ragozzine, D., et al. 2010, Science, 330, 51 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  46. Howard, A. W., Marcy, G. W., Bryson, S. T., et al. 2012, ApJS, 201, 15 [NASA ADS] [CrossRef] [Google Scholar]
  47. Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398 [NASA ADS] [CrossRef] [Google Scholar]
  48. Izidoro, A., Morbidelli, A., & Raymond, S. N. 2014, ApJ, 794, 11 [NASA ADS] [CrossRef] [Google Scholar]
  49. Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, MNRAS, 470, 1750 [NASA ADS] [CrossRef] [Google Scholar]
  50. Johansen, A., Blum, J., Tanaka, H., et al. 2014, Protostars and Planets VI, eds. H, Beuther, R S., Klessen, C. P., Dullemond, & T., Henning (Tucson: University of Arizona), 547 [Google Scholar]
  51. Jontof-Hutter,D., Lissauer, J. J., Rowe, J. F., & Fabrycky, D. C. 2014, ApJ, 785, 15 [NASA ADS] [CrossRef] [Google Scholar]
  52. Jontof-Hutter,D., Rowe, J. F., Lissauer, J. J., Fabrycky, D. C., & Ford, E. B. 2015, Nature, 522, 321 [NASA ADS] [CrossRef] [Google Scholar]
  53. Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131 [NASA ADS] [CrossRef] [Google Scholar]
  54. Kopparapu, R. k., Wolf, E. T., Haqq-Misra, J., et al. 2016, ApJ, 819, 84 [Google Scholar]
  55. Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [Google Scholar]
  56. Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Lissauer, J. J. 2007, ApJ 660, L149 [NASA ADS] [CrossRef] [Google Scholar]
  58. Lithwick, Y., Xie, J., & Wu, Y. 2012, ApJ, 761, 122 [NASA ADS] [CrossRef] [Google Scholar]
  59. Lodders, K. 2003, ApJ, 591, 1220 [NASA ADS] [CrossRef] [Google Scholar]
  60. Luger, R., Lustig-Yaeger, J., & Agol, E. 2017a, ApJ, 851, 94 [NASA ADS] [CrossRef] [Google Scholar]
  61. Luger, R., Sestovic, M., Kruse, E., et al. 2017b, New Astron., 1, 0129 [Google Scholar]
  62. Mandel, K., & Agol, E. 2002, ApJ, 580, L171 [NASA ADS] [CrossRef] [Google Scholar]
  63. Matsumoto, Y., Nagasawa, M., & Ida, S. 2012, Icarus, 221, 624 [NASA ADS] [CrossRef] [Google Scholar]
  64. Mayor, M., Marmier, M., Lovis, C., et al. 2011, ArXiv e-prints [arXiv:1109.2497] [Google Scholar]
  65. Mills, S. M., Fabrycky, D. C., Migaszewski, C., et al. 2016, Nature, 533, 509 [NASA ADS] [CrossRef] [Google Scholar]
  66. Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418 [NASA ADS] [CrossRef] [Google Scholar]
  67. Neron de Surgy, O., & Laskar, J. 1997, A&A, 318, 975 [NASA ADS] [Google Scholar]
  68. Nesvorný, D., & Vokrouhlický, D. 2014, ApJ, 790, 58 [NASA ADS] [CrossRef] [Google Scholar]
  69. Ogihara, M., & Ida, S. 2009, ApJ, 699, 824 [NASA ADS] [CrossRef] [Google Scholar]
  70. Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Ormel, C. W., Liu, B., & Schoonenberg, D. 2017, A&A, 604, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Quarles, B., Quintana, E. V., Lopez, E., Schlieder, J. E., & Barclay, T. 2017, ApJ, 842, L5 [NASA ADS] [CrossRef] [Google Scholar]
  73. Raymond, S. N., Mandell, A. M., & Sigurdsson, S. 2006, Science, 313, 1413 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  74. Raymond, S. N., Scalo, J., & Meadows, V. S. 2007, ApJ, 669, 606 [Google Scholar]
  75. Raymond, S. N., Barnes, R., & Mandell, A. M. 2008, MNRAS, 384, 663 [NASA ADS] [CrossRef] [Google Scholar]
  76. Rogers, L. A. 2015, ApJ, 801, 41 [NASA ADS] [CrossRef] [Google Scholar]
  77. Schneider, J., Dedieu, C., Le Sidaner, P., Savalle, R., & Zolotukhin, I. 2011, A&A, 532, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Schoonenberg,D., & Ormel, C. W. 2017, A&A, 602, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Seager, S., & Mallén-Ornelas, G. 2003, ApJ, 585, 1038 [NASA ADS] [CrossRef] [Google Scholar]
  80. Seager, S., Kuchner, M., Hier-Majumder, C., & Militzer, B. 2007, ApJ, 669, 1279 [NASA ADS] [CrossRef] [Google Scholar]
  81. Southworth, J. 2011, MNRAS, 417, 2166 [NASA ADS] [CrossRef] [Google Scholar]
  82. Spencer, J. R., Rathbun, J. A., Travis, L. D., et al. 2000, Science, 288, 1198 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  83. Storn, R., & Price, K. 1997, J. Global Optim., 11, 341 [Google Scholar]
  84. Tamayo, D., Rein, H., Petrovich, C., & Murray, N. 2017, ApJ, 840, L19 [NASA ADS] [CrossRef] [Google Scholar]
  85. Terquem, C., & Papaloizou, J. C. B. 2007, ApJ, 654, 1110 [NASA ADS] [CrossRef] [Google Scholar]
  86. Turbet, M., Bolmont, E., Leconte, J., et al. 2018, A&A, 612, A86 [CrossRef] [EDP Sciences] [Google Scholar]
  87. Unterborn, C. T., Desch, S. J., Hinkel, N., & Lorenzo, A. 2018, Nat. Astron., 2, 297 [NASA ADS] [CrossRef] [Google Scholar]
  88. Van Grootel, V., Fernandes, C. S., Gillon, M., et al. 2018, ApJ, 853, 30 [NASA ADS] [CrossRef] [Google Scholar]
  89. Vazan, A., Kovetz, A., Podolak, M., & Helled, R. 2013, MNRAS, 434, 3283 [NASA ADS] [CrossRef] [Google Scholar]
  90. Vrugt, J. A., Ter Braak, C., Diks, C., et al. 2009, Int. J. Nonlinear Sci. Numer. Simul., 10, 273 [Google Scholar]
  91. Waldmann, I. P., Tinetti, G., Rocchetto, M., et al. 2015, ApJ, 802, 107 [NASA ADS] [CrossRef] [Google Scholar]
  92. Wang, S., Wu,D.-H., Barclay, T., & Laughlin, G. P. 2017, ArXiv e-prints [arXiv:1704.04290] [Google Scholar]
  93. Ward, W. R. 1997, Icarus, 126, 261 [NASA ADS] [CrossRef] [Google Scholar]
  94. Wheatley, P. J., Louden, T., Bourrier, V., et al. 2017, MNRAS, 465, L74 [NASA ADS] [CrossRef] [Google Scholar]
  95. Wisdom, J., & Hernandez, D. M. 2015, MNRAS, 453, 3015 [NASA ADS] [CrossRef] [Google Scholar]
  96. Wordsworth, R. D., Forget, F., Selsis, F., et al. 2010, A&A, 522, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
1

This equation is only valid for i = 90°. A discussion for the case i≠90° is given in Gimenez & Garcia-Pelayo (1983).

2

Tests show that setting higher initial values of e, up to 0.2, does not change the results. Additionally, having higher eccentricities make such a packed planetary system very likely to be dynamically unstable. In the long term stability analysis in Sect. 6.2, the eccentricities remain below 0.025.

All Tables

Table 1

Number of observations, number of degrees of freedom Ndof = Nobservations − Nparameters, χ2 from Eq. (8) and reduced χ2 = χ2∕Ndof, for all planets separately and for all planets together.

In the text
Table 2

Updated masses, radii (Delrez et al. 2018), correlation coefficients cm,r of masses and radii as well as the densities and the surface gravity of all seven planets.

In the text
Table 3

Median and standard deviation σ of the semi-major axis a, eccentricity e, argument of periapsis ω, and mean anomaly M of the seven planets.

In the text
Table 4

Planetary characteristics of the TRAPPIST-1 planets derived from our 1D numerical climate simulations.

In the text
Table A.1

Planet b.

In the text
Table A.1

continued.

In the text
Table B.1

Planet c.

In the text
Table B.1

continued.

In the text
Table C.1

Planet d.

In the text
Table D.1

Planet e

In the text
Table E.1

Planet f.

In the text
Table F.1

Planet g.

In the text
Table G.1

Planet h.

In the text

All Figures

thumbnail Fig. 1

Folded short-cadence lightcurves extracted from K2 data, corrected for their corrected for their TTVs. Relative fluxes are shown for each planet and are shifted vertically for clarity. The short-cadence data are binned to produce the coloured points, with five cadences per point. The white points are further binned, with ten points taken from the folded curve per bin.
In the text
thumbnail Fig. 2

Calculated TTVs for all planets and for 1000 different MCMC samples (grey lines). Measured transit times with the corresponding uncertainties are indicated by coloured symbols, according to the used telescope. A detailed list of all transits is given in Appendix. The differences between the solutions reflect the distribution shown in Fig. 4 for τ = 1. The two panels at the bottom show a zoomed region for planet b and c.
In the text
thumbnail Fig. 3

Residuals of the TTVs shown in Fig. 2. The transit data index corresponds to the column index of Tables A.1–G.1.
In the text
thumbnail Fig. 4

Posterior probability distribution of the mass and eccentricities of all seven planets, for different MCMC sampling “temperatures” τ, assuming a stellar mass of M⋆ = 0.09 M⊙ (Van Grootel et al. 2018). The contours correspond to significance levels of 68%, 95% and 99% for τ = 1. Sampling “temperatures” greater than one are used to cover a larger parameter space.
In the text
thumbnail Fig. 5

Left panels: evolution of the masses of all planets of the entire ensemble of 100 walker chains in blue, and the evolution of a single chain in green. Right panels: auto-correlation function of the masses for the ensemble and a single chain.
In the text
thumbnail Fig. 6

Posterior probability distribution between all pairs of parameters, the masses m, the semi-major axes a, the eccentricities e, the arguments of perihelion ω, and the mean longitudes l = Ω + ω + M of all seven planets. It is showing strong correlations between many pairs of parameters, especially between m and a and between e and ω.
In the text
thumbnail Fig. 7

Temporal evolution of the eccentricities of all planets. The eccentricities are well constrained to values below 0.015 and show a regular behaviour for 2 Myr. After that, the systems show irregularities, caused by small uncertainties in the initial conditions. We also show the relative energy error and the relative angular momentum error of the integrator.
In the text
thumbnail Fig. 8

Temporal evolution of ω, showing a fast precession and strong mutual orbital perturbations between the planets.
In the text
thumbnail Fig. 9

Temporal evolution of the three body resonant angle ϕi = pλi−1 − (p + q)λi + qλi+1, where the values of p and q for consecutive triples of planets are (2,3), (1,2), 2,3), (1,2) and (1,1) (Luger et al. 2017b). After 2 Myr, the systems show irregularities, caused by small uncertainties in the initial conditions.
In the text
thumbnail Fig. 10

Mass–radius diagram for the TRAPPIST-1 planets, Earth, and Venus. Curves trace idealised compositions of rocky and water-rich interiors (surface temperature fixed to 200 K). Median values are highlighted by a black dot. Coloured contours correspond to significance levels of 68% and 95% for each planet. The interiors are calculated with model II of Dorn et al. (2017). Rocky interiors are composed of Fe, Si, Mg, and O, assuming different bulk ratios of Fe/Mg and Mg/Si. U17 refers to Unterborn et al. (2018). Equilibrium temperatures for each planet are indicated by the coloured contours.
In the text
thumbnail Fig. 11

Ratio between the mass of a planet or moon and that of its host as a function of orbital period. We represent the known population of exoplanets (from exoplanet.eu, Schneider et al. 2011) from three different detection methods. The Solar system is also highlighted. The TRAPPIST-1 system, like the Galilean moons of Jupiter, share a similar parameter space with the sub-Neptune and super-Earth population. Orbital periods are used to reveal the various sub-populations.
In the text
thumbnail Fig. 12

Masses and the incident fluxes received by the TRAPPIST-1 planets, compared to other exoplanets found by the transit method (yellow), via the radial-velocity technique (blue), and to the terrestrial worlds of the Solar system (grey). We highlight two ranges of interest in mass and incident flux. All objects contained in the exoplanet.eu catalogue (Schneider et al. 2011), found by the RV, or the transit method, are included in this figure.
In the text
thumbnail Fig. 13

Densities, and the incident fluxes received by the TRAPPIST-1 planets (red), compared to the Solar System’s telluric planets and the Moon (grey). Other exoplanets with reported uncertainties on mass and radius smaller than 100% are shown in yellow and originate from the TEPCAT catalogue (Southworth 2011).
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.