Toward Consistent Modeling of Atmospheric Chemistry and Dynamics in Exoplanets: Validation and Generalization of the Chemical Relaxation Method

The study of exoplanets has evolved from detection to characterization, thanks to the advent of cutting-edge observational techniques. The spectra of exo-atmospheres provide us with valuable clues about the atmospheric chemistry and thermal structure. Diagnosing and interpreting these spectra to obtain chemical compositions is now at the forefront of exo-atmospheric research.

The simplest assumption is to build a model in chemical equilibrium, where the molecular composition for a given elemental abundance only depends on local, basic parameters (pressure and temperature) independent of the reaction pathways. However, equilibrium chemistry only holds in the hot (T ≳ 2000 K) or deep (P ≳ 100 bar) parts of the atmosphere. Processes like ultraviolet (UV) irradiation and atmospheric dynamics drive the chemical composition in the observable atmosphere away from equilibrium. These disequilibrium processes commonly dominate the observable parts of atmospheres in the solar system.

Chemical kinetics models, which incorporate a chemical network of hundreds to thousands of reactions, are needed to study disequilibrium chemistry (e.g., Moses et al. 2011; Kopparapu et al. 2012; Venot et al. 2012; Grassi et al. 2014; Hu & Seager 2014; Rimmer & Helling 2016; Tsai et al. 2017). Even with such a number of reactions, chemical networks are greatly reduced compared to Nature and restricted to measurements in specific temperature ranges. In other words, there does not exist a one-fits-all chemical network in practice. The accuracy of a network is determined by whether it includes the relevant reactions and if the input rate coefficients are reliable. Unfortunately, these chemical kinetics models are expensive to run, especially when one desires to couple chemistry to three-dimensional atmospheric motion. In one dimension, eddy diffusion is used to mimic large-scale atmospheric circulation, convection, turbulence, etc. The three-dimensional atmospheric circulation patterns of tidally locked, highly irradiated exoplanets are demonstrably more complicated (e.g., Showman et al. 2009; Dobbs-Dixon & Agol 2013, and see Heng & Showman 2015 for a review). In order to correctly interpret the transmission spectra of hot Jupiters, there is a need to develop three-dimensional general circulation models (GCMs) that couple the atmospheric dynamics and chemistry. Intermediate steps have already been taken in this direction by, e.g., Agúndez et al. (2012, 2014), who coupled a chemical kinetics code to a simplified model of atmospheric dynamics (constant solid-body rotation mimicking a uniform equatorial jet). There is clearly a need to build on studies like these, but a brute-force coupling between a three-dimensional solver of the fluid equations and a chemical kinetics code with a network of hundreds to thousands of reactions is computationally challenging, even without considering radiative transfer (which is needed to include photochemistry).

Another approach is to simplify the chemical scheme. Conceptually, the interaction between atmospheric motion and chemistry is a comparison between two timescales: the dynamical versus chemical timescales. The simplest approach is to use the timescale argument, or quenching approximation. It assumes that the deep atmosphere is in chemical equilibrium, because the chemical timescale is much shorter than the dynamical timescale. There is a location within the atmosphere known as the quench point, where the timescales are equal. Above the quench point, the chemical abundances are assumed to be well-mixed and frozen to their equilibrium values at the quench point. In technical parlance, the process is referred to as "transport-induced quenching." The quenching approximation has been used to understand the over-abundance of carbon monoxide (CO) in the upper troposphere of Jupiter (Prinn & Barshay 1977). Visscher (2012) compared the chemical timescale for converting methane (CH₄) to CO to the orbital timescales of highly eccentric exoplanets to study the interaction between the evolving thermal structures and the chemistry. While the quenching approximation can be used to build the bulk of our intuition, it has been shown that it should be applied with caution to know when the abundance of a molecule is controlled by the disequilibrium abundance of another parent molecule, e.g., acetylene being controlled by methane (Tsai et al. 2017). Additionally, the quenching approximation contains the ambiguity of having to specify an appropriate length scale, which is not known from first principles (Smith 1998).

A more realistic approach is chemical relaxation. It was pioneered in the exo-atmospheres literature by Cooper & Showman (2006). Instead of simplifying the treatment of atmospheric dynamics, chemical relaxation takes the approach of replacing the chemical network with a single source/sink term that depends on the chemical timescale. Cooper & Showman coupled chemical relaxation to a simplified GCM to study quenching of CO, and suggested that most of the carbon is locked up in CO in HD 209458b due to transport. However, a shortcoming of their study is the assumption of a single rate-limiting reaction for the interconversion and probably the underestimation of the timescale of CO.

A prerequisite for implementing chemical relaxation is to be able to compute the chemical timescale over a broad range of temperatures and pressures. Initially, we had hoped to seek a universal fit for the chemical timescale across temperature, pressure, and metallicity, motivated by the work of Zahnle & Marley (2014). Figure 1 shows our attempt at fitting an Arrhenius-like function to the chemical timescale for CH₄–CO conversion. We find the timescales explored in Zahnle & Marley (2014) are likely restricted to a narrow temperature–pressure range. It is apparent that such a simple approach fails to fit the timescale across the range of values of temperature and pressure needed for us to implement chemical relaxation.

Zoom In Zoom Out Reset image size

Generally, the chemical timescale is a function that draws upon different rate-limiting reactions in different temperature and pressure regimes. A challenge with the chemical relaxation method is to find these rate-limiting reactions. In doing so, we develop a simple method to identify the dominant pathway and the associated rate-limiting step (RLS) in the chemical network for given values of the temperature, pressure, and elemental abundances. The pathway analysis can be used for comparing chemical kinetics calculations performed by different groups using different chemical networks, as it allows one to identify the key reactions that essentially control the output, and assess if key reactions are missing. It is useful both from a practical point of view and for developing physical intuition (e.g., Turányi 1990).

Ever since the work of Cooper & Showman (2006), only Drummond et al. (2018) have recently implemented the relaxation method in the Meteorological Office Unified Model and improved it by consistently coupling to radiative transfer. However, to our knowledge, chemical relaxation has surprisingly never been validated. By "validation," we mean that the accuracy of chemical relaxation should be demonstrated to a factor of a few, rather than to ∼1% accuracy (or better), given the existing uncertainties in the rate coefficients and the approximate nature of the approach. In addition, transit radii are proportional to the logarithms of the chemical abundances; such a factor-of-several validation suffices for studying atmospheric chemistry in hot Jupiters. In the current study, we perform this validation step in one dimension, since it is a test of the ability of the chemical relaxation scheme to mimic the full chemical network, rather than one of the complicated three-dimensional geometry. In a future work, we will aim to couple the chemical relaxation scheme to a three-dimensional GCM, but here we restrict ourselves to validating the scheme on a factor-of-several basis. Such a strategy follows the well-established hierarchical approach of constructing climate models (Held 2005). We do not consider photochemistry for the current study.

In Section 2, we provide the background theory on how to compute the chemical timescale. In Section 3, we describe our methodology for identifying the rate-limiting chemical reactions. In Section 4, we describe how we compute the chemical timescales. In Section 5, we compare our chemical network to that of Moses and Venot using our pathway analysis tool. In Section 6, we validate the chemical relaxation method for three model atmospheres. In Section 7, we summarize our findings and list opportunities for future work.

To compute the chemical timescale for use in the chemical relaxation method, we need to first derive a simplified expression for it that depends only on the local conditions of temperature and pressure. This allows us to evaluate the rate of change of the abundances efficiently.

Consider a chemical species subject to production (${ \mathcal P }$) and loss (${ \mathcal L }$) through a network of chemical reactions without any disequilibrium process (e.g., transport). The rate of change of its volume number density (n) is

$$\begin{eqnarray}&&\displaystyle \frac{{dn}}{{dt}}={ \mathcal P }-n{ \mathcal L }.\end{eqnarray} \tag{ 1 }$$

Equation (1) is written in such a way that ${ \mathcal P }$ and ${ \mathcal L }$ do not depend on n, as production only depends on other species and loss depends linearly on the number density of the reactant for a typical bimolecular reaction. Since ${ \mathcal P }$ and ${ \mathcal L }$ depend on the number densities of other species, the ensemble of Equation (1) for every species forms a system of coupled differential equations and the calculation involves inverting a large matrix (e.g., Hu et al. 2012; Tsai et al. 2017). If one wishes to implement chemical kinetics into the dynamical core of a GCM, then one needs to include a separate Euler equation for every species in the chemical network.

The relaxation method rewrites Equation (1) as (Smith 1998; Cooper & Showman 2006)

$$\begin{eqnarray}&&\displaystyle \frac{{dn}}{{dt}}=-\displaystyle \frac{n-{n}_{\mathrm{EQ}}}{{\tau }_{\mathrm{chem}}},\end{eqnarray} \tag{ 2 }$$

replacing the production and loss terms with a source/sink term that relaxes n to n_EQ by the chemical timescale (τ_chem). When the abundance is greater (less) than its equilibrium value, then it decreases (increases) and works its way toward equilibrium. The chemical timescale is effectively determined by the employed chemical network. When coupled to a GCM, whether the species in question attains chemical equilibrium depends on the competition between atmospheric dynamics and chemistry via their corresponding timescales. Chemical relaxation is analogous to the Newtonian relaxation method used as a substitute for radiative transfer (e.g., Held & Suarez 1994) or the treatment of condensation where the supersaturated gas is relaxed to the saturated vapor number density on the condensation timescale (Hu et al. 2012).

