Atmospheres as windows into sub-Neptune interiors: coupled chemistry and structure of hydrogen–silane–water envelopes

ABSTRACT

Sub-Neptune exoplanets are commonly hypothesized to consist of a silicate-rich magma ocean topped by a hydrogen-rich atmosphere. Previous work studying the outgassing of silicate material has demonstrated that such atmosphere–interior interactions can affect the atmosphere’s overall structure and extent. However, these models only considered SiO in an atmosphere of hydrogen gas, without considering chemical reactions between them. Here, we couple calculations of the chemical equilibrium between H, Si, and O species with an atmospheric structure model. We find that substantial amounts of silane, SiH₄, and water, H₂O, are produced by the interaction between the silicate-rich interior and hydrogen-rich atmosphere. These species extend high into the atmosphere, though their abundance is greatest at the hottest, deepest regions. For example, for a 4 M_⊕ planet with an equilibrium temperature of 1000 K, a base temperature of 5000 K, and a 0.1 M_⊕ hydrogen envelope, silicon species and water can comprise 30 per cent of the atmosphere by number at the bottom of the atmosphere. Due to this abundance enhancement, we find that convection is inhibited at temperatures ≳2500 K. This temperature is lower, implying that the resultant non-convective region is thicker, than was found in previous models that did not account for atmospheric chemistry. Our findings show that significant endogenous water is produced by magma–hydrogen interactions alone, without the need to accrete ice-rich material. We discuss the observability of the signatures of atmosphere–interior interaction and directions for future work, including condensate lofting and more complex chemical networks.

convection, planets and satellites: atmospheres, planets and satellites: composition, planets and satellites: gaseous planets, planets and satellites: interiors

1 INTRODUCTION

Exoplanet surveys have revealed that planets with radii between 1 and 4 Earth radii and orbital periods shorter than 100 d are the most intrinsically common type of planet yet observed (e.g. Fressin et al. 2013). Precise mass and radius measurements indicate that these planets are distributed bimodally in radius and bulk density (Weiss & Marcy 2014; Fulton et al. 2017). These measurements separate the small planet population into super-Earths, smaller planets consistent with bulk-Earth composition, and larger sub-Neptunes, which must contain some low-density material to explain their measured radii and densities.

Typically, sub-Neptunes are assumed to contain some combination of terrestrial rock and metal, icy material, and/or hydrogen gas (e.g. Rogers & Seager 2010; Dorn et al. 2017; Zeng et al. 2019). As hydrogen has the lowest density of these three materials, an envelope containing a small amount of it, of the order of 1 per cent of a planet’s total mass, can greatly increase a planet’s size and thus decrease its bulk density. It is therefore possible to model most sub-Neptunes as Earth-like cores with hydrogen envelopes of a few per cent of the planet’s mass (e.g. Lopez & Fortney 2014). However, density measurements do not rule out compositions with large mass fractions in ices, such as water, with correspondingly smaller hydrogen atmospheres (e.g. Zeng et al. 2019; Luque & Pallé 2022).

Due to their ubiquity, a large number of models of sub-Neptune atmospheric evolution and composition have been put forward. Crucially, neither the radii nor masses of these envelopes are expected to remain constant in time. Rather, sub-Neptune atmospheres shrink in extent as the planets cool into space (Lopez & Fortney 2014), and they can be susceptible to atmospheric stripping. This stripping, which can be due to either photoevaporative (Owen & Jackson 2012; Owen & Wu 2017) or core-powered mass-loss mechanisms (Ginzburg, Schlichting & Sari 2016; Gupta & Schlichting 2019), is thought to have stripped some sub-Neptunes entirely, turning them into super-Earths and forming the observed radius valley. The sub-Neptunes that remain were able to mostly resist this stripping. The attributes of the observed radius valley are best matched by mass-loss models if the cores of sub-Neptunes are mostly rocky, with little ice (Gupta & Schlichting 2019; Rogers & Owen 2021).

However, these and other models typically assume that each compositional constituent of the planet is contained in its own layer. Such a structure may be the simplest model to first order, but it may not accurately describe the interactions between these constituents at the high temperatures and pressures expected within sub-Neptunes. Awareness of mixing between layers previously modelled as separate is gaining traction across planetary science. In the Solar system, Jupiter and Saturn show evidence for non-discrete cores that have blended with their metallic H surroundings (Wahl et al. 2017; Mankovich & Fuller 2021). Similar mixing between water and hydrogen has been proposed in the ice giants (Bailey & Stevenson 2021), and water and rock may be miscible at the temperature–pressure conditions of sub-Neptune interiors (Vazan, Sari & Kessel 2022).

Another blurred line between layers arises from the interaction between a hydrogen atmosphere and the potential rocky, silicate-dominated core beneath it. Recent work has shown that at chemical equilibrium, significant silicate vapour is stable in the gas at the base of young sub-Neptune atmospheres (Schlichting & Young 2022), where temperatures can exceed 5000 K (Ginzburg et al. 2016). As silicate vapour will decline in abundance with decreasing temperatures, this implies a compositional gradient deep within sub-Neptunes. This gradient in composition, and thus in mean molecular weight, has been demonstrated to inhibit convection (Markham, Guillot & Stevenson 2022; Misener & Schlichting 2022). The inhibition occurs because a deeper parcel is heavier than one higher up, which overcomes its thermal buoyancy. Such an effect has long been known to apply at the conditions within the Solar system gas and ice giants, where the condensable considered is usually water vapour (Guillot 1995; Leconte et al. 2017; Markham & Stevenson 2021). However, it has only recently been applied to sub-Neptune planets with magma oceans. In addition to arguments from chemical equilibrium, gradients in Si abundance in a hydrogen-dominated atmosphere, and therefore a non-convective region, could also form as a consequence of the accretion of pebbles during formation (Brouwers & Ormel 2020; Ormel, Vazan & Brouwers 2021), though the subsequent evolution of such structures remains unclear.

These studies mark initial forays into understanding the interiors of this ubiquitous class of planet, sub-Neptunes composed of hydrogen and silicate. Much work remains to be done to further quantify the effects of interaction between the interior and atmosphere. In particular, both Misener & Schlichting (2022) and Markham et al. (2022) assume that the gas released into the atmosphere is pure SiO vapour. However, oxidized SiO is not inert in a hydrogen-dominated background gas. Rather, it will react with the hydrogen, producing water (H₂O) and silane (SiH₄). These species will alter the impact of magma condensation on the overall atmospheric structure. It will also change the observable signatures of interior–atmosphere interactions we expect. In this paper, we construct a coupled atmospheric–chemical model that captures the interactions we expect vapourizing silicate magma to have with a hydrogen atmosphere.

2 MODEL

In this section, we detail our atmospheric and chemical model of a sub-Neptune consisting of a silicate-dominated interior and a hydrogen-dominated atmosphere.

2.1 Chemistry

In Misener & Schlichting (2022), to quantify the partial pressure of rock vapour in the atmosphere above a pure SiO₂ magma ocean, the authors use an exponential relationship between SiO and temperature appropriate for congruent evaporation of SiO₂ melt, taken from Visscher & Fegley (2013). The evaporation reaction can be written as

$$\begin{eqnarray} \mathrm{SiO}_{2,{\rm l}} \rightleftharpoons \mathrm{SiO} + 0.5 \mathrm{O}_{2}, \end{eqnarray}$$

(R1)

where on the left-hand side, the l subscript denotes a liquid species; here and going forward, all chemical species without this subscript are gaseous. Reaction (R1) can be quantified using an equilibrium constant:

$$\begin{eqnarray} K_{\mathrm{eq,~R1}} = \frac{P_\mathrm{SiO} P_\mathrm{O_2}^{0.5}}{a_\mathrm{SiO_2}}, \end{eqnarray}$$

(1)

