11institutetext: Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117
Heidelberg, Germany
11email<EMAIL_ADDRESS>22institutetext: Aix-Marseille Université, CNRS, CNES,
Institut Origines, LAM, Marseille, France 33institutetext: Institut
universitaire de France (IUF), France

# Interior-atmosphere modelling to assess the observability of rocky planets
with JWST

L. Acuña 1122 M. Deleuil 2233 O. Mousis 2233

(Received 20 December 2022; accepted 24 April 2023)

###### Abstract

Context. Super-Earths present compositions dominated by refractory materials.
However, there is a degeneracy in their interior structure between a planet
with no atmosphere and a small Fe content, and a planet with a thin atmosphere
and a higher core mass fraction. To break this degeneracy, atmospheric
characterization observations are required.

Aims. We present a self-consistent interior-atmosphere model to constrain the
volatile mass fraction, surface pressure, and temperature of rocky planets
with water and CO2 atmospheres. These parameters obtained in our analysis can
then be used to predict observations in emission spectroscopy and photometry
with JWST, which can determine the presence of an atmosphere, and if present,
its composition.

Methods. We couple a 1D interior model with a supercritical water layer with
an atmospheric model. To obtain the bolometric emission and Bond albedo for an
atmosphere in radiative-convective equilibrium, we present the k-uncorrelated
approximation for fast computations within our retrieval on planetary mass,
radius and host stellar abundances. For the generation of emission spectra, we
use our k-correlated atmospheric model. An adaptive Markov chain Monte Carlo
(MCMC) is used for an efficient sampling of the parameter space at low
volatile mass fractions.

Results. We show how to use our modelling approach to predict observations
with JWST for TRAPPIST-1 c and 55 Cancri e, which have been proposed in Cycle
1. TRAPPIST-1 c’s most likely scenario is a bare surface, although the
presence of an atmosphere cannot be ruled out. If the emission in the MIRI
F1500 filter is 731 ppm or higher, there would be a water-rich atmosphere. For
fluxes between 730 and 400 ppm, no atmosphere is present, while low emission
fluxes (300 ppm) indicate a CO2-dominated atmosphere. In the case of 55 Cancri
e, a combined spectrum with NIRCam and MIRI LRS may present high uncertainties
at wavelengths between 3 and 3.7 $\mu$m. However, this does not affect the
identification of H2O and CO2 because they do not present spectral features in
this wavelength range.

###### Key Words.:

Planets and satellites: atmospheres – Planets and satellites: interiors –
Planets and satellites: composition – Planets and satellites: individual:
TRAPPIST-1 d – Planets and satellites: individual: 55 Cnc e – Methods:
statistical – Methods: numerical

##  1 Introduction

Low-mass exoplanets ($M<20\ M_{\oplus}$) have two different sub-populations
based on their radius and density: super-Earths and sub-Neptunes. Super-Earths
have radii of $R=1.3\ R_{\oplus}$, while the radius of sub-Neptunes
corresponds to $R=2.4\ R_{\oplus}$ (Fulton et al., 2017; Fulton & Petigura,
2018). If we compare these radii with planet interior and evolution models,
super-Earths are mostly composed of Fe and Si-bearing rocks, whereas sub-
Neptunes have a significant volatile (H/He, water) content. Despite having an
idea of the main component for these planets from their mass and radius data,
and interior structure models, we do not know their exact interior composition
due to degeneracies.

In the case of super-Earths, we still have the question of whether the planet
could have a thin atmosphere, or if it presents a bare rock surface.
Atmospheres containing H/He are discarded, since a very small fraction of H/He
entails a minimum radius of $\simeq$ 1.6 $R_{\oplus}$ (Lopez & Fortney, 2014).
Therefore, an atmosphere composed of water formed from ice pebbles accreted
beyond or in the vicinity of the water ice line (Mousis et al., 2019;
Krissansen-Totton et al., 2021; Kimura & Ikoma, 2022), a secondary atmosphere
built-up by outgassing (Ortenzi et al., 2020; Baumeister et al., 2021; Liggins
et al., 2022), or a silicate atmosphere (Zahnle et al., 2009) are the most
likely scenarios for super-Earths. This produces a degeneracy in the internal
structure of super-Earths and Earth-sized planets: the same planetary mass and
radius can be explained with a planet with no atmosphere and a low-Fe content
rocky bulk (Madhusudhan, 2012; Dorn et al., 2017), or a planet with a thin
atmosphere and a core mass fraction similar to that of Earth (32% of core mass
fraction).

This degeneracy in interior structure can only be broken with the support of
atmospheric characterization data. The presence of an atmosphere has been
confirmed in the hot super-Earth $\pi$ Mensae c, whose detection of C II ions
indicates atmospheric escape of a high-molecular atmosphere (García Muñoz et
al., 2021). Phase curves have also been used to determine the existence of a
silicate atmosphere in K2-141 b (Zieba et al., 2022), as well as transmission
spectroscopy for the terrestrial planet LHS 3844 b (Diamond-Lowe et al.,
2020). Moreover, Kreidberg et al. (2019) use the phase curves to confirm the
absence of an atmosphere, as well as to constrain which material constitutes
the planetary surface.

JWST (Gardner et al., 2006) will observe several super-Earths to confirm the
presence of an atmosphere, or even narrow their possible atmospheric
compositions. In this study, we present a self-consistent interior-atmosphere
model, Marseille’s Super-Earth Interior model (MSEI), to perform retrievals
from estimated mass, radius, and stellar host abundances. As a result, we
determine the posterior distribution functions (PDFs) of the atmospheric mass,
and surface pressure and temperature of water and CO2-dominated atmospheres.
These atmospheric parameters obtained from our retrieval analysis can be used
as input for an atmospheric model to produce spectra. We set an example of
this with our k-correlated atmospheric model, MSEIRADTRAN, to generate
emission spectra to predict observations with JWST with MIRI photometric
filters, and NIRCam and MIRI Low-Resolution Spectrometer (MIRI LRS). Combined
interior and atmospheric models have been used to constrain the water mass
fractions (WMF) of rocky planets as WMF ¡ $10^{-3}$ (Agol et al., 2021).
However, the consistent exploration of the parameter space in the region close
to WMF = 0 is necessary to obtain the PDFs of the water mass fraction and the
surface pressure accurately. For this reason, we employ an adaptive Markov
chain Monte Carlo (MCMC) (Director et al., 2017) to explore the low surface
pressures for possible water and CO2 atmospheres in rocky planets.

We recall the basics of our interior model, MSEI, in Sect. 2. We explain the
updates implemented in our atmosphere model, MSEIRADTRAN, with respect to
similar previous k-correlated models (Marcq et al., 2017; Pluriel et al.,
2019) in Sect. 3, as well as introduce the k-uncorrelated approximation for
fast computations of the bolometric flux and Bond albedo within our MCMC
framework. In Sect. 4, we detail the implementation of the adaptive MCMC, and
show and example of the retrieval with it and a non-adaptive one. We assess
the observability of two planets that have been proposed for observations in
Cycle 1 of JWST with our model: TRAPPIST-1 c (Gillon et al., 2016; Grimm et
al., 2018) and 55 Cancri e (Ehrenreich et al., 2012; Bourrier et al., 2018).
The planet and instrument parameters that we use as input for our interior-
atmospheric analyses and Pandexo (Batalha et al., 2020) to predict
uncertainties in JWST observations are summarised in Sect. 5. We present our
results and conclusions in Sect. 6 and Sect. 8, respectively.

##  2 Interior structure model

We remind the reader the fundamental principles on which our interior
structure model is based. The input of the interior structure model are the
total mass and two compositional parameters, which are the core mass fraction
(CMF), and the water mass fraction (WMF). The CMF and the WMF are defined as
the mass of the Fe-rich core and the hysdrosphere divided by the total
planetary mass, respectively. In the 1D interior model, the planetary radius,
$r$, is represented by a one-dimensional grid. Along this grid, the pressure,
$P(r)$, the temperature, $T(r)$, the gravity acceleration, $g(r)$, and the
density, $\rho(r)$, are calculated at each point. These four variables are
obtained by solving the equation that defines it. The pressure is computed by
integrating the equation of hydrostatic equilibrium (see Eq. 1), while the
temperature requires integrating the adiabatic gradient profile (Eq. 2). In
low-mass planets, the opacity in the deep interior is large enough for the
radiative temperature gradient to be higher than the adiabatic gradient,
making the layers unstable against convection according to the Schwarzschild
criterion. In Eq. 2, $\gamma$ and $\phi$ correspond to the Grüneisen and
seismic parameters, respectively. The former describes the behaviour of the
temperature in a crystal to its density. The latter parameter provides the
speed at which seismic waves propagate in the same crystalline structure.
Their formal definitions are shown in Eq. 3, where we can see that the seismic
parameter is related to the slope of the density at constant pressure, while
the Grüneisen parameter depends on the derivative of the pressure with respect
to the internal energy, $E$. The acceleration of gravity is obtained by
solving the integral that results from Gauss’s theorem (Eq. 4), where $G$ is
the gravitational constant, and $m$ corresponds to the mass at a given radius,
$r$.

$\dfrac{dP}{dr}=-\rho g$ (1)

$\dfrac{dT}{dr}=-g\dfrac{\gamma T}{\phi}$ (2)

$\begin{cases}\phi=\dfrac{dP}{d\rho}\\\ \gamma=V\
\left(\dfrac{dP}{dE}\right)_{V}\end{cases}$ (3)

$\dfrac{dg}{dr}=4\pi G\rho-\dfrac{2Gm}{r^{3}}$ (4)

The density, $\rho(r)$, is computed with the Equation of State (EOS), which
provides the density as a function of temperature and pressure. The interior
structure model is divided three separate layers: a Fe-rich core, a mantle
rich in silicates, and a water layer. We use a different EOS to calculate the
density for each of these layers. We adopt the Vinet EOS (Vinet et al., 1989),
with a thermal correction, for the core and the mantle. More details about
this EOS and its reference parameter values for the core and mantle can be
found in Brugger et al. (2016, 2017). For the hydrosphere, we use the EOS and
specific internal energy of Mazevet et al. (2019) for supercritical and plasma
phases of water, which is valid within the pressure and temperature regime
($P>300$ bar, $T>700$ K) covered by our interior structure model. We discuss
the validity ranges of different water EOS for this regime in Acuña et al.
(2021), while a detailed comparison of different EOS for high-pressure and
high-temperature water and its effects on the total radius of the planet can
be found in Aguichine et al. (2021).

The final input for our interior structure model are the surface temperature
and pressure. These are the boundary conditions, together with the
gravitational acceleration at the centre of the planet, whose value is zero,
$g(r=0)$ = 0. Finally, the mass of each planetary layer is obtained by
integrating the equation of conservation of mass (Eq. 5). The total planetary
mass is the sum of the individual mass of the layers. When the total input
mass and the initial boundary conditions are met, the model reaches
convergence.

$\dfrac{dm}{dr}=4\pi r^{2}\rho$ (5)

###  2.1 Interior-atmosphere coupling

