ANALYTICAL MODELS OF EXOPLANETARY ATMOSPHERES. III. GASEOUS C–H–O–N CHEMISTRY WITH NINE MOLECULES

Atmospheric chemistry is an indispensible ingredient in the study of exoplanetary atmospheres, as it teaches the practitioner how and when to be surprised. For example, if hydrogen-dominated atmospheres are in chemical equilibrium, then we expect the dominant carbon carriers to be methane and carbon monoxide at low and high temperatures, respectively. To date, the contributions to the atmospheric chemistry literature have mostly taken the form of numerical calculations using equilibrium chemistry and chemical kinetics codes (e.g., Burrows & Sharp 1999; Lodders & Fegley 2002; Line et al. 2011; Moses et al. 2011, 2013a, 2013b; Visscher & Moses 2011; Kopparapu et al. 2012; Madhusudhan 2012; Agúndez et al. 2014; Hu & Seager 2014; Hu et al. 2015; Venot et al. 2015). A complementary approach, which is standard in the astrophysical literature, is to develop analytical models (Heng & Lyons 2016; Heng et al. 2016). The current study is the third in a series of papers devoted to constructing analytical models for exoplanetary atmospheres to aid in the development of intuition, following Heng & Workman (2014; for shallow-water fluid dynamics) and Heng et al. (2014; for two-stream radiative transfer).

Specifically, Heng et al. (2016) and Heng & Lyons (2016) have previously derived solutions for the relative abundances of molecules for purely gaseous chemistry and in C–H–O (carbon–hydrogen–oxygen) systems. Here, we present novel, generalized analytical solutions for purely gaseous C–H–O–N (carbon–hydrogen–oxygen-nitrogen) systems with six and nine molecules. In Section 2, we concisely describe the theoretical setup, which we use to consider six and nine molecules in Sections 3 and 4, respectively. We present our results are in Section 5 and conclude in Section 6. A Python script that implements our analytical formulae may be found at http://github.com/exoclime/VULCAN.

2.1. Net Chemical Reactions

Our C–H–O–N network contains six net reactions,

$$\begin{eqnarray}{\mathrm{CH}}_{4}+{{\rm{H}}}_{2}{\rm{O}} & \leftrightarrows & \mathrm{CO}+3{{\rm{H}}}_{2},\\ {\mathrm{CO}}_{2}+{{\rm{H}}}_{2} & \leftrightarrows & \mathrm{CO}+{{\rm{H}}}_{2}{\rm{O}},\\ 2{\mathrm{CH}}_{4} & \leftrightarrows & {{\rm{C}}}_{2}{{\rm{H}}}_{2}+3{{\rm{H}}}_{2},\\ {{\rm{C}}}_{2}{{\rm{H}}}_{4} & \leftrightarrows & {{\rm{C}}}_{2}{{\rm{H}}}_{2}+{{\rm{H}}}_{2},\\ 2{\mathrm{NH}}_{3} & \leftrightarrows & {{\rm{N}}}_{2}+3{{\rm{H}}}_{2},\\ {\mathrm{NH}}_{3}+{\mathrm{CH}}_{4} & \leftrightarrows & \mathrm{HCN}+3{{\rm{H}}}_{2}.\end{eqnarray} \tag{ 1 }$$

The last three reactions, involving nitrogen, were informed by Burrows & Sharp (1999), Lodders & Fegley (2002), and Moses et al. (2011). Excluding molecular hydrogen, there are nine molecules in total.

2.2. Normalized Equilibrium Constants

For each net reaction, there is a corresponding equilibrium constant. In a departure from Heng et al. (2016) and Heng & Lyons (2016), we write the normalized equilibrium constants without primes as superscripts. Obeying the order in Equation (1), they are

$$\begin{eqnarray}K & = & \displaystyle \frac{{\tilde{n}}_{\mathrm{CO}}}{{\tilde{n}}_{{\mathrm{CH}}_{4}}{\tilde{n}}_{{{\rm{H}}}_{2}{\rm{O}}}}={\left(\displaystyle \frac{{P}_{0}}{P}\right)}^{2}\exp \left(-\displaystyle \frac{{\rm{\Delta }}{\tilde{G}}_{\mathrm{0,1}}}{{{ \mathcal R }}_{\mathrm{univ}}T}\right),\\ {K}_{2} & = & \displaystyle \frac{{\tilde{n}}_{\mathrm{CO}}{\tilde{n}}_{{{\rm{H}}}_{2}{\rm{O}}}}{{\tilde{n}}_{{\mathrm{CO}}_{2}}}=\exp \left(-\displaystyle \frac{{\rm{\Delta }}{\tilde{G}}_{\mathrm{0,2}}}{{{ \mathcal R }}_{\mathrm{univ}}T}\right),\\ {K}_{3} & = & \displaystyle \frac{{\tilde{n}}_{{{\rm{C}}}_{2}{{\rm{H}}}_{2}}}{{\tilde{n}}_{{\mathrm{CH}}_{4}}^{2}}={\left(\displaystyle \frac{{P}_{0}}{P}\right)}^{2}\exp \left(-\displaystyle \frac{{\rm{\Delta }}{\tilde{G}}_{\mathrm{0,3}}}{{{ \mathcal R }}_{\mathrm{univ}}T}\right),\\ {K}_{4} & = & \displaystyle \frac{{\tilde{n}}_{{{\rm{C}}}_{2}{{\rm{H}}}_{2}}}{{\tilde{n}}_{{{\rm{C}}}_{2}{{\rm{H}}}_{4}}}=\displaystyle \frac{{P}_{0}}{P}\exp \left(-\displaystyle \frac{{\rm{\Delta }}{\tilde{G}}_{\mathrm{0,4}}}{{{ \mathcal R }}_{\mathrm{univ}}T}\right),\\ {K}_{5} & = & \displaystyle \frac{{\tilde{n}}_{{{\rm{N}}}_{2}}}{{\tilde{n}}_{{\mathrm{NH}}_{3}}^{2}}={\left(\displaystyle \frac{{P}_{0}}{P}\right)}^{2}\exp \left(-\displaystyle \frac{{\rm{\Delta }}{\tilde{G}}_{\mathrm{0,5}}}{{{ \mathcal R }}_{\mathrm{univ}}T}\right),\\ {K}_{6} & = & \displaystyle \frac{{\tilde{n}}_{\mathrm{HCN}}}{{\tilde{n}}_{{\mathrm{NH}}_{3}{\tilde{n}}_{{\mathrm{CH}}_{4}}}}={\left(\displaystyle \frac{{P}_{0}}{P}\right)}^{2}\exp \left(-\displaystyle \frac{{\rm{\Delta }}{\tilde{G}}_{\mathrm{0,6}}}{{{ \mathcal R }}_{\mathrm{univ}}T}\right),\end{eqnarray} \tag{ 2 }$$

