A chemical kinetics code for modelling exoplanetary atmospheres

TL;DR

This paper introduces LEVI, a new chemical kinetics code for modeling exoplanet atmospheres, focusing on hot gas giants, and demonstrates its ability to relate atmospheric composition to physical conditions, especially C/O ratios.

Contribution

LEVI is a novel chemical kinetics code that models key atmospheric chemistry in exoplanets, enabling better interpretation of spectroscopic data.

Findings

Atmospheric C/O ratio strongly influences molecular detections.

Models are consistent with observed species in HD 209458b's atmosphere.

High N/O ratios can affect the C/O threshold for certain molecules.

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 super-solar N/O can decrease the C/O required for HCN and NH3 detection. Using a new 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 $\sim$ 1. A C/O $\gtrsim$ 0.9 (1.6 times solar) was required to meet the minimum reported value for HCN, while a C/O $\lesssim$ 1 (1.8 times solar) was required to fit the nominal H2O abundance.

Figures

Tables

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

Model	Network Basis	Diffusion	Photochemistry	Elements	Reactions	Molecules
Liang et al. (2003)	From Laboratory Measurements 222These rates were measured for a low temperature Jovian Model.	Yes	Yes	H/C/O	253	Unknown
Zahnle et al. (2009)	Zahnle et al. (1995)	Yes	Yes	H/C/O/N/S	507	49
Moses et al. (2011)	Gladstone et al. (1996)	Yes	Yes	H/C/O/N	1600	90
Venot et al. (2012)	From Combustion Studies	Yes	Yes	H/C/O/N	1918	104
Hu et al. (2012)	New	Yes	Yes	H/C/O/N/S	800	111
Rimmer & Helling (2016)	Stand2015	Yes	Yes	H/C/O/N/Metals	2980	162
Drummond et al. (2016)	Venot et al. (2012)	Yes	Yes	H/C/O/N	1918	104
Tsai et al. (2017)	Glassman et al. (2015)	Yes	No	H/C/O	300	29
This work	Stand2018	Yes	Yes	H/C/O/N	2000	150

Equations

$$\frac{\partial n_i}{\partial t}={\mathcal{P}}_i-{\mathcal{L}}_i-\frac{\partial {\Phi }_i}{\partial z},$$

$${\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{dT}{dz}\right),$$

$$D_i=\frac{3}{8n_t(\frac{1}{2}[d_i+d_t]{)}^2}{\left(\frac{k_{b}T(m_i+m_t)}{2\pi m_i{m}_t}\right)}^{1/2},$$

$$n_t=\frac{P(z)}{k_{b}T(z)},$$

$${\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 }}).$$

$$F(\lambda ,z)=F_0(\lambda )e^{-\tau (\lambda ,z)/{\mu }_0},$$

$$\tau (\lambda ,z)={\int }_z^{\infty }{n}_i(z^{\prime }){\sigma }_i(\lambda )d{z}^{\prime },$$

$${\sigma }_r(\lambda )=C_r\frac{8{\pi }^3(m_r(\lambda {)}^2-1{)}^2}{3{\lambda }^4{N}_s},$$

$$A+B\to C+D,$$

$$k=\alpha {\left(\frac{T}{300}\right)}^{\beta }\exp \left(-\frac{\gamma }{T}\right),$$

$$A+M\to C+D+M,$$

$$A+B+M\to C+M,$$

$$k_0$$

$$k_{\infty }$$

$$k=\frac{k_0[M]}{1+\frac{k_0[M]}{k_{\infty }}},$$

$$\frac{k_f}{k_r}=K_c{\left(\frac{k_{b}T}{P_0}\right)}^{-\Delta \nu },$$

$$K_c=\exp \left(-\frac{\Delta G^o}{\mathcal{R}T}\right),$$

$$\Delta G^o=\Delta H^o-T\Delta S^o,$$

$$\frac{H^o}{\mathcal{R}T}$$

$$\frac{S^o}{\mathcal{R}}$$

$$k=t_f\int {\sigma }_a(\lambda )F(\lambda ,z)d\lambda$$

$$\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},$$

$$\begin{gathered} {\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} \\ +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], \end{gathered}$$

