© The Authors 2024

1 Introduction

Large-scale melting of the silicate mantles of terrestrial planets is common during the late phases of their growth (Stevenson 1981; Wood et al. 2006; Elkins-Tanton 2012; Chao et al. 2021). These episodes of magma ocean (MO) give rise to extensive chemical equilibration between the silicate reservoir of the planet and its exterior, allowing the outgassing of volatile elements that are originally present in its building blocks (Abe & Matsui 1988; Elkins-Tanton 2008), or the ingassing of nebular material (Olson & Sharp 2019). The variability in solubility exhibited by the different volatile species leads to different outgassing patterns for different elements (Elkins-Tanton 2008; Lebrun et al. 2013; Bower et al. 2019; Nikolaou et al. 2019; Lichtenberg et al. 2021), further altered by redox-sensitive speciation (Katyal et al. 2020; Bower et al. 2022; Gaillard et al. 2022). For modest atmospheric masses (≪1% of the planetary mass) that are typical of rocky planets, the bulk redox state of the MO-atmosphere system is set by core-mantle equilibration during the planetary differentiation and is quantified by the oxygen fugacity ($f_{{\mathrm{O}}_2}$). In the Solar System, the terrestrial MO likely outgassed a neutral to oxidized atmosphere (Hirschmann 2012; Sossi et al. 2020; Deng et al. 2020; Kuwahara et al. 2023) due to the significant $f_{{\mathrm{O}}_2}$ gradient between the MO bottom, where $f_{{\mathrm{O}}_2}$ is buffered by equilibration with Fe alloy, and its surface, which in turn buffers the redox state of the atmosphere. Conversely, smaller bodies have been theorized to outgas at more reducing conditions (Deng et al. 2020; Armstrong et al. 2019) because they lack such a pressure-induced negative gradient and because the $f_{{\mathrm{O}}_2}$ pressure-gradient is positive at shallow depths.

Beyond the Solar System, rocky exoplanets likely orbit most stars in the Milky Way, and the catalog of observations is growing. With the current limitations of rocky exoplanet observations, the radius of rocky exoplanets ranges to about 1.6 Earth radii (the lower bound of the so-called radius valley; Fulton et al. 2017), corresponding to approximately 6 Earth masses. The outgassed atmospheres during MO episodes on exoplanets probably exhibit a diversity that reflects their sizes, compositions and chemistries (e.g., redox state). The fate of highly volatile elements that are essential for life (H–C–N–O–S) is strongly influenced by the MO outgassing because once they are delivered to the atmosphere, these elements are prone to escaping to space, especially during the pre-main sequence of the host star, when it emits a strong X-UV flux that can drive hydrodynamic escape (Lammer et al. 2018). Understanding MO outgassing is thus key to any theory of planetary habitability. Furthermore, atmospheric characterization by ongoing and upcoming missions will soon give access to the chemical composition of atmospheres of terrestrial exoplanets. The origin of these atmospheres can be diverse (nebular accretion, volcanic outgassing, delivery of volatile-rich cometary materials, etc.; Krisanssen-Totton et al. 2020; Kite & Barnett 2020), and MO outgassing provides a potentially important source mechanism that is particularly relevant for young systems. Therefore, observations need to be accompanied by the most advanced models in order to understand the planetary context of the studied exoplanets. Several studies have investigated the composition of MO-outgassed volatile atmospheres relevant to early Earth or Mars (Elkins-Tanton 2008; Lebrun et al. 2013; Salvador et al. 2017; Nikolaou et al. 2019; Lichtenberg et al. 2021), as well as the effect of the redox state of the planet (Katyal et al. 2020; Bower et al. 2022; Gaillard et al. 2022), but a systematic study of the planetary parameters relevant for the diversity of exoplanets, in particular, linking these parameters with the oxidation state of the atmosphere, is still needed.

While some outgassing trends can be derived from simple physical principles, the chemical speciation of the atmosphere can significantly complicate the picture. Hence, in this work, we use an MO outgassing chemical model to calculate the composition of outgassed H–C–N–O–S atmospheres on MO exoplanets as a function of their mass and of the extent of their MO (more specifically, the molten silicate mass fraction, i.e., the ratio of the mass of the MO to the total mass of silicate, which we denote as Φ), and of their redox state. We characterize the extent of the outgassing of each element and the gaseous speciation in the atmosphere. These results are used to infer the possible gas-depletion patterns of exoplanets (i.e., the relative depletion in one volatile element compared to another in the bulk planet) and to predict observables in terms of alterations in the mass–radius relations for MO exoplanets.

2 Model

We first determine the mass–radius relations of a set of planets with various masses and Φ. This allows us to calculate the mass of the MO, its bottom pressure, and its profile (all necessary quantities for the chemical model). We then compute the surface $f_{{\mathrm{O}}_2}$ based on the core-mantle equilibrium $f_{{\mathrm{O}}_2}$ ($f_{{\mathrm{O}}_2,\mathrm{e}\mathrm{q}}$) at the bottom of the MO, which we vary as another free parameter. Finally, we use a chemical equilibrium model coupled to volatile mass conservation to calculate the atmospheric pressure and composition of these planets.

2.1 Mass-radius relation of refractory reservoirs

There exists a wealth of literature on rocky exoplanets mass–radius relations (Valencia et al. 2006; Dorn et al. 2015; Zeng et al. 2016, 2019; Otegi et al. 2020; Agol et al. 2021; Unterborn et al. 2016, 2023). However, to the authors’ knowledge, only Dorn & Lichtenberg (2022) took into consideration the density difference between solid and (hot) liquid silicate in their mass–radius relations. We considered a set of planets of masses (M_p) between 0.1 and 6 Earth masses with a molten silicate mass fraction (Φ) between 0.1 and 1. For each set of parameter values, we used burnman (Myhill et al. 2021) to generate a three-layer planet composed of pure fcc iron (using the corresponding equation of state of Saxena & Erikson 2015), a bridgmanite solid lower mantle (using the corresponding equation of state from de Koker et al. 2013), and a liquid MgSiO₃ MO (using the corresponding equation of state from de Koker et al. 2013). While the smallest planets (approximate mass of Mars) are unlikely to have stable bridgmanite, and the melting curves we used to compute the depth of the MO are matched for a natural peridotite composition, we kept the same equation of state in order to isolate the effect of the outgassed atmosphere on the mass–radius relations discussed below. For all these planets, we adjusted the core-mantle boundary, MO bottom (when different), the planetary radii (R_CMB, R_MO, and R_p, respectively), and the potential temperature of the MO (T_MO), to recover (1) the desired total mass and molten silicate mass fraction, (2) an Earth-like iron core mass fraction (0.32), and (3) to match the solidus with the MO adiabat at the MO bottom. We used the solidus from Fiquet et al. (2010) for Earth- and sub-Earth-sized planets and that from Stixrude (2014) for super-Earths. Criteria 1-3 are fulfilled with a precision of 1%. All reservoirs were assumed to be adiabatic, and the thermal boundary layers between them were ignored. The resulting geometries are indicated in Table B.1, and the temperature profiles are represented in Fig. 1a.

2.2 MO redox state

