A chemical kinetics code for modelling exoplanetary atmospheres (original) (raw)

Journal Article

,

Institute of Astronomy, University of Cambridge, Cambridge CB3 0HA, UK

Search for other works by this author on:

,

Institute of Astronomy, University of Cambridge, Cambridge CB3 0HA, UK

Cambridge Earth Sciences, University of Cambridge, Cambridge CB2 3EQ, UK

Search for other works by this author on:

,

Institute of Astronomy, University of Cambridge, Cambridge CB3 0HA, UK

Search for other works by this author on:

Cambridge Earth Sciences, University of Cambridge, Cambridge CB2 3EQ, UK

MRC Laboratory of Molecular Biology, Cambridge CB2 OQH, UK

Cavendish Astrophysics, University of Cambridge, Cambridge CB3 0HE, UK

Search for other works by this author on:

Received:

25 January 2019

Revision received:

03 May 2019

• PDF
• Split View
• • Article contents
  • Figures & tables
  • Video
  • Audio
  • Supplementary Data
• Cite

Cite

Richard Hobbs, Oliver Shorttle, Nikku Madhusudhan, Paul Rimmer, A chemical kinetics code for modelling exoplanetary atmospheres, Monthly Notices of the Royal Astronomical Society, Volume 487, Issue 2, August 2019, Pages 2242–2261, https://doi.org/10.1093/mnras/stz1333
Close

Navbar Search Filter Mobile Enter search term Search

Abstract

Chemical compositions of exoplanets can provide key insights into their physical processes and formation and evolutionary histories. Atmospheric spectroscopy provides a direct avenue to probe exoplanetary compositions. However, whether obtained in transit or thermal emission, spectroscopic observations probe limited pressure windows of planetary atmospheres and are directly sensitive to only a limited set of spectroscopically active species. It is therefore critical to have chemical models that can relate retrieved atmospheric compositions to an atmosphere’s bulk physical and chemical state. To this end, we present Levi a new chemical kinetics code for modelling exoplanetary atmospheres. Levi calculates the gas phase hydrogen, oxygen, carbon, and nitrogen chemistry in planetary atmospheres. Here, we focus on hot gas giants. Applying Levi, we investigate how variations in bulk C/O and N/O affects the observable atmospheric chemistry in hot Jupiters. For typical hot Jupiters, we demonstrate the strong sensitivity of molecular detections to the atmospheric C/O. Molecular detections are conversely less sensitive to the atmospheric N/O ratio, although highly supersolar N/O can decrease the C/O required for HCN and NH3 detection. Using a new pressure–temperature (P–T) profile for HD 209458b without a thermal inversion, we evaluate recently reported detection’s of CO, H2O, and HCN in its day-side atmosphere. We find that our models are consistent with the detected species, albeit with a narrow compositional window around C/O ∼ 1. A C/O ≳ 0.9 (1.6 times solar) was required to meet the minimum reported value for HCN, while a C/O ≲ 1 (1.8 times solar) was required to fit the nominal H2O abundance.

1 INTRODUCTION

Identifying and characterizing features of exoplanets is the observational cornerstone to answering some of the biggest questions in planetary astrophysics. Within the field of giant planets, major questions remain unresolved regarding their formation, composition, evolution during atmospheric escape, and how significant migration is in their history (Beaugé & Nesvorný 2012; Mordasini et al. 2016; Nelson, Ford & Rasio 2017; Dawson & Johnson 2018; McGuire et al. 2018). Both the location of where giant planets form and possible later migration through the protoplanetary disk may leave signs in the compositions of their atmospheres (e.g. Öberg, Murray-Clay & Bergin 2011; Moses et al. 2013; Madhusudhan, Amin & Kennedy 2014a; Mordasini et al. 2016). Tracing giant planet migration history is also important for understanding terrestrial planet formation and survivability. This question is clearest for hot Jupiter systems, where the presumed inward migration of a Jupiter mass object may have significantly disrupted habitable zone planets (Del Popolo, Ercan & Yeşilyurt 2005; Hands & Alexander 2016; Sánchez, de Elía & Darriba 2018).

A key tool for characterizing an exoplanet’s atmosphere is its spectrum, observed via transits or thermal emission. The molecules in a planet’s atmosphere absorb and emit light at particular frequencies, which can imprint a signature of the planet’s atmosphere on the light received at Earth from the planet–star system. Thus, it is possible to infer atmospheric chemistry based upon the features of the planet’s transmission or emission spectrum (Seager & Sasselov 2000; Brown 2001; Deming et al. 2005; Grillmair et al. 2008; Seager & Deming 2010; Madhusudhan et al. 2016).

However, spectra only contain emission or absorption features for the spectrally active species present in a planet’s atmosphere, and therefore provide only a partial view of atmospheric chemistry and dynamics. In addition, the observed transmission or emission spectrum comes from a narrow band of pressure ranges, generally between 1 and 10−3 bar, thus lacking sensitivity to chemistry in the rest of the atmosphere. Yet, with vertical mixing and diffusion transporting species through the atmosphere, chemistry at higher and lower pressures can contribute to atmospheric species in the observed pressure window. It is also possible for layers of clouds or hazes to block observations, resulting in muted spectral features that can imply incorrect abundances. This means that models are required to complete the picture of physical and chemical processes occurring in a planet.

Atmospheric models can fill in details of atmospheric composition and structure that retrievals are not directly sensitive to, enabling physically and chemically informed interpolation and extrapolation of the retrieval observations. Improving the estimates of bulk atmospheric composition from retrievals is also then possible, by comparing the model atmosphere at different C/O and N/O ratios to the retrieved atmospheric abundances. Since the composition of planetary atmospheres may reflect when and where in the prototplanetary disk a planet formed (Öberg et al. 2011; Madhusudhan et al. 2014a; Cleeves et al. 2018), these elemental ratios may ultimately provide a key insight into the history of planet formation.

To construct the chemical model, a network of atmospheric chemical reactions is compiled and used as input to a model that calculates the production and loss of all species in the atmosphere. This produces a profile of chemical abundances for an atmosphere in thermochemical equilibrium. However, a planetary atmosphere may not be in chemical equilibrium over its entire extent. It is expected that the deeper, hotter parts of the atmosphere (P ≳ 0.1–1 bar) approach thermochemical equilibrium; however, higher up (at lower pressures) various non-equilibrium processes can become significant. Therefore, to more accurately model atmospheric chemistry, it is necessary to include the non-equilibrium processes that can dominate at lower pressures, such as vertical mixing and photochemistry (Cooper & Showman 2006; Moses et al. 2011, 2013).

We report a new code for one-dimensional (1D) modelling of chemistry in exoplanetary atmospheres, named Levi.1 The chemical processes that it can model include thermochemical equilibrium, eddy-diffusion, molecular diffusion, thermal diffusion, and photochemistry. Through inputs in the form of pressure–temperature (P–T) profiles, vertical mixing profiles, stellar properties, and a network of possible chemical reactions, the model can predict the mixing ratios of chemical species in a planet’s atmosphere. Currently, the code uses a network that includes molecules comprised of H, C, O, N, or He. The code is built to be modular in the chemical processes and reactions it can model, and the species it includes. We use this functionality to better understand the effects any of these components of the model can have on the chemistry of the atmosphere. Here, we focus on applying the code to the hydrogen-rich atmospheres of hot Jupiters, although it can be generalized to apply to the atmospheres of terrestrial planets.

Several models of exoplanetary atmospheric chemistry have been previously reported in the literature. While some focus purely on thermochemical equilibrium (Burrows & Sharp 1999; Lodders & Fegley 2002; Blecic, Harrington & Bowman 2016), others have also sought to encapsulate the non-equilibrium chemistry occurring (i.e. Liang et al. 2003; Zahnle et al. 2009; Moses et al. 2011; Hu et al. 2012; Miller-Ricci Kempton, Zahnle & Fortney 2012; Venot et al. 2012; Rimmer & Helling 2016; Tsai et al. 2017). Table 1 shows a comparison of the processes included in a selection of chemical kinetics models. All these models work similarly, using a P–T and eddy-diffusion (Kzz) profile as inputs to their codes to calculate abundance profiles of the atmosphere, with most of the differences between their predictions arising in their chemical network. Liang et al. (2003) first reported a photochemical model with H–C–O chemistry for a select set of species in highly irradiated atmospheres. Zahnle et al. (2009) investigated the effects of photochemistry on sulphur bearing molecules in hot Jupiters, and later warm Jupiters (Zahnle et al. 2016). This model was also adopted and furthered in Miller-Ricci Kempton et al. (2012). Moses et al. (2011) and Venot et al. (2012) both produced codes that focused on H–C–N–O photochemistry on hot Jupiters. Agúndez et al. (2014) adapted the work of Venot et al. (2012) to report a two-dimensional model for hot Jupiters. Hu et al. (2012) created a model for the atmosphere of super-Earths, as well as incorporating the additional non-equilibrium effects that could occur on terrestrial planets. Rimmer & Helling (2016) focused on the effects that ions produced by lightning could have on the chemistry of hot Jupiters, as well as including a selection of metals in their network. Drummond et al. (2016) reported a model that self-consistently solved for both the P–T profile and chemical abundance profiles. Tsai et al. (2017) use a reduced chemical network to reproduce the results of Moses et al. (2011) and Rimmer & Helling (2016) for several hot Jupiters, and analyse the accuracy of the quenching approximation. Ultimately, all of these approaches, including the one presented in this work, represent a compromise between efficiency and complexity, with full 3D coupled chemical and dynamical simulations being required to capture all aspects of atmospheric chemistry.

Table 1.

A comparison of the main features of a selection of non-equilibrium atmospheric chemistry codes in the literature.

_a_These rates were measured for a low temperature Jovian Model.

Table 1.

A comparison of the main features of a selection of non-equilibrium atmospheric chemistry codes in the literature.

_a_These rates were measured for a low temperature Jovian Model.

We demonstrate Levi by applying it to two of the most studied hot Jupiters, HD 209458b and HD 189733b. HD 209458b and HD 189733b are excellent candidates for benchmarking Levi, as they are some of the best characterized planets (e.g. Seager & Sasselov 2000; Moses et al. 2013; Liu & Showman 2013; Mayne et al. 2014; Schwarz et al. 2015; Amundsen et al. 2016; Brewer, Fischer & Madhusudhan 2017), and several previous chemical kinetic models have been applied to them. Between the two planets, a preliminary parameter space can be explored, comparing the differences between hot Jupiters with and without thermal inversions, the influence of strong UV radiation, and the effects of an increase in a planet’s effective temperature. In particular, HD 189733b is a very good candidate to study the effect of photochemistry, because of both the high stellar UV it receives and the relatively low atmospheric temperature. HD 209458b was thought to have a temperature inversion (Knutson et al. 2008), although that is now debated (Diamond-Lowe et al. 2014; Schwarz et al. 2015). The temperature inversion is being kept in most of this work, both to match Moses et al. (2011) for benchmarking and to better investigate the effects of high temperature at high and low pressures. We do, however, also use a P–T profile for HD 209458b that lacks a thermal inversion when comparing our model to a recent work that provides evidence of the presence of number of molecular species in its atmosphere (Hawker et al. 2018).

