Mass–Radius Relationships for Irradiated Ocean Planets

Our model iteratively solves the equations describing the interior of a planet:

$\begin{eqnarray}&&\displaystyle \frac{{dg}}{{dr}}=4\pi G\rho -\displaystyle \frac{2{Gm}}{{r}^{3}},\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{dr}}=-\rho g,\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&\displaystyle \frac{{dT}}{{dr}}=-g\gamma T\displaystyle \frac{d\rho }{{dP}},\end{eqnarray} \tag{ 3 }$

$\begin{eqnarray}&&P=f(\rho ,T),\end{eqnarray} \tag{ 4 }$

which correspond to the Gauss theorem, hydrostatic equilibrium, adiabatic profile with use of the Adams–Williamson equation, and the EOS of the considered medium, respectively. g, P, T, and ρ are gravity, pressure, temperature, and density profiles, respectively. m is the mass encapsulated within the radius r, G is the gravitational constant, and γ is the Grüneisen parameter. The Grüneisen parameter is key to computing the thermal profile of the planet, and the literature sometimes refers to the adiabatic gradient instead, expressed as follows (Kippenhahn et al. 2012; Mazevet et al. 2019; Haldemann et al. 2020):

$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mathrm{ad}}={\left(\displaystyle \frac{\partial \mathrm{ln}T}{\partial \mathrm{ln}P}\right)}_{S}=\gamma \displaystyle \frac{P}{\rho }\displaystyle \frac{1}{{c}^{2}},\end{eqnarray} \tag{ 5 }$

where S is the entropy, and c the speed of sound.

The interior model can display up to five distinct layers, depending on the planet's characteristics:

• 1.  
  a core made of metallic Fe and FeS alloy;
• 2.  
  a lower mantle made of bridgmanite and periclase;
• 3.  
  an upper mantle made of olivine and enstatite;
• 4.  
  an ice VII phase;
• 5.  
  a hydrosphere covering the whole fluid region of H₂O.

The Vinet EOS (Vinet et al. 1989) with thermal Debye correction is used for all solid phases:

$\begin{eqnarray}&&\begin{array}{l}P\left(\rho ,T\right)=3{K}_{0}\left[{\left(\displaystyle \frac{\rho }{{\rho }_{0}}\right)}^{\tfrac{2}{3}}-{\left(\displaystyle \frac{\rho }{{\rho }_{0}}\right)}^{\tfrac{1}{3}}\right]\\ \,\times \,\exp \left\{\displaystyle \frac{3}{2}\left({K}_{0}^{{\prime} }-1\right)\left[1-{\left(\displaystyle \frac{\rho }{{\rho }_{0}}\right)}^{-\tfrac{1}{3}}\right]\right\}\\ \,+{\rm{\Delta }}P,\end{array}\end{eqnarray} \tag{ 6 }$

with

$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}P=9\displaystyle \frac{\gamma \rho R}{{M}_{\mathrm{mol}}{\theta }^{3}}\\ \,\times \left[{T}^{4}{\displaystyle \int }_{0}^{\tfrac{\theta }{T}}\displaystyle \frac{{t}^{3}}{{e}^{t}-1}{dt}-{T}_{0}^{4}{\displaystyle \int }_{0}^{\tfrac{\theta }{{T}_{0}}}\displaystyle \frac{{t}^{3}}{{e}^{t}-1}{dt}\right],\end{array}\end{eqnarray} \tag{ 7 }$

where $\theta ={\theta }_{0}{\left(\tfrac{\rho }{{\rho }_{0}}\right)}^{\gamma }$, $\gamma ={\gamma }_{0}{\left(\tfrac{\rho }{{\rho }_{0}}\right)}^{-q}$, R is the ideal gas constant, and M_mol is the molar mass of the considered material. All quantities with an index 0 are reference parameters obtained by a fit on experimental data, given in Table 1. The EOS used by Mousis et al. (2020) to solve Equation (4) is the one formulated by Duan & Zhang (2006), valid up to 10 GPa and 2573.15 K.

Table 1. List of Thermodynamic and Compositional Parameters used in the Interior Model

Layer		Core		Lower Mantle				Upper Mantle
Phases		Iron-rich phase		Perovskite		Periclase		Olivine		Enstatite
Composition (%)		100		79.5		20.5		41		59
Components		Fe	FeS	FeSiO3	MgSiO3	FeO	MgO	Fe2SiO4	Mg2SiO4	Fe2Si2O6	Mg2Si2O6
Composition (%)		87	13	10	90	10	90	10	90	10	90
Molar mass (g mol−1)	Mmol	55.8457	87.9117	131.9294	100.3887	71.8451	40.3044	203.7745	140.6931	263.8588	200.7774
Reference density (kg m−3)	${\rho }_{0}$	8340	4900	5178	4108	5864	3584	4404	3222	4014	3215
Reference temperature (K)	T0	300		300		300		300		300
Reference bulk modulus (GPa)	K0	135		254.7		157		128		105.8
Pressure derivative of bulk modulus	${K}_{0}^{{\prime} }$	6		4.3		4		4.3		8.5
Reference Debye temperature (K)	${\theta }_{0}$	474		736		936		757		710
Reference Grüneisen parameter	${\gamma }_{0}$	1.36		2.23		1.45		1.11		1.009
Adiabatic power exponent	q	0.91		1.83		3		0.54		1

