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

The atmospheres of small exoplanets likely derive from a combination of geochemical outgassing and primordial gases left over from formation. Secondary atmospheres, such as those of Earth, Mars and Venus, are sourced by outgassing. Persistent outgassing into long-lived, primordial, hydrogen-helium envelopes produces hybrid atmospheres of which there are no examples in the Solar System. We construct a unified theoretical framework for calculating the outgassing chemistry of both secondary and hybrid atmospheres, where the input parameters are the surface pressure, oxidation and sulfidation states of the mantle, as well as the primordial atmospheric hydrogen, helium and nitrogen content. Non-ideal gases (quantified by the fugacity coefficient) and non-ideal mixing of gaseous components (quantified by the activity coefficient) are considered. Both secondary and hybrid atmospheres exhibit a rich diversity of chemistries, including hydrogen-dominated atmospheres. The abundance ratio of carbon dioxide to carbon monoxide serves as a powerful diagnostic for the oxygen fugacity of the mantle, which may conceivably be constrained by James Webb Space Telescope spectra in the near future. Methane-dominated atmospheres are difficult to produce and require specific conditions: atmospheric surface pressures exceeding $\sim 10$ bar, a reduced (poorly oxidised) mantle and diminished magma temperatures (compared to modern Earth). Future work should include photochemistry in these calculations and clarify the general role of atmospheric escape. Exoplanet science should quantify the relationship between the mass and oxygen fugacity for a sample of super Earths and sub-Neptunes; such an empirical relationship already exists for Solar System bodies.

1. INTRODUCTION 1.1.Motivation I: unified framework for outgassing theory 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, e.g., 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.
Calculations of outgassing chemistry is a mature topic in the geochemical literature (e.g., French 1966;Ohmoto & Kerrick 1977;Holloway 1981Holloway , 1987;;Connolly & Cesare 1993;Huizenga 2011).However, studies tend to focus on Earth-centric 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.

Motivation II: puzzle of super Earths and sub-Neptunes
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.6R ⊕ (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 photo-evaporation or core-powered mass loss (Rogers et al. 2021).Furthermore, Kite et al. (2019) proposed that considerable H 2 dissolution into magma oceans stalls sub-Neptune growth via H 2 accretion, leading to a "radius cliff" (Figure 1 in their study) and thus an over-abundance of sub-Neptunes in the population statistics.
The atmospheric chemistry of super Earths and sub-Neptunes remain 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)?Mod-els 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 arrive, it is imperative that a hierarchy of theoretical models for predicting outgassing chemistry is available.
Figure 1.Schematic depicting secondary versus hybrid atmospheres.Secondary atmospheres are fully sourced by geochemical outgassing, while hybrid atmospheres derive from outgassing into a primordial hydrogen-helium envelope left over from the process of formation and evolution.

Motivation III: non-ideal effects
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 behaviours 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 ability, we include these non-ideal 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 non-existent for arbitrary mixtures.

Structure of paper
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 carbonhydrogen-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 methanedominated 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 generalisation of the pressure accounting for non-ideal-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., page 78 of DeVoe 2020).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., page 114 of Denbigh 1981).
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), where G is the specific Gibbs free energy (Gibbs free energy per unit mass), G 0 ≡ G(P 0 ), P 0 is the reference pressure (usually set to 1 atm or 1 bar), R is the specific gas constant and T is the temperature.The task is to generalise Ψ = P/P 0 for gas mixtures with non-ideal-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.

Fugacity
For a pure gas (one species), 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, ρ = P/RT and one obtains equation ( 1) and Ψ = P/P 0 .For non-ideal 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), 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.

