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/00046361/201732233  
Published online  04 June 2018 
The nature of the TRAPPIST1 exoplanets
^{1}
University of Bern, Center for Space and Habitability,
Gesellschaftsstrasse 6,
3012
Bern, Switzerland
email: brice.demory@csh.unibe.ch
^{2}
Space Sciences, Technologies and Astrophysics Research (STAR) Institute, Université de Liège,
Allée du 6 août 19C,
4000
Liège, Belgium
^{3}
Astronomy Department, University of Washington,
Seattle,
WA
98195, USA
^{4}
NASA Astrobiology Institute’s Virtual Planetary Laboratory,
Seattle,
WA
98195, USA
^{5}
Cavendish Laboratory,
J J Thomson Avenue,
Cambridge
CB3 0HE, UK
^{6}
Institute of Astronomy,
Madingley Road,
Cambridge
CB3 0HA, UK
^{7}
School of Physics & Astronomy, University of Birmingham,
Edgbaston,
Birmingham
B15 2TT, UK
^{8}
Laboratoire de Météorologie Dynamique, IPSL, Sorbonne Universités, UPMC Univ. Paris 06, CNRS,
4 place Jussieu,
75005
Paris, France
^{9}
Université Paris Diderot, AIM, Sorbonne Paris Cité, CEA, CNRS,
91191
GifsurYvette, France
^{10}
Laboratoire d’astrophysique de Bordeaux, Univ. Bordeaux, CNRS,
B18N, Allée Geoffroy SaintHilaire,
33615
Pessac, France
^{11}
Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology,
77 Massachusetts Avenue,
Cambridge,
MA
02139, USA
^{12}
Center for Astrophysics and Space Science, University of California San Diego,
La Jolla,
CA
92093, USA
^{13}
IPAC, Mail Code 3146, Calif. Inst. of Technology,
1200 E California Blvd,
Pasadena,
CA
91125, USA
^{14}
Department of Astronomy and Astrophysics, Univ. of Chicago,
5640 S Ellis Ave,
Chicago,
IL
60637, USA
^{15}
Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology,
77 Massachusetts Ave.,
Cambridge,
MA
02139, USA
^{16}
NASA Johnson Space Center,
2101 NASA Parkway,
Houston,
TX
77058, USA
^{17}
Guggenheim Fellow
^{18}
Institut d’Astrophysique de Paris,
98 bis Boulevard Arago,
75014
Paris, France
^{19}
University of Zurich, Institute of Computational Sciences,
Winterthurerstrasse 190,
8057
Zurich, Switzerland
Received:
3
November
2017
Accepted:
21
January
2018
Context. The TRAPPIST1 system hosts seven Earthsized, temperate exoplanets orbiting an ultracool dwarf star. As such, it represents a remarkable setting to study the formation and evolution of terrestrial planets that formed in the same protoplanetary disk. While the sizes of the TRAPPIST1 planets are all known to better than 5% precision, their densities have significant uncertainties (between 28% and 95%) because of poor constraints on the planet’s masses.
Aims. The goal of this paper is to improve our knowledge of the TRAPPIST1 planetary masses and densities using transittiming variations (TTVs). The complexity of the TTV inversion problem is known to be particularly acute in multiplanetary systems (convergence issues, degeneracies and size of the parameter space), especially for resonant chain systems such as TRAPPIST1.
Methods. To overcome these challenges, we have used a novel method that employs a genetic algorithm coupled to a full Nbody integrator that we applied to a set of 284 individual transit timings. This approach enables us to efficiently explore the parameter space and to derive reliable masses and densities from TTVs for all seven planets.
Results. Our new masses result in a five to eightfold improvement on the planetary density uncertainties, with precisions ranging from 5% to 12%. These updated values provide new insights into the bulk structure of the TRAPPIST1 planets. We find that TRAPPIST1 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 less than 5%.
Key words: methods: numerical / planets and satellites: detection / planets and satellites: individual: TRAPPIST1
© ESO 2018
1 Introduction
The TRAPPIST1 system, which harbours seven Earthsize exoplanets orbiting an ultracool star (Gillon et al. 2017), represents a fascinating setting to study the formation and evolution of tightly packed small planet systems. While the TRAPPIST1 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 longterm stability of the system. So far, most exoplanet masses have been measured using the radialvelocity technique. But because of the TRAPPIST1 faintness (V = 19), precise constraints on Earthmass 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 TRAPPIST1 system shows dynamical similarities to Kepler90 system (Cabrera et al. 2014) which contains also seven planets and resonances conditions between pairs of them.
Planetary masses published in the TRAPPIST1 discovery paper (Gillon et al. 2017) were bearing conservative uncertainties because the different techniques used by the authors suggested a nonmonotonous parameter space with the absence of a single global minimum. Subsequent studies have adequately invoked the requirement for longterm 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 TRAPPIST1 (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 loworder symplectic coordinates 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 Nbody 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 transittiming 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 TRAPPIST1 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 shortcadence photometry
For the purpose of this analysis we have included the transit timings derived from the K2 photometry (Howell et al. 2014), which observed TRAPPIST1 during Campaign 12 (Luger et al. 2017b). We detail in the following the data reduction of this dataset. We used the K2’s pipelinecalibrated shortcadence target pixel files (TPF) that includes the correct timestamps. The K2 TPF TRAPPIST1 (EPIC ID 246199087) aperture is a 9 × 10 postage stamp centred on the target star, with 1minute 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 tophat 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 Gaussianprocesses (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 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}, σ) 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 maximumlikelihood 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 intransit points and flares). Out of the entire dataset, we discarded only two transits because of a low signaltonoise ratio at the following BJD_{TDB}: 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 TRAPPIST1b, 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 midtime centre and impact parameter for all planets to the published values (Gillon et al. 2017). We computed the quadratic limbdarkening coefficients u_{1} and u_{2} 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 TRAPPIST1 planets. We report the median and 1sigma credible intervals of the posterior distribution functions for the 124 K2 transit timings in Tables A.1–G.1. The resulting K2 shortcadence stacked lightcurves are shown in Fig. 1.
Fig. 1 Folded shortcadence 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 shortcadence 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 TRAPPIST1 system are challenging due to the 7planet, Laplace resonancechain 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 Nbody 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 multidimensional 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, semimajor 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 coordinate. The midtransit 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)
The quantities x_{i} and y_{i} are the astrocentric coordinates of the planet i.
We useda prechecker in the integrator to determine if a transit is expected to happen during the next time step. This will happen if the value of g_{i} moves from a negative value to a positive value during the next time step, and if z_{i} > 0. In addition, we used the conditions g∕ġ < 3.5 ⋅ dt and r_{sky} < (R_{star} + R_{i}) + v_{i}⋅dt to refine the prechecker, where r_{sky} = (x_{i}, y_{i}) is the radial coordinate on the sky plane, v_{i} is the norm of the velocity, and R_{i} the radius ofplanet i; and R_{⋆} the radius of the star.
Since the integration is performed with a symplectic integrator, the coordinates 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 midtransit time. But when the time step is small enough, this error can be safely neglected. Also, in highly eccentric orbits, the described prechecker 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 fourthorder integrator scheme with a time step of 0.08 days. When the prechecker has detected a transit candidate, then all planets are integrated with a Bulirsh–Stoer direct Nbody method for a time step and the Eq. (5) is iterated until δt is smaller than a tolerance value. A transit is reported if r_{sky} < R_{star} + R_{i}.
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 x_{i}. 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, π(x_{p})∕π(x_{i})). When the proposal is accepted, then x_{i} is replaced by x_{p} in the next generation; otherwise the state remains unchanged. In each 30^{th} generation, we set γ = 0.98 to allow jumps between multimodal solutions (Braak 2006). In addition, we set γ = 0.01 and x _{i} = x_{l≠i} during the burnin phase to eliminate outliers. Alternatively, we also tested the affine invariant ensemble walker MCMC method (Goodman & Weare 2010; ForemanMackey et al. 2013), which yields a comparable performance.
For the probability density function, π(x_{i}) we used (8)
where the running index t refers to the transit epoch, T_{calc} are the calculated mid transit times, T_{obs} are the observed midtransit 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 soobtained 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 P_{i}, the time of the first transit T_{i}, the start time of the simulation t, the mean anomaly at the first transit , the mean anomaly M_{i}, the eccentricity e_{i}, 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 semimajor axis. We used the square root of the eccentricity in the parameters x_{i} and y_{i} 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 threedimensional 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 m_{i} ∈[0, 6 × 10^{−6}M_{⊙}], e_{i} ∈ [0, 0.05]^{2}, ω_{i} ∈ [0, 2π] and M_{i} ∈ [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 P_{i} 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 semimajor axes and mean anomalies.
3.3.1 Substep 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 nonlinear 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 semimajor 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 substep optimisation scheme to find the optimal values for the semimajor axis a and the mean anomaly M. The substep 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 x_{1}, x_{2}, and x_{3} as follows:
where x means either a or M, and y_{j} means the values of χ^{2} at locations x_{j}. The cost of the described substep optimisation scheme is, that three times as many walkers are needed to generate the values y_{1}, y_{2}, and y_{3}, and each DEMCMC step has to be followed by 14 substeps 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 substep sampling for 100 chains are shown in Fig. 5, which shows efficient convergence of the chains.
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 planeparallel geometry, and alternates between drifts and universal Kepler steps between pairs of planets. A drift in Cartesian coordinatesis defined as an update of some or all positions assuming constant velocities. The initial conditions are constructed with Jacobi coordinates, and the integration uses Cartesian coordinates in a centerofmass frame. The transit times are found by tracking the projected sky position and velocity, and finding when the dot product changes sign, using Eqs. 46. 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 massratio posteriors and , along with the well constrained stellar density ρ_{⋆} determined photometrically from the Spitzer dataset (Seager & MallénOrnelas 2003). The planet density then is (JontofHutter 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 reanalysed our initial Spitzer and groundbased 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.
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. 
Updated masses, radii (Delrez et al. 2018), correlation coefficients c_{m,r} of masses and radii as well as the densities and the surface gravity of all seven planets.
Median and standard deviation σ of the semimajor axis a, eccentricity e, argument of periapsis ω, and mean anomaly M of the seven planets.
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 TRAPPIST1 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 timescales.
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 TRAPPIST1 planets. A combination of the 1) highprecision Spitzer photometry, 2) K2’s 80day long quasicontinuous coverage and the 3) resolved highfrequency component of the TTV pattern known as “chopping” (Holman et al. 2010, Deck & Agol 2015 all contribute to mitigate the masseccentricity 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 Longterm 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 threebody 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 verylong 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 MercuryT 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.
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 TRAPPIST1 planets
Our improved masses and densities show that TRAPPIST1 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 silicatemantle 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 TRAPPIST1 planets with theoretical mass–radiusrelations calculated with published models (Dorn et al. 2017). We estimate the probability p_{volatiles} of each planet to be volatilerich 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 volatilerich layers (with p_{volatiles} of at least 0.96, 0.99, 0.66, 1, and 0.71, respectively). Volatilerich layers could comprise atmosphere, oceans, and/or ice layers. In contrast, planets c and e may be rocky (with p_{volatiles} of at least 0.24 and <0.01, respectively).
The comparison of masses and radii with the idealised interior endmember of m_{water}∕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.
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: autocorrelation function of the masses for the ensemble and a single chain. 
7.2 Limits for the possible atmospherescenarios
We use the LMDG onedimensional cloudfree numerical climate model (Wordsworth et al. 2010)to simulate the vertical temperature profiles of the seven TRAPPIST1 planets. Calculations are performed using a synthetic spectrum of TRAPPIST1 and molecular spectroscopic properties from Turbet et al. (2018). We assume atmospheric compositions that range from pure H_{2}, H_{2} –CH_{4}, H_{2} –H_{2}O to pure CO_{2}. 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 10^{3} bar, see Table 4), we decompose the thermal structure of the atmosphere into 500 logspaced layers in altitude coordinates. 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 P_{surf} can be estimated by integrating the hydrostatic equation, which yields: (16)
where P_{transit} is the pressure at the transit radius R_{transit}, R_{core} 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 m_{H} 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 Earthlike composition) without unrealistically large quantities of gas.
For the colder, low density planets (f, g, and h), explaining the radius with only a CO_{2} atmosphere is difficult due to the small pressure scale height and the fact that CO_{2} 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 10^{1} –10^{4} bar, which could explain its relatively low density (assuming P_{transit} = 20 mbar). As such, TRAPPIST1 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 H_{2} dominated atmospheresthicker than a few bars. Within these assumptions and considering the expected intense atmospheric escape around TRAPPIST1 (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.
Fig. 6 Posterior probability distribution between all pairs of parameters, the masses m, the semimajor 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 Massratios, densities and irradiation
The mass ofa celestial object is its most fundamental property. We now compare the masses of the TRAPPIST1 planets and place them into wider context. Exoplanet discoveries are heavily biased towards single Sunlike stars. TRAPPIST1 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 massratio versus period diagram comparing the TRAPPIST1 system to other exoplanets (where the orbital period is used to separate various subpopulations). We also push the comparison to include the planets of the Solar system, and the principal moons of Jupiter. TRAPPIST1’s planets cover massratio a range 10^{−4}–10^{−5}, which is shared with the subNeptune and superEarth exoplanets population that orbits Sunlike stars. The similarity may suggest a similar formation mechanism, or at minimum a comparable scaling in protoplanetary disk mass. Notably, subNeptunes and superEarths are the most abundant planet types for Sunlike 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 TRAPPIST1 planets plays an important role in their evolution. It is thus also insightful to compare the TRAPPIST1 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 TRAPPIST1d & 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) & Kepler138b (0.64, 2.33) (JontofHutter 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 TRAPPIST1 planets with their irradiation level in Fig. 13. We note that the density trend vs. irradiation increases until TRAPPIST1e, where it peaks and then decreases for the outer planets. One object stands out of that pattern, TRAPPIST1d, which interestingly, has a similar bulk density as the Moon. These aspects will be particularly relevant to better understand the formation pathway of the TRAPPIST1 system.
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. 
Fig. 8 Temporal evolution of ω, showing a fast precession and strong mutual orbital perturbations between the planets. 
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 shortlived radionuclides (Grimm & McSween 1993), water loss during giant impacts between embryos (Genda & Abe 2005), and heating during the very rapid growth expected around lowmass 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 ~Earthmass 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 solidstogas 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 inwardsdrifting 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 selflimiting. 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 TRAPPIST1 (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 TRAPPIST1 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 diskdriven 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 TRAPPIST1 can be interpreted as a system that underwent a relatively slow migration creating a longlived resonant system.
Fig. 10 Mass–radius diagram for the TRAPPIST1 planets, Earth, and Venus. Curves trace idealised compositions of rocky and waterrich 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. 
Planetary characteristics of the TRAPPIST1 planets derived from our 1D numerical climate simulations.
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 TRAPPIST1 system, like the Galilean moons of Jupiter, share a similar parameter space with the subNeptune and superEarth population. Orbital periods are used to reveal the various subpopulations. 
Fig. 12 Masses and the incident fluxes received by the TRAPPIST1 planets, compared to other exoplanets found by the transit method (yellow), via the radialvelocity 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. 
Fig. 13 Densities, and the incident fluxes received by the TRAPPIST1 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 TRAPPIST1 system to constrain the masses, densities of the seven Earthsize planets found earlier this year. Our purposebuilt TTV code enables an extensive exploration of the parameter space combined to a full nbody 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 TRAPPIST1 planets are the following:

The TRAPPIST1 planets display densities ranging from 0.6 to 1.0 ρ_{⊕}.

TRAPPIST1 c and e likely have largely rocky interiors.

TRAPPIST1 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%.

TRAPPIST1 d, f, g and h are unlikely to have an enriched atmosphere (e.g. CO_{2}) above a bare core (assuming a standard Earthlike composition) without invoking unrealistically large quantities of gas.

TRAPPIST1 b is the only planet above the runaway greenhouse limit that seems to require volatiles, with pressures of water vapour of the order of 10^{1}–10^{4} 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, Earthsized 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 TRAPPIST1 system. B.O.D. acknowledges support from the Swiss National Science Foundation (PP00P2163967). 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 WalloniaBrussels 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 ANR13BS050003002 (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 AST1615315, 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
Planet b.
continued.
Appendix B Timings – planet c
Planet c.
continued.
Appendix C: Timings – planet d
Planet d.
Appendix D Timings – planet e
Planet e
Appendix E Timings – planet f
Planet f.
Appendix F Timings – planet g
Planet g.
Appendix G: Timings – planet h
Planet h.
References
 Agol, E., & Fabrycky, D. 2017, Handbook of Exoplanets, eds. H. J., Deeg & J. A., Belmonte (Basel: Springer), 7 [Google Scholar]
 Agol, E., Steffen, J., Sari, R., & Clarkson, W. 2005, MNRAS, 359, 567 [NASA ADS] [CrossRef] [Google Scholar]
 Aigrain, S., Hodgkin, S. T., Irwin, M. J., Lewis, J. R., & Roberts, S. J. 2015, MNRAS, 447, 2880 [NASA ADS] [CrossRef] [Google Scholar]
 Aigrain, S., Parviainen, H., & Pope, B. J. S. 2016, MNRAS, 459, 2408 [NASA ADS] [Google Scholar]
 Ambikasaran, S., ForemanMackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2016, IEEE Trans. Pattern. Anal., 38, 252 [CrossRef] [Google Scholar]
 Armitage, P. J., Eisner, J. A., & Simon, J. B. 2016, ApJ, 828, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 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]
 Bolmont, E., Selsis, F., Owen, J. E., et al. 2017, MNRAS, 464, 3728 [Google Scholar]
 Bouchet, J., Mazevet, S., Morard, G., Guyot, F., & Musella, R. 2013, Phys. Rev. B, 87, 094102 [NASA ADS] [CrossRef] [Google Scholar]
 Bourrier, V., Ehrenreich, D., Wheatley, P. J., et al. 2017a, A&A, 599, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bourrier, V., de Wit, J., Bolmont, E., et al. 2017b, AJ, 154, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Braak, C. J. F. T. 2006, Stat. Comput., 16, 239 [Google Scholar]
 Cabrera, J., Csizmadia, S., Lehmann, H., et al. 2014, ApJ, 781, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Canup, R. M., & Ward, W. R. 2006, Nature, 441, 834 [NASA ADS] [CrossRef] [Google Scholar]
 Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chambers, J. E. 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Claret, A., & Bloemen, S. 2011, A&A, 529, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Connolly, J. 2009, Geochem. Geophys. Geosyst., 10 [Google Scholar]
 Cossou, C., Raymond, S. N., Hersant, F., & Pierens, A. 2014, A&A, 569, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deck, K. M., & Agol, E. 2015, ApJ, 802, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Deck, K. M., Agol, E., Holman, M. J., & Nesvorný, D. 2014, ApJ, 787, 132 [CrossRef] [Google Scholar]
 Delrez, L., Gillon, M., Triaud, A. H. M. J., et al. 2018, MNRAS, 475, 3577 [NASA ADS] [CrossRef] [Google Scholar]
 de Wit, J., Wakeford, H. R., Gillon, M., et al. 2016, Nature, 537, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Dittmann, J. A., Irwin, J. M., Charbonneau, D., et al. 2017, Nature, 544, 333 [NASA ADS] [CrossRef] [Google Scholar]
 Dorn, C., Venturini, J., Khan, A., et al. 2017, A&A, 597, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dra̧żkowska, J., & Dullemond, C. P. 2014, A&A, 572, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fabrycky, D. C. 2010, NonKeplerian Dynamics of Exoplanets, ed. S., Seager (Tuscon: University of Arizona Press), 217 [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [CrossRef] [Google Scholar]
 Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Genda, H., & Abe, Y. 2005, Nature, 433, 842 [NASA ADS] [CrossRef] [Google Scholar]
 Gillon, M., Demory, B.O., Benneke, B., et al. 2012, A&A, 539, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gillon, M., Jehin, E., Lederer, S. M., et al. 2016, Nature, 533, 221 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Gillon, M., Triaud, A. H. M. J., Demory, B.O., et al. 2017, Nature, 542, 456 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Gimenez, A., & GarciaPelayo, J. M. 1983, Ap&SS, 92, 203 [NASA ADS] [CrossRef] [Google Scholar]
 Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65 [Google Scholar]
 Goździewski,K., Migaszewski, C., Panichi, F., & Szuszkiewicz, E. 2016, MNRAS, 455, L104 [Google Scholar]
 Grimm, R. E., & McSween, H. Y. 1993, Science, 259, 653 [NASA ADS] [Google Scholar]
 Grimm, S. L.,& Stadel, J. G. 2014, ApJ, 796, 23 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 Hernandez, D. M., & Bertschinger, E. 2015, MNRAS, 452, 1934 [NASA ADS] [CrossRef] [Google Scholar]
 Holman, M. J. 2005, Science, 307, 1288 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Holman, M. J., Fabrycky, D. C., Ragozzine, D., et al. 2010, Science, 330, 51 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Howard, A. W., Marcy, G. W., Bryson, S. T., et al. 2012, ApJS, 201, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398 [NASA ADS] [CrossRef] [Google Scholar]
 Izidoro, A., Morbidelli, A., & Raymond, S. N. 2014, ApJ, 794, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, MNRAS, 470, 1750 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 JontofHutter,D., Lissauer, J. J., Rowe, J. F., & Fabrycky, D. C. 2014, ApJ, 785, 15 [NASA ADS] [CrossRef] [Google Scholar]
 JontofHutter,D., Rowe, J. F., Lissauer, J. J., Fabrycky, D. C., & Ford, E. B. 2015, Nature, 522, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Kopparapu, R. k., Wolf, E. T., HaqqMisra, J., et al. 2016, ApJ, 819, 84 [Google Scholar]
 Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [Google Scholar]
 Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lissauer, J. J. 2007, ApJ 660, L149 [NASA ADS] [CrossRef] [Google Scholar]
 Lithwick, Y., Xie, J., & Wu, Y. 2012, ApJ, 761, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Lodders, K. 2003, ApJ, 591, 1220 [NASA ADS] [CrossRef] [Google Scholar]
 Luger, R., LustigYaeger, J., & Agol, E. 2017a, ApJ, 851, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Luger, R., Sestovic, M., Kruse, E., et al. 2017b, New Astron., 1, 0129 [Google Scholar]
 Mandel, K., & Agol, E. 2002, ApJ, 580, L171 [NASA ADS] [CrossRef] [Google Scholar]
 Matsumoto, Y., Nagasawa, M., & Ida, S. 2012, Icarus, 221, 624 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., Marmier, M., Lovis, C., et al. 2011, ArXiv eprints [arXiv:1109.2497] [Google Scholar]
 Mills, S. M., Fabrycky, D. C., Migaszewski, C., et al. 2016, Nature, 533, 509 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418 [NASA ADS] [CrossRef] [Google Scholar]
 Neron de Surgy, O., & Laskar, J. 1997, A&A, 318, 975 [NASA ADS] [Google Scholar]
 Nesvorný, D., & Vokrouhlický, D. 2014, ApJ, 790, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Ogihara, M., & Ida, S. 2009, ApJ, 699, 824 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., Liu, B., & Schoonenberg, D. 2017, A&A, 604, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Quarles, B., Quintana, E. V., Lopez, E., Schlieder, J. E., & Barclay, T. 2017, ApJ, 842, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Raymond, S. N., Mandell, A. M., & Sigurdsson, S. 2006, Science, 313, 1413 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Raymond, S. N., Scalo, J., & Meadows, V. S. 2007, ApJ, 669, 606 [Google Scholar]
 Raymond, S. N., Barnes, R., & Mandell, A. M. 2008, MNRAS, 384, 663 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, L. A. 2015, ApJ, 801, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, J., Dedieu, C., Le Sidaner, P., Savalle, R., & Zolotukhin, I. 2011, A&A, 532, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schoonenberg,D., & Ormel, C. W. 2017, A&A, 602, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seager, S., & MallénOrnelas, G. 2003, ApJ, 585, 1038 [NASA ADS] [CrossRef] [Google Scholar]
 Seager, S., Kuchner, M., HierMajumder, C., & Militzer, B. 2007, ApJ, 669, 1279 [NASA ADS] [CrossRef] [Google Scholar]
 Southworth, J. 2011, MNRAS, 417, 2166 [NASA ADS] [CrossRef] [Google Scholar]
 Spencer, J. R., Rathbun, J. A., Travis, L. D., et al. 2000, Science, 288, 1198 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Storn, R., & Price, K. 1997, J. Global Optim., 11, 341 [Google Scholar]
 Tamayo, D., Rein, H., Petrovich, C., & Murray, N. 2017, ApJ, 840, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Terquem, C., & Papaloizou, J. C. B. 2007, ApJ, 654, 1110 [NASA ADS] [CrossRef] [Google Scholar]
 Turbet, M., Bolmont, E., Leconte, J., et al. 2018, A&A, 612, A86 [CrossRef] [EDP Sciences] [Google Scholar]
 Unterborn, C. T., Desch, S. J., Hinkel, N., & Lorenzo, A. 2018, Nat. Astron., 2, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Van Grootel, V., Fernandes, C. S., Gillon, M., et al. 2018, ApJ, 853, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Vazan, A., Kovetz, A., Podolak, M., & Helled, R. 2013, MNRAS, 434, 3283 [NASA ADS] [CrossRef] [Google Scholar]
 Vrugt, J. A., Ter Braak, C., Diks, C., et al. 2009, Int. J. Nonlinear Sci. Numer. Simul., 10, 273 [Google Scholar]
 Waldmann, I. P., Tinetti, G., Rocchetto, M., et al. 2015, ApJ, 802, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, S., Wu,D.H., Barclay, T., & Laughlin, G. P. 2017, ArXiv eprints [arXiv:1704.04290] [Google Scholar]
 Ward, W. R. 1997, Icarus, 126, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Wheatley, P. J., Louden, T., Bourrier, V., et al. 2017, MNRAS, 465, L74 [NASA ADS] [CrossRef] [Google Scholar]
 Wisdom, J., & Hernandez, D. M. 2015, MNRAS, 453, 3015 [NASA ADS] [CrossRef] [Google Scholar]
 Wordsworth, R. D., Forget, F., Selsis, F., et al. 2010, A&A, 522, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
This equation is only valid for i = 90°. A discussion for the case i≠90° is given in Gimenez & GarciaPelayo (1983).
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
Updated masses, radii (Delrez et al. 2018), correlation coefficients c_{m,r} of masses and radii as well as the densities and the surface gravity of all seven planets.
Median and standard deviation σ of the semimajor axis a, eccentricity e, argument of periapsis ω, and mean anomaly M of the seven planets.
Planetary characteristics of the TRAPPIST1 planets derived from our 1D numerical climate simulations.
All Figures
Fig. 1 Folded shortcadence 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 shortcadence 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 
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 
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 
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 
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: autocorrelation function of the masses for the ensemble and a single chain. 

In the text 
Fig. 6 Posterior probability distribution between all pairs of parameters, the masses m, the semimajor 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 
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 
Fig. 8 Temporal evolution of ω, showing a fast precession and strong mutual orbital perturbations between the planets. 

In the text 
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 
Fig. 10 Mass–radius diagram for the TRAPPIST1 planets, Earth, and Venus. Curves trace idealised compositions of rocky and waterrich 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 
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 TRAPPIST1 system, like the Galilean moons of Jupiter, share a similar parameter space with the subNeptune and superEarth population. Orbital periods are used to reveal the various subpopulations. 

In the text 
Fig. 12 Masses and the incident fluxes received by the TRAPPIST1 planets, compared to other exoplanets found by the transit method (yellow), via the radialvelocity 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 
Fig. 13 Densities, and the incident fluxes received by the TRAPPIST1 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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.