The chemical timescale is conventionally expressed without justification as

$$\begin{eqnarray}&&{\tau }_{\mathrm{chem}}=\displaystyle \frac{n}{| {dn}/{dt}| }=\displaystyle \frac{{n}_{\mathrm{EQ}}}{{n}_{\mathrm{EQ}}^{{\prime} }{{ \mathcal L }}_{\mathrm{EQ}}^{{\prime} }},\end{eqnarray} \tag{ 3 }$$

where the second equality assumes that the species in question has a number density that is equal to its chemical equilibrium value. Furthermore, ${n}_{\mathrm{EQ}}^{{\prime} }{{ \mathcal L }}_{\mathrm{EQ}}^{{\prime} }$ refers to the loss rate determined by the RLS involving other species but not n, denoted with primes. (The ambiguity comes from using the equilibrium abundance for the numerator but not the denominator in Equation (3) since dn/dt vanishes in chemical equilibrium.)

We now wish to demonstrate that Equation (3) may be derived from Equations (1) and (2). First, consider the situation when n ≪ n_EQ. This implies that the loss of the species being considered is negligible compared to production, which means ${dn}/{dt}\approx { \mathcal P }$. If we assume that all of the other species in the network are close to chemical equilibrium (as most of the intermediate species are fast-reacting radicals), then we can further write ${ \mathcal P }\approx {{ \mathcal P }}_{\mathrm{EQ}}^{{\prime} }$ because we recall that ${ \mathcal P }$ does not depend on n. For the RLS, we can write ${{ \mathcal P }}_{\mathrm{EQ}}^{{\prime} }={n}_{\mathrm{EQ}}^{{\prime} }{{ \mathcal L }}_{\mathrm{EQ}}^{{\prime} }$. By inserting this expression into Equations (1) and (2), we obtain Equation (3). In the opposite limit of n ≫ n_EQ, the loss of the species dominates production and we have ${dn}/{dt}\,=-n{ \mathcal L }\approx -{n}_{\mathrm{EQ}}^{{\prime} }{{ \mathcal L }}_{\mathrm{EQ}}^{{\prime} }$. This again leads to Equation (3) from Equations (1) and (2). Since the chemical timescale expression is approximately correct for both limits, it is expected to work at order-of-magnitude accuracy at least. This expectation will be confirmed by full numerical calculations of chemical kinetics.

It is common for the chemical conversion of one species to another not to occur in one step. Rather, it takes multiple steps to surmount the energy barrier via the breaking or forming of chemical bonds. These sequences of reactions form a pathway, and the chemical timescales associated with each step in the pathway may differ by many orders of magnitude. The efficiency of a pathway is bottlenecked by its slowest reaction. Specifically, the RLS is defined as the slowest reaction along the fastest pathway. It informs the effective loss rate, ${n}_{\mathrm{EQ}}^{{\prime} }{{ \mathcal L }}_{\mathrm{EQ}}^{{\prime} }$, in Equation (3). Thus, computing the chemical timescale involves identifying the RLS. For example, Zahnle & Marley (2014) have remarked how the conversion of CO to CH₄ may be visualized as the reduction of the bond between C and O from a triple bond to a double bond to a single bond and eventually splitting C from O, in three steps: first between CO and formaldehyde (H₂CO), second between formaldehyde and methanol (CH₃OH), and finally between methanol and methane. We build upon and extend the diagram in Figure 2, where the temperature-and-pressure-dependent pathways and RLSs are included, as explained in the following subsections.

Zoom In Zoom Out Reset image size

3.1. CH₄–CO Interconversion

3.1.1. Identifying the Rate-limiting Steps from Full Chemical Kinetics

In addition to the RLS being a function of temperature and pressure, extracting this information to identify it is not straightforward because the possible number of paths grows exponentially with increasing number of species, making it difficult to track them all. Our approach is to develop a tool using Dijkstra's algorithm to find the shortest path and the associated RLS; see Appendix B for more details.

Figure 3 shows a survey of the different RLSs as functions of temperature and pressure and for protosolar elemental abundances from Lodders (2009) (C/H = 2.776 × 10⁻⁴, O/H =6.062 × 10⁻⁴, N/H = 8.185 × 10⁻⁵, He/H = 9.69 × 10⁻²). Unsurprisingly, there does not exist a single RLS for the range of temperatures (500–3000 K) and pressures (0.1 mbar–1 kbar) considered. The different reactions, labeled (R1)–(R10), are

$$\begin{eqnarray}&&\mathrm{OH}+{\mathrm{CH}}_{3}+{\rm{M}}\to {\mathrm{CH}}_{3}\mathrm{OH}+{\rm{M}}\end{eqnarray} \tag{ R1 }$$

$$\begin{eqnarray}&&\mathrm{OH}+{\mathrm{CH}}_{3}\to {\mathrm{CH}}_{2}\mathrm{OH}+{\rm{H}}\end{eqnarray} \tag{ R2 }$$

$$\begin{eqnarray}&&{\mathrm{CH}}_{2}\mathrm{OH}+{\rm{M}}\to {\rm{H}}+{{\rm{H}}}_{2}\mathrm{CO}+{\rm{M}}\end{eqnarray} \tag{ R3 }$$

$$\begin{eqnarray}&&{\mathrm{CH}}_{3}+{\rm{O}}\to {{\rm{H}}}_{2}\mathrm{CO}+{\rm{H}}\end{eqnarray} \tag{ R4 }$$

$$\begin{eqnarray}&&\mathrm{OH}+{\rm{C}}\to \mathrm{CO}+{\rm{H}}\end{eqnarray} \tag{ R5 }$$

$$\begin{eqnarray}&&{\rm{H}}+{\mathrm{CH}}_{4}\to {\mathrm{CH}}_{3}+{{\rm{H}}}_{2}\end{eqnarray} \tag{ R6 }$$

$$\begin{eqnarray}&&{\mathrm{CH}}_{3}\mathrm{OH}+{\rm{H}}\to {\mathrm{CH}}_{3}{\rm{O}}+{{\rm{H}}}_{2}\end{eqnarray} \tag{ R7 }$$

$$\begin{eqnarray}&&{{\rm{C}}}_{2}{{\rm{H}}}_{2}+{\rm{O}}\to {\mathrm{CH}}_{2}+\mathrm{CO}\end{eqnarray} \tag{ R8 }$$

$$\begin{eqnarray}&&{\mathrm{CH}}_{3}+{{\rm{H}}}_{2}{\rm{O}}\to {\mathrm{CH}}_{3}\mathrm{OH}+{\rm{H}}\end{eqnarray} \tag{ R9 }$$

$$\begin{eqnarray}&&{{\rm{C}}}_{2}{{\rm{H}}}_{2}+\mathrm{OH}\to {\mathrm{CH}}_{3}+\mathrm{CO}\end{eqnarray} \tag{ R10 }$$

where M refers to any third body (i.e., the total number density of the gas). Reactions (R8) and (R10) are only relevant when C/O > 1, as we will describe shortly in the example of C/O = 2.

Zoom In Zoom Out Reset image size

3.1.2. Grouping of Reactions Into Three Schemes

It is possible to understand CH₄–CO interconversion as consisting of three schemes (at least, for solar-like elemental abundances). As temperature increases, the scheme moves from (C) to (A), as higher kinetic energy allows more ambitious steps:

• (A)  
  though the progressive dehydrogenation of CH₄ into C, followed by oxidization into CO (blue arrows in Figure 2);
• (B)  
  via H₂CO (formaldehyde) from directly oxidizing CH₃ as described by (R4);
• (C)  
  via intermediate species like CH₂OH, CH₃OH, or CH₃O from oxidizing CH₃ (through the molecules shown in green in Figure 2).

At high temperatures and low pressures (the magenta region in Figure 3), scheme (A) is favored because it requires a high abundance of atomic hydrogen, produced mainly by thermal decomposition of H₂. Scheme (B) sits in the transition between scheme (A) and (C). Scheme (C) covers the broadest range of temperature and pressure and contains the pathways previously identified by, e.g., Yung et al. (1988), Bézard et al. (2002), Moses et al. (2011), and Visscher (2012). As shown in Figure 2, CH₄ is first converted to CH₃ before being oxidized by OH or H₂O to form CH₂OH or CH₃OH/CH₃O depending on the temperature and pressure. These singly bonded (C–O) intermediate species make forming the double bond (C=O) in H₂CO easier than directly from C and O. H₂CO goes on to efficiently produce HCO and finally the triple-bonded structure of CO. At high enough temperatures (T ≳ 2000 K), there is sufficient energy to directly form the double bond between C and O into H₂CO via reaction (R4) without passing through the intermediate species, which is scheme (B). This is similar to the RLS, ${\mathrm{CH}}_{3}+\mathrm{OH}\to {{\rm{H}}}_{2}+{{\rm{H}}}_{2}\mathrm{CO}$, initially suggested by Prinn & Barshay (1977), to explain the quenched CO found in Jupiter. At even higher temperatures, if the pressure remains low enough then scheme (A) dominates. Molecular hydrogen is dissociated into atomic hydrogen, which in turn promotes the dehydrogenation of methane. Eventually, the accumulated C is present at high enough abundances that allow for its oxidization into CO.

A typical pathway of scheme (A) is:

