Interior Characterization in Multiplanetary Systems: TRAPPIST-1

The interior model g(·) assumes planets to be composed of a pure iron core, a silicate mantle comprising the oxides Na₂O–CaO–FeO–MgO–Al₂O₃–SiO₂, a pure water layer, and an isothermal gas layer. Those model parameters ${\boldsymbol{m}}$ that determine the rocky and the water layers are ${r}_{\mathrm{core}}$, ${r}_{\mathrm{core}+\mathrm{mantle}}$, $\mathrm{Fe}/{\mathrm{Si}}_{\mathrm{mantle}}$, $\mathrm{Mg}/{\mathrm{Si}}_{\mathrm{mantle}}$, and m_water.

The structural model for the interior uses self-consistent thermodynamics for core, mantle, high-pressure ice, and water ocean. For the core density profile, we use the equation-of-state (EOS) fit of iron in the hcp (hexagonal close-packed) structure provided by Bouchet et al. (2013) on ab initio molecular dynamics simulations. For the silicate mantle, we compute equilibrium mineralogy and density as a function of pressure, temperature, and bulk composition by minimizing Gibbs free energy (Connolly 2009). For the water layers, we follow Vazan et al. (2013) using a quotidian equation of state (QEOS), and above 44.3 GPa, we use the tabulated EOS from Seager et al. (2007) that is derived from density functional theory (DFT) simulations. Depending on pressure and temperature, the water can be in solid, liquid, or supercritical phase. Water vapor is excluded, since we impose the condition that if water is present, then there must be a gas layer on top that imposes a pressure at least as high as the vapor pressure of water. We assume an adiabatic temperature profile within each layer of core, mantle, and solid water. The surface temperature of the water layer is set equal to the temperature of the bottom of the gas layer. Also, temperature is assumed to be continuous across other layer boundaries.

For the gas layer, we test two different models: an isothermal model (model I) and a fully adiabatic model (model II). Model I assumes a thin, isothermal atmosphere in hydrostatic equilibrium and ideal gas behavior, which is calculated using the scale height model. Those model parameters of ${\boldsymbol{m}}$ that parameterize the gas layer in model I and that we aim to constrain are the pressure at the bottom of the gas layer ${P}_{\mathrm{batm}}$, the mean molecular weight μ, and the mean temperature (parameterized by α; see below). The thickness of the opaque gas layer d_atm in model I is given by

$$\begin{eqnarray}&&{d}_{\mathrm{atm}}=H\mathrm{ln}\displaystyle \frac{{P}_{\mathrm{batm}}}{{P}_{\mathrm{out}}},\end{eqnarray} \tag{ 1 }$$

where the amount of opaque scale heights H is determined by the ratio of ${P}_{\mathrm{batm}}$ and ${P}_{\mathrm{out}}$. ${P}_{\mathrm{out}}$ is the pressure level at the optical photosphere for a transit geometry that we fix to 20 mbar (Fortney et al. 2007). The scale height H is the increase in altitude for which the pressure drops by a factor of e and can be expressed by

$$\begin{eqnarray}&&H=\displaystyle \frac{{T}_{\mathrm{atm}}{R}^{* }}{{g}_{\mathrm{batm}}\mu },\end{eqnarray} \tag{ 2 }$$

where g_batm and T_atm are gravity at the bottom of the atmosphere and mean atmospheric temperature, respectively. R* is the universal gas constant (8.3144598 J mol⁻¹ K⁻¹), and μ is the mean molecular weight. The mass of the atmosphere m_atm is directly related to the pressure P_batm as

$$\begin{eqnarray}&&{m}_{\mathrm{atm}}=4\pi {P}_{\mathrm{batm}}\displaystyle \frac{{r}_{\mathrm{batm}}^{2}}{{g}_{\mathrm{batm}}},\end{eqnarray} \tag{ 3 }$$

where r_batm is the radius at the bottom of the atmosphere.

The atmosphere's constant temperature is defined as

$$\begin{eqnarray}&&{T}_{\mathrm{atm}}=\alpha {T}_{\mathrm{star}}\sqrt{\displaystyle \frac{{R}_{\mathrm{star}}}{2a}},\end{eqnarray} \tag{ 4 }$$

where R_star and T_star are the radius and effective temperature of the host star, respectively, and a is the semimajor axis. The factor α accounts for possible cooling and heating of the atmosphere and can vary between 0.5 and α_max. The upper bound α_max is because there is a physical limit to the amount of warming by greenhouse gases. We approximate α_max for a moist (water-saturated) atmosphere (see Appendix A; Dorn et al. 2017).

The assumption of a constant mean temperature in model I might affect the interior predictions. In order to test how sensitive our results are to this assumption, we propose a second atmosphere model (model II) that assumes a convective adiabat. For model II, the temperature at the bottom of the gas layer, T_batm, is calculated as

$$\begin{eqnarray}&&{T}_{\mathrm{batm}}={T}_{\mathrm{out}}{\left(\displaystyle \frac{{P}_{\mathrm{batm}}}{{P}_{\mathrm{out}}}\right)}^{\kappa },\end{eqnarray} \tag{ 5 }$$

where T_out is calculated with Equation (4), while restricting α to values between 0.5 and 1, which is equivalent to a range of albedos from 0 to 0.94. The exponent κ is equal to 2/(2 + n), with n being the number of degrees of freedom (which we vary randomly between 5 and 6 for diatomic or triatomic gases, thus neglecting vibrational modes). Following a convective adiabat, the thickness of the opaque atmosphere is set to

$$\begin{eqnarray}&&{d}_{\mathrm{atm}}=\displaystyle \frac{{c}_{{\rm{p}}}}{{g}_{\mathrm{batm}}}({T}_{\mathrm{batm}}-{T}_{\mathrm{atm}}),\end{eqnarray} \tag{ 6 }$$

where the heat capacity is ${c}_{{\rm{p}}}=(n/2+1){R}^{* }/\mu$.

In summary, the entire interiors are parameterized by the full set of interior parameters ${\boldsymbol{m}}$: ${r}_{\mathrm{core}}$, ${r}_{\mathrm{core}+\mathrm{mantle}}$, $\mathrm{Fe}/{\mathrm{Si}}_{\mathrm{mantle}}$, $\mathrm{Mg}/{\mathrm{Si}}_{\mathrm{mantle}}$, m_water, ${P}_{\mathrm{batm}}$, α, μ, and for model II also n. We refer to Dorn et al. (2017) for more details on the interior structure model.