In Section 2, the model and network being used for this code are described. Section 3 contains a validation against previously constructed codes and an examination of how the abundances of species in the atmosphere vary from thermochemical equilibrium due to vertical transport and photochemistry. Section 4 contains an exploration of the C/O and N/O parameter space. This is used to investigate how the bulk chemistry of an atmosphere can affect the mixing ratios of species within it, and to compare our model’s predictions to evidence of species detected on HD 209458b. Section 5 contains a short summary and discussion of the results of this work.

2 THE ATMOSPHERIC MODEL

Our model considers a 1D plane-parallel atmosphere in hydrostatic equilibrium, with various processes affecting the chemistry. Using a given P–T profile, considerations of hydrostatic equilibrium, and the ideal gas law, the thermodynamic parameters of the atmosphere are determined. In this manner, we follow the standard approach of existing chemical models (see Table 1) in using predetermined profiles of both the temperature and vertical mixing strength, without self-consistently solving the radiative transfer equation.

We model chemical diffusion, photochemistry, and, as a limiting case, chemical equilibrium. Thermochemical equilibrium occurs once the rate of production and loss rate of every molecule due to chemical reactions is the same. Including diffusion, in the form of eddy, molecular, or thermal diffusion, can lead to species being lifted from deep in the atmosphere into higher, cooler layers, or can cause large, heavy molecules to sink out of the upper atmosphere. Species that get lifted to higher atmosphere layers can have their abundance frozen in above thermochemical equilibrium values, since the time it takes for them to deplete back down to their equilibrium abundance is significantly longer than their resupply time from vertically mixing. This process is known as quenching. The eddy diffusion used is an approximation for atmospheric motion, so that the complex micro-physics of atmospheric mixing can be treated as a simple macroscopic process. In the upper atmosphere, high energy photons can cause some molecules to disassociate, and new species to form from the breakdown products. The rate of these reactions is determined by the UV flux received by the planet.

Levi solves a set of coupled differential equations that govern the chemical kinetics (described in Section 2.1.1) to obtain a steady-state solution for the chemical composition of the atmosphere. By dividing up the atmosphere into vertically distributed plane-parallel layers, these equations can be described in just one dimension. Each equation describes the time evolution of one species in the atmosphere due to production or loss from chemical interactions, or fluxes from other atmosphere layers. Other factors that could influence the evolution of the atmosphere, such as evaporation, atmospheric escape, and geological outgassing are not currently implemented, but, in the context of giant planets, these processes are in general not operating and/or are not chemically significant (Moses et al. 2011).

A steady-state solution to these equations is found by time-stepping the code using a second-order Rosenbrock solver (Rosenbrock 1963), from an initial guess for the chemical abundances, until the production and loss rates of every species in every atmosphere layer are balanced. For a closed system, this will produce a solution that is stable over large time periods.

The code is designed to be highly configurable, allowing modelling of a wide range of exoplanets. This includes atmospheric parameters such as elemental composition, the P–T profile of the atmosphere, and the eddy diffusion coefficient. The stellar characteristics, such as the stellar spectrum and the planet–star separation, and the physical characteristics of the planet such as its surface gravity can also be changed as required to model different examples of hot Jupiters. The gravitational acceleration used throughout the atmosphere is an approximation of the average gravity, taken at the surface of the planet. For Jupiter-like planets, the ‘surface’ is taken to be at |$\mathrm{1\, bar}$|⁠. It is also possible to selectively exclude a subset of reactions, such as those that consider the effects of photochemistry, or the reactions of ions, to better understand the effects those reactions can have on the atmosphere. In addition, any reactions containing an undesired chemical species can be easily stripped away, to both explore the atmospheric response to the presence of particular species and to improve the computational runtime.

2.1 Model set-up

Here, we describe the equations that govern the model being used, as well as the chosen boundary and initial conditions for the model.

2.1.1 Chemical processes

The model used considers contribution from multiple processes, which can be described by a series of coupled 1D continuity equations,

\begin{eqnarray*} \frac{\partial n_{i}}{\partial t} = \mathcal {P}_{i} - \mathcal {L}_{i} - \frac{\partial \Phi _{i}}{\partial z}, \end{eqnarray*}

(1)

where ni (m−3) is the number density of species i, with i = 1, ..., Ni, with Ni being the total number of species. |$\mathcal {P}_{i}$| (m−3 s−1) and |$\mathcal {L}_{i}$| (m−3 s−1) are the production and loss rates of the species i. ∂t (s) and ∂z (m) are the infinitesimal time-step and altitude step, respectively. Φ_i_ (m−2 s−1) is the upward vertical flux of the species, given by

\begin{eqnarray*} \Phi _{i} = -(K_{zz}+D_i)n_{t}\frac{\partial X_{i}}{\partial z} + D_i n_i\left(\frac{1}{H_0} - \frac{1}{H} - \frac{\alpha _{T,i}}{T}\frac{{\rm d}T}{{\rm d}z}\right), \end{eqnarray*}

(2)

where Xi and nt (m−3) are the mixing ratio and total number density of molecules such that ni = Xi nt. The eddy-diffusion coefficient, Kzz (m2 s−1), approximates the rate of vertical transport and Di (m2 s−1) is the molecular diffusion coefficient of species i. H_0 (m) is the mean scale height, H (m) is the molecular scale height, T (K) is the temperature, and α_T, i is the thermal diffusion factor. The equation for the molecular diffusion coefficient is adopted from Chapman & Cowling (1970), where it is defined as

\begin{eqnarray*} D_i = \frac{3}{8n_{t}(\frac{1}{2}[d_i + d_t])^2}\left(\frac{k_bT(m_i+m_t)}{2\pi m_i m_t}\right)^{1/2}, \end{eqnarray*}

(3)

where mi (kg) and di (m) are the mass and diameter of species i, respectively, and mt (kg) and dt (m) are the average mass and average diameter of all species in the atmosphere, respectively. The thermal diffusion factor is estimated from experimental and theoretical data from Chapman & Cowling (1970): for atomic hydrogen, α_T, i_ ≈ −0.1(1 − ni/nt); for He,α_T, i_ ≈ 0.145(1 − ni/nt); for H2, α_T, i_ ≈ −0.38; and for all other molecules α_T, i_ ≈ 0.25. Under the assumption that the atmosphere can be approximated to an ideal gas, the total number density, nt (m−3), is

\begin{eqnarray*} n_{t} = \frac{P(z)}{k_b\ T(z)}, \end{eqnarray*}

(4)

where P (Pa) and T (K) are the altitude dependent pressure and temperature, and _k_b is the Boltzmann constant.

The production and loss rates of species i can be calculated by considering all the reactions that include said species, (e.g. |$n_i {+} n_{i^{\prime }} \rightleftharpoons n_{i^{\prime \prime }} {+} n_{i^{\prime \prime \prime }}$|⁠), then finding the rate at which ni is produced (⁠|$k_2 n_{i^{\prime \prime }} n_{i^{\prime \prime \prime }}$|⁠) and lost (⁠|$k_1 n_i n_{i^{\prime }}$|⁠) in each reaction, where _k_1 and _k_2 are the rate of the forward and reverse reactions, respectively, before summing over all such reactions such that

\begin{eqnarray*} \mathcal {P}_{i} - \mathcal {L}_{i} = \sum (k_2 n_{i^{\prime \prime \prime }} n_{i^{\prime \prime }} - k_1 n_i n_{i^{\prime }}). \end{eqnarray*}

(5)

2.1.2 Atmospheric profile

The physical characteristics of the atmospheres Levi simulates have a profound impact on their chemistry. This is realized through a P–T profile that describes the atmosphere being investigated. Currently, the P–T profiles are not calculated self-consistently, but instead use pre-calculated profiles as an input. This simplification is typical for chemical kinetics models, since implementing radiative transfer to allow iterative recalculation of the P–T profile would significantly increase the runtime, in addition to the possibility of converging to a non-physical P–T profile and atmosphere chemistry. There have been other works investigating how having a self-consistent P–T profile may affect the results of chemical kinetics models. Drummond et al. (2016) found that for strong disequilibrium chemistry, self-consistent P–T profiles can vary by up to 100K. For the two test cases that we study here, the hot Jupiters HD 189733b and HD 209458b, the P–T profiles are given in Fig. 1(a). We also consider an isothermal profile for reference.

Figure 1.

Plot (a) are the P–T profiles for the average day-side of HD 209458b and HD 189733b assuming a solar composition atmosphere, adopted from Moses et al. (2011). At pressures below 10−7 bar, an approximation for the thermosphere is used from Garcia Munoz (2007). Plot (b) displays the Kzz profiles that are being used for each planet, again drawn from Moses et al. (2011).

The P–T and Kzz profiles for HD 209458b and HD 189733b that are used in most of this work have been adopted from Moses et al. (2011). The deepest layer of the atmosphere is chosen to be |$\mathrm{10^{3}\, bar}$|⁠, a pressure at which all species are expected to be in chemical equilibrium, while the highest atmosphere layer (⁠|$\mathrm{10^{-8}\, bar}$|⁠) should be sufficient to cover the photochemical destruction and production region for all relevant molecules. For pressures below 10−7 bars, an approximation for the thermosphere is used based on the work of Garcia Munoz (2007).

2.1.3 Actinic flux

Many species in the atmosphere, upon absorbing a sufficiently energetic photon, can be photodissociated into a variety of species, including radicals. Since the rate of photodissociation is proportional to the number density of UV photons, it is necessary to define the actinic flux, the spherically integrated flux in any atmosphere layer due to irradiation from the host star. This can be calculated as

\begin{eqnarray*} F(\lambda ,z) = F_0(\lambda)\ {\rm e}^{-\tau (\lambda ,z)/\mu _0}, \end{eqnarray*}

(6)

where _F_0 (photons m−2 s−1 nm−1) is the flux at the top of the atmosphere (where by definition τ = 0), μ0 is angle of the path of sunlight, chosen to be 57.3° (Hu et al. 2012). τ is the total optical depth, as a function of wavelength, integrated over all the atmosphere that the light passes through, defined as

\begin{eqnarray*} \tau (\lambda ,z) = \int _{z}^{\infty } n_i(z^{\prime }) \sigma _i(\lambda) {\rm d}z^{\prime } , \end{eqnarray*}

(7)

where σ_i_ is the wavelength dependent cross-section of species. σ_i_ is found from the sum of the molecular absorption cross-section, σ_a_(λ), and the Rayleigh scattering cross-section (Liou 2002),

\begin{eqnarray*} \sigma _{\rm r}(\lambda) = C_{\rm r} \frac{8\pi ^3(m_{\rm r}(\lambda)^2-1)^2}{3\lambda ^4N_{\rm s}}, \end{eqnarray*}