The surface pressure for the interior model depends on the mass of the
atmospheric mass on top of the outermost interface of the interior model. For
envelopes whose bottom pressure is greater or equal to $P$ = 300 bar, the
interior model’s surface pressure is set constant to 300 bar, which is the
interface at which the interior and the atmosphere are coupled. Then the
supercritical water layer extends from this interface to the boundary between
the hydrosphere and the silicate mantle at higher pressures. For atmospheres
whose surface pressure is less than 300 bar, the interior and atmosphere are
coupled at the atmosphere-mantle interface, having the water envelope in
vapour phase only. The water mass fraction takes into account the mass of the
atmosphere, $M_{atm}$. The atmospheric mass is calculated as in Eq. 6, where
$P_{base}$ is the pressure at the base of the atmosphere (atmosphere-interior
interface), $R_{bulk}$ is the radius comprised between the centre of the
planet and the base of the atmosphere, and $g_{surf}$ is the acceleration of
gravity at this interface. The coupling interface between the interior and the
atmosphere models at a maximum pressure of 300 bar is sufficiently close to
the critical point ($P$ = 220 bar) of water to prevent the atmospheric model
from taking over pressures at which convection dominates over radiation. The
EOS we use for the interior (Mazevet et al., 2019) and the atmosphere
(Haldemann et al., 2020) are based on the IAPWS-95 EOS. Wagner & Pruß (2002)
report that the IAPWS-95 EOS presents unsatisfactory features in a small
pressure and temperature range around the critical point. To prevent
discontinuities in the adiabat and the density between the interior and the
atmosphere for planets whose adiabat passes through this area, we set the
coupling interface at 300 bar, and not at $P_{crit}$ = 220 bar.

$M_{atm}=\dfrac{P_{base}4\pi R_{bulk}^{2}}{g_{surf}}$ (6)

The atmospheric model calculates the outgoing longwave radiation (OLR) and the
Bond albedo, $A_{B}$, given as a function of bulk mass, radius and temperature
at the bottom of the atmosphere. If an atmosphere is in radiative equilibrium,
its absorbed flux, $F_{abs}$ must be equal to its emitted radiation, which is
the OLR. The absorbed flux depends on the Bond albedo via Eqs. 7 and 8, where
$\sigma$ is the Stefan-Boltzmann constant, and $T_{eq}$ is the planetary
equilibrium temperature. This requires to know the semi-major axis of the
planet, $a_{d}$, as well as the stellar radius and effective temperature,
$R_{\star}$ and $T_{\star}$, respectively.

$F_{abs}=\sigma\ T_{\mathrm{eq}}^{4}.$ (7)

$T_{\mathrm{eq}}=(1-A_{B})^{0.25}\left(0.5\frac{R_{\star}}{a_{d}}\right)^{0.5}T_{\star},$
(8)

For a constant planetary mass and radius, we find the temperature at the base
of the atmosphere by solving OLR$(T_{base})-F_{abs}(T_{base})=0$ with a root-
finding method, such as the bisection method. Then this is the input boundary
condition for the interior structure model. The radius calculated by the
interior structure model (from the centre of the planet up to the base of the
atmosphere) is an input for the atmospheric model, while the temperature at
the bottom of the atmosphere is an input for both the interior and the
atmospheric model. Therefore, the self-consistent coupling of both models is
not straightforward, and requires an iterative algorithm that checks that
convergence is reached for the total radius and surface temperature. The total
radius are computed as the sum of the bulk radius calculated by the interior
model, and the atmospheric thickness obtained by the atmospheric model. We
refer the reader to Acuña et al. (2021) for a detailed description of this
algorithm.

##  3 Atmospheric model

The interior-atmosphere coupling presented in our previous work (Mousis et
al., 2020; Acuña et al., 2021, 2022) was done by using grids of data generated
by the atmospheric model of Pluriel et al. (2019). These grids provide the
OLR, Bond albedo and atmospheric thickness for a given set of mass, radius and
surface temperature assuming a constant surface pressure. However, the use of
these grids presents the following disadvantage: the grids do not enable us to
generate emission spectra that could be used to simulate observations.
Therefore, we develop our own atmospheric model, MSEIRADTRAN. We started the
development of MSEIRADTRAN by modifying the atmosphere model presented in
Marcq et al. (2017)111http://marcq.page.latmos.ipsl.fr/radconv1d.html to
include up-to-date opacity and EOS data. In the following, we recall the basic
structure and principles of MSEIRADTRAN and the atmospheric models presented
in Marcq et al. (2017); Pluriel et al. (2019).

We consider two scenarios for the composition of the envelope: water-dominated
envelopes (99% water plus 1% CO2) and CO2-dominated envelopes (99% CO2 plus 1%
water). Including a wider variety of relative mass fractions between water and
CO2 in our models only increases the degeneracies between atmospheric mass and
atmospheric composition. Therefore, we only consider the two end-members to
assess the observability of water and CO2 spectral features with JWST. We do
not model 100% pure water or CO2 atmospheres because these pure compositions
are very unlikely due to outgassing and atmospheric escape (Krissansen-Totton
et al., 2021). In addition, to make our comparison between the uncorrelated-k
approximation and the k-correlated atmospheric model of Pluriel et al. (2019)
consistent (see Section 3.3), we use the same exact compositions of 99%:1%
instead of 100% pure water or CO2.

The 1D atmospheric model first proposes a pressure-temperature (PT) profile.
This profile consists of a near-surface, dry convective layer, followed by a
wet convective region where condensation takes place, and a isothermal
mesosphere on top. If the surface temperature is cold enough to allow for
condensation of water, the dry troposphere will not exist. For the isothermal
mesosphere, we assume a constant temperature of 200 K (Marcq, 2012; Marcq et
al., 2017). The OLR is not very dependent on the temperature of an upper
mesosphere (Kasting, 1988). In addition, we do not take into account
mesospheric stellar heating, which could increase significantly the
temperature of the mesosphere. Therefore, adopting a low mesospheric
temperature yields similar thermal profiles to self-consistent atmospheric
calculations (Lupu et al., 2014). The one-dimensional grid that represents the
pressure contains 512 computational layers. The adiabatic gradient used to
calculate the temperature in each of these points in the convective regions
depends on whether it is located in the dry or wet convective layer. The
details of the computation of the wet and dry adiabatic gradients are recalled
in Sect. 3.1.

The calculations of the emission spectrum and the Bond albedo are performed by
bands. We divide the spectrum from 0 to 10100 cm-1 (equivalent to $\geq 1\
\mu$m in wavelength) in 36 bands to obtain the OLR, similar to Pluriel et al.
(2019). For each band, we calculate the total optical depth in each
computational layer, which has four different contributions. These
contributions are the optical depth due to collision-induced absorption (CIA),
and line opacity (see Sect. 3.4), Rayleigh scattering, and clouds. We treat
Rayleigh scattering as in Pluriel et al. (2019), where the Rayleigh scattering
opacity is related to wavelength, $\lambda$, following Eq. 9. The parameters
$\kappa_{0}$ and $\lambda_{0}$ are adopted from Kopparapu et al. (2013) and
Sneep & Ubachs (2005) for H2O, and CO2, respectively. The opacity of clouds is
considered for the atmospheric layers where condensation takes place. Similar
to Marcq et al. (2017) and Pluriel et al. (2019), the cloud opacity is
parameterised after Kasting (1988), who assumes a cloud opacity proportional
to the extinction coefficient, $Q_{ext}$ (see Eq. 10). The dependence of the
extinction coefficient on wavelength (Eq. 11) is similar to that of water
clouds on Earth (Kasting, 1988; Marcq et al., 2017).

$\kappa_{Rayleigh}(\lambda)=\kappa_{0}\left(\dfrac{\lambda_{0}}{\lambda}\right)^{4}$
(9)

$\kappa_{clouds}(\lambda)=130\ Q_{ext}(\lambda)$ (10)

$Q_{ext}=\begin{cases}1&\lambda\leq 20\ \mu m\\\
3.26\cdot\lambda^{-0.4}&\lambda>20\ \mu m\end{cases}$ (11)

The total optical depth (Eq. 12), together with the PT profile, are the input
for the radiative transfer equation solver, DISORT (Stamnes et al., 2017).
DISORT obtains the emitted upward flux at the top of the atmosphere (TOA). The
TOA flux is calculated for all 36 bands, which are then summed to obtain the
bolometric, wavelength-integrated TOA flux, or OLR. Immediately after the OLR
is computed, we start the calculation of reflection spectra and the Bond
albedo, described in Sect. 3.5.

$\tau_{total}=\tau_{line}+\tau_{Rayleigh}+\tau_{CIA}+\tau_{clouds}$ (12)

###  3.1 Atmospheric EOS

The thermal structure of the atmosphere is divided in two main layers. The
near-surface layer is adiabatic, which means convection takes place, while the
layer on top is a mesosphere with constant temperature. This is set to 200 K,
which is representative of the cool temperatures that hot low-mass planets
present in their mesospheres (Lupu et al., 2014; Leconte et al., 2013). The
adiabatic layer is divided in two sub-layers: a dry convective layer and a wet
convective one. Condensation may take place depending on the value of the
pressure compared to the water saturation pressure. In the two scenarios we
consider for the composition of the envelope, 1) 99% water and 1% CO2 and 2)
1% water and 99% CO2, water is a condensable species, while CO2 is a non-
condensable gas. To obtain the temperature in a computational layer, $i$, we
consider two approximations. The first one is that the change in temperature
and pressure within an atmospheric layer is small enough to approximate
$\left(\dfrac{\partial T}{\partial
P}\right)_{S}\simeq\dfrac{T_{i-1}-T_{i}}{P_{i-1}-P_{i}}$. The second
approximation is $\Delta P\sim P_{i}\ \Delta ln(P)$, because $\dfrac{d\
ln(P)}{dP}\sim\dfrac{\Delta ln(P)}{\Delta P}=\dfrac{1}{P_{i}}$. Then the
temperature of a given computational layer, $T_{i}$ can be calculated as a
function of the pressure in that layer, $P_{i}$, and the derivative of the
temperature with respect to pressure at constant entropy, $(\partial T\partial
P)_{S}$, as seen in Eq. 13.

$T_{i}=T_{i-1}-\left(\dfrac{\partial T}{\partial P}\right)_{S}\ P_{i}\ \Delta
ln(P_{i})$ (13)

If the pressure of water vapour is below the gas saturation pressure,
$P_{v}<P_{sat}(T)$, or if its temperature is above the temperature of the
critical point of water, $T>T_{crit}$, we are under dry convection (Marcq,
2012; Marcq et al., 2017). We calculate the derivative $(dT/dP)_{S}$ in the
dry case following Marcq et al. (2017) (Eq. 14). $\rho_{v}$ and $\rho_{c}$ are
the densities of water vapour and CO2, respectively; and $C_{p,v}$ and
$C_{p,c}$, their heat capacities. $V_{v}=1/\rho_{v}$ is the specific volume of
water vapour.

$\left(\dfrac{\partial T}{\partial P}\right)_{S,\ dry}=\dfrac{\rho_{v}\ T\
(\partial V_{v}/\partial T)_{P}}{\rho_{v}\ C_{p,v}+\rho_{c}\ C_{p,c}}$ (14)