Note. Thermodynamic data are summarized in Table 1 of Sotin et al. (2007), and compositional data are the results of model calibration for Earth by Stacey (2005) and Sotin et al. (2007).

Download table as:  ASCIITypeset image

All thermodynamic and compositional parameters of mineral layers are taken equal to those of Earth (Stacey 2005; Sotin et al. 2007, 2010), and summarized in Table 1. We refer the reader to Brugger et al. (2017) to get all the computational details.

Mousis et al. (2020) computed the Grüneisen parameter for water via a bilinear interpolation in a grid generated from the python library of the IAPWS formulation, ³ computing the Grüneisen parameter in the form $\gamma =f(\rho ,T)$ with ρ and T varying in the 316–2500 kg m⁻³ and 650–10,000 K ranges, respectively. An important issue is that the density range is very limited, since this quantity can easily vary from ∼10 kg m⁻³ at the planetary surface to ∼5000 kg m⁻³ at the center of a 100% water planet of 1 M_⊕, implying that the computation of γ is erroneous at the top and at the bottom of the hydrosphere. A solution for overcoming this limitation is provided in Section 3.2.

Apart from compositional inputs, the main physical inputs of the model are the core mass fraction (CMF) x_core and WMF ${x}_{{{\rm{H}}}_{2}{\rm{O}}}$; the mantle mass fraction is then ${x}_{\mathrm{mantle}}\,=1-{x}_{\mathrm{core}}-{x}_{{{\rm{H}}}_{2}{\rm{O}}}$. Pressure and temperature profiles are integrated from outside, and require the inputs of the boundary pressure P_b and boundary temperature T_b. Finally, the model also requires the input of the planet's mass M_b (subscript b denotes the mass encapsulated within the boundary of the interior model, excluding the contribution of any potential atmosphere). Once defined, these input parameters allow for the computation of the planet's internal structure and associated boundary radius. In the case of the Earth (${x}_{\mathrm{core}}=0.325$, ${x}_{{{\rm{H}}}_{2}{\rm{O}}}=0.0005$, ${M}_{{\rm{b}}}=1$ ${M}_{\oplus }$), the model computes a radius R_b equal to 0.992 ${R}_{\oplus }$, which is less than 1% of error, indicating that errors from the model are negligible compared to errors on measurements. In the following, subscript b refers to quantities at the boundary between the interior model and the atmospheric model, such as bulk mass M_b, radius R_b, gravity g_b, pressure P_b and temperature T_b.

The atmospheric model generates the properties of a 1D spherical atmosphere of H₂O by integrating the thermodynamic profiles from bottom to top. The model takes as inputs the planet's mass and radius, as well as the thermodynamic conditions at its bottom. We choose to connect the atmospheric model with the interior model at a pressure P = ${P}_{{\rm{b}}}=300\,\mathrm{bar}$ (slightly above the critical pressure ${P}_{\mathrm{crit}}=220.67$ bar) and at a temperature T = T_b. $(P,T,\rho )$ profiles are then integrated upward via the prescription from Kasting (1988) in the case of an adiabat at hydrostatic equilibrium. Once the temperature reaches the top temperature of the atmospheric layer, here set to ${T}_{\mathrm{top}}=200$ K, an isothermal radiative mesosphere at $T={T}_{\mathrm{top}}$ is assumed. Figure 1 shows several (P,T) profiles representing the whole hydrospheres of planets (under Mazevet et al. (2019, hereafter Ma19+) parameterization, see Section 3.2) for masses and irradiation temperatures in the range 1–20 ${M}_{\oplus }$ and ${T}_{\mathrm{irr}}\,=$ 400–1300 K (see Equation (10)), respectively.

Zoom In Zoom Out Reset image size

The atmosphere transition radius is controlled by the altitude of the top of the H₂O clouds, corresponding to the top of the moist convective layer, assumed to be at a pressure ${P}_{\mathrm{top}}=0.1\,\mathrm{Pa}$. We choose this limit as the observable transiting radius, assuming that results are similar for cloudy and cloud-free atmospheres (Turbet et al. 2019, 2020). The EOS is taken from the NBS/NRC steam tables (Haar et al. 1984), implying the atmosphere is not treated as an ideal gas. The discontinuities in (P, T) profiles occurring for ${T}_{\mathrm{irr}}=1300$ K are due to the limited range of these tables, but the height of this region (P = 100–300 bar) is negligible compared to the thickness of the atmosphere. Increasing the T_top temperature will impact the final structure of the atmosphere, decreasing both the thickness of the atmosphere and the interior. Numerical tests with T_top varying from 200 K to ${T}_{\mathrm{skin}}={T}_{\mathrm{eff}}/{2}^{0.25}$ decrease the final radius of the planet by, at most, ∼200 km for the cases considered in this study. It corresponds to a difference of 2% in radius at most, but this difference is mainly below 1%.

Shortwave and thermal fluxes are then computed using four-stream approximation. Gaseous (line and continuum) absorptions are computed using the k-correlated method on 38 spectral bands in the thermal infrared, and 36 in the visible domain. Absorption coefficients are exactly the same as those in Leconte et al. (2013) and Turbet et al. (2019), which includes several databases specifically designed for H₂O-dominated atmospheres. Rayleigh opacity is also included. This method computes the total outgoing long-wave radiation (OLR, in W m⁻²) of the planet that gives the temperature that the planet would have if it was a blackbody:

$\begin{eqnarray}&&{T}_{{\rm{p}}}={\left(\displaystyle \frac{\mathrm{OLR}}{{\sigma }_{\mathrm{sb}}}\right)}^{1/4},\end{eqnarray} \tag{ 8 }$

with ${\sigma }_{\mathrm{sb}}$ being the Stefan–Boltzmann constant. In order to quantify the irradiation of the planet by its host star, we define the irradiance temperature

$\begin{eqnarray}&&{T}_{\mathrm{irr}}={T}_{\mathrm{eff}}\sqrt{\displaystyle \frac{{R}_{\star }}{2a}},\end{eqnarray} \tag{ 9 }$

where T_eff and ${R}_{\star }$ are the host star's effective temperature and radius, respectively, and a is the semimajor axis of the planet. The atmospheric model computes the Bond albedo from the atmosphere's reflectance (Pluriel et al. 2019) assuming a G-type star linking both temperatures:

$\begin{eqnarray}&&{T}_{\mathrm{irr}}={\left(\displaystyle \frac{\mathrm{OLR}}{(1-A){\sigma }_{\mathrm{sb}}}\right)}^{1/4}=\displaystyle \frac{{T}_{{\rm{p}}}}{{\left(1-A\right)}^{1/4}}.\end{eqnarray} \tag{ 10 }$

The literature often approximates T_irr to the equilibrium temperature T_eq, which is the temperature the planet would have for an albedo A = 0 (all the incoming heat is absorbed and re-emitted by the planet). Since it is the observable quantity, our results will be presented in term of T_irr. Equation (10) assumes that the planet is in radiative equilibrium with its host star. Any heating source in the planet's interior would add an additional term in the radiative equilibrium of the planet with its host star, increasing the effective temperature of the planet for the same received irradiation (Nettelmann et al. 2011). In this work, we model the structure of planets that have either no interior heating source, or that had time to cool off.

For a given planetary mass, boundary radius and irradiation temperature, the atmospheric thickness and boundary temperature are retrieved from the atmospheric model. The latter is then used to compute the interior structure, and the former is taken into account to compute the total (transiting) radius.

The choice of the EOS is critical, as it strongly impacts the estimate of the mass–radius relationships. Three EOS are then considered in this study:

• 1.  
  EOS from the latest revision of the IAPWS-95 formulation from Wagner & Pruß (2002) ⁴ (hereafter WP02). This reference EOS gives an analytical expression of the specific Helmholtz free energy $f(\rho ,T)$. Any thermodynamic quantity (pressure, heat capacity, internal energy, entropy, etc.) can be computed by taking the right derivative of f, and those quantities have analytical expressions.
• 2.  
  EOS from Duan & Zhang (2006) (hereafter DZ06). This EOS is corrected around the critical point, and gives an analytical expression for pressure as function of density and temperature $P(\rho ,T)$.
• 3.  
  EOS from Ma19. This formulation was developed for planetary interiors by extending the IAPWS-95 EOS with ingredients from statistical physics allowing a transition to plasma and superionic states. The authors created a Fortran implementation ⁵ that computes pressure, specific Helmholtz free energy, specific internal energy and specific heat capacity for a given couple $(\rho ,T)$.

The validity ranges of the different EOSs, which rely on the availability of experimental data, are given in Table 2. Extended ranges proposed by WP02 and DZ06 are also indicated because the mathematical expressions of their EOSs allow for extrapolations beyond the corresponding validity ranges. However, they become invalid when phase transition occurs (e.g., dissociation of water). Other EOSs exist in the literature, covering various regions of the phase diagram of water, or are being used for specific purposes. Our choice of EOSs among others is discussed in Section 6.

Table 2. Validity Ranges of the Different EOSs

EOS	Valid	Extended
WP02	$P\lt 1$ GPa	$P\lt 100$ GPa
	$T\lt 1273$ K	$T\lt 5000$ K
DZ06	$P\lt 10$ GPa	$P\lt 35$ GPa a
	$T\lt 2573.15$ K	$T\lt 2800$ K a
Ma19	$\rho \lt 100\times {10}^{3}\,\mathrm{kg}$ m−3	not specified
	$T\lt {\rm{100,000}}$ K	not specified

Note.

^a Limit given by Duan et al. (1996).

Download table as:  ASCIITypeset image

Figure 2 shows the $P(\rho )$ profiles derived from the considered EOSs at different temperatures. All EOSs present minor differences in their validity range, regardless of the considered temperature. Ma19 find that WP02 overestimates the pressure beyond its extended range. For a given pressure in a planet's interior, this would underestimate the density, and then overestimate the total radius of the planet. A more pronounced deviation is visible for DZ06 above its validity range. Around the critical point ($\rho \sim 350\,\mathrm{kg}$ m⁻³, mostly visible at 650 K), WP02 is closer to DZ06, compared to Ma19, as expected. In the low-density limit, all EOSs behave following the ideal gas law $P\propto \rho T$, which has a characteristic slope of 1 in a log–log scale.

Zoom In Zoom Out Reset image size

The Grüneisen parameter γ, already introduced in Equation (3), has many definitions. For solids, it gives the rate of change in phonon frequencies ${\omega }_{i}$ relative to a change in volume V (Grüneisen 1912):

