© The Authors 2023

1 Introduction

There are more than 350 ultrahot exoplanets (with equilibrium temperature of > 1500 K) discovered as of July 2022. About half of them have radii of <5 Earth radii and have masses that are consistent with either pure rocky planets or sub-Neptunes. At such temperatures and assuming a silicate-based mantle, their surfaces should be partially or entirely covered by a magma ocean. Such planets are of particular interest for future spectroscopic observations because convection in the mantle should be sufficiently rapid (Solomatov 2000) such that their atmospheres are readily replenished by chemical exchange with the underlying molten surface. As such, any spectroscopic constraints on the nature of their atmospheres may provide insights into the composition of their molten mantles. However, the search for atmospheres on ultrahot rocky worlds has so far brought mostly negative results. Spectroscopic observation of LHS 3844b (with an equilibrium temperature ~1000 K) yielded results that are compatible with the absence of an atmosphere (Kreidberg et al. 2019). The detection of an atmosphere around 55 CNC e is debated. Although a hotspot offset has been reported in the planet’s phase curve (Demory et al. 2016) - which can be interpreted as the consequence of a thick and opaque atmosphere –, detailed spectroscopic observations of 55 CNC e have not unambiguously detected any gaseous species (Bourrier et al. 2018). Some hints of the possible presence of HCN were reported in Tsiaras et al. (2016), but this would seem unlikely on chemical equilibrium grounds given that hydrogen itself has never been detected (see e.g., Ehrenreich et al. 2012).

Recently, Zieba et al. (2022) reported the possible indirect detection of an atmosphere around the ultrahot planet (T_eq ~ 2000 K) K2-141b by comparing Spitzer to synthetic spectra. However, the resolution of Spitzer spectra is too low to permit any clear identification of the molecular species in the partial pressures expected above such a rocky planet. Therefore, the detection of any sort of atmosphere around a lava planet remains elusive, despite growing international research efforts towards this goal. For the sake of completeness, we mention that a low-mean-density atmosphere containing HCN and CH₄ has been detected around the rocky world GJ 1132, but it is a “warm” planet (T ~ 530 K), and may not have a surface magma ocean (Swain et al. 2021).

Interestingly, ultrahot rocky planets tend to have high average densities. For example, 55 CNC e has an estimated density of ~5.9 g cm⁻³ (Crida et al. 2018; Mercier et al. 2022). When comparing the density of exoplanets with their equilibrium temperatures (Fig. 1 bottom), strongly irradiated rocky planets with equilibrium temperatures of T_eq > 1700 K have, on average, a much higher density than those planets with T_eq below 1700 K. More specifically, the population of terrestrial planets with densities < 5 g/cm³ almost completely disappears at T_eq larger than 1700 K (Fig. 1), and all observed ultrahot rocky planets smaller than 3 Earth radii have densities ranging from 5 to 11 g/cm⁻³ (with only one remarkable exception of a planet with T_eq ~ 2400 K and ρ ~ 3 g/cm⁻³). Otegi et al. (2020) demonstrate that volatile-rich and rocky exoplanet populations may be separated on the basis of their density, with a transition around 3 g/cm⁻³. Ultrahot mini-Neptunes and terrestrial planets may have a thin- to moderately massive atmosphere (less than 50% of the planet’s mass) above a molten mantle and an iron-rich core, as sketched out in Fig. 2. The broad positive trend of density versus equilibrium temperature might reflect a decreasing mantle/core mass ratio with increasing surface temperature. This trend could either result from efficient iron condensation (relative to silicates) in the hottest region of the protoplanetary disk during the formation of first solids (see e.g., Pignatale et al. 2017; Johansen & Dorn 2022) or arise from collisional (O’Neill & Palme 2008; Denman et al. 2020) or evaporative (Young et al. 2019; Charnoz et al. 2021; Jäggi et al. 2021) mass loss. However, most terrestrial planets are thought to have moved inward during their formation due to efficient migration, and so highly irradiated magma ocean exoplanets could have formed further away at much lower temperatures and with much higher hydrogen budgets than they have today.

In the present paper, we wish to investigate the possibility that the mantles of rocky planets in a magma ocean state evaporate at high temperatures in the presence of H₂ gas. Depending on the efficiency of this process, it might lead to a decrease in the mantle mass with increasing temperature, and thus an increase in planet density, given that the iron-rich core remains isolated from the evaporation occurring at the surface. Although the general treatment is not new, and evaporation from a magma ocean has been studied under varying circumstances (see e.g., Schaefer & Fegley 2010; Lupu et al. 2014; Schaefer et al. 2012; Ito et al. 2015; Ito & Ikoma 2021; Charnoz et al. 2021; Jäggi et al. 2021). It was found that the vapor above a magma ocean with a composition like that of Earth’s mantle (also called bulk silicate Earth, or BSE) has, on an anhydrous basis (i.e., neglecting H, C, Cl, S, N, and other volatiles) very low pressure (≪1 bar in general) and is dominated by Na and K for T < 2500 K, and by SiO for T > 3000 K (Schaefer & Fegley 2010; Schaefer et al. 2012; Ito et al. 2015; Ito & Ikoma 2021). The escape of a mineral atmosphere was studied in detail by Ito & Ikoma (2021), and was found to be very inefficient. A mineral-based atmosphere released by a magma ocean is rich in sodium, potassium, oxygen, iron, magnesium, and silicon (Schaefer & Fegley 2010). Schaefer & Fegley (2004) show that the dominant metallic gas above a silicate melt should be SiO for T > 2700 K and Na for lower temperature. Ito & Ikoma (2021) show that this mineral vapor dissociates to simple atomic and ionic compounds in the upper atmosphere. Because Na and K are efficient atmospheric coolants (Ito & Ikoma 2021), most of the energy transferred to the planet through irradiation by the star is re-emitted to space, making thermal escape inefficient with resulting heating efficiency in the range of ~5 × 10⁻⁴ to ~5 × 10⁻³, compared to 0.14−0.4 for hydrogen-dominated atmosphere (Valencia et al. 2010). Ito & Ikoma (2021) therefore concluded that atmospheric escape would have only a marginal influence on the composition of ultrahot rocky planets covered by a magma ocean. Nonetheless, Ito & Ikoma (2021) note that, at 3000 K, most of the Na and K content of the planet’s mantle should escape during the lifetime of the star. However, as Na and K may only represent a very small proportion of the mantles of such planets (together they comprise <1 mol% in the BSE for example), even if the entire planetary budget of Na and K were to escape, their loss would be unlikely to significantly affect either the mass or the radius of the rocky planet. For the atmospheric escape of the magma to have a significant impact on the mass or radius of the planet, major components, such as O, Mg, or Si, must be released in the vapor in very high quantities and then lost to space.

Here, we revisit the idea that atmospheric escape in the presence of hydrogen can engender a significant loss of the mantle of a rocky planet. The specificity of our approach lies in our consideration of the presence of a hydrogen envelope above the magma ocean and in computation of the consequences this has on the atmosphere’s speciation (which is modeled as a mixture between the H captured from the protoplanetary disk and species evaporated from the magma ocean) and atmospheric loss. To date, several studies have investigated the dissolution of H₂ or H₂O in a silicate magma ocean, and the effect on the planet’s mass and radius atmospheric loss (see e.g., Kite et al. 2020; Kite & Schaefer 2021), but neither the change in atmospheric speciation nor the identities and partial pressures of the major rock-forming elements have been investigated.

To our knowledge, the only other attempt to characterize the effect of H on the composition of a hot rocky planet atmosphere was by Fegley et al. (2016), who computed (based on the thermodynamic properties of hydroxide- and halide gas species) the compositions of metal-bearing atmospheres above a BSE magma for various partial pressures of H₂O, but did not consider the case of a captured H₂ envelope above a magma ocean. Although the effect of H₂ has been examined in its influence on the escape of K from rocky planets (Erkaev et al. 2022), these authors neglect H solubility in the magma ocean and do not provide a self-consistent model of atmospheric speciation in equilibrium with the magma. Lichtenberg et al. (2021); Dorn & Lichtenberg (2021) considered the effect of dissolved H₂ and H₂O on the structure of an exoplanet, but did not compute the effect of H on the equilibrium chemistry of the atmosphere. On the other hand, Sossi et al. (2020b, 2023); Bower et al. (2022); Gaillard et al. (2022) determine the equilibrium speciation of an atmosphere in equilibrium with magma oceans for a range of H/C, fO₂, and H masses, and argue that oxidized conditions promote the dissolution of H₂O in the magma ocean. However, these studies did not consider the influence of H₂O or H₂ on the speciation of rock-forming elements in the gas phase.

Because H₂ is the dominant constituent of protoplanetary disks, it may be captured during the planet formation process, even around growing rocky planets (e.g., Olson & Sharp 2019). The ubiquity of H may therefore play a major role in influencing the chemistry of atmospheres degassing from rocky planets. We concentrate our efforts on examining the effect of H on the gaseous species present above silicate melts in the system Na-K-Mg-Al-Fe-Si-O-H. We find that the evaporation of metal-bearing species from a BSE-like magma ocean is strongly enhanced by the presence of hydrogen. It has been shown that a H₂-rich atmosphere can be lost rapidly for ultrashort period planets, in less than 10 Myr (see e.g., Lopez 2017; Erkaev et al. 2022). Therefore, the effect of a captured H₂ envelope on atmospheric chemistry and its subsequent loss merits furtherinvestigation. In particular, we find that even minor amounts of H (with respect to the planet’s mass) may have a strong influence on the atmosphere composition.