$$\begin{gathered} \frac{\Delta n_{i,j}}{\Delta t}={\mathcal{P}}_{i,j}-{\mathcal{L}}_{i,j}+\frac{1}{\Delta z_{j-1/2}}[(k_{i,j+1/2}^{-}\frac{n_{t,j+1/2}}{n_{t,j+1}})n_{i,j+1}- \\ (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}})n_{i,j}+(k_{i,j-1/2}^{+}\frac{n_{t,j-1/2}}{n_{t,j-1}})n_{i,j-1}], \end{gathered}$$

$$\begin{gathered} k_{j+1/2}^{\pm }=\frac{K_{zz,j+1/2}+D_{i,j+1/2}}{\Delta z_j}\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{gathered}$$

$$P(z)=P_0{e}^{\frac{-m(z)gz}{k_{b}T(z)}},$$

$$\frac{dn}{dt}=f(n),$$

$$n_{k+1}$$

$$(\text{I}-\gamma \Delta t\text{J})g_1$$

$$(\text{I}-\gamma \Delta t\text{J})g_2$$

$$\text{J}\equiv \frac{\partial f}{\partial n}{|}_{n_k},$$

$$J=\left[\begin{array}{ccccccccc}{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\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮\\ 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{array}\right]$$

Full text

††pubyear: 2019††pagerange: A chemical kinetics code for modelling exoplanetary atmospheres–A chemical kinetics code for modelling exoplanetary atmospheres

1 Introduction

Identifying and characterising 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 et al. 2017; Dawson & Johnson 2018; McGuire et al. 2018). Both the location of where giant planets form and possible later migration through the proto-planetary disk may leave signs in the compositions of their atmospheres (e.g. Öberg et al. 2011; Moses et al. 2013; Madhusudhan et al. 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 et al. 2005; Hands & Alexander 2016; Sánchez et al. 2018).

A key tool for characterising 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$ - $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 \ceC/\ceO and \ceN/\ceO ratios to the retrieved atmospheric abundances. Since the composition of planetary atmospheres may reflect when and where in a disk planets 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 $\gtrsim$ 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 photo-chemistry (Cooper & Showman 2006; Moses et al. 2011; Moses et al. 2013).