where ${P}_{0}=1\,\mathrm{bar}$ is the reference pressure and ${{ \mathcal R }}_{\mathrm{univ}}\ =8.3144621$ J K⁻¹ mol⁻¹ is the universal gas constant. For a molecule X, we have defined ${\tilde{n}}_{{\rm{X}}}\equiv {n}_{{\rm{X}}}/{n}_{{{\rm{H}}}_{2}}$, where ${n}_{{\rm{X}}}$ denotes the number density. We call ${\tilde{n}}_{{\rm{X}}}$ the "mixing ratios."

Appendix A lists the Gibbs free energies of formation for C₂H₄, N₂, NH₃, and HCN. Appendix A also lists the Gibbs free energies associated with the last three net reactions (denoted by ${\rm{\Delta }}{\tilde{G}}_{0,i}$ for the ith reaction). The Gibbs free energies for the other molecules and the first three net reactions have previously been stated in Tables 1 and 2 of Heng & Lyons (2016), respectively.

2.3. Particle Conservation Equations

By counting the number of atoms sequestered in each molecule, we obtain

$$\begin{eqnarray}&&{n}_{{\mathrm{CH}}_{4}}+{n}_{\mathrm{CO}}+{n}_{{\mathrm{CO}}_{2}}+2{n}_{{{\rm{C}}}_{2}{{\rm{H}}}_{2}}+2{n}_{{{\rm{C}}}_{2}{{\rm{H}}}_{4}}\\ &&\quad +{n}_{\mathrm{HCN}}={n}_{{\rm{C}}},\\ &&4{n}_{{\mathrm{CH}}_{4}}+2{n}_{{{\rm{H}}}_{2}{\rm{O}}}+2{n}_{{{\rm{H}}}_{2}}+2{n}_{{{\rm{C}}}_{2}{{\rm{H}}}_{2}}+4{n}_{{{\rm{C}}}_{2}{{\rm{H}}}_{4}}\\ &&\quad +{n}_{\mathrm{HCN}}+3{n}_{{\mathrm{NH}}_{3}}={n}_{{\rm{H}}},\\ &&{n}_{{{\rm{H}}}_{2}{\rm{O}}}+{n}_{\mathrm{CO}}+2{n}_{{\mathrm{CO}}_{2}}={n}_{{\rm{O}}},\\ &&2{n}_{{{\rm{N}}}_{2}}+{n}_{{\mathrm{NH}}_{3}}+{n}_{\mathrm{HCN}}={n}_{{\rm{N}}}.\end{eqnarray} \tag{ 3 }$$

The number densities of atomic carbon, hydrogen, oxygen, and nitrogen are given by ${n}_{{\rm{C}}}$, ${n}_{{\rm{H}}}$, ${n}_{{\rm{O}}}$ and ${n}_{{\rm{N}}}$, respectively.

To develop our current analytical method, we first ignore CO₂, C₂H₂, and C₂H₄ in our analysis. Heng & Lyons (2016) have shown that CO₂ is subdominant compared to CO and H₂O, unless the metallicity is orders of magnitude higher than solar abundance.

If we make the simplification that $2{n}_{{{\rm{H}}}_{2}}={n}_{{\rm{H}}}$ (hydrogen-dominated atmospheres) and render the number densities dimensionless, then we end up with

$$\begin{eqnarray}&&{\tilde{n}}_{{\mathrm{CH}}_{4}}+{\tilde{n}}_{\mathrm{CO}}+{\tilde{n}}_{\mathrm{HCN}}=2{\tilde{n}}_{{\rm{C}}},\\ &&{\tilde{n}}_{{{\rm{H}}}_{2}{\rm{O}}}+{\tilde{n}}_{\mathrm{CO}}=2{\tilde{n}}_{{\rm{O}}},\\ &&2{\tilde{n}}_{{{\rm{N}}}_{2}}+{\tilde{n}}_{{\mathrm{NH}}_{3}}+{\tilde{n}}_{\mathrm{HCN}}=2{\tilde{n}}_{{\rm{N}}}.\end{eqnarray} \tag{ 4 }$$

