Atmospheric Chemistry of Secondary and Hybrid Atmospheres of Super Earths and Sub-Neptunes

Does the solar system mislead us? The planets of our solar system come in two flavors: gas/ice giants and rocky bodies. The gas and ice giants have primary atmospheres with compositions that largely reflect the chemistry of the primordial nebula out of which they formed. The terrestrial planets have secondary atmospheres with compositions that are largely determined by outgassing from their mantles and geochemical cycles (which cycle volatiles between the atmosphere, surface, and interior). Hybrid atmospheres, where the influences of both channels are comparable, do not exist in our solar system. Since we expect to hunt for biosignature gases in secondary atmospheres (e.g., Seager et al. 2013) that considerably deviate from Earth-like conditions, for example, surface pressure and mantle oxygen fugacity, it becomes imperative to develop a theory of outgassing chemistry that is capable of exploring the diverse set of physical and chemical conditions anticipated in exoplanets, since the outgassed species may be false positives for biosignatures. In other words, if biosignature hunting on exoplanets proceeds via astronomical remote sensing then biosignatures are spectral anomalies relative to a geochemical background.

Calculation of outgassing chemistry is a mature topic in the geochemical literature (e.g., French 1966; Ohmoto & Kerrick 1977; Holloway 1981, 1987; Connolly & Cesare 1993; Huizenga 2011). However, studies tend to focus on Earthcentric conditions (e.g., Iacono-Marziano et al. 2012; Gaillard & Scaillet 2014; Gaillard et al. 2022). Furthermore, the calculations are typically performed using computer codes that include dozens of species and hundreds, if not thousands, of chemical reactions (e.g., Schaefer et al. 2012; Fegley et al. 2016; Schaefer & Fegley 2017), which hinders understanding from first principles. In the astrophysical and climate literatures (Held 2005), it is common practice to construct a hierarchy of theoretical models such that the relative transparency of simpler models (e.g., French 1966) may be used to sweep parameter space and inform the detailed explorations of more complex models.

To the best of our knowledge, a simple, unified theoretical framework for computing the outgassing chemistry of both secondary and hybrid atmospheres does not exist in the exoplanet literature. Part of the challenge lies in curating a consistent set of thermodynamic definitions, notations, and experimental data in calculations of mixed-phase equilibrium chemistry. It is thus the goal of the present study to construct such a framework, which will inform future explorations of outgassing calculations for predicting the chemistry of secondary and hybrid atmospheres.

The seemingly clean dichotomy of planet types in the solar system has been broken by the discovery of super Earths and sub-Neptunes, which are exoplanets with a continuum of radii between 1 and 4 times the radius of Earth (Howard et al. 2012; Fressin et al. 2013; Petigura et al. 2013). Such exoplanets are common (see recent review by Bean et al. 2021). Exoplanets with radii above 1.5–1.6 R_⊕ (Rogers 2015; Fulton et al. 2017) appear to have puffy atmospheres dominated by hydrogen and/or helium (Owen 2019)—these are the sub-Neptunes. Below this critical radius, the bulk densities of these exoplanets are consistent with a predominantly rocky body—these are the super Earths. There is an ongoing debate on whether this population of exoplanets is sculpted by photoevaporation or core-powered mass loss (Rogers et al. 2021). Furthermore, Kite et al. (2019) proposed that considerable H₂ dissolution into magma oceans stalls sub-Neptune growth via H₂ accretion, leading to a "radius cliff" (Figure 1 in their study) and thus a paucity of sub-Neptunes with radii larger than thrice that of Earth's.

Zoom In Zoom Out Reset image size

The atmospheric chemistry of super Earths and sub-Neptunes remains largely unknown as current observations with the Hubble Space Telescope provide only loose constraints (e.g., Fisher & Heng 2018; Benneke et al. 2019; Mikal-Evans et al. 2021, and see Bean et al. 2021 for a summary). These exoplanets will be the subject of intense scrutiny by the James Webb Space Telescope (JWST) in the near future. Do super Earths have secondary atmospheres and are some of them hydrogen and/or helium dominated (Hu et al. 2015)? Are some of these atmospheres hybrid between primary and secondary atmospheres (Figure 1)? Models for predicting the chemistry of hybrid atmospheres are still in their infancy (e.g., Kite et al. 2020). As the anticipated flood of JWST data on these objects arrives, it is imperative that a hierarchy of theoretical models for predicting outgassing chemistry be available.

Almost without exception in the current literature, theoretical models of exoplanetary atmospheres assume ideal gases with constituent atoms and molecules that are ideally mixed. Departures from such ideal behaviors are quantified by the fugacity and activity coefficients, respectively. While ideal mixing may be a reasonable approximation at surface pressures ≲1000 bar (Figure 2), the assumption of an ideal gas may be inaccurate for simple molecules (Figure 3). Kite et al. (2019) has previously elucidated the relevance of the nonideal gas behavior of molecular hydrogen for sub-Neptunes. To the best of our current knowledge, we include these nonideal effects via the fugacity coefficient (for departures from an ideal gas) and activity coefficient (for departures from ideal mixing of gaseous components); our efforts are limited by the currently available experimental data on these coefficients. In particular, activity coefficients are specific to the mixture of molecules being considered and are difficult to obtain or nonexistent for arbitrary mixtures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Section 2 describes our formalism, which places secondary and hybrid atmospheres on an equal theoretical footing. It elucidates the governing equations, clarifies the thermodynamic definitions (including some confusion over the activity coefficient), states the numerical solution method, and discusses the relevant range of values for the atmospheric surface pressure and melt temperature. We examine both carbon–hydrogen–oxygen (C–H–O) and carbon–hydrogen–oxygen–nitrogen–sulfur (C–H–O–N–S) systems in turn. In Section 3, we review and curate the thermodynamic data that are critical for performing the outgassing calculations. Section 4 is devoted to a comprehensive exploration of parameter space for both secondary and hybrid atmospheres. In doing so, we discover the difficulty of producing methane-dominated atmospheres, which motivates a more in-depth investigation. In Section 5, we compare our current study to previous work, explore its observational implications, and describe opportunities for future work.

While none of the thermodynamics described here is novel, there are multiple definitions of the fugacity and activity present in the literature that are sometimes difficult to reconcile. Conceptually, the fugacity is the generalization of the pressure accounting for nonideal gas effects occurring at high pressures (e.g., Figure 3). The activity describes nonideal mixing in a gaseous system with multiple species (e.g., Figure 2). The assumptions of an ideal gas and ideal mixing is based on the simplification that the constituent gas molecules possess only kinetic, and not potential, energy (e.g., DeVoe 2020, p. 78). This assumes that intermolecular forces and their contribution to the Gibbs free energy are negligible (Holloway 1987). By definition, an ideal gas assumes ideal mixing of its components, implying that the fugacity and activity coefficients, which we will define shortly, are unity (e.g., Denbigh 1981, p. 114).

From the first law of thermodynamics, one may derive for an ideal gas with a single species (e.g., Denbigh 1981; Heng et al. 2016)

$$\begin{eqnarray}&&G\,=\,{G}_{0}+{ \mathcal R }T\mathrm{ln}{\rm{\Psi }},\end{eqnarray} \tag{ 1 }$$

where G is the specific Gibbs free energy (Gibbs free energy per unit mass), G₀ ≡ G(P₀), P₀ is the reference pressure (usually set to 1 atm or 1 bar), ${ \mathcal R }$ is the specific gas constant, and T is the temperature. The task is to generalize Ψ = P/P₀ for gas mixtures with nonideal gas behavior.

As far as possible, we respect the established notation in the geochemical literature. The fugacity is commonly denoted as f. Holloway (1977) and DeVoe (2020) use ϕ and γ for the fugacity coefficient and activity coefficient, respectively. In the current study, we follow this convention.

2.1.1. Fugacity

For a pure gas (one species),

$$\begin{eqnarray}&&G\,=\,{G}_{0}+{\int }_{{P}_{0}}^{P}\displaystyle \frac{1}{\rho }\,{dP},\end{eqnarray} \tag{ 2 }$$

where ρ is the mass density and 1/ρ is the volume per unit mass. We call the second term after the equality the "volume integral term." For an ideal gas, $\rho =P/{ \mathcal R }T$ and one obtains Equation (1) and Ψ = P/P₀. For nonideal equations of state, it is not possible to write down a general closed-form expression for Ψ. Instead, one defines a quantity known as the fugacity (Lewis 1901),

$$\begin{eqnarray}&&f\,=\,{\phi }_{i}P,\end{eqnarray} \tag{ 3 }$$

where ϕ_i is the fugacity coefficient of the pure gas composed of species i. Fugacity coefficients are specific to the chemical species being considered. For an ideal gas, ϕ_i = 1.

The generalization of Equation (1) is done by analogy with no underlying mathematical rigor. Essentially, one does the ad hoc substitution (Equation (29) of Lewis 1901),

$$\begin{eqnarray}&&\displaystyle \frac{P}{{P}_{0}}\to \displaystyle \frac{f}{{f}_{0}},\end{eqnarray} \tag{ 4 }$$

where f₀ = ϕ_i,0 P₀ and ϕ_i,0 ≡ ϕ_i (P₀).

2.1.2. Activity

At a pressure P and temperature T, the Gibbs free energy of a gas component i in a gaseous mixture is (Denbigh 1981, p. 115),

$$\begin{eqnarray}&&{G}_{\mathrm{mix}}=G+{ \mathcal R }T\mathrm{ln}{X}_{i},\end{eqnarray} \tag{ 5 }$$

where G is given by Equation (2) and X_i is the volume mixing ratio (relative abundance by number) of the ith gaseous species. One now defines the activity as the generalization of X_i ,

$$\begin{eqnarray}&&{a}_{i}={\gamma }_{i}{X}_{i},\end{eqnarray} \tag{ 6 }$$

where γ_i is the activity coefficient of the ith gaseous species. Activity coefficients are specific to the chemical species and mixture of species being considered. For ideal mixing, γ_i = 1.

Again, by analogy, the specific Gibbs free energy of the ith gaseous species in a gas mixture with nonideal mixing of its constituents is (Denbigh 1981, p. 287)

$$\begin{eqnarray}&&{G}_{\mathrm{mix}}=G+{ \mathcal R }T\mathrm{ln}{a}_{i}.\end{eqnarray} \tag{ 7 }$$

2.1.3. Equilibrium Constant

As we will see later in Sections 2.2 and 2.3, the equilibrium constant K_eq is constructed from the Ψ of the reactants and products in a chemical reaction. It is related to the Gibbs free energy of the reaction ΔG_r by (DeVoe 2020, p. 352)

$$\begin{eqnarray}&&{K}_{\mathrm{eq}}\,=\,{e}^{-\tfrac{{\rm{\Delta }}{G}_{r}}{{ \mathcal R }T}}.\end{eqnarray} \tag{ 8 }$$

For gases, one may write

$$\begin{eqnarray}&&{G}_{\mathrm{mix}}={G}_{0}+{ \mathcal R }T\mathrm{ln}\left(\displaystyle \frac{{a}_{i}f}{{f}_{0}}\right),\end{eqnarray} \tag{ 9 }$$

such that

$$\begin{eqnarray}&&{\rm{\Psi }}=\displaystyle \frac{{a}_{i}f}{{f}_{0}}=\displaystyle \frac{{f}_{i}}{{P}_{0}}.\end{eqnarray} \tag{ 10 }$$

We have written the fugacity of the ith gaseous species as f_i = a_i f. Since P₀ is typically chosen to be 1 bar, it is reasonable to set ϕ_i,0 = 1 (a gas at 1 bar behaves like an ideal gas) and thus f₀ = ϕ_i,0 P₀ = P₀. Given that pure-gas fugacity f is a generalization of total pressure P and that activity a_i is a generalization of volume mixing ratio X_i , the fugacity of the ith gaseous species f_i = a_i f can be deemed a generalization of partial pressure. It is worth emphasizing that Ψ is not the activity. Some confusion stems from the fact that the reduced expression of Ψ = f_i /P₀ is sometimes ⁶ referred to as the activity (see also Denbigh 1981, p. 287). This apparent ambiguity originates from the mathematical freedom to combine different terms into the argument of the natural logarithm.

