© The Authors 2023

1 Introduction

Low-mass exoplanets (M < 20 M_⊕) 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_⊕, while the radii of sub-Neptunes correspond to R = 2.4 R_⊕ (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 such a planet could have a thin atmosphere or a bare rock surface. Atmospheres containing H/He have been discarded since a very small fraction of H/He entails a minimum radius of ≃1.6 R_⊕ (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 variety in the possible atmospheric composition produces a degeneracy in the internal structure of super-Earths and Earth-sized planets, as the same planetary mass and radius can be explained by 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 (CMF) similar to that of Earth (CMF = 32%).

This degeneracy in interior structure can only be broken with the support of atmospheric characterisation data. The presence of an atmosphere has been confirmed in the hot super-Earth π Mensae c, whose detected C II ions indicate 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), and transmission spectroscopy has been used 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 as well as the surface pressure and temperature of water- and CO₂-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 application with our k-correlated atmospheric model, MSEIRADTRAN, to generate emission spectra to predict observations with JWST with the Mid-Infrared Instrument (MIRI) photometric filters and NIR-Cam and MIRI Low-Resolution Spectrometer (MIRI LRS). Combined interior and atmospheric models have been used to constrain the water mass fractions (WMFs) of rocky planets as WMF < 10⁻³ (Agol et al. 2021). However, the consistent exploration of the parameter space in the region close to WMF = 0 is necessary to accurately obtain the PDFs of the WMF and the surface pressure. 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 CO₂ atmospheres in rocky planets.

We describe 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. In Sect. 4, we detail the implementation of the adaptive MCMC, and we show an example of the retrieval with it as well as with a non-adaptive MCMC. With our model, we assess the observability of two planets that have been proposed for observations in Cycle 1 of JWST: TRAPPIST-1 c (Gillon et al. 2016; Grimm et al. 2018) and 55 Cancri e (Ehrenreich et al. 2012; Bourrier et al. 2018). In Sect. 5, we summarise the planet and instrument parameters we use as input for our interior-atmospheric analyses and Pandexo (Batalha et al. 2020) to predict uncertainties in JWST observations. We present our results and conclusions in Sects. 6 and 8, respectively.

2 Interior structure model

2.1 Physical model and equation of state (EOS)

In this section, we review 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: the CMF and WMF. The CMF is defined as the mass of the Fe-rich core divided by the total planetary mass, while the WMF is the mass of the hydrosphere divided by the total planetary mass. In the 1D interior model, the planetary radius, r, is represented by a 1D grid. Along this grid, the pressure, P(r), the temperature, T(r), the gravity acceleration, g(r), and the density, ρ(r), are calculated at each point. These four variables were obtained by solving the corresponding equation that defines its behaviour. The pressure was computed by integrating the equation of hydrostatic equilibrium (see Eq. (1)), while the temperature required integrating the adiabatic gradient profile (Eq. (2)). In low-mass planets, the opacity in their deep interior is high enough for the radiative temperature gradient to be greater than the adiabatic gradient, making the layers unstable against convection, according to the Schwarzschild criterion. In Eq. (2), γ and ϕ correspond to the Grüneisen and seismic parameters, respectively. The former describes the behaviour of the temperature in a crystal relative 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 the seismic parameter can be seen as being 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 was 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: (1) (2) (3) (4)

The density, p(r), was computed with the equation of state (EOS), which provides the density as a function of temperature and pressure. The interior structure model was divided into three separate layers: an Fe-rich core, a mantle rich in silicates, and a water layer. We used a different EOS to calculate the density for each of these layers. We adopted 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 used 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 their effects on the total radius of the planet can be found in Aguichine et al. (2021).

The final input for our interior structure model were the surface temperature and pressure. Together with the gravitational acceleration at the centre of the planet, whose value is zero, g(r = 0) = 0, these are the boundary conditions. Finally, the mass of each planetary layer was 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: (5)

2.2 Interior-atmosphere coupling

The surface pressure for the interior model depends on the atmospheric mass on top of the outermost interface of the interior model. For envelopes whose bottom pressure is greater than 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 WMF takes the mass of the atmosphere, M_atm, into account. The atmospheric mass is calculated as shown in Eq. (6), where P_base is the pressure at the base of the atmosphere (atmosphere-interior interface), R_bulk is the radius from the centre of the planet to 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 used 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, not at P_crit = 220 bar: (6)

The atmospheric model calculates the outgoing long-wave 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 σ is the Stefan–Boltzmann constant, and T_eq is the planetary equilibrium temperature. This requires knowledge of the semi-major axis of the planet, a_d, as well as the stellar radius and effective temperature, R_★ and T_★, respectively: (7) (8)

For a constant planetary mass and radius, the temperature at the base of the atmosphere can be found by solving OLR(T_base) − F_abs(T_base) = 0 with a root-finding method, such as the bisection method. Then, this root 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 was computed as the sum of the bulk radius calculated by the interior model, and the atmospheric thickness was 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 when 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 developed our own atmospheric model, MSEIRADTRAN. We started the development of MSEIRADTRAN by modifying the atmosphere model presented in Marcq et al. (2017)¹ to include up-to-date opacity and EOS data. In the following, we summarise the basic structure and principles of MSEIRADTRAN and the atmospheric models presented in Marcq et al. (2017); Pluriel et al. (2019).

We considered two scenarios for the composition of the envelope: water-dominated envelopes (99% water plus 1% CO₂) and CO₂-dominated envelopes (99% CO₂ plus 1% water). Including a wider variety of relative mass fractions between water and CO₂ in our models would only increase the degeneracies between atmospheric mass and atmospheric composition. Therefore, we only considered the two end-members to assess the observability of water and CO₂ spectral features with JWST. We did not model 100% pure water or CO₂ atmospheres because such pure compositions are very unlikely due to outgassing and atmospheric escape (Krissansen-Totton et al. 2021). In addition, to make our comparison between MSEIRADTRAN and the atmospheric model of Pluriel et al. (2019) consistent (see Sect. 3.5), we used the same exact compositions of 99%:1% instead of 100% pure water or CO₂.

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 an isothermal mesosphere on top. If the surface temperature is cold enough to allow for condensation of water, the dry troposphere does not exist. For the isothermal mesosphere, we assumed 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 did not take into account mesospheric stellar heating, which could significantly increase 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 1D 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 presented in Sect. 3.1.

The calculations of the emission spectrum and the Bond albedo were performed by bands. We divided the spectrum from 0 to 10 100 cm⁻¹ (equivalent to ≥1 µm in wavelength) into 36 bands to obtain the OLR, similar to Pluriel et al. (2019). For each band, we calculated 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.3), Rayleigh scattering, and clouds. We treated Rayleigh scattering as was done in Pluriel et al. (2019), where the Rayleigh scattering opacity is related to wavelength, λ, following Eq. (9). The parameters κ₀ and λ₀ were adopted from Kopparapu et al. (2013) and Sneep & Ubachs (2005) for H₂O, and CO₂, respectively. The opacity of clouds was considered for the atmospheric layers where condensation takes place. Similar to Marcq et al. (2017) and Pluriel et al. (2019), the cloud opacity was parameterised after Kasting (1988), who assumed 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): (9) (10) (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 was calculated for all 36 bands, which were then summed to obtain the bolometric, wavelength-integrated TOA flux, or OLR. Immediately after the OLR was computed, we started the calculation of the reflection spectrum and the Bond albedo, described in Sect. 3.4: (12)

3.1 Atmospheric equation of state

The thermal structure of the atmosphere is divided into two main layers. The near-surface layer is adiabatic, which means convection takes place, while the layer on top is a mesosphere with a constant temperature. We set this temperature 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 into two sub-layers: a dry convective layer and a wet convective layer. Condensation may take place depending on the value of the pressure compared to the water saturation pressure. We considered the following envelope compositions: (1) 99% water and 1% CO₂ and (2) 1% water and 99% CO₂. Water is a condensable species, while CO₂ is a non-condensable gas. To obtain the temperature in a computational layer, i, we considered two approximations. The first is that the change in temperature and pressure within an atmospheric layer is small enough to approximate . The second approximation is ∆P ~ P_i ∆ln(P) because . The temperature of a given computational layer, T_i, was then 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, (∂T∂P)_S, as seen in Eq. (13): (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 > Tcri, we are under dry convection (Marcq 2012; Marcq et al. 2017). We calculated the derivative (dT/dP)_S in the dry case following (Marcq et al. 2017; Eq. (14)). The densities of water vapour and CO₂ are ρ_v and ρ_c, and C_p,v and C_p,c are their heat capacities. The specific volume of water vapour is indicated by V_v = 1/ρ_v: (14)

When the atmospheric pressure reaches the water saturation curve, P = P_sat(T), water vapour starts to condense 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 is provided in Eq. (15) (Marcq et al. 2017). We note that CO₂ is the only non-condensable gas. The molecular weight of carbon dioxide is M_c, and C_v,c is the specific heat capacity at a constant volume of CO₂. The ideal gas constant is R, while α_v is the mixing ratio of the water vapour density relative to CO₂, α_v = ρ_v /ρ_c. Its derivative, ∂ ln(α_v)/∂ ln(T) (Eq. (16)), needed to be calculated as well (Kasting 1988; Marcq 2012; Marcq et al. 2017) (15) (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 used 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– 1273 K in temperature and up to 1 GPa in pressure. Haldemann et al. (2020) transitioned at 1200 K to an EOS that is valid at low pressures and high temperatures. This second EOS is the Chemical Equilibrium with Applications (CEA) 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 the gas phase. Figure 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, κ_ad (Eq. (17)). The reduction of the dry adiabatic coefficient at T = 1000–2500 K is due to thermal dissociation, whereas the decrease at higher temperatures (T ≥ 6000 K) is caused by thermal ionisation (Haldemann et al. 2020): (17)

We assumed that CO₂ is an ideal gas to calculate its density, ρ_c. As we treat CO₂ as an ideal gas, we can calculate its heat capacity C_v,c = C_p,c – R/M_c. We calculated the specific heat of CO₂ by using Shomate’s equation (Eq. (18)), whose tabulated coefficients A to E are provided by Chase (1998). This is part of the National Institute of Standards and Technology (NIST) chemistry WebBook database², and it is valid for temperatures up to 6000 K: (18)

Finally, we determined 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 ρ_total,i is the total mass density at altitude z_i, ρ_total,i = ρ_v,i + ρ_c,i. This expression was derived by approximating (dP/dz) ~ ∆P/z_i−1 – z_i. The gravity acceleration was computed at each point of the 1D grid with the equation for hydrostatic equilibrium. This is noted explicitly in Eq. (19) by the labelling of the gravity at the computational layer i as g_i. The final atmospheric thickness was obtained by evaluating z at the transit pressure, 20 mbar: (19)

Fig. 1 Dry adiabatic coefficient for water, κ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 correspond to κad = 0.35, 0.25, and 0.15, respectively.

3.2 K-correlated method implementation

We employed the k-correlated method (Liou 1980) to compute the emission spectrum. We discretised the dependence of the opacity on wavenumber, v, by dividing the spectrum into bands. The spectral transmittance of a spectral band b, 𝒯_b, is defined as the exponential sum of the opacity of the band, κ^b, times the column density m, which is only dependent on pressure and temperature (Sanchez-Lavega 2011). The exponential sum was performed over G-points, which are the abscissa values chosen to discretise the cumulative probability function of the opacity, G(κ). 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 . The total number of G-points is indicated by N_G: (20) (21)

In the case of an arbitrary set of G-points, the weights are equal to w_i = ∆G_i, which are the widths of the bins in the 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 integrated 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, , was obtained by invoking the radiative solver for a total optical depth whose line optical depth is (Eq. (22)): (22)

Following Malik et al. (2017, 2019), we defined our discrete G-points as the roots of the Legendre polynomial, G_LG,i (Eq. (23)). The corresponding weights, w_i, are the Legendre-Gaussian (LG) weights associated with the Legendre polynomial of the N_Gth order, (Eq. (24)). The LG weights were calculated from the ith root of the NG th order Legendre polynomial, y_i, as well as from the polynomial’s derivative, (Abramowitz & Stegun 1972; Malik et al. 2017). We chose to have N_G = 16 G-points: (23) (24)

To generate emission spectra and assess the observability with JWST, we used the original spectral resolution of the opacity k-table data (see Sect. 3.3) in our k-correlated model. This spectral resolution corresponds to R = 200 to 300 in the spectral range λ = 1 to 20 µm.

3.3 Opacity data

The total optical depth computed in the atmospheric model includes contributions from CIA and line absorption. The CIA absorption is particularly important in dense gases, such as steam and CO₂ at high pressures, especially if the line opacity is weak (Pluriel et al. 2019). For the atmospheric compositions we considered in this work, namely 99% H₂O:1% CO₂ and 99% CO₂:1% H₂O, we required CIA data for collisions between these two gases and their self-induced absorptions. We adopted CIA absorption data for H₂O–CO₂ and H₂O–H₂O collisions provided by Ma & Tipping (1992) and Tran et al. (2018)³, respectively. The CO₂–CO₂ CIA opacities are 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 H₂O–H₂O CIA table covers the complete spectral range where we calculated both our emission and reflection spectra, while the H₂O–CO₂ CIA table covers the bands with wavelength ≥1 µm, which corresponds to the bands necessary for the emission spectrum only. For the bands whose wavelengths are outside the spectral range of the CIA table, we assumed a constant CIA opacity value equal to the opacity at the limit band of the table.

Grimm & Heng (2015) provide a database⁴ of pre-calculated opacity k-tables for different species and line lists. For waterand CO₂-dominated atmospheres, we adopted 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 the HITEMP maximum temperature is 4000 K. Following the procedure described in Leconte (2021), we binned the k-correlated opacities to the same spectral bins of Marcq et al. (2017) and Pluriel et al. (2019). We calculated the k-coefficients for our water–CO₂ mixture by assuming that the spectral features of the individual gases are correlated (Malik et al. 2017). In the correlated approximation, the mixed opacity was estimated as indicated in Eq. (25), whereχ_j is the mixing ratio by mass of the jth gas and κ_j,i is the k-coefficient of the jth gas evaluated at the G_i point. The mixing ratio by mass is defined as , where is the volume mixing ratio of the jth species, MW_j is its molecular weight, and µ is the mean molecular weight of the mixture: (25)

3.4 Reflection spectra and Bond albedo

Once the bolometric OLR was obtained, we initiated the calculation of the reflectivity in 30 bands, from 5 to 0.29 µm, to obtain the Bond albedo (Pluriel et al. 2019). The bands for which we calculated both the emission flux and the reflectivity (from 1 to 5 µm) could not have the two quantities calculated simultaneously since DISORT requires different input settings to calculate them. For the emission, we assumed 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), θ in Eq. (26). This is the angle that the incident light forms with the normal of the incident surface. Once we obtained the dependence of the reflectivity with SZA, we could average it as indicated in Eq. (26) (Simonelli & Veverka 1988). To integrate Eq. (26), we evaluated the reflectivity at ten different SZA values. We assumed four streams for DISORT in both the calculation of the OLR and the reflectivity: (26)

After averaging the reflectivity over SZA, we obtained the reflection spectrum, which is the dependence of the albedo as a function of the wavenumber. To obtain the bolometric Bond albedo, we integrated Eq. (27) (Pluriel et al. 2019). The variable A_B(ν) is the reflectivity as a function of the wavenumber, B_ν(T_★) is Planck’s function for a temperature equal to the effective temperature of the host star, and σ is the Stefan-Boltzmann constant: (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 assumed the Henyey-Greenstein phase function for clouds. DISORT requires the calculation of the Legendre moments of the combined phase function, which we estimated as the weighted average of the moments of the two individual phase functions (Liou 1980; Boucher 1998). The weights were calculated as the ratio of the optical depth due to Rayleigh scattering or clouds divided by the total optical depth, τ_Rayleigh + τ_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).

As input, DISORT also requires the single scattering albedo of each atmospheric layer. The single scattering albedo is defined as the ratio of scattering efficiency to total extinction efficiency. The total extinction is a sum of both extinction by scattering and extinction by absorption. Therefore, a single scattering albedo of one 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 estimated the combined single scattering albedo from gas (Rayleigh) and clouds with their weighted average. The single scattering albedo due to Rayleigh scattering was calculated as the Rayleigh optical depth divided by the total optical depth, τ_clear + τ_Rayleigh. The clear optical depth is the sum of the line and CIA optical depths (Sect. 3.3). For fast computations of the Bond albedo within our interior–atmosphere model in retrievals, we used a grey model for the line opacity. The line opacity is constant with wavelength and is 0.01 m² kg⁻¹ for H₂O. 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 m² kg⁻¹ is representative of the opacity of water in the 8–20 µm spectral window at the Standard Reference Point (Ingersoll 1969). The single scattering albedo due to clouds was calculated as the ratio of the clouds optical depth divided by τ_clear + τ_cloud, times the cloud single scattering albedo defined in Kasting (1988): (28)

3.5 Comparison to previous models

To compare the effect of the temperature at the interior– atmosphere boundary on the total radius of planets with water-rich envelopes, we computed two sets of mass-radius relationships (Fig. 2, left panel). The first set was obtained by coupling the interior model with our k-correlated model, and the second was obtained with that of Pluriel et al. (2019). The difference in radius between the two models is less than 1% at all masses and water contents. The difference in temperature between the two models is within 50 K for WMFs of 1% and 20% (Fig. 2, right panel). For WMF = 70% and masses below 8 M_⊕, the difference in temperature between the two models can reach up to 130 K. This discrepancy is caused by differences in the OLR between k-correlated models that use different opacity data. Nonetheless, planets with M < 8 M_⊕ are unlikely to accrete WMFs above 50% (Miguel et al. 2020; Kimura & Ikoma 2022). This means that differences in opacity data in atmospheric models are unlikely to affect mass-radius relations and interior structure retrievals of detected exoplanets. We tested approximations to atmospheric models, including grey models, and find that differences in interface temperature greater than 150 K can produce changes in radius of more than 1% (not shown).

Fig. 2 Comparison of the radius and interior–atmosphere boundary temperature between our fc-correlated model (MSEIRADTRAN) and that of Pluriel et al. (2019). Left panel: mass-radius relationships for a planet with a water-dominated atmosphere orbiting a Sun-like star at ad = 0.05 AU. Dashed lines indicate the total radius calculated by MSEIRADTRAN, while the solid line corresponds to the interior radius, which comprises the core, mantle, and supercritical water (SW) hydrosphere. 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.

4 Markov chain Monte Carlo

4.1 Introduction to MCMC

The MCMC Bayesian algorithm described in Dorn et al. (2015) was later adapted by Acuña et al. (2021) to our forward interior-atmosphere model. In this work, we use it to perform the retrievals. 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 . The available data are the total mass M, the total radius R, and the Fe/Si abundance, d = {M_obs,R_obs,Fe/Si_obs}, whose observational errors are σ(M_obs), σ(R_obs), σ(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 = {M_obs,R_obs}. The prior information consists of a Gaussian distribution centred on the mean value of the observed mass with a standard deviation equal to the observational uncertainty. For the CMF and WMF, we considered uniform distributions ranging from zero to one as priors. The MCMC scheme starts by first drawing a value for each of the model parameters from their prior distributions, which we denote as . The interior model calculates the planetary radius and Fe/Si abundance that correspond to these model parameters, which is g(m₁) = {R₁, M₁,Fe/Si₁}. We then computed the likelihood of this model (Eqs. (29) and (30)), and we drew a new model from the prior distributions, m₂: (29) (30)

The log-likelihoods, l(m_i | d) = log(L(m_i | d)), of both models were used to estimate the acceptance probability (Eq. (31)). Consecutively, a random number was drawn from a uniform distribution between zero and one. If P_accept was greater than this random number, m₂ was then accepted, and the chain moved to this set of model parameters, starting the following chain n + 1. Otherwise, the chain remained in m₁, and a different set of model parameters was proposed, m₃. The accepted models were stored, and the values of their parameters composed the PDF that would enable us to estimate their mean and uncertainties: (31)

4.2 Adaptive Markov chain Monte Carlo

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

We compared the non-adaptive and adaptive MCMC for one planet, TOI-220 b (Hoyer et al. 2021). We considered as input data, the total mass and radius as well as the Fe/Si mole ratio, which was calculated with the stellar abundances of the host star. We did not establish maximum limits for the CMF and the WMF. The planet 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 discrepancy in uncertainty indicates that the non-adaptive MCMC is not as effective as the adaptive MCMC at sampling all possible pairs that could reproduce the mass and radius data. As a consequence, the uncertainties of the WMF are 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.

In Fig. 3, 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 WMFs 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 better sample 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% for the non-adaptive case within the same time.

5 Planetary and observation parameters

Our MCMC analysis required as input the planetary mass, radius, and Fe/Si mole ratio for the two planets we considered as test cases 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 were needed, as seen in Eq. (8). The values we adopted and their references are given in Table 2.

To simulate observations with JWST in photometry, we assumed the atmospheric parameters retrieved in our adaptive MCMC analysis and generated emission spectra with their respective temperature-pressure profiles. Consecutively, we binned the emission spectrum using the response functions of each of the MIRI photometry filters⁶ (Glasse et al. 2015; Piette et al. 2022). The mean flux, 〈f_λ〉, of an emission spectrum, f(λ), observed with a filter with transmission function R(λ), is defined in Eq. (32) (Stolker et al. 2020). We considered a random Gaussian noise of 100 ppm for each filter in order to derive the uncertainties of the mean flux (Lustig-Yaeger et al. 2019; Piette et al. 2022): (32)

For the observation of the emission spectrum of 55 Cancri e, we used Pandexo (Batalha et al. 2020) to simulate the expected noise. Our input included the stellar effective temperature as well as the stellar and planet radius (see Table 2). Additional input parameters can be found in the database accessible by Pandexo and ExoMast, which are shown in Table 3. We adopted observation and instrumentation variables from Hu et al. (2021).

Fig. 3 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 with a dashed black line. The mantle mass fraction (MMF) is defined as MMF = 1 – CMF -WMF.

6 Interior composition and simulated spectra

6.1 TRAPPIST-1 c

TRAPPIST-1 c has been proposed to be observed in thermal emission with MIRI in JWST Cycle 1 (Kreidberg et al. 2021). It will be observed with the F 1500W filter during four eclipses, which is the filter centred at λ = 15 μm. We analysed TRAPPIST-1 c with our adaptive MCMC for a water- and a CO₂ -dominated atmosphere. Figure 4 shows the marginal posterior distributions in ID 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⁻³. 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⁻⁵ estimate can only be shown with theoretical forward models (see their Fig. 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 when we do not consider the stellar Fe/Si constraint, and WMF_adap = 3.3±3.3 × 10⁻⁵ for the stellar Fe/Si scenario, while the non-adaptive value is . This difference 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 were explored, and thus more models are accepted in the posterior distribution, which becomes wider with a larger standard deviation than the non-adaptive posterior distribution. The corresponding surface pressure derived by the adaptive MCMC is P_surf = 40±40 bar for a water-dominated envelope (Fig. 4).

We can conclude that TRAPPIST-1 c could have an H₂O atmosphere of up to ≃80 bar of surface pressure, an atmosphere with a less volatile species (CO₂, O₂), or no atmosphere at all. The analyses performed with the k-correlated atmospheric model indicate that an H₂O atmosphere in TRAPPIST-1 c would have a surface temperature between 1000 and 2000 K and an atmospheric thickness of 150–250 km. The posterior distributions of the surface temperature and atmospheric thickness show a bimodal shape in Fig. 4 due to the inclusion of the stellar Fe/Si constraint in our retrieval. The peak at low values of z_atm corresponds to low CMFs (≃0.20), while the second peak corresponds to a CMF ≃ 0.40. If we only use the mass and radius as observable constraints in our retrieval, the first peak disappears and the PDF becomes Gaussian.

For a CO₂-dominated atmosphere, we retrieved a CMF similar to the water case. For the volatile mass fraction, the PDF of the surface pressure is approximately uniform (see Fig. 4, right panel). The molecular weight of CO₂ 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 CO₂ yield 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 CO₂ envelope can accommodate a more massive atmosphere for TRAPPIST-1 c than the water models, making it impossible to constrain the surface pressure of a CO₂-dominated atmosphere from the mass and radius (and stellar Fe/Si) alone. We ran a retrieval analysis with a different sampler, emcee (Foreman-Mackey et al. 2013), with a log-uniform prior for the surface pressure, and obtained a similar PDF for the CO₂-rich atmosphere.

We assumed the atmospheric parameters retrieved in our adaptive MCMC analysis and generated emission spectra with their respective temperature-pressure profiles, as explained in Sect. 5. Figure 5 shows the complete emission spectra and mean filter fluxes for TRAPPIST-1 c. Both our interior structure retrievals and evolution models (Krissansen-Totton et al. 2021) predict that TRAPPIST-1 c very likely has a bare surface. Therefore, we also considered the possibility of a bare surface on TRAPPIST-1 c, given the high probability obtained in our MCMC analysis for a volatile mass fraction equal to zero. Hu et al. (2012) obtained the emission spectra of bare terrestrial surfaces for different minerals. We estimated the brightness temperature for the irradiation conditions of TRAPPIST-1 c from the results of Hu et al. (2012) for two minerals. These minerals are a metal-rich surface and a granitoid one since they are the two surfaces with the highest and lowest emission for the same irradiation conditions, respectively. We approximated the emission spectrum of these surfaces to that of a black body with a temperature equal to the estimated brightness temperature. These two brightness temperatures are 443 K and 346 K (black and grey in Fig. 5, respectively). For very low emission fluxes (<200ppm), TRAPPIST-1 c would present a CO₂-dominated atmosphere. For fluxes greater than 550 ppm in the F1500W filter, TRAPPIST-1 c would have no atmosphere and an emission that corresponds to the bare surface with the lowest albedo.

Zieba et al. (2023) find that the emission of TRAPPIST-1 c in the F1500W filter is 421 ± 94 ppm. This value is within the range 200 and 500 ppm, which presents a degeneracy between bare surfaces with high albedos, such as granitoid, feldspathic, and clay (Hu et al. 2012), and thin atmospheres with little or no CO₂. Our water model shows an emission flux slightly above the CO₂ model in the 15 μm band. This may be due to the presence of clouds, which is considered in our ID atmospheric model. Nonetheless, self-consistent 3D Global Circulation Models (GCMs) of water-dominated envelopes with surface pressures consistent with our 1er estimate show a higher emission flux (Lincowski et al. 2023, Turbet et al., in prep.). In this case, O₂– and H₂O-dominated atmospheres could be compatible with the measured emission of TRAPPIST-1 c (presenting a degeneracy between the two) and with a bare surface with a high albedo. To break this degeneracy, follow-up photometric observations are needed at longer wavelengths (i.e. filters F1800W and F2100W) to distinguish between a bare rock and O₂ and H₂O atmospheres. We observe that for the near-infrared filters (λ = 5.60 to 11.30 μm), the models have very similar fluxes that are compatible within the uncertainties, which makes distinguishing between the different scenarios in these wavelengths impossible.

Greene et al. (2023) and Ih et al. (2023) find that the emission flux of TRAPPIST-1 b is mostly consistent with a bare surface. This has implications for the Fe bulk content of all planets in this planetary system. In our analysis of TRAPPIST-1 c without a constraint on the Fe/Si mole ratio, we obtain that its CMF = 0.26 ± 0.08. This value is in agreement within the uncertainties with the CMF derived for both planets b and c using the previous version of our modelling framework (Acuña et al. 2021). The uncertainties in this work are larger than those of the previous analysis due to the adaptive sampling in the MCMC (see Sect. 4). This CMF constraint from planet b breaks the degeneracy between the water and Fe content for the outer planets of TRAPPIST-1, supporting our WMF results for the outer planets of the system in scenario 2 of Acuña et al. (2021).

Fig. 4 PDF corner plot of the atmospheric parameters of TRAPPIST-1 c for two different atmospheric compositions. The atmospheric parameters include the surface temperature (Tsurf), atmospheric thickness (zatm), Bond albedo, and the surface pressure (Psurf). The two envelope compositions we consider are H2O (left panel) and CO2 (right panel).

Fig. 5 Simulated emission spectra of TRAPPIST-1 c with MIRI photometric filter mean fluxes for water and CO2 atmospheres. The spectra were obtained with the high-resolution, k-correlated MSEIRADTRAN model. We show spectra for pure (100%) water and CO2 compositions. Colour boxes indicate the wavelength ranges of the different filters.

6.2 55 Cancri e

The super-Earth 55 Cancri e is also in close orbit (P = 0.66 days) to a bright star. This exoplanet has had several interior and atmospheric hypotheses proposed. Madhusudhan (2012) explored a carbon-rich interior given the high C/O ratio found for the 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. However, a classical Fe-rich core and a silicate mantle are compatible with a volatile envelope rich in secondary atmosphere species. Furthermore, the absence of an H-/He-dominated envelope seems likely to be 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. The latter suggests that 55 Cancri e could have an optically thick atmosphere with a low mean molecular weight (Hammond & Pierrehumbert 2017). The possibility of a fully H₂O–dominated atmosphere was discarded since it would 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 CO₂, CO, and N₂, among other compounds.

Hu et al. (2021) propose observing 55 Cancri e in emission spectroscopy, combining the NIRCam F444W filter (3–5 µm) and MIRI LrS (5–14 µm). We summarize the 1σ confidence intervals of the interior and atmospheric parameters obtained in our retrievals for 55 Cancri e in Table 4. We observe that a water-dominated atmosphere reproduces the observed data well, with WMFs of up to 1.8%, indicating the possibility of a thick envelope with P_surf > 300 bar.

We did not model the scenario of a CO₂-dominated envelope for 55 Cancri e. The reason is twofold: (1) the CO₂ envelope is not extended enough to match the density of 55 Cancri e and yields a denser interior than the data show, and (2) at temperatures higher than 4000 K, CO, not CO₂, would be the dominant species in a C-rich atmosphere. This changes the emission of the atmosphere, as CO is a different absorber than CO₂. 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 CO₂ and would yield 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 would have a smaller thickness than a CO₂ atmosphere due to heavier molecular weights under similar atmospheric surface conditions. Adding silicate absorbers decreases the total planetary radius in H/He envelopes (Misener & Schlichting 2022). 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 matches the low planetary density for CMFs below 13% (Table 4), which is indicative of a bulk interior less dense and less Fe-rich than that of Earth (CMF_⊕ = 32%). An alternative bulk composition for 55 Cancri e in the volatile-poor scenario would be a carbon-rich mantle, as suggested by Madhusudhan (2012).

In most of the applications of our modelling framework, we considered a mesospheric temperature of T_top = 200 K. 55 Cancri e is extremely irradiated and likely has a higher mesospheric temperature. We repeated our analysis with a mesospheric temperature of T_top = 1000 K. We find no difference in the interior parameters compared to the low Ttop analysis. We considered the high top temperature case to generate the emission spectrum at a higher resolution. Figure 6 shows the complete predicted emission spectrum of 55 Cancri e from 3 to 14 µm. The water line at 4.3–4.4 µm can be seen. At wavelengths above 3.7 µm, the noise level is low (Fig. 6), which makes the spectral features in this range easy to identify with JWST in the high-molecular weight atmosphere scenario of 55 Cancri e. For comparison, in Fig. 6 we also show the black body emission at 1952 and 1785 K, which correspond to the irradiation temperatures typically associated with bare rocks (Hu et al. 2012).

Fig. 6 Predicted emission spectrum with the k-correlated, highresolution MSEIRADTRAN model for a water-rich atmosphere in 55 Cancri e with NIRCam and MIRI LRS. We assumed a mesosphere temperature of Ttop = 1000 K. We show the black body emission at different irradiation temperatures for comparison.

7 Discussion

In this section, we discuss the processes that were not included in our model but that 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_⊕) and a radius of ≃1.1 R_⊕, a dry mantle model and a wet magma ocean model estimate a WMF = 3 × 10⁻³ and 8 × 10⁻³, respectively. This difference decreases with a greater planet mass and a lower WMF. Including dissolution of silicates and a magma ocean would increase our WMF estimates by a factor of less than 2.7 (see Fig. 4 in Dorn & Lichtenberg 2021). Since our MCMC realisations span at least two orders of magnitude (from 10⁻⁶ to 10⁻⁴), 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 on 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_⊕, the difference is negligible, and it would yield 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 layers would form as well as their exact temperature values. We tested how changes in the mesospheric temperature may impact the OLR and the atmospheric thickness with our atmospheric model. We performed a test with a new mesospheric temperature of 1000 K and compared it to our default temperature of 200 K. The surface temperature at which the atmosphere is in radiative-convective equilibrium is similar in the two 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, such layers are more likely to form in atmospheres composed of H/He and silicates than pure water or CO₂ envelopes (Misener & Schlichting 2022). This means that our envelope mass fraction 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 Fig. 7) and not at the near-surface.

In our interior analysis, we considered a grey model for the calculation of the albedo (Sect. 3.4). The estimated grey opacity for water is based on that of Earth (Nakajima et al. 1992). Nonetheless, this parameter may be different for highly irradiated planets. Our interior–atmosphere models yield albedos between 0 and 0.30 in the water envelope case for TRAPPIST-1 c and 55 Cancri e. The change in surface temperature induced by this variation in the albedo is below 30 K. This difference in surface temperature produces changes in radius of less than 1% (Sect. 3.5). Planets with a magma ocean underneath their envelope may present hazes and aerosols (Kempton et al. 2023) whose albedo is higher than that of Earth-like water clouds. The effects of different cloud properties on interior modelling has been explored for gas giants by Poser (2019), who find that if optically thick clouds are high up in the upper atmosphere, they have a negligible effect on the inference of metal content (water or rock).

To calculate the density of CO₂, we used 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) carried out a comparison between the ideal EOS and non-ideal EOS for carbon dioxide. They find that CO₂ starts to behave as a non-ideal gas at ≃7 MPa, which corresponds to 70 bar. The non-ideal EOS yields a higher density than the ideal EOS (Böttcher et al. 2012). This means that the atmospheric thickness of a CO₂ envelope with a nonideal 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 CO₂-dominated atmosphere, which further supports that CO is more likely to constitute the atmosphere of this planet instead. For TRAPPIST-1 c, a higher density for the CO₂ envelope in comparison to the water envelope strengthens the degeneracy with the bare rock scenario. The implementation of a 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 focus 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 CO₂ atmospheres when given their observed mass, radius, and their host stellar abundances. We coupled the interior and the atmosphere self-consistently to obtain the boundary conditions at the top of a supercritical water layer or a silicate mantle given the irradiation conditions of a low-mass planet. This was done by calculating the bolometric emission flux and the Bond albedo in order to compute the flux emitted and absorbed by an atmosphere in radiative-convective equilibrium. We used a 1D k-correlated atmospheric model with updated opacity and EOS data for the computation of the bolometric emission within our interior-atmospheric model. We also demonstrated that using a constant step size when sampling the prior distribution in a MCMC scheme is not efficient for exploring the parameter space in interior modelling. This constant maximum step size 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 for rocky Earth-sized planets and super-Earths, whose WMFs are close to zero but nonetheless remain important for determining surface pressure.

Moreover, we used the surface pressure and temperature conditions retrieved with our interior–atmosphere model to generate emission spectra with our k-correlated atmospheric model, MSEIRADTRAN. We computed 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 characterisation proposals. The particular parameter that interior models can provide self-consistently for atmosphere models are the surface pressure and temperature, which are usually chosen arbitrarily in order to generate spectra. Emission spectra are more sensitive to the choice of surface temperature and thermal structure than transmission spectra.

We showcased 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 emission photometry and spectroscopy observations with JWST. We binned our emission spectra according to the response functions of the MIRI filters in order to predict emission fluxes for TRAPPIST-1 c in 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 that it lacks an atmosphere (WMF = 0). Nonetheless, the presence of a secondary atmosphere cannot be ruled out. In this scenario, TRAPPIST-1 c could have an H₂O-dominated atmosphere of up to 80 bar of surface pressure. The density of a CO₂- or O₂-rich envelope is not low enough to place any constraints on the surface pressure from mass and radius data alone. We have presented emission flux estimates for the filter centred at 15 µm, F1500W, which can be compared with observations (Zieba et al. 2023). Moreover, we discuss that a bare surface on TRAPPIST-1 b (Greene et al. 2023; Ih et al. 2023) has implications for the Fe content of all planets in the system. The Fe/Si mole ratios for planets b and c as retrieved in Acuña et al. (2021) and revisited in this work can be used to constrain the Fe contents of the other planets in the TRAPPIST-1 system. This breaks the degeneracy between the WMF and CMF for the outer planets in TRAPPIST-1, supporting the WMFs obtained in scenario 2 of Acuña et al. (2021). For 55 Cancri e, a massive water envelope with more than 300 bar of surface pressure is necessary to fit its low density and an Fe content similar to its stellar host or an Earth-like core simultaneously. We determined 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 µm. However, this part of the spectrum does not contain any spectral lines of water or CO₂, which are essential for determining the abundances in the envelope.

In our modelling approach, we have considered water- and CO₂-dominated atmospheres (99% H₂O and 1% CO₂ and vice versa). 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, according to observations and models (Madhusudhan et al. 2020; Bézard et al. 2022; 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 aim 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 in order to prepare proposals for JWST and future atmospheric characterisation facilities, such as Ariel (Tinetti et al. 2018). Our model can also be used within retrieval frameworks to simultaneously interpret mass, radius, and upcoming JWST emission spectral data to break degeneracies in exoplanet compositions.

Acknowledgements

M.D. and O.M. acknowledge support from CNES. We thank the anonymous referee whose comments helped to improve and clarify this manuscript. We acknowledge Emmanuel Marcq and Jérémy Leconte for their exchange on the k-correlated method. L.A. thanks Paul Mollière and Laura Kreidberg for useful discussions on emission spectra.

References

All Tables

All Figures

	Fig. 1 Dry adiabatic coefficient for water, κ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 correspond to κad = 0.35, 0.25, and 0.15, respectively.
In the text

Fig. 2 Comparison of the radius and interior–atmosphere boundary temperature between our fc-correlated model (MSEIRADTRAN) and that of Pluriel et al. (2019). Left panel: mass-radius relationships for a planet with a water-dominated atmosphere orbiting a Sun-like star at ad = 0.05 AU. Dashed lines indicate the total radius calculated by MSEIRADTRAN, while the solid line corresponds to the interior radius, which comprises the core, mantle, and supercritical water (SW) hydrosphere. 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.

In the text

	Fig. 3 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 with a dashed black line. The mantle mass fraction (MMF) is defined as MMF = 1 – CMF -WMF.
In the text

	Fig. 4 PDF corner plot of the atmospheric parameters of TRAPPIST-1 c for two different atmospheric compositions. The atmospheric parameters include the surface temperature (Tsurf), atmospheric thickness (zatm), Bond albedo, and the surface pressure (Psurf). The two envelope compositions we consider are H2O (left panel) and CO2 (right panel).
In the text

	Fig. 5 Simulated emission spectra of TRAPPIST-1 c with MIRI photometric filter mean fluxes for water and CO2 atmospheres. The spectra were obtained with the high-resolution, k-correlated MSEIRADTRAN model. We show spectra for pure (100%) water and CO2 compositions. Colour boxes indicate the wavelength ranges of the different filters.
In the text

	Fig. 6 Predicted emission spectrum with the k-correlated, highresolution MSEIRADTRAN model for a water-rich atmosphere in 55 Cancri e with NIRCam and MIRI LRS. We assumed a mesosphere temperature of Ttop = 1000 K. We show the black body emission at different irradiation temperatures for comparison.
In the text