where P_i denotes the partial pressure of species i and |$a_\mathrm{SiO_2}$| is the activity of SiO₂ in the melt, which we take to be 1; i.e. we assume the silicate melt is fully SiO₂. This and following equilibrium constants are found by taking the difference between the Gibbs free energies of the products and reactants. All Gibbs free energy values are calculated from the NIST Chemistry WebBook parametrizations of enthalpy and entropy, based on data from the JANAF tables (Chase 1998). For SiO_2,l, the NIST data extend to a maximum temperature of 4500 K; in this work, we extrapolate the NIST fit to 5000 K, following Schlichting & Young (2022). We also extrapolate the NIST fit below 1996 K, where solid SiO_2,s should be the predominant species. We verify that using the NIST values for the Gibbs free energy of the solid form does not alter our results. The equilibrium constant is a strong function of temperature, where higher temperatures produce relatively more evaporation and thus higher partial pressures of the gaseous species. We demonstrate the temperature dependence of K_eq,R1, as well as the equilibrium constants defined in equations (2) and (3), in Appendix  A. Application of this equilibrium constant produces good agreement with the Visscher & Fegley (2013) equation, though there are slight differences due to the different thermodynamic values used.

However, both species on the right-hand side of reaction (R1) are oxidized, and therefore neither is expected to be inert in the presence of H₂ gas, a strong reducer. Rather, each product should interact with the background hydrogen. Oxygen combines with hydrogen to produce water:

$$\begin{eqnarray} 0.5 \mathrm{O}_2 + \mathrm{H}_{2} \rightleftharpoons \mathrm{H}_2\mathrm{O}, \end{eqnarray}$$

(R2)

while SiO reacts with hydrogen to produce silane, SiH₄:

$$\begin{eqnarray} \mathrm{SiO} + 3 \mathrm{H}_2 \rightleftharpoons \mathrm{SiH}_4 + \mathrm{H}_2\mathrm{O}. \end{eqnarray}$$

(R3)

The production of silane has been suggested to occur deep within Jupiter, at low concentrations (e.g. Fegley & Lodders 1994), due to silane being more thermochemically favourable than SiO at the relevant temperature–pressure conditions (Visscher, Lodders & Fegley 2010b). More recently, the same reaction has been suggested to occur in sub-Neptune atmospheres that interface with magma oceans (Markham et al. 2022), where plentiful oxidized silicate is presumed to exist. Moreover, diamond-anvil cell experiments indicate that mixtures of SiO₂ and H₂ fluid produce silane and water in the fluid at 2 × 10⁴ bar and 1700 K (Shinozaki et al. 2014), and similar silane production was observed in a MgSiO₃–H₂ mixture at 3.6 × 10⁴ bar and 2000 K (Shinozaki et al. 2016). These pressures and temperatures are relevant both to the inner mantle of Earth and to the sub-Neptune conditions we consider here (see Fig. 2).

Figure 1.

Partial pressures of the chemical species present in the atmosphere as a function of radius, in core radii R_c, for a 4 M_⊕ planet with an equilibrium temperature of 1000 K, a base temperature of 5000 K, and an atmospheric hydrogen mass fraction of 2.5  per cent. The species we consider are H₂ (grey), H₂O (blue dashed), SiH₄ (yellow), SiO (red solid), and O₂ (pink). For comparison, the SiO abundance derived from the congruent evaporation equation in Visscher & Fegley (2013) is shown as a red dotted line. The outer vertical black line at 1.19R_c represents the outer radiative convective boundary, while the inner black line at 1.04R_c represents the inner transition within which convection is inhibited.

Figure 2.

Example of sub-Neptune envelope structure. Panel (a) shows the temperature T in kelvin, (b) the total pressure P in bar, and (c) the mean molecular weight μ in proton masses m_p, as functions of radius, in core radii R_c. The model is the same one as shown in Fig. 1: a 4 M_⊕ planet with an equilibrium temperature of 1000 K, a base temperature of 5000 K, and an atmospheric hydrogen mass fraction of 2.5 per cent. In the top panel, the dot represents the outer radiative–convective boundary, while the square represents the inner transition inside of which convection is inhibited.

The relative abundances of these species in the atmosphere are determined by the equilibrium constants for the reactions, respectively:

$$\begin{eqnarray} K_{\mathrm{eq,~R2}} = \frac{P_\mathrm{H_2O}}{P_\mathrm{H_2} P_\mathrm{O_2}^{0.5}} \end{eqnarray}$$

(2)

and

$$\begin{eqnarray} K_{\mathrm{eq,~R3}} = \frac{P_\mathrm{SiH_4} P_\mathrm{H_2O}}{P_\mathrm{SiO} P_\mathrm{H_2}^3}. \end{eqnarray}$$

(3)

As with equation (1), these equilibrium constants are calculated from Gibbs free energies taken from NIST and are functions of temperature.

To calculate the atmospheric composition at each level of the atmosphere, we use the temperature, T, and total pressure, P, derived from the atmospheric structure equations detailed in Section 2.2. We then apply the fact that in our model, the changes in atmospheric composition are driven by the evaporation and condensation of SiO₂. Therefore, we can use the stoichiometric relationship that for every Si atom at a given level in the atmosphere, there must be two O atoms. More quantitatively, we can calculate the ‘partial pressures’ of each element, which both convert through common factors to the numbers at each pressure level,

$$\begin{eqnarray} \sum P_\mathrm{Si} \equiv P_\mathrm{SiO} + P_\mathrm{SiH_4} \end{eqnarray}$$

(4)

and

$$\begin{eqnarray} \sum P_\mathrm{O} \equiv P_\mathrm{SiO} + 2 P_\mathrm{O_2} + P_\mathrm{H_2O}. \end{eqnarray}$$

(5)

We can then mandate that

$$\begin{eqnarray} \sum P_\mathrm{O} = 2 \sum P_\mathrm{Si} \end{eqnarray}$$

(6)

at each level of the atmosphere. Equations (1), (2), (3), and (6) are sufficient to solve for the partial pressures of each component. We detail our analytic method in Appendix  B.

2.2 Atmospheric structure

This atmospheric structure model builds on Misener & Schlichting (2022), which added consideration of compositional gradients and moist convection to the model of Misener & Schlichting (2021). Here, as in both of those works, we model the outer atmosphere as being in thermal equilibrium with the incident flux from the star and thus isothermal at the planet’s equilibrium temperature, T_eq (e.g. Lee & Chiang 2015; Ginzburg et al. 2016).

The atmosphere transitions to convective at the outer radiative–convective boundary (R_rcb). The radial pressure gradient is found following hydrostatic equilibrium:

$$\begin{eqnarray} \frac{\partial P}{\partial R} = -\frac{G M_\mathrm{c}}{R^2} \frac{\mu P}{k_\mathrm{B} T}, \end{eqnarray}$$

(7)

where G is the gravitational constant, M_c is the planet mass, k_B is the Boltzmann constant, and μ is the local mean molecular weight. The mean molecular weight is a function of temperature and pressure, and is calculated using chemical equilibrium, as detailed in the previous section.

The temperature profile in the convective region follows a moist adiabat:

$$\begin{eqnarray} \frac{\partial \ln T}{\partial \ln P} = \frac{k_\mathrm{B}}{\mu } \cfrac{1+\cfrac{P_\mathrm{Si}}{P_\mathrm{H}} \cfrac{\partial \ln P_\mathrm{Si}}{\partial T}}{c_\mathrm{p} +\cfrac{P_\mathrm{Si}}{P_\mathrm{H}} \cfrac{k_\mathrm{B}}{\mu } T^2 \bigg (\cfrac{\partial \ln P_\mathrm{Si}}{\partial T}\bigg)^2}, \end{eqnarray}$$

(8)

where P_H is the partial pressure of hydrogen and P_Si is the combined partial pressures of Si-bearing species (e.g. Leconte et al. 2017). The determination of these partial pressures via chemical equilibrium was described in Section 2.1. The specific heat capacity c_p is a mass-weighted average of the specific heat capacities of each component:

$$\begin{eqnarray} c_\mathrm{p} = \sum _i q_i c_{\mathrm{p}, i}, \end{eqnarray}$$

(9)

where q_i is the mass mixing ratio, defined as q_i ≡ μ_iP_i/(μP), with μ_i being the molecular weight of species i. c_p,i, the heat capacity of each species, is given by

$$\begin{eqnarray} c_{\mathrm{p}, i} = \frac{k_\mathrm{B}}{\mu _i}\frac{\gamma _i}{\gamma _i-1}. \end{eqnarray}$$