For a species in its solid phase, the fugacity is not explicitly stated and the Gibbs free energy is instead,

$$\begin{eqnarray}&&G\,=\,{G}^{{\prime} }+{ \mathcal R }T\mathrm{ln}{a}_{s},\end{eqnarray} \tag{ 11 }$$

where the quantity ${G}^{{\prime} }$ depends on the pressure P (and is not referenced to P₀, unlike for gases) and the activity associated with the solid is a_s . We simply choose Ψ = a_s while paying attention to the definition of ${G}^{{\prime} }$. The exact expression for ${G}^{{\prime} }$ is provided in Section 3. In the current study, we consider only the activity associated with graphite.

Consider the standard net chemical reaction for converting methane to carbon monoxide (e.g., Burrows & Sharp 1999; Lodders & Fegley 2002; Heng & Lyons 2016),

$$\begin{eqnarray}&&{\mathrm{CH}}_{4}+{{\rm{H}}}_{2}{\rm{O}}\leftrightarrows \mathrm{CO}+3{{\rm{H}}}_{2}.\end{eqnarray} \tag{ 12 }$$

In the limit of an ideal gas (ϕ_i = 1) with ideal mixing (γ_i = 1), we have Ψ = P_i /P₀.

We may write down the Gibbs free energies associated with the products and reactants,

$$\begin{eqnarray}\begin{array}{rcl}{G}_{{\mathrm{CH}}_{4}} & = & {G}_{{\mathrm{CH}}_{4},0}+{ \mathcal R }T\mathrm{ln}{{\rm{\Psi }}}_{{\mathrm{CH}}_{4}},\\ {G}_{{{\rm{H}}}_{2}{\rm{O}}} & = & {G}_{{{\rm{H}}}_{2}{\rm{O}},0}+{ \mathcal R }T\mathrm{ln}{{\rm{\Psi }}}_{{{\rm{H}}}_{2}{\rm{O}}},\\ {G}_{\mathrm{CO}} & = & {G}_{\mathrm{CO},0}+{ \mathcal R }T\mathrm{ln}{{\rm{\Psi }}}_{\mathrm{CO}},\\ {G}_{{{\rm{H}}}_{2}} & = & {G}_{{{\rm{H}}}_{2},0}+{ \mathcal R }T\mathrm{ln}{{\rm{\Psi }}}_{{{\rm{H}}}_{2}}.\end{array}\end{eqnarray} \tag{ 13 }$$

The differences in Gibbs free energies are

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}G & = & {G}_{\mathrm{CO}}+3{G}_{{{\rm{H}}}_{2}}-{G}_{{\mathrm{CH}}_{4}}-{G}_{{{\rm{H}}}_{2}{\rm{O}}},\\ {\rm{\Delta }}{G}_{0} & = & {G}_{\mathrm{CO},0}+3{G}_{{{\rm{H}}}_{2},0}-{G}_{{\mathrm{CH}}_{4},0}-{G}_{{{\rm{H}}}_{2}{\rm{O}},0}.\end{array}\end{eqnarray} \tag{ 14 }$$

It follows that

$$\begin{eqnarray}&&{\rm{\Delta }}G-{\rm{\Delta }}{G}_{0}={ \mathcal R }T\mathrm{ln}\left(\displaystyle \frac{{{\rm{\Psi }}}_{\mathrm{CO}}{{\rm{\Psi }}}_{{{\rm{H}}}_{2}}^{3}}{{{\rm{\Psi }}}_{{\mathrm{CH}}_{4}}{{\rm{\Psi }}}_{{{\rm{H}}}_{2}{\rm{O}}}}\right).\end{eqnarray} \tag{ 15 }$$

Following the reasoning by Heng et al. (2016), the system attains chemical equilibrium by adjusting to ΔG = 0; ΔG₀ is interpreted as the change of Gibbs free energy of formation at the reference pressure P₀. By comparing the preceding equation with Equation (8), the equilibrium constant is

$$\begin{eqnarray}&&{K}_{\mathrm{eq}}=\displaystyle \frac{{{\rm{\Psi }}}_{\mathrm{CO}}{{\rm{\Psi }}}_{{{\rm{H}}}_{2}}^{3}}{{{\rm{\Psi }}}_{{\mathrm{CH}}_{4}}{{\rm{\Psi }}}_{{{\rm{H}}}_{2}{\rm{O}}}},\end{eqnarray} \tag{ 16 }$$

in agreement with Equations (6) and (26) of Heng & Lyons (2016). The Gibbs free energy of the reaction is ΔG_r = ΔG₀. Generally, the approach of constructing K_eq from Ψ is consistent with Equation (13) of Heng et al. (2016).

French (1966) considered the system of net chemical reactions,

$$\begin{eqnarray}\begin{array}{rcl}{\rm{C}}+{{\rm{O}}}_{2} & \leftrightarrows & {\mathrm{CO}}_{2},\\ {\rm{C}}+0.5{{\rm{O}}}_{2} & \leftrightarrows & \mathrm{CO},\\ {\rm{C}}+2{{\rm{H}}}_{2} & \leftrightarrows & {\mathrm{CH}}_{4},\\ {{\rm{H}}}_{2}+0.5{{\rm{O}}}_{2} & \leftrightarrows & {{\rm{H}}}_{2}{\rm{O}}.\end{array}\end{eqnarray} \tag{ 17 }$$

In this system, carbon is present in its solid phase as graphite. French (1966) assumed ideal mixing (γ_i = 1).

If we focus on the first net reaction, the Gibbs free energies are

$$\begin{eqnarray}\begin{array}{rcl}{G}_{{\rm{C}}} & = & {G}_{{\rm{C}}}^{{\prime} }+{ \mathcal R }T\mathrm{ln}{a}_{{\rm{C}}},\\ {G}_{{{\rm{O}}}_{2}} & = & {G}_{{{\rm{O}}}_{2},0}+{ \mathcal R }T\mathrm{ln}\left(\displaystyle \frac{{f}_{{{\rm{O}}}_{2}}}{{P}_{0}}\right),\\ {G}_{{\mathrm{CO}}_{2}} & = & {G}_{{\mathrm{CO}}_{2},0}+{ \mathcal R }T\mathrm{ln}\left(\displaystyle \frac{{f}_{{\mathrm{CO}}_{2}}}{{P}_{0}}\right).\end{array}\end{eqnarray} \tag{ 18 }$$

The critical detail is that ${G}_{{\rm{C}}}^{{\prime} }(P,T)$ depends on the pressure P, whereas ${G}_{{{\rm{O}}}_{2},0}({P}_{0})$ and ${G}_{{\mathrm{CO}}_{2},0}({P}_{0})$ depend only on the reference pressure P₀.

The difference in Gibbs free energies is

$$\begin{eqnarray}&&{\rm{\Delta }}G={G}_{{\mathrm{CO}}_{2}}-{G}_{{{\rm{O}}}_{2}}-{G}_{{\rm{C}}}.\end{eqnarray} \tag{ 19 }$$

If we again argue that ΔG = 0 when chemical equilibrium is attained (Heng et al. 2016), then the expression relating the equilibrium constant to the Gibbs free energy of the reaction follows,

$$\begin{eqnarray}&&{K}_{\mathrm{eq},1}=\displaystyle \frac{{f}_{{\mathrm{CO}}_{2}}}{{f}_{{{\rm{O}}}_{2}}{a}_{{\rm{C}}}}={e}^{-\tfrac{{\rm{\Delta }}{G}_{r}}{{ \mathcal R }T}},\end{eqnarray} \tag{ 20 }$$

where the Gibbs free energy of the reaction is

$$\begin{eqnarray}&&{\rm{\Delta }}{G}_{r}={G}_{{\mathrm{CO}}_{2},0}-{G}_{{{\rm{O}}}_{2},0}-{G}_{{\rm{C}}}^{{\prime} }.\end{eqnarray} \tag{ 21 }$$

In addition to the usual Gibbs free energies of formation associated with the gaseous species, which are defined at the reference pressure P₀, one considers the Gibbs free energies of formation associated with graphite, which are defined at the pressure of interest P. If we further set a_C = 1, Equation (1) of French (1966) is recovered.

For the other three net chemical reactions in Equation (1), the equilibrium constants are

$$\begin{eqnarray}\begin{array}{rcl}{K}_{\mathrm{eq},2} & = & \displaystyle \frac{{f}_{\mathrm{CO}}}{{\left({f}_{{{\rm{O}}}_{2}}{P}_{0}\right)}^{1/2}{a}_{{\rm{C}}}},\\ {K}_{\mathrm{eq},3} & = & \displaystyle \frac{{f}_{{\mathrm{CH}}_{4}}{P}_{0}}{{f}_{{{\rm{H}}}_{2}}^{2}{a}_{{\rm{C}}}},\\ {K}_{\mathrm{eq},4} & = & \displaystyle \frac{{f}_{{{\rm{H}}}_{2}{\rm{O}}}{P}_{0}^{1/2}}{{f}_{{{\rm{H}}}_{2}}{f}_{{{\rm{O}}}_{2}}^{1/2}}.\end{array}\end{eqnarray} \tag{ 22 }$$

If one again sets a_C = 1 and omits writing the factors of P₀ explicitly (by setting them to unity), then one recovers Equations (2) to (4) of French (1966).

The oxygen fugacity (${f}_{{{\rm{O}}}_{2}}$) quantifies how reduced or oxidized the mantle of an exoplanet is. It may be interpreted as the equivalent partial pressure of oxygen even if the oxygen is locked up in solids in the form of minerals. A representative reaction for iron in the core of the Earth and iron oxide (wüstite) in its mantle is (Wade & Wood 2005)

$$\begin{eqnarray}&&\mathrm{Fe}+0.5{{\rm{O}}}_{2}\leftrightarrows \mathrm{FeO}.\end{eqnarray} \tag{ 23 }$$

If factors of P₀ are ignored, then the oxygen fugacity buffered by this reaction has the following equilibrium constant,

$$\begin{eqnarray}&&{K}_{\mathrm{eq}}=\displaystyle \frac{{a}_{\mathrm{FeO}}}{{a}_{\mathrm{Fe}}{f}_{{{\rm{O}}}_{2}}^{1/2}}.\end{eqnarray} \tag{ 24 }$$

By using Equation (8) and performing a change of base of the logarithm (from base 2 to 10), we obtain

$$\begin{eqnarray}&&\mathrm{log}{f}_{{{\rm{O}}}_{2}}=\displaystyle \frac{2}{\mathrm{ln}10}\displaystyle \frac{{\rm{\Delta }}{G}_{r}}{{ \mathcal R }T}+2\mathrm{log}\left(\displaystyle \frac{{a}_{\mathrm{FeO}}}{{a}_{\mathrm{Fe}}}\right).\end{eqnarray} \tag{ 25 }$$

An idealization known as the iron–wüstite (IW) buffer is defined when iron and wüstite occur in their pure forms, such that their activities may be set to unity (Wade & Wood 2005),

$$\begin{eqnarray}&&\mathrm{IW}\equiv \mathrm{log}{f}_{{{\rm{O}}}_{2}}^{\mathrm{IW}}=\displaystyle \frac{2}{\mathrm{ln}10}\displaystyle \frac{{\rm{\Delta }}{G}_{r}}{{ \mathcal R }T},\end{eqnarray} \tag{ 26 }$$