We define the elemental abundances as ${\tilde{n}}_{{\rm{C}}}\equiv {n}_{{\rm{C}}}/{n}_{{\rm{H}}}$, ${\tilde{n}}_{{\rm{O}}}\equiv {n}_{{\rm{O}}}/{n}_{{\rm{H}}}$, and ${\tilde{n}}_{{\rm{N}}}\equiv {n}_{{\rm{N}}}/{n}_{{\rm{H}}}$. The goal is to decouple this system of nonlinear equations such that one obtains a polynomial equation describing only one of the molecules. This requires that we first rewrite some of the mixing ratios in terms of only ${\tilde{n}}_{\mathrm{CO}}$ and ${\tilde{n}}_{{{\rm{H}}}_{2}{\rm{O}}}$,

$$\begin{eqnarray}{\tilde{n}}_{{\mathrm{CH}}_{4}} & = & \displaystyle \frac{{\tilde{n}}_{\mathrm{CO}}}{K{\tilde{n}}_{{{\rm{H}}}_{2}{\rm{O}}}},\\ {\tilde{n}}_{{\mathrm{NH}}_{3}} & = & \displaystyle \frac{K{\tilde{n}}_{\mathrm{CO}}^{2}-{D}_{1}{\tilde{n}}_{\mathrm{CO}}+4K{\tilde{n}}_{{\rm{C}}}{\tilde{n}}_{{\rm{O}}}}{{K}_{6}{\tilde{n}}_{\mathrm{CO}}},\\ {\tilde{n}}_{\mathrm{HCN}} & = & \displaystyle \frac{{K}_{6}{\tilde{n}}_{{\mathrm{NH}}_{3}}{\tilde{n}}_{\mathrm{CO}}}{K{\tilde{n}}_{{{\rm{H}}}_{2}{\rm{O}}}},\\ {\tilde{n}}_{{{\rm{N}}}_{2}} & = & {K}_{5}{\tilde{n}}_{{\mathrm{NH}}_{3}}^{2},\end{eqnarray} \tag{ 5 }$$

where we have defined

$$\begin{eqnarray}&&{D}_{1}\equiv 1+2K({\tilde{n}}_{{\rm{C}}}+{\tilde{n}}_{{\rm{O}}}).\end{eqnarray} \tag{ 6 }$$

We then use

$$\begin{eqnarray}&&{\tilde{n}}_{{{\rm{H}}}_{2}{\rm{O}}}=2{\tilde{n}}_{{\rm{O}}}-{\tilde{n}}_{\mathrm{CO}}\end{eqnarray} \tag{ 7 }$$

to eliminate the mixing ratio of water.

By substituting these expressions into the equation involving molecular nitrogen, we obtain a quintic equation for the mixing ratio of CO,

$$\begin{eqnarray}&&\displaystyle \sum _{i=0}^{5}{A}_{i}{\tilde{n}}_{\mathrm{CO}}^{i}=0.\end{eqnarray} \tag{ 8 }$$

The coefficients of this quintic equation are

$$\begin{eqnarray}{A}_{0} & = & 256{K}^{3}{K}_{5}{\tilde{n}}_{{\rm{O}}}^{3}{\tilde{n}}_{{\rm{C}}}^{2},\\ {A}_{1} & = & 32{K}^{2}{\tilde{n}}_{{\rm{O}}}^{2}{\tilde{n}}_{{\rm{C}}}[{K}_{6}-4{K}_{5}({D}_{1}+K{\tilde{n}}_{{\rm{C}}})],\\ {A}_{2} & = & 16{{KK}}_{5}{\tilde{n}}_{{\rm{O}}}(8{K}^{2}{\tilde{n}}_{{\rm{O}}}{\tilde{n}}_{{\rm{C}}}+{D}_{1}^{2}+4{{KD}}_{1}{\tilde{n}}_{{\rm{C}}})\\ & & +\,8{{KK}}_{6}{\tilde{n}}_{{\rm{O}}}[2{K}_{6}({\tilde{n}}_{{\rm{C}}}-{\tilde{n}}_{{\rm{N}}})-2K{\tilde{n}}_{{\rm{C}}}-{D}_{1}],\\ {A}_{3} & = & -8{{KK}}_{5}(4{{KD}}_{1}{\tilde{n}}_{{\rm{O}}}+8{K}^{2}{\tilde{n}}_{{\rm{O}}}{\tilde{n}}_{{\rm{C}}}+{D}_{1}^{2})\\ & & +4{{KK}}_{6}(2K{\tilde{n}}_{{\rm{O}}}+{D}_{1})+4{K}_{6}^{2}(2K{\tilde{n}}_{{\rm{N}}}-{D}_{1}),\\ {A}_{4} & = & 16{K}^{2}{K}_{5}(K{\tilde{n}}_{{\rm{O}}}+{D}_{1})+4{{KK}}_{6}({K}_{6}-K),\\ {A}_{5} & = & -\,8{K}^{3}{K}_{5}.\end{eqnarray} \tag{ 9 }$$

This derivation demonstrates that it is possible to decouple a C–H–O–N system and obtain an equation in terms of only ${\tilde{n}}_{\mathrm{CO}}$.

We now add CO₂, C₂H₂, and C₂H₄ back into the analysis and use the method developed in the previous section. The particle conservation equations become