When the atmospheric pressure reaches the water saturation curve,
$P=P_{sat}(T)$, water vapour starts to condense out and clouds form. Since the
phase change requires energy in the form of latent heat, the wet adiabatic
coefficient is different from the dry adiabatic one (Eq. 14). The expression
for the derivative $(dT/dP)_{S}$ in the wet case with CO2 as the only non-
condensable gas is provided in Eq. 15) (Marcq et al., 2017). $M_{c}$ is the
molecular weight of carbon dioxide, and $C_{v,c}$ is the specific heat
capacity at constant volume of CO2. $R$ is the ideal gas constant, while
$\alpha_{v}$ is the mixing ratio of water vapour density relative to CO2,
$\alpha_{v}=\rho_{v}/\rho_{c}$. Its derivative, $\partial\
ln(\alpha_{v})/\partial\ ln(T)$ (Eq. 16), needs to be calculated as well
(Kasting, 1988; Marcq, 2012; Marcq et al., 2017).

$\left(\dfrac{\partial T}{\partial P}\right)_{S,\
wet}=\dfrac{1}{(dP_{sat}/dT)+\rho_{c}R/M_{c}(1+\partial\
ln(\rho_{v})/\partial\ ln(T)-\partial\ ln(\alpha_{v})/\partial\ ln(T))}$

(15)

$\dfrac{\partial\ ln(\alpha_{v})}{\partial\ ln(T)}=\dfrac{R/M_{c}(\partial\
ln(\rho_{v})/\partial\ ln(T))-C_{v,c}-\alpha_{v}(\partial s_{v}/\partial\
ln(T))}{\alpha_{v}(s_{v}-s_{c})+R/M_{c}}$ (16)

The density and heat capacity of water were previously obtained by using the
steam tables of Haar et al. (1984). These tables treat water as a non-ideal
gas, although they are not valid for $T>2500$ K. Therefore, for temperatures
higher than 2500 K, we use the EOS tables from Haldemann et al. (2020) to
calculate the thermodynamic properties of water. These tables are a
compilation of different EOSs, where each EOS is applied in its validity
region of the water phase diagram. There are two EOSs that are used in the
region relevant for the atmospheres of low-mass, highly-irradiated planets.
The first EOS is the IAPWS95 (Wagner & Pruß, 2002), whose validity range for
the high-pressure supercritical regime corresponds to 251 to 1273 K in
temperature, and up to 1 GPa in pressure. Haldemann et al. (2020) transition
to an EOS valid at low pressures and high temperatures at 1200 K. This second
EOS is the CEA (Chemical Equilibrium with Applications) package (Gordon, 1994;
McBride, 1996). This package incorporates the effects of single ionisation and
thermal dissociation, which are processes that occur only at high temperatures
in gas phase. Fig. 1 shows the dry adiabatic coefficient as a function of
pressure and temperature in the region of the water phase diagram relevant for
hot planetary atmospheres. The temperature derivative $(dT/dP)_{S}$ is closely
related to the adiabatic coefficient, $\kappa_{ad}$ (Eq. 17). The reduction of
the dry adiabatic coefficient at $T=$ 1000 to 2500 K is due to thermal
dissociation, whereas the decrease at higher temperatures ($T\geq$ 6000 K) is
caused by thermal ionisation (Haldemann et al., 2020).

$\left(\dfrac{\partial T}{\partial P}\right)_{S}=\dfrac{T}{P}\
\kappa_{ad}(P,T)$ (17)

Figure 1: Dry adiabatic coefficient for water, $\kappa_{ad}$, as a function of
pressure and temperature. The displayed region covers the cold and hot gas
phase of water up to pressures close to the critical point. Solid, dashed, and
dotted contours corresponds to $\kappa_{ad}$ = 0.35, 0.25 and 0.15,
respectively.

We assume that CO2 is an ideal gas to calculate its density, $\rho_{c}$. As we
treat CO2 as an ideal gas, we can calculate its heat capacity
$C_{v,c}=C_{p,c}-R/M_{c}$. We calculate the specific heat of CO2 by using
Shomate’s equation (Eq. 18), whose tabulated coefficients $A$ to $E$ are
provided by Chase (1998). This is part of the NIST (National Institute of
Standards and Technology) chemistry webbook
database222https://webbook.nist.gov/cgi/cbook.cgi?ID=C124389&Mask=1, and it is
valid for temperatures up to 6000 K.

$C_{p,c}(T)=A+B\ T+C\ T^{2}+D\ T^{3}+E/T^{2}$ (18)

Finally, we determine the atmospheric thickness, $z_{atm}$ under the
assumption of hydrostatic equilibrium (Eq. 1). In the case of the atmospheric
model, the pressure is integrated over altitude, $z$, instead of internal
radius, $r$. The altitude of the computational atmospheric layer $i$ is shown
in Eq. 19, where $\rho_{total,i}$ is the total mass density at altitude
$z_{i}$, $\rho_{total,i}=\rho_{v,i}+\rho_{c,i}$. This expression is derived by
approximating $(dP/dz)\sim\Delta P/z_{i-1}-z_{i}$. The gravity acceleration is
computed at each point of the 1D grid with the equation for hydrostatic
equilibrium. This is noted explicitly in Eq. 19 by labelling the gravity at
the computational layer $i$ as $g_{i}$. The final atmospheric thickness is
obtained by evaluating $z$ at the transit pressure, 20 mbar.

$z_{i}=z_{i-1}+\dfrac{P_{i}\ \Delta ln(P)}{g_{i}\ \rho_{total,i}}$ (19)

###  3.2 K-correlated method implementation

We employ the k-correlated method (Liou, 1980) to compute the emission
spectrum. We discretise the dependence of the opacity on wavenumber, $\nu$, by
dividing the spectrum in bands. The spectral transmittance of a spectral band
$b$, $\mathcal{T}_{b}$, is defined as the exponential sum of the opacity of
the band, $\kappa^{b}$, times the column density $m$, which is only dependent
on pressure and temperature (Sanchez-Lavega, 2011). The exponential sum is
performed over $G$-points, which are the abscissa values chosen to discretise
the cumulative probability function of the opacity, $G(\kappa)$. Each
$G$-point, $G_{i}$, has an associated weight in the exponential sum in Eq. 20,
$w_{i}$. The discrete opacity value is not only dependent on the spectral
band, but also on the $G$-point, which is noted explicitly in Eq. 20 as
$\kappa^{b}_{i}$. $N_{G}$ corresponds to the total number of $G$-points.

$\mathcal{T}^{b}(m)=\int_{0}^{1}e^{-\kappa^{\nu}(G)\ m}\
dG\simeq\sum_{i=1}^{N_{G}}e^{-\kappa^{b}_{i}\ m(P,T)}\ w_{i}$ (20)

$F^{\uparrow}_{TOA,\ b}=\int_{0}^{1}F^{b}(G)\
dG=\sum_{i=1}^{N_{G}}F^{b}_{G_{i}}\ w_{i}$ (21)

In the case of an arbitrary set of $G$-points, the weights are equal to
$w_{i}=\Delta G_{i}$, which are the widths of the bins in $G$-space. In each
atmospheric layer, the pressure and temperature are considered constant. Under
this condition and within each spectral band, we can exchange wavenumber with
$G$ (Mollière, 2017). Consequently, we integrate the upward flux over $G$ to
obtain the emitted flux within each band (Eq. 21). The upward top flux per bin
and per $G$-point, $F^{b}_{G_{i}}$, is obtained by invoking the radiative
solver for a total optical depth whose line optical depth is
$\tau_{line}=\tau^{b}_{i}$ (Eq. 22).

$\tau^{b}_{i}=\kappa^{b}_{i}\ \rho\ \Delta z$ (22)

Following Malik et al. (2017, 2019), we define our discrete $G$-points as the
roots of the Legendre polynomial, $G_{LG,i}$ (Eq. 23). The corresponding
weigths, $w_{i}$, are the Legendre-Gaussian (LG) weights associated to the
Legendre polynomial of $N_{G}$th order, $P_{N_{G}}$ (Eq. 24). The LG weights
are calculated from the $i$th root of the $N_{G}$th order Legendre polynomial,
$y_{i}$, as well as from the polynomial’s derivative, $P^{\prime}_{N_{G}}$
(Abramowitz & Stegun, 1972; Malik et al., 2017). We choose to have $N_{G}=$ 16
$G$-points.

$G_{LG,i}=\dfrac{(1+y_{i})}{2}$ (23)

$w_{LG,i}=\dfrac{2}{\left(1-y_{i}^{2}\right)P^{\prime}_{N_{G}}\left(y_{i}\right)^{2}}$
(24)

To generate emission spectra to assess the observability with JWST, we use the
original spectral resolution of the opacity k-table data (see Sect. 3.4) in
our k-correlated model, which corresponds to $R=200$ to 300 in the spectral
range $\lambda$ = 1 to 20 $\mu$m.

###  3.3 K-uncorrelated approximation

The k-correlated method requires solving the radiative transfer equation 16
times per band, with a different line optical depth (see Eq. 22). The
radiative transfer solver is the most computationally expensive subroutine in
our 1D atmospheric model, so we adopt a new approximation that we refer to as
the k-uncorrelated method, which require less calls to DISORT per spectral
bin. The difference between the uncorrelated-k and the correlated-k methods is
where the iteration over $G$-points occurs. In the k-correlated method, the
spectral lines are assumed to be correlated from an atmospheric layer to the
next one, making it necessary to propagate the fluxes through the atmosphere
for all $G$-points before performing the integral in Eq. 20. In the
uncorrelated-k approach, the iteration over $G$-points ends before inputting
the optical depth to the RT solver. The final line optical depth that is the
input for the RT solver in the k-uncorrelated method assumes that all the
contributions from the $G$-points are included before the sum in Eq. 20, being
calculated as $\tau^{b}=-ln(\mathcal{T}^{b})$. Therefore, the radiative solver
is invoked once per spectral bin. This reduces the computation time of the
emission flux in one spectral band from 0.40 seconds to 0.22 seconds. In
addition, we reduce the spectral resolution to $R=\lambda/\Delta\lambda=10$
for the emission spectrum and Bond albedo, which is fast enough to perform a
high number of interior-atmosphere models in a retrieval analysis (see Sect.
4) within a few days.

Figure 2: Comparison of the radius and interior-atmosphere boundary
temperature between a k-correlated model and the k-uncorrelated approximation.
Left panel: Mass-radius relationships for a planet with a water-dominated
atmosphere orbiting a Sun-like star at $a_{d}=0.05$ AU. Dashed lines indicate
the total radius calculated by the k-uncorrelated version of MSEIRADTRAN,
while the solid line corresponds to the interior radius, which comprises the
core, mantle, and supercritical water (SW). Triangles and circles indicate the
total radius and the interior radius obtained when the interior model is
coupled with the atmospheric model of Pluriel et al. (2019), respectively.
Right panel: Temperature at the 300 bar interface as a function of planetary
mass.

We use the version of our atmospheric model with the k-uncorrelated
approximation to generate the OLR and Bond albedo necessary to calculate the
boundary conditions for the interior model. To compare the effect of the
k-uncorrelated approximation on the total radius of planets with water-rich
envelopes, we compute two sets of mass-radius relationships (Fig. 2, upper
panel). The first set is obtained by coupling the interior model with our
k-uncorrelated model (dotted lines), and the second one with the k-correlated
model of Pluriel et al. (2019). The k-uncorrelated approximation causes
differences of less than 1% in radius compared to the k-correlated model of
Pluriel et al. (2019). This difference stems from the difference in the
temperature at the interior-atmosphere interface ($P$ = 300 bar) between
MSEIRADTRAN and Pluriel et al. (2019). The difference in temperature is
between 200 to 230 K (Fig. 2, lower panel). Note that the constant surface
temperature at masses below 5 $M_{\oplus}$ in the Pluriel et al. (2019) model
are due to the limit in the data grid we use, having no physical
interpretation. In Fig. 3 we show the OLR and Bond albedo as a function of
bottom temperature for our atmospheric model and that of Pluriel et al.
(2019). The trends of the Bond albedo are very similar for both models,
obtaining a similar absorbed flux (dashed lines). Consequently, the
discrepancy in bottom temperature is caused by the difference in OLR between
the two models, which is approximately a factor of 2.