where we follow the convention in the geochemical literature of abbreviating $\mathrm{log}{f}_{{{\rm{O}}}_{2}}^{\mathrm{IW}}$ as "IW." The oxygen fugacity may span 50 orders of magnitude in value (Frost 1991), and it is both inconvenient and nonintuitive to quote its absolute value (in bar). Rather, geochemists favor quoting the oxygen fugacity relative to simple ⁷ buffers such as IW. From a computational standpoint, IW is a convenient reference point that allows $\mathrm{log}{f}_{{{\rm{O}}}_{2}}$ to be reported in terms of intuitive, order-of-unity numbers.

It is possible to use Equation (25) to estimate the oxygen fugacity associated with the core–mantle boundary or core formation of the Earth. By interpreting the activities as relative abundances by number, we may approximate a_FeO ≈ 0.08 (mantle) and a_Fe ≈ 0.8 (core) (Wade & Wood 2005) such that $\mathrm{log}({a}_{\mathrm{FeO}}/{a}_{\mathrm{Fe}})=-1$. Equation (25) then yields $\mathrm{log}{f}_{{{\rm{O}}}_{2}}=2{\rm{\Delta }}{G}_{r}/{ \mathcal R }T\mathrm{ln}10-2$. Using Equation (26), we obtain $\mathrm{log}{f}_{{{\rm{O}}}_{2}}=\mathrm{IW}-2$ for the core–mantle boundary of the Earth or the mantle during Earth's core formation (Wade & Wood 2005). Physically, the iron oxide in the mantle and the iron in the core jointly buffer the oxygen fugacity. It is worth noting that the oxygen fugacity of the Earth varies both temporally and spatially. For modern Earth, the upper mantle has an estimated oxygen fugacity of IW + 3 to IW + 7 (Frost & McCammon 2008). Another buffer that is commonly used is the FMQ buffer, where F, M, and Q represent fayalite (Fe₂SiO₄), magnetite (Fe₃O₄), and quartz (SiO₂), respectively (Frost 1991). IW + 5 converts to about FMQ (Frost 1991; Frost & McCammon 2008). Thus, even for modern Earth the oxygen fugacity varies by about 9 orders of magnitude from the core–mantle boundary to the upper mantle. Temporally, it is believed that the oxygen fugacity of the upper mantle increased to its current value by at least 3.6 billion years ago (Delano 2001).

Since ΔG_r generally depends on temperature and pressure, it is unsurprising that IW has the same dependences. To obtain an absolute value for IW, Ballhaus et al. (1991) provides a convenient, empirical fitting function for IW based on the experimental data from O'Neill (1987),

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{IW} & = & 14.07-28784{\left(\displaystyle \frac{T}{1\,{\rm{K}}}\right)}^{-1}-2.04\mathrm{log}\left(\displaystyle \frac{T}{1\,{\rm{K}}}\right)\\ & & +0.053\left(\displaystyle \frac{P}{1\mathrm{bar}}\right){\left(\displaystyle \frac{T}{1\,{\rm{K}}}\right)}^{-1}+3\times {10}^{-6}\left(\displaystyle \frac{P}{1\mathrm{bar}}\right),\end{array}\end{eqnarray} \tag{ 27 }$$

which is valid for 900 ≤ T/K ≤ 1420. The pressure dependence is derived by volume integration. It is understood that IW, as represented above, has physical units of the logarithm (base 10) of pressure in bar.

Analogously, a sulfur fugacity buffer may be written down based on pyrrhotite (FeS) and pyrite (FeS₂) (Ferry & Baumgartner 1987),

$$\begin{eqnarray}&&\mathrm{FeS}+0.5{{\rm{S}}}_{2}\leftrightarrows {\mathrm{FeS}}_{2},\end{eqnarray} \tag{ 28 }$$

which we abbreviate as "PP" for convenience. Froese & Gunter (1976) provides an empirical fitting function based on the experimental data from Toulmin & Barton (1964),

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{PP}\equiv \mathrm{log}{f}_{{{\rm{S}}}_{2}}^{\mathrm{PP}} & = & -16073{\left(\displaystyle \frac{T}{1\,{\rm{K}}}\right)}^{-1}+15.74\\ & & +0.06\left[\left(\displaystyle \frac{P}{1\mathrm{bar}}\right)-1\right]{\left(\displaystyle \frac{T}{1\,{\rm{K}}}\right)}^{-1},\end{array}\end{eqnarray} \tag{ 29 }$$

which is valid for 598 ≤ T/K ≤ 1016. The pressure dependence is derived by volume integration. It is again understood that PP, as represented above, has physical units of the logarithm (base 10) of pressure in bar.

Compared to the oxygen fugacity, it is worth noting that data on the sulfur fugacity in different geological settings are sparse, although estimates span an impressive range of values from $\mathrm{log}{f}_{{{\rm{S}}}_{2}}\approx -20$ (PP − 8) in submarine hydrothermal vents (Keith et al. 2014), $\mathrm{log}{f}_{{{\rm{S}}}_{2}}\approx -5$ to −2 (PP − 3 to PP) associated with metamorphic degassing (Poulson & Ohmoto 1989), to $\mathrm{log}{f}_{{{\rm{S}}}_{2}}\approx 4$ (PP + 1) at high pressures (∼10⁴ bar) associated with the sulfide saturation conditions of basaltic melts (Jégo & Dasgupta 2014, sample G280).

The sulfur and oxygen fugacities may be related if the FeO content of the melt is known (Bockrath et al. 2004), a caveat that is also noted in Section 3.1 of Liggins et al. (2020). Since we do not explicitly model the melt composition, this additional layer of complexity is left to future work and we regard ${f}_{{{\rm{S}}}_{2}}$ and ${f}_{{{\rm{O}}}_{2}}$ as independent input parameters.

2.5.1. Net Chemical Reactions

Using the French (1966) model as a baseline, we consider the following expanded set of net chemical reactions (Holloway 1981; Burrows & Sharp 1999; Lodders & Fegley 2002; Moses et al. 2011; Gaillard & Scaillet 2014),

$$\begin{eqnarray}\begin{array}{rcl}{\rm{C}}+{{\rm{O}}}_{2} & \leftrightarrows & {\mathrm{CO}}_{2},\\ {\rm{C}}+0.5{{\rm{O}}}_{2} & \leftrightarrows & \mathrm{CO},\\ {\rm{C}}+2{{\rm{H}}}_{2} & \leftrightarrows & {\mathrm{CH}}_{4},\\ {{\rm{H}}}_{2}+0.5{{\rm{O}}}_{2} & \leftrightarrows & {{\rm{H}}}_{2}{\rm{O}},\\ {{\rm{N}}}_{2}+3{{\rm{H}}}_{2} & \leftrightarrows & 2{\mathrm{NH}}_{3},\\ {\mathrm{NH}}_{3}+{\mathrm{CH}}_{4} & \leftrightarrows & \mathrm{HCN}+3{{\rm{H}}}_{2},\\ 0.5{{\rm{S}}}_{2}+{{\rm{O}}}_{2} & \leftrightarrows & {\mathrm{SO}}_{2},\\ 0.5{{\rm{S}}}_{2}+{{\rm{H}}}_{2}{\rm{O}} & \leftrightarrows & {{\rm{H}}}_{2}{\rm{S}}+0.5{{\rm{O}}}_{2}.\end{array}\end{eqnarray} \tag{ 30 }$$

Similar to the oxygen fugacity, we treat the sulfur fugacity ${f}_{{{\rm{S}}}_{2}}$ as an input parameter (e.g., Holloway 1981). The nitrogen content of the atmosphere is parameterized via the partial pressure of nitrogen ${P}_{{{\rm{N}}}_{2}}$. Assuming these two additional input parameters allows the number of unknowns and equations to be equal.

2.5.2. Equilibrium Constants

Modern thermodynamic data allow us to consider nonideal mixing (γ_i ≠ 1) and nonideal gas behavior (ϕ_i ≠ 1), unlike for French (1966). For compactness of notation, we define

$$\begin{eqnarray}&&{\alpha }_{i}\equiv {\gamma }_{i}{\phi }_{i}.\end{eqnarray} \tag{ 31 }$$

This motivates us to write down more generalized expressions for the equilibrium constants,

$$\begin{eqnarray}\begin{array}{rcl}{K}_{\mathrm{eq},1} & = & \displaystyle \frac{{\alpha }_{{\mathrm{CO}}_{2}}{P}_{{\mathrm{CO}}_{2}}}{{\gamma }_{{{\rm{O}}}_{2}}{f}_{{{\rm{O}}}_{2}}{a}_{{\rm{C}}}},\\ {K}_{\mathrm{eq},2} & = & \displaystyle \frac{{\alpha }_{\mathrm{CO}}{P}_{\mathrm{CO}}}{{\left({\gamma }_{{{\rm{O}}}_{2}}{f}_{{{\rm{O}}}_{2}}{P}_{0}\right)}^{1/2}{a}_{{\rm{C}}}},\\ {K}_{\mathrm{eq},3} & = & \displaystyle \frac{{\alpha }_{{\mathrm{CH}}_{4}}{P}_{{\mathrm{CH}}_{4}}{P}_{0}}{{\left({\alpha }_{{{\rm{H}}}_{2}}{P}_{{{\rm{H}}}_{2}}\right)}^{2}{a}_{{\rm{C}}}},\\ {K}_{\mathrm{eq},4} & = & \displaystyle \frac{{\alpha }_{{{\rm{H}}}_{2}{\rm{O}}}{P}_{{{\rm{H}}}_{2}{\rm{O}}}{P}_{0}^{1/2}}{{\alpha }_{{{\rm{H}}}_{2}}{P}_{{{\rm{H}}}_{2}}{\left({\gamma }_{{{\rm{O}}}_{2}}{f}_{{{\rm{O}}}_{2}}\right)}^{1/2}}.\\ {K}_{\mathrm{eq},5} & = & \displaystyle \frac{{\left({\alpha }_{{\mathrm{NH}}_{3}}{P}_{{\mathrm{NH}}_{3}}{P}_{0}\right)}^{2}}{{\alpha }_{{{\rm{N}}}_{2}}{P}_{{{\rm{N}}}_{2}}{\left({\alpha }_{{{\rm{H}}}_{2}}{P}_{{{\rm{H}}}_{2}}\right)}^{3}},\\ {K}_{\mathrm{eq},6} & = & \displaystyle \frac{{\alpha }_{\mathrm{HCN}}{P}_{\mathrm{HCN}}{\left({\alpha }_{{{\rm{H}}}_{2}}{P}_{{{\rm{H}}}_{2}}\right)}^{3}}{{\alpha }_{{\mathrm{NH}}_{3}}{\alpha }_{{\mathrm{CH}}_{4}}{P}_{{\mathrm{NH}}_{3}}{P}_{{\mathrm{CH}}_{4}}{P}_{0}^{2}},\\ {K}_{\mathrm{eq},7} & = & \displaystyle \frac{{\alpha }_{{\mathrm{SO}}_{2}}{P}_{{\mathrm{SO}}_{2}}{P}_{0}^{1/2}}{{\gamma }_{{{\rm{O}}}_{2}}{f}_{{{\rm{O}}}_{2}}{\left({\gamma }_{{{\rm{S}}}_{2}}{f}_{{{\rm{S}}}_{2}}\right)}^{1/2}},\\ {K}_{\mathrm{eq},8} & = & \displaystyle \frac{{\alpha }_{{{\rm{H}}}_{2}{\rm{S}}}{P}_{{{\rm{H}}}_{2}{\rm{S}}}{\left({\gamma }_{{{\rm{O}}}_{2}}{f}_{{{\rm{O}}}_{2}}\right)}^{1/2}}{{\alpha }_{{{\rm{H}}}_{2}{\rm{O}}}{P}_{{{\rm{H}}}_{2}{\rm{O}}}{\left({\gamma }_{{{\rm{S}}}_{2}}{f}_{{{\rm{S}}}_{2}}\right)}^{1/2}}.\end{array}\end{eqnarray} \tag{ 32 }$$