$$\begin{eqnarray*}&&{\mathrm{CH}}_{4}+{\rm{H}}\to {\mathrm{CH}}_{3}+{{\rm{H}}}_{2}\\ &&{\mathrm{CH}}_{3}+{\rm{H}}\to {\mathrm{CH}}_{2}+{{\rm{H}}}_{2}\\ &&{\mathrm{CH}}_{2}+{\rm{H}}\to \mathrm{CH}+{{\rm{H}}}_{2}\\ &&\mathrm{CH}+{\rm{H}}\to {\rm{C}}+{{\rm{H}}}_{2}\\ &&{\rm{C}}+\mathrm{OH}\to \mathrm{CO}+{\rm{H}}\\ &&{\rm{H}}+{{\rm{H}}}_{2}{\rm{O}}\to \mathrm{OH}+{{\rm{H}}}_{2}\\ &&\displaystyle \frac{2({{\rm{H}}}_{2}+{\rm{M}}\to 2{\rm{H}}+{\rm{M}})}{\mathrm{net}:{\mathrm{CH}}_{4}+{{\rm{H}}}_{2}{\rm{O}}\to \mathrm{CO}+3{{\rm{H}}}_{2}}\end{eqnarray*}$$

Scheme (B), which does not involve intermediate species like ${\mathrm{CH}}_{3}\mathrm{OH}$ or ${\mathrm{CH}}_{2}\mathrm{OH}$, goes through the pathway:

$$\begin{eqnarray*}&&{\mathrm{CH}}_{4}+{\rm{H}}\to {\mathrm{CH}}_{3}+{{\rm{H}}}_{2}\\ &&{\mathrm{CH}}_{3}+{\rm{O}}\to {{\rm{H}}}_{2}\mathrm{CO}+{\rm{H}}\\ &&{{\rm{H}}}_{2}\mathrm{CO}+{\rm{H}}\to \mathrm{HCO}+{{\rm{H}}}_{2}\\ &&\mathrm{HCO}+{\rm{H}}\to \mathrm{CO}+{{\rm{H}}}_{2}\\ &&\mathrm{OH}+{\rm{H}}\to {\rm{O}}+{{\rm{H}}}_{2}\\ &&{\rm{H}}+{{\rm{H}}}_{2}{\rm{O}}\to \mathrm{OH}+{{\rm{H}}}_{2}\\ &&\displaystyle \frac{2({{\rm{H}}}_{2}+{\rm{M}}\to 2{\rm{H}}+{\rm{M}})}{\mathrm{net}:{\mathrm{CH}}_{4}+{{\rm{H}}}_{2}{\rm{O}}\to \mathrm{CO}+3{{\rm{H}}}_{2}}\end{eqnarray*}$$

Scheme (C) includes several main pathways for T ≲ 2000 K or P ≳ 1 bar. One example of scheme (C) is:

$$\begin{eqnarray*}&&{\mathrm{CH}}_{4}+{\rm{H}}\to {\mathrm{CH}}_{3}+{{\rm{H}}}_{2}\\ &&{\mathrm{CH}}_{3}+\mathrm{OH}+{\rm{M}}\to {\mathrm{CH}}_{3}\mathrm{OH}+{\rm{M}}\\ &&{\mathrm{CH}}_{3}\mathrm{OH}+{\rm{H}}\to {\mathrm{CH}}_{3}{\rm{O}}+{{\rm{H}}}_{2}\\ &&{\mathrm{CH}}_{3}{\rm{O}}+{\rm{M}}\to {{\rm{H}}}_{2}\mathrm{CO}+{\rm{H}}+{\rm{M}}\\ &&{{\rm{H}}}_{2}\mathrm{CO}+{\rm{H}}\to \mathrm{HCO}+{{\rm{H}}}_{2}\\ &&\mathrm{HCO}+{\rm{M}}\to {\rm{H}}+\mathrm{CO}+{\rm{M}}\\ &&{\rm{H}}+{{\rm{H}}}_{2}{\rm{O}}\to \mathrm{OH}+{{\rm{H}}}_{2}\\ &&\displaystyle \frac{{{\rm{H}}}_{2}+{\rm{M}}\to 2{\rm{H}}+{\rm{M}}}{\mathrm{net}:{\mathrm{CH}}_{4}+{{\rm{H}}}_{2}{\rm{O}}\to \mathrm{CO}+3{{\rm{H}}}_{2}}.\end{eqnarray*}$$

More examples of pathways belonging to scheme (C) may be found in Moses et al. (2011), Rimmer & Helling (2016), and Tsai et al. (2017).

3.1.3. Comparison with Previous Work

Previous work studying CO–CH₄ interconversion has typically assumed one or two RLSs, because their intentions were to estimate the location of the quench point (Yung et al. 1988; Lodders & Fegley 2002; Madhusudhan & Seager 2011; Visscher & Moses 2011). For the purpose of implementing chemical relaxation, this is insufficient. In Figure 4, we demonstrate this by comparing the calculations of chemical timescales associated with CH₄ and CO with those from Cooper & Showman (2006) and Visscher (2012). The discrepancies between our calculations and those from Visscher (2012) are mainly present at high temperatures and low pressures. There are significant discrepancies between our calculations and the approximate ones of Cooper & Showman (2006) due to a different RLS assumed by them as explained in the following paragraph. We use Equation (16), instead of (19), of Cooper & Showman (2006), because the latter is an approximation that is valid only when T ≲ 2000 K and the mole fraction of H₂ is close to 1. Instead of using Equations (1)–(5) from Cooper & Showman (2006), we perform the full chemical equilibrium calculations. We also include chemical timescales associated with NH₃ and N₂ for completeness.

Zoom In Zoom Out Reset image size

Cooper & Showman (2006) consider a single CO–${\mathrm{CH}}_{4}$ pathway

$$\begin{eqnarray}&&{{\rm{H}}}_{2}\mathrm{CO}+{\rm{H}}+{\rm{M}}\to {\mathrm{CH}}_{3}{\rm{O}}+{\rm{M}}\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&{\mathrm{CH}}_{3}{\rm{O}}+{{\rm{H}}}_{2}\to {\mathrm{CH}}_{3}\mathrm{OH}+{\rm{H}}\end{eqnarray} \tag{ 4b }$$

$$\begin{eqnarray}&&{\mathrm{CH}}_{3}\mathrm{OH}+{\rm{H}}\to {\mathrm{CH}}_{3}+{{\rm{H}}}_{2}{\rm{O}}\end{eqnarray} \tag{ 4c }$$

proposed by Yung et al. (1988), where CO reacts with hydrogen to form ${{\rm{H}}}_{2}\mathrm{CO}$ (formaldehyde) and goes through ${\mathrm{CH}}_{3}{\rm{O}}$ (methoxide) and ${\mathrm{CH}}_{3}\mathrm{OH}$ (methanol) to get to ${\mathrm{CH}}_{3}$ (methyl). They suggested (4a), which is involved in breaking the C=O bond, as being the RLS. Cooper & Showman (2006) adopt the rate constants from Page et al. (1989) and Bézard et al. (2002) for the low- and high-pressure limits, respectively. We verify that the rate coefficient from Page et al. (1989) has a value that is similar to what we use in our chemical network and is not the source of the discrepancies between our calculations and those of Cooper & Showman (2006).

Rather, the discrepancies stem from the rate coefficient associated with (4b), which had not been measured experimentally. Yung et al. (1988) assumed this reaction to be relatively fast, based on comparison with other similar reactions (see their Appendix A). Motivated by the importance of CH₃OH kinetics (see Visscher et al. 2010 for details), Moses et al. (2011) performed ab initio calculations for the rate coefficients. According to their rate coefficients, (4b) always reacts slower than (4c) and thus should be the RLS⁵ We have chosen to use the rate coefficients of Moses et al. (2011).

Visscher (2012) identified (R1) and (R2) as being the RLSs and adopt the rate coefficients from Jasper et al. (2007). We have also taken the rate coefficients for (R1) from Jasper et al. (2007), but the reverse rate coefficient for (R2) from Tsang (1987). The differences between these rate coefficients are within a factor of 2. The two RLSs in Visscher (2012) control the most relevant temperature–pressure regions for hot/warm Jupiters (∼1000–1500 K) and their timescale agrees well with our calculation until entering the high-temperature and low-pressure regime.

3.2. NH₃–N₂ Interconversion

The timescales of nitrogen species are less constrained than that of CH₄–CO interconversion. HCN participates when T ≳ 1000 K and complicates the interconversion. We find it not straightforward to quantify the contribution of HCN, and the real timescales deviate from those simply considering NH₃–N₂ interconversion.

The RLSs for NH₃–N₂ interconversion are

$$\begin{eqnarray}&&{\mathrm{NH}}_{2}+{\mathrm{NH}}_{2}\to {{\rm{N}}}_{2}{{\rm{H}}}_{2}+{{\rm{H}}}_{2}\end{eqnarray} \tag{ R11 }$$

$$\begin{eqnarray}&&{{\rm{N}}}_{2}{{\rm{H}}}_{3}+{\rm{M}}\to {{\rm{N}}}_{2}{{\rm{H}}}_{2}+{\rm{H}}+{\rm{M}}\end{eqnarray} \tag{ R12 }$$

$$\begin{eqnarray}&&\mathrm{NH}+{\mathrm{NH}}_{2}\to {{\rm{N}}}_{2}{{\rm{H}}}_{2}+{\rm{H}}\end{eqnarray} \tag{ R13 }$$

$$\begin{eqnarray}&&\mathrm{NO}+{\mathrm{NH}}_{2}\to {{\rm{N}}}_{2}+{{\rm{H}}}_{2}{\rm{O}}\end{eqnarray} \tag{ R14 }$$