The process resulting in an MO (e.g., giant impact) is generally associated with an event of alloy-silicate equilibration (the alloy originating from the core of the impactor in the case of a giant impact). This equilibrium buffers the redox state of the MO at a certain value of $f_{{\mathrm{O}}_2,\mathrm{e}\mathrm{q}}$ depending on the compositions of the alloy and silicate melts (Rubie et al. 2015). This value, expressed as ΔIW_eq, that is, the difference in log-units from the iron-wüstite (IW) buffer at the same P-T conditions, was treated as a free parameter and was varied between 0 and −5, roughly encompassing values estimated for the Solar System (Rubie et al. 2015). $f_{{\mathrm{O}}_2,\mathrm{e}\mathrm{q}}$ sets the Fe³⁺/∑Fe ratio in the MO, where the equilibrium between ferrous and ferric iron oxides is the redox buffer after scavenging of Fe alloy in the core (Hirschmann 2012, 2022; Armstrong et al. 2019; Deng et al. 2020; Sossi et al. 2020). We calculated this ratio considering that our chosen $f_{{\mathrm{O}}_2,\mathrm{e}\mathrm{q}}$ was reached at the bottom of the MO where we assumed that the alloy-silicate equilibrium took place. From the obtained Fe³⁺ /∑Fe (assumed to be homogeneous in the MO) and the MO P-T profile (Fig. 1a), we computed the $f_{{\mathrm{O}}_2}$ profile through the MO using the parameterization from Hirschmann (2022) and Deng et al. (2020) (represented in Fig. 1b). The $f_{{\mathrm{O}}_2}$ reached at the surface is the value at which the atmosphere is buffered.

It has been shown that the slope of the $f_{{\mathrm{O}}_2}$ profile in the MO changes sign between 5 and 10 GPa, that is, the MO becomes oxidized with increasing depth at shallow levels (P <~ 6 GPa; O’Neill et al. 2006; Zhang et al. 2017) and then starts to be reduced with increasing depth at deeper levels (P >~ 6 GPa; Hirschmann 2012; Deng et al. 2020; Sossi et al. 2020; Kuwahara et al. 2023). A consequence of this is that as the planet becomes larger and/or the MO becomes deeper, beyond 5 and 10 GPa, the surface of the MO becomes increasingly oxidized (Fig. 1). None of the existing parameterizations is calibrated at pressures higher than 55 GPa (Deng et al. 2020). However, the pressure gradient of the $f_{{\mathrm{O}}_2}$ may cease to exist at large depths (Fig. 1b), and we therefore considered all MOs with a bottom deeper than 55 GPa to have a surface $f_{{\mathrm{O}}_2}$ that corresponds to the highest value reached for MOs shallower than 55 GPa, which is ~1.5 log-units above the bottom value ΔIW_eq.

We kept a constant core-mass ratio while varying the redox state at which the core and the MO equilibrate. However, these two quantities are not independent: a larger core generally implies a more reduced mantle (Rubie et al. 2015). Nevertheless, a straightforward relation only holds for a single-stage core-mantle equilibration, which is a highly unrealistic scenario. In addition, other processes (e.g., an evolution in the composition of the accreting material) can induce a deviation in the redox state of the MO from that imposed by equilibration with the bulk core. Therefore, bearing this caveat in mind, we chose to vary the redox state of the MO as a free parameter while keeping a constant, Earth-like core mass ratio, similar to that which is most frequently represented in rocky planets of the Solar System.

Fig. 1 Vertical profiles of temperature (a) and $f_{{\mathrm{O}}_2}$ (noted as the log-unit difference from the IW buffer), (b) through the MO for all Mp and Φ, for ΔIWeq = −2. The thick dashed and dotted black lines in panel a correspond to the solidus from Fiquet et al. (2010; used for Earth and sub-Earth-mass planets), and the solidus from Stixrude (2014; used for super-Earth-mass planets), respectively. The solidus defines the temperature reached at the bottom of the MO. The thin dotted black line corresponds to the maximum pressure in panel b. The solid vertical line in panel b corresponds to ΔIWeq = −2, which is the value reached at the bottom of the MO. The profile of ΔIW for cases with ΔIWeq = 0 (ΔIWeq = −5) are offset by 2 (−3) log-units on the x-axis of panel b.

2.3 MO atmosphere equilibrium

The chemical model computes the partial pressures of the various gaseous species in the atmosphere and the elemental solubilities based on the mass of the MO, the radius, and the gravity at the surface of the planet, a set of T–dependent equilibrium constants (evaluated at the potential temperature of the MO), and the bulk masses of H, C, N, and S present in the MO-atmosphere system. It accounts for mass conservation of all four volatile elements, present under 11 different gaseous species in the atmosphere (H₂, H₂O, CO, CO₂, CH₄, N₂, NH₃, HCN, SH₂, S₂, and SO₂), dissolution equilibria between gas and silicate melt, and the seven following gas-gas redox equilibria: $${\mathrm{H}}_2+1/2{\mathrm{O}}_2={\mathrm{H}}_2\mathrm{O}$$(1) $$\mathrm{C}\mathrm{O}+1/2{\mathrm{O}}_2={\mathrm{C}\mathrm{O}}_2$$(2) $${\mathrm{C}\mathrm{H}}_4+{\mathrm{O}}_2=\mathrm{C}\mathrm{O}+2{\mathrm{H}}_2\mathrm{O}$$(3) $${\mathrm{N}}_2+3{\mathrm{H}}_2=2{\mathrm{N}\mathrm{H}}_3$$(4) $${\mathrm{H}}_2\mathrm{O}+\mathrm{H}\mathrm{C}\mathrm{N}=\mathrm{C}\mathrm{O}+{\mathrm{N}\mathrm{H}}_3$$(5) $${\mathrm{S}}_2+2{\mathrm{O}}_2=2{\mathrm{S}\mathrm{O}}_2$$(6) $$2{\mathrm{H}}_2\mathrm{S}={\mathrm{S}}_2+2{\mathrm{H}}_2\text{. }$$(7)

The system of equations that is solved consists of the seven chemical equilibria above, and the mass conservation of H, C, N and S. We note that O mass conservation is omitted and replaced with the buffering of $f_{{\mathrm{O}}_2}$ by the MO, as described in the previous section. The general form of the chemical equilibria equations is $${\displaystyle \prod _s}{f}_{\mathrm{s}}^{c_{\mathrm{s}}}=K(T),$$(8)

where f_s is the fugacity of the gaseous species s at the surface, c_s is the signed stoichiometric coefficient (positive for products, negative for reactants) of the gaseous species s in the equilibrium under consideration, K is the equilibrium constant of the equilibrium under consideration, and T is the potential surface of the MO (assumed to coincide with the surface temperature). The fugacities are calculated as the partial pressures divided by a reference pressure of one bar. The equilibrium constants of Eqs. (1)–(7) are calculated from each species’ T-dependent formation Gibbs-free energy tabulated from the JANAF tables (Chase et al. 1985).