$\begin{eqnarray}&&{\gamma }_{i}=-{\left(\displaystyle \frac{\partial \mathrm{ln}{\omega }_{i}}{\partial \mathrm{ln}V}\right)}_{T}.\end{eqnarray} \tag{ 12 }$

By averaging over all lattice frequencies, it is possible to obtain a thermodynamic definition (using the internal energy U and entropy S) of the Grüneisen parameter (Arp et al. 1984) via the following expression:

$\begin{eqnarray}&&\gamma =V{\left(\displaystyle \frac{\partial P}{\partial U}\right)}_{V}=\displaystyle \frac{V}{{C}_{V}}{\left(\displaystyle \frac{\partial P}{\partial T}\right)}_{V}=\displaystyle \frac{\rho }{T}{\left(\displaystyle \frac{\partial T}{\partial \rho }\right)}_{S}.\end{eqnarray} \tag{ 13 }$

γ relates a pressure (or density) variation to a temperature change. Although initially defined for solids, the meaning of γ holds for fluids. In planetary interiors, adiabatic heat exchange is mostly driven by convective heat transfer (Stacey & Hodgkinson 2019). At planetary scales, the Grüneisen parameter can thus be used for both solids and fluids. From identities in Equation (13), γ can be expressed using other thermodynamic constants such as the thermal expansion coefficient α, the isothermal bulk modulus K_T , and the specific isochoric heat capacity c_V :

$\begin{eqnarray}&&\gamma =\displaystyle \frac{\alpha {K}_{T}}{\rho {c}_{V}}.\end{eqnarray} \tag{ 14 }$

γ is assumed to be temperature independent in the solid phase, and its value is fitted from experimental data, taking into account small density variations. In this study, we use the Helmholtz free energy F given in Wagner & Pruß (WP02) and Ma19. In the IAPWS-95 release, the specific Helmholtz free energy f in its dimensionless form ϕ is divided into its ideal part (superscript ◦) and a residual (superscript r) via the following expression:

$\begin{eqnarray}&&\displaystyle \frac{f(\rho ,T)}{{RT}}=\phi (\delta ,\tau )=\phi ^\circ (\delta ,\tau )+{\phi }^{{\rm{r}}}(\delta ,\tau ),\end{eqnarray} \tag{ 15 }$

with $\delta =\rho /{\rho }_{c}$ and $\tau ={T}_{c}/T$, ${\rho }_{c}$ and T_c being the supercritical density and temperature, respectively. After defining the derivatives of the ideal and residual parts:

$\begin{eqnarray}&&{\phi }_{{mn}}^\circ =\displaystyle \frac{{\partial }^{m+n}\phi ^\circ (\tau ,\delta )}{\partial {\tau }^{m}\partial {\delta }^{n}},\end{eqnarray} \tag{ 16 }$

$\begin{eqnarray}&&{\phi }_{{mn}}^{{\rm{r}}}=\displaystyle \frac{{\partial }^{m+n}{\phi }^{{\rm{r}}}(\tau ,\delta )}{\partial {\tau }^{m}\partial {\delta }^{n}},\end{eqnarray} \tag{ 17 }$

where integers m and n define the order of the derivative with respect to τ and δ, respectively. From these expressions, one can derive:

$\begin{eqnarray}&&{\gamma }_{-}=-\displaystyle \frac{1+\delta {\phi }_{01}^{{\rm{r}}}-\delta \tau {\phi }_{11}^{{\rm{r}}}}{{\tau }^{2}\left({\phi }_{20}^{^\circ }+{\phi }_{20}^{{\rm{r}}}\right)},\end{eqnarray} \tag{ 18 }$

where ${\gamma }_{-}$ is the formulation of the Grüneisen parameter computed following the approach of WP02.

The fortran implementation of Ma19 computes $F(\rho ,T)$, along with other useful quantities such as ${\chi }_{T}={\left(\tfrac{\partial \mathrm{ln}P}{\partial \mathrm{ln}T}\right)}_{V}$ and the specific isochoric heat capacity c_V . In this case, the Grüneisen parameter is expressed as

$\begin{eqnarray}&&{\gamma }_{+}=\displaystyle \frac{P(\rho ,T){\chi }_{T}(\rho ,T)}{\rho {c}_{V}T},\end{eqnarray} \tag{ 19 }$

where ${\gamma }_{+}$ is the formulation of the Grüneisen parameter derived from the quantities calculated via the approach of Ma19.

In the case of an ideal gas, one can derive the theoretical value $\gamma =\tfrac{2}{l}$, where l is the number of degrees of freedom for a given molecule. For H₂O, $\gamma \simeq \tfrac{1}{3}$, since l = 6 (3 rotational and 3 vibrational degrees of freedom).

The Grüneisen parameter is crucial for computing the adiabatic temperature gradient inside a planet's interior. However, because temperature a has low impact on EOSs used in the solid phase, it is possible to assume isothermal layers in interior models when thermodynamic data are lacking, and generate internal structures close to reality (Zeng et al. 2019). In the case of fluids (here H₂O), temperature rises sharply with depth. This strongly impacts the EOS and leads to different phase changes that are not visible in the case of isothermal profiles.