$$\begin{eqnarray}&&{\rm{N}}+\mathrm{NO}\to {{\rm{N}}}_{2}+{\rm{O}}\end{eqnarray} \tag{ R15 }$$

$$\begin{eqnarray}&&{{\rm{N}}}_{2}{{\rm{H}}}_{4}+{\rm{H}}\to {{\rm{N}}}_{2}{{\rm{H}}}_{3}+{{\rm{H}}}_{2}.\end{eqnarray} \tag{ R16 }$$

Figure 5 visualizes the network of RLSs for nitrogen chemistry, while Table 1 lists the RLSs of NH₃–N₂, across temperature and pressure.

Zoom In Zoom Out Reset image size

Table 1.  ${\mathrm{NH}}_{3}$ $\leftrightarrow$ ${{\rm{N}}}_{2}$ Rate-limiting Reactions

T (K)	500	1000	1500	2000	2500
P (bar)
10−4	R11	R14	R15	R15	R15
10−3	R11	R14	R13	R15	R15
10−2	R11	R13	R13	R15	R15
10−1	R11	R12	R13	R13	R15
1	R12	R12	R13	R13	R13
10	R11	R12	R13	R13	R13
100	R11	R12	R16	R13	R13

Download table as:  ASCIITypeset image

NH₃–N₂ interconversion can be divided into two schemes, depending on whether ${{\rm{N}}}_{2}$ is formed from ${{\rm{N}}}_{2}{\rm{H}}$ or NO. At high pressures (the parameter space occupied by (R11)–(R13) in Table 1), ${{\rm{N}}}_{2}$ is mainly formed by the dissociation of ${{\rm{N}}}_{2}{\rm{H}}$, with a pathway such as

$$\begin{eqnarray*}&&2({\mathrm{NH}}_{3}+{\rm{H}}\to {\mathrm{NH}}_{2}+{{\rm{H}}}_{2})\\ &&{\mathrm{NH}}_{2}+{\mathrm{NH}}_{2}\to {{\rm{N}}}_{2}{{\rm{H}}}_{2}+{{\rm{H}}}_{2}\\ &&{{\rm{N}}}_{2}{{\rm{H}}}_{2}+{\rm{H}}\to {{\rm{N}}}_{2}{\rm{H}}+{{\rm{H}}}_{2}\\ &&{{\rm{N}}}_{2}{\rm{H}}+{\rm{M}}\to {{\rm{N}}}_{2}+{\rm{H}}+{\rm{M}}\\ &&\displaystyle \frac{{{\rm{H}}}_{2}+{\rm{M}}\to 2{\rm{H}}\,+\,{\rm{M}}}{\mathrm{net}:2{\mathrm{NH}}_{3}\to {{\rm{N}}}_{2}+3{{\rm{H}}}_{2}},\end{eqnarray*}$$

where the second reaction (R11) is the RLS. As temperature increases, (R11) is replaced by (R13), or a channel through ${{\rm{N}}}_{2}{{\rm{H}}}_{3}$ and ${{\rm{N}}}_{2}{{\rm{H}}}_{4}$ (R12). This pathway is similar to that identified for HD 209458b by Moses et al. (2011).

At low pressures (the parameter space occupied by (R14) and (R15) in Table 1), H₂ is attacked by the more abundant free atomic O and produces OH, which in turn forms NO with N via ${\rm{N}}+\mathrm{OH}$ $\to$ $\mathrm{NO}+{\rm{H}}$. NO then reacts with N or NH₂, depending on the temperature, to produce N₂. This step involves forming the N ≡ N bond and is usually the RLS. The pathway becomes

$$\begin{eqnarray*}&&2({\mathrm{NH}}_{3}+{\rm{H}}\to {\mathrm{NH}}_{2}+{{\rm{H}}}_{2})\\ &&2({\mathrm{NH}}_{2}+{\rm{H}}\to \mathrm{NH}+{{\rm{H}}}_{2})\\ &&2(\mathrm{NH}+{\rm{H}}\to {\rm{N}}+{{\rm{H}}}_{2})\\ &&{\rm{N}}+\mathrm{NO}\to {{\rm{N}}}_{2}+{\rm{O}}\\ &&{\rm{O}}+{{\rm{H}}}_{2}\to \mathrm{OH}+{\rm{H}}\\ &&{\rm{N}}+\mathrm{OH}\to \mathrm{NO}+{\rm{H}}\\ &&\displaystyle \frac{2({{\rm{H}}}_{2}+M\to 2{\rm{H}}+M)}{\mathrm{net}:2{\mathrm{NH}}_{3}\to {{\rm{N}}}_{2}+3{{\rm{H}}}_{2}},\end{eqnarray*}$$

where (R15) is the RLS.

The full expressions for the chemical timescales, including those involving methane, carbon monoxide, water, ammonia, and molecular nitrogen, are stated in Appendix A with the rate-limiting reactions and rate coefficients listed in Table 4. In the following, we explain the reasoning behind their construction.

4.1. Revisiting CH₄–CO Interconversion

The study of interconversion between methane and carbon monoxide has a long and rich history. In the current study, we focus on identifying the pathways in reducing atmospheres, while the same steps can be applied to other types of atmospheres. In this subsection, our goal is to provide an analytical expression for computing the timescale associated with CH₄–CO interconversion as a function of the reactions (R1)–(R10). Table 2 shows the RLSs for solar metallicity, while Table 3 shows them for C/O = 2. Reactions (R8) and (R10) are only relevant for C/O = 2, as previously mentioned.

Table 2.  ${\mathrm{CH}}_{4}$ $\leftrightarrow$ $\mathrm{CO}$ Rate-limiting Reactions (Solar Abundance)

T (K)	500	1000	1500	2000	2500
P (bar)
10−4	R9	R2	R3	R5	R6
10−3	R9	R2	R3	R4	R5
10−2	R9	R2	R2	R4	R5
10−1	R9	R1	R2	R3	R4
1	R9	R1	R2	R3	R3
10	R9	R1	R1	R3	R3
100	R9	R1	R1	R7	R7

Download table as:  ASCIITypeset image

Table 3.  ${\mathrm{CH}}_{4}$ $\leftrightarrow$ $\mathrm{CO}$ Rate-limiting Reactions (C/O = 2)

T (K)	500	1000	1500	2000	2500
P (bar)
10−4	R9	R2	R8	R8	R5
10−3	R9	R2	R8	R8	R8
10−2	R9	R2	R10	R8	R8
10−1	R9	R1	R2	R8	R8
1	R9	R1	R2	R3	R8
10	R9	R1	R1	R3	R8
100	R9	R1	R7	R7	R7

Download table as:  ASCIITypeset image

Under differing conditions of temperature and pressure, the various RLSs can either collaborate or compete with one another. We find that a useful analogy for understanding the pathways is to visualize them as the resistors in an electrical circuit, built in series or in parallel. Using such an analogy, we can construct an analytical expression for the chemical timescale that consists of a network of reactions working either in series or in parallel. The upper diagram of Figure 6 visualizes how such an analogous electrical circuit would look like for CH₄–CO interconversion. We identify the series and parallel RLSs and group those that operate in similar temperatures and pressures together. For example, (R2) and (R3) are in series but (R2) and (R4) are in parallel. Depending on the temperature and pressure, either (R1), (R9) or the group of (R2)–(R4) is in control.

Zoom In Zoom Out Reset image size

Mathematically, a pair of reactions in parallel can be expressed as an operation that takes the maximum of the two rate coefficients. For a pair of reactions in series, the operation instead takes the minimum of the two rate coefficients. The groups that operate in particular temperatures and pressures are then added together for simplicity. In this way, the series of relationships between the reactions can be expressed as

$$\begin{eqnarray}{\tau }_{{\mathrm{CH}}_{4}} & = & \displaystyle \frac{[{\mathrm{CH}}_{4}]}{{\rm{r}}1+\max (\min ({\rm{r}}2,{\rm{r}}3),{\rm{r}}4)+{\rm{r}}9+\min ({\rm{r}}5,{\rm{r}}6)}\\ & & +{\tau }_{{{\rm{H}}}_{2}}\times \displaystyle \frac{3[\mathrm{CO}]}{[{{\rm{H}}}_{2}]},\end{eqnarray} \tag{ 5 }$$

where [X] represents the equilibrium number density of species X in cm⁻³, rx stands for the reaction rate of (Rx), and the factor $\tfrac{3[\mathrm{CO}]}{[{{\rm{H}}}_{2}]}$ stems from the amount of hydrogen that participates in the ${\mathrm{CH}}_{4}$–CO interconversion (three ${{\rm{H}}}_{2}$ for every CO according to the net reaction). The second term is usually orders of magnitude smaller than the first term except at high temperatures and low pressures, where the dissociation and recombination between H and H₂ become important and the relatively slower conversion of hydrogen starts to bottleneck the process. We demonstrate how methane is controlled by hydrogen in Figure 7 where we manually vary the rate of H–H₂ dissociation/recombination.

Zoom In Zoom Out Reset image size

We have also replaced (R7) with (R1) for simplicity because the rates of both reactions are very similar at high pressures. The same formula is applied to τ_CO except one needs to replace the numerator with [CO], since the interconversion goes both ways.

4.2. Water

H₂O reacts efficiently with atomic hydrogen via

$$\begin{eqnarray}&&{{\rm{H}}}_{2}{\rm{O}}+{\rm{H}}\to \mathrm{OH}+{{\rm{H}}}_{2}\end{eqnarray} \tag{ R17 }$$