$$\begin{eqnarray}&&{\tilde{n}}_{{\mathrm{CH}}_{4}}+{\tilde{n}}_{\mathrm{CO}}+{\tilde{n}}_{{\mathrm{CO}}_{2}}+{\tilde{n}}_{\mathrm{HCN}}\\ &&\quad +2{\tilde{n}}_{{{\rm{C}}}_{2}{{\rm{H}}}_{2}}+2{\tilde{n}}_{{{\rm{C}}}_{2}{{\rm{H}}}_{4}}=2{\tilde{n}}_{{\rm{C}}},\\ &&{\tilde{n}}_{{{\rm{H}}}_{2}{\rm{O}}}+{\tilde{n}}_{\mathrm{CO}}+2{\tilde{n}}_{{\mathrm{CO}}_{2}}=2{\tilde{n}}_{{\rm{O}}},\\ &&2{\tilde{n}}_{{{\rm{N}}}_{2}}+{\tilde{n}}_{{\mathrm{NH}}_{3}}+{\tilde{n}}_{\mathrm{HCN}}=2{\tilde{n}}_{{\rm{N}}}.\end{eqnarray} \tag{ 10 }$$

The various mixing ratios are now described by

$$\begin{eqnarray}{\tilde{n}}_{{{\rm{H}}}_{2}{\rm{O}}} & = & \displaystyle \frac{{D}_{5}}{{D}_{4}},\\ {\tilde{n}}_{{\mathrm{CO}}_{2}} & = & \displaystyle \frac{{\tilde{n}}_{\mathrm{CO}}{\tilde{n}}_{{{\rm{H}}}_{2}{\rm{O}}}}{{K}_{2}},\\ {\tilde{n}}_{{\mathrm{NH}}_{3}} & = & \displaystyle \frac{{D}_{2}}{{{KK}}_{6}{D}_{4}^{2}{D}_{5}{\tilde{n}}_{\mathrm{CO}}},\\ {\tilde{n}}_{\mathrm{HCN}} & = & \displaystyle \frac{{K}_{6}{D}_{4}{\tilde{n}}_{{\mathrm{NH}}_{3}}{\tilde{n}}_{\mathrm{CO}}}{{{KD}}_{5}}.\end{eqnarray} \tag{ 11 }$$

The expressions for ${\tilde{n}}_{{\mathrm{CH}}_{4}}$ and ${\tilde{n}}_{{{\rm{N}}}_{2}}$ are the same as those given in Equation (5). To make the algebra tractable, we have defined

$$\begin{eqnarray}{D}_{2} & \equiv & -\,{{KD}}_{4}^{2}{D}_{5}{\tilde{n}}_{\mathrm{CO}}+{K}^{2}{D}_{4}{D}_{5}^{2}(2{\tilde{n}}_{{\rm{C}}}-{\tilde{n}}_{\mathrm{CO}})\\ & & -2{K}_{3}{D}_{3}{D}_{4}^{3}{\tilde{n}}_{\mathrm{CO}}^{2}-\displaystyle \frac{{K}^{2}{D}_{5}^{3}{\tilde{n}}_{\mathrm{CO}}}{{K}_{2}},\\ {D}_{3} & \equiv & 1+\displaystyle \frac{1}{{K}_{4}},\\ {D}_{4} & \equiv & 1+\displaystyle \frac{2{\tilde{n}}_{\mathrm{CO}}}{{K}_{2}},\\ {D}_{5} & \equiv & 2{\tilde{n}}_{{\rm{O}}}-{\tilde{n}}_{\mathrm{CO}}.\end{eqnarray} \tag{ 12 }$$

If we write D₂ as

$$\begin{eqnarray}&&{D}_{2}=\displaystyle \sum _{i=0}^{5}{F}_{i}{\tilde{n}}_{\mathrm{CO}}^{i},\end{eqnarray} \tag{ 13 }$$

then one may show that the coefficients are

$$\begin{eqnarray}{F}_{0} & = & 8{K}^{2}{\tilde{n}}_{{\rm{O}}}^{2}{\tilde{n}}_{{\rm{C}}},\\ {F}_{1} & = & 2K{\tilde{n}}_{{\rm{O}}}\left\{-1+2K\left[2{\tilde{n}}_{{\rm{C}}}\left(\displaystyle \frac{2{\tilde{n}}_{{\rm{O}}}}{{K}_{2}}-1\right)-{\tilde{n}}_{{\rm{O}}}\right]\right.\\ & & \left.-\displaystyle \frac{4K{\tilde{n}}_{{\rm{O}}}^{2}}{{K}_{2}}\right\},\\ {F}_{2} & = & K\left(1-\displaystyle \frac{8{\tilde{n}}_{{\rm{O}}}}{{K}_{2}}\right)-2{K}_{3}{D}_{3}\\ & & +2{K}^{2}\left[{\tilde{n}}_{{\rm{C}}}\left(1-\displaystyle \frac{8{\tilde{n}}_{{\rm{O}}}}{{K}_{2}}\right)+2{\tilde{n}}_{{\rm{O}}}\left(1+\displaystyle \frac{{\tilde{n}}_{{\rm{O}}}}{{K}_{2}}\right)\right],\\ {F}_{3} & = & \displaystyle \frac{4K}{{K}_{2}}\left(1-\displaystyle \frac{2{\tilde{n}}_{{\rm{O}}}}{{K}_{2}}\right)-\displaystyle \frac{12{K}_{3}{D}_{3}}{{K}_{2}}\\ & & +{K}^{2}\left[\displaystyle \frac{2(2{\tilde{n}}_{{\rm{C}}}+{\tilde{n}}_{{\rm{O}}})}{{K}_{2}}-1\right],\\ {F}_{4} & = & \displaystyle \frac{K}{{K}_{2}}\left(\displaystyle \frac{4}{{K}_{2}}-K\right)-\displaystyle \frac{24{K}_{3}{D}_{3}}{{K}_{2}^{2}},\\ {F}_{5} & = & -\,\displaystyle \frac{16{K}_{3}{D}_{3}}{{K}_{2}^{3}}.\end{eqnarray} \tag{ 14 }$$