We conclude that the k-uncorrelated MSEI may be used to calculate the radius
of water-rich planets within our MCMC retrieval framework, since fast
computations of the OLR and Bond albedo are necessary. The k-uncorrelated
approach is an alternative to grey models. We show a detailed comparison of
this approximation with a grey model and our k-correlated model in Appendix A
for the OLR and the Bond albedo. This approximation should not be used to
compute spectra at higher resolution since the planetary emission differs by a
factor of 2 compared to a k-correlated method. This entails a significant
difference in predictions of the emission spectrum for observations.

Figure 3: Outgoing Longwave Radiation (OLR, solid lines) and absorbed flux
(dashed lines) as a function of bottom atmospheric temperature, calculated
with our atmospheric model, MSEIRADTRAN (red), and Pluriel et al. (2019)
(blue). Dotted lines indicate the temperature at which the OLR and the
absorbed flux are equal, meaning the atmosphere is in radiative-convective
equilibrium.

###  3.4 Opacity data

The total optical depth computed in the atmospheric model include
contributions from collision-induced absorption (CIA) and line absorption. CIA
absorption is particularly important in dense gases, such as steam and CO2 at
high pressures, specially if the line opacity is weak (Pluriel et al., 2019).
For the atmospheric compositions we consider in this work, namely 99% H2O:1%
CO2, and 99% CO2:1% H2O, we require CIA data for collisions between these two
gases, and their self-induced absorptions. We adopt CIA absorption data for
H2O-CO2 and H2O-H2O collisions provided by Ma & Tipping (1992) and Tran et al.
(2018)333https://www.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/continuum_data/,
respectively. CO2-CO2 CIA opacities are read from a look-up table obtained by
Bézard et al. (2011) and Marcq et al. (2008), which is also used in the
atmospheric model by Marcq et al. (2017). Our H2O-H2O CIA table covers the
complete spectral range where we calculate both our emission and reflection
spectra, while the H2O-CO2 CIA table covers the bands with wavelength $\geq$ 1
$\mu$m, which corresponds to the bands necessary for the emission spectrum
only. For the bands whose wavelength is outside the spectral range of the CIA
table, we assume a constant CIA opacity value equal to the opacity at the
limit band of the table.

Grimm & Heng (2015) provide a
database444https://chaldene.unibe.ch/data/Opacity3/ of pre-calculated opacity
k-tables for different species and line lists. For water and CO2-dominated
atmospheres, we adopt the POKAZATEL (Polyansky et al., 2018) and HITEMP2010
(Rothman et al., 2010) opacity data, respectively. POKAZATEL presents the
widest validity range in temperature for water in planetary atmospheres, with
a maximum temperature of 5000 K, while HITEMP maximum temperature is 4000 K.
Following the procedure described in Leconte (2021), we bin the k-correlated
opacities to the same spectral bins of Marcq et al. (2017) and Pluriel et al.
(2019). We calculate the k-coefficients for our water-CO2 mixture by assuming
that the spectral features of the individual gases are correlated (Malik et
al., 2017). In the correlated approximation, the mixed opacity is estimated as
indicated in Eq. 25, where $\chi_{j}$ is the mixing ratio by mass of the $j$th
gas, $\kappa_{j,i}$ is the k-coefficient of the $j$th gas evaluated at the
$G_{i}$ point. The mixing ratio by mass is defined as $\chi_{j}=\dfrac{X_{j}\
MW_{j}}{\mu}$, where $X_{j}=\dfrac{P_{j}}{P}$ is the volume mixing ratio of
the $j$th species, $MW_{j}$ is its molecular weight, and $\mu$ is the mean
molecular weight of the mixture.

$\kappa_{mix,i}=\sum_{j=1}^{N_{gases}}\chi_{j}\ \kappa_{j,i}$ (25)

###  3.5 Reflection spectra and Bond albedo

Once the bolometric OLR is obtained, we initiate the calculation of the
reflectivity in 30 bands, from 5 to 0.29 $\mu$m, to obtain the Bond albedo
(Pluriel et al., 2019). The bands for which we calculate both the emission
flux and the reflectivity (from 1 to 5 $\mu$m) cannot have the two quantities
calculated simultaneously since DISORT requires different input settings to
calculate them. For the emission, we assume zero illumination from the top of
the atmosphere, as well as an upward flux that forms 90 degrees with the
surface of the planet, which corresponds to a polar angle equal to zero.
DISORT calculates the reflectivity of the atmosphere as a function of incident
beam angle, which corresponds to the solar zenith angle (SZA), $\theta$ in Eq.
26. This is the angle that the incident light forms with the normal of the
incident surface. Once we obtain the dependence of the reflectivity with SZA,
we can average it as indicated in Eq. 26 (Simonelli & Veverka, 1988). To
integrate Eq. 26, we evaluate the reflectivity at 10 different SZA values. We
assume 4 streams for DISORT in both the calculation of the OLR and the
reflectivity.

$A_{B}(\nu)=2\int_{0}^{\pi/2}A_{B}(\nu,\theta)\ cos(\theta)\ sin(\theta)\
d\theta$ (26)

After averaging the reflectivity over SZA, we obtain the reflection spectrum,
which is the dependence of the albedo as a function of wavenumber. To obtain
the bolometric Bond albedo, we integrate Eq. 27 (Pluriel et al., 2019).
$A_{B}(\nu)$ is the reflectivity as a function of wavenumber,
$B_{\nu}(T_{\star})$ is Planck’s function for a temperature equal to the
effective temperature of the host star, and $\sigma$ is the Stefan-Boltzmann
constant.

$A_{Bond,\ bol}=\dfrac{\pi\ \int_{0}^{\infty}A_{B}(\nu)\ B_{\nu}(T_{\star})\
d\nu}{\sigma T_{\star}^{4}}$ (27)

The Bond albedo is a parameter particularly sensitive to the choice of phase
function. For atmospheric layers that present clouds, the gas contributes to
scattering with a Rayleigh phase function, while we assume the Henyey-
Greenstein phase funtion for clouds. DISORT requires the calculation of the
Legendre moments of the combined phase function, which we estimate as the
weighted average of the moments of the two individual phase functions (Liou,
1980; Boucher, 1998). The weights are calculated as the ratio of the optical
depth due to Rayleigh scattering or clouds divided by the total optical depth,
$\tau_{Rayleigh}+\tau_{clouds}$, for Rayleigh and Henyey-Greenstein phase
function moments, respectively. For clear atmospheric layers, the only
contribution to scattering is Rayleigh scattering due to the gas, so the total
phase function moment corresponds to that of Rayleigh scattering (Marcq et
al., 2017).

DISORT also requires as input the single scattering albedo of each atmospheric
layer. The single scattering albedo is defined as the ratio of scattering
efficiency to total extinction effi- ciency. The total extinction is a sum of
both extinction by scattering and extinction by absorption. Therefore, a
single scattering albedo of 1 indicates that all extinction is due to
scattering, whereas a value of zero means that absorption dominates. Similar
to the moments of the phase function, we estimate the combined single
scattering albedo from gas (Rayleigh) and clouds with their weighted average.
The single scattering albedo due to Rayleigh scattering is calculated as the
Rayleigh optical depth divided by the total optical depth,
$\tau_{clear}+\tau_{Rayleigh}$. The clear optical depth is the sum of the line
and CIA optical depths (Sect. 3.4). The single scattering albedo due to clouds
is calculated as the ratio of the clouds optical depth divided by
$\tau_{clear}+\tau_{cloud}$, times the cloud single scattering albedo defined
in Kasting (1988):

$\varpi_{0}=\begin{cases}1&\lambda\leq 2\ \mu m\\\
1.24\cdot\lambda^{-0.32}&\lambda>2\ \mu m\\\ \end{cases}$ (28)

##  4 Markov chain Monte Carlo (MCMC)

We use the MCMC Bayesian algorithm described in Dorn et al. (2015) and later
adapted by Acuña et al. (2021) to our forward interior-atmosphere model. We
recall that the model parameters are the planetary mass, $M$, the CMF,
$x_{core}$, and the WMF, $x_{H2O}$. Therefore, one single model is determined
by these three parameters as m = $\left\\{M,x_{core},x_{H2O}\right\\}$. The
available data are the total mass $M$, the total radius $R$, and the Fe/Si
abundance, d = $\left\\{M_{obs},R_{obs},Fe/Si_{obs}\right\\}$, whose
observational errors are
$\sigma(M_{obs}),\sigma(R_{obs}),\sigma(Fe/Si_{obs})$, respectively. When the
Fe/Si mole ratio is not considered in the inverse problem, the data is reduced
to only the total planetary mass and radius, d =
$\left\\{M_{obs},R_{obs}\right\\}$. The prior information consists on a
Gaussian distribution centered on the mean value of the observed mass, with a
standard deviation equal to the observational uncertainty. For the CMF and
WMF, we consider as priors uniform distributions ranging from 0 to 1. The MCMC
scheme first starts by drawing a value for each of the model parameters from
their prior distributions, which we denote as
m${}_{1}=\left\\{M_{1},x_{core,1},x_{H2O,1}\right\\}$. The interior model
calculates the planetary radius and Fe/Si abundance that corresponds to these
model parameters, which is
g(m${}_{1})=\left\\{R_{1},M_{1},Fe/Si_{1}\right\\}$. The likelihood of this
model is computed (Eqs. 29 and 30) and a new model is drawn from the prior
distributions, $\textbf{m}_{2}$.

$L(\textbf{m}_{i}\mid\textbf{d})=C\ exp\
\Biggl{(}-\frac{1}{2}\Biggl{[}\left(\frac{(R_{i}-R_{obs})}{\sigma(R_{obs})}\right)^{2}+\left(\frac{(M_{i}-M_{obs})}{\sigma(M_{obs})}\right)^{2}\\\
+\left(\frac{(Fe/Si_{i}-Fe/Si_{obs})}{\sigma(Fe/Si_{obs})}\right)^{2}\Biggr{]}\Biggr{)},$
(29)

$C=\dfrac{1}{(2\pi)^{3/2}\left[\sigma^{2}(M_{obs})\cdot\sigma^{2}(R_{obs})\cdot\sigma^{2}(Fe/Si_{obs})\right]^{1/2}}$
(30)

The log-likelihoods,
$l(\textbf{m}_{i}\mid\textbf{d})=log(L(\textbf{m}_{i}\mid\textbf{d}))$, of
both models are used to estimate the acceptance probability (Eq. 31).
Consecutively, a random number is drawn from a uniform distribution between 0
and 1. If $P_{accept}$ is greater than this random number, m2 is accepted and
the chain moves to this set of model parameters, starting the following chain
$n+1$. Otherwise, the chain remains in m1 and a different set of model
parameters is proposed, m3. The accepted models are stored, and values of
their parameters conform the PDF that will enable us to estimate their mean
and uncertainties.