in a hydrogen-rich atmosphere. Being the major oxygen carrier, water prevailingly participates in various reaction pathways, e.g., CO $\leftrightarrow$ CO₂ and H₂ $\leftrightarrow$ 2H. For solar metallicity, carbon is the richest heavy element next to oxygen. We find the timescale of water is effectively determined by the interconversion rate of CH₄–CO,

$$\begin{eqnarray}{\tau }_{{{\rm{H}}}_{2}{\rm{O}}} & = & \displaystyle \frac{[{{\rm{H}}}_{2}{\rm{O}}]}{{\rm{r}}1+\max (\min ({\rm{r}}2,{\rm{r}}3),{\rm{r}}4)+{\rm{r}}9+\min ({\rm{r}}5,{\rm{r}}6)}\\ & & +{\tau }_{{{\rm{H}}}_{2}}\times \displaystyle \frac{3[\mathrm{CO}]}{[{{\rm{H}}}_{2}]}.\end{eqnarray} \tag{ 6 }$$

4.3. Carbon Dioxide

Carbon dioxide is produced through the relatively fast scheme (Line et al. 2010; Moses et al. 2011)

$$\begin{eqnarray*}&&{{\rm{H}}}_{2}{\rm{O}}+{\rm{H}}\to \mathrm{OH}+{{\rm{H}}}_{2}\\ &&\displaystyle \frac{\mathrm{CO}+\mathrm{OH}\to {\mathrm{CO}}_{2}+{\rm{H}}}{\mathrm{net}:\mathrm{CO}+{{\rm{H}}}_{2}{\rm{O}}\to {\mathrm{CO}}_{2}+{{\rm{H}}}_{2}}.\end{eqnarray*}$$

The chemical timescale is simply

$$\begin{eqnarray}&&{\tau }_{{\mathrm{CO}}_{2}}=\displaystyle \frac{[{\mathrm{CO}}_{2}]}{{k}_{{\mathrm{CO}}_{2}}[\mathrm{CO}][\mathrm{OH}]},\end{eqnarray} \tag{ 7 }$$

where ${k}_{{\mathrm{CO}}_{2}}$ is the rate coefficient for $\mathrm{CO}+\mathrm{OH}\to {\mathrm{CO}}_{2}+{\rm{H}}$. Due to the fast conversion, CO₂ still maintains pseudo-equilibrium with CO and H₂O after the latter two are quenched, before CO₂ reaches its own quench point (see the discussion in Section 3.1 of Moses et al. 2011). Owing to this coupling, instead of relaxing CO₂ to its equilibrium value, we find it correct to relax CO₂ toward the pseudo-equilibrium value determined by the (possibly quenched) CO and H₂O, as expressed in Equation (18) below, similar to the treatment in Equation (43) of Zahnle & Marley (2014).

4.4. Ammonia and Molecular Nitrogen

The chemical timescale of ammonia is approximately determined by the ${\mathrm{NH}}_{3}$–${{\rm{N}}}_{2}$ conversion. Omitting HCN, we group the RLSs with similar operating temperatures and pressures and construct the effective electrical circuit in the lower diagram of Figure 6.

For ammonia production, we have

$$\begin{eqnarray}{\tau }_{{\mathrm{NH}}_{3}} & = & \displaystyle \frac{1}{2}\left(\displaystyle \frac{[{\mathrm{NH}}_{3}]}{\max ({\rm{r}}11,{\rm{r}}12)+{\rm{r}}13+{\rm{r}}14+{\rm{r}}15}+{\tau }_{{{\rm{H}}}_{2}}\right.\\ & & \left.\times \displaystyle \frac{3[{{\rm{N}}}_{2}]}{[{{\rm{H}}}_{2}]}\right),\end{eqnarray} \tag{ 8 }$$

where the factor $\tfrac{3[{{\rm{N}}}_{2}]}{[{{\rm{H}}}_{2}]}$ is again the amount of hydrogen that participates in the ${\mathrm{NH}}_{3}$–${{\rm{N}}}_{2}$ interconversion, as in the ${\mathrm{CH}}_{4}$–CO interconversion limited by H–H₂ interconversion at high temperatures and low pressures, and the factor of 1/2 comes from the fact that the net reaction converts two NH₃ molecules to one N₂ molecule.

Following the same steps, the chemical timescale of molecular nitrogen is expressed as

$$\begin{eqnarray}&&{\tau }_{{{\rm{N}}}_{2}}=\displaystyle \frac{[{{\rm{N}}}_{2}]}{\max ({\rm{r}}11,{\rm{r}}12)+{\rm{r}}13+{\rm{r}}14+{\rm{r}}15}+{\tau }_{{{\rm{H}}}_{2}}\times \displaystyle \frac{3[{{\rm{N}}}_{2}]}{[{{\rm{H}}}_{2}]}.\end{eqnarray} \tag{ 9 }$$

4.5. The Effects of Metallicity

4.6. The Effects of C/O

The carbon-to-oxygen ratio (C/O) is a crucial factor in controlling the atmospheric chemistry and thermal structure (Madhusudhan 2012; Moses et al. 2013a; Venot et al. 2015; Rocchetto et al. 2016). Equilibrium chemistry is sensitive to C/O and undergoes a qualitative transition at C/O = 1. We explore C/O values ranging from 0.1 to 2. The reason to limit ourselves at C/O = 2 is that C/O much larger than unity is considered unlikely as the surplus carbon tends to condense and form graphite (Moses et al. 2013b). Comparing Tables 2 and 3, the major change as C/O exceeds unity is that carbon takes the route through C₂H₂ at high temperature. The typical pathway in a hot, carbon-rich atmosphere is

$$\begin{eqnarray*}&&{\mathrm{CH}}_{4}+{\rm{H}}\to {\mathrm{CH}}_{3}+{{\rm{H}}}_{2}\\ &&{\mathrm{CH}}_{3}+{\mathrm{CH}}_{2}\to {{\rm{C}}}_{2}{{\rm{H}}}_{4}+{\rm{H}}\\ &&{{\rm{C}}}_{2}{{\rm{H}}}_{4}+{\rm{H}}\to {{\rm{C}}}_{2}{{\rm{H}}}_{3}+{{\rm{H}}}_{2}\\ &&{{\rm{C}}}_{2}{{\rm{H}}}_{3}+M\to {{\rm{C}}}_{2}{{\rm{H}}}_{2}+{\rm{H}}+M\\ &&{{\rm{C}}}_{2}{{\rm{H}}}_{2}+{\rm{O}}\to {\mathrm{CH}}_{2}+\mathrm{CO}\\ &&\mathrm{OH}+{\rm{H}}\to {\rm{O}}+{{\rm{H}}}_{2}\\ &&{\rm{H}}+{{\rm{H}}}_{2}{\rm{O}}\to \mathrm{OH}+{{\rm{H}}}_{2}\\ &&\displaystyle \frac{{{\rm{H}}}_{2}+M\to 2{\rm{H}}+M}{\mathrm{net}:{\mathrm{CH}}_{4}+{{\rm{H}}}_{2}{\rm{O}}\to \mathrm{CO}+3{{\rm{H}}}_{2}},\end{eqnarray*}$$

where ${{\rm{C}}}_{2}{{\rm{H}}}_{2}+{\rm{O}}$ $\to$ ${\mathrm{CH}}_{2}+\mathrm{CO}$ (R8) is the RLS in the above pathway. The corresponding timescale for C/O between 1 and 2 is

$$\begin{eqnarray}&&{\tau }_{{\mathrm{CH}}_{4}}=\displaystyle \frac{[{\mathrm{CH}}_{4}]}{{\rm{r}}1+\min ({\rm{r}}2,{\rm{r}}3)+{\rm{r}}9+\max ({\rm{r}}8,{\rm{r}}10)}.\end{eqnarray} \tag{ 10 }$$

In this scheme, CH₄–CO conversion is no longer limited by hydrogen dissociation.

We recognize that the chemical pathways taken depend entirely on the network one is using. With a different network, the rate-limiting chemical reactions need to be re-identified. The expressions for chemical timescales can be worked out following the same steps in Section 4. Here, we compare our network with two others that have included high-temperature chemical kinetics and been applied to hot exoplanets. We first compare our network with that of Moses et al. (2011); our two networks are naturally similar, because we have taken the values of certain rate coefficients from that study. We then compare our network with that of Venot et al. (2012). In Figure 8, we use our pathway analysis tool to compare the route taken by each network (since the reaction rates are calculated with equilibrium composition, the forward rate equals the reverse rate; there is no directionality in the pathways, i.e., CH₄ $\to$ CO takes the same pathway as CO $\to$ CH₄).

Zoom In Zoom Out Reset image size

At solar abundance and temperatures less than 2000 K, the network of Moses et al. (2011) shows almost exactly the same CH₄–CO pathways as ours. Due to a different choice of the rate coefficients associated with (R2), part of the parameter space occupied by (R3) is replaced by (R2) in their network for T ≳ 2000 K and P ∼ 10 bar. (R7) also extends to lower pressures in their network. At high temperatures and low pressures in scheme (A), their network experiences the same dehydrogenation process, but instead of (R4) and (R5), their network chooses a faster path through water: $\mathrm{CH}+{{\rm{H}}}_{2}{\rm{O}}$ $\to$H2CO + H (which is not included in our network). Yet, in this regime, the timescale is limited by hydrogen dissociation and yields the same timescale of CH₄–CO interconversion as ours does. At C/O = 2, carbon forms abundant ${{\rm{C}}}_{2}{{\rm{H}}}_{2}$ and then gets oxidized to CO in a similar way, but except via (R10), it takes ${{\rm{C}}}_{2}{{\rm{H}}}_{2}+\mathrm{OH}$ $\to$ ${{\rm{H}}}_{2}\mathrm{CCO}+{\rm{H}}$ or ${{\rm{C}}}_{2}{{\rm{H}}}_{2}+{\rm{O}}$ $\to$ HCCO + H in their network. ${{\rm{H}}}_{2}\mathrm{CCO}$ or HCCO then proceeds to be split into CO by H. In general, our network is consistent with that of Moses et al. (2011), as suggested by the comparison in Tsai et al. (2017).