Activity
At a pressure P and temperature T , the Gibbs free energy of a gas component i in a gaseous mixture is (page 115 of Denbigh 1981), where G is given by equation ( 2) and X i is the volume mixing ratio (relative abundance by number) of the i-th gaseous species.One now defines the activity as the generalisation of where γ i is the activity coefficient of the i-th 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 i-th gaseous species in a gas mixture with non-ideal mixing of its constituents is (page 287 of Denbigh 1981)

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 (page 352 of DeVoe 2020) For gases, one may write We have written the fugacity of the i-th gaseous species as f i = a i f .Since P 0 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 0 = ϕ i,0 P 0 = P 0 .Given that pure-gas fugacity f is a generalisation of total pressure P and that activity a i is a generalisation of volume mixing ratio X i , the fugacity of the i-th gaseous species f i = a i f can be deemed a generalisation 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 0 is sometimes1 referred to as the activity (see also page 287 of Denbigh 1981).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, where the quantity G ′ depends on the pressure P (and is not referenced to P 0 , 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 ′ .The exact expression for G ′ is provided in Section 3. In the current study, we consider only the activity associated with graphite.

Constructing the equilibrium constant from Ψ in gaseous C-H-O system
Consider the standard net chemical reaction for converting methane to carbon monoxide (e.g., Burrows & Sharp 1999;Lodders & Fegley 2002;Heng & Lyons 2016), In the limit of an ideal gas (ϕ i = 1) with ideal mixing (γ i = 1), we have Ψ = P i /P 0 .We may write down the Gibbs free energies associated with the products and reactants, The differences in Gibbs free energies are It follows that Following the reasoning by Heng et al. (2016), the system attains chemical equilibrium by adjusting to ∆G = 0. ∆G 0 is interpreted as the change of Gibbs free energy of formation at the reference pressure P 0 .By comparing the preceding equation with equation ( 8), the equilibrium constant is in agreement with equations ( 6) and ( 26) of Heng & Lyons (2016).The Gibbs free energy of the reaction is ∆G r = ∆G 0 .Generally, the approach of constructing K eq from Ψ is consistent with equation (13) of Heng et al. (2016).
2.3.Constructing equilibrium constants in mixed-phase C-H-O system of French (1966) French (1966) considered the system of net chemical reactions, 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 The critical detail is that G ′ C (P, T ) depends on the pressure P , whereas G O2,0 (P 0 ) and G CO2,0 (P 0 ) depend only on the reference pressure P 0 .
The differences in Gibbs free energies is 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, where the Gibbs free energy of the reaction is In addition to the usual Gibbs free energies of formation associated with the gaseous species, which are defined at the reference pressure P 0 , one considers the Gibbs free energies of formation associated with graphite, which is 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 ( 17), the equilibrium constants are If one again sets a C = 1 and omits writing the factors of P 0 explicitly (by setting them to unity), then one recovers equations (2) to (4) of French (1966).

Review: oxygen and sulfur fugacities
The oxygen fugacity (f O2 ) quantifies how reduced or oxidised 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) Fe + 0.5O 2 ⇆ FeO. (23) If factors of P 0 are ignored, then the oxygen fugacity buffered by this reaction has the following equilibrium constant, By using equation ( 8) and performing a change of base of the logarithm (from base 2 to 10), we obtain An idealisation 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), where we follow the convention in the geochemical literature of abbreviating log f IW O2 as "IW".The oxygen fugacity may span 50 orders of magnitude in value (Frost 1991) and it is both inconvenient and non-intuitive to quote its absolute value (in bar).Rather, geochemists favor quoting the oxygen fugacity relative to simple2 buffers such as IW.From a computational standpoint, IW is a convenient reference point that allows log f O2 to be reported in terms of intuitive, order-ofunity 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 log (a FeO /a Fe ) = −1.Equation (25) then yields log f O2 = 2∆G r /RT ln 10 − 2. Using equation (26), we obtain log f O2 = 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 2 SiO 4 ), magnetite (Fe 3 O 4 ) and quartz (SiO 2 ), 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), 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 2 ) (Ferry & Baumgartner 1987), which we abbreviate as "PP" for convenience.Froese & Gunter (1976) provides an empirical fitting function based on the experimental data from Toulmin & Barton (1964), 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.
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 S2 and f O2 as independent input parameters.
2.5.Generalised model for mixed-phase C-H-O-N-S system 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), (30) Similar to the oxygen fugacity, we treat the sulfur fugacity f S2 as an input parameter (e.g., Holloway 1981).The nitrogen content of the atmosphere is parametrised via the partial pressure of nitrogen P N2 .Assuming these two additional input parameters allows the number of unknowns and equations to be equal.

Equilibrium constants
Modern thermodynamic data allow us to consider nonideal mixing (γ i ̸ = 1) and non-ideal gas behavior (ϕ i ̸ = 1), unlike for French (1966).For compactness of notation, we define This motivates us to write down more generalised expressions for the equilibrium constants, 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 non-ideal mixing (expressed through γ i ) from non-ideal-gas equations of state (expressed through ϕ i ) by explicitly writing f i = ϕ i P i .

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, which depend on the parameters f O2 , f S2 and P N2 .We note that P HCN ∝ P 1/2 H2 .The total pressure is given by P =P H2 + P He + P O2 + P S2 + P N2 + P CO2 + P CO + P CH4 + P H2O + P NH3 + P HCN + P SO2 + P H2S . (34)

Solution method
The solution method depends on whether one is solving for a secondary or hybrid atmosphere.Following French (1966), if one considers a secondary atmosphere then the partial pressure of molecular hydrogen (P H2 ) is a quantity that one solves for.In the original C-H-O system considered by French (1966), one solves a quadratic equation for P H2 .When we include the additional molecular species beyond what French (1966) considered, we need to numerically solve equation ( 34), which provides a quartic equation in x where x ≡ P 1/2 H2 .The input parameters are T , P , f O2 , f S2 and P N2 .
For a hybrid atmosphere, one imagines a primordial envelope of molecular hydrogen that is the consequence of the formation and evolutionary history of the exoplanet.In the absence of such a robust theory of formation and evolution, we prescribe the value of P H2 but solve for P .The input parameters are T , P H2 , f O2 , f S2 and P N2 .This is done by iteration, where one makes a first guess for P and updates the contributions to it and the equilibrium constants (which depend on P ) until convergence is attained.Physically, one is solving for a system in which one has a mantle with fixed oxygen and sulfur fugacities interacting with a primordial envelope of hydrogen and helium.
Our implementation is as follows: • Choose values for T and P (secondary atmospheres) or P H2 (hybrid atmospheres).
• For all species other than O 2 and S 2 , calculate γ i based on guessed or iterated atmospheric compositions.Since data are unavailable for N-and S-bearing molecules, we set ϕ i = γ i = 1 for these species (see discussion in Section 3.6).Iterate between γ i and atmosphere composition (partial pressures P i ) until convergence attains.
The partial pressure P He is an additional input parameter that considers the presence of helium as an inert gas contributing to the atmospheric pressure.When solving for a secondary atmosphere, P is replaced by P − P He , which alters the solution for P H2 .When solving for a hybrid atmosphere, the assumed relative content of hydrogen versus helium affects the partial pressures of various gases directly.The inclusion of helium is motivated by the work of Hu et al. (2015), who postulated that atmospheric escape over ∼ 0.1 Gyr timescales may evolve primordial H 2 -He atmospheres to He-dominated ones.However, as the parameter space is already very broad, its influence and implications are deferred to a future study and we set P He = 0 throughout this study.
The partial pressure of molecular nitrogen P N2 plays a similar role with the key difference that it participates in the net chemical reactions, while helium does not.
It is clear from the above solution procedure that our outgassing model is zero-dimensional, meaning that it does not resolve the respective temperature and pressure gradients in the atmosphere and the melt-producing interior, and hence is a geochemical "box" model.Moreover, it's noteworthy that for the "box" of melt-bearing interior, gas dissolution into the melts is neglected, which differs from previous studies (Gaillard & Scaillet 2014;Ortenzi et al. 2020;Bower et al. 2022;Gaillard et al. 2022).We intentionally avoid considering gas solubilities in this study for two reasons.First, consideration of gas solubility requires prescription of the bulk volatile budget of an exoplanet (in absolute mass terms), which is subject to large uncertainties.This is also the reason why Ortenzi et al. (2020), Bower et al. (2022) and Gaillard et al. (2022) explored a range of initial volatile inventories.In other words, including gas solubility into the calculations in turn introduces more free parameters that are difficult to assign values to.Second, the method ignoring gas solubility and thus without volatile budget constraints corresponds to global Gibbs energy minimization, whereas the method that includes volatile budgets corresponds to "constrained" Gibbs energy minimization; the difference is that between phase diagram and pseudosection as elaborated in Powell et al. (1998).Within the hierarchical modeling approach (Held 2005), our model is one step simpler than models that include gas dissolution into melts.Nonetheless, in Appendix A, we provide an example using the hydrogen-oxygen (H-O) subsystem to illustrate how gas dissolution into melt can be incorporated into the current modeling framework.