$P_{accept}=min\left\\{1,e^{(l(\textbf{m}_{new,i}\mid\textbf{d})-l(\textbf{m}_{old}\mid\textbf{d}))}\right\\}$
(31)

###  4.1 Adaptive MCMC

In our initial implementation of the MCMC (Acuña et al., 2021), the random
walker uses a uniform distribution to choose the next state where it is going
to move in the parameter space of the CMF and WMF. This is called a naive walk
(Mosegaard & Tarantola, 1995), in which all points in the parameter space have
a probability of being chosen proportional to their number of neighbours. This
poses the following problem: for the states whose CMF or WMF are close to 0 or
1, they are less likely to be sampled in the random walk, because they have
less neighbours than the central values. A model with WMF = 1 is not physical,
although many highly-irradiated rocky planets might present low-mass
atmospheres, which correspond to a WMF close to 0. To compensate the lower
probability of being chosen in the limiting states of the prior, we use an
adaptive step size in the walker. This consists on adapting the maximum size
of the perturbation used to generate a new model instead of using a fixed
value everywhere in the parameter space. This adaptive step size will decrease
in the limiting areas of the prior (i.e low WMF states) and have its greatest
value in the centre of the prior (WMF = 0.5). The self-adjusting step size is
carried by a transformation of the parameter space, which ranges from
exponential to spherical transforms (Chaudhry et al., 2021). In this work, we
choose to implement the self-adjusting logit transform (SALT), proposed by
Director et al. (2017). The SALT transform is publicly available in the
SALTSampler R package555https://rdrr.io/cran/SALTSampler/man/SALTSampler-
package.html, which eases its implementation in Python for our own model.

We compare the non-adaptive and adaptive MCMC for one planet, TOI-220 b (Hoyer
et al., 2021). We consider as input data the total mass and radius, as well as
the Fe/Si mole ratio, which has been calculated with the stellar abundances of
the host star. No maximum limits have been established for the CMF or the WMF.
TOI-220 b has an equilibrium temperature of 806 K, which means that it is
strongly irradiated and could present steam and supercritical phases. Table 1
presents the input data, and the retrieved parameters of the non-adaptive and
adaptive MCMCs. All three agree within uncertainties for mass, radius and
Fe/Si. The uncertainties of the mass and radius in the non-adaptive MCMC are
smaller than the input data. This difference in uncertainties is significant
in the case of the total mass. This indicates that the non-adaptive MCMC is
not as effective as the adaptive MCMC at sampling all possible
$\left\\{x_{core},x_{H2O}\right\\}$ pairs that could reproduce the mass and
radius data. As a consequence, the uncertainties of the WMF are being
underestimated in the non-adaptive MCMC, while the adaptive MCMC produces a
greater confidence interval for the WMF, and retrieves the exact uncertainties
of the mass and radius.

Table 1: TOI-220 b MCMC input (Data), and output mean values and 1$\sigma$ uncertainties for the non-adaptive and adaptive MCMCs. | Data (Hoyer et al., 2021) | Non-adaptive | Adaptive
---|---|---|---
$M$ [$M_{\oplus}$] | 13.8$\pm$1.0 | 13.8$\pm$0.7 | 13.7$\pm$1.0
$R$ [$R_{\oplus}$] | 3.03$\pm$0.15 | 3.06$\pm$0.12 | 2.98$\pm$0.15
Fe/Si | 0.65$\pm$0.09 | 0.64$\pm$0.11 | 0.64$\pm$0.10
$x_{core}$ |  | 0.08$\pm$0.03 | 0.09$\pm$0.03
$x_{H2O}$ |  | 0.62$\pm$0.10 | 0.58$\pm$0.14

In Fig. 4, we show the sampled 2D PDFs for the CMF and the WMF in the ternary
diagram. In addition to the same area of the ternary diagram as the non-
adaptive algorithm, the adaptive MCMC explores an area at lower WMF along the
Fe/Si = 0.65 isoline, going down to WMF = 0.10 in the driest simulations. This
is a consequence of the ability of the adaptive MCMC to sample better the
extremes of the prior distribution of the WMF, in comparison to the non-
adaptive MCMC. Furthermore, the acceptance rate is also improved in the
adaptive case, having an acceptance rate of 53% in comparison to the original
acceptance rate of 35% of the non-adaptive case within the same time.

Figure 4: Sampled 2D PDFs of the CMF and WMF in the ternary diagram for
TOI-220 b for non-adaptive (red) and adaptive (blue) MCMCs. The mean value of
the input Fe/Si mole ratio is indicated in a dashed black line. The mantle
mass fraction (MMF) is defined as MMF = 1 - CMF - WMF.

##  5 Planetary and observation parameters

Our MCMC analysis requires as input the planetary mass, radius and Fe/Si mole
ratio for the two planets we consider as test-case in this work, TRAPPIST-1 c
and 55 Cancri e. Their values and references are shown in Table 2. To
determine the surface temperature at which the atmosphere is in radiative-
convective equilibrium, the stellar effective temperature, stellar radius, and
semi-major axis are needed, as seen in Eq. 8. The values we adopt and their
references are given in Table 2.

Table 2: Planetary parameters for TRAPPIST-1 c and 55 Cancri e: masses, radii,
Fe/Si mole ratios, semi-major axes, and host stellar effective temperatures
and radii.

|  $M$ [$M_{\oplus}$] |  $R$ [$R_{\oplus}$] | Fe/Si |  $a_{d}$ [AU] |  $T_{\star}$ [K] |  $R_{\star}$ [$R_{\odot}$] | References
---|---|---|---|---|---|---|---
TRAPPIST-1 c | 1.308$\pm$0.056 | 1.097${}^{+0.014}_{-0.012}$ | 0.76$\pm$0.12 | 1.58$\times\ 10^{-2}$ | 2566 | 0.119 | 1, 2
55 Cancri e | 7.99${}^{+0.32}_{-0.33}$ | 1.875$\pm$0.029 | 0.60$\pm$0.14 | 1.54$\times\ 10^{-2}$ | 5172 | 0.943 | 3, 4

(1) Agol et al. (2021); (2) Unterborn et al. (2018); (3) Bourrier et al.
(2018); (4) Luck (2016)

To simulate the observations with JWST in photometry, we assume the
atmospheric parameters retrieved in our adaptive MCMC analysis and generate
emission spectra with their respective temperature-pressure profiles.
Consecutively, we bin the emission spectrum using the response functions of
each of the MIRI photometry filters666http://svo2.cab.inta-
csic.es/svo/theory/fps3/index.php?id=JWST (Glasse et al., 2015; Piette et al.,
2022). The mean flux, $\langle f_{\lambda}\rangle$, of an emission spectrum,
$f(\lambda)$, observed with a filter with transmission function $R(\lambda)$,
is defined in Eq. 32 (Stolker et al., 2020). We consider random Gaussian noise
of 100 ppm for each filter to derive the uncertainties of the mean flux
(Lustig-Yaeger et al., 2019; Piette et al., 2022).

$\langle f_{\lambda}\rangle=\dfrac{\int f(\lambda)\ R(\lambda)\ d\lambda}{\int
R(\lambda)\ d\lambda}$ (32)

For the observation of the emission spectrum of 55 Cancri e, we use Pandexo
(Batalha et al., 2020) to simulate the expected noise. Our input includes the
stellar effective temperature, and stellar and planet radius (see Table 2).
Additional input parameters are found in the database accessible by Pandexo
and ExoMast, which are shown in Table 3. We adopt observation and
instrumentation variables from Hu et al. (2021). We adopt stellar spectrum
templates generated by petitRADTRANS (Mollière et al., 2019) to calculate the
planet-to-star flux ratio with our planetary emission spectra.

Table 3: Input parameters for Pandexo to simulate observations of the emission
spectrum of 55 Cancri e with JWST’s MIRI Low Resolution Spectroscopy (LRS) and
NIRCam instruments.

Parameter | Value
---|---
Star |  |
Metallicity, $log$[Fe/H] | 0.35
Gravity, $log\ g$ [cgs] | 4.43
J Magnitude | 4.59
Planet |  |
Transit duration [d] | 0.0647
Observation |  |
Baseline [h] | 3.2
Number of eclipses | 2
Instrument | MIRI LRS | NIRCam
Mode | Slitless | F444W, subgrism 64
Saturation limit [full well] | 80%
Constant minimum noise | 100 ppm

##  6 Interior composition and simulated spectra

###  6.1 TRAPPIST-1 c

(a)

(b)

Figure 5: PDF corner plot of the atmospheric parameters of TRAPPIST-1 c for
two different atmospheric compositions. The atmospheric parameters include the
surface temperature ($T_{surf}$), atmospheric thickness ($z_{atm}$), Bond
albedo, and the surface pressure ($P_{surf}$). The two envelope compositions
we consider are H2O (left panel) and CO2 (right panel). Figure 6: Simulated
emission spectra of TRAPPIST-1 c with MIRI photometric filter mean fluxes for
water and CO2 atmospheres, obtained with the high-resolution, k-correlated
MSEIRADTRAN model. Color boxes indicate the wavelength ranges of the different
filters.

TRAPPIST-1 c has been proposed to be observed in thermal emission with the
Mid-infrared instrument (MIRI) in JWST Cycle 1 (Kreidberg et al., 2021). It
will be observed with the F1500W filter during 4 eclipses, which is the filter
centered at $\lambda$ = 15 $\mu$m. We analyse TRAPPIST-1 c with our adaptive
MCMC with a water and a CO2-dominated atmosphere. Fig. 5 shows the marginal
posterior distributions (PDF) in 1D and 2D for this analysis. The CMF agrees
well with the value obtained previously with our non-adaptive MCMC in Acuña et
al. (2021), and with the possible CMFs considered in Agol et al. (2021). Agol
et al. (2021) constrain the WMF with their MCMC and interior-atmosphere model
with an upper limit of WMF = $10^{-3}$. This upper limit is two orders of
magnitude larger than the mean of our PDF, showing a clear improvement in the
resolution of the MCMC in the low-surface pressure region of the parameter
space. In their case, a maximum WMF = $10^{-5}$ estimate can only be shown
with theoretical forward models (see their figure 18). In comparison with our
own previous work, the distribution of the WMF derived by the adaptive MCMC is
wider than the non-adaptive one, with mean and standard deviation values of
WMFadap = (11.0$\pm$5.6)$\times\ 10^{-6}$, and WMFnon-adap =
(0.0${}^{+2.7}_{-0.0}$)$\times\ 10^{-6}$, respectively. This is because the
adaptive MCMC presents a higher acceptance rate due to a more effective
exploration of the parameter space at WMF close to zero. More models in this
region are explored, and thus these are accepted in the posterior
distribution, which becomes wider with larger standard deviation than the non-
adaptive posterior distribution. The corresponding surface pressure derived by
the adaptive MCMC is $P_{surf}$ = 15$\pm$7 bar for a water-dominated envelope.
This 1$\sigma$ confidence interval must be taken carefully since the PDF of
the surface pressure does not present a Gaussian distribution shape.