We report a new code for 1-D modelling of chemistry in exoplanetary atmospheres, named Levi. 111Named after the Latin Levis, meaning light. 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 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 \ceH, \ceC, \ceO, \ceN or \ceHe. 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 generalised 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 et al. 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 et al. 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. These models all work similarly, using a pressure - temperature (P-T) and eddy-diffusion ($K_{zz}$) 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 \ceH-\ceC-\ceO chemistry for a select set of species in highly irradiated atmospheres. Zahnle et al. (2009) investigated the effects of photochemistry on sulfur 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 \ceH-\ceC-\ceN-\ceO photochemistry on hot Jupiters. Agúndez et al. (2014) adapted the work of Venot et al. (2012) to report a 2-D 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.

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 characterised planets (e.g., Seager & Sasselov 2000; Moses et al. 2013; Lewis et al. 2013; Liu & Showman 2013; Mayne et al. 2014; Schwarz et al. 2015; Amundsen et al. 2016; Brewer et al. 2017), and several previous chemical kinetics 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 it’s 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 \ceC/\ceO and \ceN/\ceO 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 1-D 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 pre-determined 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 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 break-down 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 pressure-temperature (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 further explore the atmospheric response to the presence of particular species, and to improve the computational runtime.

2.1 Model Setup

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 one-dimensional continuity equations,

[TABLE]

where $n_{i}$ (${\mathrm{m}}^{-3}$) is the number density of species $i$, with $i=1,...,N_{i}$, with $N_{i}$ being the total number of species. $\mathcal{P}_{i}$ (${\mathrm{m}}^{-3}\text{\,}{\mathrm{s}}^{-1}$) and $\mathcal{L}_{i}$ (${\mathrm{m}}^{-3}\text{\,}{\mathrm{s}}^{-1}$) are the production and loss rates of the species $i$. $\partial t$ ($\mathrm{s}$) and $\partial z$ ($\mathrm{m}$) are the infinitesimal time step and altitude step respectively. $\Phi_{i}$ (${\mathrm{m}}^{-2}\text{\,}{\mathrm{s}}^{-1}$) is the upward vertical flux of the species, given by,

[TABLE]

where $X_{i}$ and $n_{t}$ (${\mathrm{m}}^{-3}$) are the mixing ratio and total number density of molecules such that $n_{i}=X_{i}n_{t}$. The eddy-diffusion coefficient, $K_{zz}$ (${\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$), approximates the rate of vertical transport and $D_{i}$ (${\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$) is the molecular diffusion coefficient of species i. $H_{0}$ ($\mathrm{m}$) is the mean scale height, $H$ ($\mathrm{m}$) is the molecular scale height, $T$ (K) is the temperature, and $\alpha_{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,

[TABLE]

where $m_{i}$ (kg) and $d_{i}$ (m) are the mass and diameter of species $i$, respectively, and $m_{t}$ (kg) and $d_{t}$ (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, $\alpha_{T,i}\approx-0.1(1-n_{i}/n_{t})$; for \ceHe, $\alpha_{T,i}\approx 0.145(1-n_{i}/n_{t})$; for \ceH2, $\alpha_{T,i}\approx-0.38$; and for all other molecules $\alpha_{T,i}\approx 0.25$. Under the assumption that the atmosphere can be approximated to an ideal gas, the total number density, $n_{t}$ (${\mathrm{m}}^{-3}$), is

[TABLE]

where $P$ ($\mathrm{Pa}$) and $T$ ($\mathrm{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 which include said species, (e.g. $n_{i}\ce{+}n_{i^{\prime}}\rightleftharpoons n_{i^{\prime\prime}}\ce{+}n_{i^{\prime\prime\prime}}$), then finding the rate at which $n_{i}$ 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,

[TABLE]

2.1.2 Atmospheric Profile

The physical characteristics of the atmospheres Levi simulates have a profound impact on their chemistry. This is realised through a pressure-temperature (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 Figure 1(a). We also consider an isothermal profile for reference.

The P-T and $K_{zz}$ 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 photo-dissociated into a variety of species, including radicals. Since the rate of photo-dissociation 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,

[TABLE]

where $F_{0}$ (photons ${\mathrm{m}}^{-2}\text{\,}{\mathrm{s}}^{-1}$ $\mathrm{nm^{-1}}$) is the flux at the top of the atmosphere (where by definition $\tau=0$), $\mu_{0}$ is angle of the path of sunlight, chosen to be $57.3^{o}$ (Hu et al. 2012). $\tau$ is the total optical depth, as a function of wavelength, integrated over all the atmosphere that the light passes through, defined as,

[TABLE]

where $\sigma_{i}$ is the wavelength dependent cross-section of species. $\sigma_{i}$ is found from the sum of the molecular absorption cross-section, $\sigma_{a}(\lambda)$, and the Rayleigh scattering cross section (Liou 2002),

[TABLE]

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 \ceCO2, 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 \ceH2.

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_{\oplus}$) and the planet and its host star ($d_{r}$), $(d_{\oplus}/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 Figure 2.

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 \ceH/\ceC/\ceN/\ceO/\ceHe 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,

[TABLE]

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

[TABLE]

where $k$ (${\mathrm{m}}^{3}\text{\,}{\mathrm{s}}^{-1}$) is the rate coefficient, $T$ (K) is the temperature and $\alpha$ (${\mathrm{m}}^{3}\text{\,}{\mathrm{s}}^{-1}$), $\beta$ and $\gamma$ (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,

[TABLE]

or combination reactions,

[TABLE]

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,

[TABLE]

with $k_{0}$ ($\mathrm{m}^{6}\,\mathrm{s}^{-1}$ for combination or $\mathrm{m}^{3}\,\mathrm{s}^{-1}$ for decomposition) as the low-pressure rate constant and $k_{\infty}$ ($\mathrm{m}^{3}\,\mathrm{s}^{-1}$ for combination or $\mathrm{s}^{-1}$ for decomposition) as the high-pressure rate constant. These can be combined to form an effective two-body reaction rate,

[TABLE]

where $[M]$ (${\mathrm{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,

[TABLE]

in which $k_{f}$ and $k_{r}$ are the rate of the forward and reverse reactions, $k_{b}$ is the Boltzmann constant, $\Delta\nu$ 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,

[TABLE]

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

[TABLE]

with $\Delta H^{o}$ and $\Delta S^{o}$ as 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,

[TABLE]

The values for these polynomials were drawn from Burcat 333http://garfield.chem.elte.hu/Burcat/burcat.html 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 photo-dissociation for photochemically active molecules is also calculated. The molecular absorption cross-sections, $\sigma_{a}$ (${\mathrm{m}}^{2}$), for these molecules were drawn primarily from PhIDRates 444http://phidrates.space.swri.edu/, apart from \ceC4H2 and \ceC4H4, which were not in this database. These two molecules used the cross-sections from the MPI-Mainz-UV-VIS Spectral Atlas of Gaseous Molecules 555www.uv-vis-spectral-atlas-mainz.org. 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 \ceCO2 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 emphasises the need for more high temperature cross-sections to accurately compute the photo-chemistry occuring 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 photo-dissociation rate is then,

[TABLE]

with $F(\lambda,z)$ as the Actinic Flux as described in 2.1.3, and $t_{f}$ 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. Photo-dissociation 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. Photo-dissociation 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 discretised into a form that can be solved numerically.

2.3.1 Discretisation

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

[TABLE]

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,…,$N_{z}$, with $N_{z}$ 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 $\Delta z$. Equation 2 can be similarly discretised to

[TABLE]

in which any $Y_{j\pm 1/2}$ has been approximated to $(Y_{j\pm 1}+Y_{j})/2$. $m_{t,j}$ is the mean mass of atmosphere layer j, and $m_{i}$ is the mass of species i, $g$ is the gravitational acceleration, and $k_{b}$ is the Boltzmann constant. Now, it is possible to combine Equation 20 and 21 to form a set of ordinary differential equations with time as the independent variable,

[TABLE]

in which

[TABLE]

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

[TABLE]

with $m$ as the mean mass of the molecules in an atmosphere layer, found by $\sum\limits_{i}m_{i}n_{i}$, where $m_{i}$ 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 $\mathrm{m}\text{\,}{\mathrm{s}}^{-2}$ for an isothermal profile, 21.4 $\mathrm{m}\text{\,}{\mathrm{s}}^{-2}$ for HD 189733b and 9.36 $\mathrm{m}\text{\,}{\mathrm{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 which have sub-microsecond timescales and diffusion which can occur over timescales of hundreds of years. The general differential equation is,

[TABLE]

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,

[TABLE]

where $\Delta t$ is the size of the time-step, $\gamma$ 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

[TABLE]

is the Jacobian matrix evaluated at $n_{k}$.

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 $N_{i}N_{z}$, 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 $N_{i}$ that describe the interactions between species in the same atmosphere layer. The off-diagonal blocks, A and C, also matrices of size $N_{i}$, represent the effects of adjacent atmosphere layers via diffusion. The subscripts for the A, B and C matrices show which atmosphere layer they represent.

[TABLE]

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

[TABLE]

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 generalisation to block tridiagonal by treating each matrix element in the inversion method as a full matrix. As a result, only the three $N_{i}$ size matrices are ever kept in the memory, which is much less expensive than trying to invert the complete $N_{i}N_{z}$ 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, $\epsilon$ is sufficiently small.

The truncation error is found by calculating the difference between second-order solution and the first-order solution, $n^{1}_{k+1}=n_{k}+\Delta tg_{1}$, such that,

[TABLE]

Since the timescales 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,

[TABLE]

2.4 Initial and Boundary Conditions

2.4.1 Initial Conditions

In the present 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 \ceC/\ceO and \ceN/\ceO 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_{\ce{H}}=1.707$, $X_{\ce{He}}=0.145$, $X_{\ce{C}}=4.594\times 10^{-4}$, $X_{\ce{N}}=1.154\times 10^{-4}$ and $X_{\ce{O}}=8.359\times 10^{-4}$. When calculating solar abundances all \ceH is assumed to be in the form of \ceH2, thus leading to a $X_{\ce{H}}$ that is greater than one. This gives a \ceC/\ceO ratio of 0.55 and a \ceN/\ceO ratio of 0.138.

When choosing the initial conditions, conservation of mass is ensured by summing over the product of the species present ($n_{i}$) and the number of atoms of $x$ contained by $n_{i}$ ($A_{i,x}$),

[TABLE]

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 Figure 3, with the initial species sets being: (1) \ceH2, \ceCO, \ceOH and \ceNH3; (2) \ceH2, \ceCH4, \ceH2O and \ceN2; and (3) atomic \ceH, \ceC, \ceO, and \ceN. 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 initialise 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. Plots b), c) and d) of Figure 3 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.

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 (Linsky et al. 2010; Murray-Clay et al. 2009). For a boundary in chemical equilibrium, the assumption is that layers deeper in the atmosphere have a chemical timescale much shorter than the dynamical timescale, 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 Figure 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 \ceN2. To better quantify this change, the right plot in Figure 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 percent, however closer to the lower boundary it can increase up to 10%. 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% 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.

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 prevoius codes (Moses et al. 2011 and 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 2.2.

3.1 Validation

We chose HD 209458b as one of the best candidates for benchmarking Levi, due to both the well characterised 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) (Figure 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 \ceH2, \ceH and \ceH2O (Figure 5a), \ceCO, \ceCH4 and \ceCO2 (Figure 5b), and \ceN2, \ceNH3 and \ceHCN (Figure 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.

Looking at the comparison plots in Figure 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% of the oxygen abundance to be sequestered in silicates.

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

[TABLE]

The observed difference in the quenching location of \ceCH4 arises for similar reasons. In the network used by Levi, there are several pathways that convert \ceCH4 into \ceCO, one of which has a highly uncertain reaction rate (Rimmer & Helling 2016), and can explain the nearly 1 dex and 0.5 dex higher quenching abundance of \ceCH4 seen in Moses et al. (2011) and Venot et al. (2012), respectively, than in Levi. The large differences in abundance seen in \ceHCN 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 \ceHCN.

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 vs. Diffusion vs. 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 Figure 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 1-D 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).

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 \ceCH4 and \ceH2O dominate over \ceCO, and \ceNH3 approaches the abundance of \ceN2, while higher in the atmosphere this trend is reversed (Heng & Tsai 2016). When solving for equilibrium analytically, this same trend is found. In the observable region of the atmosphere, approximately $1$ - $10^{-4}$ bar for emission spectrography, \ceCO is the dominant carbon species, with most of the remaining oxygen ending up in \ceH2O, and \ceN2 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 Figure 6. Notable for comparison later on are the abundances of \ceHCN and \ceNH3 in the observable region of the atmosphere; in this equilibrium scenario they are very minor constituents, far below detectable limits (\ceHCN or \ceNH3 abundances at least 1% of the \ceH2O abundance) for solar composition (MacDonald & Madhusudhan 2017b).

3.2.2 Diffusion

The inclusion of diffuion, the dashed lines in Figure 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.

Figure 6 shows a number of species being affected by quenching in observable regions of the atmosphere, including \ceCH4, \ceHCN and \ceNH3. \ceNH3 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 timescale compared to the mixing timescale, 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 the planet’s upper atmosphere not having a thermal inversion below the thermosphere, so the chemical timescale monotonically increases with decreasing pressure. On the other hand, HD 209458b does have a thermal inversion, so the chemical timescale 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 Figure 6 show how the inclusion of photochemistry affects the atmosphere. The photolysis of \ceCO, \ceH2O, \ceN2, \ceNH3 and \ceCH4 drive most of the photochemical reactions that occur in the upper atmosphere.

The destruction of \ceH2O by a UV photon into \ceOH and \ceH, as deep as 0.1 bar, sets up the reaction that turns \ceH2 into \ceH:

[TABLE]

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 Figure 6a and 6b, \ceH2O 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 \ceCO can produce \ceC and \ceO radicals that can also contribute to the destruction of \ceH2. Like \ceH2O, \ceCO is also being replenished quickly, in this case by a reaction between \ceCH4 and \ceH2O:

[TABLE]

However, as Figure 6d shows, on HD 189733b the formation of \ceCO is not efficient enough to prevent loss of \ceCO above $10^{-5}$ bar and so some of the carbon can end up in other species, e.g., \ceCO2 and \ceHCN. On HD 209458b, the temperature, and thus the reaction rates, are high enough to keep replenishing \ceCO until pressures as low as $10^{-7}$ bar.

The \ceC radicals from the photolysis of \ceCO often end up producing \ceCH4:

[TABLE]

which leads to the increase of methane seen on HD 209458b. \ceCH4 is, in turn, susceptible to destruction by \ceH radicals, thus leading to an overall decrease in methane if sufficient \ceH is produced by the photo-dissociation of \ceH2O, as discussed earlier, and as can be seen at $10^{-5}$ bar on HD 189733b.

\ce

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

\ce

C2H2 is an important byproduct of photochemistry. The carbon atoms released from dissociation of \ceCH4 or \ceCO are reprocessed to form \ceC2H2. This is the result of a highly efficient pathway:

[TABLE]

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

Photo-chemistry can also speed up the reaction converting \ceNH3 to \ceHCN,

[TABLE]

causing a significant increase in the amount of \ceHCN. On HD 209458b this proceeds to the point where \ceHCN 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 \ceNH3 seen in the upper atmosphere of HD 189733b, compared to models without photochemistry. Although loss processes exist for \ceHCN, the products are normally not very stable, and often form \ceHCN upon their destruction. The result is that destruction reactions are inefficient in regulating \ceHCN’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 \ceC/\ceO ratio and \ceN/\ceO ratio, and consider how it may effect the detectability of certain species in exoplanet atmospheres. Apart from where stated otherwise, these models use the P-T profile and $K_{zz}$ 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 \ceC/\ceO and \ceN/\ceO ratios that would best fit with new evidence of multiple species on HD 209458b.

4.1 The \ceC/\ceO and \ceN/\ceO 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 \ceC/\ceO ratios substantially different from solar values (Madhusudhan et al. 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 it’s past (Öberg et al. 2011, Madhusudhan et al. 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 \ceC/\ceO and \ceN/\ceO 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 \ceC/\ceO ratio varying between 0.055 and 5.5, and an \ceN/\ceO 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 $K_{zz}$ profile of HD 189733b from Figures 1(a) and 1(b).

In Figure 7 and 8 maps of the abundance of \ceC and \ceO species are shown as a function of atmospheric composition. None of the species show any strong dependence on the amount of \ceN in the atmosphere. In hydrogen dominated atmospheres, the abundances of \ceH2O, \ceCH4 and \ceCO are primarily determined by the equation:

[TABLE]

which at 100 mbar in the atmosphere of HD 189733b favours the formation of CO. At solar \ceC, a \ceC/\ceO ratio of 0.55, this means that \ceCO is the primary carbon carrier and \ceH2O contains most of the remaining oxygen. Since there is nearly twice as much oxygen as carbon in the atmosphere, \ceCO and \ceH2O have very similar abundances. For sub-solar carbon, there is less \ceCO, and thus less oxygen bound to \ceCO, resulting in an increase of \ceH2O. Once \ceC/\ceO > 1 (1.8$\times$ solar carbon), oxygen, not carbon, is now the limiting factor in producing \ceCO, 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 \ceCH4, resulting in a sharp increase in its abundance. \ceCO2 is also affected indirectly by equation 39; while there is excess oxygen not bound to \ceCO, its abundance is directly related to the amount of carbon, however once \ceC/\ceO > 1, then the majority of the oxygen is bound to \ceCO, and the abundance of \ceCO2 drops off rapidly. Larger hydrocarbons can also be very strongly affected by the variation in carbon, with the abundance of \ceC2H2 changing by more than eight orders of magnitude for only two orders of magnitude change in carbon.

In Figure 9 the variation of the three major nitrogen species within the parameter space is shown. The abundance of \ceN2 has no dependence on carbon, while both \ceNH3 and \ceHCN do. This is because of the fast reaction,

[TABLE]

thus causing the abundance of \ceNH3 and \ceHCN to have a dependence on the amount of carbon in the atmosphere. This dependence is particularly strong when \ceC/\ceO > 1, since there can now be an excess of \ceCH4, without it being converted into \ceCO. More \ceCH4 therefore leads to more \ceNH3 being converted into \ceHCN. As expected, all the nitrogen species have a strong dependence on the amount of nitrogen in the atmosphere, though \ceHCN 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 \ceC/\ceO ratio are seen in Madhusudhan (2012).

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

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 it’s 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 \ceCO, \ceH2O and \ceHCN in the atmosphere of HD 209458b (Hawker et al. 2018). Hawker et al. (2018) find the best fit mixing abundance of \ceHCN to be $10^{-5}$, with a minimum mixing ratio of $10^{-6.5}$. They also find \ceH2O and \ceCO mixing ratios consistent with previous values of approximately $10^{-5\pm 0.5}$ and $10^{-4\pm 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 (Figure 11). The same $K_{zz}$ profile and UV flux for HD 209458b from earlier are used, seen in Figure 1(b) and Figure 2 respectively.

In Figure 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 \ceHCN or \ceCH4, are present in greater quantities. For \ceHCN and \ceCH4 this results in an increase in abundance of nearly two orders of magnitude at $10^{-3}$ bar.

In Figure 13 the abundance profiles for \ceCO, \ceHCN and \ceH2O are displayed for a range of \ceC/\ceO 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 \ceCO provides little information since across all \ceC/\ceO ratios being modelled it stays at approximately $10^{-3}$, consistent with the expected value. \ceHCN provides more information, with a \ceC/\ceO $\gtrsim$ 0.9 required to pass the minimum estimated mixing ratio of $10^{-6.5}$, and a \ceC/\ceO = 1.2 is required for \ceHCN to match the best fit abundance of $10^{-5}$. However, \ceH2O provides an upper limit of \ceC/\ceO $\approx$ 1; higher \ceC/\ceO than this quickly depletes \ceH2O below $10^{-6}$. As such, it is clear that alterations to the \ceC/\ceO ratio by itself may not be sufficient to explain the observed abundances.

In Figure 13 we also included two models in which we increased the \ceN/\ceO ratio by two times (\ceN/\ceO = 0.28) and ten times (\ceN/\ceO=1.4) the solar value, to discover the effect this has on the abundance of these molecules. For these models, we kept \ceC/\ceO = 0.9. Increasing the amount of nitrogen in the atmosphere by two times solar increases the \ceHCN 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 \ceHCN since the amount of free carbon is the limiting factor beyond this. The \ceN/\ceO ratio has no significant effect on the abundance of either \ceCO or \ceH2O.

Through comparison of our model and the expected molecular abundances of HD 209458b we can produce some initial constraints for the \ceC/\ceO and \ceN/\ceO ratio of this hot Jupiter: a 0.9 $\lesssim$ \ceC/\ceO $\gtrsim$ 1, and preferentially \ceN/\ceO > 1. However, varying \ceC/\ceO and \ceN/\ceO is not sufficient to match the estimated best fit abundance of \ceHCN, 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 $K_{zz}$ 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 one-dimensional 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 \ceH, \ceC, \ceO and \ceN chemistry, using a chemical network that contains over 1000 reactions and 150 species. Through inputs in the form of a pressure - temperature 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 \ceHCN nor \ceNH3 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 \ceNH3 and \ceHCN to potentially detectable levels.

The influence of the parameter space of the \ceC/\ceO and \ceN/\ceO 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 \ceC or \ceO were strongly affected by variations in the \ceC/\ceO ratio, while other molecules were weakly affected, if at all. In general, only species that contain nitrogen are affected by changes to the \ceN/\ceO ratio. We also looked at how the fractions \ceNH3/\ceH2O and \ceHCN/\ceH2O 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 \ceNH3 or \ceHCN, the quantity of carbon in the atmosphere is much more important than the amount of nitrogen. At a solar \ceN/\ceO ratio, sub-solar carbon would make detecting either of these molecules near impossible, however with at least 5 times solar nitrogen in the atmosphere, \ceNH3 and \ceHCN 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 \ceCO, \ceH2O and \ceHCN with best fit mixing ratios of approximately $10^{-4\pm 0.5}$, $10^{-5\pm 0.5}$ and $10^{-5}$ respectively, and with a minimum mixing ratio of \ceHCN at $10^{-6.5}$. We applied our model to HD 209458b, using a new P-T profile without a thermal inversion, to discover which \ceC/\ceO and \ceN/\ceO ratios could best fit with these observations. A \ceC/\ceO > 1 resulted in a significant depletion of \ceH2O. To obtain at least the minimum mixing ratio of \ceHCN, a \ceC/\ceO $\gtrsim$ 0.9 was required. A large \ceN/\ceO ratio, ten times solar or more, can help increase the \ceHCN abundance closer to the best fit value, but it results in much smaller increases compared to increasing the \ceC/\ceO ratio. Further testing with both these parameters and others are outside the scope of this work, and have been left for future investigations.

Acknowledgements

R.H. and N.M. acknowledge support from the UK Science and Technology Facilities Council (STFC). PBR acknowledges funding from the Simons Foundation SCOL awards 599634.