The overall CH₄–CO timescale in Venot et al. (2012) is shorter than ours and that of Moses et al. (2011). We find the two key reactions that make the RLSs and essentially the timescales different are (R9) and

$$\begin{eqnarray}&&{\mathrm{CH}}_{2}\mathrm{OH}+{\rm{M}}\to {{\rm{H}}}_{2}\mathrm{CO}+{\rm{H}}+{\rm{M}}.\end{eqnarray} \tag{ R18 }$$

Venot et al. (2012) includes faster rate coefficients for both (R9) and (R18). Their rate coefficient for (R9) is based on the work of Hidaka et al. (1989), which has been suggested as overestimating the rate; see the discussion in Visscher (2012) and Moses (2014). This rate coefficient is significantly higher (∼10 orders of magnitude) than the ab initio calculation in Moses et al. (2011), which is also used in our network.

At lower temperatures (T ≲ 1000 K), Venot et al. (2012) takes the same pathway through (R9), forming ${\mathrm{CH}}_{3}\mathrm{OH}$ from ${\mathrm{CH}}_{3}$ and ${{\rm{H}}}_{2}{\rm{O}}$. However, with a faster rate it never bottlenecks the pathway and controls the timescale. Their RLSs are instead the reactions involving forming or destroying ${{\rm{H}}}_{2}\mathrm{CO}$, e.g., ${\mathrm{CH}}_{3}{\rm{O}}$ + M $\to$ ${{\rm{H}}}_{2}\mathrm{CO}+{\rm{H}}$ + M and ${{\rm{H}}}_{2}\mathrm{CO}+{\rm{H}}$ $\to$ $\mathrm{HCO}+{{\rm{H}}}_{2}$. Similarly, due to the faster ${\mathrm{CH}}_{3}$–${\mathrm{CH}}_{3}\mathrm{OH}$ channel via (R9), for 1000 K ≲ T ≲ 1500 K, Venot et al. (2012) exhibits pathways close to ours, except that our (R1) is replaced by ${\mathrm{CH}}_{2}\mathrm{OH}+{{\rm{H}}}_{2}$ $\to$ ${\mathrm{CH}}_{3}\mathrm{OH}+{\rm{H}}$ as the RLS. At higher temperature where T ≳ 1500 K, the differences are mainly attributed to (R18). Venot et al. (2012) uses the rate from Greenhill et al. (1986), validated for 600–1000 K, while this work and Moses et al. (2011) use the rate from Cribb et al. (1992), validated for 1900–2700 K. The former is about two orders of magnitude larger than the latter in this temperature range. Consquently, the fast ${\mathrm{CH}}_{2}\mathrm{OH}$–${{\rm{H}}}_{2}\mathrm{CO}$ interconversion in Venot et al. (2012) again never limits the pathway (e.g., the right pathway in Figure 8). Their RLS remains ${\mathrm{CH}}_{2}\mathrm{OH}+{{\rm{H}}}_{2}$ $\to$ ${\mathrm{CH}}_{3}\mathrm{OH}+{\rm{H}}$ for high pressures and switches to ${\mathrm{CH}}_{3}+\mathrm{OH}$ $\to$ ${\mathrm{CH}}_{2}\mathrm{OH}+{\rm{H}}$ for low pressures.

In conclusion, our pathway analysis tool is useful in identifying the key reactions for a given network, which allows us to diagnose the divergent behaviors of different networks. By isolating the rate coefficients of the key reactions involved, we hope to motivate future laboratory and/or theoretical studies that will resolve these discrepancies.

We are finally ready to validate the chemical relaxation method to a factor of several, having assembled the necessary ingredients.

6.1. Setup

The goal of any chemical calculation, either from full kinetics or any simplified method, is to provide the rate of change of every species locally. The applicability of the relaxation method should not depend on the format or complexity of atmospheric dynamics. Therefore, before including the chemical relaxation method in a GCM, we evaluate whether and to what extent the relaxation method can replace full chemical kinetics with a one-dimensional full kinetics model, where eddy diffusion (K_zz) is used to represent vertical mixing. For w (vertical velocity) ∼1 km s⁻¹ and H ∼ 100 km, we have ${K}_{{zz}}\lesssim {wH}\sim {10}^{12}$ cm² s⁻¹. To exaggerate the effects of vertical mixing, we explore a range of K_zz values up to 10¹⁵ cm² s⁻¹.

We run the same atmospheric conditions with the relaxation method and VULCAN, a full kinetics model with a C–H–O network including 29 species with up to two carbon atoms and about 300 forward and reverse reactions (Tsai et al. 2017). For nitrogen chemistry, an updated N–C–H–O network is implemented including 53 species and about 600 forward and reverse reactions. The chemical equilibrium abundances are calculated using the FastChem (Stock et al. 2018) code. We perform our chemistry calculations over three temperature–pressure profiles constructed using the analytical formula in Heng et al. (2014) to represent cool, warm, and hot atmospheres (Figure 9). These profiles are meant to mimic GJ 1214b-, HD 189733b-, and WASP-18b-like atmospheres. With the atmospheres we tested, the computational time for integrating one step using the chemical relaxation method was about 5 ms, compared to about 0.5 s using VULCAN. The computational speed can be increased by about 100 times using the chemical relaxation method.

Zoom In Zoom Out Reset image size

6.2. Hot Atmospheres (WASP-18b-like)

At solar metallicity and high temperatures, CO, H₂O, and N₂ are the dominant molecules in chemical equilibrium (Moses et al. 2011; Madhusudhan 2012; Heng & Tsai 2016). In Figure 10, it is therefore unsurprising that the mixing ratios of the first two molecules (not showing N₂) are insensitive to, or nearly independent of, pressure. For ${K}_{{zz}}\lesssim {10}^{13}$ cm² s⁻¹, the mixing ratios of CO and H₂O essentially track their chemical -equilibrium values closely, deviating only with stronger mixing and/or lower pressures (≲1 mbar).

Zoom In Zoom Out Reset image size

The mixing ratios of CH₄ and NH₃ exhibit a much larger range of values as they both drop off significantly with increasing altitude: 13 and 7 orders of magnitude for CH₄ and NH₃, respectively, over the range of pressure examined (0.1 mbar to 1 kbar). Over these broad ranges, the chemical relaxation method performs fairly well, exhibiting an accuracy of within an order of magnitude for the most part. With strong vertical mixing, CH₄ and NH₃ begin to quench at about 10 bar and 100 bar, respectively. At about 1 mbar, hydrogen dissociation/recombination slows down the interconversion and sets the second quench level. If the effect of hydrogen dissociation is neglected, the estimated chemical timescale will be too short and the prediction of CH₄ and NH₃ will be too close to chemical equilibrium.

6.3. Warm Atmospheres (HD 189733b-like)

Figure 11 shows the mixing ratios of CO, CO₂, CH₄, H₂O, and NH₃ for a HD 189733b-like atmosphere. The chemical relaxation method performs with an accuracy of better than a factor of 2 in most parts across a broad range of pressures and mixing ratio values, with ${\mathrm{NH}}_{3}$ being least accurate due to the error in estimating its timescale. For comparison, we show the chemical relaxation calculations performed using our implementation of the method of Cooper & Showman (2006) as dotted–dashed curves. Recall that Cooper & Showman use essentially a shorter single chemical timescale of CO (tied to a single RLS), whereas the main goal of the present study is to use a set of rate-limiting chemical reactions depending on the temperature and pressure conditions to obtain more accurate timescales.

Zoom In Zoom Out Reset image size

Cooper & Showman (2006) and Drummond et al. (2018) only calculate CO using the relaxation method and relate ${\mathrm{CH}}_{4}$ (and other species) through mass balance, assuming that all the carbon is locked in either ${\mathrm{CH}}_{4}$ and CO; hence the mixing ratios ${X}_{{\mathrm{CH}}_{4}}$ and X_CO are conserved:

$$\begin{eqnarray}&&{X}_{{\mathrm{CH}}_{4}}+{X}_{\mathrm{CO}}={\rm{C}}/({{\rm{H}}}_{2}+\mathrm{He}),\end{eqnarray} \tag{ 11 }$$