Relevant range of temperatures and pressures for outgassing
The input pressure P is interpreted as the surface pressure of the atmosphere (Gaillard & Scaillet 2014).We consider P = 1 mbar to 10 4 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 2 P/g.It follows that the ratio of the atmospheric mass to the total mass of the exoplanet is (35) It is thus not unreasonable to regard sub-Neptunes as exoplanets with surface pressures ∼ 10 4 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 (Hirchmann 2012; Keppler & Golabek 2019;Sossi et al. 2020;Gaillard et al. 2022;Bower 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 4 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 7a 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 4 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 non-negligible 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.The study of Gaillard & Scaillet (2014) chose T = 1573 K, which we adopt and approximate as 1600 K to emphasise 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.Note: Hi,0 and Si,0 are stated for T0 = 298.15K and P0 = 1 bar.For SO2, N2, NH3 and HCN, Hi,0 and Si,0 data are directly taken from the JANAF database, and Aj,i (j = 1, 2, 3, 4) are from regression of isobaric (at 1 bar) heat capacity data from the JANAF database.All other data from Holland & Powell (1998).
For both gaseous and solid phases, the Gibbs free energy is given by equation (A1) of Powell (1978), where G i = G i (T, P ) and G i,0 ≡ G i (T, P 0 ).For convenience, we write the integrand as V i .The Gibbs free energy at the reference pressure P 0 may be further expanded as G i,0 = H i (T, P 0 )−T S i (T, P 0 ), where H i (T, P 0 ) and S i (T, P 0 ) are the enthalpy and entropy at the reference pressure, respectively.If we further expand H i (T, P 0 ) and S i (T, P 0 ) about a reference temperature T 0 , then we obtain where we have defined H i,0 ≡ H i (T 0 , P 0 ) and S i,0 ≡ S i (T 0 , P 0 ) for convenience.We note that where C P,i (T, P 0 ) 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) 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 3 .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 non-ideal behavior to be computed.

Heat capacity at constant pressure
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, which are valid for T ≤ 2273 K. Table 1 summarises 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 2 , O 2 , CO, CO 2 , H 2 O, CH 4 , H 2 S and S 2 , all the thermodynamic data are from the extended dataset by Holland & Powell (1998, 2011).For other S-and N-bearing species, i.e., SO 2 , N 2 , NH 3 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 from the JANAF database.
By substituting equation ( 40) into (39), we obtain The preceding equation 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 non-ideally.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 ′ 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 .

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

Equation of state
For non-ideal gases, we use the Compensated-Redlich-Kwong (CORK) equation of state (Holland & Powell 1991) for H 2 O, CO 2 , CH 4 , CO and H 2 .Specifically, this refers to V i (T, P ).For nitrogen-and sulfur-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., page 815 of Symonds & Reed 1993) and G i ≈ G i,0 is obtained.

Fugacity coefficient
With the equations of state of non-ideal gases in hand, one may obtain the fugacity coefficient of a pure gas by numerically evaluating the volume integral (page 186 of DeVoe 2020), where we have set f 0 = P 0 = 1 bar.Since f = ϕ i P , one obtains ϕ i of species i numerically.

Activity coefficient
Holland & Powell (2003) provide fitting functions of activity coefficients (γ i ) for CO-CO 2 -CH 4 -H 2 O-H 2 mixtures, which are valid for 723 ≤ T / K ≤ 2073 and 500 ≤ P/ bar ≤ 4 × 10 4 .In practice, data on activity coefficients are sparser than for fugacity coefficients.Activity data for sulfur-bearing species in the C-H-O-S system are partially available (Evans et al. 2010), but they are unknown when nitrogen-bearing species (N 2 , HCN and NH 3 ) are considered.Since a full treatment of non-ideality 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 non-ideal 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 non-ideality is considered for all the species except for SO 2 .We follow this approach in the full C-O-H-S-N system to account for the full non-ideality 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 non-ideal 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).

RESULTS
Now that we have fully elucidated our theoretical framework for computing the atmospheric chemistry of secondary and hybrid atmospheres, we next explore their anticipated chemical diversity.We first examine carbon-hydrogenoxygen (C-H-O) systems followed by carbon-hydrogenoxygen-nitrogen-sulfur (C-H-O-N-S) chemical systems.The goal is to elucidate what is chemically possible when only considering an atmosphere-melt system in chemical equilibrium.It is understood that, for realistic comparisons to observed systems, one needs to account for how photochemistry may alter the various chemical abundances of some irradiated atmospheres, which is out of scope of the present