Again, the goal is to obtain a single equation for ${\tilde{n}}_{\mathrm{CO}}$ by substituting all of these expressions into the equation involving molecular nitrogen,

$$\begin{eqnarray}&&2{K}_{5}{D}_{2}^{2}+{{KK}}_{6}{D}_{4}^{2}{D}_{5}{D}_{2}{\tilde{n}}_{\mathrm{CO}}\\ &&\quad -2{K}^{2}{K}_{6}^{2}{D}_{4}^{4}{D}_{5}^{2}{\tilde{n}}_{{\rm{N}}}{\tilde{n}}_{\mathrm{CO}}^{2}+{K}_{6}^{2}{D}_{4}^{3}{D}_{2}{\tilde{n}}_{\mathrm{CO}}^{2}=0.\end{eqnarray} \tag{ 15 }$$

Evaluating ${D}_{2}^{2}$ is particularly tedious (see Appendix C). It is worth nothing that D₁ and D₃ are essentially numbers that only depend on the normalized equilibrium constants, while D₂, D₄, and D₅ are functions of ${\tilde{n}}_{\mathrm{CO}}$. By "opening up" these terms in the preceding expression, we may re-express it as

$$\begin{eqnarray}&&\displaystyle \sum _{i=0}^{10}{A}_{i}{\tilde{n}}_{\mathrm{CO}}^{i}=0.\end{eqnarray} \tag{ 16 }$$

The coefficients of this decic equation are given in Appendix C.

We define the solar abundance of elements as being ${\tilde{n}}_{{\rm{O}}}=5\times {10}^{-4}$, ${\tilde{n}}_{{\rm{C}}}=2.5\times {10}^{-4}$, and ${\tilde{n}}_{{\rm{N}}}={10}^{-4}$. This implies that the carbon-to-oxygen and nitrogen-to-oxygen ratios are, respectively,

$$\begin{eqnarray}&&{\rm{C}}/{\rm{O}}\equiv \displaystyle \frac{{\tilde{n}}_{{\rm{C}}}}{{\tilde{n}}_{{\rm{O}}}}=0.5\,\mathrm{and}\,{\rm{N}}/{\rm{O}}\equiv \displaystyle \frac{{\tilde{n}}_{{\rm{N}}}}{{\tilde{n}}_{{\rm{O}}}}=0.2.\end{eqnarray} \tag{ 17 }$$

Unless otherwise stated, these are our default parameter values.

5.1. Benchmarking

We validate our analytical formulae by comparing them to calculations done using the TEA code (Blecic et al. 2016), which performs Gibbs free energy minimization. We emphasize that we use TEA at its full capability and include species beyond the 9 we have considered in our analytical formulae: H, C, O, H₂, O₂, CO, CO₂, CH₄, H₂O, CH, CH₂, CH₃, C₂H₂, OH, H₂CO, HCO, C₂, C₂H, C₂H₄, N₂, NH₃, HCN. In attempting to solve the decic equation in (16), we find that the procedure is sometimes numerically unstable, because the values of the various A_i coefficients may vary by many orders of magnitude. We emphasize that this is an issue of implementation and not of theory. In practice, it is sufficient to obtain the mixing ratio of CO by solving the quintic equation in (8), which is numerically stable. We use the polyroots routine in Python. The other mixing ratios are then obtained using Equations (2), (5), and (11).

In Figures 1–3, we represent calculations from the TEA code as circles overplotted on our calculations, which are shown as curves. Using the TEA calculations as a reference, we calculate the errors associated with our analytical formulae. For Figure 1, we find that the errors are $\sim 1 \%$, except for C₂H₂ and C₂H₄ with ${\rm{C}}/{\rm{O}}=1$ (middle panel), where they are $\sim 10 \%$, but the increased inaccuracy is due to our use of the quintic, rather than the decic, equation. For Figures 2 and 3, we find that the errors are $\sim 1 \%$ or smaller.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.2. Basic Trends

Figure 1 shows the trends associated with the mixing ratios versus temperature for solar-abundance, ${\rm{C}}/{\rm{O}}=1$, and nitrogen-rich atmospheres. The trend of CH₄ and CO being the dominant carbon carriers at low and high temperatures, respectively, persists even in the presence of nitrogen (Madhusudhan 2012). Analogously, NH₃ and N₂ are the dominant nitrogen carriers at low and high temperatures, respectively (Burrows & Sharp 1999; Lodders & Fegley 2002; Moses et al. 2011). HCN closely tracks the rise of CO with temperature. In a ${\rm{C}}/{\rm{O}}=1$ environment, HCN inhibits the formation of CH₄ and C₂H₂, as we can see from comparing the solutions derived in the current study versus the nitrogen-free solutions of Heng & Lyons (2016). Even when N/O is increased tenfold from 0.2 to 2, the trends produced resemble those of the solar-abundance case.