(10)

where γ_i is the adiabatic index of a species. For the diatomic molecules O₂, H₂, and SiO, we assume γ = 7/5, while for H₂O we take γ = 4/3 and for SiH₄ we take γ = 1.3, following that of similarly structured methane. We use fixed values of each species’ heat capacity for simplicity; we verify that temperature-dependent values from NIST (Chase 1998) deviate by less than a factor of 2 from these constant values over the adiabatic region, which only marginally alters the atmospheric structure we obtain. Generally, we find that the atmospheric composition in the adiabatic regime is hydrogen-dominated, so the heat capacity is very near that of hydrogen alone. We assume the hydrogen remains fully diatomic throughout the atmosphere. Hydrogen dissociation is expected at sufficiently high pressures and temperatures. Molecular dynamics simulations indicate that hydrogen likely remains diatomic throughout most of our atmospheric structure, though it may start to dissociate at the highest temperatures and pressures we consider (e.g. Tamblyn & Bonev 2010; French et al. 2012; Soubiran et al. 2017). If hydrogen begins to dissociate, the effective adiabatic index of the atmosphere would be lower due to the energy required for breaking the H–H bonds (Lee & Chiang 2015), similar to the effect of latent heat of vapourization, making the adiabatic temperature gradient shallower.

To construct the atmospheric structure, we begin from the outer radiative–convective boundary, where R = R_rcb and P(R_rcb) are chosen, and T(R_rcb) = T_eq. We increment in small pressure steps: P_new = P + ΔP. The new radius is thus |$R_\mathrm{new} = R + \Delta P/(\partial P/\partial R)$|⁠, and the new temperature, |$T_\mathrm{new} = T + \Delta P (\partial T/\partial P)$|⁠, employing equations (7) and (8), respectively.

As described in Misener & Schlichting (2022), the change in composition of the atmosphere with temperature due to condensation causes a mean molecular weight gradient. In a hydrogen-dominated atmosphere, these condensation effects cause the mean molecular weight to increase with temperature, stabilizing the gas against convection. This is in contrast to an Earth-like atmosphere, in which the major condensable, water, is lighter than the background air. In previous work, in which one condensable species was considered, the tipping point beyond which convection is inhibited could be quantified by a critical mass mixing ratio q_crit (Guillot 1995; Leconte et al. 2017; Markham et al. 2022; Misener & Schlichting 2022):

$$\begin{eqnarray} q_\mathrm{crit} = \cfrac{1}{\bigg (1-\cfrac{\mu _\mathrm{H}}{\mu _\mathrm{cond}}\bigg) \cfrac{\partial \ln P_\mathrm{cond}}{\partial \ln T}}. \end{eqnarray}$$

(11)

In an atmosphere with multiple species changing abundance, this relationship is in principle more complex. However, the inhibition of convection can be quantified by calculating a ‘convection criterion’ for each non-hydrogen species χ_i:

$$\begin{eqnarray} \chi _i = q_i \bigg (1-\frac{\mu _\mathrm{H}}{\mu _{i}}\bigg) \frac{\partial \ln P_{i}}{\partial \ln T}, \end{eqnarray}$$

(12)

where q_i is the mass-mixing ratio of species i. Convection is inhibited if ∑_iχ_i ≥ 1. In the case of a single species, this method yields the same results as equation (11). We derive the fact that the overall stability criterion is the sum of the criteria for each species in Appendix  C.

2.2.1 Non-convective region

In regions where the criterion is fulfilled, convection is inhibited, no matter how super-adiabatic the temperature profile becomes. In these regions, heat must be transported by either radiation or conduction. At the typical temperatures and pressures in the deep non-convective regions of sub-Neptunes we consider in this work, conduction may be competitive with radiation in transporting heat (e.g. Vazan & Helled 2020; Misener & Schlichting 2022).

We quantify the competition between conduction and radiation by comparing the conductivity, |$\lambda _\mathrm{cond} = L/(4 \pi r^2)/(\partial T/\partial r)_\mathrm{cond}$|⁠, to the equivalent radiative heat transport term, |$\lambda _\mathrm{rad} \equiv L/(4 \pi r^2)/(\partial T/\partial r)_\mathrm{rad} = 16 \sigma T^3/(3 \kappa \rho)$|⁠, where σ is the Stefan–Boltzmann constant and κ is the local Rosseland mean opacity. Conduction dominates radiation when Λ ≡ λ_cond/λ_rad > 1, or, written in scaling form:

$$\begin{eqnarray} \Lambda &\approx \bigg (\frac{\lambda _\mathrm{cond}}{7 \times 10^4~\mathrm{erg~s^{-1} cm^{-1} K^{-1}}}\bigg) \bigg (\frac{\kappa }{10^3~\mathrm{cm^2 g^{-1}}}\bigg) \\ & \bigg (\frac{\rho }{1~\mathrm{g~cm^{-3}}}\bigg) \bigg (\frac{T}{6000~\mathrm{K}}\bigg)^{-3} \\ &\gt 1 \mathrm{.} \end{eqnarray}$$

(13)

From equation (13), it is apparent that larger conductivities, opacities, and gas densities favour conduction over radiation.

Due to the widely varying atmospheric compositions, the atmospheric Rosseland mean opacity κ is difficult to determine without a detailed radiative model, which is beyond the scope of this work. Therefore, following Misener & Schlichting (2022) and Markham et al. (2022), we use the Freedman et al. (2014) relation for solar metallicity gas throughout the atmosphere. We extend these opacities to pressures beyond their asserted validity, so we acknowledge that the opacity of the interior is a source of uncertainty in our model. We discuss the relevance of different opacity choices in Section 3.

As with the opacities, the conductivities in these regions are uncertain and depend on the material properties of the atmosphere, which are not entirely clear for the exotic mixtures we encounter. Therefore, we employ a simple approach based on the behaviour of pure H/He, as described in Misener & Schlichting (2022). To calculate the thermal conductivity, we use the electrical conductivity scaling with temperature and pressure from McWilliams et al. (2016), which is based on experimental results for pure hydrogen. We then convert these to thermal conductivities using the Wiedemann–Franz law. At relatively low temperatures and pressures, where the electrical conductivity is low, we use a minimum value of λ_cond = 2 × 10⁵ erg s⁻¹ cm⁻¹ K⁻¹, appropriate for the nucleic contribution over a broad range of relevant temperatures and pressures (French et al. 2012) and consistent with approximations used in previous work modelling Earth-like silicate in planetary interiors (e.g. Stevenson, Spohn & Schubert 1983; Vazan & Helled 2020).

We can incorporate both conduction and radiation by adding the thermal conductivity to the equivalent radiative term to make an ‘effective conductivity,’ λ_eff ≡ λ_cond + λ_rad. An alternative, but equivalent, framing is to employ an effective opacity, κ_eff (Vazan & Helled 2020):

$$\begin{eqnarray} \kappa _\mathrm{eff} = \cfrac{1}{\cfrac{1}{\kappa } + \cfrac{1}{\kappa _\mathrm{cond}}}, \end{eqnarray}$$

(14)

where κ_cond, the ‘conductive opacity,’ is given by

$$\begin{eqnarray} \kappa _\mathrm{cond} = \frac{16 \sigma T^3}{3 \rho \lambda _\mathrm{cond}}. \end{eqnarray}$$

(15)

The temperature gradient of the non-convective region is then determined by the energy flux, L, as well as the local pressure, temperature, and the effective opacity κ_eff:

$$\begin{eqnarray} \frac{\partial \ln T}{\partial \ln P} = \frac{3 \kappa _\mathrm{eff} P L}{64\pi G M_\mathrm{c} \sigma T^4}. \end{eqnarray}$$

(16)

In steady state, the energy flux that must be transported across the radiative boundary is equal to the radiative flux of the planet into space, i.e. the luminosity at the radiative–convective boundary:

$$\begin{eqnarray} L = \frac{\gamma -1}{\gamma }\frac{64\pi G M_\mathrm{c} \sigma T_\mathrm{eq}^4}{3 \kappa _\mathrm{rcb} P_\mathrm{rcb}}, \end{eqnarray}$$