Each computation for the interior model can be performed by using any of the three EOSs (WP02, DZ06, Ma19) to solve Equation (4), and the WP02 or Ma19 EOS to solve Equation (3) (i.e., computing γ with the EOS from WP02 or Ma19). In the following, we will use the name of the EOS used to solve Equation (4), and add + or − depending on the EOS used to compute the Grüneisen parameter, ${\gamma }_{+}$ (Ma19) or ${\gamma }_{-}$ (WP02) respectively. For example, Ma19- indicates that the Ma19 EOS was used to solve Equation (4), and that the WP02 approach was used to solve Equation (3).

Figure 3 shows the values of ${\gamma }_{+}$ and ${\gamma }_{-}$ in the H₂O phase diagram. Since γ is integrated to obtain the temperature gradient, a small difference leads to different paths in the (P,T) plane. The indiscernability between the WP02- and the Ma19- profiles shows that the internal structure (and thus mass–radius relationships) is more impacted by the temperature profile than the difference in the EOS in the case of a water layer. The difference in temperature between the Ma19+ and Ma19-/WP02- profiles is as high as ∼2000 K, which also results in a difference of ∼200 kg m⁻³ in density at the center of the planet.

Zoom In Zoom Out Reset image size

Atmospheric properties (OLR, albedo, mass and thickness) are all quantities that evolve smoothly. To enable a smooth connection between the two models, we implemented a trilinear interpolation module that can estimate atmospheric properties for a planet whose physical parameters g_b, M_b, and T_b are in the 3–30 m s⁻², 0.2–20 ${M}_{\oplus }$, and 750–4500 K ranges, respectively. This allows us to correct the slight deviations from nodes of the grid, and trilinear interpolation ensures that properties computed at a node are exactly those at the node, which would not be the case if a polynomial fit was performed on data. Details of the connection between the two models are given in Appendix A.

Figure 4 shows T_irr as a function of T_p for a set of fixed g_b and M_b. Due to a strong greenhouse (or blanketing) effect from the steam atmosphere, most cases lead to ${T}_{{\rm{b}}}\gt 2000$ K. As previously stated, this consequence discards any EOS that does not hold for such high temperatures. A second observation is that at low temperatures, one input an irradiation temperature T_irr that can correspond to two different planet temperatures T_p (and atmospheric properties). Since our work focuses on highly irradiated exoplanets, we will only investigate cases with ${T}_{\mathrm{irr}}\,\gt \,400$ K to bypass this degeneracy.

Zoom In Zoom Out Reset image size

One widely known process of atmospheric escape is the Jeans escape. Gas molecules have a velocity distribution given by the Maxwell–Boltzmann distribution, which displays an infinite extension in the velocity space, meaning that some particles have velocities greater than the escape velocity. By integrating this distribution, one can derive the Jeans particle flux (particles per time unit per surface unit) escaping the atmosphere at the exobase (Jeans 1925):

$\begin{eqnarray}&&{{\rm{\Phi }}}_{J}=\displaystyle \frac{{n}_{{\rm{e}}}{v}_{\mathrm{esc}}}{2\sqrt{\pi }}\displaystyle \frac{1}{\sqrt{\lambda }}\left(1+\lambda \right){e}^{-\lambda },\end{eqnarray} \tag{ 20 }$

where n_e is the particle number density at the top of the atmosphere (exobase) and ${v}_{\mathrm{esc}}=\sqrt{2{g}_{{\rm{b}}}{R}_{{\rm{p}}}}$ is the escape velocity (we assume ${R}_{{\rm{p}}}\simeq {R}_{{\rm{b}}}$ and ${M}_{{\rm{p}}}\simeq {M}_{{\rm{b}}}$). $\lambda ={\left(\tfrac{{v}_{\mathrm{esc}}}{{v}_{\mathrm{th}}}\right)}^{2}$ is the escape parameter, with ${v}_{\mathrm{th}}=\sqrt{2{R}_{g}{T}_{{\rm{e}}}/\mu }$ the average thermal velocity of molecules of mean molar mass μ at the exobase temperature T_exo, and R_g is the ideal gas constant.

We wish to provide an estimate of the physical characteristics of the planets that would lose more than a fraction ${x}_{\mathrm{lost}}=0.1$ of water content over a typical timescale of ${\rm{\Delta }}t=1\,\mathrm{Gyr}$. This condition is met when

$\begin{eqnarray}&&4\pi {R}_{{\rm{p}}}^{2}\displaystyle \frac{\mu }{{{ \mathcal N }}_{{\rm{A}}}}{{\rm{\Phi }}}_{J}\geqslant \displaystyle \frac{{x}_{\mathrm{lost}}{M}_{{\rm{p}}}}{{\rm{\Delta }}t},\end{eqnarray} \tag{ 21 }$

with ${{ \mathcal N }}_{{\rm{A}}}$ being the Avogadro number. Solving Equation (21) with Earth's properties (${n}_{{\rm{e}}}=\tfrac{{P}_{\mathrm{top}}{{ \mathcal N }}_{{\rm{A}}}}{{{RgT}}_{\mathrm{exo}}}\sim {10}^{19}$, ${R}_{{\rm{p}}}={R}_{\oplus }$, ${M}_{{\rm{p}}}\,={M}_{\oplus }$) yields $\lambda \leqslant 100$. Due to the exponential term, the result is not very sensitive to changes in parameters, including the exact location of the exobase. Assuming ${T}_{\mathrm{exo}}={T}_{\mathrm{irr}}$, this condition can be rewritten as