The general form of the chemical equilibrium equations (derived in Appendix C) is $$M_{\mathrm{e}}={\displaystyle \sum _s}\left[(M_{\mathrm{M}\mathrm{O}}{S}_s(p_{\mathrm{s}},P,f_{{\mathrm{O}}_2})+\frac{4\pi R_{\mathrm{p}}^2}{g}{p}_{\mathrm{s}}\frac{{\mu }_{\mathrm{s}}}{{\mu }_{\mathrm{a}\mathrm{t}\mathrm{m}}})\frac{{\mu }_e{\lambda }_{\mathrm{s}}^e}{{\mu }_{\mathrm{s}}}\right],$$(9)

where M_e is the bulk mass of element e, M_MO is the mass of the MO; S_s, p_s, and μ_s are the solubility law, partial pressure at the surface, and molecular weight of species s, respectively; P and $f_{{\mathrm{O}}_2}$ are the total pressure and the oxygen fugacity at the surface; R_p and g are the planetary radius and surface gravity, respectively; μ_atm is the average molecular weight in the atmosphere; μ_e is the atomic weight of element e; and ${\lambda }_{\mathrm{s}}^e$ is the number of atoms of e in one molecule of s. In Eq. (9), the first term in the innermost brackets corresponds to the mass of species s dissolved in the MO, and the second term corresponds to the mass present in the atmosphere. The term in the outermost bracket converts the mass of species s into the mass of element e, and the total is summed over all e-bearing species. We benchmarked our out-gassing model against other similar models from the literature (Bower et al. 2022; Gaillard et al. 2022) for the cases presented in Bower et al. (2022; see Appendix D).

In order to maintain a manageable parameter space, we restricted our study to a fiducial case with bulk contents of 100 ppm of the MO mass for each element. The volatile delivery during planet formation is not well understood, and the observed fractionation between these elements could be inherited from the building blocks (e.g. Piani et al. 2020; Grewal et al. 2021) or might postdate the MO phase (Albarede 2009). In addition, H ingassing from the protoplanetary disks (for MOs on planets that are still embedded in the primordial nebula) could also alter the H budget, with an increasing effect of planetary mass (Olson & Sharp 2019; Young et al. 2023). Stökl et al. (2016) calculated that in a solar nebula-like disk, a Mars-sized planet could attract an H₂ envelope of ~ 100 ppm of its bulk mass by the time of disk dispersal (1-10 Myr), which corresponds to our fiducial value for H. Larger planets are likely to accrete more nebular gas (Stökl et al. 2016 found that super-Earths are likely to trigger a gravitational runaway phase, but their model provided an upper bound on accretion because it neglected erosion processes). We therefore extended our parameter space to a higher H budget (1000 ppm) in order to investigate this possibility. Finally, as rocky bodies in the Solar System exhibit a consistent N depletion compared to C and H (Marty 2012; Bergin et al. 2015; Dasgupta & Grewal 2019), we also tested cases in which the N budget was decreased to 10 ppm.

2.4 Volatile solubilities

Unless stated otherwise, the solubilities were set by gas-melt equilibria. The corresponding species-specific solubility laws, their references, and their calibration domains are reported in Table A.1: We used Henry law coefficients from Hirschmann et al. (2012) for H₂ and those of Hirschmann (2016) for CO, CO₂, and for CH₄ (C solubility under “oxidized”, “reduced”, and “very reduced conditions” in their Table 2, respectively), and the recent parameterization for water solubility in peridotitic melts by Sossi et al. (2023). We neglected dissolution of HCN, NH₃, SH₂, and SO₂ and used the $f_{{\mathrm{O}}_2}$-dependent solubility laws from Dasgupta et al. (2022) for N₂ and from Gaillard et al. (2022) for S₂. Hence, the elemental solubilities of all elements are sensitive to $f_{{\mathrm{O}}_2}$: H and C because of the different solubilities of their redox end-member carriers, and N and S through a direct $f_{{\mathrm{O}}_2}$-sensitive solubility. It is important to note that these various parameterizations have been calibrated by laboratory experiments in a parameter space that sometimes does not extend to actual MO-like conditions. In particular, the upper bound in temperature at which these experiments are performed corresponds to the lower range of the temperatures envisaged here (see Table B.1). Hence, high-temperature experiments, although challenging, are needed to calibrate models such as the one used here, and to constrain the influence of temperature on volatile solubility. Melt composition is also a limitation: Most experiments (except for Sossi et al. 2023) were performed on basaltic composition, which is relevant for volcanic outgassing, but is not accurate for bulk rocky reservoirs.

It is possible that a given volatile element, rather than being either dissolved in the MO or outgassed to the atmosphere, is also incorporated in a condensed phase. Based on the solubilities, this possibility is highest for graphite/diamond precipitation. To account for the possibility of graphite precipitation, we calculated the CO₂ and CO concentrations at graphite saturation based on the P − T − $f_{{\mathrm{O}}_2}$ profiles in the MO (Fig. 1) using the parameterizations from Eguchi & Dasgupta (2018) and Yoshioka et al. (2019), respectively. Assuming a well-mixed MO, the concentration of each species is constant throughout the MO, and graphite will precipitate if the C concentration exceeds the minimum value of the respective saturation profiles. If this occured, we decreased the bulk C mass by 1% and computed the new chemical equilibrium. We iterated until the CO₂ and CO concentrations in the MO were both lower than or equal to their minimum saturation values in the MO. The CO–CO₂ content of the MO was thus buffered at the graphite saturation at the bottom of the MO due to the decreasing solubility of these species with pressure (at equilibrium with graphite). As these parameterizations have only been calibrated up to a few GPa and ~1500 K and yield virtually zero solubility when extrapolated at higher pressures, we only tested the shallowest MO cases (Φ = 0.1, M_p ≤ 0.5 Earth mass).

3 Results

3.1 Physics of outgassing

The influence of the planetary mass and molten silicate mass fraction (Φ) on outgassing is best understood for a simplified case, where a single volatile element is present as a single gaseous species, obeying a Henrian solubility law (i.e., the molar concentration of volatile in the MO is proportional to the partial pressure of the species, which itself is equal to the total atmospheric pressure P). In a single-species atmosphere, considering a Henrian solubility law, the ingassed and outgassed masses are $$M_{\mathrm{i}\mathrm{n}}=M_{\mathrm{M}\mathrm{O}}\alpha P,$$(10) $$M_{\text{out }}=\frac{4\pi R_{\mathrm{p}}^2}{g}P=\frac{4\pi R_{\mathrm{p}}^4}{G{M}_{\mathrm{p}}}P\text{, }$$(11)