(17)

where ‘rcb’ subscripts represent values calculated at the conditions of the radiative–convective boundary. Here, the atmosphere is dominated by hydrogen, so we use the adiabatic index for pure H₂ gas and opacity relations of Freedman et al. (2014) for solar composition gas.

3 RESULTS

In this section, we present the atmospheric profile and chemical abundances we obtain for our fiducial planet, a 4 M_⊕ planet with an equilibrium temperature of 1000 K. The atmosphere has a total hydrogen mass of 2.5  per cent the core’s total mass. These values are all typical of sub-Neptunes (e.g. Lopez & Fortney 2014). This planet has a base temperature of 5000 K, which is a reasonable value early in the planet’s evolution (e.g. Ginzburg et al. 2016; Misener & Schlichting 2022). We begin by examining the region interior to the outer radiative–convective boundary in Section 3.1, where the chemistry has the largest effect on the atmospheric structure. We then describe the chemical equilibrium of the outer radiative region in Section 3.2, which has implications for the observability of these interior–atmosphere interactions.

3.1 Inner atmosphere

In this section, we focus on the region interior to the outer radiative–convective boundary, which is at ∼1.19R_c in our fiducial model. This radius and corresponding pressure P(R_rcb) are found by iterating the atmospheric profile until the desired planet characteristics, i.e. the hydrogen mass and base temperature T(R_c), are achieved, following the method of Misener & Schlichting (2022). In Fig. 1, we show the partial pressures of the species we consider, namely H₂ (grey), H₂O (blue dashed), SiH₄ (yellow), SiO (red solid), and O₂ (pink), as functions of radius. In Fig. 2, we show aspects of the overall atmospheric profile, namely the temperature, pressure, and mean molecular weight as functions of radius, and in Fig. 3, we present an alternative view of the same model as in Fig. 1, but in pressure–temperature space.

Figure 3.

Partial pressures of the chemical species present in the atmosphere as a function of temperature. The line meanings remain the same as those in Fig. 1, as do the physical parameters: a 4 M_⊕ planet with an equilibrium temperature of 1000 K, a base temperature of 5000 K, and an atmospheric hydrogen mass fraction of 2.5 per cent. The black line at 2300 K represents the point at which convection becomes inhibited at hotter temperatures.

The atmosphere is hydrogen-dominated by number at all radii (i.e. |$P_\mathrm{H_2} \gt P_i$| for all other species i). The dominance of hydrogen by number is demonstrated in Fig. 4, which shows the number fraction of the species we consider throughout the atmosphere. The abundances of all the secondary species we consider, namely H₂O, SiH₄, SiO, and O₂, increase in abundance with depth as the temperature increases. Water and silane have partial pressures a factor of ∼10⁻⁷ lower than that of H₂ at the outer radiative–convective boundary, while SiO is lower still, with a number fraction of ∼10⁻¹⁴, and O₂ many orders of magnitude less than this. However, at the magma-atmosphere interface, r = R_c, secondary species abundances are much higher, with SiO, H₂O, and SiH₄ together making up nearly 40 per cent of the atmosphere by number and overtaking H₂ by mass, as shown in Fig. 5. Accordingly, the mean molecular weight, displayed in panel (c) of Fig. 2, remains near 2 m_p, that of H₂, throughout most of the atmosphere but rises sharply to larger than 8 m_p near the inner edge.

Figure 4.

The number fraction of the species we consider as a function of radius, in core radii R_c. The line meanings remain the same as those in Fig. 1, as do the physical parameters: a 4 M_⊕ planet with an equilibrium temperature of 1000 K, a base temperature of 5000 K, and an atmospheric hydrogen mass fraction of 2.5 per cent. The number fraction of oxygen remains less than 10⁻⁴ and so is not shown. The atmosphere remains hydrogen dominated by number.

Figure 5.

The mass fraction of the different species as a function of radius, in core radii R_c. The line meanings remain the same as those in Fig. 1, as do the physical parameters: a 4 M_⊕ planet with an equilibrium temperature of 1000 K, a base temperature of 5000 K, and an atmospheric hydrogen mass fraction of 2.5  per cent. The mass fraction of oxygen remains less than 10⁻⁴ and so is not shown. Near the base of the atmosphere, when T ∼ 5000 K, the atmosphere becomes dominated in mass by water and SiO.

We show the number ratios of H₂O to SiO (in purple) and SiH₄ to SiO (in orange) in Fig. 6. Throughout most of the atmosphere, SiH₄ dominates over SiO, except at the very base for our chosen parameters, and H₂O dominates SiO in the entire atmosphere. This dominance of reduced species compared to oxidized ones is due to the presence of hydrogen gas as the background species. In contrast to inert background gases such as nitrogen, hydrogen is highly reactive and will act to reduce the outgassed magma ocean species, per reactions (R2) and (R3). We find these equations are typically driven towards the products at the temperatures and hydrogen gas fractions we consider, favouring the production of water and silane as the primary oxygen and silicon carriers in the atmosphere, respectively. This is in contrast to previous work on the atmosphere–interior interaction, such as Misener & Schlichting (2022) and Markham et al. (2022), which considered only SiO as the main carrier of the silicon and oxygen taken up from the underlying magma ocean.

Figure 6.

Number ratios of H₂O to SiO (in purple) and SiH₄ to SiO (in orange) as a function of radius, in core radii R_c. The model is the same one as shown in Fig. 1: a 4 M_⊕ planet with an equilibrium temperature of 1000 K, a base temperature of 5000 K, and an atmospheric hydrogen mass fraction of 2.5 per cent. The more reduced species, H₂O and SiH₄, dominate at nearly all points in the atmosphere, except near the hot base.

Due to the numbers of silicon and oxygen atoms remaining in a fixed ratio, as described by equation (6), the partial pressures of water and silane are similar, and remain in lockstep as the temperatures increase. At sufficiently high temperatures deep in the atmosphere, SiH₄ production is no longer favoured over SiO, and the abundance of SiO approaches and, for the conditions we consider, just surpasses that of SiH₄, as shown most clearly in Fig. 4. This behaviour appears qualitatively similar to the CH₄–CO transition, extensively studied in the context of exoplanet observations (e.g. Burrows & Sharp 1999; Visscher, Moses & Saslow 2010a; Fortney et al. 2020), though that transition occurs at much lower temperatures (see also Section 3.2).

Another notable result shown by Fig. 1 is that the overall abundance of silicon-bearing species is much higher than was previously found in Misener & Schlichting (2022). In that work, the authors used a vapour pressure formula appropriate for congruent evaporation of silicate magma into a vacuum taken from Visscher & Fegley (2013). This value for rock vapour is denoted in Fig. 1 with a red dotted line. However, Misener & Schlichting (2022) ignored the chemical effects of the non-inert background gas, H₂. When we account for the hydrogen chemistry, we find that the production of silane and water via reaction (R3) draws out more SiO from the interior in order to continue to satisfy the relationship with K_eq,R3, which is solely a function of temperature, in equation (3). In order to maintain the equality of equation (1), the oxygen partial pressure decreases to the values in Fig. 1. These values are much lower than the congruent ratio of 0.5 mole of O₂ for every mole of SiO.

Meanwhile, water is the most abundant product of magma–hydrogen interaction in our model, despite no water being present in the system initially. Water as a fundamental by-product of silicate-hydrogen chemistry has previously been found in the context of sub-Neptunes (Schlichting & Young 2022; Zilinskas et al. 2023) and early Earth (Young, Shahar & Schlichting 2023). This result shows that the presence of water in a sub-Neptune atmosphere does not necessarily indicate it formed with substantial ices, as has been suggested for some sub-Neptunes (e.g. Zeng et al. 2019; Madhusudhan et al. 2020; Venturini et al. 2020; Emsenhuber et al. 2021). We investigate factors that could alter the abundance of water in a sub-Neptune in Section 4.2.

The increased abundances of condensable species alter the overall atmospheric structure compared to the previous SiO-only model of Misener & Schlichting (2022). These increased abundances increase the value of the convective criterion, and therefore lower the temperature at which the atmosphere transitions from convective to radiative to ∼2300 K, as marked in Fig. 2 by the square, a temperature significantly lower than the typical transition temperature of ∼4000 K found in Misener & Schlichting (2022) for similar planet parameters.