A WMF of zero is also compatible with the density of TRAPPIST-1 c.
Consequently, we can conclude that TRAPPIST-1 c could have a H2O atmosphere of
up to $\simeq$ 25 bar of surface pressure, or no atmosphere at all. The
analyses performed with the k-correlated atmospheric model and the
k-uncorrelated approximation agree that a H2O atmosphere in TRAPPIST-1 c would
have a surface temperature between 1000 and 1500 K, and an atmospheric
thickness of 150 to 250 km. The retrieved Bond albedos differ by 0.05 due to
the use of different atmospheric opacity data.

For a CO2-dominated atmosphere, we retrieve a similar CMF to the water case,
although the volatile mass fraction increases by one order of magnitude, VMF =
(2.49$\pm$2.07)$\times\ 10^{-5}$. The molecular weight of CO2 is higher than
that of water vapour, producing a more compressed atmosphere for a similar
surface pressure and temperature. In addition, the radiative properties (i.e
opacity) of CO2 yields a lower surface temperature for the same irradiation
conditions in comparison to a water-dominated envelope, which contributes to a
lower atmospheric thickness. As a consequence, the models with a CO2 envelope
can accommodate a more massive atmosphere for TRAPPIST-1 c than the water
models, retrieving a surface pressure of $P_{surf}$ = 35$\pm$29 bar. The
surface temperature and atmospheric thickness are $T_{surf}$ = 807$\pm$102 K,
and $z_{atm}$ = 63$\pm$12 km, respectively.

We assume the atmospheric parameters retrieved in our adaptive MCMC analysis
and generate emission spectra with their respective temperature-pressure
profiles as explained in Sect. 5. Fig. 6 shows the complete emission spectra
and mean filter fluxes for TRAPPIST-1 c. We observe that for the nIR filters
($\lambda=5.60$ to $11.30\ \mu$m), both atmosphere models have very similar
fluxes that are compatible within uncertainties, which makes it not possible
to distinguish between the two compositions in these wavelengths. For the
filters F1500, F1800 and F2100, the mean flux uncertainties of the water and
the CO2 atmospheres do not overlap, allowing the different compositions to be
distinguished. Therefore, observing TRAPPIST-1 c in emission with filter
F1500, as proposed by Kreidberg et al. (2021), is well-adapted to
differentiate between a water or a CO2-dominated atmosphere. We also consider
the possibility of a bare surface in TRAPPIST-1 c, given the high probability
obtained in our MCMC analysis for a volatile mass fraction equal to zero. Hu
et al. (2012) obtain the emission spectra of bare terrestrial surfaces for
different minerals. We estimate the brightness temperature for the irradiation
conditions of TRAPPIST-1 c for two minerals, a metal-rich surface and a
granitoid one, since these are the two surfaces with the highest and lowest
emission for the same irradiation conditions, respectively. We approximate the
emission spectrum of these surfaces to that of a black body with a temperature
equal to the estimated brightness temperature. The emission flux in the F1500
filter for a water atmosphere is slightly higher than of the metal-rich
surface. If TRAPPIST-1 c emission flux in this filter is 731 ppm or higher, it
is indicative of a thin, water-rich atmosphere. For very low emission fluxes
($\simeq$ 300 ppm), TRAPPIST-1 c would present a CO2-dominated atmosphere. For
fluxes between 730 to 400 ppm, TRAPPIST-1 c would have no atmosphere, with an
emission that corresponds to a bare surface. The surfaces with the lowest
emission are granitoid, feldspathic, or clay (Hu et al., 2012).

###  6.2 55 Cancri e

Table 4: MCMC retrieved mean value and 1$\sigma$ uncertainties for observable (mass, radius and Fe/Si mole ratio), compositional (core and volatile mass fractions) and atmospheric (surface pressure and temperature, atmospheric thickness and Bond albedo) parameters of 55 Cancri e. | Data | H2O | CO2
---|---|---|---
$M$ [$M_{\oplus}$] | 7.99${}^{+0.32}_{-0.33}$ | 7.99${}^{+0.29}_{-0.32}$ | 8.17$\pm$0.29
$R$ [$R_{\oplus}$] | 1.875$\pm$0.029 | 1.877$\pm$0.020 | 1.850$\pm$0.021
Fe/Si | 0.60$\pm$0.14 | 0.62$\pm$0.14 | 0.48$\pm$0.13
CMF |  | 0.20$\pm$0.05 | 0.15$\pm$0.05
VMF |  | $(6.7^{+7.4}_{-5.9})\times 10^{-5}$ | $(5.1^{+5.9}_{-4.6})\times 10^{-5}$
$P_{surf}$ [bar] |  | 209$\pm$93 | ¿ 300
$T_{surf}$ [K] |  | 4161$\pm$199 | 4035$\pm$597 (at 300 bar)
$z_{atm}$ [km] |  | 522$\pm$46 | 152$\pm$33 (from 300 bar)
Bond albedo |  | 0.191$\pm$0.001 | 0.351$\pm$0.004

55 Cancri e is a super-Earth in a close orbit (P = 0.66 days) to a bright
star, for which several interior and atmospheric hypothesis have been
proposed. Madhusudhan (2012) explored a carbon-rich interior given the high
C/O ratio found for its host star, showing that in this case the planetary
bulk density would be lower than that of a silicate-rich mantle planet, such
as Earth. They concluded that a volatile layer would not be necessary to
account for its density. On the other hand, a classical Fe-rich core and a
silicate mantle are compatible with a volatile envelope rich in secondary
atmosphere species. Furthermore, the absence of a H/He-dominated envelope
seems likely due to the lack of hydrogen and helium emission and absorption
lines in the spectrum (Ehrenreich et al., 2012; Zhang et al., 2021). The
presence of a secondary atmosphere is supported by phase curve data from the
Spitzer Space Telescope (Demory et al., 2011; Angelo & Hu, 2017), and 3D GCM
modelling that suggests an optically thick atmosphere with a low mean
molecular weight (Hammond & Pierrehumbert, 2017). A fully H2O-dominated
atmosphere has been discarded, since it will require the presence of water and
hydrogen simultaneously in the atmosphere due to water dissociation. Possible
compositions for the atmosphere of 55 Cancri e are a mixture of silicate
compounds (Keles et al., 2022), such as HCN, detected by Tsiaras et al.
(2016), with traces of water (detected by Esteves et al., 2017), or CO2, CO,
and N2, among other compounds.

Hu et al. (2021) have proposed to observe 55 Cancri e in emission spectroscopy
combining NIRCam F444W filter (3-5 $\mu$m), and MIRI’s Low Resolution
Spectrograph (MIRI LRS; 5-14 $\mu$m). In Table 4, we observe that a water-
dominated atmosphere reproduces well the observed data, with a surface
pressure higher than 200 bar, whereas the CO2 envelope is not extended enough
to match the density of 55 Cancri e, yielding a more dense interior. At
temperatures higher than 4000 K, CO2 would not be the dominant species in a
C-rich atmosphere, but CO. This changes the emission of the atmosphere as CO
is a different absorber from CO2. A CO-rich atmosphere could also explain the
low-density of 55 Cancri e in this scenario, since CO has a lower molecular
weight than CO2, yielding a larger atmospheric scale height. When H/He is not
included in the interior modelling, water as a trace species is necessary to
explain the low density of 55 Cancri e, since a purely dry silicate atmosphere
will have a smaller thickness than a CO2 atmosphere due to their heavier
molecular weights under similar atmospheric surface conditions. Adding
silicate absorbers decreases the total planetary radius in H/He envelopes
(Misener & Schlichting, 2022a). However, more modelling work is necessary to
explore the effect of silicates in atmospheres that have lost their primordial
H/He. A planet with no volatiles could match the low planetary density if the
core and the mantle where less dense than that of an Earth-like interior,
pointing to a carbon-rich mantle as suggested by Madhusudhan (2012). The
emission spectrum in this scenario would be that of a bare surface, requiring
terrestrial surface models such as the ones presented in Hu et al. (2012).

Figure 7: Predicted emission spectrum with the k-correlated, high-resolution
MSEIRADTRAN model for a water-rich atmosphere in 55 Cancri e with NIRCam and
MIRI LRS. The blue error bars are computed by Pandexo by using the same input
Hu et al. (2021) provide in their proposal, which we display in Table 3.

Fig. 7 shows the complete predicted emission spectrum of 55 Cancri e from 3 to
14 $\mu$m. We notice that for wavelengths below $\simeq$ 3.7 $\mu$m the noise
is quite high, although in the range $\lambda=$ 3 to 3.7 $\mu$m there are no
water or CO2 spectral lines. Nonetheless, the rest of the spectral coverage of
the proposed observations have low noise level, which makes the spectral
features of water easy to identify with JWST in the high-molecular weight
atmosphere scenario of 55 Cancri e.

##  7 Discussion

In this section, we discuss what processes that are not included in our model
may affect our estimates on the volatile mass fraction, such as solubility in
a magma ocean, and near-surface isothermal layers.

TRAPPIST-1 c and 55 Cancri e are warm enough to hold a magma ocean on their
surface. The dissolution of silicates and the density of the mantle can change
between a dry mantle model and that of a wet magma ocean (Dorn & Lichtenberg,
2021). For a planet less massive than TRAPPIST-1 c (0.8 $M_{\oplus}$) and a
radius of $\simeq 1.1\ R_{\oplus}$, a dry mantle model and a wet magma ocean
model estimate a WMF = 3 $\times 10^{-3}$ and 8 $\times 10^{-3}$,
respectively. This difference decreases with higher planet mass and lower WMF.
Including dissolution of silicates and a magma ocean would increase our WMF
estimates by less than a factor of 2.7 (see figure 4 in Dorn & Lichtenberg,
2021). Since our MCMC realisations span at least 2 orders of magnitude (from
$10^{-6}$ to $10^{-4}$), the wet magma ocean model would only shift our WMF
distribution by less than one order of magnitude. The maximum surface pressure
of TRAPPIST-1 c would go from 25 bar to 75 bar at most. A retrieval with a wet
magma ocean would not discard a bare surface in TRAPPIST-1 c, since the radius
of a wet melt-solid interior is less than that of a dry solid interior,
leaving more room for an atmosphere above. For 55 Cancri e, which has a mass
of $8\ M_{\oplus}$, the difference is negligible, yielding similar WMF
estimates if we considered a wet magma ocean surface.

In our atmospheric model, we prescribe the atmosphere’s thermal structure as a
near surface dry convective layer, followed by a wet convective layer and an
isothermal mesosphere. A self-consistent treatment of the shortwave radiation,
together with an iterative scheme on the temperature profile would enable us
to compute the regions of the atmosphere where radiative layer would form, and
their exact temperature value. We test with our atmospheric model how changes
in the stratospheric temperature may impact the OLR and the atmospheric
thickness. We perform this test for three stratospheric temperatures: 200 K,
500 K, and 1000 K. The surface temperature at which the atmosphere is in
radiative-convective equilibrium is similar in the three cases, yielding
similar atmospheric thicknesses. Therefore, our envelope mass fraction
estimates are robust against different upper radiative temperatures. In
contrast, near-surface radiative layers may decrease the thickness of the
atmosphere compared to a convective atmosphere. However, these are more likely
to form in atmospheres composed of H/He and silicates than pure water or CO2
envelopes (Misener & Schlichting, 2022b). This means that our envelope mass
fractions estimates are lower limits compared to those that would be obtained
with an atmosphere that presents near-surface radiative layers. Moreover,
Vazan et al. (2022) find that water envelopes with silicates only develop
radiative layers at low pressures (100-10 bar, see their figure 7), and not in
the near-surface.