where G is the universal gravity constant, and α is the Henrian solubility constant. For a fixed Φ, we have M_MO ∝ M_p, and the out-to-ingassed mass ratio is therefore proportional to $R_{\mathrm{p}}^4/M_{\mathrm{p}}^2$, which is a decreasing function of R_p for any realistic mass–radius relation (e.g., for a homogeneous sphere where $M_{\mathrm{p}}\propto R_{\mathrm{p}}^3$, or also for more realistic cases: $M_{\mathrm{p}}\propto R_{\mathrm{p}}^{1/0.27}$ (Valencia et al. 2006) or $M_{\mathrm{p}}\propto R_{\mathrm{p}}^{1/0.29}$ (Otegi et al. 2020).

In order to ensure mass conservation, P increases as M_p increases (at a fixed Φ and bulk volatile concentration). Summing Eqs. (10) and (11), we can write $$P=\frac{M_{\text{bulk }}}{M_{\text{MO}}\alpha +\frac{4\pi R_p^2}{g}},$$(12)

where M_bulk = M_in + M_out is the bulk mass of volatiles in the MO+atmosphere system. Taking the limit of Eq. (12) for low M_p yields $P\sim M_{\text{bulk}}g/\left(4\pi R_{\mathrm{p}}^2\right)$, which is the expression for the pressure of an atmosphere of mass M_bulk, thus corresponding to a fully outgassed case. Conversely, taking the limit of Eq. (12) for high M_p yields: P ~ M_bulk/(αM_MO), that is, αP ~ M_bulk/M_mo, which is the expression of the solubility (αP) when volatiles are completely ingassed (replacing M_in with M_bulk in Eq. (10)).

If there are multiple species in the atmosphere, Eq. (12) becomes $$p=\frac{M_{\text{bulk }}}{M_{\mathrm{M}\mathrm{O}}\alpha +\frac{4\pi R_{\mathrm{p}}^2\mu /{\mu }_{\mathrm{a}\mathrm{t}\mathrm{m}}}{g}},$$(13)

where p is the partial pressure of the species under consideration, μ is its molecular mass, and μ_atm is the mean molecular mass of the atmosphere. If μ_atm decreases (e.g., because a light species is added), then p decreases as well. This it does not mean that the species is being ingassed, however: as p decreases, the ingassed mass (M_MOαp) decreases as well. On the contrary, the positive effect of the decrease of μ_atm on the outgassed mass dominates the negative effect of the decrease of p: the species is being outgassed although its partial pressure decreases. The opposite occurs when μ_atm, and addition of heavy species leads to net ingassing.

Finally, keeping planetary parameters (mass, gravity, and radius) constant, the outgassed reservoir only varies with P, while the ingassed reservoir varies as Φ × P. The ingassed-to-outgassed mass fraction thus varies as Φ. Again, the atmospheric pressure has to increase so that the in- and outgassed mass fractions add up to one.

We verified these trends by monitoring the outgassed mass fraction of each element upon varying the planetary mass (Fig. 2). C-bearing species have very low solubilities compared to most other species (except for the most reduced cases, where graphite precipitates). They remain mostly in the atmosphere, although their outgassed fraction decreases slightly at high planetary masses and molten silicate mass fraction. In the most reduced cases (for 0.1 and 0.5 M_E with Φ = 0.1; Fig. 2a) graphite precipitates, significantly decreasing the outgassed mass fraction. While we only accounted for graphite saturation in the smallest planets with shallow MOs, this probably also applies to other reduced cases, and our result could overestimate the amount of outgassed carbon in other reduced cases. Unfortunately, the conditions for graphite precipitation at higher P − T remain unconstrained. N solubility follows a similar trend (although slightly less pronounced). However, at a low planetary mass, it features a decrease that is not explained by these physical considerations but by the $f_{{\mathrm{O}}_2}$–dependence on the N solubility (see Sec. 3.2). Because the elemental solubility of H is offset by the very high solubility of water, its outgassed mass fraction is always much smaller. However, it also follows a clear trend of a decreasing outgassed mass fraction with increasing planetary mass. The water solubility does not follow a Henrian behavior, but still induces the same trend, as observed by Sossi et al. (2023). Finally, while this trend is also followed by S for shallow and/or reduced MOs (Fig. 2a), the decreasing solubility of S with increasing $f_{{\mathrm{O}}_2}$ (similar to N) reverses this trend as ΔIW_eq and Φ increase. In general, decreasing outgassed fractions with increasing planetary mass suggest that the volatile budget of super-Earths is much less prone to loss due to MO outgassing and subsequent atmospheric escape than smaller planets, in particular, for H (Sossi et al. 2023). Conversely, small planets such as Mars outgas a large fraction of their H, which is then prone to escaping due to the low binding energy of these planets (see Sec. 4.2).

Fig. 2 Outgassed mass fraction of all volatiles as a function of planetary mass for a shallow MO (a), a half-mantle MO (b), and a whole-mantle MO (c). The curve colors and style correspond to element and ΔIWeq, respectively.

3.2 Chemistry of outgassing

In this section, we show how variations in the redox conditions provide another important control of the extent of outgassing and the composition of the outgassed atmosphere for a variety of rocky planets in terms of mass. The case of S already illustrated how the decreasing solubility with increasing $f_{{\mathrm{O}}_2}$ could dominate the physical controls on outgassing for the deepest MOs and yield an increase in the outgassed mass fraction of S (from ~0.2 weight percent, wt%, of the bulk S budget for Mars-sized planets up to 50 wt% for 6 Earth masses super Earths; Fig. 2c and supplementary materials).

Figure 3 represents the atmospheric composition over the complete parameter space. Two trends are clearly observed: first, the large fraction of H₂ at low planetary mass and/or low Φ. This effect, already hinted at in the previous section, is due to the low molecular mass of H₂, which amplifies the influence of the planetary mass on its outgassing, in particular, at low M_p. Second, at constant ΔIW_eq, the mixing ratios of the oxidized species systematically increase with increasing planetary mass and Φ. In MO deeper than 6 GPa (i.e., in this study, all but the shallowest MO on the smallest planet), an increasing depth of the MO causes its surface to become more oxidized (Hirschmann 2012; Armstrong et al. 2019; Deng et al. 2020; Kuwahara et al. 2023), thereby affecting the relative abundances of species in the same redox system. This second effect accentuates the first effect to promote H₂ in the atmospheres of small planets and shallow MOs. Over the whole M_p − Φ space, the surface ΔIW varies by 3 log-units for each value of ΔIW_eq (see Fig. 1b and supplementary materials), leading to large variations in the atmospheric speciation.

The element H is present in the most gaseous species (6 out of the 11). As a major species in the atmosphere, however, it is essentially carried by H₂ and H₂O, and it is marginally carried by CH₄. The H₂−H₂O system is dominated by the effect of the low molecular mass of H₂ already described and by the high (and non-Henrian) solubility of water. The T-dependence of equilibrium (1) should promote H₂O at low temperature (the equilibrium constant increases by ~12 orders of magnitude between the highest and the lowest MO surface temperature investigated here), but it is the opposite of what we observe. This is due to the decrease in surface $f_{{\mathrm{O}}_2}$ between these same cases (Fig. 1b), by more than 12 orders of magnitude, the decrease in surface ΔIW adds up to the T–dependence on the IW buffer, which dominates the variation of the equilibrium constant. From this equilibrium, for an Earth-sized planet with Φ ~ 0.5 (corresponding to an MO surface temperature of 2265 K), the fugacities of these two species coincide around IW. At higher oxygen fugacity, water dominates and allows for a very large ingassing of H. At lower oxygen fugacity (see supplementary materials), abundant H₂ in the atmosphere allows for the marginal presence (subpercent level) of other hydrogenated species, such as CH₄, HCN, NH₃,and SH₂.

CO is always a major gas in the atmosphere. The importance of CO₂ increases with ΔIW_eq, balancing CO around a surface $f_{{\mathrm{O}}_2}$ of IW+1 (for an Earth-sized planet and Φ ~ 0.5). CO₂ has a higher solubility in silicate melts than CO, thus altering out-gassing of C with increasing $f_{{\mathrm{O}}_2}$ (i.e., increasing M_p, Φ, and/or ΔIW_eq). Under reducing conditions, the high H₂ fugacity allows for some C to be transferred to CH₄ and HCN, but always as minor or trace species. As for the H system, the increase in the equilibrium constant of Eq. (2) with decreasing temperature is balanced by that of $f_{{\mathrm{O}}_2}$, preventing promotion of CO₂ over CO on small planets with a colder MO surface. While graphite saturation in the most reduced cases significantly alters C outgassing, the resulting atmospheric composition is only slightly altered compared to the same cases when this process is not accounted for (although the C/N ratio in the atmosphere is impacted; see Fig. 4). Other species (in particular, H₂) are indirectly affected. They maintain comparable molar fractions while the total pressure decreases when C outgassing is suppressed by graphite precipitation.

N is consistently present as N₂, but other N carriers such as NH₃ and HCN are only marginally present. While the N solubility strongly depends on $f_{{\mathrm{O}}_2}$ and becomes highly soluble under reducing conditions, the N ingassed mass fraction increases only under the most reduced conditions (ΔIW_eq = −5, Φ ≤ 0.5 and M_p ≤ 1 Earth mass; Figs. 2a and b). Over the remaining parameter space, the physical effect of ingassing with increasing planetary mass and molten silicate mass fraction dominates the decreasing solubility of N. This is consistent with the results of Gaillard et al. (2022), who found N solubility to increase significantly around IW-4. This value is reached (at the surface) only in the most extreme cases of our parameter space (i.e., Mars-size planets, or planets with a shallow MO). The N redox system is not as T-sensitive as the H and C ones.

As already pointed out, due to its high solubility at low $f_{{\mathrm{O}}_2}$ (Gaillard et al. 2022), S is mostly absent in the atmosphere except in the most oxidized cases, in which it is present as the oxidized end-member SO₂, one of the major species of our oxidized super-Earth cases. Here, the increase in the equilibrium (7) constant with increasing T acts in the same direction as the increase in the $f_{{\mathrm{O}}_2}$ on large planets, promoting speciation of S as SO₂ and thereby decreasing its effective solubility. In more reduced cases, the presence of H promotes SH₂ over S₂, which hardly exceeds the ppm level in the atmosphere. S can also precipitate as sulfide under reduced conditions (Blanchard et al. 2021). We tested for it in the same way as we tested for graphite saturation (but in the whole parameter space, sulfide saturation parameterization being calibrated over wider P-T conditions), and did not find sulfide precipitation because the S concentration we considered was moderate.

Over the parameter space investigated, lava planets whose oxidation states have been buffered by an event of core formation exhibit a diversity of atmospheric composition, dominated by H for small planets and by C (reduced CO or oxidized CO₂) for super-Earths. Nitrogen (as N₂) is consistently present, while the presence of of S in the atmosphere is diagnostic of oxidizing conditions. In general, oxidized atmospheres exhibit a larger diversity of major species than reduced atmospheres, which are dominated either by CO and H₂ or by H₂ alone.

Fig. 3 Atmospheric composition (expressed as the molar fraction of each gas in the atmosphere) for exoplanets from 0.1 to 6 Earth masses (ME, columns) and an MO mass from 0.1 to 1 times the planetary silicate mass (Φ, rows) for all three $f_{{\mathrm{O}}_2,\mathrm{e}\mathrm{q}}$, investigated (expressed as ΔIW). The numbers indicate the total surface pressure in bar (the number in parentheses indicates the total surface pressure for cases without graphite precipitation). All cases have bulk silicate-equivalent abundances of 100 ppm of H, C, N, and S. The amount of some species (e.g., NH3, HCN, SH2, or S2) is always too small to be noticeable in these charts.

Fig. 4 Atmospheric C/N (top row), C/H (middle row), and C/S (bottom row) normalized to the bulk value for all cases with Φ = 0.1 (left column), Φ = 0.5 (middle column), and Φ = 1 (right column). The blue, green, and orange lines correspond to the fiducial, H-rich, and N-depleted cases, respectively. The line shades correspond to the MO redox state (more transparent lines correspond to more oxidized cases, ranging from ΔIWeq = −5 to 0). The dashed lines (only in the left column) correspond to cases when graphite saturation was ignored. A bulk-similar atmosphere is plotted along y = 1, a C-enriched atmosphere (compared to N, H or S) above this line, and an or S-enriched atmosphere (compared to C) is plotted below. The y-axis is linear in the top row and logarithmic in the middle and bottom rows.

3.3 Influence of the bulk volatile budget

The previous results apply for the fiducial case of planets with contents of 100 ppm per mass of MO of each volatile element (H, C, N, and S). We investigated how the observed trends are affected when these budgets were varied. It is extremely difficut to vary the initial budget in volatile elements because it depends on a wide range of processes, from parent body processing in the planetary building blocks to its accretion history. To narrow our parameter space, we focused on two alternative scenarios: 1) an H-enriched case (1000 ppm H), accounting for the possibility of nebular gas capture during the lifetime of the protoplanetary disk (in particular, by massive super-Earths), and 2) an N-depleted case (10 ppm N), as this pattern is a general feature of rocky and icy bodies in the Solar System (Bergin et al. 2015).