We note the convention in the geochemical literature (the use of Equation (10); e.g., Ferry & Baumgartner 1987) is to absorb the activity coefficient (γ_i ) into the fugacity (f_i ). By departing from this practice, we separate the effects of nonideal mixing (expressed through γ_i ) from nonideal gas equations of state (expressed through ϕ_i ) by explicitly writing f_i = ϕ_i P_i .

2.5.3. Partial and Total Pressures

By rewriting the expressions for the equilibrium constants, one may write down expressions for the partial pressures of the various species,

$$\begin{eqnarray}\begin{array}{rcl}{P}_{{\mathrm{CO}}_{2}} & = & \displaystyle \frac{{a}_{{\rm{C}}}{K}_{\mathrm{eq},1}{\gamma }_{{{\rm{O}}}_{2}}{f}_{{{\rm{O}}}_{2}}}{{\alpha }_{{\mathrm{CO}}_{2}}},\\ {P}_{\mathrm{CO}} & = & \displaystyle \frac{{a}_{{\rm{C}}}{K}_{\mathrm{eq},2}{\left({\gamma }_{{{\rm{O}}}_{2}}{f}_{{{\rm{O}}}_{2}}{P}_{0}\right)}^{1/2}}{{\alpha }_{\mathrm{CO}}},\\ {P}_{{\mathrm{CH}}_{4}} & = & \displaystyle \frac{{a}_{{\rm{C}}}{K}_{\mathrm{eq},3}{\left({\alpha }_{{{\rm{H}}}_{2}}{P}_{{{\rm{H}}}_{2}}\right)}^{2}}{{\alpha }_{{\mathrm{CH}}_{4}}{P}_{0}},\\ {P}_{{{\rm{H}}}_{2}{\rm{O}}} & = & \displaystyle \frac{{K}_{\mathrm{eq},4}{\alpha }_{{{\rm{H}}}_{2}}{P}_{{{\rm{H}}}_{2}}{\left({\gamma }_{{{\rm{O}}}_{2}}{f}_{{{\rm{O}}}_{2}}\right)}^{1/2}}{{\alpha }_{{{\rm{H}}}_{2}{\rm{O}}}{P}_{0}^{1/2}},\\ {P}_{{\mathrm{NH}}_{3}} & = & \displaystyle \frac{{K}_{\mathrm{eq},5}^{1/2}{\left({\alpha }_{{{\rm{N}}}_{2}}{P}_{{{\rm{N}}}_{2}}\right)}^{1/2}{\left({\alpha }_{{{\rm{H}}}_{2}}{P}_{{{\rm{H}}}_{2}}\right)}^{3/2}}{{\alpha }_{{\mathrm{NH}}_{3}}{P}_{0}},\\ {P}_{\mathrm{HCN}} & = & \displaystyle \frac{{K}_{\mathrm{eq},6}{\alpha }_{{\mathrm{NH}}_{3}}{\alpha }_{{\mathrm{CH}}_{4}}{P}_{{\mathrm{NH}}_{3}}{P}_{{\mathrm{CH}}_{4}}{P}_{0}^{2}}{{\alpha }_{\mathrm{HCN}}{\left({\alpha }_{{{\rm{H}}}_{2}}{P}_{{{\rm{H}}}_{2}}\right)}^{3}},\\ {P}_{{\mathrm{SO}}_{2}} & = & \displaystyle \frac{{K}_{\mathrm{eq},7}{\gamma }_{{{\rm{O}}}_{2}}{f}_{{{\rm{O}}}_{2}}{\left({\gamma }_{{{\rm{S}}}_{2}}{f}_{{{\rm{S}}}_{2}}\right)}^{1/2}}{{\alpha }_{{\mathrm{SO}}_{2}}{P}_{0}^{1/2}},\\ {P}_{{{\rm{H}}}_{2}{\rm{S}}} & = & \displaystyle \frac{{K}_{\mathrm{eq},8}{\alpha }_{{{\rm{H}}}_{2}{\rm{O}}}{P}_{{{\rm{H}}}_{2}{\rm{O}}}{\left({\gamma }_{{{\rm{S}}}_{2}}{f}_{{{\rm{S}}}_{2}}\right)}^{1/2}}{{\alpha }_{{{\rm{H}}}_{2}{\rm{S}}}{\left({\gamma }_{{{\rm{O}}}_{2}}{f}_{{{\rm{O}}}_{2}}\right)}^{1/2}},\end{array}\end{eqnarray} \tag{ 33 }$$

which depend on the parameters ${f}_{{{\rm{O}}}_{2}}$, ${f}_{{{\rm{S}}}_{2}}$, and ${P}_{{{\rm{N}}}_{2}}$. We note that ${P}_{\mathrm{HCN}}\propto {P}_{{{\rm{H}}}_{2}}^{1/2}$.

The total pressure is given by

$$\begin{eqnarray}\begin{array}{rcl}P & = & {P}_{{{\rm{H}}}_{2}}+{P}_{\mathrm{He}}+{P}_{{{\rm{O}}}_{2}}+{P}_{{{\rm{S}}}_{2}}+{P}_{{{\rm{N}}}_{2}}+{P}_{{\mathrm{CO}}_{2}}+{P}_{\mathrm{CO}}\\ & & +{P}_{{\mathrm{CH}}_{4}}+{P}_{{{\rm{H}}}_{2}{\rm{O}}}+{P}_{{\mathrm{NH}}_{3}}+{P}_{\mathrm{HCN}}+{P}_{{\mathrm{SO}}_{2}}+{P}_{{{\rm{H}}}_{2}{\rm{S}}}.\end{array}\end{eqnarray} \tag{ 34 }$$

2.5.4. Solution Method

The input pressure P is interpreted as the surface pressure of the atmosphere (Gaillard & Scaillet 2014). We consider P = 1 mbar to 10⁴ bar. The lower limit is arbitrary and empirically inspired by the Martian atmosphere. The upper limit is motivated by the estimate that a hydrogen-dominated atmosphere with ∼1% of the mass of the exoplanet is massive enough to double the radius of its core (Owen 2019). Denoting its radius by R and its surface gravity by g, the mass of the atmosphere of an exoplanet is M_atm = 4π R² P/g. It follows that the ratio of the atmospheric mass to the total mass of the exoplanet is

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{\mathrm{atm}}}{M}=0.026\,\left(\displaystyle \frac{P}{{10}^{4}\mathrm{bar}}\right){\left(\displaystyle \frac{R}{1.6\,{R}_{\oplus }}\displaystyle \frac{\rho }{2.0{\rm{g}}\ {\mathrm{cm}}^{-3}}\right)}^{-2}.\end{eqnarray} \tag{ 35 }$$

It is thus not unreasonable to regard sub-Neptunes as exoplanets with surface pressures ∼10⁴ bar.

The input temperature T is interpreted as that of the melt from which the secondary atmospheres are outgassed. We use the term "melt" to refer to liquids at both low and high temperatures; the term "magma" refers to a high-temperature melt. Its range of values is bracketed by two plausible melting scenarios associated with rocky exoplanets. The first scenario is a fully molten magma ocean (Hirschmann 2012; Keppler & Golabek 2019; Sossi et al. 2020; Bower et al. 2022; Gaillard et al. 2022). The second scenario considers volcanic outgassing from a largely solidified body (Symonds & Reed 1993; Gaillard & Scaillet 2014; Liggins et al. 2020; Höning et al. 2021). The lower temperature bound that we adopt is T = 873 K, which corresponds roughly to the minimum partial melting temperature (solidus) of wet granite under the pressure range of 1 mbar to 10⁴ bar (Huang & Wyllie 1973). To fully melt a "dry" peridotite, the temperature (liquidus) is about 1973 K at a pressure of 1 bar. Using Figure 7(a) of Takahashi et al. (1993), we determine that this temperature value will not change much at 1 mbar and scale its value to about 2073 K at 10⁴ bar. It is known that volatiles in a melt can reduce the liquidus by about 180–300 K (Hort 1998), which motivates us to set the upper temperature bound for melts to 1873 K to accommodate the scenario of a planetary-scale magma ocean. Temperatures higher than 1873 K will likely involve nonnegligible silicate vaporization (e.g., Schaefer & Fegley 2003; Schaefer et al. 2012; Herbort et al. 2020), which is beyond the scope of this study.

Motivated by the possibility that melt compositions may be different from peridotitic (e.g., basalt), Gaillard et al. (2022) chose T = 1773 K. Gaillard & Scaillet (2014) chose T = 1573 K, which we adopt and approximate as 1600 K to emphasize the somewhat ad hoc choice of melt temperature in the absence of a full interior-mantle model. In plots where we have to fix the temperature and vary other parameters, we adopt T = 1600 K, an intermediate value between the lower (873 K) and upper (1873 K) bounds of melt temperatures.

For both gaseous and solid phases, the Gibbs free energy is given by Equation (A1) of Powell (1978),

$$\begin{eqnarray}&&{G}_{i}={G}_{i,0}+{\int }_{{P}_{0}}^{P}{\left.\displaystyle \frac{\partial {G}_{i}}{\partial P}\right|}_{T}\,{dP},\end{eqnarray} \tag{ 36 }$$

where G_i = G_i (T, P) and G_i,0 ≡ G_i (T, P₀). For convenience, we write the integrand as V_i . The Gibbs free energy at the reference pressure P₀ may be further expanded as G_i,0 = H_i (T, P₀) − TS_i (T, P₀), where H_i (T, P₀) and S_i (T, P₀) are the enthalpy and entropy at the reference pressure, respectively. If we further expand H_i (T, P₀) and S_i (T, P₀) about a reference temperature T₀, then we obtain

$$\begin{eqnarray}\begin{array}{rcl}{G}_{i,0} & = & {H}_{i,0}-{{TS}}_{i,0}+{\int }_{{T}_{0}}^{T}\displaystyle \frac{\partial {H}_{i}(T,{P}_{0})}{\partial T}\,{dT}\\ & & -T{\int }_{{T}_{0}}^{T}\displaystyle \frac{\partial {S}_{i}(T,{P}_{0})}{\partial T}\,{dT},\end{array}\end{eqnarray} \tag{ 37 }$$

where we have defined H_i,0 ≡ H_i (T₀, P₀) and S_i,0 ≡ S_i (T₀, P₀) for convenience. We note that

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {H}_{i}(T,{P}_{0})}{\partial T} & = & {C}_{P,i}(T,{P}_{0}),\\ \displaystyle \frac{\partial {S}_{i}(T,{P}_{0})}{\partial T} & = & \displaystyle \frac{{C}_{P,i}(T,{P}_{0})}{T},\end{array}\end{eqnarray} \tag{ 38 }$$

where C_P,i (T, P₀) is the isobaric heat capacity of species i at the reference pressure. It follows that the full expression for the Gibbs free energy is (Powell 1978)

$$\begin{eqnarray}\begin{array}{rcl}{G}_{i} & = & {H}_{i,0}-{{TS}}_{i,0}+{\int }_{{T}_{0}}^{T}{C}_{P,i}(T,{P}_{0})\,{dT}\\ & & -T{\int }_{{T}_{0}}^{T}\displaystyle \frac{{C}_{P,i}(T,{P}_{0})}{T}\,{dT}+{\int }_{{P}_{0}}^{P}{V}_{i}\,{dP}.\end{array}\end{eqnarray} \tag{ 39 }$$

In the evaluation of net chemical reactions, it is not the absolute energy levels that are relevant. Rather, it is the difference in energies between the reactants and the products that inform the equilibrium constants. It is analogous to how potential energies are always defined relative to a reference value. In Heng et al. (2016) and Heng & Lyons (2016), the change in G_i,0 is interpreted as the Gibbs free energy of formation. ⁸ In other words, the reference level of G_i,0 is irrelevant because we are interested only in its difference between the reactants and products of the net chemical reaction. Similarly, since the reference level is irrelevant we are free to interpret H_i,0 as the enthalpy of formation (e.g., Appendix A of Powell 1978), as long as we are cognizant of the fact that we ultimately wish to compute relative, rather than absolute, energies. In the current study, our approach is to use H_i,0 and S_i,0 to compute G_i using Equation (39). We consider this approach to be more general as it is applicable beyond gaseous phases and explicitly allows for nonideal behavior to be computed.