To calculate the density of CO2, we use the ideal gas EOS, in contrast to a
non-ideal EOS for water. In the following we discuss how the use of an ideal
EOS for carbon dioxide may affect our results. Böttcher et al. (2012) carry
out a comparison between the ideal EOS and non-ideal EOS for carbon dioxide.
They find that CO2 starts to behave as a non-ideal gas at $\simeq$ 7 MPa which
corresponds to 70 bar. In the case of TRAPPIST-1 c, the maximum surface
pressure for a CO2-dominated atmosphere is 25 bar in our retrievals, so we can
safely say that the CO2 gas in this case still has an ideal behaviour. The
non-ideal EOS yields a higher density than the ideal EOS (Böttcher et al.,
2012). This means that for 55 Cancri e, the atmospheric thickness of a CO2
envelope with a non-ideal EOS would be even lower than that calculated with
the ideal EOS. This would make it even harder to match the current density of
55 Cancri e with a CO2-dominated atmosphere, which further supports that CO is
more likely to constitute the atmosphere of this planet instead. The
implementation of non-ideal EOS, such as the SESAME EOS 5210 for carbon
dioxide (Lyon, 1992), and the Podolak et al. (2023) EOS for CO, will be the
scope of future work.

##  8 Conclusions

In this work, we present a self-consistent model built to estimate the
internal compositions and structures of low-mass planets with water and CO2
atmospheres given their observed mass, radius, and their host stellar
abundances. We couple self-consistently the interior and the atmosphere to
obtain the boundary conditions at the top of a supercritical water layer or a
silicate mantle given the irradiation conditions of the planet. This is done
by calculating the bolometric emission flux and the Bond albedo, to compute
the flux emitted and absorbed by an atmosphere in radiative-convective
equilibrium. For fast computations of the radiative quantities within our MCMC
retrieval approach, we present the k-uncorrelated approximation. We show that
the k-uncorrelated approximation presents sufficient accuracy to obtain the
radius of an irradiated water planet within a 1% error with k-correlated
models. This error increases with higher water mass fraction (WMF ¿ 70%) and
lower planetary mass. We also demonstrate that using a constant step size when
sampling the prior distribution in a MCMC scheme is not efficient in exploring
the parameter space in interior modelling. This causes an underestimation of
the uncertainties of the compositional parameters. Therefore, it is necessary
to use an adaptive MCMC when performing retrieval with interior models,
especially for planets whose compositional parameters can reach the maximum or
minimum possible values. This is the case of rocky Earth-sized planets and
super-Earths, whose WMFs are close to zero, but are nonetheless important to
determine their surface pressure.

Moreover, we use the surface pressure and temperature conditions retrieved
with our interior-atmosphere model to generate emission spectra with our
k-correlated atmospheric model, MSEIRADTRAN. We compute emission spectra to
show how the output of our interior-atmosphere model can be used to predict
the input necessary for atmosphere models, and help prepare atmospheric
characterization proposals. The particular parameter that interior models can
provide self-consistently for atmosphere models is the surface pressure, which
is usually chosen arbitrarily to generate spectra. Emission spectra are more
sensitive to the choice of surface pressure and temperature than transmission
spectra.

We showcase how to use interior and atmospheric modelling simultaneously to
predict observations for two rocky planets, TRAPPIST-1 c and 55 Cancri e,
which have been proposed for observations in emission photometry and
spectroscopy with JWST. We bin our emission spectra according to the response
functions of the MIRI filters to predict emission fluxes for TRAPPIST-1 c for
different scenarios, while for 55 Cancri e we input our emission spectra to
Pandexo to predict observational uncertainties.

The most likely scenario for TRAPPIST-1 c is the lack of an atmosphere (WMF =
0). Nonetheless, the presence of a secondary atmosphere cannot be ruled out.
In this scenario, TRAPPIST-1 c could have a H2O-dominated atmosphere of up to
25 bar of surface pressure. A CO2-rich envelope would have higher surface
pressures, and therefore it would be more massive than the water case, with a
maximum surface pressure of 64 bar. We present emission flux estimates for the
filter centered at 15 $\mu$m, F1500, that can be compared with observations by
Kreidberg et al. (2021) to determine if TRAPPIST-1 c has a bare surface. For
55 Cancri e, a massive envelope with more than 300 bar of surface pressure
that contains water is necessary to fit its low density. We determine that a
combined spectrum with NIRCam and MIRI LRS, as proposed by Hu et al. (2021),
may present a high noise level at wavelengths between 3 and 3.7 $\mu$m.
However, this part of the spectrum does not contain any spectral lines of
water or CO2, which are essential to determine the abundances in the envelope.

In our modelling approach, we have considered water and CO2-dominated
atmospheres (99% H2O:1% CO2, and viceversa). However, the atmospheres of low-
mass planets are more diverse than these two compositional scenarios. The
atmospheric compositions of sub-Neptunes are proving to be a mixture of H/He,
water and other compounds via observations and models (Madhusudhan et al.,
2020; Bézard et al., 2020; Guzmán-Mesa et al., 2022), while super-Earths can
have more exotic atmospheric compositions, such as mineral atmospheres (Keles
et al., 2022). Therefore, the scope of future work will be to include more
gases in the atmospheric model, as well as the calculation of transmission
spectra in addition to the existing implementation of emission and reflection
spectra. Our interior-atmosphere model, MSEI, serves as a precedent to develop
models with more diverse envelope compositions to prepare proposals for JWST
and future atmospheric characterization facilities, such as Ariel (Tinetti et
al., 2018). Our model can also be used within retrieval frameworks to
interpret mass, radius and upcoming JWST emission spectral data simultaneously
to break degeneracies in exoplanet composition.

###### Acknowledgements.

M.D. and O.M. acknowledge support from CNES. We acknowledge the anonymous
referee whose comments helped improve and clarify this manuscript.