Outgassing couplings between the various elements (and their bulk contents) chiefly occur via an alteration of the mean molecular mass of the atmosphere, which modifies the mass distribution between in the MO and the atmosphere. Other couplings, chemical in origin, could in theory involve species bearing more than one element (not counting O), such as CH₄, NH₃, HCN or SH₂. However, the redox and temperature conditions and the high solubility of sulfide strongly disfavor the first three and the last one, respectively.

Figure 4 represents three elemental ratios (C/N, C/H, and C/S) in the atmosphere, normalized by their bulk value in the MO + atmosphere system. In N-poor cases, N depletion is more pronounced in the outgassed atmosphere than in the bulk MO + atmosphere system (compare the orange and blue lines in the top row of Fig. 4). In the most reduced cases, the enhanced solubility of N accentuates this effect, which is made obvious when graphite precipitation in the shallowest MOs on the smallest planets is ignored (dashed lines in the top left panel of Fig. 4. The bulk elemental solubility of N increases by a factor 2 between the fiducial case and the N-depleted case, which we attribute to the decrease in mean molecular mass of the atmosphere caused by shortage of N₂ (see supplementary materials), and which results in ingassing. These trends affect the N system more strongly due to its decreased abundance, and a divergence from the fiducial case on other elemental ratios is hardly noticeable. The nonmonotonic change in the atmospheric-to-bulk C/N ratio for Φ = 1 in the most reduced case (amplified in the N-depleted case) is due to nonlinearities in the chemical system (e.g., the link between surface $f_{{\mathrm{O}}_2}$ and planet size). In the most oxidized cases, C is increasingly speciated as CO₂, thereby increasing its effective solubility and decreasing the atmosphere-to-bulk C/N ratio (see the top right panel in Fig. 4).

Interestingly, increasing the H budget suppresses graphite precipitation in the most reduced cases, enforcing a high atmospheric C/N ratio. The N-depleted case is due to the decrease in the mean molecular mass of the atmosphere, this time not due to depletion of heavy species (N₂), but to the addition of light ones (H₂). Addition of H also significantly reduces atmospheric C/H and C/S, which causes both ratios to become closer to the bulk value (graphite precipitation can decrease the atmospheric C/H below the bulk value in the most reduced cases; see the middle left panel in Fig. 4). Here again, the addition of H lowers the mean molecular mass in the atmosphere (which always remains below 9 amu in the H-rich cases, but reaches 26 amu in the other cases; see supplementary materials), thereby promoting outgassing of all species and thus causing the atmosphere to become more similar in bulk.

4 Discussion

We use the results presented above to derive updated mass–radius relations that take the complex chemistry of an outgassed atmosphere into account. We then discuss the fate of these outgassed atmospheres because potential atmospheric escape can alter the chances of observations as well as the volatile content patterns on exoplanets.

4.1 Mass–radius relations

Mass-radius relations are crucial theoretical tools for constraining the properties of exoplanets. They are used to infer the internal structure of a planet based on its bulk density. Mass-radius relations for rocky exoplanets have mostly been studied in the light of their refractory composition (e.g. Valencia et al. 2006; Dorn et al. 2015; Unterborn et al. 2016, 2023; Zeng et al. 2016, 2019; Agol et al. 2021; Acuña et al. 2023). Exoplanet observations by the transit method measure the transit radius, which, if considered as the actual surface radius of the solid or condensed part of the planet, can lead to a significant error on the planetary structure. Recently, Turbet et al. (2020) showed that the presence of water as steam rather than in a condensed phase could significantly alter the mass–radius relations of terrestrial planets, which are more irradiated than the runaway greenhouse limit, yielding a much larger transit radii. The resulting thick steam atmosphere maintains a high surface temperature, and when the planet is hot enough to melt its silicate envelope, Dorn & Lichtenberg (2022) showed that the dissolution of water in the MO could further alter the influence of the steam atmosphere on the mass–radius relation and decrease the transit radius.

While previous studies have focused on steam atmospheres, we showed that outgassing of endogeneous volatiles from an MO can follow very different chemical patterns depending on the planet size, the molten silicate mass fraction, and the redox state (Fig. 3). In order to assess the influence of these outgassed atmospheres on the mass–radius relations for terrestrial exoplanets, we calculated their transit radii, shown in Figs. 5a–c (for the fiducial volatile contents). The transit radii were calculated as a function of wavelength using petitRADTRANS (Molliere et al. 2019) and were averaged between 0.5 and 1 μm, which is similar to the wavelength range of most photometric missions, e.g., the CHaracterising ExOPlanets Satellite (CHEOPS), the Transiting Exoplanets Survey Satellite (TESS), the Kepler satellite, or the upcoming PLAnetary Transits and Oscillations of stars satellite (PLATO). As observed in Sec. 3.1, small planets exhibit the strongest relative outgassing, in particular of H₂, yielding a large atmospheric vertical extension and thus an observed transit radius that is up to twice as large as the actual planetary radius (Fig. 5). For example, a 0.1 Earth-mass planet could appear to be as large as an Earth-sized planet (Figs. 5a and d). Interestingly, graphite precipitation on a planet like this, although decreasing the surface pressure, also decreases the mean molar mass of the atmosphere, yielding larger vertical extension (see Table B.1).

Changing the volatile budgets by either decreasing N or increasing H also affects the mean molecular mass of the atmosphere. While the former only has a negligible effect, the latter (increasing the H budget by one order of magnitude) leads to significantly inflated super-Earths, whose bulk densities are similar to a 100% H₂O composition. No realistic transit radii for H-rich reduced small planets could be obtained because the assumption of hydrostatic equilibrium breaks down in these conditions of surface pressure, temperature, and gravity. Furthermore, we note that our mass–radius relations for the refractory part of the planet yield similar densities for whole-mantle MO planets (Φ = 1) and 100 % (solid) MgSiO₃ planets from Zeng et al. (2016; Fig. 5c), which is partly due to the presence of an MO (less dense than cold mantle), and partly due to potential differences between the equations of state from burnman and those used in Zeng et al. (2016). While this is likely an underestimate of the bulk density, it is still dominated by the effect of the atmospheres for small planets and for H-enriched super-Earths.

The very low apparent density of small reduced exoplanets (an overestimate by 100% of the radius corresponding to an eightfold underestimate of the density) is a striking feature that might be mistaken for a nebular accreted primordial atmosphere, however. Nevertheless, as an MO-causing giant impact can occur much later than disk dissipation (Raymond et al. 2004), and considering the likely limited survival time of H₂-rich atmospheres around small planets, an outgassed origin would provide an explanation in case a small, low-density planet were observed around a relatively young star (a few 100 Myr old). The transit depth (estimated as (R_tr/R_st)², where R_st is the stellar radius) of the most reduced 0.1 M_E planets is 80 ppm around a Sunlike star. This is similar to the detection limit of the upcoming PLATO mission (PLATO Team 2017). Lupu et al. (2014) and Bonati et al. (2019) made the case for the search of transient lava worlds in young systems, which our result helps to understand.

Fig. 5 Mass–radius relations of bare silicate core radii (dashed lines) and with the outgassed atmosphere as a function of ΔIWeq (colored contours) for Φ = 0.1 (a), Φ = 0.5 (b), and Φ = 1 (c), and relative sizes (to scale) of the solid, liquid, and gaseous part of the planets for the fiducial volatile contents and ΔIWeq = −5.

4.2 Fate of secondary atmospheres

After they are outgassed, volatiles are prone to atmospheric escape. H₂-rich inflated atmospheres are particularly sensitive to hydrodynamic escape, where heating of the upper atmosphere results in a pressure gradient that destabilizes hydrostatic equilibrium and induces an outward flow. The gas is advected outside of the gravitational zone of influence of the planet (the Bondi or Hill sphere), and is lost. An upper limit for the mass-loss rate is provided by the energy-limited regime, where the entirety of the stellar extreme-UV (XUV) flux received by the planet is used to lift the atmosphere (Watson et al. 1981), $$F_{{\mathrm{H}}_2}=\frac{\pi R_{\mathrm{X}\mathrm{U}\mathrm{V}}^2{f}_{\mathrm{X}\mathrm{U}\mathrm{V}}ϵR_{\mathrm{p}}}{G{M}_{\mathrm{p}}},$$(14)

where $F_{{\mathrm{H}}_2}$ is the H₂ escape flow (in kg s⁻¹), R_XUV is the XUV absorption radius, which we took to be equal to the transit radius calculated in the previous section, f_XUV is the incoming XUV flux (in W m⁻²), ϵ is an efficiency coefficient, which we took to be equal to 0.3 following Katyal et al. (2020), and G is the universal gravitational constant. Assuming energy-limited H₂ escape, Fig. 6 represents the lifetime of an outgassed atmosphere as a function the planetary mass and the incoming XUV flux for the most reduced cases (ΔIW = −5) and Φ = 0.1. While Mars-sized planets lose their H₂-dominated atmospheres in less than 1 Myr under most instellation conditions for young stars, H₂-dominated secondary atmospheres on larger planets that also exhibit anomalously large transit radii can be retained over several million years for early-Earth-like XUV instellations. Furthermore, the presence of heavier species (already at the percent level) reduces the efficiency of hydrodynamic H loss by several orders of magnitude on Mars- (Yoshida & Huramoto 2020) and Earth-sized planets (Yoshida & Kuramoto 2021), provided that these species are infrared emitters. In particular, Yoshida & Huramoto (2020) showed that an H₂ atmosphere of 6 bar on Mars with 7% CO–CH₄ (comparable to our most reduced case) can be retained for ~8 Myr due to the radiative cooling of C-bearing species that compete with XUV heating and decrease the escape efficiency.

Escape of H₂-rich atmospheres can carry diluted heavier species and trigger loss of elements other than H (Kite & Barnett 2020). This occurs for species lighter than the crossover mass (Hunten et al. 1987), which is defined as $${\mu }_{\text{crossover }}={\mu }_{{\mathrm{H}}_2}+\frac{kT{F}_{{\mathrm{H}}_2}}{bg{X}_{{\mathrm{H}}_2}}\text{, }$$(15)

with k the Boltzmann constant, T the temperature in the upper atmosphere (set to the stratospheric temperature, 200 K), b = 2.2 × 10²¹ m⁻¹.s⁻¹ the diffusion parameter (Hunten et al. 1987), and $X_{{\mathrm{H}}_2}$ the molar mixing ratio of H₂ in the atmosphere. In Fig. 6, we also indicate by contours where the crossover mass matches the molecular mass of H₂O, CO, or N₂, and CO₂. Above these lines, the respective species escapes along with H₂, and hydrodynamic escape affects the bulk volatile content of their constitutive elements. However, we stress again that most of these species (except for N₂) are strong infrared emitters, can cool the upper atmosphere, and can significantly decrease the escape rate. Furthermore, it is questionable whether the most massive super-Earths, whose atmosphere is not dominated by H₂ (Fig. 3, third last row), are in the hydrodynamic escape regime. Only an accurate modeling of hydrodynamic escape of MO atmospheres, which is beyond the scope of this study, can be conclusive on the escape of the various volatile elements during this phase.

Massive loss of H₂ can also alter the redox state of the planet because H loss corresponds to a net oxidation (Sharp et al. 2013; Wordsworth et al. 2018; Kite & Schaefer 2021). In this case, a new $f_{{\mathrm{O}}_2}$ forcing sets in, offsetting the FeO–Fe₂O₃ buffer toward a higher $f_{{\mathrm{O}}_2}$. Sharp et al. (2013) calculated that the loss of 1.3 Earth water-ocean-equivalent of H (corresponding to 455 ppm H₂O of the BSE mass, or 50 ppm H) would raise the $f_{{\mathrm{O}}_2}$ of Earth’s mantle from IW to FMQ, which is an increase of ~4 log-units. Assuming that we can extrapolate the trend to lower $f_{{\mathrm{O}}_2}$ values, this corresponds to a crossing of the complete $f_{{\mathrm{O}}_2}$ space we investigated upon losing half of the H endowment. This could raise the initially most reduced cases to the most oxidized ones as a result of H escape. With increasing $f_{{\mathrm{O}}_2}$, H₂ is converted into H₂O, which is heavier and dissolves more readily in the magma, thus shutting down atmosphere escape. It is therefore unlikely that even the smallest reduced exoplanets become fully dessicated.

Fig. 6 Lifetime (in Myr) of the MO-outgassed atmosphere as a function of planetary mass (x-axis), and incoming XUV flux (LXUV, y-axis), normalized by the solar XUV flux received by the Earth 4.5 Ga ($L_{\mathrm{X}\mathrm{U}\mathrm{V}}^{\mathrm{s}\mathrm{u}\mathrm{n}}$(t = 4.5 Ga)) from Ribas et al. (2005), for the cases most prone to hydrodynamic escape (ΔIW = −5 and Φ = 0.1). The loci in which the crossover mass matches the molecular mass of H2O (black), CO or N2 (purple), and CO2 (white) are also represented.

5 Conclusions

By investigating the atmospheric composition and extent of MO outgassing on exoplanets from 0.1 to 6 Earth masses, with molten silicate mass fractions from 0.1 (shallow MO) to 1 (whole-mantle MO), and with core-mantle equilibration $f_{{\mathrm{O}}_2}$ between IW-5 and IW (yielding to MO-atmosphere eaquilibrium between IW-5.5 and IW+1.5), we reached the conclusions that we list below:

• 1.

  The extent of outgassing decreases with increasing planetary mass and MO depth. The low molecular mass of H₂ amplifies this mechanism, leading to strong outgassing of H (as H₂) on small planets and shallow MOs;

• 2.

  The MO and atmosphere oxidation increases with increasing planetary mass and MO depth, increasing H and C retention in MO, but decreasing N and S retention;

• 3.

  C and N are mostly in the atmosphere, while H and S are mostly dissolved, except under the most oxidized conditions, where SO₂ is outgassed, and the most reduced conditions, where atmospheres become consistently H₂ rich. Graphite saturation in the most reduced cases can also drastically decrease C outgassing;

• 4.

  Extensive H outgassing on reduced small planets leads to a high vertical extension of their atmospheres, significantly affecting their mass–radius relations. This effect is amplified by limited C outgassing due to graphite saturation;

• 5.

  Extensive H outgassing on small planets makes them prone to H₂ escape, while volatiles are efficiently retained in the oxidized atmospheres with high molecular mass of super Earths (in addition to their increased gravity).

These results should inform future atmosphere-retrieval models applied to potential young lava worlds. In particular, several of the species predicted by our model, exhibiting strong absorption features in the infrared, offer an interesting target for spectroscopic studies.

Acknowledgements

The authors thank the anonymous reviewer for their constructive comments that helped improve the manuscript. M.M. thanks Tatsuya Yoshida for insightful discussions on hydrodynamic escape. This work received funding from the NASA grant 80NSSC18K0828.

Appendix A Solubility laws

Table A.1 lists the solubilitiy parameterizations we used, as well as their temperature, pressure, and $f_{{\mathrm{O}}_2}$ calibration ranges, when available. In addition to the experiments presented in these studies, the parameterizations we derived rely on a collection of data from the literature. The ranges indicated in Table A.1 reflect this as accurately as possible. We refer to the studies cited in these references for more information.

Appendix B Mass-radius relations

Table B.1 lists the mass–radius relations for the silicate cores of planets that were calculated as described in Section 4.1, and the transit radii and atmospheric pressures.

Appendix C Derivation of the mass-conservation equation

Let e be an element of atomic mass μ_e whose mass is conserved, and s be a gaseous species of molecular mass μ_s, composed of (among others) ${\lambda }_{\mathrm{s}}^e$ atoms of e. If p_s is the partial pressure of s in the atmosphere and P and μ_atm are the total atmospheric pressure and mean molecular mass, respectively, then the molar fraction of s in the atmosphere is p_s/P and its mass fraction is p_s/P × μ_s/μ_atm. Since the total mass of the atmosphere (under the assumption of hydrostatic equilibrium) is its the mass of s in the atmosphere is $4\pi R_{\mathrm{p}}^2/g$, then the mass of s in the atmosphere is $$M_s^{\mathrm{a}\mathrm{t}\mathrm{m}}=\frac{4\pi R_{\mathrm{p}}^2}{g}{p}_{\mathrm{s}}\frac{{\mu }_{\mathrm{s}}}{{\mu }_{\mathrm{a}\mathrm{t}\mathrm{m}}}.$$(C.1)

Let S_s be the solubility law of species s (in general, a function of p_s, P, and $f_{{\mathrm{O}}_2}$. Assuming that the atmosphere and the MO are in equilibrium (see Salvador & Samuel (2023) for a discussion of the validity of this assumption)), the mass of s dissolved in the MO is simply $$M_s^{\mathrm{M}\mathrm{O}}=M_{\mathrm{M}\mathrm{O}}{S}_{\mathrm{s}}(p_{\mathrm{s}},P,f_{{\mathrm{O}}_2}).$$(C.2)

The total mass of s in the system is the sum of C.2 and C.1, and the mass of e carried by molecules of s is the mass of s multiplied by $\left({\mu }_{\mathrm{e}}/{\mu }_{\mathrm{s}}\right){\lambda }_{\mathrm{s}}^e$. Summing over all species (whose masses are not conserved because of chemical reactions), we obtain Equation 9).