To cast Equation (39) in a more useful form, we need an expression for the heat capacity at constant pressure. Based on experimental data, the study of Holland & Powell (1998) provides empirical fitting functions,

$$\begin{eqnarray}&&{C}_{P,i}\left(T,{P}_{0}\right)={A}_{1,i}+{A}_{2,i}T+{A}_{3,i}{T}^{-2}+{A}_{4,i}{T}^{-1/2},\end{eqnarray} \tag{ 40 }$$

which are valid for T ≤ 2273 K. Table 1 summarizes the fitting coefficients (A_j,i , where j = 1, 2, 3, 4) for the gases considered in the current study, as well as for graphite. For H₂, O₂, CO, CO₂, H₂O, CH₄, H₂S, and S₂, all the thermodynamic data are from the extended data set by Holland & Powell (1998, 2011). For other S- and N-bearing species, that is, SO₂, N₂, NH₃, and HCN, we derive the fitting coefficients by performing regression to the isobaric heat capacity data at 1 bar from the JANAF database; the enthalpy and entropy data are directly taken from the JANAF database.

By substituting Equations (40) into (39), we obtain

$$\begin{eqnarray}\begin{array}{rcl}{G}_{i} & = & {H}_{i,0}-{{TS}}_{i,0}\\ & & +{A}_{1,i}(T-{T}_{0})+\displaystyle \frac{{A}_{2,i}}{2}\left({T}^{2}-{T}_{0}^{2}\right)\\ & & -{A}_{3,i}\left(\displaystyle \frac{1}{T}-\displaystyle \frac{1}{{T}_{0}}\right)+2{A}_{4,i}\left(\sqrt{T}-\sqrt{{T}_{0}}\right)\\ & & -{A}_{1,i}T\mathrm{ln}\left(\displaystyle \frac{T}{{T}_{0}}\right)-{A}_{2,i}T\left(T-{T}_{0}\right)\\ & & +\displaystyle \frac{{A}_{3,i}T}{2}\left(\displaystyle \frac{1}{{T}^{2}}-\displaystyle \frac{1}{{T}_{0}^{2}}\right)+2{A}_{4,i}T\left(\displaystyle \frac{1}{\sqrt{T}}-\displaystyle \frac{1}{\sqrt{{T}_{0}}}\right)\\ & & +{\int }_{{P}_{0}}^{P}{V}_{i}\,{dP}.\end{array}\end{eqnarray} \tag{ 41 }$$

Equation (41) elucidates the various ingredients needed to numerically evaluate G_i : the entropy and enthalpy of formation (first line) and equation of state (last line).

For a pure gas, the conventional practice is to separate out the volume integral term in Equation (41), because this term depends on whether the gas behaves ideally or nonideally. Subsequently, in the computation of the equilibrium constant only G_i,0 enters ΔG_r , whereas the volume integral becomes part of the equilibrium constant. For solids, the ${G}^{{\prime} }$ term in Equation (11) is exactly given by G_i in Equation (41). Thus, it is G_i , rather than G_i,0, that enters ΔG_r .

The entropy (S_i,0) and enthalpy of formation (H_i,0) are taken from Holland & Powell (1998), except for the N- and S-bearing species, for which we obtain these quantities from the JANAF database. These quantities are stated for the reference temperature of T₀ = 298.15 K and the reference pressure of P₀ = 1 bar; they are tabulated in Table 1.

For nonideal gases, we use the compensated-Redlich–Kwong (CORK) equation of state (Holland & Powell 1991) for H₂O, CO₂, CH₄, CO, and H₂. Specifically, this refers to V_i (T, P). For N- and S-bearing species, we assume ideal gases due to the absence of data. For solids, we use the equation of state for graphite provided by Holland & Powell (1998).

In practice, we note that solids have much smaller molar volume and are less compressible than gases, implying that the volume integral is sometimes neglected under lower pressures (e.g., Symonds & Reed 1993, p. 815) and G_i ≈ G_i,0 is obtained.

With the equations of state of nonideal gases in hand, one may obtain the fugacity coefficient of a pure gas by numerically evaluating the volume integral (DeVoe 2020, p. 186),

$$\begin{eqnarray}&&f={P}_{0}{e}^{\tfrac{1}{{ \mathcal R }T}{\int }_{{P}_{0}}^{P}V\,{dP}},\end{eqnarray} \tag{ 42 }$$

where we have set f₀ = P₀ = 1 bar. Since f = ϕ_i P, one obtains ϕ_i of species i numerically.

Holland & Powell (2003) provide fitting functions of activity coefficients (γ_i ) for CO–CO₂–CH₄–H₂O–H₂ mixtures, which are valid for 723 ≤ T/K ≤ 2073 and 500 ≤ P/bar ≤ 4 × 10⁴. In practice, data on activity coefficients are sparser than for fugacity coefficients. Activity data for S-bearing species in the C–H–O–S system are partially available (Evans et al. 2010), but they are unknown when N-bearing species (N₂, HCN, and NH₃) are considered. Since a full treatment of nonideality requires a complete knowledge of every pair-wise interaction among C–H–O–N–S-bearing species, like in the C–H–O system, existing studies choose to either ignore nonideal mixing (Sun & Lee 2022) or only partially consider ideal mixing. For example, Ague et al. (2022) simulated the C–O–H–S subsystem where full nonideality is considered for all the species except for SO₂. We follow this approach in the full C–H–O–N–S system to account for the full nonideality of the C–O–H subsystem (Holland & Powell 2003), whereas for S- and N-bearing species we assume full ideality (γ_i = ϕ_i = 1) due to lack of nonideal thermodynamic parameters. Such an approximation remains to be validated by and thus calls for future experimental measurements of activity coefficient data, especially for N- and S-bearing species in the full C–H–O–N–S system (e.g., Kite et al. 2020).

To build initial intuition for the chemistry of secondary and hybrid atmospheres, we consider the simplest nontrivial chemical system that consists of carbon (C), hydrogen (H), and oxygen (O). We will see that key trends emerge that carry over to systems that include nitrogen and sulfur.

4.1.1. Secondary Atmospheres

Figures 4, 5, and 6 examine the trends in the relative abundances of H₂, H₂O, CO, CO₂, and CH₄ as functions of the atmospheric surface pressure (P), oxygen fugacity of the mantle (${f}_{{{\rm{O}}}_{2}}$), and melt temperature (T), respectively. An atmosphere deriving from a reduced mantle with low carbon content is H₂ and H₂O dominated, not unlike a gas giant exoplanet (top left panel of Figure 4). As the mantle becomes increasingly oxidized (moving from the first to the third column of Figure 4), two important trends appear:

• 1.  
  The atmosphere transitions from being dominated by H₂ to being dominated by H₂O;
• 2.  
  The major carbon carrier transitions from being CO to CO₂.

These trends persist even when the carbon content of the mantle becomes high (bottom row of Figure 4). Methane appears when the carbon content is high, but its abundance is never comparable to those of the other species unless the surface pressure is high (Figure 4(d)). Oxidized mantles further suppress the appearance of methane (Figures 4(d)–(f)). Figure 5 confirms the preceding trends. The competition between CO and CO₂ motivates the use of their relative abundances as a diagnostic for the oxygen fugacity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 demonstrates that the exact melt temperature assumed needs to be considered carefully, because the abundances of CO, CO₂, and CH₄ vary by orders of magnitude across a relatively modest range of T values. An important feature of Figures 4, 5, and 6 is that the carbon-to-oxygen (C/O) ratio of the gaseous phase varies by orders of magnitude. Nonideal effects appear to be minimal, unless large pressures are considered.

4.1.2. Hybrid Atmospheres

The qualitative trends and lessons learned from our study of C–H–O secondary atmospheres carry over when we examine hybrid atmospheres in Figures 7, 8, and 9. A key difference is that, when the partial pressure of atmospheric molecular hydrogen (${P}_{{{\rm{H}}}_{2}}$) becomes large (e.g., Figures 8(b) and (c)), the production of CO and CO₂ is suppressed in favor of H₂O. It is again emphasized that modest changes in the melt temperatures lead to order-of-magnitude changes in the relative abundances of the gaseous species (Figure 9).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The addition of nitrogen (N) and sulfur (S) to the chemical system results in five more gaseous species: hydrogen sulfide (H₂S), sulfur dioxide (SO₂), ammonia (NH₃), molecular nitrogen (N₂), and hydrogen cyanide (HCN). This doubles the number of species compared to C–H–O systems.

For illustration, we choose ${P}_{{{\rm{N}}}_{2}}/P=0.1$ (for secondary atmospheres) and ${P}_{{{\rm{N}}}_{2}}/{P}_{{{\rm{H}}}_{2}}=0.1$ (for hybrid atmospheres), although the chemical trends are also explored for ${P}_{{{\rm{N}}}_{2}}/P$ and ${P}_{{{\rm{N}}}_{2}}/{P}_{{{\rm{H}}}_{2}}$ in the range of 10⁻⁴–1 (Table 3 and Figure B5). While these choices are arbitrary, it is better than choosing an absolute value for ${P}_{{{\rm{N}}}_{2}}$ as we wish to avoid situations where the trends associated with N-bearing species appear artificially weak because an arbitrarily low value of the nitrogen partial pressure was chosen.

Table 3. Explored Parameter Ranges ^{a , b} and Key Findings for the C–H–O–N–S System

Parameter	Range for Secondary Atmosphere	Range for Hybrid Atmosphere
$\mathrm{log}{f}_{{{\rm{S}}}_{2}}$	Between PP−10 and PP	Between PP−10 and PP
${P}_{{{\rm{N}}}_{2}}/P$	10−4–1	not applicable
${P}_{{{\rm{N}}}_{2}}/{P}_{{{\rm{H}}}_{2}}$	Not applicable	10−4–1
Key Findings	Qualitative trends of the C–H–O system carries over; the H2S/SO2 ratio is sensitive to surface pressure and ${f}_{{{\rm{O}}}_{2}};$ the NH2/HCN ratio is sensitive to surface pressure but not to ${f}_{{{\rm{O}}}_{2}};$ the H2S/SO2 ratio is not sensitive to ${f}_{{{\rm{S}}}_{2}}$ but their abundances are; likewise, the NH2/HCN ratio is not sensitive to ${P}_{{{\rm{N}}}_{2}}$ but their abundances are.	Same as the secondary atmospheres

Notes.

^a The range of a_C, $\mathrm{log}{f}_{{{\rm{O}}}_{2}}$, T, and P or ${P}_{{{\rm{H}}}_{2}}$ are the same as for the C–H–O system (Table 2). ^b Some ranges are changed to allow finding physically realistic solutions that are trend revealing, e.g., the $\mathrm{log}{f}_{{{\rm{O}}}_{2}}$ range in Figure 11(d).

Download table as:  ASCIITypeset image

As already noted in Section 2.4, the sulfur fugacity associated with the Earth spans an enormous range of values. Our choices of values for the sulfur fugacity are to allow us to visually display trends in an illustrative manner. Similar to the strategy of parameter value selection in the C–H–O system, we choose $\mathrm{log}{f}_{{{\rm{S}}}_{2}}=$ PP −10 if a value of $\mathrm{log}{f}_{{{\rm{S}}}_{2}}$ needs to be fixed (e.g., Figure 10) and separately explore the chemical trends when $\mathrm{log}{f}_{{{\rm{S}}}_{2}}$ varies between PP −10 and PP (Table 3 and Figure B4).

Zoom In Zoom Out Reset image size

4.2.1. Secondary Atmospheres