## References

  * Abramowitz & Stegun (1972) Abramowitz, M. & Stegun, I. A. 1972, Handbook of Mathematical Functions
  * Acuña et al. (2021) Acuña, L., Deleuil, Magali, Mousis, Olivier, et al. 2021, A&A, 647, A53
  * Acuña et al. (2022) Acuña, L., Lopez, T. A., Morel, T., et al. 2022, A&A, 660, A102
  * Agol et al. (2021) Agol, E., Dorn, C., Grimm, S. L., et al. 2021, The Planetary Science Journal, 2, 1
  * Aguichine et al. (2021) Aguichine, A., Mousis, O., Deleuil, M., & Marcq, E. 2021, The Astrophysical Journal, 914, 84
  * Angelo & Hu (2017) Angelo, I. & Hu, R. 2017, The Astronomical Journal, 154, 232
  * Batalha et al. (2020) Batalha, N., Fix, M., Fraine, J., et al. 2020, natashabatalha/PandExo: Release 1.5
  * Baumeister et al. (2021) Baumeister, P., Tosi, N., Grenfell, J. L., & MacKenzie, J. 2021, in Plato Mission Conference 2021. Presentations and posters of the online PLATO Mission Conference 2021, 87
  * Bézard et al. (2020) Bézard, B., Charnay, B., & Blain, D. 2020, arXiv e-prints, arXiv:2011.10424
  * Bézard et al. (2011) Bézard, B., Fedorova, A., Bertaux, J.-L., Rodin, A., & Korablev, O. 2011, Icarus, 216, 173
  * Boucher (1998) Boucher, O. 1998, Journal of the Atmospheric Sciences, 55, 128
  * Bourrier et al. (2018) Bourrier, Dumusque, X., Dorn, C., et al. 2018, A&A, 619, A1
  * Brugger et al. (2017) Brugger, B., Mousis, O., Deleuil, M., & Deschamps, F. 2017, The Astrophysical Journal, 850, 93
  * Brugger et al. (2016) Brugger, B., Mousis, O., Deleuil, M., & Lunine, J. I. 2016, The Astrophysical Journal, 831, L16
  * Böttcher et al. (2012) Böttcher, N., Taron, J., Kolditz, O., et al. 2012, International Association of Hydrological Sciences-AISH publication. ModelCARE. Conference (Leipzig 2011-09)
  * Chase (1998) Chase, M. W., J. 1998, J. Phys. Chem. Ref. Data, Monograph 9.
  * Chaudhry et al. (2021) Chaudhry, S., Lautzenheiser, D., & Ghosh, K. 2021, arXiv e-prints, arXiv:2110.10840
  * Demory et al. (2011) Demory, B. O., Gillon, M., Deming, D., et al. 2011, A&A, 533, A114
  * Diamond-Lowe et al. (2020) Diamond-Lowe, H., Charbonneau, D., Malik, M., Kempton, E. M. R., & Beletsky, Y. 2020, AJ, 160, 188
  * Director et al. (2017) Director, H. M., Gattiker, J., Lawrence, E., & Wiel, S. V. 2017, Journal of Statistical Computation and Simulation, 87, 3521
  * Dorn et al. (2015) Dorn, C., Khan, A., Heng, K., et al. 2015, A&A, 577, A83
  * Dorn & Lichtenberg (2021) Dorn, C. & Lichtenberg, T. 2021, ApJ, 922, L4
  * Dorn et al. (2017) Dorn, C., Venturini, J., Khan, A., et al. 2017, A&A, 597, A37
  * Ehrenreich et al. (2012) Ehrenreich, Bourrier, V., Bonfils, X., et al. 2012, A&A, 547, A18
  * Esteves et al. (2017) Esteves, L. J., de Mooij, E. J. W., Jayawardhana, R., Watson, C., & de Kok, R. 2017, The Astronomical Journal, 153, 268
  * Fulton & Petigura (2018) Fulton, B. J. & Petigura, E. A. 2018, AJ, 156, 264
  * Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, The Astronomical Journal, 154, 109
  * García Muñoz et al. (2021) García Muñoz, A., Fossati, L., Youngblood, A., et al. 2021, ApJ, 907, L36
  * Gardner et al. (2006) Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485
  * Gillon et al. (2016) Gillon, M., Jehin, E., Lederer, S. M., et al. 2016, Nature, 533, 221
  * Glasse et al. (2015) Glasse, A., Rieke, G. H., Bauwens, E., et al. 2015, Publications of the Astronomical Society of the Pacific, 127, 686
  * Gordon (1994) Gordon, S. 1994, Tech. rep. NASA Lewis Research Center.
  * Grimm et al. (2018) Grimm, S. L., Demory, Brice-Olivier, Gillon, Michaël, et al. 2018, A&A, 613, A68
  * Grimm & Heng (2015) Grimm, S. L. & Heng, K. 2015, The Astrophysical Journal, 808, 182
  * Guzmán-Mesa et al. (2022) Guzmán-Mesa, A., Kitzmann, D., Mordasini, C., & Heng, K. 2022, MNRAS, 513, 4015
  * Haar et al. (1984) Haar, L., Gallagher, J. S., & Kell, G. S. 1984
  * Haldemann et al. (2020) Haldemann, J., Alibert, Y., Mordasini, C., & Benz, W. 2020, A&A, 643, A105
  * Hammond & Pierrehumbert (2017) Hammond, M. & Pierrehumbert, R. T. 2017, The Astrophysical Journal, 849, 152
  * Hoyer et al. (2021) Hoyer, S., Gandolfi, D., Armstrong, D. J., et al. 2021, Monthly Notices of the Royal Astronomical Society, 505, 3361
  * Hu et al. (2021) Hu, R., Brandeker, A., Damiano, M., et al. 2021, Determining the Atmospheric Composition of the Super-Earth 55 Cancri e, JWST Proposal. Cycle 1, ID. #1952
  * Hu et al. (2012) Hu, R., Ehlmann, B. L., & Seager, S. 2012, The Astrophysical Journal, 752, 7
  * Ingersoll (1969) Ingersoll, A. P. 1969, Journal of Atmospheric Sciences, 26, 1191
  * Kasting (1988) Kasting, J. F. 1988, Icarus, 74, 472
  * Keles et al. (2022) Keles, E., Mallonn, M., Kitzmann, D., et al. 2022, MNRAS, 513, 1544
  * Kimura & Ikoma (2022) Kimura, T. & Ikoma, M. 2022, Nature Astronomy, 6, 1296
  * Kopparapu et al. (2013) Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, The Astrophysical Journal, 765, 131
  * Kreidberg et al. (2021) Kreidberg, L., Agol, E., Bolmont, E., et al. 2021, Hot Take on a Cool World: Does Trappist-1c Have an Atmosphere?, JWST Proposal. Cycle 1, ID. #2304
  * Kreidberg et al. (2019) Kreidberg, L., Koll, D. D. B., Morley, C., et al. 2019, Nature, 573, 87
  * Krissansen-Totton et al. (2021) Krissansen-Totton, J., Galloway, M. L., Wogan, N., Dhaliwal, J. K., & Fortney, J. J. 2021, ApJ, 913, 107
  * Leconte (2021) Leconte, J. 2021, A&A, 645, A20
  * Leconte et al. (2013) Leconte, J., Forget, F., Charnay, B., Wordsworth, R., & Pottier, A. 2013, Nature, 504, 268
  * Liggins et al. (2022) Liggins, P., Jordan, S., Rimmer, P. B., & Shorttle, O. 2022, Journal of Geophysical Research (Planets), 127, e07123
  * Liou (1980) Liou, K. N. 1980, An introduction to atmospheric radiation.
  * Lopez & Fortney (2014) Lopez, E. D. & Fortney, J. J. 2014, ApJ, 792, 1
  * Luck (2016) Luck, R. E. 2016, The Astronomical Journal, 153, 21
  * Lupu et al. (2014) Lupu, R. E., Zahnle, K., Marley, M. S., et al. 2014, The Astrophysical Journal, 784, 27
  * Lustig-Yaeger et al. (2019) Lustig-Yaeger, J., Meadows, V. S., & Lincowski, A. P. 2019, The Astronomical Journal, 158, 27
  * Lyon (1992) Lyon, S. P. 1992, Los Alamos National Laboratory report LA-UR-92-3407
  * Ma & Tipping (1992) Ma, Q. & Tipping, R. H. 1992, The Journal of Chemical Physics, 96, 8655
  * Madhusudhan (2012) Madhusudhan, N. 2012, The Astrophysical Journal, 758, 36
  * Madhusudhan et al. (2020) Madhusudhan, N., Nixon, M. C., Welbanks, L., Piette, A. A. A., & Booth, R. A. 2020, ApJ, 891, L7
  * Malik et al. (2017) Malik, M., Grosheintz, L., Mendonça, J. M., et al. 2017, The Astrophysical Journal, 153, 56
  * Malik et al. (2019) Malik, M., Kitzmann, D., Mendonça, J. M., et al. 2019, The Astrophysical Journal, 157, 170
  * Marcq (2012) Marcq, E. 2012, Journal of Geophysical Research: Planets, 117
  * Marcq et al. (2008) Marcq, E., Bézard, B., Drossart, P., et al. 2008, Journal of Geophysical Research (Planets), 113, E00B07
  * Marcq et al. (2017) Marcq, E., Salvador, A., Massol, H., & Davaille, A. 2017, Journal of Geophysical Research (Planets), 122, 1539
  * Mazevet et al. (2019) Mazevet, S., Licari, A., Chabrier, G., & Potekhin, A. Y. 2019, A&A, 621, A128
  * McBride (1996) McBride, B. J. G. 1996, Tech. rep. NASA Lewis Research Center.
  * Misener & Schlichting (2022a) Misener, W. & Schlichting, H. E. 2022a, MNRAS, 514, 6025
  * Misener & Schlichting (2022b) Misener, W. & Schlichting, H. E. 2022b, MNRAS, 514, 6025
  * Mollière et al. (2019) Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67
  * Mollière (2017) Mollière, P. M. 2017, PhD thesis, Heidelberg University
  * Mosegaard & Tarantola (1995) Mosegaard, K. & Tarantola, A. 1995, Journal of Geophysical Research (Solid Earth), 100, 12,431
  * Mousis et al. (2020) Mousis, O., Deleuil, M., Aguichine, A., et al. 2020, The Astrophysical Journal Letters, 896, L22
  * Mousis et al. (2019) Mousis, O., Ronnet, T., & Lunine, J. I. 2019, ApJ, 875, 9
  * Nakajima et al. (1992) Nakajima, S., Hayashi, Y.-Y., & Abe, Y. 1992, Journal of Atmospheric Sciences, 49, 2256
  * Ortenzi et al. (2020) Ortenzi, G., Noack, L., Sohl, F., et al. 2020, Scientific Reports, 10, 10907
  * Piette et al. (2022) Piette, A. A. A., Madhusudhan, N., & Mandell, A. M. 2022, MNRAS, 511, 2565
  * Pluriel et al. (2019) Pluriel, W., Marcq, E., & Turbet, M. 2019, Icarus, 317, 583
  * Podolak et al. (2023) Podolak, M., Levi, A., Vazan, A., & Malamud, U. 2023, Icarus, 394, 115424
  * Polyansky et al. (2018) Polyansky, O. L., Kyuberis, A. A., Zobov, N. F., et al. 2018, Monthly Notices of the Royal Astronomical Society, 480, 2597
  * Rothman et al. (2010) Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010, Journal of Quantitative Spectroscopy and Radiative Transfer, 111, 2139
  * Sanchez-Lavega (2011) Sanchez-Lavega, A. 2011, Contemporary Physics, 52, 487
  * Simonelli & Veverka (1988) Simonelli, D. P. & Veverka, J. 1988, Icarus, 74, 240
  * Sneep & Ubachs (2005) Sneep, M. & Ubachs, W. 2005, Journal of Quantitative Spectroscopy and Radiative Transfer, 92, 293
  * Stamnes et al. (2017) Stamnes, K., Tsay, S. C., Jayaweera, K., et al. 2017, DISORT: DIScrete Ordinate Radiative Transfer
  * Stolker et al. (2020) Stolker, Quanz, S. P., Todorov, K. O., et al. 2020, A&A, 635, A182
  * Tinetti et al. (2018) Tinetti, G., Drossart, P., Eccleston, P., et al. 2018, Experimental Astronomy, 46, 135
  * Tran et al. (2018) Tran, H., Turbet, M., Chelin, P., & Landsheere, X. 2018, Icarus, 306, 116
  * Tsiaras et al. (2016) Tsiaras, A., Rocchetto, M., Waldmann, I. P., et al. 2016, The Astrophysical Journal, 820, 99
  * Unterborn et al. (2018) Unterborn, C. T., Desch, S. J., Hinkel, N. R., & Lorenzo, A. 2018, Nature Astronomy, 2, 297
  * Vazan et al. (2022) Vazan, A., Sari, R., & Kessel, R. 2022, ApJ, 926, 150
  * Vinet et al. (1989) Vinet, P., Rose, J. H., Ferrante, J., & Smith, J. R. 1989, Journal of Physics: Condensed Matter, 1, 1941
  * Wagner & Pruß (2002) Wagner, W. & Pruß, A. 2002, Journal of Physical and Chemical Reference Data, 31, 387
  * Zahnle et al. (2009) Zahnle, K., Marley, M. S., Freedman, R. S., Lodders, K., & Fortney, J. J. 2009, ApJ, 701, L20
  * Zhang et al. (2021) Zhang, M., Knutson, H. A., Wang, L., et al. 2021, The Astronomical Journal, 161, 181
  * Zieba et al. (2022) Zieba, S., Zilinskas, M., Kreidberg, L., et al. 2022, arXiv e-prints, arXiv:2203.00370

##  Appendix A K-uncorrelated approximation

We compare the behaviour of the k-uncorrelated approximation with the
k-correlated model and a grey model. In the grey model, the line opacity is
constant with wavelength, being 0.01 m2/kg for H2O. This grey opacity is
benchmarked with non-grey atmospheric models (Nakajima et al. 1992; Marcq et
al. 2017). In the case of water, a grey opacity of 0.01 m2/kg is
representative of the opacity of water in the 8-20 $\mu$m spectral window at
the Standard Reference Point (Ingersoll 1969).

For our comparison, we consider a water-dominated composition, consisting of a
99% water and 1% CO2 in mass. Figs. 8 and 9 show the dependence of the
outgoing longwave radiation (OLR), and the Bond albedo, as functions of the
temperature at the bottom of the atmosphere. In Fig. 8, we consider the case
of a planet with bulk mass $M=1.308\ M_{\oplus}$, and bulk radius $R=1.075\
R_{\oplus}$, which correspond to the mass and radius comprised between the
centre and the top of the interior model for TRAPPIST-1 c, as estimated in our
MCMC simulations. The atmospheric surface pressure is set to 1 bar, while the
stellar effective temperature and radius are equal to those of TRAPPIST-1 c.
We can see that the k-uncorrelated approximation increases with higher surface
temperatures, similar to the grey model. The Bond albedo agrees within less
than 1% between the k-uncorrelated approximation and the k-correlated model.
In addition, in Fig. 9 we examine the case of an atmosphere with a surface
pressure of 300 bar, with an underlying bulk mass of $M=7.99\ M_{\oplus}$, and
bulk radius $R=1.80\ R_{\oplus}$, which have been obtained by our retrieval on
55 Cancri e, with its corresponding stellar host parameters for the
calculation of the Bond albedo. Similar to the less massive case of TRAPPIST-1
c, the OLR in k-uncorrelated approximation increases with surface temperature,
while the Bond albedo agrees within less than 5% in the surface temperature
range relevant to 55 Cancri e (3900 to 4400 K).

Given this comparison between the k-uncorrelated approximation, and the grey
and k-correlated models in OLR and Bond albedo, the final surface temperature
for an atmosphere in radiative-convective equilibrium is 150 to 200 K in both
cases. As shown in Fig. 2, this difference in the boundary temperature for the
interior model produces differences in radius of 1% approximately. Therefore,
using a grey or k-correlated model instead of the k-uncorrelated approximation
would leave our conclusions on the interior structure of TRAPPIST-1 c and 55
Cancri e unchanged.

Figure 8: Outgoing longwave radiation (OLR) and Bond albedo as function of the
surface temperature for a grey model, our k-correlated model, and the
k-uncorrelated approximation, assuming the water-dominated case of TRAPPIST-1
c (see text). Figure 9: Outgoing longwave radiation (OLR) and Bond albedo as
function of the surface temperature for a grey model, our k-correlated model,
and the k-uncorrelated approximation, assuming the water-dominated case of 55
Cancri e (see text).