where ${\rm{C}}$/(${{\rm{H}}}_{2}$+$\mathrm{He}$) is the ratio of carbon atoms to molecular hydrogen and helium (i.e., the bulk gas), determined by solar metallicity. We emphasize that this mass balance relation is only valid when (1) the system is in or close to chemical equilibrium, (2) the temperature is not too high such that all hydrogen remains in molecular form, and (3) ${\mathrm{CH}}_{4}$ and CO have close abundances. The elemental abundance (${\rm{C}}$/(H₂+$\mathrm{He}$) in this case) is a local property, and can be violated by disequilibrium processes.⁶ We demonstrate this in the top right panel in Figure 11, with the mixing ratios of ${\mathrm{CH}}_{4}$ calculated from Equation (11) using X_CO obtained from chemical relaxation (the top left panel in Figure 11) as dotted–dashed curves. The mass balance approach overpredicts the quenching of ${\mathrm{CH}}_{4}$. Furthermore, as CO is around three orders of magnitude more abundant than ${\mathrm{CH}}_{4}$ in this case, which requires the estimation of CO to be as accurate as around three decimal places for the mass balance approach to work. Unfortunately, this precision is not attainable with the relaxation method or any other kinetics models. Therefore, the mass balance approach can only be applied to a system in chemical equilibrium, not to the relaxation method.

For CO₂, we show for comparison the calculation of relaxing CO₂ to the equilibrium abundance of CO2 (dotted–dashed curves) versus its pseudo-equilibrium value as determined by the quenched abundances of CO and H₂O (solid curves). The difference in accuracy is substantial: ∼10% versus an order of magnitude. Neglecting this effect will lead to the inaccurate prediction that the mixing ratio of CO₂ is close to equilibrium.

For illustration, we show the mixing ratios of CH₄ when the metallicity is 100× solar, as well as those of H₂O and CO₂ when C/O = 2. In the latter case, H₂O loses its dominance to CH₄ and its mixing ratio becomes sensitive to the strength of vertical mixing. In these cases, our more general treatment of chemical relaxation allows the accuracy of the method to remain about the same as for the solar metallicity case.

6.4. Cool Atmospheres (GJ 1214b-like)

Figure 12 shows the mixing ratios of CO, CO₂, CH₄, and N₂ for a GJ 1214b-like atmosphere. In this range of temperatures, the chemical relaxation method is highly accurate (∼10%). For CO, the difference from using the single-RLS timescale of Cooper & Showman (2006) increases significantly to a few orders of magnitudes. However, we note that photochemistry will potentially influence the abundances in these cooler atmospheres, which is not taken into account in this work.

Zoom In Zoom Out Reset image size

7.1. Summary

7.2. Opportunities for Future Work

S.-M.T. and K.H. acknowledge partial financial support from the PlanetS National Center of Competence in Research (NCCR), the Center for Space and Habitability, the Swiss National Science Foundation and the Swiss-based MERAC Foundation, as well as useful conversations with Julie Moses and Paul Rimmer.

Software: Python,⁷ SciPy,⁸ NumPy⁹ (van der Walt & Varoquaux 2011), Matplotlib¹⁰ (Hunter 2007).

[X] represents the number density of species X in chemical equilibrium and M refers to any third body.

For C/O < = 1:

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\tau }_{{\mathrm{CH}}_{4}}=\displaystyle \frac{[{\mathrm{CH}}_{4}]}{{k}_{1}[{\mathrm{CH}}_{3}][\mathrm{OH}][M]+\max (\min ({k}_{2}[{\mathrm{CH}}_{2}\mathrm{OH}][{\rm{H}}],{k}_{3}[{\mathrm{CH}}_{2}\mathrm{OH}][M]),{k}_{4}[{\mathrm{CH}}_{3}][{\rm{O}}]+{k}_{9}[{\mathrm{CH}}_{3}\mathrm{OH}][{\rm{H}}]+\min ({k}_{5}[{\rm{C}}][\mathrm{OH}],{k}_{6}[{\mathrm{CH}}_{4}][{\rm{H}}])}+{\tau }_{{{\rm{H}}}_{2}}\times \displaystyle \frac{3[\mathrm{CO}]}{[{{\rm{H}}}_{2}]}\\ {\tau }_{{{\rm{H}}}_{2}}=\displaystyle \frac{[{{\rm{H}}}_{2}]}{{k}_{{\rm{H}}}[{\rm{H}}][{\rm{H}}][M]}\end{array}\right.\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\tau }_{\mathrm{CO}}=\displaystyle \frac{[\mathrm{CO}]}{{k}_{1}[{\mathrm{CH}}_{3}][\mathrm{OH}][M]+\max (\min ({k}_{2}[{\mathrm{CH}}_{2}\mathrm{OH}][{\rm{H}}],{k}_{3}[{\mathrm{CH}}_{2}\mathrm{OH}][M]),{k}_{4}[{\mathrm{CH}}_{3}][{\rm{O}}]+{k}_{9}[{\mathrm{CH}}_{3}\mathrm{OH}][{\rm{H}}]+\min ({k}_{5}[{\rm{C}}][\mathrm{OH}],{k}_{6}[{\mathrm{CH}}_{4}][{\rm{H}}])}+{\tau }_{{{\rm{H}}}_{2}}\times \displaystyle \frac{3[\mathrm{CO}]}{[{{\rm{H}}}_{2}]}.\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\tau }_{{{\rm{H}}}_{2}{\rm{O}}}=\displaystyle \frac{[{{\rm{H}}}_{2}{\rm{O}}]}{{k}_{1}[{\mathrm{CH}}_{3}][\mathrm{OH}][M]+\max (\min ({k}_{2}[{\mathrm{CH}}_{2}\mathrm{OH}][{\rm{H}}],{k}_{3}[{\mathrm{CH}}_{2}\mathrm{OH}][M]),{k}_{4}[{\mathrm{CH}}_{3}][{\rm{O}}]+{k}_{9}[{\mathrm{CH}}_{3}\mathrm{OH}][{\rm{H}}]+\min ({k}_{5}[{\rm{C}}][\mathrm{OH}],{k}_{6}[{\mathrm{CH}}_{4}][{\rm{H}}])}+{\tau }_{{{\rm{H}}}_{2}}\times \displaystyle \frac{3[\mathrm{CO}]}{[{{\rm{H}}}_{2}]}\end{eqnarray} \tag{ 14 }$$

For C/O > 1:

$$\begin{eqnarray}&&{\tau }_{{\mathrm{CH}}_{4}}=\displaystyle \frac{[{\mathrm{CH}}_{4}]}{{k}_{1}[{\mathrm{CH}}_{3}][\mathrm{OH}][M]+\min ({k}_{2}[{\mathrm{CH}}_{2}\mathrm{OH}][{\rm{H}}],{k}_{3}[{\mathrm{CH}}_{2}\mathrm{OH}][M])+{k}_{9}[{\mathrm{CH}}_{3}\mathrm{OH}][{\rm{H}}]+\max ({k}_{8}[{{\rm{C}}}_{2}{{\rm{H}}}_{2}][{\rm{O}}],{k}_{10}[{{\rm{C}}}_{2}{{\rm{H}}}_{2}][\mathrm{OH}])}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\tau }_{\mathrm{CO}}=\displaystyle \frac{[\mathrm{CO}]}{{k}_{1}[{\mathrm{CH}}_{3}][\mathrm{OH}][M]+\min ({k}_{2}[{\mathrm{CH}}_{2}\mathrm{OH}][{\rm{H}}],{k}_{3}[{\mathrm{CH}}_{2}\mathrm{OH}][M])+{k}_{9}[{\mathrm{CH}}_{3}\mathrm{OH}][{\rm{H}}]+\max ({k}_{8}[{{\rm{C}}}_{2}{{\rm{H}}}_{2}][{\rm{O}}],{k}_{10}[{{\rm{C}}}_{2}{{\rm{H}}}_{2}][\mathrm{OH}])}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\tau }_{{{\rm{H}}}_{2}{\rm{O}}}=\displaystyle \frac{[{{\rm{H}}}_{2}{\rm{O}}]}{{k}_{1}[{\mathrm{CH}}_{3}][\mathrm{OH}][M]+\min ({k}_{2}[{\mathrm{CH}}_{2}\mathrm{OH}][{\rm{H}}],{k}_{3}[{\mathrm{CH}}_{2}\mathrm{OH}][M])+{k}_{9}[{\mathrm{CH}}_{3}\mathrm{OH}][{\rm{H}}]+\max ({k}_{8}[{{\rm{C}}}_{2}{{\rm{H}}}_{2}][{\rm{O}}],{k}_{10}[{{\rm{C}}}_{2}{{\rm{H}}}_{2}][\mathrm{OH}])}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}{[{\mathrm{CO}}_{2}]}_{\mathrm{pseudo}-\mathrm{eq}} & = & \displaystyle \frac{[\mathrm{CO}][{{\rm{H}}}_{2}{\rm{O}}]{[{\rm{H}}]}_{\mathrm{eq}}}{{[\mathrm{CO}]}_{\mathrm{eq}}{[{{\rm{H}}}_{2}{\rm{O}}]}_{\mathrm{eq}}[{\rm{H}}]}\\ {\tau }_{{\mathrm{CO}}_{2}} & = & \displaystyle \frac{[{\mathrm{CO}}_{2}]}{{k}_{{\mathrm{CO}}_{2}}[\mathrm{CO}][\mathrm{OH}]}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{\tau }_{{\mathrm{NH}}_{3}}=\displaystyle \frac{1}{2}\left(\displaystyle \frac{[{\mathrm{NH}}_{3}]}{\max ({k}_{11}[{\mathrm{NH}}_{2}][{\mathrm{NH}}_{2}],{k}_{12}[{{\rm{N}}}_{2}{{\rm{H}}}_{3}][M])+{k}_{13}[\mathrm{NH}][{\mathrm{NH}}_{2}]+{k}_{14}[\mathrm{NO}][{\mathrm{NH}}_{2}])+{k}_{15}[{\rm{N}}][\mathrm{NO}]}+{\tau }_{{{\rm{H}}}_{2}}\times \displaystyle \frac{3[{{\rm{N}}}_{2}]}{[{{\rm{H}}}_{2}]}\right)\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\tau }_{{{\rm{N}}}_{2}}=\displaystyle \frac{[{{\rm{N}}}_{2}]}{\max ({k}_{11}[{\mathrm{NH}}_{2}][{\mathrm{NH}}_{2}],{k}_{12}[{{\rm{N}}}_{2}{{\rm{H}}}_{3}][M])+{k}_{13}[\mathrm{NH}][{\mathrm{NH}}_{2}]+{k}_{14}[\mathrm{NO}][{\mathrm{NH}}_{2}])+{k}_{15}[{\rm{N}}][\mathrm{NO}]}+{\tau }_{{{\rm{H}}}_{2}}\times \displaystyle \frac{3[{{\rm{N}}}_{2}]}{[{{\rm{H}}}_{2}]}.\end{eqnarray} \tag{ 20 }$$