$\begin{eqnarray}&&{R}_{{\rm{p}}}\gt \displaystyle \frac{1}{\lambda }\displaystyle \frac{G\mu }{{R}_{g}{T}_{\mathrm{irr}}}{M}_{{\rm{p}}},\end{eqnarray} \tag{ 22 }$

with G being the gravitational constant. This estimate is consistent with today's composition of planets of the solar system (see Figure 8). Equation (22) gives an indication of the properties of the planets that are subject to H₂ or H₂O escape, implying that their atmospheres should be dominated by heavier molecules (H₂O, CO₂, O₂, CH₄, etc.) or be rocky planets, respectively.

Hydrodynamic escape, also referred to as hydrodynamic blowoff, occurs when upper layers of the atmosphere are heated by intercepting the high-energy irradiation (far UV, extreme UV and X-ray fluxes, the sum of which is often called the XUV flux) from the host star. This heating induces an upward flow of gas, leading to mass loss at a rate (Erkaev et al. 2007; Owen & Wu 2013)

$\begin{eqnarray}&&\dot{M}=\epsilon \displaystyle \frac{{L}_{\mathrm{XUV}}{R}_{{\rm{p}}}^{3}}{{{GM}}_{{\rm{p}}}{\left(2a\right)}^{2}},\end{eqnarray} \tag{ 23 }$

where L_XUV is the host star XUV luminosity, a is the planet's orbital distance and is a conversion factor between incident irradiation energy and mechanical blowoff energy. Note that Equation (23) is only true in the energy-limited case. Heating occurs by absorption of high-energy photons by molecules which are dissociated in the upper atmosphere, meaning that blowoff can be limited by (i) the number of photons as one photon breaks one molecule, and (ii) the recombination time as a dissociated molecule may recombine before being able to absorb the XUV irradiation again. Boundaries between these regimes have been explored by Owen & Alvarez (2016), who showed that the sub-Neptune population undergoes mostly energy-limited mass loss, validating the use of Equation (23) in our case.

For our estimate, we use the X-ray and UV luminosities obtained by fits on observational data for M to F type stars by Sanz-Forcada et al. (2011):

$\begin{eqnarray}&&{L}_{\mathrm{EUV}}={10}^{3.8}{L}_{{\rm{X}}}^{0.86},\end{eqnarray} \tag{ 24 }$

$\begin{eqnarray}\begin{array}{rcl}{L}_{{\rm{X}}} & = & 6.3\times {10}^{-4}{L}_{\star },\quad \tau \lt {\tau }_{\mathrm{sat}}\\ & = & 1.89\times {10}^{21}{\tau }^{-1.55},\quad \tau \gt {\tau }_{\mathrm{sat}}\end{array}\end{eqnarray} \tag{ 25 }$

where τ is the host star age in gigayears and ${\tau }_{\mathrm{sat}}\,=5.72\times {10}^{15}{L}_{\star }^{-0.65}$ (Sanz-Forcada et al. 2011). To estimate the XUV luminosity, the star's bolometric luminosity is assumed constant, a hypothesis supported by the stellar evolution tracks of Baraffe et al. (2015). Integrating the XUV luminosity in the saturation regime ($0\lt \tau \lt {\tau }_{\mathrm{sat}}$) and beyond gives the finite quantity ${E}_{\mathrm{XUV}}={\int }_{0}^{+\infty }{L}_{\mathrm{XUV}}\,{dt}=1.8\times {10}^{39}$ W for a solar type star.

Again, we look for planets that could lose more than 10% of their mass over a 1 Gyr period due to atmospheric blowoff:

$\begin{eqnarray}&&\epsilon \displaystyle \frac{{E}_{\mathrm{XUV}}{R}_{{\rm{p}}}^{3}}{{{GM}}_{{\rm{p}}}{\left(2a\right)}^{2}}\geqslant {x}_{\mathrm{lost}}{M}_{{\rm{p}}}.\end{eqnarray} \tag{ 26 }$

Combining Equation (9) and Stefan–Boltzmann's law ${L}_{\star }\,=\,4\pi {R}_{\star }^{2}{\sigma }_{\mathrm{sb}}{T}_{\mathrm{eff}}^{4}$ gives

$\begin{eqnarray}&&{\left(2a\right)}^{2}=\displaystyle \frac{1}{{T}_{\mathrm{irr}}^{4}}\displaystyle \frac{{L}_{\star }}{4\pi {\sigma }_{\mathrm{sb}}}.\end{eqnarray} \tag{ 27 }$

Substituting this expression in Equation (26) yields the condition:

$\begin{eqnarray}&&{R}_{{\rm{p}}}\geqslant {M}_{{\rm{p}}}^{\tfrac{2}{3}}{\left(\displaystyle \frac{{x}_{\mathrm{lost}}G}{\epsilon 4\pi {\sigma }_{\mathrm{sb}}{T}_{\mathrm{irr}}^{4}{E}_{\mathrm{XUV}}}\right)}^{\tfrac{1}{3}}.\end{eqnarray} \tag{ 28 }$

This condition only gives an indication of the planets that are subject to substantial hydrodynamic escape. All arbitrary quantities such as $\epsilon \simeq 0.1$ (Owen & Jackson 2012; Bolmont et al. 2017) and E_XUV are affected by a power of 1/3, resulting in a low dependency on the chosen values.