A transition at lower temperatures means the atmosphere is non-convective for a larger range of temperatures in the deep interior of sub-Neptune planets, as can be seen in Fig. 2. The gas density at the transition point is also lower. Since the opacities and conductivities we expect are lower at lower densities, the temperature–pressure profiles are not as steep in the non-convective region as was found in Misener & Schlichting (2022). We compare the opacities we use in Fig. 7. We find that radiation and conduction are comparably efficient in transporting energy deep in sub-Neptune atmospheres, as found in Misener & Schlichting (2022), with conduction dominating as the temperatures and pressures increase. As discussed in Section 2.2.1, the opacity we use is extrapolated beyond the pressure regime fit in Freedman et al. (2014). However, it is apparent in Fig. 7 that the effective opacity is close to the conductive opacity in the interior, with the radiative opacity being larger. If the opacity were larger than the solar value we use (∼10³ cm² g⁻¹), as might be expected in a hot, high-metallicity region, heat transport would be even more conduction dominated than we find here, with little effect on the structure we obtain. We therefore conclude that our results are insensitive to the exact value of the opacity in the interior so long as it is higher than the conductive opacity. This is similar to the conclusion reached regarding opacities in Misener & Schlichting (2022). Fig. 7 also confirms that radiation dominates conduction at the outer radiative–convective boundary.

Figure 7.

Comparison of the radiative opacity κ_rad (in blue), the conductive opacity κ_cond (in orange), and the effective opacity κ_eff (in green), as a function of radius, in core radii R_c. The model is the same one as shown in Fig. 1: a 4 M_⊕ planet with an equilibrium temperature of 1000 K, a base temperature of 5000 K, and an atmospheric hydrogen mass fraction of 2.5 per cent. The black line represents the inner non-convective boundary. Radiation dominates at the outer radiative–convective boundary, as expected, while conduction becomes more important than radiation in the interior, though by a factor of 10 or less.

3.2 Outer atmosphere

We now examine the implications of this interior chemistry and structure on the outer atmosphere, the region accessible to spectroscopic observations. We model the outer atmosphere as isothermal at T = T_eq, which for our model planet is set to 1000 K. In this region, the total pressure falls off exponentially. We assume chemical equilibrium is maintained throughout the atmosphere and that the vapour is saturated in SiO₂, with the SiO₂ activity always equal to one. In particular, by assuming chemical equilibrium we implicitly assume that the region we probe is below the homopause, above which each species follows its own scale height, and that mixing is always faster than the chemical kinetic time-scale, i.e. the reactions we consider are never quenched. In Fig. 8, we show the partial pressures of the chemical species we consider as a function of planet radius over the whole atmosphere. The interior region, r < R_rcb ≈ 1.19R_c, contains the same abundances presented in Fig. 1.

Figure 8.

Partial pressures of the chemical species present in the atmosphere as a function of radius, in core radii R_c, through the entire atmosphere out to P = 10⁻⁶ bar (r ∼ 3R_c), including the outer radiative region. The profile is for a 4 M_⊕ planet with an equilibrium temperature of 1000 K, a base temperature of 5000 K, and an atmospheric hydrogen mass fraction of 2.5 per cent. The species we consider are H₂ (grey), H₂O (blue dashed), SiH₄ (yellow), SiO (red solid), and O₂ (pink). The outer vertical black line at 1.19 R_c represents the outer radiative convective boundary.

Fig. 8 shows sharp kinks in the abundances of the species we consider as the temperature ceases declining. Intriguingly, not only do the abundances change with altitude in the isothermal region, but the relative proportions of the different species do as well. The changes in relative abundances is shown more clearly in Fig. 9, which displays the number fraction of each species as a function of the total pressure in the isothermal region only. This total pressure is very nearly the pressure of H₂, as evidenced by its number fraction being nearly 1 in this region. The relative abundance of SiO increases throughout the isothermal region, while the relative abundance of SiH₄ decreases. The two become equal near 10⁻¹ bar: at lower pressures (i.e. higher altitudes) than this, SiO is the dominant silicon-bearing species, rather than SiH₄.

Figure 9.

Number fraction of each chemical species as a function of total pressure in the outer isothermal region. The temperature is fixed at 1000 K. The number fraction of O₂ remains less than 10⁻²⁰ and so is not shown. SiO becomes the dominant silicon-bearing species over SiH₄ at ∼10⁻¹ bar.

The increasing dominance of SiO over SiH₄ can be understood by examining equation (3). At constant temperature, K_eq,R3 is constant. As |$P_\mathrm{H_2}$| decreases with altitude, the partial pressure of SiO must increase to maintain the equilibrium. In other words, the reaction no longer so strongly favours the production of SiH₄. Due to the fixed Si:O atomic number ratio, equation (6), the abundance of water first decreases with decreasing pressure, tied to the abundance of SiH₄, then increases in lockstep with SiO once the latter becomes dominant.

This transition from SiH₄ to SiO as the dominant silicon-bearing species occurs at an observationally relevant pressure for the planet parameters studied here. Therefore, the relative SiO/SiH₄ abundance could serve as a probe of the overall chemistry of the interior and atmosphere, similarly to the well-studied CO–CH₄ transition, which occurs due to a similar reaction (e.g. Burrows & Sharp 1999; Visscher et al. 2010a; Fortney et al. 2020). Future work will elucidate how this transition point varies across the sub-Neptune parameter space and in time, as well as with more complex interior compositions. It will also further constrain the absolute abundances we expect. We note that our simple model assumes the initial atmosphere comprises pure hydrogen, such that all heavier species in our results are due to outgassing from the magma interior. This results in elemental abundances in the upper atmosphere which are sub-solar, and which do not match solar ratios. Specifically, the solar abundance of silicon is 3.3 × 10⁻⁵N_H, and that of oxygen is 5.7 × 10⁻⁴N_H, where N_H is the hydrogen abundance (Lodders 2021), implying an oxygen-to-silicon ratio of N_O/N_Si = 17.3. Our abundances and oxygen-to-silicon ratio are lower than these values throughout the isothermal region. We further discuss the effects of different compositions of core and atmosphere on our results in Section 4.2.

4 DISCUSSION AND FUTURE DIRECTIONS

In Section 3, we demonstrated that silane and water are expected to be by-products of silicate–hydrogen interactions in sub-Neptunes with underlying magma oceans. In this section, we discuss the observability of these species with current and future facilities. We then discuss how other species beyond the simplified chemical network we consider here could alter the atmospheric profiles we obtain, how these new profiles could alter the evolution in time of these atmospheres, and how our assumptions about opacities could impact these results.

4.1 Observability

The two most abundant species in our model besides hydrogen are water and silane. Water vapour has numerous strong absorption bands in the infrared, which JWST is well suited to exploit. Due to its potential as a tracer of planet formation and its importance for life, water vapour has been extensively searched for in the atmospheres of sub-Neptunes in previous campaigns with Hubble and JWST. Water features have been detected in Hubble observations of sub-Neptunes (Benneke et al. 2019a, b) and in the atmospheres of giant exoplanets using JWST (e.g. Alderson et al. 2023; Rustamkulov et al. 2023), with observations of smaller planets underway. Our results show that detections of water vapour in sub-Neptune atmospheres are not necessarily diagnostic of formation beyond the snow line. Instead, it could be produced endogenously via magma–hydrogen interactions.

Silane also has features in the mid-infrared, with the largest cross-section at ∼4.5 μm (Owens et al. 2017) according to the ExoMol data base (Tennyson et al. 2016). This feature is squarely within the wavelength capabilities of JWST’s NIRSpec and PRISM instruments, as evidenced by the detections of features at similar wavelengths, such as CO₂ in the hot Jupiter WASP-39 b (JWST Transiting Exoplanet Community Early Release Science Team 2023). Silane has previously been considered in small exoplanet atmospheres due to its potential as a biosignature photosynthetic product in reducing conditions, though it was deemed unlikely to occur (Seager, Bains & Hu 2013). In the case studied in this work, silane would instead be a product of magma–atmosphere interactions.