The prior distributions of the interior parameters are listed in Table 1. The priors are chosen conservatively. The cubic uniform priors on ${r}_{\mathrm{core}}$ and ${r}_{\mathrm{core}+\mathrm{mantle}}$ reflect equal weighing of masses for both core and mantle. Prior bounds on $\mathrm{Fe}/{\mathrm{Si}}_{\mathrm{mantle}}$ and $\mathrm{Mg}/{\mathrm{Si}}_{\mathrm{mantle}}$ are determined by the abundance proxies. Since iron is distributed between core and mantle, the abundance constraint only sets an upper bound on $\mathrm{Fe}/{\mathrm{Si}}_{\mathrm{mantle}}$, thus making iron-free mantles possible. For the atmospheric models I and II, a ln-uniform prior is set for ${P}_{\mathrm{batm}}$. P_batm,max is estimated by Equation (24) from Ginzburg et al. (2016), who investigated super-Earth atmospheres. A uniform prior in μ equally favors various dominating species (e.g., N₂, CO₂, H₂O), which seems appropriate for anticipated terrestrial-like atmospheres. The factor α incorporates possible cooling (model I and II) and heating (model I) of the atmosphere; it can vary linearly between 0.5 and α_max for the isothermal model I and between 0.5 and 1 for the adiabatic model II. The upper bound α_max is a physical limit to the amount of warming by greenhouse gases and is described in Appendix A. The lower bound of 0.5 is equivalent to an albedo of 0.94, which is similar to the highest albedo estimate among solar system bodies (0.96 for Eris; Sicardy et al. 2011).

Table 1.  Prior Ranges

Parameter	Prior Range	Distribution
${r}_{\mathrm{core}}$	(0.01–1) rmantle	Uniform in ${r}_{\mathrm{core}}^{3}$
$\mathrm{Fe}/{\mathrm{Si}}_{\mathrm{mantle}}$	0–$\mathrm{Fe}/{\mathrm{Si}}_{\mathrm{star}}$	Uniform
$\mathrm{Mg}/{\mathrm{Si}}_{\mathrm{mantle}}$	$\mathrm{Mg}/{\mathrm{Si}}_{\mathrm{star}}$	Gaussian
${r}_{\mathrm{core}+\mathrm{mantle}}$	(0.01–1) R	Uniform in ${r}_{\mathrm{mantle}}^{3}$
mwater	0–0.98 M	Uniform
${P}_{\mathrm{batm}}$	20 mbar–Pbatm,max	Uniform in ln(${P}_{\mathrm{batm}}$)
α (model I)	0.5–αmax	Uniform
α (model II)	0.5–1.	Uniform
μ	2.3–50.0	Uniform in 1/μ
n (model II)	5–6	Uniform

Download table as:  ASCIITypeset image

The data that we use to characterize the interiors of the planets comprise planetary mass and radius, stellar irradiation, and constraints on refractory element abundances. Planetary radii are taken from Delrez et al. (2018), masses and semimajor axes are taken from Grimm et al. (2018), and stellar radius and luminosity are given by Gillon et al. (2017). These data are based on several hundred planet transits that were observed between 2015 September and 2017 March. In addition to the transits presented in earlier studies (de Wit et al. 2016; Gillon et al. 2016, 2017; Luger et al. 2017), Grimm et al. (2018) included additional observations from Spitzer Space Telescope, Kepler, and K2.

For TRAPPIST-1, there is a correlation between planetary radius and mass for the individual planets, as well as a correlation between all seven planetary masses and between all seven planetary radii. In order to better distinguish between the two kinds of correlation, we will use the following nomenclature. Unless otherwise mentioned, the terms correlated data and uncorrelated data refer to the correlation between mass and radius of individual planets, while the terms dependent and independent refer to the correlation between planetary data of the different planets. Planetary masses and radii are shown in Figure 1 and listed in Table 2, where we also list the corresponding correlation c_m,r. Figure 2 shows the interdependency of planetary masses as derived from Grimm et al. (2018).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the constraints on refractory element abundances we have compiled two different scenarios. In the first scenario (U), we use the bulk abundance constraints as suggested by Unterborn et al. (2017), who analyzed F-G-K stars of similar metallicity to TRAPPIST-1 and estimated their statistical median to Fe/Mg = 1.72 ± 0.46 in mass ratios, while fixing the Mg/Si mass ratio to 0.87. Here we adopted these estimates: Fe/Si = 1.49 ± 0.4, Mg/Si = 0.89 ± 0.3. For comparison, the solar mass ratio estimates are Fe/Si = 1.69 and Mg/Si = 0.89 (Lodders et al. 2009). In a second scenario (A), we do not impose any constraints on bulk abundances that are based on stellar estimates. Instead, we analyze the possible range of refractory elements for the densest planet, TRAPPIST-1e, while only using its data of mass, radius, and stellar irradiation. The estimated ranges of Mg/Si_T1e and Fe/Si_T1e are subsequently used as input constraints under the premise that planet e is analogous in its relative refractory element budget to the rocky interiors of all planets. This excludes a priori the possibility of planet e being a Mercury-type planet.

In total, we discuss five different data scenarios, labeled U, A, UCM, UHM, and UHMR, which we analyze in detail. All data scenarios are summarized in Table 3. Scenarios U and A differ in terms of the abundance constraints used; in scenario UCM, we additionally account for the fact that planetary masses of different planets depend on each other; for UHM and UHMR, we use hypothetical reduced uncertainties on planetary masses and radii. All scenarios help to determine the value of the different data with respect to making interior predictions.

Our characterization scheme involves two main steps. The first step comprises the interior characterization of individual planets b–h using the generalized Bayesian inference scheme of Dorn et al. (2017). The outputs are posterior distributions for each of the planets that do not depend on each other. In a second step, we account for the interdependency of the different planetary data. In order to do so, we have developed a resampling scheme that yields posterior distributions for each of the planets that depend on each other. Thereby, we take all available information of the mass–radius data into account.

In the following, we first introduce the formulation of Tarantola & Valette (1982) on inference problems, and then we describe each of the two main steps. Our scheme is summarized in Figure 3.

Zoom In Zoom Out Reset image size

3.4.1. Bayes's Theorem and the Tarantola–Valette Formulation

Formally, the posterior distribution based on Bayes's theorem for a fixed model parameterization and conditional on data is given by

$$\begin{eqnarray}&&p({\boldsymbol{m}}| {\boldsymbol{d}})=\displaystyle \frac{{\text{}}p({\boldsymbol{m}})L({\boldsymbol{m}}| {\boldsymbol{d}})}{p({\boldsymbol{d}})}\,,\end{eqnarray} \tag{ 7 }$$

where p(${\boldsymbol{m}}$) represents prior information on model parameter ${\boldsymbol{m}}$, $L({\boldsymbol{m}}| {\boldsymbol{d}})$ is the likelihood function, and p(${\boldsymbol{d}}$) signifies the evidence and is here a normalizing constant. The likelihood function is defined as