The nature of escaping particles is not considered in Equation (28), meaning the computed quantity is the total lost mass. Bolmont et al. (2017) developed a method to quantify the hydrodynamic outflow r_F (how many atoms of oxygen leave for each hydrogen atom). Based on their work, we compute ${r}_{{\rm{F}}}\sim 0.2$, indicating substantial loss of both H and O, with an accumulation of O₂. Mass loss of water content and accumulation of O₂ have several implications for the habitability of exoplanets (Ribas et al. 2016; Schaefer et al. 2016). The power laws for mass is 1 and $\tfrac{2}{3}$ in the cases of Equations (22) and (28), respectively. This implies that hydrodynamic escape is more efficient for less dense planets. In contrast, Jeans escape is dominant in the case of denser planets. The power law for T_irr is −1 and $-\tfrac{4}{3}$ in the cases of Equations (22) and (28), respectively, implying that hydrodynamic escape will take over Jeans escape at higher irradiation temperatures. As shown in Figure 5, Equations (22) and (28) leave a window for planets that lost their H₂ reservoir but kept heavier volatiles from which they formed (Zeng et al. 2019). This result highlights the consistency between the possible existence of irradiated ocean planets and atmospheric escape.

Zoom In Zoom Out Reset image size

Figure 5 presents computed mass–radius relationships for the ${\gamma }_{+}$ parameterization and assumes Earthlike properties for the rocky part (see Table 1). As predicted from the shape of EOSs curves, WP02 and DZ06 EOSs underestimate the density and thus produce larger planets. This effect is accentuated for more-massive planets with a larger amount of water, corresponding to cases where water pressure reaches the highest values. The radius is also overestimated for low-mass planets, because the hydrosphere becomes extended due to the low gravity, implying that a slight underestimation of the density can still lead to a substantial difference in radius. These results show the incontestable asset of the EOS developed by Ma19, and rule out the possibility of using the WP02 or DZ06 EOSs to produce reliable mass–radius relationships for planets with substantial amounts of water. To remain in the true validity ranges of the WP02 or DZ06 EOSs, one should consider a few percent of water content at most in the planet.

As discussed in Section 3.2, ${\gamma }_{+}$ is always lower than ${\gamma }_{-}$ in Earth-sized planets fully made of water. As a result, (P,T) profiles for ${\gamma }_{+}$ parameterizations are steeper than for ${\gamma }_{-}$ parameterizations (see Figure 3), meaning the interior is colder for ${\gamma }_{+}$. In turn, colder planets will be denser and thus smaller. The impact of the choice between ${\gamma }_{+}$ and ${\gamma }_{-}$ is shown in Figure 6, where the relative difference on the radius between Ma19+ and Ma19- parameterizations is presented. In all cases, the relative difference between the models is 10% at most.

Zoom In Zoom Out Reset image size

As the mass of a planet increases, its gravity becomes more important, and its hydrosphere (interior structure and atmosphere) consequently thinner. Thinner hydrospheres, especially in the case of massive planets, lead to smaller relative differences in radii. Moreover, values of ${\gamma }_{+}$ and ${\gamma }_{-}$ become closer (and even equal) in the 10¹–10² GPa pressure range (see Figure 3), thus reducing the radii differences between the Ma19+ and Ma19- parameterizations even more significantly.

The value of γ increases when the (P,T) curves of a hydrosphere approach the liquid–ice VII transition, which leads to a more important temperature gradient that prevents the formation of high-pressure ices. This observation is in major disagreement with models assuming isothermal hydrospheres (Valencia et al. 2006; Seager et al. 2007; Valencia et al. 2007; Zeng & Sasselov 2013; Brugger et al. 2017; Zeng et al. 2019), a hypothesis often justified by assuming that temperature has a secondary impact on EOSs, which remains a valid statement for solid phases but not in the case of the hydrosphere. A correct treatment of the temperature gradient (Mousis et al. 2020) leads to the presence of high-temperature phases for H₂O (ionic, superionic, plasma), which are more dilated, significantly impacting the mass–radius relationships.

In the following, we use ${\gamma }_{+}$ and Ma19 to compute the mass–radius relationships. Indeed, the pressure and temperature ranges in the hydrospheres of sub-Neptune-like planets lie well in the region for which the Ma19 formulation was developed. Also, due to the blanketing effect of the atmosphere, even the coldest planets irradiated at ${T}_{\mathrm{irr}}\,=\,400$ K have a temperature of more than 2000 K at the 300 bar interface (see Figure 4), which corresponds to the pressure at which the atmospheric and internal models are connected. This interface is already located well above the range of validity of ${\gamma }_{-}$.

Mass–radius relationships only provide an estimate of a possible exoplanet's composition. A more precise assessment is achieved via the use of compositional ternary diagrams. For a given planet's mass and irradiation temperature, such a diagram shows the radius as a function of the planet's WMF and CMF. Possible compositions as thus retrieved from the contour at the level of the planet's measured radius. Computations presented here use only the central value of the mass of each planet, thus not taking into account the measurement error on the planet's mass.