We assume throughout this work that the condensate, in our case liquid SiO₂, is present in small enough quantities that it does not affect the atmospheric properties, such as, e.g. its opacity or heat capacity. However, some liquid must remain suspended in the gas. Lofting of condensates could alter the atmospheric profile, due to the liquid’s larger specific heat (Graham et al. 2021). Such condensate retention would also by definition produce clouds, which we do not explicitly model but which may have dramatic effects on the observability of these features. Therefore, such lofting should be investigated further, despite its reliance on complex microphysical processes. Possibly helpful analogues include more massive bodies such as hot Jupiters and brown dwarfs, in which silicate clouds have been previously studied (e.g. Burningham et al. 2021; Gao & Powell 2021).

4.2 Other species and compositions

We considered a chemical network including three reactions, in order to demonstrate the potential importance of these reactions on the atmospheric structure. Chemical networks involving more species and reactions have been employed in the study of sub-Neptune atmospheres in general, as well as magma–atmosphere interactions in various planetary contexts (e.g. Moses et al. 2013; Schlichting & Young 2022; Zilinskas et al. 2023), though these works did not couple the products of silicate–hydrogen reactions with atmospheric structure in the sub-Neptune regime as we have here. In this section, we highlight possible extensions of the chemistry we consider that could influence the structure and evolution of sub-Neptune planets.

First, we considered an interior composed of pure SiO₂. However, including a more realistic, Earth-like composition would alter the products outgassed into the atmosphere. For example, if the rocky interior were approximated as pure MgSiO₃ instead, we would expect the evaporation products to be atomic Mg, SiO, and O₂, in equal molar proportions. This increase in the proportion of oxygen to silicon would likely lead to increased water production compared to our model. However, experimental results indicate that silicon may preferentially dissolve into H₂ compared to magnesium (Shinozaki et al. 2013, 2016), potentially affecting, e.g. the Mg/Si ratio we would predict. We also do not consider that more complex chemistry in the melts will lower the activities of the relevant melt species and thus affect the equilibrium vapour pressures of SiO and O₂, which may in turn alter the physics of convection inhibition. This example illustrates that interior composition can alter the atmospheric composition of sub-Neptunes, but that such interactions may be complex.

We do not expect an ice layer to form at the base of the atmosphere. Extrapolations indicate water and hydrogen are likely miscible at these temperatures and pressures, though experiments are needed to confirm this (Bailey & Stevenson 2021). Additionally, solid phases of water are not stable at the temperatures and pressures at the base of the atmosphere (e.g. Madhusudhan et al. 2020). In fact, the temperatures and pressures are above the critical temperature and pressure of pure water, though the bulk behaviour depends on the properties of the hydrogen–water–silane mixture, which is poorly constrained, not those of water alone (Markham et al. 2022). The critical temperature of the silicate interior is likely higher than the 5000 K base temperature we assume (e.g. Xiao & Stixrude 2018), indicating the overall mixture may not be supercritical, but more modelling of such mixtures is needed. If the mixture is supercritical, condensation of silicate vapour can no longer occur, and so convection would no longer be inhibited in the supercritical region of the envelope (Markham et al. 2022; Pierrehumbert 2023).

While not yet supercritical, the incompressibility of the melt and gases at high pressures can lead to changes in the Gibbs free energies of formation not captured in our equilibrium model, which assumes ideal gases and melt behaviour. The equation of state of silicate melt is relatively well constrained (e.g. De Koker & Stixrude 2009), allowing calculation of its molar volume as a function of pressure and therefore its change in chemical potential (e.g. Schlichting & Young 2022). We find the change in chemical potential of the melt at the base of our atmosphere to be approximately 230 kJ mol⁻¹, enough to significantly alter the equilibrium state. However, this change in chemical potential may be compensated by non-ideal behaviour in the product gaseous species on the other side of reaction (R1). Specifically, an increase in the fugacity coefficient of SiO of approximately 40–100 at the conditions of the base of the atmosphere (P ∼ 10⁵ bar, T ∼ 5000 K) would balance the PV effects on the chemical potential of the melt. Unfortunately, the equations of state, and thus the fugacities, of the silicon vapour species we consider are not well constrained. The observed behaviour of other species, such as water vapour, indicates that such fugacities are plausible at the pressure conditions we consider (e.g. Otsuka & Karato 2011). However, their precise determination is beyond the scope of this work.

We also ignore ingassing reactions, such as the solubility of hydrogen and water in the interior. Such reactions may be important in driving the interior composition and long-term atmospheric evolution of sub-Neptunes and super-Earths, although the relevant solubilities are highly uncertain (e.g. Chachan & Stevenson 2018; Olson & Sharp 2019; Dorn & Lichtenberg 2021; Schlichting & Young 2022).

Finally, we assume that the entire atmosphere, i.e. a given mass of hydrogen gas, has chemically equilibrated with the underlying magma ocean, similar to previous work on magma–hydrogen interactions (e.g. Markham et al. 2022; Schlichting & Young 2022; Zilinskas et al. 2023). This chemical processing leads to the subsolar abundances we predict in the outer atmosphere: most of the oxygen and silicon is segregated to the deep atmosphere where it is more thermodynamically favourable. The result is an atmosphere that is overall strongly supersolar in e.g. oxygen abundance – O makes up 2 per cent of the total atmosphere by number – but with strong variations with depth. The likelihood of reaching full chemical equilibrium depends on the mixing efficiency in the atmosphere during and after formation. If there is sufficient transport between atmospheric layers over the age of the planet, then we would expect chemical equilibrium to be achieved throughout the whole planet system, as modelled here. Conversely, higher elemental abundances of Si and O in the outer atmosphere than predicted here could indicate that chemical equilibrium between the outer and inner layers of the planet has not been reached, and therefore that the interior and atmosphere are poorly coupled chemically (e.g. Zilinskas et al. 2023). Enhancement in these species could come from a higher metallicity in the initial accreted gas, although some of these enhanced metals would have condensed into refractory materials. Another source could be the ablation of accreting solid material; however, ablation is thought to occur in deeper regions than those probed by spectroscopic observations (e.g. Brouwers & Ormel 2020).

4.3 Time evolution

Atmospheric loss processes are expected to affect a significant portion of sub-Neptunes, so the interplay between such processes and any evolution of the atmospheric composition should be considered carefully. We find that compositional gradients at depth tend to shrink the overall planet radius, if everything else is kept constant. A smaller radius tends to inhibit atmospheric loss, as the atmosphere must be removed from deeper within the gravitational potential well. However, atmospheric accretion dictates that planets’ initial radiative–convective boundaries are close to the Bondi radius (e.g. Ginzburg et al. 2016) and deep non-convective layers slow thermal contraction over time (Misener & Schlichting 2022). Therefore, whether the overall effect of the compositional gradients examined here inhibits or furthers atmospheric loss remains to be determined, highlighting the need for fully self-consistent accretion and loss models.

Another possibility is that hydrogen could be preferentially lost, increasing the molecular weight of the atmosphere (e.g. Malsky et al. 2023). The two major loss processes thought to shape exoplanet demographics, core-powered mass-loss and photoevaporation, are hydrodynamic. Their strong winds are sufficient to drag along heavier species such as those we find could be present in this study (e.g. Misener & Schlichting 2021). However, due to the gradient in molecular weight, a wind escaping from the top of the atmosphere is relatively less enriched in outgassed species than the mean overall atmosphere, which could produce effective fractionation. Additionally, as hydrogen is depleted, the chemistry of the atmosphere could change, with the magma–hydrogen interaction producing more oxidized species (e.g. Zilinskas et al. 2023). A model that simultaneously evolves the thermal state, atmospheric composition, and atmospheric mass could investigate these potential interplays further.

5 CONCLUSIONS