$$\begin{eqnarray}\begin{array}{rcl}L({\boldsymbol{m}}| {\boldsymbol{d}}) & = & \displaystyle \frac{1}{{(2\pi )}^{N/2}{(\det {\rm{\Sigma }})}^{1/2}}\\ & & \times \ \exp \left(-\displaystyle \frac{1}{2}{(g({\boldsymbol{m}})-{\boldsymbol{d}})}^{T}{{\rm{\Sigma }}}^{-1}(g({\boldsymbol{m}})-{\boldsymbol{d}})\right),\end{array}\end{eqnarray} \tag{ 8 }$$

where Σ is the data covariance matrix and det(Σ) denotes the determinant of Σ. For uncorrelated data errors, Σ is a diagonal matrix. Here the correlation between planetary mass and radius is prescribed by the off-diagonal entries of Σ.

There is a different formulation of inference problems by Tarantola & Valette (1982) that we will use in the following. In this formulation, the posterior p(${\boldsymbol{d}}$, ${\boldsymbol{m}}$) is expressed in the joint space of data and model parameters (${ \mathcal D }$ and ${ \mathcal M }$) as

$$\begin{eqnarray}&&p({\boldsymbol{d}},{\boldsymbol{m}})\equiv (\rho \wedge \theta )({\boldsymbol{d}},{\boldsymbol{m}})=\displaystyle \frac{\rho ({\boldsymbol{d}},{\boldsymbol{m}})\theta ({\boldsymbol{d}},{\boldsymbol{m}})}{z({\boldsymbol{d}},{\boldsymbol{m}})}\,,\end{eqnarray} \tag{ 9 }$$

where the prior ρ, the forward model density θ, and the null information density z are also defined in the joint space ${ \mathcal D }\times { \mathcal M }$. The null information is a normalizing constant. The posterior of the model parameters evaluated in the model space is then

$$\begin{eqnarray}&&{p}_{m}({\boldsymbol{m}})={\int }_{{ \mathcal D }}p({\boldsymbol{d}},{\boldsymbol{m}})d{\boldsymbol{d}}.\end{eqnarray} \tag{ 10 }$$

In Appendix B, we show that the Tarantola–Valette formulation is equivalent to Bayes's theorem under simplified conditions.

3.4.2. Bayesian Inference of Individual Planets

3.4.3. Resampling Scheme to Account for the Interdependency of Data from Different Planets

The above-described Bayesian inference analysis yields posterior samples in the model and data space (${ \mathcal M }$ and ${ \mathcal D }$), p(${\boldsymbol{d}}$, ${\boldsymbol{m}}$). Note that the posterior data p_d(${\boldsymbol{d}}$) can be different from the prior data ρ_d(${\boldsymbol{d}}$). The data prior is the distribution of the noise-contaminated data, which is inherited by the given distribution of observational uncertainties. The posterior data are recomputed data, taking into account constraints from the model ${\boldsymbol{m}}$. Prior data and posterior data are therefore not necessarily identical. The aim of our resampling is to preserve the differences between posterior and prior (which is added by model ${\boldsymbol{m}}$), while adding the information of the interdependency of prior data but not the prior data themselves (since these have been already used in the Bayesian inference).

The posteriors p(${\boldsymbol{d}}$, ${\boldsymbol{m}}$) are obtained for each planet individually. Consequently, the obtained data posteriors for a specific planet evaluated in ${ \mathcal D }$, p_d(${\boldsymbol{d}}$), do not depend on the posteriors of other planets, but they are independent. The interdependency of planetary masses as described by Grimm et al. (2018) quantifies the dependencies of data priors between the planets. Formally, we can say that ${\rho }_{{d}_{{ij}}}({{\boldsymbol{d}}}_{{ij}})\ne {\rho }_{{d}_{i}}({{\boldsymbol{d}}}_{i}){\rho }_{{d}_{j}}({{\boldsymbol{d}}}_{j})$, where i, j are generic planet indices i, j ∈ {b, c, d, e, f, g, h}. In the case in which the individual data priors are independent from each other, we can write the data posteriors as

$$\begin{eqnarray}&&{p}_{{d}_{{ij}}}^{\mathrm{indep}}({\boldsymbol{d}})={L}_{{m}_{i}}({{\boldsymbol{d}}}_{i}){L}_{{m}_{j}}({{\boldsymbol{d}}}_{j}){\rho }_{{d}_{i}}({{\boldsymbol{d}}}_{i}){\rho }_{{d}_{j}}({{\boldsymbol{d}}}_{j})\,,\end{eqnarray} \tag{ 11 }$$

while for dependent data priors the posteriors are

$$\begin{eqnarray}\begin{array}{rcl}{p}_{{d}_{{ij}}}^{\mathrm{dep}}({\boldsymbol{d}}) & = & {L}_{{m}_{{ij}}}({{\boldsymbol{d}}}_{i},{{\boldsymbol{d}}}_{j}){\rho }_{{d}_{{ij}}}({{\boldsymbol{d}}}_{i},{{\boldsymbol{d}}}_{j})\\ & = & \ {L}_{{m}_{i}}({{\boldsymbol{d}}}_{i}){L}_{{m}_{j}}({{\boldsymbol{d}}}_{j}){\rho }_{{d}_{{ij}}}({{\boldsymbol{d}}}_{i},{{\boldsymbol{d}}}_{j}).\end{array}\end{eqnarray} \tag{ 12 }$$

If ${L}_{{m}_{{ij}}}({{\boldsymbol{d}}}_{i},{{\boldsymbol{d}}}_{j})={L}_{{m}_{i}}({{\boldsymbol{d}}}_{i}){L}_{{m}_{j}}({{\boldsymbol{d}}}_{j})$, then

$$\begin{eqnarray}&&{p}_{{d}_{{ij}}}^{\mathrm{dep}}({{\boldsymbol{d}}}_{{ij}})={p}_{{d}_{{ij}}}^{\mathrm{indep}}({{\boldsymbol{d}}}_{{ij}})\displaystyle \frac{{\rho }_{{d}_{{ij}}}({{\boldsymbol{d}}}_{i},{{\boldsymbol{d}}}_{j})}{{\rho }_{{d}_{i}}({{\boldsymbol{d}}}_{i}){\rho }_{{d}_{j}}({{\boldsymbol{d}}}_{j})}.\end{eqnarray} \tag{ 13 }$$

The condition that ${L}_{{m}_{{ij}}}({{\boldsymbol{d}}}_{i},{{\boldsymbol{d}}}_{j})={L}_{{m}_{i}}({{\boldsymbol{d}}}_{i}){L}_{{m}_{j}}({{\boldsymbol{d}}}_{j})$ can be understood from the following. L(d) is a measure of the likelihood of the specific data set d, based solely on the prior information about the models m = g⁻¹(d) that map into d. If g is bijective (the inverse problem has a unique best-fitting solution), ${L}_{m}({\boldsymbol{d}})={\rho }_{m}({g}^{-1}({\boldsymbol{d}}))$. If the problem is nonunique, the right-hand side of this equation expresses an integral/sum over all models that map into d. For simplicity, consider a case where our data consist of two (possibly dependent) parts, d₁ and d₂. In this case, the data likelihood can be expressed as