Possible compositions of the three planets of the GJ 9827 system are shown in Figure 7, based on the planets' parameter measurements made by Rice et al. (2019). Planet b exhibits an Earthlike interior without the need to invoke a significant steam atmosphere. The presence of a thick steam atmosphere is quite consistent with the low-density measurements made for planet c, with a water content ranging from 1% to 8%. Physical properties (mass, radius, and temperature) of planet c lead to important Jeans escape (with our criterion in Equation (21), see Figure 8), suggesting the absence of H₂ and He in the atmosphere. Moreover, planet c is unlikely to accrete a substantial amount of H₂ and He due to its low mass. Although planet d is consistent with a Jupiter-like interior due to their similar bulk densities, again, its high irradiation temperature suggests the presence of a H₂–He free atmosphere. Isochrones used by Rice et al. (2019) fix a lower limit of 5 Gyr on the age of the system, which makes an H₂–He atmosphere less likely as Jeans escape would remove them. Applying our model to the current measurements yields a WMF in the 5%–30% range for planet d. These results are summarized in Table 3.

Zoom In Zoom Out Reset image size

Table 3. Planetary Parameters of the GJ 9827 System Used as Inputs for the Model, and Estimated WMFs using Ternary Diagrams (Figure 7)

Planet	b	c	d
${M}_{{\rm{p}}}\,({M}_{\oplus })$	4.91 ± 0.49	0.84 ± 0.66	4.04 ± 0.83
${R}_{{\rm{p}}}\,({R}_{\oplus })$	1.58 ± 0.03	1.24 ± 0.03	2.02 ± 0.05
Tirr (K)	1184	820	686
WMF (%)	0–5	1–5	5–30

Download table as:  ASCIITypeset image

Ternary diagrams presented here do not take into account the uncertainty on each planet's mass, and were computed for the central value only. If a planet's mass is slightly higher (respectively lower), its density increases (respectively decreases), while the estimated WMF diminishes (respectively grows). This implies that the mass and radius of a planet must be measured with extreme accuracy to constrain the WMF properly. Additional constraints can be applied from observational data such as the stellar elemental ratios (Fe/Si, Mg/Si) that could help in constraining the core to mantle mass ratio (Brugger et al. 2017), and methods such as Markov Chain Monte Carlo can be performed to simultaneously determine all parameters (Acuña et al. 2021).

Figure 8 represents the computed mass–radius relationships for WMFs of 0.2, 0.5, and 1. In this figure, the condition for substantial atmospheric loss due to Jeans escape is derived by solving Equation (21) for each planet. One already known effect is that steam atmospheres are very extended (Mousis et al. 2020), allowing the computation of compositions without invoking small H₂–He envelopes (1%–5% by mass). The second effect is heating due to the adiabatic gradient, which decreases the density, and then increases the radius. In the 10–20 ${M}_{\oplus }$ range, the radius of a planet with a WMF of 50% made of liquid H₂O is equal to that of a planet with a WMF of 20% constituted of supercritical H₂O. Also, the radius of a planet fully made of liquid H₂O is equivalent to that of a planet with half of its mass constituted of supercritical H₂O. This shows how important the error on the computation of WMF can be, depending on the physical assumptions made. In the figures presented in Mousis et al. (2020), where the DZ06 EOS was used, the model was able to match Neptune's mass (17 ${M}_{\oplus }$) and radius (3.88 ${R}_{\oplus }$) with a 95% H₂O interior at 300 K. With the Ma19 EOS, a 100% water planet presents a radius of 3.25 ${R}_{\oplus }$ at ${T}_{\mathrm{irr}}=400$ K and 3.6 ${R}_{\oplus }$ at ${T}_{\mathrm{irr}}=1300$ K.

Zoom In Zoom Out Reset image size

All produced mass–radius relationships are very well approximated by an equation of the form

$\begin{eqnarray}&&\mathrm{log}{R}_{{\rm{p}}}=a\mathrm{log}{M}_{{\rm{p}}}+b+\exp \left(-d(\mathrm{log}{M}_{{\rm{p}}}+c)\right),\end{eqnarray} \tag{ 29 }$

where log denotes the decimal logarithm, and R_p and M_p are normalized to Earth units. a, b, c, and d are coefficients obtained by fits, and have one value for each composition $({x}_{\mathrm{core}},{x}_{{{\rm{H}}}_{2}{\rm{O}}})$ and each temperature T_irr. For each fitted curve, we define the mean absolute error between the data and fit as

$\begin{eqnarray}&&\mathrm{MAE}=\displaystyle \frac{1}{N}\sum _{i=1}^{N}\left|\displaystyle \frac{{R}_{{\rm{p}},\mathrm{model}}-{R}_{{\rm{p}},\mathrm{fit}}}{{R}_{{\rm{p}},\mathrm{model}}}\right|.\end{eqnarray} \tag{ 30 }$

Values of the MAE are 0.01%–1% for all fits, indicating a good accuracy. The largest deviation between one point $({M}_{{\rm{p}}},{R}_{{\rm{p}}})$ and the fitted curve is of 2.3%, meaning the deviation between data and fit can be neglected. Fitted coefficients vary smoothly with respect to the three parameters $({x}_{\mathrm{core}},{x}_{{{\rm{H}}}_{2}{\rm{O}}},{T}_{\mathrm{irr}})$, allowing a good interpolation of the intermediate values. The produced grid uses the compositional parameters for the core and mantle calibrated for Earth (see Table 1), and data may be different if Fe/Si or Mg/Si ratios are different.