Figures 10, 11, 12, 13, and 14 examine the trends in the relative abundances of H₂, H₂O, CO, CO₂, CH₄, H₂S, SO₂, NH₃, N₂, and HCN as functions of the atmospheric surface pressure (P), oxygen fugacity of the mantle (${f}_{{{\rm{O}}}_{2}}$), melt temperature (T), sulfur fugacity of the mantle (${f}_{{{\rm{S}}}_{2}}$), and atmospheric partial pressure of molecular nitrogen (${P}_{{{\rm{N}}}_{2}}$), respectively. The qualitative trends and lessons learned from examining C–H–O chemical systems carry over. Additionally, the following trends emerge:

• 1.  
  The competition between H₂S and SO₂ depends sensitively on the surface pressure (Figure 10) and oxygen fugacity (Figure 11).
• 2.  
  Likewise, the competition between NH₃ and HCN is sensitive to the surface pressure (Figure 10). neither NH₃ nor HCN are oxygen carriers, their abundances are somewhat insensitive to the oxygen fugacity (Figure 11).
• 3.  
  The absolute abundances of H₂S and SO₂ are sensitive to the sulfur fugacity, but the ratio of their abundances is not (Figure 13).
• 4.  
  Likewise, the absolute abundances of NH₃ and HCN are sensitive to the partial pressure of nitrogen, but the ratio of their abundances is not (Figure 14).

As before, varying the melt temperature by a factor ∼2 results in order-of-magnitude variations in the molecular abundances (Figure 12). The competition between H₂S and SO₂ motivates the use of their relative abundances as a supporting diagnostic for the oxygen fugacity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2.2. Hybrid Atmospheres

In our explorations of secondary and hybrid atmospheres, within both the C–H–O and C–H–O–N–S chemical systems, we noticed that methane-dominated atmospheres are somewhat rare. From a theoretical standpoint, it seems difficult to make CH₄-dominated atmospheres. The trends elucidated in Figures 4–14 and B1–B5 suggest that high surface pressures and reduced mantles are needed. Motivated by these trends, we perform a set of Monte Carlo calculations, where we set a_C = 1 (i.e., graphite in presence) and sample the remaining parameters over the following ranges: 10⁻³ ≤ P/bar ≤ 10⁴ (for secondary atmospheres only), ${10}^{-3}\,\leqslant {P}_{{{\rm{H}}}_{2}}/\mathrm{bar}\leqslant {10}^{4}$ (for hybrid atmospheres only), 873 ≤ T/K ≤ 1873, $\mathrm{IW}-6\leqslant \mathrm{log}{f}_{{{\rm{O}}}_{2}}\leqslant \mathrm{IW}+6$, $\mathrm{PP}-10\leqslant \mathrm{log}{f}_{{{\rm{S}}}_{2}}\,\leqslant \mathrm{PP}$, ${10}^{-2}\leqslant {P}_{{{\rm{N}}}_{2}}/P\leqslant 1$ (for secondary atmospheres only), and ${10}^{-2}\leqslant {P}_{{{\rm{N}}}_{2}}/{P}_{{{\rm{H}}}_{2}}\leqslant 1$ (for hybrid atmospheres only). A methane-dominated atmosphere is identified using the criterion

$$\begin{eqnarray}&&{X}_{{\mathrm{CH}}_{4}}\gt \displaystyle \sum _{i}{X}_{i},\end{eqnarray} \tag{ 43 }$$

where the summation is performed over all of the molecular species besides methane.

Figures 15 and 16 demonstrate that the production of methane-dominated atmospheres requires somewhat low melt temperatures (T ≲ 1000 K), reduced mantles (${f}_{{{\rm{O}}}_{2}}\lesssim \mathrm{IW}$), and high atmospheric surface pressures exceeding ∼10 bar. The same calculations performed with C–H–O systems produced essentially the same outcomes (not shown). These rather specific conditions imply that, if a methane-dominated atmosphere is discovered by JWST, some of its atmospheric conditions will be strongly constrained if the methane is indeed sourced by geochemical outgassing (Thompson et al. 2022).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In terms of its scientific intentions, the closest study to compare this study to is Liggins et al. (2020), who used what they termed a "magma degassing code" to study whether early Earth and Mars could have sustained an atmosphere with nonnegligible amounts of molecular hydrogen. As in the present study, Liggins et al. (2020) prescribed or parameterized the oxygen fugacity and total pressure in mixed-phase, chemical-equilibrium, carbon–hydrogen–oxygen–sulfur (C–H–O-S) calculations. Unlike the present study, Liggins et al. (2020) explicitly considered outgassing fluxes to be balanced by atmospheric escape. The former is generally difficult to calculate from first principles, and Liggins et al. (2020) parameterized its value relative to that of modern Earth. Upon specifying the total surface pressure as an input parameter, Liggins et al. (2020) estimated the atmospheric escape efficiency needed to satisfy the corresponding atmospheric abundance of H₂ and the assumed escape flux. In that study, it remains the case that atmospheric escape involves complex radiative transfer processes and chemistry that are difficult to model from first principles and are instead parameterized by a "fudge factor."

Liggins et al. (2020) concluded that ∼1% H₂ abundance (by number) is possible for early Earth, but ∼10% H₂ abundance is unlikely given the plausible space of parameters explored. For early Mars, about 2%–8% of H₂ (by number) is plausible but only if the magmas are water-rich. Liggins et al. (2020) did not extend these explorations to exoplanets and their atmospheric chemical diversity.

In terms of methodology, previous studies closest to the current one are Herbort et al. (2020), Woitke et al. (2021), and Ortenzi et al. (2020). First, some form of volatile budget constraints has to be prescribed in all three studies, namely, various crustal or chondritic bulk compositions assumed in Herbort et al. (2020), a large range of combinations of hydrogen, carbon, oxygen, and nitrogen abundances in Woitke et al. (2021), and initial H₂O and CO₂ budget in Ortenzi et al. (2020). As reasoned earlier and suggested in Appendix A, volatile inventories and/or bulk compositions on exoplanets are subject to large uncertainties, and assimilating these constraints into phase equilibrium modeling of exoplanetary atmospheres introduces more uncertainties; this is why Woitke et al. (2021) conducted an exhaustive exploration of C–O–H–N abundances. On the other hand, our approach here circumvents the necessity of bookkeeping elemental abundances and thus enables transparency and simplicity amenable to analytical insights (e.g., graphite instability). Second, Herbort et al. (2020) surveyed the temperature range 600–3500 K, and one focus of the study is on the water lockup into hydrous minerals whose solid solution, however, remains unaccounted for. In addition to outgassing, Ortenzi et al. (2020) also modeled thermal evolution of rocky interiors, which enables simulation of temporal evolution of outgassing fluxes for rocky planets with stagnant-lid convection, but their approach is parameterized and methane is completely excluded from the model. The lower limit of the temperature range investigated in Woitke et al. (2021) is below 600 K, which, according to the melt temperature range 873–1873 K discussed in Section 2.6, corresponds to nonmagmatic outgassing. Given the disparate model details and implementation, the three studies are comparable to the current one with respect to methodology rather than model results.

As for methane-rich atmospheres, Krissansen-Totton et al. (2018) proposed disequilibrium coexistence of abundant CH₄ and CO₂ as a potential biosignature for exoplanets, which could be further corroborated by absence of abundant CO. Our results can confirm that abundant CH₄ and CO₂ in coexistence represents chemical disequilibrium in the melt temperature range of 873–1873 K. However, as shown by Woitke et al. (2021), at temperatures below 600 K, outgassing, likely metamorphic, can indeed produce abundant coexisting CH₄ and CO₂ in chemical equilibrium, rendering simultaneous detection of CH₄ and CO₂ a false-positive biosignature for exoplanets. Thompson et al. (2022) explored the range of abiotic CH₄ fluxes in the solar system and argues that the short photochemical lifetime of atmospheric CH₄ requires replenishing CH₄ fluxes higher than those from the abiotic outgassing, thus necessitating biological CH₄ fluxes. For an exoplanet, its planetary context and astrophysical environment needs to be carefully considered to gauge the magnitude of abiotic CH₄ fluxes and then the likelihood of CH₄ detection as a true biosignature. Since the current framework is incapable of computing outgassing fluxes, it is an opportunity for future work. We note that whether methane and its coexistence with abundant carbon dioxide may be regarded as a biosignature remains highly context dependent.

The current study has demonstrated that secondary and hybrid atmospheres are expected to exhibit a rich diversity of chemistries. Is it possible to constrain some of these input parameters, such as the oxygen fugacity of the mantle, from the measured spectra of these atmospheres? Generally, the posterior distributions of chemical abundances may be extracted from measured spectra via the technique of Bayesian atmospheric retrieval (see Barstow & Heng 2020 for a recent review). It is expected that abundance ratios may be extracted at a higher precision than absolute abundances. Therefore, we are motivated to explore the theoretical relationships between the various input parameters and the ratio of abundances of simple molecules.

To accomplish this task, we perform Monte Carlo calculations that provide random realizations of ensembles of C–H–O–N–S chemical models. We consider secondary and hybrid atmospheres separately. In the absence of a robust theory for how these atmospheres formed and evolved, we sample each parameter uniformly either in a linear or logarithmic sense. While it has been described elsewhere in the current study, we state the sampled ranges of values of the parameters here for the convenience of the reader (see Table 4 for a summary):

• 1.  
  Melt temperatures (T) are sampled uniformly from 873 to 1873 K.
• 2.  
  For the surface pressures (P) of secondary atmospheres, $\mathrm{log}(P/\mathrm{bar})$ is sampled uniformly from −3 to 4 and $\mathrm{log}\left({P}_{{{\rm{N}}}_{2}}/P\right)$ is sampled uniformly from −2 to 0. For the hydrogen partial pressures (${P}_{{{\rm{H}}}_{2}}$) of hybrid atmospheres, $\mathrm{log}({P}_{{{\rm{H}}}_{2}}/\mathrm{bar})$ is sampled uniformly from −3 to 4 and $\mathrm{log}\left({P}_{{{\rm{N}}}_{2}}/{P}_{{{\rm{H}}}_{2}}\right)$ is sampled uniformly from −2 to 0.
• 3.  
  For the oxygen fugacity (${f}_{{{\rm{O}}}_{2}}$), $\mathrm{log}{f}_{{{\rm{O}}}_{2}}$ is sampled uniformly from IW − 6 to IW + 6. In the absence of better knowledge, the sulfur fugacities are sampled in the same way as $\mathrm{log}{f}_{{{\rm{S}}}_{2}}$ from PP − 10 to PP.
• 4.  
  The carbon activity (a_C), which is a proxy for the carbon content of the mantle, is sampled uniformly as $\mathrm{log}{a}_{{\rm{C}}}$ from −8 to 0.

Table 4. Monte Carlo Sampled Parameter Space for the C–H–O–N–S System and Key Findings

Parameter	Sample Distribution and Space
$\mathrm{log}{a}_{{\rm{C}}}$	uniform a between −8 and 0
$\mathrm{log}{f}_{{{\rm{O}}}_{2}}$	uniform between IW−6 and IW+6
$\mathrm{log}{f}_{{{\rm{S}}}_{2}}$	uniform between PP−10 and PP
T	uniform between 873 and 1873 K
$\mathrm{log}(P/\mathrm{bar})$ for secondary atmospheres	uniform between −3 and 4
$\mathrm{log}({P}_{{{\rm{N}}}_{2}}/P)$ for secondary atmospheres	uniform between −2 and 0
$\mathrm{log}({P}_{{{\rm{H}}}_{2}}/\mathrm{bar})$ for hybrid atmospheres	uniform between −3 and 4
$\mathrm{log}({P}_{{{\rm{N}}}_{2}}/{P}_{{{\rm{H}}}_{2}})$ for hybrid atmospheres	uniform between −2 and 0
Key Findings	Joint ${X}_{{\mathrm{CO}}_{2}}/{X}_{\mathrm{CO}}$ and ${X}_{{{\rm{H}}}_{2}{\rm{S}}}/{X}_{{\mathrm{SO}}_{2}}$ as a strong constraint on the ${f}_{{{\rm{O}}}_{2}}$ of a rocky mantle; complemented by ${X}_{{{\rm{H}}}_{2}{\rm{O}}}/{X}_{{\mathrm{CH}}_{4}}$ and ${X}_{{\mathrm{NH}}_{3}}/{X}_{\mathrm{HCN}}$, melt carbon content and temperature may be further constrained; ${f}_{{{\rm{S}}}_{2}}$ and ${P}_{{{\rm{N}}}_{2}}$ are challenging to constrain from the four ratios.
	Methane-dominated atmospheres requires low melt temperatures (T ≲ 1000 K), reduced mantles (${f}_{{{\rm{O}}}_{2}}\lesssim \mathrm{IW}$), and high atmospheric surface pressures exceeding ∼10 bar.