Figure 2 shows the mixing ratios versus C/O. The low- and high-temperature trends have previously been elucidated, namely that carbon-rich atmospheres are water-poor and methane-rich (Madhusudhan 2012; Moses et al. 2013a; Heng et al. 2016). HCN closely tracks the abundance of CH₄ as C/O increases, suggesting that the ratio of their abundances should be a constant. Figure 3 shows that the mixing ratios are somewhat insensitive to N/O. Unsurprisingly, only the nitrogen-bearing species (N₂, NH₃ and HCN) show any dependence of their mixing ratios on N/O.

Zoom In Zoom Out Reset image size

We note that our formulae do not consider graphite formation, which is expected to occur in carbon-rich atmospheres (Moses et al. 2013b). For this reason, we urge caution when applying these formulae to ${\rm{C}}/{\rm{O}}\gt 1$ situations.

5.3. Observational Diagnostics

Figure 4 shows the ratio of abundances of various pairs of molecules. The C₂H₂/CO₂ ratio is a sensitive diagnostic for C/O (Venot et al. 2015; Heng & Lyons 2016), spanning more than 10 orders of magnitude as C/O varies from 0.1 to 10, suggesting that this ratio may be used as an observational diagnostic for inferring the value of C/O, provided a given spectrum of an exoplanetary atmosphere has the sufficient resolution and signal-to-noise for such an inference to be made via an inversion technique. The CH₄/H₂O ratio is somewhat less sensitive to C/O, but provides an additional check on the inferred value of C/O. The HCN/CH₄ ratio is essentially constant across a factor of 100 in C/O and its value depends only on the metallicity, implying that it may be used as a robust diagnostic for the metallicity of the atmosphere.

Zoom In Zoom Out Reset image size

We have developed a novel analytical method for computing the abundances of six and nine molecules in a C–H–O–N system in chemical equilibrium. Our work demonstrates a useful trick, which is that trace molecules may formally be left out of the system of nonlinear equations and computed later using the mixing ratios of other molecules. Since our formulae have been successfully validated by a Gibbs free energy minimization code, they may be used to benchmark chemical kinetics codes. Reproducing chemical equilibrium is a key test of a chemical kinetics code. Our formulae may also be used in retrieval calculations to enforce chemical equilibrium throughout the atmosphere. Such an approach may be used to test the Bayesian evidence for chemical disequilibrium when interpreting the spectrum of an exoplanetary atmosphere.

We acknowledge financial and administrative support from the Center for Space and Habitability (CSH), the PlanetS NCCR framework and the Swiss-based MERAC Foundation.