Appendix D Benchmark of the chemical equilibrium model

The chemical equilibrium model we used relies on a number of parameterizations for gas solubilities, gas-gas equilibria, and oxygen fugacity buffers. These choices, in addition to numerical implementation, contribute to potential discrepancies with other models. In order to assess these discrepancies, we reproduced the cases presented in Figures 4 and 8 of Bower et al. (2022) using our model as well as the model of Gaillard et al. (2022) (see supplementary materials). These cases (corresponding to a whole-mantle MO on an Earth-sized planet) cover a wide range of $f_{{\mathrm{O}}_2}$ (from IW-2 to IW+4 at the surface) and of H and C bulk contents (from 39 to 389 ppm, and from 0 to 1944 ppm, respectively). Bower et al. (2022) considered a C-O-H system, and we therefore used a negligible bulk N and S content (10⁻⁴ ppm). We only reproduced the initial state, as the final state (also presented in these Figures) has a volatile element content that is affected by MO crystallization, which is not reported.

We find that the results of all three models agree in terms of orders of magnitude and C/H outgassing. We note a systematic offset in the reduced-to-oxidized species ratio (e.g., H₂/H₂O for the H system) between our model and that of Gaillard et al. (2022) on one hand and the model Bower et al. (2022) on the other hand, the latter having higher ratios at the same $f_{{\mathrm{O}}_2}$. This discrepancy is likely due to the difference in the IW buffer used by Bower et al. (2022), from O’Neill & Eggins (2002), which differs by 0.5 log units from the model we used. The model of Gaillard et al. (2022) only converged for ΔIW < 0. Importantly, while this exercise illustrates agreements and discrepancies between models, the results are still only as good as the extrapolations of the parameterizations in use outside of the calibration ranges (see Table A.1).