$$\begin{eqnarray}&&{L}_{m}({\boldsymbol{d}})={L}_{m}({{\boldsymbol{d}}}_{1},{{\boldsymbol{d}}}_{2})={\rho }_{{m}_{1},{m}_{2}}({g}^{-1}({{\boldsymbol{d}}}_{1},{{\boldsymbol{d}}}_{2})).\end{eqnarray} \tag{ 14 }$$

However, when m₁ and m₂ are a priori (statistically) independent and m_i is functionally related only to d_i (for i = 1, 2), we get

$$\begin{eqnarray}&&{\rho }_{{m}_{1},{m}_{2}}({g}^{-1}({{\boldsymbol{d}}}_{1},{{\boldsymbol{d}}}_{2}))={\rho }_{{m}_{1}}({g}^{-1}({{\boldsymbol{d}}}_{1},{{\boldsymbol{d}}}_{2}))\cdot {\rho }_{{m}_{2}}({g}^{-1}({{\boldsymbol{d}}}_{1},{{\boldsymbol{d}}}_{2}))\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&=\,{\rho }_{{m}_{1}}({g}^{-1}({{\boldsymbol{d}}}_{1}))\cdot {\rho }_{{m}_{2}}({g}^{-1}({{\boldsymbol{d}}}_{2})),\end{eqnarray} \tag{ 16 }$$

leading to

$$\begin{eqnarray}&&{L}_{m}({{\boldsymbol{d}}}_{1},{{\boldsymbol{d}}}_{2})={L}_{{m}_{1}}({{\boldsymbol{d}}}_{1})\cdot {L}_{{m}_{2}}({{\boldsymbol{d}}}_{2}).\end{eqnarray} \tag{ 17 }$$

So much for the theory.

In practice, the resampling works as follows. The inference analysis of individual planets yields independent data posteriors ${p}_{{d}_{{ij}}}^{(\mathrm{indep})}({{\boldsymbol{d}}}_{i},{{\boldsymbol{d}}}_{j})={p}_{{d}_{i}}({{\boldsymbol{d}}}_{i}){p}_{{d}_{j}}({{\boldsymbol{d}}}_{j})$. For a total of K samples, we have ${d}_{i}^{1}$, ${d}_{i}^{2}$, ..., ${d}_{i}^{K}$ from ${p}_{{d}_{i}}$ and ${d}_{j}^{1}$, ${d}_{j}^{2}$, ..., ${d}_{j}^{K}$ from ${p}_{{d}_{j}}$. Thus, any pair (e.g., $({d}_{i}^{1},{d}_{j}^{1})$, $({d}_{i}^{2},{d}_{j}^{2})$, ... $({d}_{i}^{K},{d}_{j}^{K})$) is a sample from ${p}_{{d}_{i}}({{\boldsymbol{d}}}_{i}){p}_{{d}_{j}}({{\boldsymbol{d}}}_{j})$. We shall resample from ${p}_{{d}_{i}}({{\boldsymbol{d}}}_{i}){p}_{{d}_{j}}({{\boldsymbol{d}}}_{j})$ such that the new samples are drawn from ${p}_{{d}_{{ij}}}^{\mathrm{dep}}({{\boldsymbol{d}}}_{i},{{\boldsymbol{d}}}_{j})$. In order to do so, we calculate a factor Q that signifies a maximum difference between the data priors of the data-independent and data-dependent cases:

$$\begin{eqnarray}&&Q={\max }_{k}({\rho }_{{d}_{{ij}}}^{\mathrm{dep}}({{\boldsymbol{d}}}_{i}^{k},{{\boldsymbol{d}}}_{j}^{k})/{\rho }_{{d}_{i}}({{\boldsymbol{d}}}_{i}^{k}){\rho }_{{d}_{j}}({{\boldsymbol{d}}}_{j}^{k})).\end{eqnarray} \tag{ 18 }$$

For k = 1, ..., K,

• 1.  
  we generate a random number w ∈ [0,1];
• 2.  
  if $w\lt {\rho }_{{d}_{{ij}}}^{\mathrm{dep}}({{\boldsymbol{d}}}_{i}^{k},{{\boldsymbol{d}}}_{j}^{k})/(Q\times {\rho }_{{d}_{i}}({{\boldsymbol{d}}}_{i}^{k}){\rho }_{{d}_{j}}({{\boldsymbol{d}}}_{j}^{k}))$, then ${({{\boldsymbol{d}}}_{i}^{k},{{\boldsymbol{d}}}_{j}^{k})=({{\boldsymbol{d}}}_{i}^{k},{{\boldsymbol{d}}}_{j}^{k})}^{\mathrm{dep}}$; otherwise, discard the sample.

Thereby, we obtain samples ${({{\boldsymbol{d}}}_{i}^{k},{{\boldsymbol{d}}}_{j}^{k})}^{\mathrm{dep}}$ that represent the data posterior ${p}_{{d}_{{ij}}}^{\mathrm{dep}}({{\boldsymbol{d}}}_{i},{{\boldsymbol{d}}}_{j})$. Any interior model pairs ${({{\boldsymbol{m}}}_{i}^{k},{{\boldsymbol{m}}}_{j}^{k})}^{\mathrm{dep}}$ that generated ${({{\boldsymbol{d}}}_{i}^{k},{{\boldsymbol{d}}}_{j}^{k})}^{\mathrm{dep}}$ are samples from the model posterior ${p}_{{m}_{{ij}}}^{\mathrm{dep}}({{\boldsymbol{m}}}_{i},{{\boldsymbol{m}}}_{j})$.

In other words, we can understand the resampling as follows. We randomly combine samples from the analysis of individual planets i and j, which represents a set of paired samples $({d}_{i}^{k},{d}_{j}^{k})$ drawn from the data posterior under the premise that data of different planets are independent. In order to obtain the data posterior samples ${({d}_{i}^{k},{d}_{j}^{k})}^{\mathrm{dep}}$ that incorporate the dependencies of data from different planets, we reject some of the original samples $({d}_{i}^{k},{d}_{j}^{k})$. Those samples for which the ratio between dependent data prior and independent data prior is large (close to 1 or larger) are likely accepted as samples from ${({d}_{i}^{k},{d}_{j}^{k})}^{\mathrm{dep}}$, while they are rejected if the ratio is closer to 0. This procedure is done for all combinations of planets i and j.