All data have been compiled using the NIST-JANAF database (http://kinetics.nist.gov/janaf/). Note that the molar Gibbs free energy associated with N₂ and H₂ are 0 J mol⁻¹ by definition. The Gibbs free energy of formation for C₂H₄ are (in units of kJ/mol and from 500 to 3000 K, in intervals of 100 K): 80.933, 88.017, 95.467, 103.180, 111.082, 119.122, 127.259, 135.467, 143.724, 152.016, 160.331, 168.663, 177.007, 185.357, 193.712, 202.070, 210.429, 218.790, 227.152, 235.515, 243.880, 252.246, 260.615, 268.987, 277.363, 285.743. For NH₃, we have: 4.800, 15.879, 27.190, 38.662, 50.247, 61.910, 73.625, 85.373, 97.141, 108.918, 120.696, 132.469, 144.234, 155.986, 167.725, 179.447, 191.152, 202.840, 214.509, 226.160, 237.792, 249.406, 261.003, 272.581, 284.143, 295.689. For HCN, we have: 117.769, 114.393, 111.063, 107.775, 104.525, 101.308, 98.120, 94.955, 91.812, 88.687, 85.579, 82.484, 79.403, 76.333, 73.274, 70.226, 67.187, 64.158, 61.138, 58.127, 55.124, 52.130, 49.144, 46.167, 43.198, 40.237. For ${\rm{\Delta }}{\tilde{G}}_{\mathrm{0,4}}$, we have: 116.519, 103.718, 90.63, 77.354, 63.959, 50.485, 36.967, 23.421, 9.864, −3.697, −17.253, −30.802, −44.342, −57.87, −71.385, −84.888, −98.377, −111.855, −125.322, −138.777, −152.222, −165.657, −179.085, −192.504, −205.916, −219.322. For ${\rm{\Delta }}{\tilde{G}}_{\mathrm{0,5}}$, we have: −9.6, −31.758, −54.38, −77.324, −100.494, −123.82, −147.25, −170.746, −194.282, −217.836, −241.392, −264.938, −288.468, −311.972, −335.45, −358.894, −382.304, −405.68, −429.018, −452.32, −475.584, −498.812, −522.006, −545.162, −568.286, −591.378. For ${\rm{\Delta }}{\tilde{G}}_{\mathrm{0,6}}$, we have: 145.71, 121.401, 96.516, 71.228, 45.662, 19.906, −5.977, −31.942, −57.955, −83.992, −110.035, −136.073, −162.096, −188.098, −214.075, −240.023, −265.94, −291.826, −317.679, −343.5, −369.29, −395.047, −420.775, −446.472, −472.141, −497.784.

If only CO₂ is excluded, then we end up with a hexic/sextic equation for ${\tilde{n}}_{\mathrm{CO}}$ with the following coefficients,

$$\begin{eqnarray}{A}_{0} & = & 2{K}_{5}{J}_{0},\\ {A}_{1} & = & 2{K}_{5}{J}_{1}+2{{KK}}_{6}{\tilde{n}}_{{\rm{O}}}{F}_{0},\\ {A}_{2} & = & 2{K}_{5}{J}_{2}+{{KK}}_{6}(2{\tilde{n}}_{{\rm{O}}}{F}_{1}-{F}_{0})+{K}_{6}^{2}{F}_{0}-8{K}^{2}{K}_{6}^{2}{\tilde{n}}_{{\rm{O}}}^{2}{\tilde{n}}_{{\rm{N}}},\\ {A}_{3} & = & 2{K}_{5}{J}_{3}+{{KK}}_{6}(2{\tilde{n}}_{{\rm{O}}}{F}_{2}-{F}_{1})+{K}_{6}^{2}{F}_{1}+8{K}^{2}{K}_{6}^{2}{\tilde{n}}_{{\rm{O}}}{\tilde{n}}_{{\rm{N}}},\\ {A}_{4} & = & 2{K}_{5}{J}_{4}+{{KK}}_{6}(2{\tilde{n}}_{{\rm{O}}}{F}_{3}-{F}_{2})+{K}_{6}^{2}{F}_{2}-2{K}^{2}{K}_{6}^{2}{\tilde{n}}_{{\rm{N}}},\\ {A}_{5} & = & 2{K}_{5}{J}_{5}+{K}_{6}{F}_{3}({K}_{6}-K),\\ {A}_{6} & = & 2{K}_{5}{J}_{6},\end{eqnarray} \tag{ 18 }$$

where ${J}_{0}={F}_{0}^{2}$, ${J}_{1}=2{F}_{0}{F}_{1}$, ${J}_{2}=2{F}_{0}{F}_{2}+{F}_{1}^{2}$, ${J}_{3}=2{F}_{0}{F}_{3}+2{F}_{1}{F}_{2}$, ${J}_{4}=2{F}_{1}{F}_{3}+{F}_{2}^{2}$, ${J}_{5}=2{F}_{2}{F}_{3}$, ${J}_{6}={F}_{3}^{2}$, ${F}_{0}=8{K}^{2}{\tilde{n}}_{{\rm{O}}}^{2}{\tilde{n}}_{{\rm{C}}}$, ${F}_{1}=-2K{\tilde{n}}_{{\rm{O}}}[1+2K(2{\tilde{n}}_{{\rm{C}}}+{\tilde{n}}_{{\rm{O}}})]$, ${F}_{2}=K-2{K}_{3}{D}_{3}+2{K}^{2}({\tilde{n}}_{{\rm{C}}}+2{\tilde{n}}_{{\rm{O}}})$, and ${F}_{3}=-{K}^{2}$.

To render the algebra tractable, we have written

$$\begin{eqnarray}&&{D}_{2}^{2}=\displaystyle \sum _{i=0}^{10}{J}_{i}{\tilde{n}}_{\mathrm{CO}}^{i},\end{eqnarray} \tag{ 19 }$$

where the coefficients are ${J}_{0}={F}_{0}^{2}$, ${J}_{1}=2{F}_{0}{F}_{1}$, ${J}_{2}={F}_{1}^{2}\ +2{F}_{0}{F}_{2}$, ${J}_{3}=2{F}_{0}{F}_{3}+2{F}_{1}{F}_{2}$, ${J}_{4}=2{F}_{0}{F}_{4}+2{F}_{1}{F}_{3}+{F}_{2}^{2}$, ${J}_{5}\ =2{F}_{0}{F}_{5}+2{F}_{1}{F}_{4}+2{F}_{2}{F}_{3}$, ${J}_{6}=2{F}_{1}{F}_{5}+2{F}_{2}{F}_{4}+{F}_{3}^{2}$, ${J}_{7}\ =2{F}_{2}{F}_{5}+2{F}_{3}{F}_{4}$, ${J}_{8}=2{F}_{3}{F}_{5}+{F}_{4}^{2}$, ${J}_{9}=2{F}_{4}{F}_{5}$, and ${J}_{10}={F}_{5}^{2}$. The F_i coefficients are defined in Equation (14). For convenience, we also write ${C}_{2}\equiv 1/{K}_{2}$. The coefficients of Equation (16) are:

$$\begin{eqnarray}{A}_{0} & = & 2{K}_{5}{J}_{0},\\ {A}_{1} & = & 2{K}_{5}{J}_{1}+2{{KK}}_{6}{\tilde{n}}_{{\rm{O}}}{F}_{0},\\ {A}_{2} & = & 2{K}_{5}{J}_{2}+{{KK}}_{6}[2{\tilde{n}}_{{\rm{O}}}{F}_{1}+{F}_{0}(8{C}_{2}{\tilde{n}}_{{\rm{O}}}-1)]\\ & & +{K}_{6}^{2}{F}_{0}-8{K}^{2}{K}_{6}^{2}{\tilde{n}}_{{\rm{O}}}^{2}{\tilde{n}}_{{\rm{N}}},\\ {A}_{3} & = & 2{K}_{5}{J}_{3}+{{KK}}_{6}[2{\tilde{n}}_{{\rm{O}}}{F}_{2}+{F}_{1}(8{C}_{2}{\tilde{n}}_{{\rm{O}}}-1)\\ & & +4{F}_{0}{C}_{2}(2{C}_{2}{\tilde{n}}_{{\rm{O}}}-1)]+{K}_{6}^{2}({F}_{1}+6{F}_{0}{C}_{2})\\ & & +8{K}^{2}{K}_{6}^{2}{\tilde{n}}_{{\rm{O}}}{\tilde{n}}_{{\rm{N}}}(1-8{C}_{2}{\tilde{n}}_{{\rm{O}}}),\\ {A}_{4} & = & 2{K}_{5}{J}_{4}+{{KK}}_{6}[2{\tilde{n}}_{{\rm{O}}}{F}_{3}+{F}_{2}(8{C}_{2}{\tilde{n}}_{{\rm{O}}}-1)\\ & & +4{F}_{1}{C}_{2}(2{C}_{2}{\tilde{n}}_{{\rm{O}}}-1)-4{F}_{0}{C}_{2}^{2}]\\ & & +{K}_{6}^{2}({F}_{2}+6{F}_{1}{C}_{2}+12{F}_{0}{C}_{2}^{2})\\ & & -2{K}^{2}{K}_{6}^{2}{\tilde{n}}_{{\rm{N}}}(1-32{C}_{2}{\tilde{n}}_{{\rm{O}}}+96{C}_{2}^{2}{\tilde{n}}_{{\rm{O}}}^{2}),\\ {A}_{5} & = & 2{K}_{5}{J}_{5}+{{KK}}_{6}[2{\tilde{n}}_{{\rm{O}}}{F}_{4}+{F}_{3}(8{C}_{2}{\tilde{n}}_{{\rm{O}}}-1)\\ & & +4{F}_{2}{C}_{2}(2{C}_{2}{\tilde{n}}_{{\rm{O}}}-1)-4{F}_{1}{C}_{2}^{2}]\\ & & +{K}_{6}^{2}({F}_{3}+6{F}_{2}{C}_{2}+12{F}_{1}{C}_{2}^{2}+8{F}_{0}{C}_{2}^{3})\\ & & -16{K}^{2}{K}_{6}^{2}{\tilde{n}}_{{\rm{N}}}{C}_{2}(1-12{C}_{2}{\tilde{n}}_{{\rm{O}}}+16{C}_{2}^{2}{\tilde{n}}_{{\rm{O}}}^{2}),\\ {A}_{6} & = & 2{K}_{5}{J}_{6}+{{KK}}_{6}[2{\tilde{n}}_{{\rm{O}}}{F}_{5}+{F}_{4}(8{C}_{2}{\tilde{n}}_{{\rm{O}}}-1)\\ & & +4{F}_{3}{C}_{2}(2{C}_{2}{\tilde{n}}_{{\rm{O}}}-1)-4{F}_{2}{C}_{2}^{2}]\\ & & +{K}_{6}^{2}({F}_{4}+6{F}_{3}{C}_{2}+12{F}_{2}{C}_{2}^{2}+8{F}_{1}{C}_{2}^{3})\\ & & -16{K}^{2}{K}_{6}^{2}{\tilde{n}}_{{\rm{N}}}{C}_{2}^{2}(3-16{C}_{2}{\tilde{n}}_{{\rm{O}}}+8{C}_{2}^{2}{\tilde{n}}_{{\rm{O}}}^{2}),\\ {A}_{7} & = & 2{K}_{5}{J}_{7}+{{KK}}_{6}[{F}_{5}(8{C}_{2}{\tilde{n}}_{{\rm{O}}}-1)\\ & & +4{F}_{4}{C}_{2}(2{C}_{2}{\tilde{n}}_{{\rm{O}}}-1)-4{F}_{3}{C}_{2}^{2}]\\ & & +{K}_{6}^{2}{F}_{5}+2{K}_{6}^{2}{C}_{2}(3{F}_{4}+6{F}_{3}{C}_{2}+4{F}_{2}{C}_{2}^{2})\\ & & -64{K}^{2}{K}_{6}^{2}{\tilde{n}}_{{\rm{N}}}{C}_{2}^{3}(1-2{C}_{2}{\tilde{n}}_{{\rm{O}}}),\\ {A}_{8} & = & 2{K}_{5}{J}_{8}+4{{KK}}_{6}{C}_{2}[{F}_{5}(2{C}_{2}{\tilde{n}}_{{\rm{O}}}-1)-{F}_{4}{C}_{2}]\\ & & +2{K}_{6}^{2}{C}_{2}(3{F}_{5}+6{F}_{4}{C}_{2}+4{F}_{3}{C}_{2}^{2})-32{K}^{2}{K}_{6}^{2}{\tilde{n}}_{{\rm{N}}}{C}_{2}^{4},\\ {A}_{9} & = & 2{K}_{5}{J}_{9}-4{{KK}}_{6}{F}_{5}{C}_{2}^{2}+4{K}_{6}^{2}{C}_{2}^{2}(3{F}_{5}+2{F}_{4}{C}_{2}),\\ {A}_{10} & = & 2{K}_{5}{J}_{10}+8{K}_{6}^{2}{F}_{5}{C}_{2}^{3}.\end{eqnarray} \tag{ 20 }$$

From a practical standpoint, if one was implementing these expressions in a computer code, one would first code up K and K_i, followed by F_i and J_i, which would allow the construction of A_i.

We thank the developers of the Thermochemical Equilibrium Abundances (TEA) code (Blecic et al. 2016), initially developed at the University of Central Florida, Orlando, Florida, USA. The Reproducible Research Compendium and the Python code we used to produce Figure 1 are available at http://github.com/exoclime/VULCAN.