(8)

where _m_r is the real part of the refractive index of the molecule and _C_r is a corrective factor due to molecular anisotropy. _N_s is the number density at STP (1 atm, 273.15 K). For most species, _C_r is approximately unity, apart from CO2, for which _C_r is approximately 1.14 (Sneep & Ubachs 2005). The refractive index can be approximated as that of the dominant species in the atmosphere, which for the atmospheres modelled in this work is H2.

The host star of HD 209458b is type G0V, like the sun, and therefore the spectrum used was the unattenuated solar flux at the top of the Earth’s atmosphere, weighted by the relative flux received due to differing orbital radius between Earth and the sun (d⊕) and the planet and its host star (_d_r), (d⊕/_d_r)2. For HD 209458b, _d_r is approximately 0.047 au.

HD 189733b orbits a star of type K1.5V. The closest analogue that could be found for this was the spectrum of epsilon Eridani, a K2V star, for wavelengths above 115 nm (Youngblood et al. 2016), and the solar flux below 115 nm. Since HD 189733b orbits at a distance of 0.031 au, the spectrum is corrected to account for this in the same way as described above for HD 209458b.

Both these spectra are shown in Fig. 2.

Figure 2.

The unattenuated UV fluxes applied to the top of the atmospheres of HD 189733b and HD 209458b. For HD 189733b, the solar flux is used up to 115 nm and epsilon Eridani beyond that. For HD 209458b, the solar flux is used. In both cases, the flux has been corrected to account for planet–star separation.

2.2 Chemical network

Levi is able to use an external chemical network to compute the reactions in the atmosphere, and the rate that at which they occur. In this work, the chemical network being used is a subset of the Stand2018 Atmospheric Chemical Network first developed in Rimmer & Helling (2016). The version being used is a H/C/N/O/He network, containing approximately 150 different chemical species and 1000 forward reactions. To improve the generality of the code, it is possible to selectively disable certain reaction types, allowing the code to be applied to a much wider range of situations using the same network. The reactions that the model uses can be split into three general types: bimolecular, termolecular, and photochemical.

2.2.1 Bimolecular reactions

The general form of a bimolecular reaction is,

\begin{eqnarray*} A + B \rightarrow C + D, \end{eqnarray*}

(9)

for which the reaction rate is calculated using the generalized Arrhenius equation,

\begin{eqnarray*} k = \alpha \left(\frac{T}{300}\right)^{\beta }\exp {\left(-\frac{\gamma }{T}\right)}, \end{eqnarray*}

(10)

where k (m3 s−1) is the rate coefficient, T (K) is the temperature, and α (m3 s−1), β, and γ (K) are specific constants for each reaction, determined by fitting the Arrhenius equation to experimentally found or theoretically calculated values of k. Reactions of this form, and their rate constants, are tabulated in a supplementary file with the label A or RA, and are numbered 194–973 and 1040–1044.

2.2.2 Termolecular reactions

Termolecular reactions typically occur as either decomposition reactions,

\begin{eqnarray*} A + M \rightarrow C + D + M, \end{eqnarray*}

(11)

or combination reactions,

\begin{eqnarray*} A + B + M \rightarrow C + M, \end{eqnarray*}

(12)

where the M represents any third body that is involved in the reaction but is chemically unaffected as a result, often acting as either an energy source or sink for decomposition or combination reactions, respectively. To calculate the rate coefficients for these reactions, the Lindemann form is used (Lindemann et al. 1922); this form separates the rate constants into high and low pressure limits,

\begin{eqnarray*} k_0 &=& \alpha _0 \left(\frac{T}{300}\right)^{\beta _0}\exp {\left(-\frac{\gamma _0}{T}\right)}, \nonumber \\ k_{\infty } &=& \alpha _{\infty } \left(\frac{T}{300}\right)^{\beta _{\infty }}\exp {\left(-\frac{\gamma _{\infty }}{T}\right)}, \end{eqnarray*}

(13)

with _k_0 (⁠|$\mathrm{m^{6}\, s}^{-1}$| for combination or m3 s−1 for decomposition) as the low-pressure rate constant and _k_∞ (m3 s−1 for combination or s−1 for decomposition) as the high-pressure rate constant. These can be combined to form an effective two-body reaction rate,

\begin{eqnarray*} k = \frac{k_0 [M]}{1+\frac{k_0 [M]}{k_{\infty }}} , \end{eqnarray*}

(14)

where [_M_] (m−3) is the number density of the third species, where the third species can be any molecule in the atmosphere. Reactions of this form, and their rate constants, are tabulated in a supplementary file with the label B, and are numbered 1–193.

2.2.3 Reverse reactions

It is necessary to consider that the reactions described in the previous sections will also occur in the opposite direction, i.e. the previous products will react together to form the reactants. The rate at which this happens can be calculated either with the Arrhenius equation, with the constants being determined experimentally as before, or by thermodynamically reversing the reaction (Burcat & Ruscic 2005). Due to the lack of experimental data for these reverse reactions, especially at the temperatures being considered, we are using thermodynamic reversal to determine the reaction rate.

The equation used to calculate the reverse reaction rate is

\begin{eqnarray*} \frac{k_f}{k_r} = K_c \left(\frac{k_b T}{P_0}\right)^{-\Delta \nu }, \end{eqnarray*}

(15)

where _k_f and _k_r are the rate of the forward and reverse reactions, _k_b is the Boltzmann constant, Δν is the number of products minus the number of reactants, _P_0 is the standard-state pressure (which is 1 bar, under the ideal gas assumption). _K_c is an equilibrium constant determined by

\begin{eqnarray*} K_{\rm c} = \exp {\left(-\frac{\Delta G^o}{\mathcal {R}T}\right)}, \end{eqnarray*}

(16)

where |$\mathcal {R}$| is the gas constant and Δ_G_ o is the standard state gibbs free energy change on reaction, which is expressed as enthalpy and entropy changes on reaction:

\begin{eqnarray*} \Delta G^o = \Delta H^o - T \Delta S^o, \end{eqnarray*}

(17)

where Δ_H_ o and Δ_S_ o are the standard state change in enthalpy and entropy, respectively, upon reaction. The enthalpy and entropy of each species can be calculated by a series in terms of the NASA polynomials,

\begin{eqnarray*} \frac{H^o}{\mathcal {R}T} &=& a_1 + \frac{a_2}{2} T + \frac{a_3}{3}T^2 + \frac{a_4}{4}T^3 + \frac{a_5}{5}T^4 + \frac{a_6}{T}, \nonumber \\ \frac{S^o}{\mathcal {R}} &=& a_1 ln T + a_2 T + \frac{a_3}{2}T^2 + \frac{a_4}{3}T^3 + \frac{a_5}{4}T^4 + a_7. \end{eqnarray*}

(18)

The values for these polynomials were drawn from Burcat2 whenever possible, and when not possible, used a modified Benson’s method, described in the appendix of Rimmer & Helling (2016).

2.2.4 Photochemistry

The rate of photodissociation for photochemically active molecules is also calculated. The molecular absorption cross-sections, σ_a_ (m2), for these molecules were drawn primarily from PhIDRates3(Huebner et al. 1979, 1992, 2015), apart from C4H2 and C4H4, which were not in this data base. These two molecules used the cross-sections from the MPI-Mainz-UV-VIS Spectral Atlas of Gaseous Molecules4(Keller-Rudek et al. 2013). These cross-sections are often temperature dependant, but there is a significant lack of data at high temperatures, and so often the cross-sections we use are measured at room temperature. There has been some work done in identifying the UV cross-section for CO2 at high temperatures (Venot et al. 2013, 2018), in which several orders of magnitude difference were found in the absorption cross-section, resulting in significant differences to the computed chemical composition. This emphasizes the need for more high temperature cross-sections to accurately compute the photochemistry occurring in hot Jupiter atmospheres. The cross-sections were divided into 3000 bins, each 0.1 nm wide, to cover the approximately 300 nm range over which most cross-sections are defined. The photodissociation rate is then,

\begin{eqnarray*} k = t_f \int \sigma _a(\lambda) F(\lambda , z) {\rm d}\lambda \end{eqnarray*}

(19)

with F(λ, z) as the Actinic Flux as described in Section 2.1.3, and tf a dimensionless factor accounting for the fraction of time the atmosphere being modelled spends receiving stellar irradiation. This is 0 or 1 for a tidally locked planet’s night and day sides, respectively, and on average 1/2 for a planet with a day-night cycle. Photodissociation reactions are not reversed, since, while it is possible for reactions of this form to occur, they occur too slowly to have any impact on the atmosphere. Photodissociation reactions and their rate constants are tabulated in a supplementary file with the label V, numbers 974–1039.

2.3 Numerical method

Here, we describe how the governing equations presented in the previous section are discretized into a form that can be solved numerically.

2.3.1 Discretization

Equation (1) can be rewritten in a discrete form that allows it to be solved numerically,

\begin{eqnarray*} \frac{\Delta n_{i,j}}{\Delta t} = \mathcal {P}_{i,j} - \mathcal {L}_{i,j} - \frac{ \Phi _{i,j+1/2} - \Phi _{i,j-1/2}}{\Delta z}, \end{eqnarray*}

(20)

where, as previously, the subscript i refers to the _i_th species. The subscript j represents the j_th layer in the atmosphere, and takes values from j = 1,...,Nz, with Nz being the total number of atmosphere layers. The diffusion terms have subscripts j + 1/2 and j − 1/2 and thus are defined to be at the upper and lower boundary of each layer. Here, it has been assumed that only adjacent layers will contribute to the incoming fluxes, as a first-order treatment of diffusion is being used (Hu et al. 2012). The layer thickness is given by Δ_z. Equation (2) can be similarly discretized to

\begin{eqnarray*} \Phi _{i,j\pm 1/2} &=& \mp (K_{zz,j\pm 1/2}+D_{i,j\pm 1/2})\ n_{t,j\pm 1/2}\ \frac{ X_{i,j\pm 1} -X_{i,j}}{\Delta z}\nonumber \\ && + D_{i,j\pm 1/2}\ n_{i,j\pm 1/2}\left[\frac{(m_{t,j\pm 1/2} - m_i)g}{k_b T_{j\pm 1/2}} \mp \frac{\alpha _{T,i}}{T_{j\pm 1/2}}\frac{T_{j\pm 1} - T_j}{\Delta z}\right],\nonumber \\ \end{eqnarray*}

(21)

in which any Y j ± 1/2 has been approximated to (Y j ± 1 + Yj)/2. mt, j is the mean mass of atmosphere layer j, and mi is the mass of species i, g is the gravitational acceleration, and _k_b is the Boltzmann constant. Now, it is possible to combine equations (20) and (21) to form a set of ordinary differential equations with time as the independent variable:

\begin{eqnarray*} \frac{\Delta n_{i,j}}{\Delta t} &=& \mathcal {P}_{i,j} - \mathcal {L}_{i,j} + \frac{1}{\Delta z_{j-1/2}}\left[\left(k^-_{i,j+1/2}\frac{n_{t,j+1/2}}{n_{t,j+1}}\right)\ n_{i,j+1} \right.\nonumber \\ &&-\,\left(k^+_{i,j+1/2} \frac{n_{t,j+1/2}}{n_{t,j}} + k^-_{i,j-1/2}\frac{n_{t,j-1/2}}{n_{t,j}}\right)\ n_{i,j} \nonumber\\ &&\left.\,+ \left(k^+_{i,j-1/2}\frac{n_{t,j-1/2}}{n_{t,j-1}}\right)\ n_{i,j-1}\right], \end{eqnarray*}

(22)

where

\begin{eqnarray*} k^\pm _{j+1/2} &=& \frac{K_{zz,j+1/2} + D_{i,j+1/2}}{\Delta z_{j}} \nonumber \\ &\pm&\frac{D_{i,j+1/2}}{2\Delta z_j}\left[\frac{(m_{t,j}-m_i)g\ \Delta z_j}{k_b T_{j+1/2}} - \frac{\alpha _{T,i}}{T_{j+1/2}}(T_{j+1}-T_j)\right] \end{eqnarray*}

(23)

We define the atmosphere layers to be spaced equally in log pressure, and thus to find the corresponding altitude hydrostatic equilibrium is assumed,

\begin{eqnarray*} P(z) = P_0\ {\rm e}^{\frac{-m(z)gz}{k_{\rm b} T(z)}}, \end{eqnarray*}

(24)

with m as the mean mass of the molecules in an atmosphere layer, found by |$\sum \limits _i m_i n_i$|⁠, where mi is the mass of a molecule of species i. g is the surface gravity, defined for gas giants as the gravity at 1 bar, which is set to 10 m s−2 for an isothermal profile, 21.4 m s−2 for HD 189733b and 9.36 m s−2 for HD 209458b in accordance with Moses et al. (2011). The size of the atmospheric layers are recalculated at every time-step to ensure the atmosphere stays in hydrostatic equilibrium, since the mass of each layer varies as the system evolves.

2.3.2 Solver

Levi uses a semi-implicit second-order Rosenbrock solver to find a steady-state solution to equation (22), (Rosenbrock 1963). Verwer et al. (1999) show that the second-order Rosenbrock solver is stable over large range of step-sizes. This property of the solver is ideal for chemical kinetics, where step-sizes can range over many orders of magnitude due to the possibility of chemical reactions that have sub-microsecond time-scales and diffusion which can occur over time-scales of hundreds of years. The general differential equation is

\begin{eqnarray*} \frac{{\rm d} n}{{\rm d} t} = f(n), \end{eqnarray*}

(25)

where n is the dependent variable, t is the independent variable and f(n) is some general function of n. The Rosenbrock solver provides a solution to equation (25) in the form,

\begin{eqnarray*} n_{k+1} &=& n_k + 1.5 \Delta t g_1 + 0.5 \Delta t g_2, \nonumber \\ ({\bf I} - \gamma \Delta t {\bf J})g_1 &=& f(n_k), \nonumber \\ ({\bf I} - \gamma \Delta t {\bf J})g_2 &=& f(n_k + \Delta t g_1) - 2 g_1, \end{eqnarray*}

(26)

where Δ_t_ is the size of the time-step, γ is a scaler constant, suggested by Verwer et al. (1999) to be |$\gamma =1+1/\sqrt {2}$| for stability, the subscript k denotes the _k_th step in the solver, I is the identity matrix, and

\begin{eqnarray*} {\bf J} \equiv \frac{\partial f}{\partial n}|_{n_k} \end{eqnarray*}

(27)

is the Jacobian matrix evaluated at nk.

To solve equation (26), it is necessary to invert the Jacobian matrix, which is one of the most expensive parts of the calculations, and so we focused on making this as efficient as possible. The Jacobian takes the form of a tridiagonal block matrix of total size Ni Nz, the product of the total number of species, and the total number of atmosphere layers. The construction of this matrix is presented in equation (28). The diagonal blocks, B, of the Jacobian are made up of square matrices of size Ni that describe the interactions between species in the same atmosphere layer. The off-diagonal blocks, A and C, also matrices of size Ni, represent the effects of adjacent atmosphere layers via diffusion. The subscripts for the A, B, and C matrices show which atmosphere layer they represent.

\begin{eqnarray*} J = {\begin{bmatrix}B_{1} & C_{1} & 0 & 0&\dots &0&0&0& 0 \\ A_{2} & B_{2} & C_{2} & 0 & \dots &0&0&0& 0 \\ 0 & A_{3} & B_{3} & C_{3} & \dots &0&0&0& 0 \\ 0 & 0 & A_{4} & B_{4} & \dots &0&0&0& 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots &\vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \dots & B_{N_z -3} & C_{N_z - 3} & 0 & 0 \\ 0 & 0 & 0 & 0 & \dots & A_{N_z -2} & B_{N_z - 2} & C_{N_z -2} & 0 \\ 0 & 0 & 0 & 0 & \dots & 0 & A_{N_z -1} & B_{N_z - 1} & C_{N_z -1} \\ 0 & 0 & 0 & 0 & \dots & 0 & 0 & A_{N_z} & B_{N_z} \\ \end{bmatrix}} \end{eqnarray*}

(28)

By comparison of equations 5, 22, 25, and 27, it can be seen that

\begin{eqnarray*} A &=& \frac{k^-_{i,j+1/2}\ n_{t,j+1/2}}{dz_{j-1/2}\ n_{t,j+1}}, \nonumber \\ C &=& \frac{k^+_{i,j-1/2}\ n_{t,j-1/2}}{dz_{j-1/2}\ n_{t,j-1}}, \nonumber \\ B &=& -\left(\frac{k^+_{i,j+1/2}\ n_{t,j+1/2}}{dz_{j-1/2}\ n_{t,j}}+\frac{k^-_{i,j-1/2}\ n_{t,j-1/2}}{dz_{j-1/2}\ n_{t,j}}\right) \nonumber \\ && +\, \frac{\partial }{\partial n_{j}} (\sum (k_2\ n_{i^{\prime \prime }}\ n_{i^{\prime \prime \prime }} - k_1\ n_i\ n_{i^{\prime }})). \end{eqnarray*}

(29)

As a result, both A and C are purely diagonal matrices, while B is a full matrix with each position describing how one species reacts with another in the same layer.

To invert this Jacobian, we use the method of inverting a tridiagonal matrix described in Hubeny & Mihalas (2014), with a generalization to block tridiagonal by treating each matrix element in the inversion method as a full matrix. As a result, only the three Ni size matrices are ever kept in the memory, which is much less expensive than trying to invert the complete Ni Nz matrix.

2.3.3 Convergence

There are several conditions that must be satisfied before the run is considered to be complete and have reached steady state. These are that both the relative change, |$\Delta n = \frac{|n_{i,j,k+1} - n_{i,j,k}|}{n_{i,j,k}}$|⁠, and the rate of relative change, |$\frac{\Delta n}{\Delta t}$|⁠, of mixing ratios in the atmosphere are sufficiently small that running the code further would produce no significant change.

To ensure that each step is correctly approaching steady state, several conditions are checked after every time-step. If the conditions are passed then the solver continues, otherwise the latest solution is discarded and the step-size is reduced before re-running the solver. The conditions are that every species has a positive number density and that the truncation error, ϵ, is sufficiently small.

The truncation error is found by calculating the difference between the second- and first-order solutions, |$n^1_{k+1} = n_k + \Delta t g_1$|⁠, such that,

\begin{eqnarray*} \epsilon = |n_{k+1} - n^1_{k+1}|. \end{eqnarray*}

(30)

Since the time-scales being modelled can vary over many orders of magnitude, having a time-step that is self-adjusting is vital to allowing the code to have a reasonable runtime. The truncation error can thus be used as a cheap local error estimation to control the step size,

\begin{eqnarray*} \Delta t_{k+1} \sim \Delta t_{k} (1/\epsilon)^{0.5}. \end{eqnarray*}

(31)

2.4 Initial and boundary conditions

2.4.1 Initial conditions

In this work, we focus on hydrogen-rich gas giant atmospheres. Our baseline model uses solar elemental abundances from Asplund et al. (2009). In Section 4, we also explore the effects of variation in the C/O and N/O ratios. In Asplund et al. (2009), the solar abundances are given as fractions of the number of hydrogen molecules, while for Levi the elemental abundances have been transformed to the fraction of the total number of molecules (i.e. the mixing ratio), enabling easier description of non-hydrogen dominated atmospheres. The ratio of elemental atoms to the total number of molecules therefore starts as _X_H = 1.707, _X_He = 0.145, _X_C = 4.594 × 10−4, _X_N = 1.154 × 10−4, and _X_O = 8.359 × 10−4. When calculating solar abundances, all H is assumed to be in the form of H2, thus leading to a _X_H that is greater than one. This gives a C/O ratio of 0.55 and a N/O ratio of 0.138.

When choosing the initial conditions, conservation of mass is ensured by summing over the product of the species present (ni) and the number of atoms of x contained by ni (Ai, x),

\begin{eqnarray*} \sum \limits _{i} A_{i,x} n_i = X_x . \end{eqnarray*}

(32)

The initial chemistry of the atmosphere is chosen as a set of species that contain all the elements in the atmosphere, then by using equation (32) for each element to solve for the species’ initial abundance and setting the initial abundance of all other species in the network to zero. To simplify solving equation (32), we picked the initial species set to be small molecules, so that dividing the elements between them would be straightforward.

With this approach, there are still a range of viable starting conditions for the atmospheric mixing ratio and it is important to demonstrate that model convergence is not affected by this decision. A comparison of the steady-state solution emerging from a choice of three initial atmospheric species sets is shown in the first plot of Fig. 3, with the initial species sets being (1) H2, CO, OH, and NH3; (2) H2, CH4, H2O, and N2; and (3) atomic H, C, O, and N. Having the initial species set being elements in their atomic form provides a starting condition of extreme disequilibrium, compared to starting with molecules that are expected to dominate the steady-state atmosphere, as in the other two starting conditions. There is no noticeable difference between the three sets of mixing ratios, confirming that the steady-state solution is independent of the choice of atmospheric species set to initialize the calculation with. The time taken to reach this steady state, however, is significantly different between the three sets, with set (1) taking four and three times as long as sets (2) and (3), respectively, in this scenario. Figs 3(b)–(d) show the accuracy to which each of the three different initial conditions are the same. In all cases, there was less than a 10 ppm difference between any two starting sets of molecules, thus suggesting that the solutions these initial sets produce are functionally identical.

Figure 3.