We analyse the chemical equilibrium between hydrogen gas and species vapourized from the surface of a silicate magma ocean, a condition likely to be present in the depths of young sub-Neptune planets. We find that SiO and O₂ react with hydrogen gas to produce significant amounts of silane (SiH₄) and water vapour (H₂O). The resulting depletion of SiO₂ draws more magma ocean products into the atmosphere, greatly increasing the atmosphere’s silicon content compared to a model that does not account for the reducing conditions of the atmosphere. The amounts and proportions of these products likely vary depending on the chemical state of the interior. This implies that the atmospheric compositions of planets with magma oceans are a window into their interior composition, a promising prospect since such atmospheric abundances are observable with current and future telescopes. The chemical products of magma–atmosphere interaction in turn alter the atmospheric structure, inhibiting convection in the interior to a larger extent than previously found. These results imply that the presence of a magma ocean must be considered in order to understand the chemical abundances and overall atmospheric mass fractions of sub-Neptune planets.

ACKNOWLEDGEMENTS

We thank James Rogers, Namrah Habib, and the anonymous reviewer for useful comments that improved the manuscript. In this work, we use the numpy (Harris et al. 2020), matplotlib (Hunter 2007), and scipy (Virtanen et al. 2020) packages. This research has been supported in part by the Alfred P. Sloan Foundation under grant G202114194 as part of the AEThER collaboration and by NASA under grant number 80NSSC21K0392 issued through the Exoplanet Research Program.

DATA AVAILABILITY

Data are available upon request.

APPENDIX A: EQUILIBRIUM CONSTANT VALUES

In Figs A1,  A2, and  A3, we plot the equilibrium constants we use in this work as functions of temperature. As described in Section 2.1, all values are calculated from the Gibbs free energies in NIST.

Figure A1.

The equilibrium constant describing reaction (R1), as defined in equation (1), as a function of temperature.

Figure A2.

The equilibrium constant describing reaction (R2), as defined in equation (2), as a function of temperature.

Figure A3.

The equilibrium constant describing reaction (R3), as defined in equation (3), as a function of temperature.

APPENDIX B: DERIVATION OF PARTIAL PRESSURES

As stated in Section 2.1, equations (1), (2), (3), and (6), along with the total pressure P and temperature T, are sufficient to solve for the partial pressures of all components of the atmosphere at a given level, specifically |$P_\mathrm{H_2}$|⁠, |$P_\mathrm{H_2O}$|⁠, P_SiO, |$P_\mathrm{SiH_4}$|⁠, and |$P_\mathrm{O_2}$|⁠. In practice, we use the following method, though many alternative derivations are possible in principle.

We begin by assuming an oxygen partial pressure, |$P_\mathrm{O_2}$|⁠. We choose to guess this value because it varies by many orders of magnitude over a typical atmospheric profile, the most of any of the unknowns in our problem. Therefore, small changes in the other inferred partial pressures lead to large changes in the inferred |$P_\mathrm{O_2}$|⁠, making it more amenable to solve for via numerical techniques than the other values, which vary slowly and therefore lead to failures to converge. Given the assumed |$P_\mathrm{O_2}$| and the temperature, which yields values for all three equilibrium constants, it is easy to use equation (1) to solve for the partial pressure of SiO:

$$\begin{eqnarray} P_\mathrm{SiO} = K_\mathrm{eq,~R1}/P_\mathrm{O_2}^{1/2}. \end{eqnarray}$$

(B1)

Similarly, we can invert equation (2) to solve for the partial pressure of water as a function of the partial pressure of H₂, which is as yet unknown:

$$\begin{eqnarray} P_\mathrm{H_2O} = K_\mathrm{eq,~R2} P_\mathrm{H_2} P_\mathrm{O_2}^{1/2}. \end{eqnarray}$$

(B2)

Inserting equation (B2) into equation (3), we can solve for |$P_\mathrm{SiH_4}$|⁠:

$$\begin{eqnarray} P_\mathrm{SiH_4} = \frac{K_\mathrm{eq,~R3} P_\mathrm{SiO} P_\mathrm{H_2}^2}{K_\mathrm{eq,~R2} P_\mathrm{O_2}^{1/2}}, \end{eqnarray}$$

(B3)

which also depends on the unknown |$P_\mathrm{H_2}$|⁠. We can leverage the definition of P as the sum of all partial pressures to express |$P_\mathrm{H_2}$| in terms of the other species:

$$\begin{eqnarray} P_\mathrm{H_2} = P-P_\mathrm{H_2O}-P_\mathrm{SiO}-P_\mathrm{O_2}-P_\mathrm{SiH_4}. \end{eqnarray}$$

(B4)

Inserting equation (B2) into this equation and collecting terms of |$P_\mathrm{H_2}$| allows us to solve for |$P_\mathrm{H_2}$| as a function of known values and |$P_\mathrm{SiH_4}$|⁠:

$$\begin{eqnarray} P_\mathrm{H_2} = \frac{P-P_\mathrm{SiO}-P_\mathrm{O_2}-P_\mathrm{SiH_4}}{1 + K_\mathrm{eq,~R2} P_\mathrm{O_2}^{1/2}}. \end{eqnarray}$$

(B5)

Such an expression allows us to substitute |$P_\mathrm{H_2}$| in equation (B3), eliminating all other unknown variables. This yields a simple quadratic equation:

$$\begin{eqnarray} P_\mathrm{SiH_4} = \alpha (\beta -P_\mathrm{SiH_4})^2, \end{eqnarray}$$

(B6)

where

$$\begin{eqnarray} \alpha \equiv \frac{K_\mathrm{eq,~R3} P_\mathrm{SiO}}{K_\mathrm{eq,~R2} P_\mathrm{O_2}^{1/2} \big (1 + K_\mathrm{eq,~R2} P_\mathrm{O_2}^{1/2}\big)^2} \end{eqnarray}$$

(B7)

and

$$\begin{eqnarray} \beta \equiv P-P_\mathrm{SiO}-P_\mathrm{O_2} \end{eqnarray}$$

(B8)

are functions of known quantities defined for simplicity. Equation (B6) can be solved using the quadratic formula:

$$\begin{eqnarray} P_\mathrm{SiH_4} = \frac{2 \alpha \beta + 1 \pm \sqrt{4 \alpha \beta +1}}{2 \alpha }. \end{eqnarray}$$

(B9)

Mathematically, this yields two solutions; however, in all cases only the negative branch is physically reasonable (i.e. yields P_i > 0 for all species). Once |$P_\mathrm{SiH_4}$| is known, |$P_\mathrm{H_2}$| can be solved via equation (B5), and |$P_\mathrm{H_2O}$| can be solved via equation (B2). This yields a full set of partial pressures, which correctly sum to the total pressure, and conform to chemical equilibrium at a given temperature, for a guessed |$P_\mathrm{O_2}$|⁠. The only constraint not yet used is the number of atoms constraint (equation 6). We now calculate the ‘atomic partial pressures’ using equations (4) and (5), and calculate the difference ∑P_O − ∑P_Si. This difference will vary as a function of the input |$P_\mathrm{O_2}$|⁠, with one unique solution where the difference is zero. We solve for this |$P_\mathrm{O_2}$| that balances the reactions correctly numerically for each layer, using the fsolve function of the scipy.optimize package (Virtanen et al. 2020). In practice, since we solve the atmospheric structure layer by layer in small steps, the previous value of |$P_\mathrm{O_2}$| provides a good starting guess for the value in the next layer deeper, which we use to increase computational speed compared to a fully naive guess.

APPENDIX C: DERIVATION OF MULTISPECIES CONVECTION CRITERION

As discussed in Section 2.2, in this work we derive a multispecies convection criterion, which turns out to be the sum of the individual convection criteria of each species, ∑_ic_i ≥ 1, where c_i is given by equation (12). It is non-trivial that the overall convection criterion is the sum of those of each species, so we demonstrate that here. For demonstration purposes, we consider the case of two condensables, but the argument applies to n condensables equally well. Let us define a three component atmosphere, composed of a dry component with mass mixing ratio q_d and molecular weight μ_d, and two condensable species with mass mixing ratios q₁ and q₂ and molecular weights μ₁ and μ₂, respectively. By definition, q_d + q₁ + q₂ = 1. Meanwhile, the overall mean molecular weight μ is given by