Alternatively, it is in principle possible to do an inference analysis by using all data and their correlations and interdependencies between planets. However, our scheme elegantly separates the information of correlated mass and radius of individual planets from information about interdependent data of different planets. This allows us to better control each analysis step and to be more flexible with testing different data scenarios without the need to run a full inference analysis for the entire system. Furthermore, our scheme can be generally applied to other inference problems.

In Figure 4, we show the consequences of different bulk constraints on mass and radius of TRAPPIST-1e. The prior data of measured mass and radius are highlighted in blue. The use of the suggested stellar constraint of Unterborn et al. (2017) (scenario U) is not able to explain well the relatively high bulk density of TRAPPIST-1e. This is evident from the difference between the prior data distribution (blue) and the posterior distribution of scenario U (green).

Zoom In Zoom Out Reset image size

In a second scenario N, we do not impose any abundance constraints. Instead, we analyze the possible range of bulk abundances of TRAPPIST-1e, while only using the data of mass, radius, and stellar irradiation. The predicted bulk mass ratios are high, with large 1σ uncertainties: Fe/Si_T1e = 11.2 ± 5.7 and Mg/Si_T1e = 5.7 ± 3.7. The data of only mass and radius (N) allow high relative bulk iron abundances, since iron incorporated in silicates and oxides in the mantle has a limited influence on bulk density. These abundance ranges (Fe/Si_T1e and Mg/Si_T1e) are subsequently used as input constraints to analyze the other planets in scenario A.

The ratios Fe/Si_T1e and Mg/Si_T1e are much higher than those estimated for Earth and other solar terrestrial planets. However, the large uncertainties on the ratios also allow for Earth-like rocky compositions. The fact that these uncertainties estimated from mass and radius alone are high illustrates the large degeneracy of rocky interiors, i.e., very different combinations of mantle compositions and core sizes can explain a given planetary mass and radius. The suggested stellar abundance proxy of Unterborn et al. (2017), with a relatively small uncertainty of 30%, is partly incompatible with mass and radius of TRAPPIST-1e. In the discussion (Section 5), we provide two interpretations for the interior of planet e.

Our calculated interior predictions account for the variability in layer thicknesses (of core, mantle, ice/ocean, and gas), layer compositions (of mantle and gas), and thermal states. In the following, we present 1D and 2D marginalized distributions of estimated posteriors. Note that the 1σ uncertainty regions differ between a distribution that is marginalized over one or two dimensions. The uncertainty regions are generally smaller for 1D marginalized distributions. Stated uncertainty regions refer to 1D marginalized distributions.

In Figure 5 we show the predicted ranges of radius fractions for gas and water layers for the scenario U, in which we use the stellar abundance proxy as suggested by Unterborn et al. (2017). Possible radius fractions of gas envelopes range between 0% and 5%, while water layers can contribute up to 30% to the total radius. The uncertainty ranges on possible amounts of water are large. This is because of the degeneracy with gas envelope thicknesses (Figure 5), but mainly the size of the rocky interiors (Figure 6) is large. There is a strong correlation between the possible amount of water and the size of the rocky interior, as expected.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the following, we discuss how much the predicted water mass fractions depend on model assumptions and considered data.

4.2.1. Influence of the Choice of Abundance Proxy

Figure 7 illustrates how the choice of abundance constraints influences the predicted amount of water. If the stellar proxy (U) is used, rocky interiors are less dense compared to the proxy that is based on TRAPPIST-1e (A). Consequently, the predicted amount of water is smaller in scenario U, since less water can be added on top of low-density rocky interiors while still fitting mass and radius. There are differences in the median predicted water mass fraction that range from 25% up to 50% between U and A. For all but planet e, the differences are largely around 30%. For planet e, it is the discrepancy between mass and radius versus stellar proxy that leads to such large changes in m_water/M_p between the scenarios U and A. Table 4 summarizes the estimated interior parameters and their 1σ uncertainties.

Zoom In Zoom Out Reset image size

Table 4.  Interior Parameter Estimates with 1σ Uncertainties of the 1D Marginalized Posterior Distributions