Note.

^a Set a_C = 1 when investigating methane-dominated atmospheres.

Download table as:  ASCIITypeset image

Figures 17 and 18 shows the outcomes of these Monte Carlo calculations for secondary and hybrid atmospheres, respectively. As motivated by our findings reported earlier in the current study, the abundance ratio of carbon dioxide to carbon monoxide (${X}_{{\mathrm{CO}}_{2}}/{X}_{\mathrm{CO}}$) is chosen as a diagnostic for the oxidation state of the mantle (oxygen fugacity); the abundance ratio of hydrogen sulfide to sulfur dioxide (${X}_{{{\rm{H}}}_{2}{\rm{S}}}/{X}_{{\mathrm{SO}}_{2}}$) is chosen as a supporting diagnostic of the oxidation state. Since water and methane are expected to be readily detectable, their abundance ratio (${X}_{{{\rm{H}}}_{2}{\rm{O}}}/{X}_{{\mathrm{CH}}_{4}}$) is chosen as another observational diagnostic. The abundance ratio of ammonia to hydrogen cyanide (${X}_{{\mathrm{NH}}_{3}}/{X}_{\mathrm{HCN}}$), which are two important nitrogen carriers, completes the set of four observational diagnostics.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As expected, a combination of ${X}_{{\mathrm{CO}}_{2}}/{X}_{\mathrm{CO}}$ and ${X}_{{{\rm{H}}}_{2}{\rm{S}}}/{X}_{{\mathrm{SO}}_{2}}$ potentially sets a strong constraint on the oxygen fugacity of the mantle of a rocky exoplanet. Similarly, the carbon content (carbon activity) and melt temperature may be constrained using a combination of all four abundance ratios. For secondary atmospheres, it is perhaps unsurprising that the atmospheric surface pressure cannot be constrained from measuring these abundance ratios. However, the hydrogen content of hybrid atmospheres appears to be sensitive to these four abundance ratios. For both types of atmospheres, it appears that the sulfur fugacity and nitrogen content are somewhat challenging to constrain from measuring these four abundance ratios.

It is worth noting that, in cycle 1 of the JWST Guest Observer program ⁹ alone, there are approved proposals targeting about 10 small exoplanets. Therefore, it is conceivable that the mantle oxygen fugacities of some small exoplanets may be constrained in the near future.

In the study of gas giant exoplanets, astronomers often use the C/O ratio and the "metallicity" as the control parameters. The study of the C/O ratio is motivated by its varying value with distance from the star within a protoplanetary disk, due to the differing condensation temperatures of simple molecules (so-called ice lines or snow lines; Öberg et al. 2011). By measuring the C/O ratio of a gas giant exoplanet, the hope is that one may then derive its original site of formation within a protoplanetary disk (see Öberg et al. 2011 for caveats). The "metallicity" assumes that the entire set of elemental abundances, for both refractory and volatile elements, have ratios that are locked to their solar values. Only in this way may a set of elemental abundances be described by a single number. Motivated by an empirical trend in the solar system, previous studies have claimed an empirical relationship between the metallicities and the masses of exoplanets (e.g., Wakeford et al. 2017), but this assumes that the derived abundances of water translate into oxygen abundances that scale with the abundances of refractory elements in a straightforward way (Heng 2018).

Since secondary and hybrid atmospheres are at least partially sourced by geochemical outgassing, their chemistries do not straightforwardly trace their formation histories unlike in the case of gas giant exoplanets. This implies that the C/O ratio has less relevance for rocky exoplanets. When the atmosphere is sourced by outgassing, the volatile and refractory elements are partitioned between their gaseous forms (atmosphere), the melt (which the gases may dissolve into), and rocks (which are composed of a mixture of minerals). This implies that it will be complicated, if not impossible (from the viewpoint of interpreting astronomical data), to decipher the relative abundances of volatile and refractory elements, thus rendering the simplistic concept of metallicity (as defined in the preceding paragraph) suspect.

Therefore, we suggest that a better way of thinking about rocky exoplanets is in terms of the oxygen fugacity and carbon content of their mantles. In the solar system, an empirical trend exists such that more massive bodies tend to have more oxidized mantles (e.g., Wadhwa 2008; Cartier & Wood 2019). Our current work suggests that there is a clear path toward deriving the oxygen fugacities (${f}_{{{\rm{O}}}_{2}}$) of the rocky mantles of exoplanets by inferring accurate abundances of CO₂ and CO via atmospheric retrieval performed on high-quality spectra (by, e.g., the JWST), but future work needs to correct for the effects of photochemistry. An exciting prospect for exoplanet science will be to quantify the relationship between ${f}_{{{\rm{O}}}_{2}}$ and the masses of exoplanets for a sample of super Earths and sub-Neptunes.

The current study is the first to consider secondary and hybrid atmospheres within the same theoretical framework. But it will certainly not be the last, as there are many other processes to explore and other members of the model hierarchy of geochemical outgassing to construct and study. The models constructed in the current study are fundamentally 0D, meaning that the temperature–pressure profile of the mantle and atmosphere are not considered. These thermal gradients could lead to chemical abundance gradients and processes such as cold traps. Extending these calculations beyond 0D requires coupling the chemical models to radiative transfer calculations that consider relevant dynamical processes such as convection and large-scale atmospheric circulation. Another process we have neglected in the current study is the general solubility of gases in melts.

The full or ideal model would quantify the joint evolution of the interior, mantle, and atmosphere of the exoplanet, as well as the radiative influence of its host star. Its chemical geodynamics would be elucidated such that one could calculate the oxygen fugacity and carbon content of the mantle from first principles, rather than parameterize them by single numbers. In the current study, we have parameterized the carbon content via the carbon activity, from which it is difficult to directly compute the carbon content by mass. The atmospheric surface pressure and H₂ partial pressure should be outcomes of a calculation balancing the outgassing flux and atmospheric escape (e.g., Liggins et al. 2020). The latter process involves performing radiative transfer calculations alongside ionic chemistry in parts of the atmosphere where the fluid approximation no longer holds. In essence, we have strived for simplicity over sophistication in the current study by choosing to parameterize these formidable processes using two numbers (${f}_{{{\rm{O}}}_{2}}$ and P or ${P}_{{{\rm{H}}}_{2}}$).

In order to find analytical solutions that help build our intuitive understanding, we consider the simplest possible system comprising only two elements—hydrogen (H) and oxygen (O)—and only three molecular species—H₂, O₂, and H₂O. The linear algebra analysis method in Powell et al. (1998) tells us that there exists only one independent chemical reaction for this subsystem:

$$\begin{eqnarray}&&{{\rm{H}}}_{2}+0.5{{\rm{O}}}_{2}\leftrightarrows {{\rm{H}}}_{2}{\rm{O}}.\end{eqnarray} \tag{ A1 }$$

Aiming for analytical solutions, we first ignore nonideality and thus all fugacity quantities become equivalent to pressure quantities, that is, f_i = P_i . Writing out the equilibrium constant for reaction (A1),

$$\begin{eqnarray}&&K=\displaystyle \frac{{P}_{{{\rm{H}}}_{2}{\rm{O}}}}{{P}_{{{\rm{H}}}_{2}}{f}_{{{\rm{O}}}_{2}}^{\tfrac{1}{2}}},\end{eqnarray} \tag{ A2 }$$

where the equilibrium constant can be calculated in the way it was in the main text:

$$\begin{eqnarray}&&K=\exp \left(\displaystyle \frac{-{{\rm{\Delta }}}_{r,A1}G}{{RT}}\right).\end{eqnarray} \tag{ A3 }$$

One nice property of this simple H–O system is that, in Equation (A2), ${P}_{{{\rm{H}}}_{2}{\rm{O}}}$ and ${P}_{{{\rm{H}}}_{2}}$ have the same exponent of 1, which makes it convenient to convert Equation (A2) to an expression of volume mixing ratios:

$$\begin{eqnarray}&&K=\displaystyle \frac{{P}_{{{\rm{H}}}_{2}{\rm{O}}}}{{P}_{{{\rm{H}}}_{2}}{f}_{{{\rm{O}}}_{2}}^{\tfrac{1}{2}}}=\displaystyle \frac{\left({P}_{{{\rm{H}}}_{2}{\rm{O}}}/P\right)}{\left({P}_{{{\rm{H}}}_{2}}/P\right){f}_{{{\rm{O}}}_{2}}^{\tfrac{1}{2}}}=\displaystyle \frac{{X}_{{{\rm{H}}}_{2}{\rm{O}}}}{{X}_{{{\rm{H}}}_{2}}{f}_{{{\rm{O}}}_{2}}^{\tfrac{1}{2}}},\end{eqnarray} \tag{ A4 }$$

where P is the total gas pressure of the subsystem.

Another nice property of reaction (A1) is that it involves only gaseous species and thus Δ_r,A1 G and K are only temperature dependent, that is, they do not depend on total pressure P (see the discussion in Section 2.3). We can thus write K = K(T). Rearranging Equation (A4) and considering unity sum lead to

$$\begin{eqnarray}\left\{\begin{array}{ll}\displaystyle \frac{{X}_{{{\rm{H}}}_{2}{\rm{O}}}}{{X}_{{{\rm{H}}}_{2}}}=K(T){f}_{{{\rm{O}}}_{2}}^{\tfrac{1}{2}}, & \\ {X}_{{{\rm{H}}}_{2}{\rm{O}}}+{X}_{{{\rm{H}}}_{2}}=1, & \mathrm{unity}\ \mathrm{sum}.\end{array}\right.\end{eqnarray} \tag{ A5 }$$

Equation (A5) suggests that as long as temperature and oxygen fugacity (T and ${f}_{{{\rm{O}}}_{2}}$) are given, mixing ratios of H₂ and H₂O are fully determined for the simple H–O system. In other words, the solution of mixing ratios is pressure independent; the simplicity of this subsystem allows us to peel off all pressure variables in the system and thus reduce its degree of freedom by one. In spite of this, to compute partial pressures ${P}_{{{\rm{H}}}_{2}{\rm{O}}}$ and ${P}_{{{\rm{H}}}_{2}}$, we still need a total pressure P such that ${P}_{{{\rm{H}}}_{2}{\rm{O}}}={X}_{{{\rm{H}}}_{2}{\rm{O}}}P$ and ${P}_{{{\rm{H}}}_{2}}={X}_{{{\rm{H}}}_{2}}P$. This total pressure could be derived from volatile budget constraints (see below).

Now we start considering gas solubilities in the H–O system. At the outset, when accounting for volatile dissolution into molten rocks, solubility laws are generally needed that dictate how volatiles are partitioned between the atmosphere reservoir and the molten silicate reservoir. Naturally, total volatile budgets that are to be thermodynamically partitioned between the two reservoirs also need to be prescribed. With these two additional constraints, one can solve for how much mass of each volatile is distributed among the two reservoirs. For the atmosphere reservoir, if the surface area (A) and gravity (g) of a planet are further known or assumed, one can convert volatile mass M_i,atm in the atmosphere to its partial pressure P_i by P_i = M_i,atm g/A (e.g., Ortenzi et al. 2020; Bower et al. 2022).

A.2.1. Hydrogen Budget

A.2.2. Solubility Laws

With silicate melt as an infinite oxygen reservoir and the negligible dissolution of molecular hydrogen (Bower et al. 2022), only the dissolution of water into melt needs to be considered through its solubility law:

$$\begin{eqnarray}&&{X}_{{{\rm{H}}}_{2}{\rm{O}}}=\alpha {P}_{{{\rm{H}}}_{2}{\rm{O}}}^{0.5},\end{eqnarray} \tag{ A6 }$$

where ${X}_{{{\rm{H}}}_{2}{\rm{O}}}$ is the mass fraction of H₂O in the melts and α = 534 ppmw/bar^0.5 is experimentally determined for peridotitic melts (Sossi et al. 2023).

A.2.3. Mass Balance

Accompanying the two constraints from hydrogen budget and water solubility is one more equation of hydrogen mass balance:

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{{{\rm{H}}}_{2}{\rm{O}},\ \mathrm{melt}}}{{m}_{{{\rm{H}}}_{2}{\rm{O}}}}+\displaystyle \frac{{M}_{{{\rm{H}}}_{2}{\rm{O}},\ \mathrm{atm}}}{{m}_{{{\rm{H}}}_{2}{\rm{O}}}}+\displaystyle \frac{{M}_{{{\rm{H}}}_{2},\ \mathrm{atm}}}{{m}_{{{\rm{H}}}_{2}}}={{ \mathcal M }}_{{{\rm{H}}}_{2},\mathrm{total}\ \mathrm{in}\ \mathrm{mol}},\end{eqnarray} \tag{ A7 }$$