References

All Tables

All Figures

Fig. 1 Vertical profiles of temperature (a) and $f_{{\mathrm{O}}_2}$ (noted as the log-unit difference from the IW buffer), (b) through the MO for all Mp and Φ, for ΔIWeq = −2. The thick dashed and dotted black lines in panel a correspond to the solidus from Fiquet et al. (2010; used for Earth and sub-Earth-mass planets), and the solidus from Stixrude (2014; used for super-Earth-mass planets), respectively. The solidus defines the temperature reached at the bottom of the MO. The thin dotted black line corresponds to the maximum pressure in panel b. The solid vertical line in panel b corresponds to ΔIWeq = −2, which is the value reached at the bottom of the MO. The profile of ΔIW for cases with ΔIWeq = 0 (ΔIWeq = −5) are offset by 2 (−3) log-units on the x-axis of panel b.

In the text

	Fig. 2 Outgassed mass fraction of all volatiles as a function of planetary mass for a shallow MO (a), a half-mantle MO (b), and a whole-mantle MO (c). The curve colors and style correspond to element and ΔIWeq, respectively.
In the text

Fig. 3 Atmospheric composition (expressed as the molar fraction of each gas in the atmosphere) for exoplanets from 0.1 to 6 Earth masses (ME, columns) and an MO mass from 0.1 to 1 times the planetary silicate mass (Φ, rows) for all three $f_{{\mathrm{O}}_2,\mathrm{e}\mathrm{q}}$, investigated (expressed as ΔIW). The numbers indicate the total surface pressure in bar (the number in parentheses indicates the total surface pressure for cases without graphite precipitation). All cases have bulk silicate-equivalent abundances of 100 ppm of H, C, N, and S. The amount of some species (e.g., NH3, HCN, SH2, or S2) is always too small to be noticeable in these charts.