In the first part of the present paper, we aim to quantify how the presence of H modifies the vapor composition compared to a H-free case using a thermodynamic model of vapor-magma ocean equilibrium. To illustrate the physical processes at play, in Sect. 2 we present a simplified model based on two reactions that show how H modifies liquid-vapor equilibria in promoting the vaporization of metal-bearing gases due to the decrease in oxygen fugacity (= O₂ partial pressure) it engenders, as well as the stabilization of H-bearing metal molecules in the gas phase. This leads to extensive depletion of Si and other metals from the silicate liquid. In Sect. 3, we consider a more realistic case of an infinite magma ocean with BSE composition in contact with a hydrogen envelope. We compute the vapor composition as a function of temperature and the amount of captured hydrogen. At the end of this section, we also show that measurements of Na/Si, Mg/Si, and Mg/Fe and the detection of H₂O in the atmosphere of ultrahot rocky planets may be used as indicators of the presence and nature of a magma ocean.

The second part of the paper explores how the presence of H promotes atmospheric escape of heavy metallic species by increasing their concentration in the atmosphere (SiO, O, O₂, Fe, Mg, etc.). Should this process prevail over long timescales, this may lead to considerable evaporation of mantle components (relative to the Fe-Ni-rich core), and thus give rise to high-density rocky planets. We use a simple energy-limited formalism modified to take into account the inefficient loss of mineral species as reported in Ito et al. (2015) and in Ito & Ikoma (2021) to quantify the extent of atmospheric loss.