An investigation of the effect of initial conditions on atmospheric chemistry. The sets of initial species are H2, CO, OH, and NH3 (solid lines, labelled 1; H2, CH4, H2O, and N2 (dashed lines, labelled 2; and elements in their atomic form (dotted lines, labelled 3). In this example, a 1500K Isotherm was used, with no diffusion or photochemistry. Plot (a) shows all three starting conditions, while (b)–(d) show the percentage difference between two of the starting options.

2.4.2 Boundary conditions

It is necessary to specify boundary conditions at the minimum and maximum pressures chosen for the simulation. There are a variety of choices available, depending both on the type of planet and the range of pressures being explored.

For terrestrial exoplanets, the lower boundary conditions usually model any surface–atmosphere interaction, so typical boundary conditions can include permanently assigning a species an abundance at the surface to model a large reservoir, or a flux to describe surface emission and deposition. For the upper atmosphere, a flux due to atmospheric escape is common (Hu et al. 2012).

Gas-giant exoplanets have no solid surface at the lower boundary, thus the boundary conditions are either chemical equilibrium or zero flux. The upper boundary condition could be atmospheric escape or, if the upper boundary is deep enough in the atmosphere, zero flux. In general, atmospheric escape does not significantly affect the mass of the atmosphere (Murray-Clay, Chiang & Murray 2009; Linsky et al. 2010). For a boundary in chemical equilibrium, the assumption is that layers deeper in the atmosphere have a chemical time-scale much shorter than the dynamical time-scale, so the boundary would be in chemical equilibrium. A zero flux, or closed boundary, assumes that the deeper layers are chemically identical to the bottom layer in the calculation, such that there will be no chemical consequence to flow of molecules across the layer boundaries.

The effect of setting different lower boundary conditions are explored in the left plot in Fig. 4, which shows the difference between equilibrium (dashed lines) and zero flux (solid lines) boundary conditions. It can be seen that there is a slight difference in the mixing ratios as a result of these boundary condition changes, especially to N2. To better quantify this change, the right plot in Fig. 4 shows the percentage difference between zero flux and equilibrium boundary conditions. Throughout most of the atmosphere, the difference in abundance is less than one per cent; however, closer to the lower boundary it can increase up to 10 per cent. Thus, it is clear that the boundary conditions can have a much more significant impact upon atmospheric abundances than the initial conditions, although even a 10 per cent change is negligible when considering the precision with which the abundance of these molecules can be detected in exoplanet atmospheres. Since it is possible for diffusion to drive even the deepest layers slightly away from equilibrium, zero flux boundary conditions were chosen for this work to prevent a disparity between the boundary layer and the layer above. As a result, zero flux conditions were also chosen for the upper boundary to ensure conservation of mass in the atmosphere.

Figure 4.

The difference in steady-state atmospheric chemistry based on the lower boundary condition. In this example, a 1500 K isotherm was used, with 1 × 106 m s−2_Kzz_ and no photochemistry. Plot (a) shows a comparison of the atmospheric chemistry between the two boundary conditions, while (b) shows the percentage difference in the abundances of a zero-flux boundary condition compared to an equilibrium boundary condition.

3 VALIDATION AND TESTING

In this section, we apply our model to the atmospheres of several hot Jupiters. We first compare our models of HD 209458b to that of previous codes (Moses et al. 2011; Venot et al. 2012) for validation. We then investigate the effects of disequilibrium chemistry on the atmospheric chemistry of hot Jupiters and discuss why these effects occur. We also show the fastest reaction pathway the model takes for several important net reactions in the atmosphere. These pathways are found by use of a code that traces the reactions that lead from the products to the reactants, and identifying the route that has the fastest overall reaction rate. When these pathways are shown, the bracketed numbers adjacent to the reaction correspond to its position in the chemical network, as described in Section 2.2.

3.1 Validation

We chose HD 209458b as one of the best candidates for benchmarking Levi due to both the well characterized spectral observations of its atmosphere and the existence of several previous chemical models that have predicted atmospheric compositions for this planet (Moses et al. 2011; Venot et al. 2012; Agúndez et al. 2014; Rimmer & Helling 2016; Tsai et al. 2017). The results produced by this code are compared to the disequilibrium models of HD 209458b of both Moses et al. (2011) and Venot et al. (2012) – see Fig. 5. Like Levi, both other models do not self-consistently calculate the P–T profile, but rather use a fixed input. The species chosen for the comparison are H2, H, and H2O (Fig. 5a); CO, CH4, and CO2 (Fig. 5b); and N2, NH3, and HCN (Fig. 5c). These species are expected to be some of the most abundant (e.g. Swain et al. 2009; MacDonald & Madhusudhan 2017b), and several have strong spectral signatures.

Figure 5.

A comparison of models of the atmosphere of HD 209458b for a selection of major molecules in the atmosphere. The solid line is the model described in this paper, the dashed line is the model created by Moses et al. (2011), and the dotted line is the code described in Venot et al. (2012).

Looking at the comparison plots in Fig. 5, it can be seen that for most species there is no significant difference between Levi and Venot et al. (2012), or between Levi and Moses et al. (2011). In the deep atmosphere there are some slight differences in the abundances of many molecules. These differences are not significant, and are likely a consequence of slightly differing P–T profiles, thermodynamic coefficients, and choices of initial elemental ratios. For example, Moses et al. (2011) and Venot et al. (2012) model a lower-oxygen atmosphere, because they have considered 20 per cent of the oxygen abundance to be sequestered in silicates.

One of the largest differences between the models is for NH3, where there is an half an order-of-magnitude lower quenching abundance of NH3 for Levi compared to Moses et al. (2011), while Levi predicts nearly a hundred times more NH3 at 10−6 bar. Compared to Venot et al. (2012), Levi predicts a nearly identical quenching abundance of NH3, and almost 10 000 times more NH3 at 10−6 bar. These differences are due to a high degree of uncertainty in the rate of the pathway that converts NH3 to N2:

\begin{eqnarray*} \rm {H_2 + M} \,&\rightarrow& \, {\rm H + H + M} \ \rm (R5), \nonumber \\ {\rm 2(H + NH_3} \,&\rightarrow& \, \rm NH_2 + H_2) \ (R241), \nonumber \\ {\rm 2(H + NH_2} \,&\rightarrow& \, \rm HN + H_2) \ (R220), \nonumber \\ \end{eqnarray*}

\begin{eqnarray*} {\rm 2HN} \,&\rightarrow& \, {\rm N_2 + 2H \ (R383)}, \nonumber\\ {\rm {Net}{:}}\ {\rm 2NH_3} \,&\rightarrow& \, \rm N_2 + 3H_2. \end{eqnarray*}

(33)

The observed difference in the quenching location of CH4 arises for similar reasons. In the network used by Levi, there are several pathways that convert CH4 into CO, one of which has a highly uncertain reaction rate (Rimmer & Helling 2016), and can explain the nearly 1 and 0.5 dex higher quenching abundance of CH4 seen in Moses et al. (2011) and Venot et al. (2012), respectively, than in Levi. The large differences in abundance seen in HCN between the three models throughout the atmosphere are due to high degrees of uncertainty in the rate constants of many of the reactions that include HCN.

To conclude, Levi can reproduce the composition of HD 209458b predicted by several other chemical models at an acceptable degree of accuracy. While there are still some significant differences in abundances of several highly abundant species in exoplanet atmospheres, these are the result of uncertainties in rate constants and absorption cross-sections, which differ between the models. Only new experiments and calculations of the rate constants and absorption cross-sections will help resolve these differences.

3.2 Equilibrium versus diffusion versus photochemistry

In this section, the effects of the three main chemical processes that occur in atmospheres are compared and contrasted. Thermochemical equilibrium, diffusion in the form of molecular, thermal and eddy-diffusion, and photochemistry can all combine to play an important role in determining the distribution of observable species in exoplanet atmospheres, and thus should be considered when investigating the chemistry of these exoplanets. Chemical equilibrium is expected to dominate in the lower atmosphere, up to a few bars, due to the high temperatures and pressure. Diffusion effects then take over for some species between approximately 1–10−4 bars, and photochemistry dominates above this. At higher altitudes, in the thermosphere, it is sufficiently hot that most species thermally dissociate down to their base elements. We apply our code to the hot Jupiters HD 209458b and HD 189733b and discuss how and why the equilibrium and disequilibrium models diverge. We run three model scenarios, with the results presented in Fig. 6: no-diffusion, no-photochemistry (solid lines); with diffusion, but no-photochemistry (dashed lines); with diffusion and photochemistry (dotted lines). Results similar to these have also been produced by previous 1D chemical codes (e.g. Zahnle et al. 2009; Moses et al. 2011; Hu et al. 2012; Venot et al. 2012; Rimmer & Helling 2016; Tsai et al. 2017).

Figure 6.

Abundance profiles for several species on HD 290458b (left) and HD 189733b (right) for purely equilibrium chemistry (solid lines), diffusion (dashed line), and both diffusion and photochemistry (dotted line). The P–T and the Kzz profiles for these models are shown in Figs 1(a) and (b).

3.2.1 Equilibrium

In the deepest regions of the atmospheres of the planets being investigated, the high temperature and pressure may produce an atmosphere in chemical equilibrium. The relative abundances of species at chemical equilibrium depend on the P–T profile. For the investigated planets, at solar composition, we find that CH4 and H2O dominate over CO, and NH3 approaches the abundance of N2, while higher in the atmosphere this trend is reversed (Heng & Tsai 2016). When solving for equilibrium analytically, the same trend is found. In the observable region of the atmosphere, approximately 1–10−4 bar for emission spectrography, CO is the dominant carbon species, with most of the remaining oxygen ending up in H2O, and N2 is the dominant nitrogen species.

On HD 209458b, the thermosphere begins around 10−7 bar, causing most species to dissociate down to their base atomic forms. The thermosphere of HD 189733b does not begin until a higher altitude, so no thermal dissociation can be seen for this planet in Fig. 6. Notable for comparison later on are the abundances of HCN and NH3 in the observable region of the atmosphere; in this equilibrium scenario, they are very minor constituents, far below detectable limits (HCN or NH3 abundances at least 1 per cent of the H2O abundance) for solar composition (MacDonald & Madhusudhan 2017b).

3.2.2 Diffusion

The inclusion of diffusion, the dashed lines in Fig. 6, leads to some significant deviations from chemical equilibrium. At very high pressures, the temperature is high enough that the chemical time-scale of the reactions is much shorter than the dynamical time-scale, and therefore the atmosphere at these levels is still in chemical equilibrium, as predicted. However, as the pressure and temperature decrease, and reaction rates drop due to the decreased rates of molecule-molecule interactions, the rate of some reactions drop sufficiently such that the time-scale of dynamic motion is now shorter than that of chemical interaction. This results in quenching at the pressure level where the two time-scales are equal, with the species abundance being effectively frozen in at this point and this chemistry being transported higher into the atmosphere by eddy diffusion. This process leads to the possibility of significantly different abundances of some species in the observable region of the atmosphere, compared to chemical equilibrium.

Fig. 6 shows a number of species being affected by quenching in observable regions of the atmosphere, including CH4, HCN, and NH3. NH3 is of particular note as quenching increases its abundance by many orders of magnitude, on both HD 209458b and HD 189733b, to potentially detectable levels (MacDonald & Madhusudhan 2017b).

On HD 209458b, the high temperatures lead to a shorter reaction time-scale compared to the mixing time-scale, implying that the condition for chemical equilibrium to be in effect is satisfied until much lower pressures in the atmosphere, compared to HD 189733b where quenching occurs deeper in the atmosphere. HD 189733b has many species at their quenched abundance throughout the upper atmosphere. This is due to the planet’s upper atmosphere not having a thermal inversion below the thermosphere, so the chemical time-scale monotonically increases with decreasing pressure. On the other hand, HD 209458b does have a thermal inversion, so the chemical time-scale can decrease with decreasing pressure, and so the upper atmosphere may begin to return to chemical equilibrium.

3.2.3 Photochemistry

The dotted lines in Fig. 6 show how the inclusion of photochemistry affects the atmosphere. The photolysis of CO, H2O, N2, NH3, and CH4 drive most of the photochemical reactions that occur in the upper atmosphere.

The destruction of H2O by a UV photon into OH and H, as deep as 0.1 bar, sets up the reaction that turns H2 into H:

\begin{eqnarray*} {\rm H_2O + }\ h\nu \,&\rightarrow&\, \rm OH + H \ (R992), \nonumber \\ {\rm H_2 + OH} \,&\rightarrow& \,\rm H_2O + H \ (R221),\nonumber\\ {Net}{:}\ {\rm H_2 + }\ h\nu \,&\rightarrow&\, \rm 2H. \end{eqnarray*}

(34)

This reaction provides a large pool of highly reactive hydrogen radicals that can diffuse to both higher and lower altitudes to cause a chain of further reactions, drastically altering the composition of the upper atmosphere. On HD 209458b the rate of this reaction is low compared to the rate of dissociation of molecules due to the high temperatures in the upper atmosphere, and so little difference can be seen between the case with and without photochemistry. As evidenced in Figs 6(a) and (b), H2O is not permanently depleted from the atmosphere until much lower pressures due to it being replenished as fast as it is destroyed.

The photolysis of CO can produce C and O radicals that can also contribute to the destruction of H2. Like H2O, CO is also being replenished quickly, in this case by a reaction between CH4 and H2O:

\begin{eqnarray*} {\rm 2(H_2 + M} \,&\rightarrow& \, \rm 2H +M) \ (R5), \nonumber \\ {\rm H_2O + H} \,&\rightarrow& \, \rm OH + H_2 \ (R221), \nonumber \\ {\rm CH_4 + H} \,&\rightarrow& \, \rm H_2 + CH_3 \ (R251), \nonumber \\ {\rm CH_3 + OH} \,&\rightarrow& \, \rm CH_2O + H_2 \ (R479), \nonumber \\ {\rm H + CH_2O} \,&\rightarrow&\ \, \rm H_2 + CHO \ (R236), \nonumber \\ {\rm H + CHO} \,&\rightarrow& \, \rm CO + H_2 \ (R216), \nonumber\\ {Net}{:}\ {\rm CH_4 + H_2O} \,&\rightarrow& \, \rm CO + 3H_2. \end{eqnarray*}

(35)

However, as Fig. 6(d) shows, on HD 189733b the formation of CO is not efficient enough to prevent loss of CO above 10−5 bar and so some of the carbon can end up in other species, e.g. CO2 and HCN. On HD 209458b, the temperature and thus the reaction rates are high enough to keep replenishing CO until pressures as low as 10−7 bar.

The C radicals from the photolysis of CO often end up producing CH4:

\begin{eqnarray*} {\rm CO} + \ h\nu \,&\rightarrow& \, \rm C + O \ (R979), \nonumber\\ {\rm C + H_2} \,&\rightarrow& \, \rm CH + H \ (R195),\nonumber\\ {\rm CH + H_2} \,&\rightarrow& \, \rm CH_2 + H \ (R212),\nonumber \\ {\rm CH_2 + H_2} \,&\rightarrow& \, \rm CH_3 + H \ (R235),\nonumber \\ {\rm CH_3 + H_2} \,&\rightarrow& \, \rm CH_4 + H \ (R251),\nonumber \\ {Net}{:}\ {\rm CO + 4H_2} \,&\rightarrow&\, \rm CH_4 + O + 4H, \end{eqnarray*}

(36)

which leads to the increase of methane seen on HD 209458b. CH4 is, in turn, susceptible to destruction by H radicals, thus leading to an overall decrease in methane if sufficient H is produced by the photodissociation of H2O, as discussed earlier, and as can be seen at 10−5 bar on HD 189733b.

CO2 is not greatly affected by photochemistry, since its production is based upon fast reactions between H2O and CO, both of which are largely constant until high in the atmosphere. In the upper atmosphere of HD 189733b, some of the CO is transformed into CO2.

C2H2 is an important by-product of photochemistry. The carbon atoms released from dissociation of CH4 or CO are reprocessed to form C2H2. This is the result of a highly efficient pathway:

\begin{eqnarray*} {\rm CO} + \ h\nu \,&\rightarrow& \,\rm C +O \ (R979), \nonumber \\ {\rm C + H_2} \,&\rightarrow& \, \rm CH + H \ (R195), \nonumber \\ {\rm CH + H_2} \,&\rightarrow& \, \rm CH_2 + H \ (R212), \nonumber \\ {\rm CH_4} + \ h\nu \,&\rightarrow& \, \rm CH_3 +H \ (R1019), \nonumber \\ {\rm CH_3 + M} \,&\rightarrow& \, \rm CH_2 + H + M \ (R29), \nonumber \\ {\rm CH_2 + CH_2} \,&\rightarrow& \, \rm C_2H_2 + H_2 \ (R614), \\ {Net}{:}\ {\rm CO + CH_4 + 2H_2} \,&\rightarrow& \,\rm C_2H_2 + 4H, \end{eqnarray*}

(37)

which leads to many orders of magnitude more C2H2, on both HD 209458b and HD 189733b, than otherwise expected if there was no UV flux.

Photochemistry can also speed up the reaction converting NH3 to HCN,

\begin{eqnarray*} {\rm NH_3} + \ h\nu \,&\rightarrow& \, \rm NH_2 + H \ (R1013), \nonumber\\ {\rm NH_2} + \ h\nu \,&\rightarrow& \, \rm HN + H \ (R999), \nonumber\\ {\rm H + NH} \,&\rightarrow& \, \rm H + H + N \ (R204), \nonumber\\ {\rm CH_4 + M} \,&\rightarrow& \,\rm CH_3 + H + M \ (R34), \nonumber\\ {\rm CH_3 + N} \,&\rightarrow&\, \rm HCN + H_2 \ (R303), \nonumber\\ {\rm 2(H + H + M} \,&\rightarrow& \, \rm H_2 + M) \ (R5), \nonumber\\ {Net}{:}\ {\rm NH_3 + CH_4} \,&\rightarrow& \,\rm HCN + 3H_2, \end{eqnarray*}

(38)

causing a significant increase in the amount of HCN. On HD 209458b, this proceeds to the point where HCN is the second most abundant carbon and nitrogen bearing molecule in the atmosphere, and on HD 189733b, it is the most abundant carbon and nitrogen bearing molecule. This reaction also causes the large decrease of NH3 seen in the upper atmosphere of HD 189733b, compared to models without photochemistry. Although loss processes exist for HCN, the products are normally not very stable, and often form HCN upon their destruction. The result is that destruction reactions are inefficient in regulating HCN’s abundance, especially with other pathways forming it fast enough to replenish any loss.

4 EXPLORATION OF THE CHEMICAL PARAMETERS

In this section, we perform an initial exploration into the parameter space of the C/O ratio and N/O ratio, and consider how it may affect the detectability of certain species in exoplanet atmospheres. Apart from where stated otherwise, these models use the P–T profile and Kzz profile described in Section 2.1.2, the initial conditions labelled (a) and the atmospheric composition from Section 2.4.1 and the zero-flux boundary conditions from Section 2.4.2. Finally, we use a new P–T profile for HD 209458b that does not contain a thermal inversion to investigate the C/O and N/O ratios that would best fit with new evidence of multiple species on HD 209458b.

4.1 The C/O and N/O ratio

Throughout our previous models, we have used a solar composition for the atmospheres of the planets being investigated (Section 2.4.1). However, there is evidence that some hot Jupiters have C/O ratios substantially different from solar values (Madhusudhan, Burrows & Currie 2011, Stevenson et al. 2014). A planet’s composition depends on many factors, mainly arising from its formation and migration history. Thus, knowledge of the atmospheric composition of a planet can lead to insights into its past (Öberg et al. 2011; Madhusudhan, Amin & Kennedy 2014b; Mordasini et al. 2016). One way of finding an atmosphere’s composition is through consideration of the effects that a change in composition would have upon species in that atmosphere.

In this section, we independently vary the C/O and N/O ratios for a planet otherwise equivalent to HD 189733b. We investigate the effects that composition has upon the abundance of a number of species in the atmosphere, how it affects the possibility of detecting some of these species and thus how easily detectable these changes in composition are. We explore a parameter space in which the total amount of carbon and nitrogen can vary between 0.1 and 10 times solar. This corresponds to an atomic C/O ratio varying between 0.055 and 5.5, and an N/O ratio varying between 0.0138 and 1.38. The atomic fraction of oxygen is always kept constant, with the change in carbon or nitrogen abundance being accounted for by an equivalent but opposite change in the total amount of hydrogen. We pick 100 mbar as the pressure being investigated since both emission and transmission spectra are sensitive to the abundance of species at this pressure. These models use the P–T and Kzz profile of HD 189733b from Figs 1(a) and (b).

In Fig. 7 and 8 maps of the abundance of C and O species are shown as a function of atmospheric composition. None of the species show any strong dependence on the amount of N in the atmosphere. In hydrogen-dominated atmospheres, the abundances of H2O, CH4, and CO are primarily determined by the equation:

\begin{eqnarray*} {\rm CH_4 + H_2O} \rightleftharpoons {\rm CO + 3H_2}, \end{eqnarray*}

(39)

which at 100 mbar in the atmosphere of HD 189733b favours the formation of CO. At solar C, a C/O ratio of 0.55, this means that CO is the primary carbon carrier and H2O contains most of the remaining oxygen. Since there is nearly twice as much oxygen as carbon in the atmosphere, CO and H2O have very similar abundances. For sub-solar carbon, there is less CO, and thus less oxygen bound to CO, resulting in an increase of H2O. Once C/O > 1 (1.8× solar carbon), oxygen, not carbon, is now the limiting factor in producing CO, and so the abundance of water quickly drops as there is very little oxygen available to form it. Much of the excess carbon ends up in CH4, resulting in a sharp increase in its abundance. CO2 is also affected indirectly by equation (39), while there is excess oxygen not bound to CO, its abundance is directly related to the amount of carbon; however once C/O > 1, then the majority of the oxygen is bound to CO and the abundance of CO2 drops off rapidly. Larger hydrocarbons can also be very strongly affected by the variation in carbon, with the abundance of C2H2 changing by more than eight orders of magnitude for only two orders of magnitude change in carbon.

Figure 7.

Abundance variations of H2O (top), CH4 (middle), and CO (bottom) on HD 189733b at 0.1 bar for the parameter space in which the amount of nitrogen and carbon in the atmosphere varies between 0.1 and 10 times the solar amount. These models use the P–T and diffusion profile from Figs 1(a) and (b) for HD 189733b.

Figure 8.

Abundance variations of CO2 (top), C2H2 (bottom) on HD 189733b at 0.1 bar for the parameter space in which the amount of nitrogen and carbon in the atmosphere varies between 0.1 and 10 times the solar amount. These models use the P–T and diffusion profile from Figs 1(a) and (b) for HD 189733b.

In Fig. 9, the variation of the three major nitrogen species within the parameter space is shown. The abundance of N2 has no dependence on carbon, while both NH3 and HCN do. This is because of the fast reaction,

\begin{eqnarray*} {\rm H_2 + M} \,&\rightarrow&\, \rm H + H + M \ (R5), \nonumber \\ {\rm CH_4 + H} \,&\rightarrow& \, \rm CH_3 + H_2 \ (R251), \nonumber \\ {\rm H + NH_3} \,&\rightarrow& \, \rm NH_2 + H_2 \ (R241), \nonumber \\ {\rm H + NH_2} \,&\rightarrow& \, \rm HN + H_2 \ (R220), \nonumber \\ {\rm H + HN}\,&\rightarrow& \, \rm N + H + H \ (R204), \nonumber \\ {\rm N + CH_3} \,&\rightarrow& \, \rm HCN + H_2 \ (R302),\nonumber\\ {Net}{:}\ \rm {NH_3 + CH_4} \,&\rightarrow& \,\rm HCN + 3H_2, \end{eqnarray*}

(40)

thus causing the abundance of NH3 and HCN to have a dependence on the amount of carbon in the atmosphere. This dependence is particularly strong when C/O > 1, since there can now be an excess of CH4, without it being converted into CO. More CH4 therefore leads to more NH3 being converted into HCN. As expected, all the nitrogen species have a strong dependence on the amount of nitrogen in the atmosphere, though HCN has a weaker dependence than the others, due to carbon often being the limiting factor for its abundance. Many of these same trends for variation in the C/O ratio are seen in Madhusudhan (2012).

Figure 9.

Abundance variations of N2 (top), HCN (middle), and NH3 (bottom) on HD 189733b at 0.1 bar for the parameter space in which the amount of nitrogen and carbon in the atmosphere varies between 0.1 and 10 times the solar amount. These models use the P–T and diffusion profile from Figs 1(a) and (b) for HD 189733b.

In Fig. 10, the ratios of NH3 and HCN to H2O are shown, since this is what determines the observability of these species in exoplanet atmospheres. Since both NH3 and H2O decrease with increasing carbon, the overall increase in NH3/H2O is due to the more rapid decrease of H2O. Having HCN and NH3 abundances at least 1 per cent that of H2O is predicted to be needed to detect either of these molecules (MacDonald & Madhusudhan 2017b). Given this, Fig. 10 suggests that the amount of nitrogen is not very significant in the detection of these species, whereas C/O > 1 is almost always essential. However, exceptionally large amounts of nitrogen in the atmosphere can overcome this and allow detection of HCN and NH3 at slightly lower C/O ratios.

Figure 10.

Variations in the ratio of HCN and NH3 to H2O, top and bottom, respectively, on HD 189733b at 0.1 bar for the parameter space in which the amount of nitrogen and carbon in the atmosphere varies between 0.1 and 10 times the solar amount. These models use the P–T and diffusion profile from Figs 1(a) and (b) for HD 189733b.

4.2 Molecular detections on HD 209458b

The atmosphere of HD 209458b was originally believed to contain a thermal inversion (Knutson et al. 2008), and so far in this work we have modelled its atmosphere as such to enable validation of the code. More recent studies, however, have found that HD 209458b should not have a thermal inversion (Diamond-Lowe et al. 2014), and it seems valuable for us to examine how this changes the predicted abundance profiles for species in its atmosphere. In addition, recent work in the literature presents evidence of CO, H2O, and HCN in the atmosphere of HD 209458b (Hawker et al. 2018). Hawker et al. (2018) find the best-fitting mixing abundance of HCN to be 10−5, with a minimum mixing ratio of 10−6.5. They also find H2O and CO mixing ratios consistent with previous values of approximately 10−5 ± 0.5 and 10−4 ± 0.5, respectively (Madhusudhan et al. 2014a; MacDonald & Madhusudhan 2017a; Brogi & Line 2018).

In this section, we used a model P–T profile for the day-side of HD 209458b, without a thermal inversion (Gandhi & Madhusudhan 2017). For the deep adiabatic region of the atmosphere, above 100 bars, we use the same values as the earlier P–T profile for HD 209458b, translated by −100 K for a smooth connection (Fig. 11). The same Kzz profile and UV flux for HD 209458b from earlier are used, seen in Fig. 1(b) and Fig. 2, respectively.

Figure 11.

The P–T profile for the average day-side of HD 209458b, with no thermal inversion in the atmosphere.

In Fig. 12, the chemical abundance profiles for the new P–T profile of HD 209458b are compared with those from the old P–T profile. In the lower atmosphere, where the P–T profile is largely unchanged, there is very little difference between the old and new mixing ratios. In the upper atmosphere, the lower temperature of the new P–T profile, due to the lack of a thermal inversion, means that some species that would otherwise be depleted by the higher temperatures, such as HCN or CH4, are present in greater quantities. For HCN and CH4, this results in an increase in abundance of nearly two orders of magnitude at 10−3 bar.

Figure 12.

A comparison of the chemistry in the atmosphere of HD 209458b for two temperature profiles. The solid lines use a P–T profile with no thermal inversion, from Fig. 11, while the dashed lines use a P–T profile with a thermal inversion, from Fig. 1(a).

In Fig. 13, the abundance profiles for CO, HCN, and H2O are displayed for a range of C/O ratios. We can compare the mixing ratios of these species between 10−1 and 10−3 bar in our models to the estimated abundances from Hawker et al. (2018). The abundance of CO provides little information since across all C/O ratios being modelled it stays at approximately 10−3, consistent with the expected value. HCN provides more information, with a C/O ≳ 0.9 required to pass the minimum estimated mixing ratio of 10−6.5, and a C/O = 1.2 is required for HCN to match the best-fitting abundance of 10−5. However, H2O provides an upper limit of C/O ≈ 1; higher C/O than this quickly depletes H2O below 10−6. As such, it is clear that alterations to the C/O ratio by itself may not be sufficient to explain the observed abundances.

Figure 13.

The mixing ratios of CO, HCN, and H2O on HD 209458b, using the P–T profile from Fig. 11. Each of the molecules are displayed at different C/O and N/O ratios. The C/O ratio is 0.9 for the two N/O ratios shown.

In Fig. 13, we also included two models in which we increased the N/O ratio by two times (N/O = 0.28) and ten times (N/O=1.4) the solar value, to discover the effect this has on the abundance of these molecules. For these models, we kept C/O = 0.9. Increasing the amount of nitrogen in the atmosphere by two times solar increases the HCN abundance by approximately a third at 10−3 bar. Further increases to the amount of nitrogen in the atmosphere do little to the abundance of HCN since the amount of free carbon is the limiting factor beyond this. The N/O ratio has no significant effect on the abundance of either CO or H2O.

Through comparison of our model and the expected molecular abundances of HD 209458b, we can produce some initial constraints for the C/O and N/O ratio of this hot Jupiter: a 0.9 ≲ C/O ≳ 1 and preferentially N/O > 1. However, varying C/O and N/O is not sufficient to match the estimated best-fitting abundance of HCN, which lies outside the range that these two parameters alone can explore. In particular, there are a number of other important parameters that could and should also be varied, such as the Kzz strength, the atmospheric metallicity and the P–T profile. This full sweep of parameter space is, however, beyond the scope of this work, but provides a good basis for future investigations.

5 CONCLUSION

In this work, we present a new 1D diffusion and photochemistry code, named Levi, currently able to model gaseous chemistry in hot Jupiter atmospheres via solving of the continuity equation. We focus on H, C, O, and N chemistry, using a chemical network that contains over 1000 reactions and 150 species. Through inputs in the form of a P–T profile, an eddy diffusion profile, and a UV stellar flux, the code can calculate steady-state abundance profiles for all species over the desired pressure range.

Levi was validated against the disequilibrium chemistry models produced by Moses et al. (2011) and Venot et al. (2012). For HD 209458b, the abundance profiles of species that Levi produced matched closely with those of the other two models, with any differences being a consequence of the slightly differing input parameters. We discussed how the differing P–T profiles of HD 209458b and HD 189733b affected equilibrium chemistry in their atmosphere, and noted that, in chemical equilibrium, neither HCN nor NH3 would be detectable. Eddy-diffusion and photochemistry were also included, and we discussed how they affected the chemistry of the atmospheres, as well as how the chemistry of the two planets atmospheres compared to each other. Of particular, note was that the quenching caused by eddy-diffusion is capable of raising the abundance of NH3 and HCN to potentially detectable levels.

The influence of the parameter space of the C/O and N/O ratio on the abundance of various molecules at 100 mbar in the atmosphere of HD 189733b was investigated. It was found that species that contained C or O were strongly affected by variations in the C/O ratio, while other molecules were weakly affected, if at all. In general, only species that contain nitrogen are affected by changes to the N/O ratio. We also looked at how the fractions NH3/H2O and HCN/H2O varied with the atmospheric composition, an important fraction in determining whether it is possible to detect these molecules. We found that when trying to detect NH3 or HCN, the quantity of carbon in the atmosphere is much more important than the amount of nitrogen. At a solar N/O ratio, sub-solar carbon would make detecting either of these molecules near impossible; however with at least five times solar nitrogen in the atmosphere, NH3 and HCN could be detected with only 0.5 times and 0.8 times solar carbon, respectively. In an atmosphere with a deficit of nitrogen, less than two times solar carbon would make detecting either of these molecules possible.

Recent literature (Diamond-Lowe et al. 2014) has suggested that HD 209458b does not contain a thermal inversion. Other work (Madhusudhan et al. 2014a; MacDonald & Madhusudhan 2017a; Brogi & Line 2018; Hawker et al. 2018) has found evidence for CO, H2O, and HCN with best-fitting mixing ratios of approximately 10−4 ± 0.5, 10−5 ± 0.5, and 10−5, respectively, and with a minimum mixing ratio of HCN at 10−6.5. We applied our model to HD 209458b, using a new P–T profile without a thermal inversion, to discover which C/O and N/O ratios could best fit with these observations. A C/O > 1 resulted in a significant depletion of H2O. To obtain at least the minimum mixing ratio of HCN, a C/O ≳ 0.9 was required. A large N/O ratio, ten times solar or more, can help increase the HCN abundance closer to the best-fitting value, but it results in much smaller increases compared to increasing the C/O ratio. Further testing with both these parameters and others are outside the scope of this work, and have been left for future investigations.

SUPPORTING INFORMATION

Chemical_Network.txt

Included in the supplementary data is the full reaction network used in this work, along with the rate coefficients, as described in Section 2.2.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

Richard Hobbs and Nikku Madhusudhan acknowledge support from the Science and Technology Facilities Council (STFC), UK. Paul Rimmer acknowledges funding from the Simons Foundation SCOL awards 599634.

Footnotes

1

Named after the Latin Levis, meaning light.

REFERENCES

Agúndez

M.

,

Parmentier

V.

,

Venot

O.

,

Hersant

F.

,

Selsis

F.

,

2014

,

A&A

,

564

,

A73

Amundsen

D. S.

et al. .,

2016

,

A&A

,

595

,

A36

Asplund

M.

,

Grevesse

N.

,

Sauval

A. J.

,

Scott

P.

,

2009

,

ARA&A

,

47

,

481

Beaugé

C.

,

Nesvorný

D.

,

2012

,

ApJ

,

751

,

119

Blecic

J.

,

Harrington

J.

,

Bowman

M. O.

,

2016

,

ApJ

,

225

,

4

Brewer

J. M.

,

Fischer

D. A.

,

Madhusudhan

N.

,

2017

,

AJ

,

153

,

83

Brogi

M.

,

Line

M. R.

,

2018

,

AJ

,

157

,

114

Brown

T. M.

,

2001

,

ApJ

,

553

,

1006

Burcat

A.

,

Ruscic

B.

,

2005

,

Third millenium ideal gas and condensed phase thermochemical database for combustion with updates from active thermochemical tables

.

Argonne National Laboratory

,

Argonne, IL

Burrows

A.

,

Sharp

C. M.

,

1999

,

ApJ

,

512

,

843

Chapman

S.

,

Cowling

T.

,

1970

,

The Mathematical Theory of Non-uniform Gases

,

University Press

,

Cambridge

Cleeves

L. I.

,

Öberg

K. I.

,

Wilner

D. J.

,

Huang

J.

,

Loomis

R. A.

,

Andrews

S. M.

,

Guzman

V. V.

,

2018

,

ApJ

,

865

,

155

Cooper

C. S.

,

Showman

A. P.

,

2006

,

ApJ

,

649

,

1048

Dawson

R. I.

,

Johnson

J. A.

,

2018

,

ARA&A

,

56

,

175

Del Popolo

A.

,

Ercan

N.

,

Yeşilyurt

I. S.

,

2005

,

A&A

,

436

,

363

10.1051/0004-6361:20041561

Deming

D.

,

Brown

T. M.

,

Charbonneau

D.

,

Harrington

J.

,

Richardson

L. J.

,

ApJ

,

2005

,

622

,

1149

Diamond-Lowe

H.

,

Stevenson

K. B.

,

Bean

J. L.

,

Line

M. R.

,

Fortney

J. J.

,

2014

,

ApJ

,

796

,

66

Drummond

B.

,

Tremblin

P.

,

Baraffe

I.

,

Amundsen

D. S.

,

Mayne

N. J.

,

Venot

O.

,

Goyal

J.

,

2016

,

A&A

,

594

,

A69

Gandhi

S.

,

Madhusudhan

N.

,

2017

,

MNRAS

,

472

,

2334

Garcia Munoz

A.

,

2007

,

P&SS

,

55

,

1426

Gladstone

G. R.

,

Allen

M.

,

Yung

Y. L.

,

1996

,

Icar

,

119

,

1

Glassman

I.

,

Yetter

R.

,

Glumac

N.

,

2015

,

Combustion

.

Academic Press

Grillmair

C. J.

et al. .,

2008

,

Natur

,

456

,

767

Hands

T. O.

,

Alexander

R. D.

,

2016

,

MNRAS

,

456

,

4121

Hawker

G. A.

,

Madhusudhan

N.

,

Cabot

S. H. C.

,

Gandhi

S.

,

2018

,

ApJ

,

863

,

L11

Heng

K.

,

Tsai

S.-M.

,

2016

,

ApJ

,

829

,

104

Hu

R.

,

Seager

S.

,

Bains

W.

,

2012

,

ApJ

,

761

,

166

Hubeny

I.

,

Mihalas

D.

,

2014

,

Theory of Stellar Atmospheres

,

Princeton University Press

Huebner

W. F.

,

Carpenter

C. W.

,

1979

,

STIN

,

80

,

24243

Huebner

W. F.

,

Keady

J. J

, .

Lyon

S. P.

,

1992

,

Ap&SS

,

195

,

1

Huebner

W. F.

,

Mukherjee

J.

,

2015

,

P&SS

,

106

,

11

Keller-Rudek

H.

,

Moortgat

G. K.

,

Sander

R.

,

Sörensen

R.

,

2013

,

ESSD

,

5

,

365

Knutson

H.

et al. .,

2008

;

Liang

M.-C.

,

Parkinson

C. D.

,

Lee

A. Y.-T.

,

Yung

Y. L.

,

Seager

S.

,

2003

,

ApJ

,

596

,

L247

Lindemann

F. A.

,

Arrhenius

S.

,

Langmuir

I.

,

Dhar

N. R.

,

Perrin

J.

,

Lewis

W. C. McC.

,

1922

,

Trans. Faraday Soc.

,

17

,

598

Linsky

J. L.

,

Yang

H.

,

France

K.

,

Froning

C. S.

,

Green

J. C.

,

Stocke

J. T.

,

Osterman

S. N.

,

2010

,

ApJ

,

717

,

1291

Liou

K.

,

2002

,

An Introduction to Atmospheric Radiation

,

Academic Press

Liu

B.

,

Showman

A. P.

,

2013

,

ApJ

,

770

,

42

Lodders

K.

,

Fegley

B.

,

2002

,

Icar

,

155

,

393

MacDonald

R. J.

,

Madhusudhan

N.

,

2017a

,

MNRAS

,

469

,

1979

MacDonald

R. J.

,

Madhusudhan

N.

,

2017b

,

ApJ

,

850

,

L15

Madhusudhan

N.

,

2012

,

ApJ

,

758

,

36

Madhusudhan

N.

,

Burrows

A.

,

Currie

T.

,

2011

,

ApJ

,

737

,

34

Madhusudhan

N.

,

Amin

M. A.

,

Kennedy

G. M.

,

2014a

,

ApJ

,

794

,

L12

Madhusudhan

N.

,

Amin

M. A.

,

Kennedy

G. M.

,

2014b

,

ApJ

,

794

,

L12

Madhusudhan

N.

,

Agúndez

M.

,

Moses

J. I.

,

Hu

Y.

,

2016

,

Space Sci. Rev.

,

205

,

285

Mayne

N. J.

et al. .,

2014

,

A&A

,

561

,

A1

McGuire

B. A.

et al. .,

2018

, in

Murphy

E.

, ed.,

ASP Conf. Ser. Vol. 517, Science with a Next Generation Very Large Array

.

Astron. Soc. Pac

,

San Francisco

, p.

217

Miller-Ricci Kempton

E.

,

Zahnle

K.

,

Fortney

J. J.

,

2012

,

ApJ

,

745

,

3

Mordasini

C.

,

van Boekel

R.

,

Mollière

P.

,

Henning

T.

,

Benneke

B.

,

2016

,

ApJ

,

832

,

41

Moses

J. I.

et al. .,

2011

,

ApJ

,

737

,

15

Moses

J. I.

,

Madhusudhan

N.

,

Visscher

C.

,

Freedman

R. S.

,

2013

,

ApJ

,

763

,

25

Murray-Clay

R. A.

,

Chiang

E. I.

,

Murray

N.

,

2009

,

ApJ

,

693

,

23

Nelson

B. E.

,

Ford

E. B.

,

Rasio

F. A.

,

2017

,

AJ

,

154

,

106

Öberg

K. I.

,

Murray-Clay

R.

,

Bergin

E. A.

,

2011

,

ApJ

,

743

,

L16

Rimmer

P. B.

,

Helling

C.

,

2016

,

ApJS

,

224

,

9

Rosenbrock

H. H.

,

1963

,

Comput. J.

,

5

,

329

Sánchez

M. B.

,

de Elía

G. C.

,

Darriba

L. A.

,

2018

,

MNRAS

,

481

,

1281

Schwarz

H.

,

Brogi

M.

,

de Kok

R.

,

Birkby

J.

,

Snellen

I.

,

2015

,

A&A

,

576

,

A111

Seager

S.

,

Deming

D.

,

2010

,

ARA&A

,

48

,

631

Seager

S.

,

Sasselov

D. D.

,

2000

,

ApJ

,

537

,

916

Sneep

M.

,

Ubachs

W.

,

2005

,

J. Quant. Spec. Radiat. Transf.

,

92

,

293

Stevenson

K. B.

,

Bean

J. L.

,

Madhusudhan

N.

,

Harrington

J.

,

2014

,

ApJ

,

791

,

36

Swain

M. R.

et al. .,

2009

,

ApJ

,

704

,

1616

Tsai

S.-M.

,

Lyons

J. R.

,

Grosheintz

L.

,

Rimmer

P. B.

,

Kitzmann

D.

,

Heng

K.

,

2017

,

ApJ

,

228

,

20

Venot

O.

et al. .,

2013

,

A&A

,

551

,

A131

Venot

O.

et al. .,

2018

,

A&A

,

609

,

A34

Venot

O.

,

Hébrard

E.

,

Agúndez

M.

,

Dobrijevic

M.

,

Selsis

F.

,

Hersant

F.

,

Iro

N.

,

Bounaceur

R.

,

2012

,

A&A

,

546

,

A43

Verwer

J.

,

Spee

E.

,

Blom

J.

,

Hundsdorfer

W.

,

1999

,

SIAM J. Sci. Comput.

,

20

,

1456

Youngblood

A.

et al. .,

2016

,

ApJ

,

824

,

101

Zahnle

K.

,

Mac Low

M.-M.

,

Lodders

K.

,

Fegley

B.

Jr,

1995

,

GeoRL

,

22

,

1593

Zahnle

K.

,

Marley

M. S.

,

Freedman

R. S.

,

Lodders

K.

,

Fortney

J. J.

,

2009

,

ApJ

,

701

,

L20

Zahnle

K.

,

Marley

M. S.

,

Morley

C. V.

,

Moses

J. I.

,

2016

,

ApJ

,

824

,

137

© 2019 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

Supplementary data

Citations

Views

Altmetric

Metrics

Total Views 876

561 Pageviews

315 PDF Downloads

Since 5/1/2019

Citations

29 Web of Science

×

Email alerts

Citing articles via

• Latest

• Most Read

• Most Cited

More from Oxford Academic