$$\begin{eqnarray} \mu = \frac{\mu _\mathrm{d} \mu _1 \mu _2}{\mu _1 \mu _2 + (\mu _\mathrm{d} \mu _2-\mu _1 \mu _2) q_1 + (\mu _\mathrm{d} \mu _1-\mu _1 \mu _2) q_2}. \end{eqnarray}$$

(C1)

For convection to be inhibited, the density gradient in the environment must be steeper than that of a parcel moved adiabatically (e.g. Leconte et al. 2017):

$$\begin{eqnarray} \bigg (\frac {\partial \ln T}{\partial \ln P}-\frac {\partial \ln \mu }{\partial \ln P}\bigg)_\mathrm{env} \gt \bigg (\frac {\partial \ln T}{\partial \ln P}-\frac {\partial \ln \mu }{\partial \ln P}\bigg)_\mathrm{ad}. \end{eqnarray}$$

(C2)

Therefore, to assess whether convection operates in a regime with multiple species changing in abundance, we must compute the change in molecular weight with pressure |$\frac *{\partial \ln \mu }{\partial \ln P}$|⁠. The molecular weight changes with pressure for two reasons: due to changes in the abundance of species 1, or due to changes in the abundance of species 2. Quantitatively, we can express this as a sum of partial derivatives with the mixing ratio of the other species held constant:

$$\begin{eqnarray} \frac {\partial \ln \mu }{\ln P} = \frac {\partial \ln \mu }{\partial \ln P}\bigg |_{q_2}(q_1) + \frac {\partial \ln \mu }{\partial \ln P}\bigg |_{q_1}(q_2). \end{eqnarray}$$

(C3)

Here, the first term is the change in molecular weight due to changing q₁, with q₂ held fixed, and the second term is the change in molecular weight due to changing q₂, with q₁ held fixed.

Examining the first term in equation (C3), it can be expanded into the product of two partial derivatives:

$$\begin{eqnarray} \frac {\ln \mu }{\partial \ln P}\bigg |_{q_2} = \frac {\partial \ln \mu }{\partial \ln q_1}\bigg |_{q_2} \frac {\partial \ln q_1}{\partial \ln P}\bigg |_{q_2}. \end{eqnarray}$$

(C4)

The first term in the product can be computed using the definition of μ in equation (C1) and simplifies to

$$\begin{eqnarray} \frac {\partial \ln \mu }{\ln q_1}\bigg |_{q_2} = \mu q_1 \bigg (\frac{1}{\mu _\mathrm{d}}-\frac{1}{\mu _1} \bigg). \end{eqnarray}$$

(C5)

Meanwhile, the second term in equation (C4) can be expressed as the sum of two terms: the change in condensable mass mixing ratio as the total pressure changes at fixed temperature, and the change in condensable mass mixing ratio as the temperature changes due to the T–P relation:

$$\begin{eqnarray} \frac {\partial \ln q_1}{\partial \ln P}\bigg |_{q_2} = \frac {\partial \ln q_1}{\ln P}\bigg |_{T, q_2} + \frac {\partial \ln q_1}{\partial \ln T}\bigg |_{P, q_2} \frac {\ln T}{\ln P}. \end{eqnarray}$$

(C6)

The first term in this equation can be solved using the definition of the mass mixing ratio, q₁ = μ₁P₁/(μP), and reduces to

$$\begin{eqnarray} \frac {\partial \ln q_1}{\partial \ln P}\bigg |_{T, q_2} = - \frac{\mu _\mathrm{d}}{\mu }. \end{eqnarray}$$

(C7)

We note that this term is negative, because increasing the total pressure without changing the temperature leaves the partial pressure of the condensable unchanged, therefore decreasing the mass mixing ratio.

The derivative of q₁ with respect to T, which appears in the second term of equation (C6), can also be solved by inserting the definition of the mass mixing ratio, but it simplifies considerably less and introduces cross-terms dependent on q₂:

$$\begin{eqnarray} \frac {\partial \ln q_1}{\partial \ln T}\bigg |_{P, q_2} = \bigg [\frac{\mu _\mathrm{d}}{\mu } + q_2 \bigg (1-\frac{\mu _\mathrm{d}}{\mu _2}\bigg)\bigg ] \frac {\partial \ln P_1}{\partial \ln T}-q_2 \bigg (1-\frac{\mu _\mathrm{d}}{\mu _2}\bigg) \frac {\partial \ln P_2}{\partial \ln T}. \\ \end{eqnarray}$$

(C8)

Equations (C7) and (C8) can then be inserted into equation (C6), which can be substituted along with equation (C5) into equation (C4). After some algebraic manipulation, this substitution leads to an expression for the overall gradient as a function of q₁ with q₂ fixed:

$$\begin{eqnarray} \frac {\partial \ln \mu }{\partial \ln P}\bigg |_{q_2} &=& q_1 \bigg (1-\frac{\mu _\mathrm{d}}{\mu _1}\bigg) \bigg [\frac {\partial \ln P_1} {\partial \ln T}\frac {\partial \ln T}{\partial \ln P}-1\bigg ] \\ &+& \frac{\mu }{\mu _\mathrm{d}} q_1 q_2 (1-\frac{\mu _\mathrm{d}}{\mu _1}) (1-\frac{\mu _\mathrm{d}}{\mu _2}) \frac {\partial \ln T}{\partial \ln P}. \end{eqnarray}$$

(C9)

As q₁ and q₂ are completely symmetrical in these equations, the expression for |$\frac *{\partial \ln \mu }{\partial \ln P}|_{q_1}$| can be obtained by swapping the 1 and 2 subscripts in equation (C9). Upon adding these equations together per equation (C3), it is immediately clear that the second terms cancel, leaving

$$\begin{eqnarray} \frac {\partial \ln \mu }{\partial \ln P} &=& q_1 \bigg (1-\frac{\mu _\mathrm{d}}{\mu _1}\bigg)\bigg [\frac {\partial \ln P_1}{\partial \ln T}\frac {\partial \ln T}{\partial \ln P}-1\bigg ] \\ &+& q_2 \bigg (1-\frac{\mu _\mathrm{d}}{\mu _2} \bigg) \bigg [\frac {\partial \ln P_2}{\partial \ln T} \frac {\partial \ln T}{\partial \ln P}-1\bigg ] \\ &=& \alpha _1 \gamma _1 \frac {\partial \ln T}{\partial \ln P}-\alpha _1 (1-\varpi _1 q_1-\varpi _2 q_2) \\ &+& \alpha _2 \gamma _2 \frac {\partial \ln T}{\partial \ln P}-\alpha _2 (1-\varpi _1 q_1-\varpi _2 q_2) \end{eqnarray}$$

(C10)

using the notation of Leconte et al. (2017), where

$$\begin{eqnarray} \alpha _i \equiv \mu q_i \bigg (\frac{1}{\mu _\mathrm{d}}-\frac{1}{\mu _i}\bigg), \end{eqnarray}$$

(C11)

$$\begin{eqnarray} \gamma _i \equiv \frac{\mu _\mathrm{d}}{\mu }\frac {\partial \ln P_i}{\partial \ln T}, \end{eqnarray}$$

(C12)

and

$$\begin{eqnarray} \varpi _i \equiv 1-\frac{\mu _\mathrm{d}}{\mu _i}. \end{eqnarray}$$

(C13)

Finally, we can insert equation (C10) into equation (C2) to obtain the convection criterion. We note that the second and fourth terms of equation (C10) will be the same on either side of the inequality in equation (C2) and therefore cancel, eliminating all further cross terms. This leaves the inequality as

$$\begin{eqnarray} \bigg [\bigg (\frac {\ln T}{\partial \ln P}\bigg)_\mathrm{env}-\bigg (\frac {\partial \ln T}{\ln P}\bigg)_\mathrm{ad}\bigg ] (1-\alpha _1 \gamma _1-\alpha _2 \gamma _2) \gt 0. \end{eqnarray}$$

(C14)

Since the environmental temperature gradient will not be less than the adiabatic gradient, this inequality implies that

$$\begin{eqnarray} 1 \gt \alpha _1 \gamma _1 + \alpha _2 \gamma _2, \end{eqnarray}$$

(C15)

where α_iγ_i = χ_i as defined in equation (12): the inhibition criteria solely add.