Reduced Mantle
Oxidised Mantle Nominal Mantle

Carbon-rich
Figure 4. Examples of secondary atmospheres in the C-H-O chemical system, where the volume mixing ratios (relative abundances by number) of gases are shown as a function of the prescribed atmospheric surface pressure.The top and bottom rows are for low-and high-carbon content in the mantle, respectively.The first, second and third columns are for reduced, nominal and oxidised mantles, respectively.See text for specific parameter values.Regions of the plots where no curves exist are because the computed partial pressures of CO and CO2 exceed the prescribed total pressure, implying that no mathematical solutions exist.Solid and dotted curves correspond to calculations with fully non-ideal effects (see text for details) and the assumption of an ideal gas with ideally-mixed constituents, respectively.4 for secondary atmospheres in the C-H-O chemical system, but with volume mixing ratios as a function of the oxygen fugacity of the mantle.The first, second and third columns are for atmospheric surface pressures of 1 mbar (Mars-like), 1 bar (Earth-like) and 100 bar (Venus-like), respectively.The heterogeneous ranges of values for the oxygen fugacity on the horizontal axes were chosen to minimise displaying regions of parameter space where no mathematical solutions exist.study.We find that methane-dominated atmospheres are rare, which motivates a follow-up investigation into the properties required for their existence.
We present a suite of figures each consisting of a montage of 6 plots with some combination of the following parameter values: low (a C = 10 −7 ) versus high (a C = 0.1) carbon content of the mantle melt; reduced (log f O2 = IW − 3), nominal (log f O2 = IW) and oxidised (log f O2 = IW + 3) mantle melts; P = 1 mbar, 1 bar and 100 bar (inspired by Mars, Earth and Venus); P H2 = 1 mbar, 1 bar and 100 bar (in the absence of a general theory of atmospheric escape).
In plots where we have to pick a single melt temperature, we choose T = 1600 K as explained previously in Section 2.6.
In plots where we have to pick a single oxygen fugacity of the mantle, we fix log f O2 = IW for illustration unless otherwise stated in the caption.Once a combination of parameter values are selected (e.g., a C = 10 −7 , log f O2 = IW − 3, and T = 1600 K in Figure 4a), we explore the chemical trend with an extra parameter (e.g., P in Figure 4a); the various parameter ranges used for examining trends are listed in Table 2. Overall, these choices of parameter values are inspired by and extended beyond the range of Solar System rocky planets, and the overarching intention here is to illustrate qualitative trends rather than model any specific exoplanet.
The choice of these values for the carbon activity are difficult to reconcile with the absolute carbon content in mass, as we do not model mass conservation explicitly.However, some physical interpretation may be provided.A melt saturated with pure graphite corresponds to a C = 1.If a C < 1, then the graphite is undersaturated and thus dissolves into coexisting melts as some form of carbon.Since the activity is the generalisation of the relative abundance by number (volume mixing ratio for gases), one may interpret the carbon activity as the relative abundance of carbon in the melt.
The total atmospheric surface pressure (for secondary atmospheres only; P ), hydrogen partial pressure (for hybrid atmospheres only; P H2 ), helium partial pressure (P He ) and nitrogen partial pressure (P N2 ) encode our ignorance about a set of complex processes (radiative transfer, atmospheric mixing, atmospheric escape, etc), which are out of scope of the present study.
In the numerical solutions for secondary atmospheres, some combinations of parameter values for total pressure P , melt temperature T , carbon activity a C , and oxygen fugacity f O2 do not admit physically realistic solutions, as exemplified by subsequent figures.This is explicable.With the prescribed T , a C and f O2 , the partial pressures of CO and CO 2 can be straightforwardly determined through the first two net reactions in (17), and if the sum of these two partial pressures is larger than the prescribed total pressure P , it becomes impossible for the solved partial pressures to be all positive.Physically, if a C is prescribed to be 1, then the nonexistence of a realistic solution signifies graphite destablisation.This reasoning is also elucidated by equation (8) of French (1966) and is applicable to solutions for the full C-H-O-N-S systems.When numerically solving for hybrid atmospheres, the absence of a physically realistic solution is reflected by the PH 2 rise at fixed fO 2 favors H2O and suppresses CO and CO2 (Figure 8b & c) † Some ranges are changed to allow finding physically realistic solutions that are trend-revealing, e.g., the log fO 2 range in Figure 5d-f.ever-increasing total pressure from iteration to iteration (Section 2.5.4).

C-H-O Chemical Systems
To build initial intuition for the chemistry of secondary and hybrid atmospheres, we consider the simplest non-trivial 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.

Secondary atmospheres
Figures 4, 5 and 6 examine the trends in the relative abundances of H 2 , H 2 O, CO, CO 2 and CH 4 as functions of the atmospheric surface pressure (P ), oxygen fugacity of the mantle (f O2 ) and melt temperature (T ), respectively.An atmosphere deriving from a reduced mantle with low carbon content is H 2 -and H 2 O-dominated, not unlike a gas-giant exoplanet (top-left panel of Figure 4).As the mantle becomes increasingly oxidised (moving from the first to the third column of Figure 4), two important trends appear: • The atmosphere transitions from being dominated by H 2 to being dominated by H 2 O; • The major carbon carrier transitions from being CO to CO 2 .
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 4d).Oxidised mantles further suppress the appearance of methane (Figure 4d-f).Figure 5 confirms the preceding trends.The competition between CO and CO 2 motivates the use of their relative abundances as a diagnostic for the oxygen fugacity.Figure 6 demonstrates that the exact melt temperature assumed needs to be considered carefully, because the abundances of CO, CO 2 and CH 4 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.Non-ideal effects appear to be minimal, unless large pressures are considered.

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 H2 ) becomes large (e.g., Figure 8b & c), the production of CO and CO 2 is suppressed in favor of H 2 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).