Planet:	TRAPPIST-1b			TRAPPIST-1c			TRAPPIST-1d
Scenario:	U	A	UCM	U	A	UCM	U	A	UCM
${r}_{\mathrm{core}}$/${r}_{\mathrm{core}+\mathrm{mantle}}$	${0.39}_{-0.11}^{+0.09}$	${0.50}_{-0.14}^{+0.12}$	${0.39}_{-0.11}^{+0.09}$	${0.40}_{-0.11}^{+0.08}$	${0.52}_{-0.13}^{+0.11}$	${0.40}_{-0.11}^{+0.08}$	${0.39}_{-0.11}^{+0.09}$	${0.51}_{-0.14}^{+0.12}$	${0.39}_{-0.12}^{+0.08}$
${r}_{\mathrm{core}+\mathrm{mantle}}$ $/{R}_{{\rm{p}}}$	${0.84}_{-0.06}^{+0.06}$	${0.79}_{-0.07}^{+0.07}$	${0.82}_{-0.07}^{+0.08}$	${0.92}_{-0.05}^{+0.04}$	${0.86}_{-0.06}^{+0.06}$	${0.91}_{-0.05}^{+0.04}$	${0.84}_{-0.06}^{+0.05}$	${0.79}_{-0.06}^{+0.06}$	${0.84}_{-0.06}^{+0.06}$
${m}_{\mathrm{water}}/{M}_{{\rm{p}}}$	${0.15}_{-0.07}^{+0.08}$	${0.19}_{-0.08}^{+0.09}$	${0.16}_{-0.09}^{+0.10}$	${0.06}_{-0.04}^{+0.05}$	${0.10}_{-0.05}^{+0.06}$	${0.06}_{-0.04}^{+0.05}$	${0.14}_{-0.06}^{+0.08}$	${0.18}_{-0.07}^{+0.07}$	${0.14}_{-0.07}^{+0.09}$
${r}_{\mathrm{water}}/{R}_{{\rm{p}}}$	${0.15}_{-0.06}^{+0.07}$	${0.19}_{-0.07}^{+0.07}$	${0.16}_{-0.08}^{+0.07}$	${0.07}_{-0.04}^{+0.05}$	${0.12}_{-0.06}^{+0.06}$	${0.08}_{-0.04}^{+0.05}$	${0.14}_{-0.06}^{+0.06}$	${0.18}_{-0.06}^{+0.06}$	${0.13}_{-0.06}^{+0.07}$
${r}_{\mathrm{gas}}/{Rp}$	${0.01}_{-0.01}^{+0.01}$	${0.01}_{-0.01}^{+0.01}$	${0.01}_{-0.01}^{+0.01}$	${0.01}_{-0.01}^{+0.01}$	${0.01}_{-0.01}^{+0.01}$	${0.01}_{-0.01}^{+0.01}$	${0.02}_{-0.01}^{+0.02}$	${0.02}_{-0.01}^{+0.02}$	${0.02}_{-0.01}^{+0.02}$
Planet:	TRAPPIST-1e			TRAPPIST-1f			TRAPPIST-1g
Scenario:	U	A	UCM	U	A	UCM	U	A	UCM
${r}_{\mathrm{core}}$/${r}_{\mathrm{core}+\mathrm{mantle}}$	${0.45}_{-0.10}^{+0.06}$	${0.58}_{-0.11}^{+0.09}$	${0.45}_{-0.10}^{+0.07}$	${0.39}_{-0.11}^{+0.08}$	${0.52}_{-0.13}^{+0.12}$	${0.39}_{-0.11}^{+0.08}$	${0.39}_{-0.11}^{+0.08}$	${0.50}_{-0.14}^{+0.12}$	${0.39}_{-0.11}^{+0.09}$
${r}_{\mathrm{core}+\mathrm{mantle}}$ $/{R}_{{\rm{p}}}$	${0.97}_{-0.02}^{+0.02}$	${0.93}_{-0.05}^{+0.04}$	${0.97}_{-0.02}^{+0.02}$	${0.91}_{-0.03}^{+0.03}$	${0.86}_{-0.05}^{+0.05}$	${0.91}_{-0.03}^{+0.03}$	${0.86}_{-0.03}^{+0.03}$	${0.81}_{-0.05}^{+0.05}$	${0.86}_{-0.03}^{+0.03}$
${m}_{\mathrm{water}}/{M}_{{\rm{p}}}$	${0.02}_{-0.01}^{+0.02}$	${0.04}_{-0.03}^{+0.04}$	${0.02}_{-0.01}^{+0.02}$	${0.07}_{-0.03}^{+0.03}$	${0.11}_{-0.05}^{+0.05}$	${0.07}_{-0.03}^{+0.03}$	${0.13}_{-0.04}^{+0.04}$	${0.17}_{-0.06}^{+0.05}$	${0.13}_{-0.04}^{+0.04}$
${r}_{\mathrm{water}}/{R}_{{\rm{p}}}$	${0.02}_{-0.01}^{+0.02}$	${0.06}_{-0.04}^{+0.05}$	${0.02}_{-0.01}^{+0.02}$	${0.08}_{-0.03}^{+0.03}$	${0.13}_{-0.05}^{+0.05}$	${0.08}_{-0.03}^{+0.03}$	${0.13}_{-0.04}^{+0.03}$	${0.18}_{-0.06}^{+0.05}$	${0.13}_{-0.04}^{+0.04}$
${r}_{\mathrm{gas}}/{Rp}$	${0.01}_{-0.00}^{+0.01}$	${0.01}_{-0.01}^{+0.01}$	${0.01}_{-0.00}^{+0.01}$	${0.01}_{-0.01}^{+0.01}$	${0.01}_{-0.01}^{+0.01}$	${0.01}_{-0.01}^{+0.01}$	${0.01}_{-0.01}^{+0.01}$	${0.01}_{-0.01}^{+0.01}$	${0.01}_{-0.01}^{+0.01}$
Planet:	TRAPPIST-1h
Scenario:	U	A	UCM
${r}_{\mathrm{core}}$/${r}_{\mathrm{core}+\mathrm{mantle}}$	${0.40}_{-0.11}^{+0.08}$	${0.52}_{-0.14}^{+0.11}$	${0.40}_{-0.12}^{+0.08}$
${r}_{\mathrm{core}+\mathrm{mantle}}$ $/{R}_{{\rm{p}}}$	${0.89}_{-0.07}^{+0.06}$	${0.84}_{-0.07}^{+0.07}$	${0.87}_{-0.07}^{+0.06}$
${m}_{\mathrm{water}}/{M}_{{\rm{p}}}$	${0.09}_{-0.06}^{+0.09}$	${0.12}_{-0.07}^{+0.08}$	${0.11}_{-0.06}^{+0.09}$
${r}_{\mathrm{water}}/{R}_{{\rm{p}}}$	${0.10}_{-0.06}^{+0.08}$	${0.14}_{-0.07}^{+0.08}$	${0.11}_{-0.06}^{+0.08}$
${r}_{\mathrm{gas}}/{Rp}$	${0.01}_{-0.01}^{+0.02}$	${0.01}_{-0.01}^{+0.02}$	${0.01}_{-0.01}^{+0.02}$

Download table as:  ASCIITypeset image

4.2.2. Accounting for the Interdependency of Planetary Masses

Inferred planetary masses by Grimm et al. (2018) are significantly correlated between different planets (Figure 2). We account for this interdependency of planetary masses by resampling the posterior interior samples of scenario U with our new resampling scheme (Section 3.4.3) that yields the posterior distribution UCM. The difference between the posterior distribution U and UCM demonstrates the information contained in the interdependency of planetary masses. Figure 8 shows that the differences are largest for planets b and h, with up to 20% difference in both the median water mass fraction and the corresponding uncertainty (Table 4).

Zoom In Zoom Out Reset image size

The information contained in the interdependency of planetary data depends on the specific correlation between the planets. In Figure 9 we show how the distribution of water mass fraction for planet b changes while adding step by step the information on the correlation between the planetary masses. For planet b, the information gained by accounting for interdependency of masses is mainly kept in the interdependency with planets e and h. In general, how much information can be gained depends on the actual level of interdependency of planet pairs. While accounting for interdependent data, water mass fractions are shifting toward lower values, and uncertainties increase.

Zoom In Zoom Out Reset image size

Why do uncertainties on estimated water mass fractions increase? The reason is illustrated in Figure 10, showing the independent posteriors (U) and the prior ratio (${\rho }_{{d}_{{ij}}}({{\boldsymbol{d}}}_{i},{{\boldsymbol{d}}}_{j})/Q{\rho }_{{d}_{i}}({{\boldsymbol{d}}}_{i}){\rho }_{{d}_{j}}({{\boldsymbol{d}}}_{j})$ from Equation (13)), which is used to resample from the independent posteriors in order to obtain dependent posteriors. The ratio is highest at the two tails of the 2D distributions, where planetary masses are both small and large. These tails of the independent mass distribution are preferentially resampled. For a given radius, high masses are characterized by low water contents and vice versa. In consequence, the distribution of possible water mass fractions broadens.

Zoom In Zoom Out Reset image size

Why are higher water mass fractions preferentially resampled? For some of the planets, the highest possible masses cannot be well described by the stellar proxy (U). This is the case for planet e (as discussed earlier for Figure 4) but also planets c and h (Figure 10, bottom panel). It is not the case, for example, for planets f and g (Figure 10, top panel). The consequence is clear from the posteriors (U) that are shifted to smaller planetary masses compared to the prior data (for e, h, and c). During the resampling of the posterior U, only those samples survive for which the prior ratio is large. This is the case for lower masses (Figure 10, bottom panel) that are characterized by higher water contents (for a given radius).