In the text

Fig. 4 Atmospheric C/N (top row), C/H (middle row), and C/S (bottom row) normalized to the bulk value for all cases with Φ = 0.1 (left column), Φ = 0.5 (middle column), and Φ = 1 (right column). The blue, green, and orange lines correspond to the fiducial, H-rich, and N-depleted cases, respectively. The line shades correspond to the MO redox state (more transparent lines correspond to more oxidized cases, ranging from ΔIWeq = −5 to 0). The dashed lines (only in the left column) correspond to cases when graphite saturation was ignored. A bulk-similar atmosphere is plotted along y = 1, a C-enriched atmosphere (compared to N, H or S) above this line, and an or S-enriched atmosphere (compared to C) is plotted below. The y-axis is linear in the top row and logarithmic in the middle and bottom rows.

In the text

	Fig. 5 Mass–radius relations of bare silicate core radii (dashed lines) and with the outgassed atmosphere as a function of ΔIWeq (colored contours) for Φ = 0.1 (a), Φ = 0.5 (b), and Φ = 1 (c), and relative sizes (to scale) of the solid, liquid, and gaseous part of the planets for the fiducial volatile contents and ΔIWeq = −5.
In the text

Fig. 6 Lifetime (in Myr) of the MO-outgassed atmosphere as a function of planetary mass (x-axis), and incoming XUV flux (LXUV, y-axis), normalized by the solar XUV flux received by the Earth 4.5 Ga ($L_{\mathrm{X}\mathrm{U}\mathrm{V}}^{\mathrm{s}\mathrm{u}\mathrm{n}}$(t = 4.5 Ga)) from Ribas et al. (2005), for the cases most prone to hydrodynamic escape (ΔIW = −5 and Φ = 0.1). The loci in which the crossover mass matches the molecular mass of H2O (black), CO or N2 (purple), and CO2 (white) are also represented.

In the text