C-H-O-N-S Chemical Systems
The addition of nitrogen (N) and sulfur (S) to the chemical system results in 5 more gaseous species: hydrogen sulfide (H 2 S), sulfur dioxide (SO 2 ), ammonia (NH 3 ), molecular nitrogen (N 2 ) and hydrogen cyanide (HCN).This doubles the number of species compared to C-H-O systems.
For illustration, we choose P N2 /P = 0.1 (for secondary atmospheres) and P N2 /P H2 = 0.1 (for hybrid atmospheres), although the chemical trends are also explored for P N2 /P and P N2 /P H2 in the range of 10 −4 -1 (Table 3, and Figures 14 & B6 ).While these choices are arbitrary, it is better than choosing an absolute value for P N2 as we wish to avoid situations where the trends associated with nitrogen-bearing species appear artificially weak because an arbitrarily low value of the nitrogen partial pressure was chosen.
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 log f S2 = PP −10 if a value of log f S2 needs to be fixed (e.g., Figure 10), and separately explore the chemical trends when log f S2 varies between PP −10 and PP (Table 3, and Figures 13 & B5).qualitative trends of the C-H-O system carries over; the H2S/SO2 ratio is sensitive to surface pressure and fO 2 ; the NH2/HCN ratio sensitive to surface pressure but not to fO 2 ; the H2S/SO2 ratio is not sensitive to fS 2 but their abundances are; likewise, the NH2/HCN ratio is not sensitive to PN 2 but their abundances are; same as the secondary atmospheres

Secondary atmospheres
Figures 10, 11, 12, 13 and 14 examine the trends in the relative abundances of H 2 , H 2 O, CO, CO 2 , CH 4 , H 2 S, SO 2 , NH 3 , N 2 and HCN as functions of the atmospheric surface pressure (P ), oxygen fugacity of the mantle (f O2 ), melt temperature (T ), sulfur fugacity of the mantle (f S2 ) and atmospheric partial pressure of molecular nitrogen (P N2 ), respectively.The qualitative trends and lessons learned from examining C-H-O chemical systems carry over.Additionally, the following trends emerge: • The competition between H 2 S and SO 2 depends sensitively on the surface pressure (Figure 10) and oxygen fugacity (Figure 11).
• Likewise, the competition between NH 3 and HCN is sensitive to the surface pressure (Figure 10).As neither NH 3 nor HCN are oxygen carriers, their abundances are somewhat insensitive to the oxygen fugacity (Figure 11).
• The absolute abundances of H 2 S and SO 2 are sensitive to the sulfur fugacity, but the ratio of their abundances is not (Figure 13).
• Likewise, the absolute abundances of NH 3 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 2 S and SO 2 motivates the use of their relative abundances as a supporting diagnostic for the oxygen fugacity.

Hybrid atmospheres
Figures B2, B3, B4, B5 and B6 examine the trends in the relative abundances of H 2 , H 2 O, CO, CO 2 , CH 4 , H 2 S, SO 2 , NH 3 , N 2 and HCN as functions of the atmospheric partial pressure of molecular hydrogen (P H2 ), oxygen fugacity of the mantle (f O2 ), melt temperature (T ), sulfur fugacity of the mantle (f S2 ) and atmospheric partial pressure of molecular nitrogen (P N2 ), respectively.Since no new trends are uncovered beyond what we already learned from examining C-H-O systems and C-H-O-N-S secondary atmospheres, we present these calculations in Appendix B for completeness.

Methane-rich atmospheres
In our explorations of secondary and hybrid atmospheres, both within 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 4 -dominated atmospheres.The trends elucidated in Figures 4-14 and B2-B6 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 −3 ≤ P/ bar ≤ 10 4 (for secondary atmospheres only), 10 −3 ≤ P H2 / bar ≤ 10 4 (for hybrid atmospheres only), 873 ≤ T / K ≤ 1873, IW − 6 ≤ log f O2 ≤ IW + 6, PP−10 ≤ log f S2 ≤ PP, 10 −2 ≤ P N2 /P ≤ 1 (for secondary atmospheres only) and 10 −2 ≤ P N2 /P H2 ≤ 1 (for hybrid atmospheres only).A methane-dominated atmosphere is identified using the criterion, 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 O2 ≲ 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).methane-dominated atmospheres requires low melt temperatures (T ≲ 1000 K), reduced mantles (fO 2 ≲ IW) and high atmospheric surface pressures exceeding ∼ 10 bar.† set aC = 1 when investigating methane-dominated atmospheres.