The difference between U and UCM is due to not only the partial incompatibility of planetary data (mass and radius versus the stellar proxy) of planets e, h, and c but also the fact that the 2D mass distributions (Figure 2) are not perfect ellipses.

4.2.3. Influence of Atmospheric Model

Figure 7 shows that the differences in water mass fractions are small when we change from an isothermal to a fully adiabatic temperature model for the gas layer. In the isothermal model I, we have introduced a parameter α that accounts for cooling and heating of an atmosphere. However, temperatures are largely limited to 400 K, and thus, especially for the innermost planets, this model I is not able to produce temperatures as high as suggested by Grimm et al. (2018) (surface temperatures of 750–1500 K for planet b). Model II overcomes this limitation but overpredicts the temperatures in the gas layer by assuming a fully convective gas layer. In consequence, surface temperatures in model II can be very high and can range, for example, on planet b from 500 to 10,000 K for surface pressures of 1–10⁴ bars. Thus, while temperatures in model I are underestimated, they are overestimated in model II. This is illustrated in Figure 11, showing pressure–temperature profiles for an Earth-like atmosphere. Global climate models that solve for radiative transfer as presented in Turbet et al. (2018) yield much more realistic surface temperatures. However, we find that there is only limited influence due to the choice of gas layer model on the estimated mass of the underlying water layer. Median estimates generally vary by 10%–15%, expect for planets b and c with 25%. The differences are relatively small compared to other influences (Sections 4.2.1 and 5.1).

Zoom In Zoom Out Reset image size

Estimated water mass fractions are smaller for model II compared to model I. This is because the lower gas envelope temperatures in model I result in smaller scale heights and thin gas envelopes, and thus larger amounts are required to fit planet mass and radius. The solutions given the two models I and II represent end-member solutions.

4.2.4. Relevance of Data Uncertainty

Figure 12 illustrates how much our ability to constrain water mass fraction would improve by more precise estimates of mass and radius. An improved estimate on planet mass (70% of the nominal uncertainty) only has generally marginal influence on the median estimate of m_water/M_p, but it reduces the 1σ confidence interval by 5%–20%. Improved estimates on both planet mass and radius (70% of the nominal uncertainty) reduce the 1σ confidence interval by 10%–35%. The largest improvements are seen for planets e and h, while the smallest improvements are seen for planets f and g. The information gained by improved data uncertainty is small compared to the different choices of abundance proxies (Figure 7), as well as the information kept in the interdependency of planetary data (Figure 8).

Zoom In Zoom Out Reset image size

Recent announcements (Demory 2018) have indicated that previously unavailable observations for the TTV analysis of Grimm et al. (2018) provide a higher accuracy with changes in planetary masses within ±24%. In Figure 13, we show that consequences for the estimated water mass fractions are large. Planet e would not be the odd case of a super-Mercury. All planets could be consistent with an increasing water mass fraction with orbital period (within 2%–12% with respect to 1σ errors) or a uniform water mass fraction of 7% (within 1σ errors). Such large amounts of water would exclude Earth-like habitability even for the temperate planets since their rocky interiors are separated by icy layers from the water oceans in that case. In contrast, the published data estimates (Grimm et al. 2018) suggest no clear trend of water mass fraction with orbital distance as discussed below (Figure 13).

Zoom In Zoom Out Reset image size

Specifically for the TRAPPIST-1 planets, the trend of changing water mass fractions with systematic shifts in planetary masses can be described as follows. A shift of +25% in planetary mass can reduce the water mass fraction by 65% (i.e., the reference water mass fraction multiplied by a factor of 0.35). For shifts of +20%, +10%, +5%, and −15% the factor is roughly 0.4, 0.55, 0.75, and 2.6, respectively.

In general, biases from the TTV analysis are extremely difficult to quantify because they stem from limited observations. Biases from TTV can be caused by unresolved degeneracies between eccentricity, mass, and arguments of perihelion of the different planets. Limited observation time can cause the TTV analysis to prefer solutions that depart from the true solution, or the parameter search is trapped in a local minimum of the sampling parameter space. Also strong tidal effects, nearby perturbing stars, or even yet-undetected additional planets can introduce a TTV signal bias. Further TTV analysis on a larger number of transit observations is needed to increase the accuracy of planetary masses.

The degeneracies within structure and composition of rocky interiors stem from the fact that we allow for variable silicate mantle compositions and core sizes. Independent constraints on element abundances as taken from stellar proxies are highly valuable for refined rocky interior estimates (Dorn et al. 2015). Here, however, no direct estimates from the faint host star are available. Instead, we have adopted two different abundance proxies (see Table 3): a stellar proxy (scenario U) as derived from F-G-K stars of similar metallicities (Unterborn et al. 2017), and a proxy based on only mass and radius of the densest planet of the system, planet e (scenario A).

Depending on the abundance proxies used for refractory elements (U or A), the interior estimates change with significant influences on the predicted possible water budgets (Section 4.2.1). Any decision on whether scenario U or A is more likely depends on the interpretation of planet e.

If TRAPPIST-1e were a super-Mercury-type planet, it would not represent the rocky analog for all planets. In this case, scenario A would be best only to describe planet e, while all other planets would be best described by scenario UCM (in case the stellar proxy would be excluded for planet e). A super-Mercury-type interior could be explained by a high-energy giant impact (Benz et al. 1988; Stewart & Leinhardt 2009; Marcus et al. 2010). This would also allow us to explain why inferred interiors for planet e are poorest in volatiles, although its neighboring planets are volatile-rich. A mantle-stripping impact would require large impact velocities or small impact angles (Stewart & Leinhardt 2009), or an impactor larger than planet e, which does not exist in the system. It remains to be studied whether impacts are a reasonable scenario.

If TRAPPIST-1e is not a super-Mercury-type planet, there are two interpretations possible. First, TRAPPIST-1e can indeed be a rocky analog of all the planets regarding the relative ratios of rock-forming elements. In this case, all planets have compositions that are not Earth-like, implying that the iron was more abundant compared to solar compositions. In this case, all planetary interiors are best described by scenario A.

Second, let us assume that TRAPPIST-1e is indeed similar to an Earth interior and thus well represented by U, but its mass and radius estimates are affected by systematic biases (as discussed in Section 5.1). In this case, scenario U may well describe planet e, while dismissing high bulk density interiors. Thus, many mass–radius realizations are incompatible for planet e in scenario U. Given that all planetary masses (radii) are correlated to the mass (radius) of planet e, there are also many planetary masses (radii) of the other planets that are incompatible with U. In consequence, accounting for the interdependencies of planetary data improves interior estimates. If this case is true, all interiors are well described by scenario UCM.