Fig. 1 Density of known exoplanets versus their equilibrium temperatures. These data come from the online DACE database maintained by Observatory of Geneva (https://dace.unige.ch/) and are taken from the catalog of Otegi et al. (2020), which is an exoplanet catalog based on reliable, robust, and - to the greatest possible extent – accurate mass and radius measurements of transiting planets up to 120 earth masses (Otegi et al. 2020). Top: density vs. equilibrium temperature for all planets. Bottom: same as above but limited to planets with radii of smaller than 3 Earth radii. The dotted line shows a simple least-square fit to the data.

Fig. 2 Schematic illustration of the internal structure of a rocky planet as assumed in this paper, with an atmosphere containing a mixture of captured hydrogen and gaseous species derived from a molten mantle that forms a global magma ocean below the planet atmosphere. In the centre of the planet, an iron/nickel core is present. The composition of the magma ocean is assumed to be fixed, while the composition of the atmosphere is at thermodynamic equilibrium with that of the magma ocean.

2 Toy model: a SiO₂ magma ocean degassing in a H₂ atmosphere

We start with a simple chemical model to illustrate the effect of the presence of a primordial hydrogen envelope at the surface of a degassing magma ocean. For simplicity, we assume that the (infinite) magma ocean is made of SiO₂ only at temperature T and is surrounded by a H₂ atmosphere. To understand the basic ingredients, we reduce the physics to two equations (1) the vaporization of SiO_2(ℓ) (where the index (ℓ) stands for liquid and (g) for gas) in SiO_(g) and O₂, and (2) the conversion of H₂ and O₂ into H₂O. These reactions are $${\text{SiO}}_{2\left(\ell \right)}\iff {\text{SiO}}_{\left(\text{g}\right)}+1/{\text{2O}}_2,$$(1) $${\text{H}}_2\text{O}\iff {\text{H}}_2+1/{\text{2O}}_2.$$(2)

At equilibrium between the liquid and the gas, the partial pressure of SiO_(g) is given by the law of mass action at temperature T: $$P_{\text{SiO}}{}_{{}_{\left(\text{g}\right)}}{P}_{{\text{O}}_2}^{1/2}\text{=\hspace{0.17em}}a\left[{\text{SiO}}_{2\left(\ell \right)}\right]K_1\left(T\right),$$(3)

where $P_{{\text{SiO}}_{\left(\text{g}\right)}}$ is the partial pressure of SiO_(g), a[SiO_2(ℓ)] is the activity of SiO₂ in the melt (here we consider an ideal and pure liquid so that a[SiO_2(ℓ)] = 1 here) and $P_{{\text{O}}_2}$ is the partial pressure of O₂, that is the fO₂ (we consider an ideal gas, such that the fugacity, f, is equal to the partial pressure, p). Here, K₁ (T) is the equilibrium constant of this reaction at temperature T. The equilibrium constant is then related to the Gibbs-free-energy change of the evaporation reaction, K₁(T) = e^−ΔG1/RT where ΔG1 = ΔG⁰(SiO_(g)) + 1/2ΔG⁰(O₂) – ΔG⁰(SiO_2(ℓ)) and ΔG⁰ are the Gibbs free energies of formation of the different species from the elements. Similarly, for the H₂O formation equation (Eq. (2)): $$P_{{\text{O}}_2}^{1/2}{P}_{{\text{H}}_2}\text{=}{K}_2\left(T\right)P_{{\text{H}}_2\text{O}}.$$(4)

To close the system, we introduce the following limiting conditions. The first is that all oxygen atoms that go in the atmosphere come from the dissociation of SiO₂, such that the atomic O:Si ratio in the atmosphere is always equal to 2, which leads to the relation: $$P_{{\text{SiO}}_{\left(\text{g}\right)}}+P_{{\text{H}}_2\text{O}}+2P_{{\text{O}}_2}=2P_{{\text{SiO}}_{\left(\text{g}\right)}}.$$(5)

The above equation is established by considering that the total number of oxygen moles in the system is $N_O=N_{\text{SiO}}+2N_{{\text{O}}_2}+N_{{\text{H}}_2\text{O}}$, and that N_O = 2N_SiO, if we assume that all oxygen comes from the dissociation of SiO_2(ℓ). We therefore get the equality: $2N_{\text{SiO}}=N_{\text{SiO}}+2N_{{\text{O}}_2}+N_{{\text{H}}_2\text{O}}$. As partial pressures are proportional to the number of moles, one obtains $2P_{\text{SiO}}=P_{\text{SiO}}+2P_{{\text{O}}_2}+P_{{\text{H}}_2\text{O}}$. In addition, we assume that the total H content of the atmosphere is fixed and is equal to $P_{\text{H}}^0$, which is defined as: $$P_{\text{H}}^0=2P_{{\text{H}}_2\text{O}}+2P_{{\text{H}}_2},$$(6)

and is used as a proxy for the total quantity of H in the atmosphere in the remainder of the paper. Equations (3)–(6) form a close set of equations with unknowns ($P_{{\text{O}}_2}$, $P_{{\text{H}}_2\text{O}}$, $P_{{\text{SiO}}_{\left(\text{g}\right)}}$, $P_{{\text{H}}_2}$) and free parameters (T,$P_{\text{H}}^0$). We solve this simple system numerically.

The partial pressures of the various gas species as a function of the initial H content and temperature are displayed in Fig. 3. For low values of $P_{\text{H}}^0$ (<1 bar), the atmosphere is dominated by O₂ and SiO_(g) and corresponds to a pure rocky atmosphere in a hydrogen-free environment. As $P_{\text{H}}^0$ increases, H₂O begins to form, and $P_{{\text{O}}_2}$ drops sharply in response. As the system becomes increasingly oxygen-poor, the quantity of silicon in the atmosphere increases proportionally to $P_{{\text{O}}_2}^{-\left(1/2\right)}$, as given by Eq. (3). In Fig. 4, the hydrogen content of the atmosphere $P_{\text{H}}^0$ is set, while the temperature is varied. At low temperatures, H₂ predominates. However, as temperature increases, the equilibrium (2) shifts to the left, and H₂O becomes the prevailing H-bearing gas species above ~2700 K. Above ~3000 K, the H₂/H₂O ratio remains constant and SiO(g) and O₂(g) increase in tandem according to the $P_{{\text{O}}_2}^{-\left(1/2\right)}$ dependence.

We now consider two asymptotic cases. The partial pressure of SiO_(g) in the H-poor and H-dominated regimes can be easily computed as follows. If the pressure of hydrogen is extremely low, then H₂O and H₂ are mostly absent and the system reduces to the single reaction written in Eq. (1).

Therefore, the relation between the pressures of SiO_(g) and O₂ is given by Eq. (3). We set a[SiO_2(ℓ)] = 1 (pure and ideal liquid). The dissociation of SiO_2(ℓ) implies that there are always two atoms of O for 1 atom of Si in the gas. Counting the number of O and Si atoms in each species leads to the relation: $$2P_{{\text{O}}_2{}_{\left(\text{g}\right)}}+P_{{\text{SiO}}_{\left(\text{g}\right)}}=2P_{\text{SiO}}{}_{{}_{\left(\text{g}\right)}}.$$(7)

Combining Eqs. (3) and (7) and solving for $P_{{\text{SiO}}_{\left(\text{g}\right)}}$ leads to the simple solution: $$P_{{\text{SiO}}_{\left(\text{g}\right)}}\left(T\right)=2^{1/3}{K}_1{\left(T\right)}^{2/3}.$$(8)

This solution is displayed in Fig. 3 with the brown dashed line. In the opposite case, in which hydrogen becomes a major component, H₂ and H₂O control the system evolution. The system is therefore controlled by the reaction $${\text{SiO}}_{2\left(\ell \right)}+{\text{H}}_2\iff {\text{SiO}}_{\left(\text{g}\right)}+{\text{H}}_2\text{O.}$$(9)

The law of mass action gives the relationship between the partial pressures: $$P_{{\text{SiO}}_{\left(\text{g}\right)}}{P}_{{\text{H}}_2\text{O}}=K_3\left(T\right)P_{{\text{H}}_2}a\left[{\text{SiO}}_{2\left(\ell \right)}\right],$$(10)

where $$K3\left(T\right)=e^{-\frac{G^0\left({\text{SiO}}_{\left(\text{g}\right)}\right)+G^0\left({\text{H}}_2\text{O}\right)-G^0\left({\text{H}}_2\right)-G^0\left({\text{SiO}}_2{}_{\left(\ell \right)}\right)}{RT}}=K_1\left(T\right)/K_2\left(T\right).$$(11)

As all oxygen in the vapor comes from the evaporation of SiO2, the conservation of the oxygen atoms gives the relation $P_{{\text{H}}_2\text{O}}+P_{SiO}=2P_{SiO}$, which reduces to $P_{{\text{H}}_2\text{O}}=P_{SiO}$. The conservation of the hydrogen atoms gives $2P_{{\text{H}}_2\text{O}}+2P_{{\text{H}}_2}=P_{\text{H}}^0$. We set a[SiO_2(ℓ)] = 1 as usual. Solving for $P_{{\text{SiO}}_{\left(\text{g}\right)}}$ leads to a second-order polynomial whose only positive solution is $$P_{{\text{SiO}}_{\left(\text{g}\right)}}\left(P_{\text{H}}^0,T\right)=\frac{-K_3+{{\left(K_3(T\right)}^2+2K_3\left(T\right)P_{\text{H}}^0)}^{1/2}}{2}.$$(12)

This asymptotic solution is displayed as a purple dashed line in Fig. 3, and closely matches the SiO behavior for large hydrogen content. This also explains why SiO and H₂O have the same values at high H pressure. We see that at high hydrogen pressure, H₂O plays both the role of O and H reservoir, while O₂ is almost completely absent. This clearly demonstrates that, as the pressure of H increases, the content of Si in the atmosphere sharply increases also.

We conclude from this section that H plays a major in the composition of the vapor: because of the formation of H₂O, the O₂ partial pressure drops and, in turn, we observe an enhanced evaporation of SiO that scales with $P_{\text{H}}^{0^{1/2}}$. We now turn to the full model where this process acts on all metals contained in the magma ocean.

Fig. 3 Pressures of the different elements for the simple two-reaction model with only SiO(g), H2, H2O, and O2. Here, the temperature is fixed to 3000 K. SiO solution 1 and SiO solution 2 refer to asymptotic solutions with low and high hydrogen content, respectively.

Fig. 4 Pressures of the different elements for the simple two-reaction model with only SiO(g), H2, H2O, and O2. Here, the hydrogen content is fixed ($P_{\text{H}}^0=1\text{bar}$). SiO solution 1 and SiO solution 2 refer to asymptotic solutions with low and high hydrogen content, respectively.

3 Full model method calculation

Here, we describe a model of an infinite magma ocean (infinite = the liquid magma composition is fixed) in equilibrium with a H-rich primordial atmosphere. Our code is not unique; the MAGMA (Fegley & Cameron 1987) and VapoRock (Wolf et al. 2023) codes perform similar calculations. In particular, the VapoRock code takes into account nonideal interactions via a Regular solution model. For the vapor molecular composition, our approach uses a Gibbs free-energy minimization code (we use the iconic CEA/NASA code provided by Gordon & McBride 1996). We modify the CEA code to couple it to a vapor-liquid phase equilibrium model corresponding to liquid with a fixed composition. The same approach may be used with any atmospheric code (e.g., Woitke et al. 2018). We assume that the magma ocean is made of a nonideal mixture of liquid oxides, namely SiO_2(ℓ), MgO_(ℓ), FeO_(ℓ), Na₂O_(ℓ), and K₂O_(ℓ), with fixed molar abundances $X_i^l$ (where i stands for any of the liquid oxides). Each liquid oxide may vaporize into gaseous species through the limited list of reactions compiled in Table 1. The specificity of our approach compared to previous works is that we assume there is a pre-existing H envelope whose total pressure of monoatomic H is $P_{\text{H}}^0$. For now, we assume that all H stays in the vapor phase, whereas it is known that H₂ and H₂O may dissolve in the magma ocean (Hirschmann et al. 2012; Sossi et al. 2023). This simplification is not critical here as $P_{\text{H}}^0$ may be interpreted as the quantity of H contained in the atmosphere only.

3.1 The H monoatomic pressure: $P_{\text{H}}^0$

The H envelope reacts with gaseous species derived from the magma ocean to modify the atmospheric composition. Therefore, knowledge of the total amount of H in the atmosphere is useful. In order to quantify the amount of H independently of the atmosphere composition and planet mass, we introduce the “hydrogen monoatomic pressure”, $P_{\text{H}}^0$, which corresponds to the partial pressure of H if all the H of the atmosphere is in pure monoatomic form, and is computed as follows: $$P_{\text{H}}^0={\displaystyle \sum _i{v}_{\text{H}}^i{P}_i},$$(13)

where i is any molecule in the atmosphere, P_i is the partial pressure of molecule i, and $v_{\text{H}}^i$ is the stochiometric coefficient of H in molecule i. For example, in an atmosphere in which all H is in H₂, $P_{\text{H}}^0=2\times P_{{\text{H}}_2}$ where $P_{{\text{H}}_2}$ is the partial pressure of H₂.

3.2 Constitutive equations

For a given value of $P_{\text{H}}^0$ and magma temperature, the parameters we wish to determine are the molar fractions of each gaseous element in the atmosphere, namely O, Si, Mg, Fe, Na, H, and K, which we refer to collectively as $X_i^g$ (where i stands for any of the previous elements), as well as the total pressure (P_tot) of the atmosphere. For any combination of ($X_i^g$, P_tot), the CEA/NASA code can provide the molecular composition of the atmosphere at chemical equilibrium in the gas. To couple the atmospheric model to the infinite magma ocean model, we must find the only atmospheric composition that simultaneously satisfies (1) the liquid/vapor equilibrium, (2) the mass conservation of O, and (3) the mass conservation of H. The thermodynamic constraints are the following:

• Liquid/gas equilibrium states that the partial pressure of metal-bearing species evaporated from the melt (SiO, SiO₂, Na, K, Fe, Mg) depends on the atmospheric fO₂. Therefore, the following relations must be obeyed for the gas for each reaction j reported in Table 1. $$P_{{\text{vapour}}_j}=\frac{K_j\left(T\right)a\left({\text{liquid}}_j\right)}{P_{{\text{O}}_2}^{s_j}},$$(14)

  where vapour_j stands for the vapor species that bears the metal in evaporation reaction j (e.g., Na_(g) in reaction #3), liquidj is the corresponding liquid oxide (NaO_0.5 for reaction #3), s_j is the sto-chiometric coefficient (one-quarter in reaction #3), K_j(T) is the equilibrium constant of vaporisation reaction j, and a(liquid_j) is the activity of liquid #_j, that is a(liquid_j) = X_j × Γ_j, where X_j is the mole fraction of liquid_j in the magma ocean (assumed to be fixed) and Γ_j· is the activity coefficient of j in the liquid due to nonideal mixing effects. Activity coefficients, Γ_j, of all melt oxide components are interpolated using outputs of the VapoRock code (Wolf et al. 2023).

  When the atmospheric composition verifies all the relations for 1 < j < 6 given by Eq. (14), the gas is at thermodynamic equilibrium with the liquid.

• We assume that all metals and O in the atmosphere come from the vaporization of liquid oxides. This imposes a mass conservation relation between oxygen and all evaporated metals, which is governed by the stoichiometric relations of all reactions reported in Table 1. We simply get $$X_{\text{O}}^{\text{g}}=2X_{\text{Si}}^{\text{g}}+X_{\text{Mg}}^g+\frac{1}{2}{X}_{\text{Na}}^{\text{g}}+\frac{1}{2}{X}_{\text{K}}^{\text{g}}+X_{\text{Fe}}^{\text{g}},$$(15)

  where $X_i^g$ is the number of moles of atom i in the gas.

• The initial total content of H (measured as a pressure) is fixed to a constant value, $P_{\text{H}}^0$, meaning that Eq. (13) must also be verified.

  We note that we assume that Fe is only present as FeO in the melt, meaning Fe₂O₃ is ignored, as is reasonable to expect in a magma ocean in equilibrium with H₂ gas.

  The calculation is performed as follows:

  • Choose a magma composition ($X_i^l$), temperature, and a total hydrogen content ($P_{\text{H}}^0$). Here the magma composition is assumed to be BSE (Table 2).

  • Determine the atomic composition of the atmosphere ($X_i^{\text{g}}$ for all atoms i) and the total pressure P_tot. For all sets ($X_i^{\text{g}}$, P_tot), the corresponding molecular composition of the atmosphere is computed using the CEA/NASA chemical equil_gibrium code. The only solution to this problem is the ($X_i^{\text{g}}$, P_tot) set for which the molecular composition of the atmosphere simultaneously solves Eqs. (13)–(15).

  • The search is done using an iterative method, with a simplex minimization algorithm. The calculation goes on until Eqs. (13)–(15) are solved up to a relative accuracy of better than 10⁻³ for each partial pressure involved.

For the atoms considered in this paper (Si, Mg, K, Na, Fe, O, H), the following gas species are included in the CEA/NASA thermolib table (and are the same as those found in JANAF table): Fe, FeO, H, H₂, H₂O, H₂O₂, HO₂, K, K₂, K₂O, K₂O₂, K₂O₂H₂, KH, KNa, KO, KOH, Mg, Mg₂, MgH, MgO, MgOH, Na, Na₂, Na₂O, Na₂O₂, Na₂O₂H₂, NaH, NaO, NaOH, O, O₂, O₃, OH, Si, Si₂, Si₃, SiH, SiH₂, SiH₃, SiH₄, SiO, SiO₂.

3.3 Validation against other thermochemical codes

Appendix A displays the composition we find for the vapor in the absence of H. This is a case that has been extensively investigated in the past (see e.g., Visscher & Fegley 2013; Ito et al. 2015) and more recently in the VapoRock (Wolf et al. 2023) and LavAtmos codes (van Buchem et al. 2022). Our calculated partial pressures are in very good agreement with those calculated by existing codes (compare our Appendix to Fig. 4 of Wolf et al. 2023 and van Buchem et al. 2022).

3.4 Validation against molecular dynamics calculations

However, our calculations are limited to species that are included in the thermodynamic tables, such as JANAF, which may not be complete. For the sake of completeness, we compare our thermodynamic model against ab initio molecular dynamics simulations. For this, we employ the VASP package based on the planar augmented wavefunctions (Kresse & Hafner 1993) of the density functional theory. We consider pyrolite melt, with the bulk silicate Earth composition (McDonough & Sun 1995): 1/2Na₂O.2CaO.3/2Al₂O₃.4FeO.30MgO.24SiO₂. This six-component composition was used previously to model the crystallization (Caracas et al. 2019) and the structure of the magma ocean (Solomatova & Caracas 2019).

We performed simulations well inside the liquid-vapordome at 3000 K and a density of 0.4g/cm⁻³, at which neither liquid nor vapor are stable as a single phase, but coexist as a mechanical mixture. The simulations are realized within the generalized gradient approximation in the Perdew-Wang-Ernzernhofformu-lation (Perdew et al. 1996), spin-polarized, and corrected for the van der Waals interactions in the gas phase (Grimme et al. 2010).

We used five different initial configurations and ran simulations lasting between 50 and 136 picoseconds with a time step of 2 femtoseconds. For each simulation, we first build the entire interatomic connectivity matrix (Caracas et al. 2021), which allows us to define the chemical species. At 3000 K, we obtain a vapor whose main components are: SiO, FeO₂, Na, O₂, Mg, FeO, O, and SiO₂. The major difference with our thermodynamic calculation (and all other published calculations; see Appendix A) comes from the presence of FeO₂ observed in the ab initio simulation and not considered in any thermodynamic calculation. Moreover, the FeO₂ component is present at all temperatures in the simulations. This may indicate that further work is needed to reconcile thermochemical models and ab initio molecular dynamics calculations.

4 Results: atmospheric composition for a magma ocean degassing in hydrogen

4.1 Composition above the magma ocean

Our calculations are performed between T = 1800 K and 3500 K, and for total hydrogen pressures ranging from 10⁻⁶ to 10⁶ bar. At higher pressure, hydrogen could be metallic (>1 Mbar), which is beyond the scope of our paper, as we focus on terrestrial to sub-Neptune planets, and ignore giant planets. In addition, at such high pressures, the ideal-gas approximation on which the present paper is based fails. Examination of the virial coefficients for common gases indicate the assumption of ideality breaks down above ~10³ to 10⁴ bar. Refractory species such as Al and Ca may also be of interest but are not expected to comprise significant fractions of silicate atmospheres (e.g. Fegley et al. 2016); neither have they been detected in exoplanetary atmospheres. As such, we do not consider them in this work other than to define the mole fractions of melt oxides in the BSE.

Here, we present four fiducial calculations with magma oceans at 1800, 2000, 3000, and 3400 K and varying total H pressure. The total pressure of the atmosphere above the magma ocean is plotted in Fig. 5, the mass of hydrogen versus the planet mass is plotted in Fig. 6 and the compositions of the atmospheres at equilibrium with the magma oceans are plotted in Figs. 7–10. A rapid comparison of these figures shows that the transition from a pure mineral atmosphere to a hydrogenated mineral atmosphere occurs when $P_{\text{H}}^O$ becomes comparable to the mineral vapor pressure (i.e., about 10⁻⁵, 10⁻⁴, 10⁻¹, and 10 bars for T = 1800, 2000, 3000, 3400 K, respectively). These pressures can be converted to a total mass of captured H using the following relation (assuming constant surface acceleration at the planet’s surface: $g=G{M}_{\text{p}}/R_{\text{p}}^2$): $$M_{\text{H}}=\frac{4\pi R_{\text{p}}^4}{G{M}_{\text{p}}}{\displaystyle \sum _i\frac{P_i{\nu }_{\text{H}}^i{\mu }_{\text{H}}}{{\mu }_i},}$$(16)

where M_p, is the planet mass, i is any molecular species in the atmosphere, $v_{\text{H}}^i$, μ_Η, and μ_i are the stochiometric coefficient of H in molecule i, the molar mass of H, and the molar mass of molecule i, respectively. Assuming an Earth-like planet, the conversion between $P_{\text{H}}^0$ and planet mass fraction is displayed in Fig. 6.

We do see that for an Earth-like planet, amounts of H of $P_{\text{H}}^0={10}^{-5}$, 10⁻⁴, 10⁻¹, and 10 bars correspond to planet mass fractions of ~10⁻¹², 10⁻¹¹, 10⁻⁷, and 10⁻⁶, respectively. This demonstrates that a negligible amount of H can efficiently hydrogenate a mineral atmosphere. These mass fractions of H are several orders of magnitude smaller than the initial amount of captured H during the formation of the planet, which ranges from 10⁻³ % to 100% of the planet’s mass (Owen et al. 2020).

Considering the molecular composition of the atmosphere (Figs. 7–10), the same behavior is observed as in the simple two-reaction model (Sect. 2). As the H content increases (on the x axis), the partial pressures of gases derived from evaporation of the magma ocean (those containing Si, Mg, Fe, and K atoms) increase by several orders of magnitude relative to the H-absent case (compare the left- and right-hand sides of Figs. B.1 to B.4), a key result of the present paper. In Appendix B, we provide the atmospheric compositions represented with partial pressures.

In terms of molar fraction, for the T = 1800 K case (Fig. 7), for low H contents ($P_{\text{H}}^0<{10}^{-5}\text{bar}$) the atmosphere is dominated by Na, O₂ and K (in order of decreasing abundance), as in the classic rocky atmosphere case (see e.g., Schaefer et al. 2012), and the total pressure is about 10⁻⁵ bar. However, for cases in which $P_{\text{H}}^0>{10}^{-4}$ bar, the atmospheric composition becomes increasingly dominated by H₂ and H₂O as volatile species, and the dominant metal-bearing species are now Fe and Na (in order of decreasing abundance). Fe is efficiently released and increases linearly with $P_{\text{H}}^0$. For $P_{\text{H}}^0=1\text{\hspace{0.17em}bar}$, NA is 100 times more abundant that in the H-free case (see partial pressure plot in Fig. B.1), and Fr is 1000 times more abundant than in the H-free case. Interestingly, for $P_{\text{H}}^0>{10}^3$ bar, a new gas species becomes dominant (just after H₂): SiH₄ (“Silane”, a toxic gas with a strong repulsive odor) and SiH₃ (“Silanide”; 100–1000 times less abundant).

At 2000 K, the behavior is qualitatively the same as for 1800 K. While Na is the most abundant metallic species in the absence of H for T = 1800 and 2000 K, it is always overtaken by Fe and SiO for H content larger than 10 bar.

For T = 3000 K (Fig. 9), the same qualitative behavior is observed but more species are involved. At this temperature, a pure rocky atmosphere (i.e., low H) is dominated by SiO, Na, O₂, O, Fe, SiO₂, FeO, Mg, and MgO, with a total pressure of about 0.1 bar. When $P_{\text{H}}^0~0.1$ bar, the dominant volatile species become H₂ and H₂O, and the abundance of O₂ decreases. For $1<P_{\text{H}}^0<{10}^4$ bar, the dominant heavy species are SiO, Fe, Mg, and Na. For $P_{\text{H}}^0>{10}^4$ bar, SiO decreases sharply as silicon is converted into SiH₄ and SiH₃. In turn, the $P_{{\text{O}}_2}$ increases again for $P_{\text{H}}^0>{10}^4$ bar because of the liberated oxygen in this process (see Fig. B.3).

Finally, for T = 3400 K (Fig. 10), at low H content, the pure rocky atmosphere has a total pressure of about 3 bar, and is dominated by SiO, O₂, O, SiO_2(g), Na, Fe, and MgO_(g).

For $P_{\text{H}}^0>100$ bar, the atmosphere becomes dominated by H₂, H₂O, and SiO while the abundance of O₂ decreases. For $P_{\text{H}}^0>{10}^4$ bar, SiH₄ and SiH₂ become the most abundant heavy molecules. The shift of this transition to higher initial H contents is mostly due to the temperature dependence of reaction (2), in which the pH₂/pH₂O ratio decreases with increasing temperature, such that the formation of the reduced species silane and silanide are delayed to higher total H contents in hotter atmospheres.

From the above considerations, the following conclusions can be drawn:

• Adding hydrogen to the evaporated atmosphere above the magma ocean increases the partial pressure of all evaporated species owing to the decrease in fO₂. This therefore promotes efficient evaporation of the magma ocean as more material goes into the atmosphere, and will also favor spectroscopic detection as the column density of evaporated species will be higher and the scale heights larger;

• The vapor becomes strongly hydrogenated (where H₂ and H₂O become the dominant species and where the atmosphere chemistry deviates strongly for pure mineral atmosphere) when $P_{\text{H}}^0>P_{\text{RA}}$ where P_RA is the pressure of the rocky atmosphere in the absence of H. Converted into a mass of atmospheric H, this represents from 10⁻¹² to 10⁻⁶ times the planet mass only;

• The appearance of SiH₄ and SiH₃ occurs for $P_{\text{H}}^0>{10}^3$ bar in every case (representing 10⁻³ times the planet mass of hydrogen for a Earth like planet) owing to association reactions; for example, SiO + 2H₂ = SiH₄ + H₂O, in which there are three moles of gas in the reactants compared with only two in the products. At these high H pressures, SiH₄ and SiH₃ become the most abundant heavy molecules (after H₂ and H₂O).

Fig. 5 Total pressure of the gas at equilibrium with the magma ocean for different temperatures.

Fig. 6 Total pressure of monoatomic hydrogen ($P_{\text{H}}^0$) versus planet mass fraction, assuming the planet has the mass and radius of the Earth.

Fig. 7 Molar abundances of most abundant species versus total monoatomic hydrogen pressure $P_{\text{H}}^0$ for an atmosphere at equilibrium with a magma ocean at T = 1800 K. The same plot but with partial pressures is provided in Appendix B.

Fig. 8 Molar abundances of most abundant species versus total monoatomic hydrogen pressure $P_{\text{H}}^0$ for an atmosphere at equilibrium with a magma ocean at T = 2000 K. The same plot but with partial pressures is provided in Appendix B.

Fig. 9 Molar abundances of most abundant species versus total monoatomic hydrogen pressure $P_{\text{H}}^0$ for an atmosphere at equilibrium with a magma ocean at T = 3000 K. The same plot but with partial pressures is provided in Appendix B.

Fig. 10 Molar abundances of most abundant species versus total monoatomic hydrogen pressure ($P_{\text{H}}^0$) for an atmosphere at equilibrium with a magma ocean at T = 3400 K. The same plot but with partial pressures is provided in Appendix B.

4.2 Signature of a hydrogenated magma ocean?

It is interesting to speculate as to how we could detect the presence of a magma ocean beneath a hydrogen-rich atmosphere. In the absence of hydrogen, the results presented above, and in previous studies (Schaefer et al. 2012; Schaefer & Fegley 2010; Lupu et al. 2014; Ito et al. 2015; Ito & Ikoma 2021; Jäggi et al. 2021; Charnoz et al. 2021), show that Na should be released by a magma ocean for T > 1500 K and should remain the dominant species up to 2500 K. At higher temperatures (>2500 K), SiO should be the most abundant mineral species in the vapor.

When H is present - even in small quantities -, H₂O, Fe, Na, and SiO become dominant in the atmosphere. This suggests that the detection of H₂O, Si, Na, or Fe for a rocky planet with an equilibrium temperature of T_eq > 2000 K may point to the presence of a magma ocean below the atmosphere and the presence of a small amount of H. As of today, heavy species, such as Si or Mg, have never been detected in the atmosphere of a rocky or sub-Neptune planet, whereas they have been detected in the atmospheres of giant exoplanets (Zieba et al. 2022). In order to provide a framework for their detection in the future, below we report some atomic ratios that may be of interest for future observations designed to unveil the presence of a magma ocean.

The Mg/Si ratio in the atmosphere above the magma ocean is found to deviate strongly from that of the planet’s mantle (Fig. 11). For $P_{\text{H}}^0<1000$ bar, the atmospheric Mg/si ratio is almost constant at about 10% of the mantle value. For higher hydrogen content, the Mg/Si ratio decreases - by several orders of magnitude – below the values of both the star and the mantle. As Mg and Si are two major elements that condense at relatively high temperatures (~1350 K at 10⁻⁴ bar) from the nebular gas, the fractionation of these elements by condensation or evaporation during planetary accretion should be limited, at least for stars with solar-like C/O ratios (Larimer & Bartholomay 1979). Therefore, assuming that the Mg/Si ratio of the planetary mantle is equal to that of the host star, their detection on an ultrahot rocky planet may provide insights into the magma ocean composition. A caveat associated with this conclusion is that Si remains lithophilic in the planetary mantle, that is, it does not enter the core in appreciable amounts, as is the case on Mercury for example (e.g., Malavergne et al. 2010).

The Mg/Fe ratio of the atmosphere is of interest because our models show that it is independent of H content (Fig. 12). This is because pMg and pFe are proportional to aMgO and aFeO in the liquid phase, and inversely proportional to $f{\text{O}}_2^{1/2}$ (Eq. (14)). Measuring the Mg/Fe ratio in a molten exoplanet atmosphere may provide insights into the Mg/Fe mantle composition independently of the H content of the atmosphere.

On the other hand, the Na/Si ratio of the atmosphere is strongly dependent on the H content of the planet (Fig. 13). If this ratio can be measured (or deduced based on SiO partial pressure, see Wolf et al. 2023), it may provide constraints as to the H content of the atmosphere. Such a conclusion is dependent upon Na still being present in the magma ocean today despite its tendency to evaporate and escape over the lifetime of the planet (see following section, and Erkaev et al. 2022).

Fig. 11 Molar Mg/Si ratio in the vapor above a magma ocean as a function of surface temperature. Colors display the amount of hydrogen in the planet (Log($P_{\text{H}}^0$) in bars). The blue dashed line shows the solar value of Mg/Si and the black dashed line shows the BSE value.

Fig. 12 Molar Mg/Fe ratio in the vapor above a magma ocean as a function of surface temperature. Colors display the amount of hydrogen in the planet (Log($P_{\text{H}}^0$) in bars). The blue dashed line shows the solar value of Mg/Fe and the black dashed line shows the BSE value.

Fig. 13 Molar Na/Si ratio in the vapor above a magma ocean as a function of surface temperature. Colors display the amount of hydrogen in the planet (Log($P_{\text{H}}^0$) in bars). The blue dashed line shows the solar value of Na/Si and the black dashed line shows the BSE value.

5 Atmospheric escape

Now that we are able to compute the atmosphere composition above the magma ocean, we wish to compute the efficiency of escape. This is a very complex problem, and has mostly been investigated for hydrogen-rich planets. The treatment of the full problem using an MHD-hydrocode including multiple species and complex chemistry (as we do here) is well beyond the scope of the present paper. Still, to obtain an approximate estimate of the escape efficiency of such an atmosphere, we developed a model inspired by the popular energy-limited model, which we present below.

5.1 Method of calculation

We try to quantify the efficiency of atmospheric escape from planets with a long-lived magma ocean embedded in a primordial hydrogen envelope. We focus on planets very close to their star so that stellar heating is strong enough to maintain a magma ocean during the lifetime of the star. While it is well established that hot-Jupiter planets are massive enough to retain their atmospheres (see e.g., Valencia et al. 2010), Earth-sized planets should have lost their primordial H envelope in less than 10 Myr (e.g., Valencia et al. 2010; Lopez 2017; Erkaev et al. 2022). What happens in the intermediate regime of super-Earths/sub-Neptunes is comparatively poorly understood. Of particular importance is the observation of a bimodal radius distribution of super-Earths, with a well-identified valley in the radius range 1.5–2 R_⊕ (see e.g., Fulton et al. 2017; Gupta & Schlichting 2019). It has been proposed that this “valley” marks the transition from the smaller rocky planets (“super-Earths”) to planets containing a few per cent of the planet’s total mass in a H/He-rich envelope (“sub-Neptunes”; see e.g., Owen & Wu 2013; Lopez & Fortney 2014; Ginzburg et al. 2016; Ginzburg & Sari 2017). This transition may be due to rapid photoevaporation of highly irradiated H-rich atmosphere, which may not survive for planets with radii of < 1.5–2R_⊕. In this section, we do not investigate the origin of this transition; rather, we study the consequence of the presence of a small amount of hydrogen on the efficiency of magma evaporation and escape of the atmosphere produced in such a scenario.

Here, we consider the case of an atmosphere in thermodynamic equilibrium with a magma ocean of constant temperature at its base as described in the previous section. Models of strongly irradiated planets show that escape is always efficient for terrestrial-mass planets (Ito et al. 2015; Benedikt et al. 2020; Ito & Ikoma 2021) and that the particle velocities in the atmospheres of these planets far exceed those necessary for hydrostatic equilibrium. Although this conclusion is unlikely to hold for a cold atmosphere (which may indeed be hydrostatic), this assertion may be justified for ultrahot planets (T > 1800 K) with irradiated atmospheres, depending on the average molar mass of the atmosphere. For heavy species (Si, Mg), Ito & Ikoma (2021) find strong atmospheric escape (with mass flux of about 0.4 M_⊕/Gyr⁻¹) in coupled hydrodynamical and radiative simulations, for an Earth-size planet at 0.02 AU devoid of hydrogen. These authors conclude that the atmosphere is never at hydrostatic equilibrium and is in “blow-off” regime, with a Parker-wind-type transonic velocity profile. Konatham et al. (2020) argue that significant to catastrophic escape occurs if the atomic thermal velocities reach >10% of a planet’s escape velocities. This means that an Earth-like planet at a temperature of 3000 K may lose all molecules up to 100 atomic mass units, whereas a planet of 2 Earth masses may lose atoms up to gasses with molar weight up to ~25 atomic mass units. Similarly, Johnstone et al. (2019) find that even Ar (40 AMU) can escape from a strongly irradiated Earth-like planet (but these authors do not include heavier species). In the latter study, the H/N/O/C/He/Ar atomic ratios in the upper escaping atmosphere are found to vary between 0.75 and 2.5 relative to the same ratios in the lower atmosphere. As long as the atmosphere is in the hydrodynamic escape regime (blow-off), it seems that no strong atomic fractionation occurs during the escape, within a factor of a few (Johnstone et al. 2019).

In all the studies mentioned above, the gravitational stratification of species is not found because diffusive timescales are always longer than advective timescales. This is why gravitational stratification is ignored here.

We note that cooling via vibrational bands of H₂O (Yoshida et al. 2022) and atomic lines for N₂ and O₂ (Nakayama et al. 2022) has been found to be efficient enough to cool down the atmosphere and prevent escape from habitable planets under very strong XUV irradiation (up to about 1000 times the present Earth value, which is relevant for very young stars). However, in the case we consider here, that is, ultrahot and ultrashort-period exoplanets, the XUV flux may be > 10⁴ times the value for Earth today owing to the very short stellar distance (0.01–0.02 au). In such extreme cases, Yoshida et al. (2022) do find hydrodynamical escape.

A full treatment of the computation of atmospheric escape is beyond the scope of the present paper. Here, we limit our analysis to simple estimates based on the popular “energy limited” approximation, following the method of Salz et al. (2016), which provides a useful scaling law for energy-limited escape; this law is calibrated in 1D numerical simulations for irradiated, hydrogen-dominated planets. The mass flux (Kg/s) escaping from the atmosphere is approximated by: $$\dot{M}=\frac{-3{\beta }^2\eta F_{\text{XUV}}}{4KG{\rho }_{\text{p}}},$$(17)

where β = R_XUV/R_p is the ratio of the radius of absorption of the radiation from an XUV star to the physical radius of the planet (R_p), and η is a parameter standing for the “heating efficiency”, which is the fraction of input radiative energy converted into thermal expansion of the atmosphere. F_XUV is the incident stellar flux in the XUV band (in J⁻¹ s⁻¹ m⁻²) and K is a parameter accounting for the effect of stellar tides (which facilitate escape). Here, we use K = 1 for simplicity.

Numerous physical mechanisms are subsumed into η and its value is highly dependent on the stellar spectrum, the gravity of the planet, and the atmospheric composition (see e.g., Valencia et al. 2010; Salz et al. 2016; Ito & Ikoma 2021). Salz et al. (2016) find that β < 2 and depends sensitively on the planet’s gravitational potential well and the flux of the star. The following parametrization is proposed (Salz et al. 2016): $$\begin{array}{l}Lo{g}_{10}\left(\beta \right)=\text{Max}[0.,-0.185{\text{Log}}_{10}\left(\frac{{\text{Φ}}_{\text{G}}}{erg/g}\right)\hfill \\ +0.021{\text{Log}}_{10}\left(\frac{F_{\text{XUV}}}{\text{erg\hspace{0.17em}}{\text{cm}}^{-2}{\text{s}}^{-1}}\right)+2.42],\hfill \end{array}$$(18)

where Φ_G = GM_p/R_p. Numerical simulations by Salz et al. (2016), which are valid for a H₂-dominated atmosphere, show that for planets with a deep gravitational well (i.e., large Φ), escape is difficult and η may be ≪1 (giant planets), whereas for planets with shallow gravitational potential (i.e., small Φ), escape is relatively efficient (Earth- and super-Earth-mass planets). Following Salz et al. (2016) and noting $v={\text{Log}}_{10}\left(\frac{{\Phi }_{\text{G}}}{{\text{erg g}}^{-1}}\right)$ η_Η for a hydrogen-dominated atmosphere: $$Lo{g}_{10}\left({\eta }_{\text{H}}\right)=\{\begin{array}{ll}-0.5-0.44\left(\nu -12.\right)\hfill & \text{for\hspace{0.17em}}12\text{.<}\nu <13.11\hfill \\ -0.98-7.29\left(\nu -13.11\right)\hfill & \text{for\hspace{0.17em}}13\text{.11<}\nu <\mathrm{13.6.}\hfill \end{array}$$(19)

The above fit, which is valid for H₂-dominated atmospheres, may not be valid for pure mineral atmospheres. Ito & Ikoma (2021) find that, in pure mineral atmospheres, atmospheric cooling is much more efficient than in H₂-dominated atmospheres, resulting in low escape efficiency with 2 × 10⁻⁴ < η < 4 × 10⁻³, which is much lower than 0.1 < η < 0.5 for H-dominated atmospheres. To take this compositional effect into account, we replace η in Eq. (17) by: $${\eta }_{\text{eff}}={\eta }_{\text{H}}\times \left(\frac{1.-{\text{e}}^{-x}}{1.-{\text{e}}^{-1}}\right)+5\times {10}^{-4},$$(20)

where x is the atomic fraction of H in the atmosphere. Although not optimal, the above prescription ensures that H₂-dominated atmospheres escape at the efficiency given by Salz et al. (2016), whereas mineral atmospheres escape with an efficiency in the (lower) range of the values reported by Ito & Ikoma (2021).

Finally, concerning the time-dependent XUV flux at the location of the planet, we adopt the fit to XUV flux data reported for young stars by Valencia et al. (2010): $$F_{\text{XUV}}\left(\text{J\hspace{0.17em}}{\text{m}}^{-2}{\text{s}}^{-1}\right)=\{\begin{array}{l}29.7\times {10}^{-3}{\left(0.1\right)}^{-1.23}\left(\frac{1}{a^2}\right)\text{\hspace{0.17em}\hspace{0.17em}for\hspace{0.17em}t\hspace{0.17em}<0}.1\text{\hspace{0.17em}Gyr}\hfill \\ 29.7\times {10}^{-3}{\left(\frac{t}{1\text{\hspace{0.17em}Gyr}}\right)}^{-1.23}\left(\frac{1}{a^2}\right)\text{\hspace{0.17em}for\hspace{0.17em}t\hspace{0.17em}>0}.1\text{\hspace{0.17em}Gyr.}\hfill \end{array}$$(21)

The calculation is performed as follows. We start with a planet mass M_p and radius R_p with surface acceleration g, surface temperature T, and mass of H in the atmosphere M_H, corresponding to an initial pressure of monoatomic $\text{H}:P_{\text{H}}^0=M_{\text{H}}*g/\left(4\pi R_{\text{p}}^2\right)$. For a given T and $P_{\text{H}}^0$, the atmospheric molecular composition above the magma ocean is computed as described in Sect. 3. We compute the total mass loss according to Eq. (17). We assume the star is 5 Gyr old. Following the two preceding arguments regarding inefficient gravitational stratification in the blow-off regime, the total mass flux of each atom is Ṁ × X_i, where X_i is the atomic mass fraction of atom i in the escaping gas.

We do not consider here the dissolution of H₂ or H₂O in the magma ocean as this will be the subject of a future study. We note that the dissolution of H₂ and H₂O into the magma ocean may lead to the constitution of a reservoir of H and O, which may lengthen the duration of H escape (Hier-Majumder & Hirschmann 2017; Bower et al. 2022). Consequently, the timescales of H escape reported in the following section should be considered as lower bounds. In addition, the substantial solubility of H as H₂O in the magma ocean (Sossi et al. 2023) indicates that it would be largely (≥95% of its total budget) dissolved should the entire mantle of an Earth-mass planet remain molten. However, as the magma ocean cools, the solubility of H₂O exceeds the mass fraction of available liquid in which to dissolve, leading to late-stage release of significant amounts of O and H (as H₂O) into the atmosphere when the atmospheric pressure is sufficiently low (Solomatova & Caracas 2021). This may serve to increase the fO₂ relative to that initially imposed by the H₂-rich nebular gas, and thus modify atmospheric chemistry. We warn the reader that these effects are not considered here.

Fig. 14 Escape from a 1 M⊕ and R⊕ planet located at 0.02 au from its star. Top: mass of H vs. time for different initial mass fractions of hydrogen. Bottom: mass fraction of lost Si(red) and Na(blue) vs. time. The solid line shows these relations for an initial H mass fraction of 10−6 times the planet mass, and the dashed line is for an initial H mass fraction of 0.1 times the planet mass.

5.2 Results

Figure 14 displays the mass of hydrogen as a function of time for an Earth-like planet (M_p = 1 M_⊕, R_p = 1 R_⊕) located at 0.02 au, starting with different hydrogen content from 10⁻⁶ to 10⁻¹ times the planet mass. The top panel shows that most H is lost within about 100 Kyr to 10 Myr, depending on $M_{\text{H}}^0/M_{\text{p}}$. The bottom panel displays the lost mass of Na (a proxy for the moderately volatile elements) divided by its mass in the planet’s mantle magma ocean (in red; we assume to first order that the magma ocean is 50% of the planet mass). For $M_{\text{H}}^0/M_{\text{p}}={10}^{-6}$ (red solid line), Na is lost at a roughly constant rate as long as hydrogen is present and the star is <10 Myr old when its XUV flux is intense. When hydrogen has disappeared, the efficiency of Na escape increases. This result may be surprising, but can be reconciled with the fact that, for high H contents, the Na mole fraction is very low; that is, the atmosphere is dominated by SiO or Fe (see previous section). In contrast, when the abundance of H drops to zero, the mole fraction of the atmosphere comprised of H increases by three orders of magnitude, meaning the escape rate of Na increases by a proportionate factor. Over gigayear timescales, the loss of Na is dominated by loss during the last 1 Gyr and we find, at the end, that about 30% of the mantle Na content is lost.

The role of H in influencing the Si content of the planet’s mantle is more marked. The blue solid line in Fig. 14 shows Si evolution for $M_{\text{H}}^0/M_{\text{p}}={10}^{-6}$. After 5 Gyr, the planet devoid of H has lost only ~1% of its Si. Conversely, for planets starting with hydrogen contents of 10% of the mass of the entire planet, about 10% of the total budget of Si is lost after 5 Gyr. This strong increase in Si loss with H content reflects the fact that when H is abundant, SiO becomes the dominant metal-bearing gas species. In addition, the presence of H increases the efficiency of atmospheric heating (η > 0.1, see Sect. 5.1) and decreases the mean molar mass of the atmosphere. We highlight the fact that, due to the simplicity of the models, the above values may not be understood as exact solutions, but rather as order-of-magnitude estimates of the extent of escape. Indeed, the values could be somewhat underestimated owing to the conservative approach that was taken in estimating the quantities used in Sect. 5.1, though it should be mentioned that the composition of the magma ocean is assumed to remain fixed rather than evolve, as would be expected for a finite mass undergoing Rayleigh (fractional) vaporisation.

We also explored different planet masses and radii and find that the main factors limiting the amount of metals and metal-bearing oxides and hydrides able to escape is the lifetime of the H envelope (controlled by the planet’s density) and the surface temperature. For a more complete view of the effect of escape on each atmospheric constituent, Fig. 15 shows the lost mass fraction of every atom (O, Si, Mg, Fe, Na, K) for all combinations of temperature and initial mass fraction of hydrogen and for an Earth-like planet located at 0.02 au. The color of the plot indicates the mass fraction of each atom that is lost (with respect to the content in the mantle). White denotes that all atoms in the mantle are lost, purple indicates a loss of about 10%, and so on. The most volatile elements, Na and K, are completely lost for T_surf < 3500 K, which is in line with the results of Ito & Ikoma (2021). For T_surf > 3500 K, the efficiency of Na and K loss is decreased because their molar fraction in the atmosphere is supplanted by other species (Si, Fe, Mg), thus diminishing their net loss in comparison to other species. Indeed the mass flux (Eq. (17)) only depends on the F(XUV) flux and the planet mass and radius. In addition, the heating efficiency is low (<0.001), such that when temperature is >3200 K and H content is low, the species that escape are mostly those that predominate in the atmosphere, that is, mostly Si, O, and Mg. If we now turn to Si, we see that more than 10% of the planet’s mantle Si is lost for T_surf > 2600 K and for a H mass fraction of M_H/M_p > 0.1. The same result applies for oxygen. As oxygen is the most abundant species by mass (about 50% of the mass of Earth’s mantle), this implies significant removal of the planet’s mantle and an increase in average density (an effect that is not taken into account here). Concerning Fe, it is substantially lost (>10%) for temperatures >2600 K and $M_{\text{H}}^0/M_{\text{p}}>0.01$ due to its high abundance in the vapor combined with its low abundance (6.3 mol%) in the mantle.

Magnesium shows an intermediate behavior between those of Si and Na, with most efficient loss occurring for T_surf > 3200 K and $M_{\text{H}}^0/M_{\text{p}}>{10}^{-2}$. Up to 8% of Mg can be lost from the mantle in such extreme conditions. The less efficient loss of Mg stems from its lower volatility than Si and Fe for the BSE composition considered here.

Therefore, in conclusion, we do see that for extreme surface temperatures (>2600 K) and initial H mass content of >1%, a significant fraction of the mantle’s content (~ a few 10%) could be evaporated during the planet’s evolution. We note that the results of this section must be taken with care because of the simplicity of the model and the numerous simplifications we had to introduce. The energy-limited escape is known to overestimate the escape, and the escape efficiency is not well known for atmosphere heavier than solar composition. In addition, vibrational and atomic radiative cooling lines during the escape (see e.g., Yoshida et al. 2022; Nakayama et al. 2022) may act against escape. We therefore conclude that our simple model suggests that escape of heavy molecules may be made easier by the presence of H, but further investigation of this process is required to confirm this finding.

Fig. 15 Each figure shows a color plot of the mass fraction of an atom that is lost for temperature T (x axis) and initial hydrogen mass (y axis) of the planet. The planet’s solid core here is 1 M⊕ and 1 R⊕.

6 Conclusions

We have studied the composition and the escape of an atmosphere above a magma ocean planet in the presence of hydrogen. In the present study, a special effort was dedicated to computing the atmospheric composition at thermodynamical equilibrium with the magma ocean, taking into account the effect of H on the atmospheric chemistry, and its retro-action on the liquid-atmosphere equilibrium. Our main results are summarized below:

• Hydrogen present in the atmosphere has a dramatic effect on the vapor composition. For T = 2000 K and 3000 K, the atmosphere composition is strongly modified for atmospheric hydrogen contents of about 10⁻⁸ and 10⁻⁶ times the planet mass, respectively (assuming an Earth mass planet) and becomes dominated by H₂ and H₂O;

• Hydrogen lowers the O₂ partial pressure (= oxygen fugac-ity fO₂) to make H₂O, which happens in proportion to the mass of H added. As a consequence, the addition of hydrogen promotes efficient degassing of metals Na, K, Fe, Mg, and Si. These elements all evaporate according to reactions, M^x+sO_(x+s)/2(l = M^xO_x/2(g) + s/4O₂(g) in which the stoichiometric coefficient, s, is positive. Therefore, a small amount of hydrogen (comparable to the total pressure of the evaporated elements, and nearly negligible compared to the planet’s mass) greatly enhances the extraction of these elements from the magma to the vapor phase;

• The atomic abundances of heavy (i.e., rock-forming) species in the atmosphere may diverge significantly from the mantle value, and are controlled by the magma-atmosphere equilibrium. The Mg/Si ratio in the atmosphere is constant and independent of H content for T > 3000 K and may be directly related to the composition of the planet’s mantle. Conversely, the Na/Si ratio is very sensitive to H content (owing to the different dependence of the partial pressures of their major gas species on fO₂, and may be useful in constraining the presence of H). Finally, the Mg/Fe ratio is independent of H and may be linked to the mantle composition. This suggests that spectroscopic observations of magma ocean planets could provide a means to constrain the mantle composition, even in the presence of hydrogen;

• For Earth-sized planets of <0.02 AU and with a temperature above the magma ocean, that is >2000 K, most H is lost in <10 Myr, regardless of its speciation;

• For Earth-sized planets of <0.02 AU and with a temperature above the magma ocean, that is >3000 K, we find significant loss of Fe, Si, and O (>10% of their budget in the planet’s mantle). In such extreme conditions, this may lead to a decrease in the mass of the mantle, and in turn to an increase in the mean density of the planet (due to core/mantle ratio decreases). This may potentially explain the high densities of strongly irradiated rocky planets.

In conclusion, we show that the presence of relatively small fractions of H in the atmosphere may enhance the evaporation of metals from ultrahot rocky planets covered by a magma ocean, and, at very extreme temperatures, may lead to significant loss of heavy atoms such as Si, O, or Fe. Whereas no heavy species have yet been detected in rocky planetary atmospheres, we identify ratios of rock-forming elements in the vapor phase that are potentially able to reveal the presence and nature of a magma ocean, and subsequently probe the composition of the planet mantle. In contrast to some recent studies (Zilinskas et al. 2022), we do find that the presence of H strongly increases the partial pressure of major degassed species (see Figs. B.1 to B.4). This will significantly impact the resulting atmospheric spectra and (P, T) structure of the exoplanet. This will be the subject of a future study.

In the future, we will study the effect of dissolution of H₂ and H₂O on the atmospheric composition of magma ocean exo-planets. We anticipate that the effect of dissolution of H₂O will be to oxidize the atmosphere, thus lowering the reducing effect of H, which may lead to less efficient release of metallic species (boosted by H₂); it will also create an additional reservoir of H, probably resulting in an extended escape time of H. This will be studied in a future paper.

On a more technical note, we also emphasize that improvement of ab initio calculations and comparison with tabulated thermochemical data will be useful in the future, as ab initio calculations reveal the formation of a great diversity of species (sometimes not available in common data tables such as JANAF or CEA-NASA) and at higher temperatures than are typically accessible experimentally.

Acknowledgements

S.C., A.F., P.T. acknowledge financial support by LabEx UnivEarthS (ANR-10-LABX-0023 and ANR-18-IDEX-0001) and by the CNES (Centre National d’Études Spatiales). P.T. would also like to acknowledge and thank the ERC for funding this work under the Horizon 2020 program project ATMO (ID: 757858). R.C. acknowledges support from the European Research Council under EU Horizon 2020 research and innovation program (grant agreement 681818 – IMPACT to RC), the Research Council of Norway with project number 223272 and project HIDDEN 325567, and access to supercomputing facilities via eDARI stl2816 grants, PRACE RA4947 grant, Uninet2 NN9697K grant. P.S. thanks the Swiss National Science Foundation (SNSF) via an Ambizione Fellowship (180025), an Eccellenza Professorship (203668) and the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number MB22.00033, a SERI-funded ERC Starting Grant ‘2ATMO’. We thank the anonymous reviewer for their insightful comments that improved the quality of the paper.

Appendix A Composition of vapor in the absence of H

Figure A.1 shows the vapor composition at equilibrium with a liquid with BSE composition, computed with the procedure described in (Charnoz et al. 2021); it compares very well with composition reported in Ito et al. (2015) Visscher & Fegley (2013) van Buchem et al. (2022).

Fig. A.1 Molar fractions in the vapor at equilibrium with a magma with BSE composition, and without hydrogen.

Appendix B Partial pressures in the hydrogenated atmosphere

Fig. B.1 Molar abundances of most abundant species versus total monoatomic hydrogen pressure $P_{\text{H}}^0$ for an atmosphere at equilibrium with a magma ocean at T=1800 K.

Fig. B.2 Molar abundances of most abundant species versus total monoatomic hydrogen pressure $P_{\text{H}}^0$ for an atmosphere at equilibrium with a magma ocean at T=2000 K.

Fig. B.3 Molar abundances of most abundant species versus total monoatomic hydrogen pressure $P_{\text{H}}^0$ for an atmosphere at equilibrium with a magma ocean at T=3000 K.

Fig. B.4 Molar abundances of most abundant species versus total monoatomic hydrogen pressure $P_{\text{H}}^0$ for an atmosphere at equilibrium with a magma ocean at T=3400 K.

References

All Tables

All Figures

Fig. 1 Density of known exoplanets versus their equilibrium temperatures. These data come from the online DACE database maintained by Observatory of Geneva (https://dace.unige.ch/) and are taken from the catalog of Otegi et al. (2020), which is an exoplanet catalog based on reliable, robust, and - to the greatest possible extent – accurate mass and radius measurements of transiting planets up to 120 earth masses (Otegi et al. 2020). Top: density vs. equilibrium temperature for all planets. Bottom: same as above but limited to planets with radii of smaller than 3 Earth radii. The dotted line shows a simple least-square fit to the data.

In the text

Fig. 2 Schematic illustration of the internal structure of a rocky planet as assumed in this paper, with an atmosphere containing a mixture of captured hydrogen and gaseous species derived from a molten mantle that forms a global magma ocean below the planet atmosphere. In the centre of the planet, an iron/nickel core is present. The composition of the magma ocean is assumed to be fixed, while the composition of the atmosphere is at thermodynamic equilibrium with that of the magma ocean.

In the text

	Fig. 3 Pressures of the different elements for the simple two-reaction model with only SiO(g), H2, H2O, and O2. Here, the temperature is fixed to 3000 K. SiO solution 1 and SiO solution 2 refer to asymptotic solutions with low and high hydrogen content, respectively.
In the text

	Fig. 4 Pressures of the different elements for the simple two-reaction model with only SiO(g), H2, H2O, and O2. Here, the hydrogen content is fixed ($P_{\text{H}}^0=1\text{bar}$). SiO solution 1 and SiO solution 2 refer to asymptotic solutions with low and high hydrogen content, respectively.
In the text

	Fig. 5 Total pressure of the gas at equilibrium with the magma ocean for different temperatures.
In the text

	Fig. 6 Total pressure of monoatomic hydrogen ($P_{\text{H}}^0$) versus planet mass fraction, assuming the planet has the mass and radius of the Earth.
In the text

	Fig. 7 Molar abundances of most abundant species versus total monoatomic hydrogen pressure $P_{\text{H}}^0$ for an atmosphere at equilibrium with a magma ocean at T = 1800 K. The same plot but with partial pressures is provided in Appendix B.
In the text

	Fig. 8 Molar abundances of most abundant species versus total monoatomic hydrogen pressure $P_{\text{H}}^0$ for an atmosphere at equilibrium with a magma ocean at T = 2000 K. The same plot but with partial pressures is provided in Appendix B.
In the text

	Fig. 9 Molar abundances of most abundant species versus total monoatomic hydrogen pressure $P_{\text{H}}^0$ for an atmosphere at equilibrium with a magma ocean at T = 3000 K. The same plot but with partial pressures is provided in Appendix B.
In the text

	Fig. 10 Molar abundances of most abundant species versus total monoatomic hydrogen pressure ($P_{\text{H}}^0$) for an atmosphere at equilibrium with a magma ocean at T = 3400 K. The same plot but with partial pressures is provided in Appendix B.
In the text

	Fig. 11 Molar Mg/Si ratio in the vapor above a magma ocean as a function of surface temperature. Colors display the amount of hydrogen in the planet (Log($P_{\text{H}}^0$) in bars). The blue dashed line shows the solar value of Mg/Si and the black dashed line shows the BSE value.
In the text

	Fig. 12 Molar Mg/Fe ratio in the vapor above a magma ocean as a function of surface temperature. Colors display the amount of hydrogen in the planet (Log($P_{\text{H}}^0$) in bars). The blue dashed line shows the solar value of Mg/Fe and the black dashed line shows the BSE value.
In the text

	Fig. 13 Molar Na/Si ratio in the vapor above a magma ocean as a function of surface temperature. Colors display the amount of hydrogen in the planet (Log($P_{\text{H}}^0$) in bars). The blue dashed line shows the solar value of Na/Si and the black dashed line shows the BSE value.
In the text

	Fig. 14 Escape from a 1 M⊕ and R⊕ planet located at 0.02 au from its star. Top: mass of H vs. time for different initial mass fractions of hydrogen. Bottom: mass fraction of lost Si(red) and Na(blue) vs. time. The solid line shows these relations for an initial H mass fraction of 10−6 times the planet mass, and the dashed line is for an initial H mass fraction of 0.1 times the planet mass.
In the text

	Fig. 15 Each figure shows a color plot of the mass fraction of an atom that is lost for temperature T (x axis) and initial hydrogen mass (y axis) of the planet. The planet’s solid core here is 1 M⊕ and 1 R⊕.
In the text

	Fig. A.1 Molar fractions in the vapor at equilibrium with a magma with BSE composition, and without hydrogen.
In the text

	Fig. B.1 Molar abundances of most abundant species versus total monoatomic hydrogen pressure $P_{\text{H}}^0$ for an atmosphere at equilibrium with a magma ocean at T=1800 K.
In the text

	Fig. B.2 Molar abundances of most abundant species versus total monoatomic hydrogen pressure $P_{\text{H}}^0$ for an atmosphere at equilibrium with a magma ocean at T=2000 K.
In the text

	Fig. B.3 Molar abundances of most abundant species versus total monoatomic hydrogen pressure $P_{\text{H}}^0$ for an atmosphere at equilibrium with a magma ocean at T=3000 K.
In the text

	Fig. B.4 Molar abundances of most abundant species versus total monoatomic hydrogen pressure $P_{\text{H}}^0$ for an atmosphere at equilibrium with a magma ocean at T=3400 K.
In the text