Finding the chemical pathway is similar to the path finding problem in graph theory. That is, all species in the network are presented by nodes and reactions between them form the edges, weighted by the reaction rates (faster reactions form shorter connections). This is equivalent to finding the least time-consuming route from a a starting node to an end node in the network and the total time-cost of the route can be approximated by the slowest edge (the RLS). While there are several different algorithms, Lehmann (2002) for example identifies the dominant pathways that are most efficient in removing/producing a species of interest from kinetics results with temporal evolution. We are particularly interested in finding the pathway in chemical equilibrium, which fulfils the need for estimating the chemical timescales using equilibrium abundances. For this purpose, We implement Dijkstra's algorithm (Dijkstra 1959), which is easy to implement and highly efficient in finding the fastest paths in network problems (Viswanath et al. 2013).

Table 4.  Rate-limiting Reactions/Steps and Their Rate Coefficients

Index	Reaction	Rate Coefficient	Reference
R1	$\mathrm{OH}+{\mathrm{CH}}_{3}\mathop{\longrightarrow }\limits^{{\rm{M}}}{\mathrm{CH}}_{3}\mathrm{OH}$	k0 = 1.93 × 103 T−9.88 exp(−7544/T) + 5.11 × 10−11 T−6.25 exp(−1433/T) k∞ = 1.03 × 10−10 T−0.018 exp(16.74/T)	Moses et al. (2011)
R2rev	${\mathrm{CH}}_{2}\mathrm{OH}+{\rm{H}}\to \mathrm{OH}+{\mathrm{CH}}_{3}$	1.60 × 10−10	NIST (shown by Squib in the NIST database) 1987TSA471
R3	${\mathrm{CH}}_{2}\mathrm{OH}\mathop{\longrightarrow }\limits^{{\rm{M}}}{\rm{H}}+{{\rm{H}}}_{2}\mathrm{CO}$	k0 = 1.66 × 10−10 exp(−12630/T) k∞ = 3 × 109 exp(−14600/T)	NIST 1987TSA471 NIST 1975BOW343
R4	${\mathrm{CH}}_{3}+{\rm{O}}\to {{\rm{H}}}_{2}\mathrm{CO}+{\rm{H}}$	1.4 × 10−10	NIST 1992BAU/COB411-429
R5	$\mathrm{OH}+{\rm{C}}\to \mathrm{CO}+{\rm{H}}$	1.05 × 10−12 T0.5	NSRDS 67
R6	${\rm{H}}+{\mathrm{CH}}_{4}\to {\mathrm{CH}}_{3}+{{\rm{H}}}_{2}$	2.20 × 10−20 T3 exp(−4040/T)	NIST 1992BAU/COB411-429
R7	${\mathrm{CH}}_{3}\mathrm{OH}+{\rm{H}}\to {\mathrm{CH}}_{3}{\rm{O}}+{{\rm{H}}}_{2}$	6.82 × 10−20 T2.685 exp(−4643/T)	NIST 1984WAR197C
R8	${{\rm{C}}}_{2}{{\rm{H}}}_{2}+{\rm{O}}\to {\mathrm{CH}}_{2}+\mathrm{CO}$	6.78 × 10−16 T1.5 exp(−854/T)	NIST 1987CVE261
R9rev	${\mathrm{CH}}_{3}\mathrm{OH}+{\rm{H}}\to {\mathrm{CH}}_{3}+{{\rm{H}}}_{2}{\rm{O}}$	4.91 × 10−19 T2.485 exp(−10380/T)	Moses et al. (2011)
R10	${{\rm{C}}}_{2}{{\rm{H}}}_{2}+\mathrm{OH}\to {\mathrm{CH}}_{3}+\mathrm{CO}$	8.04 × 10−28 T4 exp(1010/T)	NIST 1989MIL/MEL1031-1039
R11	${\mathrm{NH}}_{2}+{\mathrm{NH}}_{2}\to {{\rm{N}}}_{2}{{\rm{H}}}_{2}+{{\rm{H}}}_{2}$	2.89 × 10−16 T1.02 exp(−5930/T)	NIST 2009KLI/HAR10241-10259
R12	${{\rm{N}}}_{2}{{\rm{H}}}_{3}+\mathop{\longrightarrow }\limits^{{\rm{M}}}{{\rm{N}}}_{2}{{\rm{H}}}_{2}+{\rm{H}}$	k0 = 3.49 × 1038 T−13.13 exp(−36825/T) k∞ = 7.95 × 1013 exp(−27463/T)	Hwang & Mebel (2003)
R13	$\mathrm{NH}+{\mathrm{NH}}_{2}\to {{\rm{N}}}_{2}{{\rm{H}}}_{2}+{\rm{H}}$	6.98 × 10−10 T−0.27 exp(39/T)	NIST 2009KLI/HAR10241-10259
R14	$\mathrm{NO}+{\mathrm{NH}}_{2}\to {{\rm{N}}}_{2}+{{\rm{H}}}_{2}{\rm{O}}$	7.9 × 10−9 T−1.1 exp(−98/T)	NIST 1994DIA/YU4034-4042
R15	${\rm{N}}+\mathrm{NO}\to {{\rm{N}}}_{2}+{\rm{O}}$	3.7 × 10−11	NIST 1992MIC/LIM3228-3234
R16	${{\rm{N}}}_{2}{{\rm{H}}}_{4}+{\rm{H}}\to {{\rm{N}}}_{2}{{\rm{H}}}_{3}+{{\rm{H}}}_{2}$	1.17 × 10−11 exp(−1260/T)	NIST 1995VAG777-790
kH	${\rm{H}}+{\rm{H}}\mathop{\longrightarrow }\limits^{{\rm{M}}}{{\rm{H}}}_{2}$	k0 = 2.7 × 10−31 T−0.6 k∞ = 3.31 × 10−6 T−1	NIST 1992BAU/COB411-429 NIST 1965JAC/GIE3688
${k}_{{\mathrm{CO}}_{2}}$	$\mathrm{CO}+\mathrm{OH}\to {\mathrm{CO}}_{2}+{\rm{H}}$	1.05 × 10−17 T1.5 exp(259/T)	NIST 1992BAU/COB411-429

Note. The notation k₀ and k_∞ are the low-pressure and high-pressure limiting rate coefficients, respectively. The rate coefficient for the termolecular or thermal dissociation reaction is expressed by Equation (28) in Tsai et al. (2017). Rev indicates the reverse reactions of Rx, which are used to derive the rate constants of Rx using the equilibrium constant calculated by the NASA polynomials (http://garfield.chem.elte.hu/Burcat/burcat.html), as described in Appendix E of Tsai et al. (2017).

Download table as:  ASCIITypeset image

Examples of comparing different chemical networks by identifying the pathways with associated RLSs are shown in Figures 7 and 13.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Below are the steps used in Dijkstra's algorithm to find the shortest path.

• 1.  
  Create a list of visited nodes (initially empty). Assign tentative distance values to all nodes: set to zero for the initial node and to infinity for all other nodes. Set the initial node as the current node.
• 2.  
  From the current node, consider all of its neighbors and calculate their tentative distances. Update the distance if the new value is smaller than the previously assigned value.
• 3.  
  Include the current node in the visited-node list. A visited node will never be checked again.
• 4.  
  Stop if the destination node has been marked visited. The pathway has been found and the longest edge (the slowest step) in the path is the RLS. Otherwise, select the unvisited node with the smallest tentative distance, set it as the new "current node," and go back to step 2.

We demonstrate the steps in the following example:

Consider the simple network with nodes a–f in Figure 14. We will try to find the shortest path between a and f.

• 1.  
  The tentative distance values are assigned to a, b, c, d, e, f as [0, inf, inf, inf, inf, inf]. The current node is a.
• 2.  
  The adjacent nodes are b, d, and e. The distance values are updated to [0, 2, inf, 1, 6, inf].
• 3.  
  The list of visited nodes is now updated to [a]. d as the unvisited node with the smallest distance value becomes the current node.
• 4.  
  The adjacent nodes are c and e. The distance values are updated to [0, 2, 10, 1, 6, inf].
• 5.  
  The list of visited nodes is now updated to [a, d]. b as the unvisited node with the smallest distance value becomes the current node.
• 6.  
  We repeat the above steps (second to fifth steps in the algorithm) until f has been included in the list of visited nodes. We then have the distance values as [0, 2, 5, 1, 6, 9]. The shortest path is obtained by going back from f and following the smallest distance value. The path is a $\to$ b $\to$ c $\to$ f.