Comparison to previous work
In terms of its scientific intentions, the closest study to compare 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 non-negligible amounts of molecular hydrogen.As in the present study, Liggins et al. (2020) prescribed or parametrised the oxygen fugacity and total pressure in mixed-phase, chemical-equilibrium, carbon-hydrogenoxygen-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) parametrised 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 2 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 parametrised by a "fudge factor".Liggins et al. (2020) concluded that ∼ 1% H 2 abundances (by number) is possible for early Earth, but ∼ 10% H 2 abundances are unlikely given the plausible space of parameters explored.For early Mars, about 2-8% of H 2 (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 2 O and CO 2 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 parameterised 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 Examples of hybrid atmospheres in the C-H-O chemical system, where the volume mixing ratios (relative abundances by number) of gases are shown as a function of the prescribed atmospheric partial pressure of molecular hydrogen.The top and bottom rows are for low-and high-carbon content in the mantle, respectively.The first, second and third columns are for reduced, nominal and oxidised mantles, respectively.See text for specific parameter values.Regions of the plots where no curves exist are because the computed partial pressures of CH4 and H2O exceed the computed total pressure, implying that no mathematical solutions exist.The oxygen fugacity for panel f is prescribed at IW+1.5 instead of that of oxidised mantle (IW+3) to minimise no-solution parameter space.Solid and dotted curves correspond to calculations with fully non-ideal effects (see text for details) and the assumption of an ideal gas with ideally-mixed constituents, respectively.The solved total pressures for ideal (P-id) and non-ideal (P-nd) cases overlap with each other until PH 2 ≳ 10 3 bar.7 for hybrid atmospheres in the C-H-O chemical system, but with volume mixing ratios as a function of the oxygen fugacity of the mantle.The first, second and third columns are for atmospheric partial pressures of molecular hydrogen of 1 mbar, 1 bar and 100 bar, respectively.The solved total pressures for ideal (P-id) and non-ideal (P-nd) cases overlap with each other until log fO ≳ IW + 2.  8 for hybrid atmospheres in the C-H-O chemical system, but with a mantle of nominal oxidation state (see text for details) and volume mixing ratios as a function of the melt temperature.The solved total pressures for ideal (P-id) and non-ideal (P-nd) cases overlap for the entire T range.
CH 4 and CO 2 as a potential biosignature for exoplanets, which could be further corroborated by absence of abundant CO.Our results can confirm that abundant CH 4 and CO 2 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 4 and CO 2 in chemical equilibrium, rendering simultaneous detection of CH 4 and CO 2 a false positive biosignature for exoplanets.Thompson et al. (2022) explored the range of abiotic CH 4 fluxes in the Solar System, and argues that the short photochemical lifetime of atmospheric CH 4 requires replenishing CH 4 fluxes higher than those from the abiotic outgassing, thus necessitating biological CH 4 fluxes.For an exoplanet, therefore, its planetary context and astrophysical environment needs to be carefully considered to gauge the magnitude of abiotic CH 4 fluxes and then the likelihood of CH 4 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 co-existence with abundant carbon dioxide may be regarded as a biosignature remains highly context-dependent.

Will telescope data allow us to constrain the properties of secondary and hybrid atmospheres?
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 pa-rameters, 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 realisations 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): • Melt temperatures (T ) are sampled uniformly from 873 to 1873 K.
• For the surface pressures (P ) of secondary atmospheres, log (P/bar) is sampled uniformly from −3 to 4 and log (P N2 /P ) is sampled uniformly from −2 to 0. For the hydrogen partial pressures (P H2 ) of hybrid atmospheres, log (P H2 /bar) is sampled uniformly from -3 to 4 and log (P N2 /P H2 ) is sampled uniformly from −2 to 0. The sulfur fugacity is arbitrarily chosen to be log fS 2 = PP − 10.The top and bottom rows are for low-and high-carbon content in the mantle, respectively.The first, second and third columns are for reduced, nominal and oxidised mantles, respectively.See text for specific parameter values.Regions of the plots where no curves exist are because the computed partial pressures of CO, CO2 and H2O exceed the prescribed total pressure, implying that no mathematical solutions exist.The oxygen fugacity for panel f is prescribed at IW+1 instead of that of oxidised mantle (IW+3) to minimise no-solution parameter space.Solid and dotted curves correspond to calculations with partially non-ideal effects (see text for details) and the assumption of an ideal gas with ideally-mixed constituents, respectively.Examples of secondary atmospheres in the C-H-O-N-S chemical system, with volume mixing ratios as a function of the oxygen fugacity of the mantle.The sulfur fugacity is arbitrarily chosen to be log fS 2 = PP − 10.The first, second and third columns are for atmospheric surface pressures of 1 mbar (Mars-like), 1 bar (Earth-like) and 100 bar (Venus-like), respectively.The heterogeneous ranges of values for the oxygen fugacity on the horizontal axes were chosen to minimise displaying regions of parameter space where no mathematical solutions exist.
• For the oxygen fugacity (f O2 ), log f O2 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 log f S2 from PP − 10 to PP.
• The carbon activity (a C ), which is a proxy for the carbon content of the mantle, is sampled uniformly as log a C from −8 to 0.
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 CO2 /X 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 H2S /X SO2 ) is chosen as a supporting diagnostic of the oxidation state.Since water and methane are expected to be readily detectable, their abundance ratio (X H2O /X CH4 ) is chosen as another observational diagnostic.The abundance ratio of ammonia to hydrogen cyanide (X NH3 /X HCN ), which are two important nitrogen carriers, complete the set of 4 observational diagnostics.
As expected, a combination of X CO2 /X CO and X H2S /X SO2 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 4 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 4 abundance ratios.For both types of atmospheres, it appears that the sulfur fugacity and nitrogen content are somewhat challenging to constrain from measuring these 4 abundance ratios.
It is worth noting that, in Cycle 1 of the JWST Guest Observer program4 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.5.3.Are we thinking about the chemistry of rocky exoplanets in a physically realistic way?
In the study of gas-giant exoplanets, astronomers often use the C/O ratio and the "metallicity" as the control parameters.The study of 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 oxidised mantles (e.g., Wadhwa 2008;Cartier & Wood 2019).Our current work suggests that there is a clear path towards deriving the oxygen fugacities (f O2 ) of the rocky mantles of exoplanets by inferring accurate abundances of CO 2 and CO via atmospheric retrieval performed on highquality spectra (by, for example, 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 O2 and the masses of exoplanets for a sample of super Earths and sub-Neptunes.

Opportunities for future work
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 zero-dimensional (0D), meaning that the temperaturepressure 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  11 for secondary atmospheres in the C-H-O-N-S chemical system, but with volume mixing ratios as a function of the melt temperature.The sulfur fugacity is arbitrarily chosen to be log fS 2 = PP−10.For display purposes, a reduced mantle (log fO 2 = IW−3) was arbitrarily chosen because it minimises the regions of parameter space with no mathematical solutions.10 for secondary atmospheres in the C-H-O-N-S chemical system, but with volume mixing ratios as a function of the sulfur fugacity of the mantle.For display purposes, different values of the oxygen fugacity and a fixed total pressure P = 100 bar were chosen because it minimises the regions of parameter space with no mathematical solutions.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 parametrise them by single numbers.In the current study, we have parametrised 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 2 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 parametrise these formidable processes using two numbers (f O2 and P or P H2 ).Outgassed atmosphere composition, partial and total pressures for the simple H-O system with T = 1600 K at the atmosphere-melt interface.(a) Variation of H2O and H2 partial pressures with oxygen fugacity.Left y-axis is for partial pressures Pi and right y-axis is for total pressure P .Dotted curves correspond to dissolution-free (α = 0) cases, whereas solid curves correspond to cases that consider H2O dissolution (α ̸ = 0).(b) Same results as in (a), but replotted to show the variation of mixing ratios (Xi) with oxygen fugacity.Note that cases with and without gas dissolution share the same trend of volume mixing ratios, despite disparate trends of total pressure (magenta curves).
Figure A1b).Moreover, due to the negligible H 2 solubility, the entire hydrogen budget (one Earth ocean) resides exclusively in the atmosphere reservoir in the form of H 2 gas.This enables us to estimate the corresponding P H2 under low oxygen fugacities, i.e., P H2 = M H2,atm g/A ≈ 30 bar, where M H2,atm is equivalent to the prescribed hydrogen budget in mass unit.The estimated P H2 ≈ 30 bar is verified by the numerical solutions for both dissolution and dissolution-free cases in Figure A1a (red curves).In contrast, under high oxygen fugacity conditions (∼IW+6), the dissolution and dissolution-free cases both feature almost pure H 2 O atmosphere (blue curve in Figure A1b).In this case, ignoring H 2 O dissolution means nearly all the hydrogen budget is oxidised to H 2 O and resides exclusively in the atmosphere.This enables estimating P H2O = M H2O,atm g/A ≈ 30 × 9 = 270 bar, where the multiplier 9 is the ratio of molecular weight between H 2 O and H 2 and it stems from the total hydrogen budget being completely oxidised to H 2 O.Such an estimate of P H2O ≈ 270 bar at high f O2 is verified by the numerical solution of the dissolution-free case in Figure A1a (dotted blue curve).When H 2 O dissolution is enabled at high f O2 , its sequestration into coexisting melts greatly reduces its atmospheric partial pressure, as suggested by Figure A1a (compare solid and dotted blue curves).Furthermore, under oxidising conditions, dissolution-induced drop of P H2O drives P H2 to even lower levels, as revealed by comparison of solid and dotted red curves in Figure A1a.Overall, because the H 2 O converted from H 2 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 A1a).
Despite the considerable effect of gas solubility on atmosphere pressures, it is worth emphasising 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 2 O in this study) solubility, as attested to by the fully overlapped curves in Figure A1b.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 2 O and H 2 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 A1b shows that the crossover of X H2O and X H2 occurs at f O2 ∼ 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 B3. Second,as shown by Figure 1 in Ortenzi et al. (2020) andFigures 2 &8 in Sossi et al. (2023), considering or not 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, e.g., CO, CO 2 , are indeed affected by dissolution behaviours (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.

Figure 2 .Figure 3 .
Figure2.Examples of activity coefficients (which we denote by γi in the present study), which quantify departures from ideal mixing of gaseous components.In addition to their dependence on temperature and pressure, activity coefficients depend on the specific mixture of molecules being considered.Here, we show γi contours from H2O-CO2 (top row) and H2O-CH4 (bottom row) binary mixtures.The top and bottom rows show activity coefficients for fixed XCO 2 = 0.05 and XCH 4 = 0.05, respectively, in order to better visualise the multi-dimensional space of possibilities.Note that a linear scale in pressure has been chosen for better visualisation of the contours.Ideal mixing occurs when γi = 1.

Figure 5 .
Figure5.Same as Figure4for secondary atmospheres in the C-H-O chemical system, but with volume mixing ratios as a function of the oxygen fugacity of the mantle.The first, second and third columns are for atmospheric surface pressures of 1 mbar (Mars-like), 1 bar (Earth-like) and 100 bar (Venus-like), respectively.The heterogeneous ranges of values for the oxygen fugacity on the horizontal axes were chosen to minimise displaying regions of parameter space where no mathematical solutions exist.

Figure 6 .
Figure6.Same as Figure5for secondary atmospheres in the C-H-O chemical system, but with volume mixing ratios as a function of the melt temperature.For display purposes, a reduced mantle (log fO 2 = IW − 3) was arbitrarily chosen because it minimises the regions of parameter space with no mathematical solutions.
Figure7.Examples of hybrid atmospheres in the C-H-O chemical system, where the volume mixing ratios (relative abundances by number) of gases are shown as a function of the prescribed atmospheric partial pressure of molecular hydrogen.The top and bottom rows are for low-and high-carbon content in the mantle, respectively.The first, second and third columns are for reduced, nominal and oxidised mantles, respectively.See text for specific parameter values.Regions of the plots where no curves exist are because the computed partial pressures of CH4 and H2O exceed the computed total pressure, implying that no mathematical solutions exist.The oxygen fugacity for panel f is prescribed at IW+1.5 instead of that of oxidised mantle (IW+3) to minimise no-solution parameter space.Solid and dotted curves correspond to calculations with fully non-ideal effects (see text for details) and the assumption of an ideal gas with ideally-mixed constituents, respectively.The solved total pressures for ideal (P-id) and non-ideal (P-nd) cases overlap with each other until PH 2 ≳ 10 3 bar.

Figure 8 .
Figure8.Same as Figure7for hybrid atmospheres in the C-H-O chemical system, but with volume mixing ratios as a function of the oxygen fugacity of the mantle.The first, second and third columns are for atmospheric partial pressures of molecular hydrogen of 1 mbar, 1 bar and 100 bar, respectively.The solved total pressures for ideal (P-id) and non-ideal (P-nd) cases overlap with each other until log fO ≳ IW + 2.

Figure 9 .
Figure9.Same as Figure8for hybrid atmospheres in the C-H-O chemical system, but with a mantle of nominal oxidation state (see text for details) and volume mixing ratios as a function of the melt temperature.The solved total pressures for ideal (P-id) and non-ideal (P-nd) cases overlap for the entire T range.

Figure 10 .
Figure10.Examples of secondary atmospheres in the C-H-O-N-S chemical system, where the volume mixing ratios (relative abundances by number) of gases are shown as a function of the prescribed atmospheric surface pressure.The sulfur fugacity is arbitrarily chosen to be log fS 2 = PP − 10.The top and bottom rows are for low-and high-carbon content in the mantle, respectively.The first, second and third columns are for reduced, nominal and oxidised mantles, respectively.See text for specific parameter values.Regions of the plots where no curves exist are because the computed partial pressures of CO, CO2 and H2O exceed the prescribed total pressure, implying that no mathematical solutions exist.The oxygen fugacity for panel f is prescribed at IW+1 instead of that of oxidised mantle (IW+3) to minimise no-solution parameter space.Solid and dotted curves correspond to calculations with partially non-ideal effects (see text for details) and the assumption of an ideal gas with ideally-mixed constituents, respectively.
Figure11.Examples of secondary atmospheres in the C-H-O-N-S chemical system, with volume mixing ratios as a function of the oxygen fugacity of the mantle.The sulfur fugacity is arbitrarily chosen to be log fS 2 = PP − 10.The first, second and third columns are for atmospheric surface pressures of 1 mbar (Mars-like), 1 bar (Earth-like) and 100 bar (Venus-like), respectively.The heterogeneous ranges of values for the oxygen fugacity on the horizontal axes were chosen to minimise displaying regions of parameter space where no mathematical solutions exist.

Figure 12 .
Figure12.Same as Figure11for secondary atmospheres in the C-H-O-N-S chemical system, but with volume mixing ratios as a function of the melt temperature.The sulfur fugacity is arbitrarily chosen to be log fS 2 = PP−10.For display purposes, a reduced mantle (log fO 2 = IW−3) was arbitrarily chosen because it minimises the regions of parameter space with no mathematical solutions.

Figure 13 .
Figure13.Same as Figure10for secondary atmospheres in the C-H-O-N-S chemical system, but with volume mixing ratios as a function of the sulfur fugacity of the mantle.For display purposes, different values of the oxygen fugacity and a fixed total pressure P = 100 bar were chosen because it minimises the regions of parameter space with no mathematical solutions.
We acknowledge partial financial support from the Chair of Theoretical Astrophysics of Extrasolar Planets at the LMU-Munich and the European Research Council (ERC) via a 2018 Consolidator Grant awarded to KH (number: 771620; acronym: EXOKLEIN).KH is especially grateful to Ralf Bender, Gudrun Niggl and Daniel Grün for their collegial support during his transition from the University of Bern to the LMU.MT and KH jointly designed the current study over dozens of discussions.MT clarified the concepts and derivations associated with thermodynamics, introduced the relevant geochemical literature to KH, curated the thermodynamic data, wrote the computer code necessary to perform the calculations, produced all of the figures and cowrote the paper.KH led the writing and structuring of the paper based on these conversations, introduced the relevant exoplanet literature to MT, derived the generalised outgassing model equations and provided general scientific direction for the paper.

Figure 16 .
Figure 16.Distributions of Monte Carlo outcomes of the atmospheric partial pressure of molecular hydrogen (panel a), melt temperature (panel b), atmospheric surface pressure (panel c), and oxygen fugacity of the mantle (panel d) corresponding to methane-dominated, hybrid atmospheres computed using the C-H-O-N-S system.
B. SUPPLEMENTARY FIGURESIn this appendix, we present for completeness the results for hybrid atmosphere simulations for the C-H-O-N-S system, which reveal no new trends beyond what we already learned from examining C-H-O systems and C-H-O-N-S secondary atmospheres.

Figure B3 .Figure B4 .
Figure B3.Same as FigureB2for hybrid atmospheres in the C-H-O-N-S chemical system, but with volume mixing ratios as a function of the oxygen fugacity of the mantle.The sulfur fugacity is arbitrarily chosen to be log fS 2 = PP − 10.The first, second and third columns are for atmospheric partial pressures of molecular hydrogen of 1 mbar, 1 bar and 100 bar, respectively.The solved total pressures for ideal (P-id) and non-ideal (P-nd) cases overlap with each other until log fO 2 ≳ IW + 2.

Table 1 .
Fitting Coefficients for Computing Heat Capacity and Gibbs Free Energy

Table 2 .
Explored Parameter Ranges † and Key Findings for the C-H-O System

Table 3 .
Explored Parameter Ranges † ‡ and Key Findings for the C-H-O-N-S System

Table 4 .
Monte Carlo Sampled Parameter Space for the C-H-O-N-S System and Key Findings /XCO and XH 2 S/XSO 2 as a strong constraint on the fO 2 of a rocky mantle; complemented by XH 2 O/XCH 4 and XNH 3 /XHCN, melt carbon content and temperature may be further constrained; fS 2and PN 2 are challenging to constrain from the four ratios.