The effect of a small amount of hydrogen in the atmosphere of ultrahot magma-ocean planets: atmospheric composition and escape

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_{\mathrm{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:

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

For a given value of $P_{\mathrm{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^{g}_{i}$ (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^{g}_{i},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 f$\mathrm{O_{2}}$. Therefore, the following relations must be obeyed for the gas for each reaction $j$ reported in Table 1.

  $$P_{vapor_{j}}=\frac{K_{j}(T)a(liquid_{j})}{P_{\mathrm{O_{2}}}^{s_{j}}},$$

  (14)

  where $vapor_{j}$ stands for the vapor species that bears the metal in evaporation reaction $j$ (e.g., $Na_{(g)}$ in reaction $\#3$), $liquid_{j}$ is the corresponding liquid oxide ($NaO_{0.5}$ for reaction $\#3$), $s_{j}$ is the stochiometric 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}\times\Gamma_{j}$, where $X_{j}$ is the mole fraction of $liquid_{j}$ in the magma ocean (assumed to be fixed) and $\Gamma_{j}$ is the activity coefficient of $j$ in the liquid due to nonideal mixing effects. Activity coefficients, $\Gamma_{j}$, of all melt oxide components are interpolated using outputs of the VapoRock code (Wolf et al., 2022).

  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^{g}_{\mathrm{O}}=2X^{g}_{\mathrm{Si}}+X^{g}_{\mathrm{Mg}}+\frac{1}{2}X^{g}_{\mathrm{Na}}+\frac{1}{2}X^{g}_{\mathrm{K}}+X^{g}_{\mathrm{Fe}},$$

  (15)

  where $X^{g}_{i}$ 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_{\mathrm{H}}^{0}$, meaning that Equation 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^{l}_{i}$), temperature, and a total hydrogen content ($P_{\mathrm{H}}^{0}$) .

• •

  Determine the atomic composition of the atmosphere ($X^{g}_{i}$ for all atoms i) and the total pressure $P_{tot}$. For all sets ($X^{g}_{i}$, $P_{tot}$), the corresponding molecular composition of the atmosphere is computed using the CEA/NASA chemical equilibrium code. The only solution to this problem is the ($X^{g}_{i}$, $P_{tot}$) set for which the molecular composition of the atmosphere simultaneously solves Equations 13 to 15.

• •

  The search is done using an iterative method, with a simplex minimization algorithm. The calculation goes on until Equations 13 to 15 are solved up to a relative accuracy of better than $10^{-3}$ 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, $\mathrm{H_{2}}$, $\mathrm{H_{2}O}$, $\mathrm{H_{2}O_{2}}$, $\mathrm{HO_{2}}$, K, $\mathrm{K_{2}}$, $\mathrm{K_{2}O}$, $\mathrm{K_{2}O_{2}}$, $\mathrm{K_{2}O_{2}H_{2}}$, $\mathrm{KH}$, $\mathrm{KNa}$, $\mathrm{KO}$, $\mathrm{KOH}$, $\mathrm{Mg}$, $\mathrm{Mg_{2}}$, $\mathrm{MgH}$, $\mathrm{MgO}$, $\mathrm{MgOH}$, $\mathrm{Na}$, $\mathrm{Na_{2}},\mathrm{Na_{2}O}$, $\mathrm{Na_{2}O_{2}}$, $\mathrm{Na_{2}O_{2}H_{2}}$, $\mathrm{NaH}$, $\mathrm{NaO}$, $\mathrm{NaOH}$, $\mathrm{O}$, $\mathrm{O_{2}}$, $\mathrm{O_{3}}$, $\mathrm{OH}$, $\mathrm{Si}$, $\mathrm{Si_{2}}$, $\mathrm{Si_{3}}$, $\mathrm{SiH}$, $\mathrm{SiH_{2}}$, $\mathrm{SiH_{3}}$, $\mathrm{SiH_{4}}$, $\mathrm{SiO}$, $\mathrm{SiO_{2}}$.

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., 2022) 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 Figure 4 of Wolf et al. (2022) and van Buchem et al. (2022)).

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–vapor dome at 3000 K and a density of 0.4 g/cm3, 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-Ernzernhof formulation (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.

Our calculations are performed between T = 1800 K and 3500 K, and for total hydrogen pressures ranging from $10^{-6}$ to $10^{6}$ 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 $\sim$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 Figure 5. Figure 6 provides the conversion between H pressure and mass of hydrogen (with respect to the planet’s mass). The composition of the atmospheres at equilibrium with the magma oceans are plotted in Figs. 7 to 10. A rapid comparison of these figures shows that the transition from a pure mineral atmosphere to a hydrogenated mineral atmosphere occurs when $P_{H}^{O}$ becomes comparable to the mineral vapor pressure (i.e., about $10^{-5}$, $10^{-4}$, $10^{-1}$, 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=GM_{p}/R_{p}^{2}$):

where $M_{p}$, is the planet mass, $i$ is any molecular species in the atmosphere, $\nu^{i}_{H}$, $\mu_{H}$, and $\mu_{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_{H}^{0}$ and planet mass fraction is displayed in Figure 6.

We do see that for an Earth-like planet, amounts of H of $P_{H}^{0}=10^{-5},10^{-4},10^{-1}$, and $10$ bars correspond to planet mass fractions of $\sim 10^{-12},10^{-11},10^{-7}$, and $10^{-6}$, 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^{-3}$% to 100$\%$ of the planet’s mass (Owen et al., 2020).

Considering the molecular composition of the atmosphere (Figs. 7, 8, 9, 10), the same behavior is observed as in the simple two-reaction model (section 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. 17 to 20), 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_{\mathrm{H}}^{0}<10^{-5}$bar) the atmosphere is dominated by Na,$\mathrm{O_{2}}$ 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^{-5}$ bar. However, for cases in which $P_{\mathrm{H}}^{0}>10^{-4}$ bar, the atmospheric composition becomes increasingly dominated by $\mathrm{H_{2}}$ and $\mathrm{H_{2}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_{H}^{0}$. For $P_{H}^{0}$=1 bar, NA is 100 times more abundant that in the H free case(see partial pressure plot in Appendix 17), and Fr is 1000 times more abundant than in the H-free case . Interestingly, for $P_{\mathrm{H}}^{0}>10^{3}$ bar, a new gas species becomes dominant (just after $\mathrm{H_{2}}$): $\mathrm{SiH_{4}}$ (”Silane”, a toxic gas with a strong repulsive odor) and $\mathrm{SiH_{3}}$ (”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_{\mathrm{H}}^{0}\sim 0.1$ bar, the dominant volatile species become $\mathrm{H_{2}}$ and $\mathrm{H_{2}O}$, and the abundance of $\mathrm{O_{2}}$ decreases. For $1<P_{\mathrm{H}}^{0}<10^{4}$ bar, the dominant heavy species are SiO, Fe, Mg, and Na. For $P_{\mathrm{H}}^{0}>10^{4}$ bar, SiO decreases sharply as silicon is converted into $\mathrm{SiH_{4}}$ and $\mathrm{SiH_{3}}$. In turn, the $P_{\mathrm{O_{2}}}$ increases again for $P_{\mathrm{H}}^{0}>10^{4}$ bar because of the liberated oxygen in this process (see Fig. 19).

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, $\mathrm{SiO}_{2(g)}$ , Na, Fe, and $\mathrm{Mg}\mathrm{O}_{(g)}$.

For $P_{\mathrm{H}}^{0}>100$ bar, the atmosphere becomes dominated by $\mathrm{H_{2}}$, $\mathrm{H_{2}O}$, and SiO while the abundance of $\mathrm{O_{2}}$ decreases. For $P_{\mathrm{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 $f\mathrm{O_{2}}$. 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 $\mathrm{H_{2}}$ and $\mathrm{H_{2}O}$ become the dominant species and where the atmosphere chemistry deviates strongly for pure mineral atmosphere) when $P_{\mathrm{H}}^{0}>P_{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^{-12}$ to $10^{-6}$ times the planet mass only.

• •

  The appearance of $\mathrm{SiH_{4}}$ and $\mathrm{SiH_{3}}$ occurs for $P_{\mathrm{H}}^{0}>10^{3}$ bar in every case (representing $10^{-3}$ 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, $\mathrm{SiH_{4}}$ and $\mathrm{SiH_{3}}$ become the most abundant heavy molecules (after $\mathrm{H_{2}}$ and $\mathrm{H_{2}O}$).

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; Schaefer et al., 2012; 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—, $\mathrm{H_{2}O}$, Fe, Na, and SiO become dominant in the atmosphere. This suggests that the detection of $\mathrm{H_{2}O}$, Si, Na, or Fe for a rocky planet with an equilibrium temperature of $T_{eq}>2000K$ 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_{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 ($\sim$ 1350 K at 10^-4 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 fO${}_{2}^{1/2}$ (cf. 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. (2022)), it may provided 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)).

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_{\oplus}$ (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_{\oplus}$. 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_{\oplus}/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 $\sim 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 $\mathrm{H_{2}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^{4}$ 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. (2015), 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:

where $\beta=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 $\eta$ 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^{2}$) 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 $\eta$ 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 $\beta$¡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):

where $\Phi_{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 $\Phi$), escape is difficult and $\eta$ may be $\ll\leavevmode\nobreak\ 1$ (giant planets), whereas for planets with shallow gravitational potential (i.e., small $\Phi$), escape is relatively efficient (Earth- and super-Earth-mass planets). Following Salz et al. (2016) and noting $\nu=Log_{10}\left(\frac{\Phi_{G}}{erg/g}\right)$ $\eta_{H}$ for a hydrogen-dominated atmosphere:

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\times 10^{-4}<\eta<4\times 10^{-3}$, which is much lower than $0.1<\eta<0.5$ for H-dominated atmospheres. To take this compositional effect into account, we replace $\eta$ in Eq. 17 by:

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):

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 H: $P_{\mathrm{H}}^{0}=M_{H}*g/(4\pi{R_{p}}^{2})$. For a given $T$ and $P_{\mathrm{H}}^{0}$ , the atmospheric molecular composition above the magma ocean is computed as described in section 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 $\dot{M}\times 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 $\mathrm{H_{2}}$ or $\mathrm{H_{2}O}$ in the magma ocean as this will be the subject of a future study. We note that the dissolution of $\mathrm{H_{2}}$ and $\mathrm{H_{2}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 $\mathrm{H_{2}O}$ in the magma ocean (Sossi et al., 2023) indicates that it would be largely ($\geq$ 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 $\mathrm{H_{2}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.

Figure 14 displays the mass of hydrogen as a function of time for an Earth-like planet ($M_{p}=1\leavevmode\nobreak\ M_{\oplus}$, $R_{p}=1\leavevmode\nobreak\ R_{\oplus}$) located at 0.02 au, starting with different hydrogen content from $10^{-6}$ to $10^{-1}$ times the planet mass. The top panel shows that most H is lost within about 100 Kyr to 10 Myr, depending on $M_{H}^{0}/M_{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_{H}^{0}/M_{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 billion years 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_{H}^{0}/M_{p}=10^{-6}$. After 5 Gyr, the planet devoid of H has lost only $\sim$ $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 ($\eta\leavevmode\nobreak\ >\leavevmode\nobreak\ 0.1$, see section 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 section 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 ¿ 2600K and $M_{H}^{0}/M_{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_{H}^{0}/M_{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 ($\sim$ 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.