The layers of water can generally be in liquid, solid, and even supercritical state. Water vapor is excluded, as we impose that any water layer requires a gaseous layer on top that imposes a pressure at least as high as the vapor pressure of water. For simplification, we use pure water and neglect other forms of ices (e.g., CO₂). For all but planet e, most of the possible planetary interiors are characterized by ice layers separating water oceans from rocky interiors. This is due to the fact that either the large amounts of water allow pressures high enough to form high-pressure ices or temperatures are sufficiently cool (especially for the outer planets f, g, and h). We note that the thermal states of the planets are unknown and significant tidal heating may allow for more extended thicknesses of liquid water. Planet e is the volatile-poor exception among the planets. Its possible tiny amounts of water and temperate conditions allow for liquid water on the surface. Also for planets c and h, the uncertainty on mass and radius allows for negligibly thin volatile layers.

In Figure 14, we show the 2σ uncertainty distribution for scenario U and all planets in comparison with mean estimates for Earth and Venus, but also Titan and icy satellites, which probably formed outside the water ice line: Europa, Ganymede, Callisto, Triton, Pluto, and Charon. The water mass fractions of the shown solar system objects range from 0% to 30% (for a review, see Nimmo & Pappalardo 2016) and are comparable to our inferred ranges of possible water mass fractions. Water mass fractions of solar system objects on the order of 1000 km or smaller can have higher water budgets. (Saturn's moon Tethys, with a radius of 1060 km, is almost entirely composed of water, with a bulk density of 0.98 g cm⁻³.) Formation models that simulate the structure and evolution of disks (Alibert & Benz 2017) predicted planets to be more volatile-rich around low-mass stars, since the snow line is located closer to the stars compared to the solar system. Indeed, the TRAPPIST-1 system is significantly more volatile-rich than the terrestrial solar system planets.

Zoom In Zoom Out Reset image size

Our results suggest an upper water mass fraction of not much larger than 30%. The amount of water loss (Bourrier et al. 2017) or water delivery (Kral et al. 2018) during the long-term evolution of the planets after formation is small compared to our predicted uncertainties on water budgets. While water budgets can be altered by 1–100 Earth oceans, the predicted water budgets range from tens to thousands of Earth oceans, with uncertainties on the scale of hundreds to thousands of Earth oceans. Therefore, our inferred water budgets largely represent the water budgets originally accreted during formation.

In comparison to the expected water mass fraction of material outside the water snow line of 50% (Lodders 2003), our maximum predicted water budgets are significantly smaller (30%). This can be due to a mixed accretion of volatile-poor and volatile-rich planetesimals from within and beyond the ice line, respectively. Sufficient mixing of volatile-poor and volatile-rich planetesimals for all planets is also reasonable given that no clear trend of water budget with orbital distance is seen as shown in Figure 13. Here we plot the inferred water mass fractions against the equilibrium temperature at zero albedo, which relates to orbital distance. The error bars on water mass fraction represent 3th–16th–84th–97th percentiles of the cumulative distribution. Neither a clear trend of water mass fraction with orbital distance nor one with planet mass could be identified. Thus, our findings suggest that the accreted material stems from regions in the disk where mixing was efficient enough to blur any possible trends of increasing water mass fraction with orbital distance, unless migration reordered the planets. Below, we discuss the implications for the formation of the planets in more detail.

The two main conclusions of our calculations are that the water content is very diverse among the planets, with no obvious correlation with either the mass or the period, and that the maximum amount of water is somewhat lower than the icy content of planetesimals usually assumed in the solar system.

The scatter in planetary water budgets is at odds with what we observe in the solar system, both in the planets and in Galilean satellites. In both cases, innermost objects tend to be dryer than outermost ones. While a uniform water budget in all planets cannot be excluded at the 2σ level (see Figure 13), this scatter is puzzling. Different origins can be invoked for this scatter. First, planets likely migrated during their formation (Alibert & Benz 2016), and therefore accreted planetesimals from different parts in the disk, some dry (inside the ice line) and some wet (beyond the ice line). In addition, the ice line itself is likely to move during the formation of planets, and this can modify the water content of solids in the disk (pebbles, but also small planetesimals). In both cases, one can therefore expect that the average water content of accreted solids (which does not necessary reflect the final water content of planets as discussed below) can vary depending on the exact formation path of the planets. What is, however, not clear is whether it is possible to end up with a water budget that shows no trend. Indeed, one would expect, from a smooth and ordered migration of planets, that the water content should monotonically increase as a function of period, as we see, for example, in the case of Galilean satellites. Breaking such a monotonic trend requires exchanging the place of some planets, which requires a strong level of dynamical rearrangement resulting from gravitational interactions. Whether such a rearrangement is possible, while at the same time ending with a flat and resonant system, remains to be investigated.

Another possibility is that the final water budget of planets is not the one of the accreted solids. Indeed, during their formation, planets accrete a gas envelope (which may be lost afterward). The interaction of the gas envelope, the protoplanetary disk, and the accreted solids can lead, under some circumstances, to the loss of some of the accreted solid material (Alibert 2017); however, see Brouwers et al. (2018). This process results from the advection of gas from the disk to the planetary atmosphere, and such a process could prevent accretion of water by forming planets. As the efficiency of advection depends on different parameters (including planetary mass, gas density in the disk, and distance to the star) in a nontrivial way, it is not clear whether such a process can destroy a preexisting water content gradient, but it is likely to increase the scatter in the final water content.

The second conclusion of our study is that the maximum water content is of the order to 30%. Interestingly enough, this value is smaller than the assumed water content in solids beyond the ice line, but at the same time similar to the maximum water content in solar system bodies up to the terrestrial mass range (see Figure 14). In fact, Ormel et al. (2017) predict moderate water budgets (∼10%) for the TRAPPIST-1 planets when planetary embryos form at the water snow line, migrate to the inner disk, and grow by dry pebbles until a critical planet mass (∼0.7 ${M}_{\oplus }$) is reached that prohibits further pebble accretion. A reduction of water budget for planets that accreted outside of the snow line can include processes such as desiccation of volatile-rich planetesimals by short-lived radionuclides (Grimm & McSween 1993; Lichtenberg et al. 2018), giant impacts between embryos (Genda & Abe 2005), and heating of planets and planetesimals during accretion and collisions, which is expected to be more efficient around low-mass stars (Lissauer 2007).

In any case, the water budget of the TRAPPIST-1 planets is an important constraint that needs to be fully considered in formation models. As a consequence, better determinations (e.g., by improving the mass, radius, or compositional constraints on refractory elements or atmospheric composition) of this quantity are key to making future progress in the understanding of the formation of the system.