where m_i is the molecular weight of species i and ${{ \mathcal M }}_{{{\rm{H}}}_{2},\mathrm{total}\ \mathrm{in}\ \mathrm{mol}}$ is the prescribed hydrogen budget. The first term is the moles of H₂ contained in H₂O dissolved in the melt, the second is the amount of H₂ contained in H₂O in the atmosphere, and the third term is the H₂ amount in the atmosphere; the H₂ dissolved in the melt is ignored due to its low solubility.

First, according to the simple relation P_i = M_i,atm g/A, we can replace M_i,atm with P_i -bearing terms:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{M}_{{{\rm{H}}}_{2}{\rm{O}},\ \mathrm{melt}}}{{m}_{{{\rm{H}}}_{2}{\rm{O}}}}+{10}^{5}\times \displaystyle \frac{{P}_{{{\rm{H}}}_{2}{\rm{O}}}A}{{{gm}}_{{{\rm{H}}}_{2}{\rm{O}}}}+{10}^{5}\times \displaystyle \frac{{P}_{{{\rm{H}}}_{2}}A}{{{gm}}_{{{\rm{H}}}_{2}}}\\ \,\,\,={{ \mathcal M }}_{{{\rm{H}}}_{2},\mathrm{total}\ \mathrm{in}\ \mathrm{mol}},\end{array}\end{eqnarray} \tag{ A8 }$$

where the prefactor 10⁵ arises from unit conversion of partial pressures from bar to Pa.

Second, the solubility law for H₂O can be substituted:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{M}_{m}\alpha {P}_{{{\rm{H}}}_{2}{\rm{O}}}^{0.5}}{{m}_{{{\rm{H}}}_{2}{\rm{O}}}}+{10}^{5}\times \displaystyle \frac{{P}_{{{\rm{H}}}_{2}{\rm{O}}}A}{{{gm}}_{{{\rm{H}}}_{2}{\rm{O}}}}+{10}^{5}\times \displaystyle \frac{{P}_{{{\rm{H}}}_{2}}A}{{{gm}}_{{{\rm{H}}}_{2}}}\\ \,\,\,={{ \mathcal M }}_{{{\rm{H}}}_{2},\mathrm{total}\ \mathrm{in}\ \mathrm{mol}},\end{array}\end{eqnarray} \tag{ A9 }$$

where M_m is the total mass of melt (molten rocks).

Third, chemical equilibrium relation ${P}_{{{\rm{H}}}_{2}{\rm{O}}}=K(T){P}_{{{\rm{H}}}_{2}}{f}_{{{\rm{O}}}_{2}}^{\tfrac{1}{2}}$ from Equation (A2) can be inserted

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{M}_{m}\alpha {K}^{\tfrac{1}{2}}{f}_{{{\rm{O}}}_{2}}^{\tfrac{1}{4}}}{{m}_{{{\rm{H}}}_{2}{\rm{O}}}}{P}_{{{\rm{H}}}_{2}}^{\tfrac{1}{2}}+{10}^{5}\times \displaystyle \frac{{{Kf}}_{{{\rm{O}}}_{2}}^{\tfrac{1}{2}}A}{{{gm}}_{{{\rm{H}}}_{2}{\rm{O}}}}{P}_{{{\rm{H}}}_{2}}+{10}^{5}\\ \,\times \displaystyle \frac{A}{{{gm}}_{{{\rm{H}}}_{2}}}{P}_{{{\rm{H}}}_{2}}={{ \mathcal M }}_{{{\rm{H}}}_{2},\mathrm{total}\ \mathrm{in}\ \mathrm{mol}},\end{array}\end{eqnarray} \tag{ A10 }$$

which, if letting $x={P}_{{{\rm{H}}}_{2}}^{\tfrac{1}{2}}$, becomes a quadratic equation in x.

A.2.4. Solution Method

The preceding derivation suggests that Equation (A10) accounts altogether for chemical speciation (A1), gas dissolution (A6), and total hydrogen budget of the simple H–O system. With a planet radius R_p and density ρ, the surface area $A=4\pi {R}_{p}^{2}$, planet mass ${M}_{p}=\rho (4/3)\pi {R}_{p}^{3}$, and surface gravity $g={{GM}}_{p}/{R}_{p}^{2}$ can be obtained. The total melt mass is M_m = X_m M_p if the mass fraction of molten rocks is denoted by X_m . For simplicity and as proof of concept, we take Earth values for R_p and ρ and assume a whole-mantle magma ocean leading to X_m ≈ 0.67. Once a total H₂ budget is further prescribed, ${P}_{{{\rm{H}}}_{2}}$ is solvable from Equation (A10) at given temperature and oxygen fugacity values. A representative temperature T = 1600 K, a hydrogen budget of one Earth ocean (7.8 × 10²² mol H₂), and a range of oxygen fugacity values are selected to observe model trend (Figure A1). With solved ${P}_{{{\rm{H}}}_{2}}$, water partial pressure ${P}_{{{\rm{H}}}_{2}{\rm{O}}}$ and total pressure P are further computed through ${P}_{{{\rm{H}}}_{2}{\rm{O}}}=K(T){P}_{{{\rm{H}}}_{2}}{f}_{{{\rm{O}}}_{2}}^{\tfrac{1}{2}}$ and $P\,={P}_{{{\rm{H}}}_{2}}+{P}_{{{\rm{H}}}_{2}{\rm{O}}}$, respectively. To calculate a case without gas dissolution, simply set the solubility coefficient α to zero in equation (A10). Note that in dissolution-free cases, the constraint from H₂ budget is fundamentally only to determine various pressures (${P}_{{{\rm{H}}}_{2}{\rm{O}}}$, ${P}_{{{\rm{H}}}_{2}}$, and P) because volume mixing ratios (${X}_{{{\rm{H}}}_{2}{\rm{O}}}$ and ${X}_{{{\rm{H}}}_{2}}$) are already determined by temperature and oxygen fugacity only (see Section A.1).

Zoom In Zoom Out Reset image size

Figure A1 shows the results with (solid curves) and without (dotted curves) considering gas dissolution into melts. Under low oxygen fugacities (∼IW−6), both the dissolution and dissolution-free cases feature almost pure H₂ atmosphere (red curve in Figure A1(b)). Moreover, due to the negligible H₂ solubility, the entire hydrogen budget (one Earth ocean) resides exclusively in the atmosphere reservoir in the form of H₂ gas. This enables us to estimate the corresponding ${P}_{{{\rm{H}}}_{2}}$ under low oxygen fugacities, that is, ${P}_{{{\rm{H}}}_{2}}={M}_{{{\rm{H}}}_{2},\mathrm{atm}}g/A\approx 30\,\mathrm{bar}$, where ${M}_{{{\rm{H}}}_{2},\mathrm{atm}}$ is equivalent to the prescribed hydrogen budget in mass unit. The estimated ${P}_{{{\rm{H}}}_{2}}\approx 30\,\mathrm{bar}$ is verified by the numerical solutions for both the dissolution and dissolution-free cases in Figure A1(a) (red curves). In contrast, under high oxygen fugacity conditions (∼IW+6), the dissolution and dissolution-free cases both feature almost pure H₂O atmosphere (blue curve in Figure A1(b)). In this case, ignoring H₂O dissolution means nearly all the hydrogen budget is oxidized to H₂O and resides exclusively in the atmosphere. This enables estimating ${P}_{{{\rm{H}}}_{2}{\rm{O}}}={M}_{{{\rm{H}}}_{2}{\rm{O}},\mathrm{atm}}g/A\approx 30\times 9=270\,\mathrm{bar}$, where the multiplier 9 is the ratio of molecular weight between H₂O and H₂ and stems from the total hydrogen budget being completely oxidized to H₂O. Such an estimate of ${P}_{{{\rm{H}}}_{2}{\rm{O}}}\approx 270\,\mathrm{bar}$ at high ${f}_{{{\rm{O}}}_{2}}$ is verified by the numerical solution of the dissolution-free case in Figure A1(a) (dotted blue curve). When H₂O dissolution is enabled at high ${f}_{{{\rm{O}}}_{2}}$, its sequestration into coexisting melts greatly reduces its atmospheric partial pressure, as suggested by Figure A1(a) (compare solid and dotted blue curves). Furthermore, under oxidizing conditions, dissolution-induced drop of ${P}_{{{\rm{H}}}_{2}{\rm{O}}}$ drives ${P}_{{{\rm{H}}}_{2}}$ to even lower levels, as revealed by comparison of solid and dotted red curves in Figure A1(a). Overall, because the H₂O converted from H₂ oxidation tends to dissolve in melts and thus decreases total pressure, the trend of total pressure variation with oxygen fugacity is flipped from dissolution-free to dissolution cases (dotted and solid magenta curves in Figure A1(a)).

Despite the considerable effect of gas solubility on atmosphere pressures, it is worth emphasizing that for the H–O system investigated here, the variational trend of atmospheric volume mixing ratios with oxygen fugacity remains immune to the effect of gas (H₂O in this study) solubility, as attested to by the fully overlapped curves in Figure A1(b). This is unsurprising because it is analytically elucidated in Section A.1 that this trend depends on temperature only; gas solubility and the associated hydrogen budget play the mere role of further constraining partial and total pressure quantities. This immunity of H₂O and H₂ mixing ratios to gas solubilities and total hydrogen budget can also be robustly confirmed by our dissolution-free results and previous studies. For example, the analysis in Figure A1(b) shows that the crossover of ${X}_{{{\rm{H}}}_{2}{\rm{O}}}$ and ${X}_{{{\rm{H}}}_{2}}$ occurs at ${f}_{{{\rm{O}}}_{2}}\sim$ IW. This crossover at ∼IW is first immune to expanding the system from H–O to C–H–O and C–H–O–N–S, as confirmed by closer inspection of Figures 5, 8, 11, and B2. Second, as shown by Figure 1 in Ortenzi et al. (2020) and Figures 2 and 8 in Sossi et al. (2023), considering or not considering gas solubilities does not alter the crossover position either.

It is worth clarifying that, despite the above invariability for the H–O system, mixing ratios of other gas species, for example, CO and CO₂, are indeed affected by dissolution behaviors (Ortenzi et al. 2020). Nevertheless, we reiterate that the invariability would not be discovered had we not followed the hierarchical modeling approach to begin with dissolution-free models.