© K. Molaverdikhani et al. 2020

1 Introduction

Ultra-hot Jupiters (UHJs) are highly irradiated, tidally locked planets with extreme day-to-night temperature differences (e.g., Delrez et al. 2017; Parmentier et al. 2018; Bell & Cowan 2018; Kreidberg et al. 2018; Komacek & Tan 2018; Lothringer et al. 2018; Lothringer & Barman 2019; Molaverdikhani et al. 2019a; Tan & Komacek 2019; Arcangeli et al. 2019; Mansfield et al. 2020; Wong et al. 2019). A day-to-night temperature difference of ≥2000 K on, for example, WASP-18b and HAT-P-7b results in different chemical regimes. The hot dayside atmosphere is predominately composed of H/He and is largely cloud-free, while the cold nightside is composed of a H₂/He atmosphere where cloud particles can form in abundance (Helling et al. 2019a,b).

The extreme insolation received by UHJs, resulting in large atmospheric scale heights, has rendered them prime targets for spectroscopic studies. The most extreme of these planets is KELT-9b, which is the hottest known exoplanet with an equilibrium temperature of ~4000 K. High-resolution HARPS-North observations of KELT-9b indicate an abundance of ions, such as Fe⁺ and Ti⁺, in its upper atmosphere (Hoeijmakers et al. 2018). The nondetection of atomic Ti, and less significant detection of Fe, point toward a highly ionized environment. This proliferation of ions is consistent with expectations from thermochemical equilibrium at these temperatures (Kitzmann et al. 2018). The somewhat cooler world MASCARA-2b/KELT-20b (T_eq ~ 2300 K) also shows signs of ions, namely Ca⁺, and Fe⁺, in addition to Na in its transmission spectrum (Casasayas-Barris et al. 2019). Successes such as these across a wide temperature regime demonstrate the observational accessibility of these atmospheric systems.

The progress in detecting spectral features in spectra of UHJs motivates detailed theoretical modeling of their atmospheres. 1D simulations, for instance, suggest that these atmospheres are in thermochemical equilibrium at the pressures probed by low-resolution spectroscopy, that is, higher than 1 mbar (Kitzmann et al. 2018). However, hydrodynamic transport processes (advection or winds), turbulent eddies, and interactions with a radiation field or charged particles (cosmic rays) can affect the local gas phase composition if they favor (or enable) a certain reaction path over its backward or forward reaction.

It has been shown that turbulent and molecular diffusion, photochemical processes, as well as interaction with free electrons may alter the chemical composition of planetary atmospheres from thermochemical equilibrium to a disequilibrium chemistry state (see, e.g., Line et al. 2010; Moses et al. 2011; Venot et al. 2012; Miguel & Kaltenegger 2013; Moses 2014; Rimmer & Helling 2016; Tsai et al. 2017; Helling & Rimmer 2019). While previous studies have typically been conducted on colder planets, Kitzmann et al. (2018) reported that the chemical state of the atmosphere of KELT-9b also deviates from thermochemical equilibrium at pressures lower than 1 mbar in a 1D setup. Moreover, according to their simulations, disequilibrium processes drive some of the atomic and molecular abundances away from their chemical equilibrium states at pressures even higher than 1 mbar. For instance, assuming an inversion in the atmosphere of KELT-9b, they report a pressure of ~0.03 bar at which or lower than which H₂, H₂O, and CO₂ abundances deviate from their thermochemical equilibrium values. CO and HCN abundances also deviate from their thermochemical equilibrium values at around a similar pressure, 0.02 bar. Transmission spectroscopy is sensitive to these pressures. This suggests that disequilibrium processes might leave observable fingerprints on the chemical composition of highly irradiated planets as well.

The 3D nature of planets and their resulting colder nightside could also enhance the effect of disequilibrium processes by extending the quenching level of major opacity species to higher pressures at the terminators and nightside. Consequently, 3D studies are needed to take into account local processes (such as disequilibrium chemistry) and connect the probed regions of the atmosphere to its global properties. Here, we begin to assess the nature of disequilibrium processes in UHJs using full 3D modeling of the governing chemical and physical processes. The chosen target of our investigation is the UHJ HAT-P-7b (T_eq ~ 2100 K). We began this exercise by studying global cloud formation and its effect on local elemental abundances in Helling et al. (2019b, hereafter Paper I). We continue our investigation by addressing the role of kinetic processes in altering the 3D gas-phase composition from thermochemical equilibrium expectations.

We note, however, that models linking kinetic gas-phase chemistry to dynamic atmosphere simulations are quite time-consuming. Relaxation schemes are therefore commonly applied to omit the long computational times required to ensure independence of the solution from the initial conditions (see, e.g., Cooper & Showman 2006; Drummond et al. 2018a,b; Mendonça et al. 2018). In such schemes, it is assumed that if the timescale for hydrodynamic transport processes becomes longer than the chemical timescale, then the local atmosphere will achieve chemical equilibrium (wherein both forward and backward reactions occur with a similar rate). This mostly occurs at hot and high-pressure atmospheric levels or in hydrodynamically inactive regions. The atmosphere will be driven out of thermochemical equilibrium if the hydrodynamic timescales (also called “mixing timescales”) become shorter than the chemical timescales. This occurs, for example, if a cell of material is transported more rapidly than its thermodynamic state can adjust to its new surrounding (e.g., Prinn & Barshay 1977; Griffith & Yelle 1999; Saumon et al. 2006; Cooper & Showman 2006; Drummond et al. 2018b). Hence, the composition of the cell may remain that of a hotter gas as long as no decay reaction is fast enough (e.g., N₂ /NH₃ and CO/CH₄). This has been demonstrated for brown dwarfs for example by Saumon et al. (2000). The effect of such a dynamically induced gas disequilibrium on phase curves has been demonstrated for example in Steinrueck et al. (2019). For instance, in wavelength bands with significant methane absorption, the amplitude of phase curves can be dramatically reduced. This effect is expected to be more significant for planets on which the nightside becomes cold enough so that under equilibrium chemistry, methane becomes the dominant nightside carbon-bearing species while CO prevails on the dayside (Steinrueck et al. 2019). We note, however, that CH₄ maybe locally destroyed by lightning (Bailey et al. 2014).

Transport-induced chemical disequilibrium mostly occurs for gaseous species at higher altitude atmospheric conditions. Here, it may lead to a chemical stratification as each species reacts according to its own timescale, and each species diffuses on its own timescale, in principle. The pressure level at which this transition between thermodynamic equilibrium and kinetic disequilibrium takes place is called the quenching point or quenching level (Smith 1998). The determination of the quenching level varies between models, which apply, for example, a constant diffusion constant for all (gas and cloud) species (Zhang & Showman 2018a,b). We further note that quenching levels can only be studied for molecules with known kinetic forward and backward rates, and with known diffusion constants.

The relative gas number densities, that is, mixing ratios, of the affected species at pressures lower than the pressure of the quenching point (i.e., higher altitudes) are believed to remain equal to the equilibrium abundances at the quenching point, or so-called quenched abundances (see, e.g., Prinn & Barshay 1977; Moses et al. 2011; Venot et al. 2018a). This uniform (quenched) mixing ratio profile has frequently been used in the atmospheric retrievals as a simplified approach to increase the computational speed (see, e.g., Madhusudhan & Seager 2011; Tsiaras et al. 2018; Pinhas et al. 2018). We propose a fast zeroth-ordermethod (based on the results of Venot et al. 2018a) to estimate the quenching levels and compare it with the gas-phase kinetic estimations.

In addition, under different circumstances and for different species, the abundance above the quenched level can deviate from a constant value in a chemical kinetic context (see, e.g., Molaverdikhani et al. 2019b). If significant, this deviation might pose a problem in the retrieved abundances, as shown by Changeat et al. (2019), for example. We discuss the validity of constant mixing ratio approximation for the case of HAT-P-7b in Sect. 3.2.

We therefore aim to quantitatively assess the importance of kinetic gas-phase chemical processes that are driven by hydrodynamic motion or by photochemistry, and map them in the case of HAT-P-7b, a UHJ. We describeour hierarchical modeling approach in Sect. 2 and present the results on the quenching levels for the atmosphere of HAT-P-7b in Sect. 3, where we add a critical evaluation of the quenching level concept. In Sect. 4 we highlight our findings on how much the abundances of H₂O, CO, and CH₄ change due to diffusion (both turbulent eddies and molecular) and photochemistry. We summarize and conclude our findings in Sects. 5 and 6.

2 Approach

We used our modeling tools as virtual laboratories, expanding our hierarchical approach of modeling environments from Paper I. In Paper I, we used 97 1D profiles extracted from the 3D general circulation model (GCM) calculated for the cloud-free HAT-P-7b (see Fig. 1 in Mansfield et al. 2018). Latitudes and longitudes are distributed with 22.5° and 15° intervals, respectively, which results in 96 1D profiles in one hemisphere. The zero-longitude reference is defined at the substellar point. An additional node at longitude = 90° and latitude = 85.5° was also chosen to probe the polar region. We then applied our kinetic cloud formation model to solve for the gas-phase composition under chemical equilibrium while consistently accounting for element abundance changes due to cloud formation. The next modeling step was to use the 97 1D (T_gas(z), p_gas (z), w(x, y, z))-profiles for HAT-P-7b (Fig. 1 in Paper I: T_gas(z) is the local gas temperature [K], p_gas(z) is the local gas pressure [bar], and w(x,y,z) is the local vertical velocity component [cm s⁻¹]; all are direct outputs of the GCM) to derive the diffusion coefficient Kzz, which assumes that all gas species have the same eddy diffusion coefficient.

In this section, we outline the two approaches we adopted to study the effect of gas-phase kinetic processes on the atmosphere of HAT-P-7b. We follow a hierarchical approach, wherein we begin in Sect. 2.1 with a simple quenching-level calculation (i.e., zeroth order), before we increase the complexity in Sect. 2.2 by considering full gas-phase chemical kinetic modeling.

2.1 Quenching levels: a zeroth-order approximation

Modeling kinetic gas chemistry is computationally expensive. It can therefore be instructive to explore fast but reasonably accurate approaches to investigate chemical quenching in 3D atmospheric models. As a zeroth-order approximation, we adapted the calculated quenching pressures by Venot et al. (2018a; presented in their Fig. 1) to determine the quenching levels for each 1D trajectory in the 3D model of HAT-P-7b. Briefly, Venot et al. (2018a) used a series of different thermal profiles (following analytical models from Parmentier & Guillot 2014) to determine the quenching level as a function of temperature and eddy diffusion coefficient (Kzz) by examining the criteria for at least one of their main species (i.e., H₂, H, H₂O, CH₄, CO, CO₂, N₂, NH₃, HCN, CH₃, and OH) to deviate from its thermochemical equilibrium abundance.

Because these functions are monotonic (the quenching pressure as a function of temperature and K_zz), their inverse functions can be easily constructed, that is, eddy diffusion coefficients at the quenching points (K) as a function of pressure and temperature. We linearly interpolated and extrapolated these K values (using radial basis functions) in the range of 1000 K < T_gas < 3000 K and 1 mbar < P_gas < 10⁶ mbar; see Fig. 1. The comparison of K with the derived eddy diffusion coefficient at any given pressure and temperature from the GCM 1D trajectories provides an estimate for whether that atmospheric level is quenched (i.e., if K_zz (P, T) > K, it is quenched).

In 1D kinetic models, the eddy diffusion coefficient (K_zz) is usually estimated by multiplying the vertical wind speed (w) by the relevant vertical length scale, for instance, the scale height (H) (see, e.g., Lewis et al. 2010; Moses et al. 2011). However, Parmentier et al. (2013) derived K_zz values from passive tracers in their GCM and found that these values were about two orders of magnitude lower than the traditional approach. Therefore, we approximated the mixing strength as follows: (1)

where H, the atmospheric scale height, is given by H = RT∕μg at any pressure level (R is the universal gas constant, T_gas is the gas temperature, μ is the mean molar mass, which is defined as the mass of the atmosphere at each level divided by the amount of substance in that level, and g is the gravitational acceleration). The average eddy diffusion coefficient at pressures higher than 1 bar shows a somewhat linear relation with pressure as log(K_zz) = 5.5–1.7log(P_gas); Fig. 2. Above this region, the eddy diffusion coefficient is less variable on average; it revolves around values of 2 × 10⁷ and 2 × 10⁶ cm²s⁻¹ on the day- and nightside, respectively. The stronger K_zz on the dayside, along with the cooler atmosphere on the nightside, indicate that between the day- and nightside lies a region with a minimum quenching pressure. This is discussed in the results section (Sect. 3.1). By taking this approach, small-scale and local variations ofK_zz might be introduced in comparison to the case of a smoothed parameterized K_zz approach. This means that small-scale variations of abundances might not be very realistic, and the results should be taken for their order-of-magnitude significance.

Although this approach is quite simple, it provides a reasonable overall picture about where we expect chemical quenching to occur. Nevertheless, a more sophisticated approach should be considered in order to quantitatively assess the importance of disequilibrium processes on UHJs. For this, we turn now to chemical kinetic modeling.

Fig. 1 Estimation of the eddy coefficient needed for quenching to occur as a function of pressure and temperature, K(p,T). Venot et al. (2018a) reported the quenching pressures for several pairs of temperature-vertical mixing strength (in their Fig. 1), and we interpolate their results to obtain a continuous K function at any chosen pressure and temperature. See Sect. 2.1 for more details.

2.2 Chemical kinetic modeling

A chemical network is essential to quantitatively trace the chemical abundances resulting from interactions between major and minor species at different pressures and temperatures. For this, we used ChemKM, a photochemical kinetic model (Molaverdikhani et al. 2019b), to perform a detailed study of gas-phase chemical kinetics and quenching levels on HAT-P-7b for CNOH-binding species. ChemKM is a modular package that calculates the vertical distribution and temporal evolution of atmospheric constituents in 1D. It employs the ODEPACK library (Hindmarsh 1983) and solves the coupled mass-continuity equations as a function of pressure for each species, (2)

where n_i is the number density (cm⁻³), Φ_i is the vertical flux (cm⁻² s⁻¹), P_i is the production rate (cm⁻³ s⁻¹), t is time, z is the vertical distance from the center of planet, and L_i is the chemical loss rate (s⁻¹) of species i.

Photochemical processes by UV and X-rays, photodissociation by galactic cosmic rays, flux attenuation due to Rayleigh scattering, condensation, ablation, outgassing, molecular diffusion, and eddy mixing could all drive the chemical composition of the atmosphere out of thermochemical equilibrium. However, here we only explore the effects of diffusion and mixing and photochemical processes by UV and X-ray because the other processes are poorly known for exoplanets.

In order to solve Eq. (2), we adopted a reduced version of the Venot et al. (2012) network that accounts for most of the dominant species expected at the atmospheric conditions of HAT-P-7b. While 97 1D profiles is not a huge number, developing a reduced network optimized for a particular case would reduce the computational time needed for studying the case. We developed the reduced network through a sensitivity analysis of the abundances of major opacity species (mainly water; see Sect. 4.5). By excluding the species in Venot et al. (2012) one by one and calculating the differences between the abundances, we obtained a list of 38 species (and their corresponding reactions). The current reduced network results in differences smaller than 10% at all pressures for major species compared with Venot et al. (2012). This is much smaller than other uncertainties in the kinetic chemical simulations and justifies order-of-magnitude studies by using this network, which is the aim of this work. Although this approach does not ensure the generality of this reduced network for any other environment, it can be used for cases that are similar to HAT-P-7b in terms of temperature range as well as metallicity and C/O ratio (both are assumed to have solar values). Comparison of our network with the recently published reduced network by Venot et al. (2019; with 30 species) shows that these networks share 23 species and only differ in some of the heavier molecules. A thorough comparison of these networks demands a further analysis.

Our reduced network contains 38 species (H₂, He, H, CO, H₂O, CH₄, N₂, NH₃, CO₂, HCN, C₂ H₂, O₂, O(3P), O(1D), C, N(4S), N(2D), OH, CH₃, ³CH₂, ¹CH₂, CH, NH₂, NH, H₂ CO, HCO, CN, C₂ H, NO, HNO, NCO, HNCO, HOCN, HCNO, C₂ H₃, C₂ H₄, C₂ H₅, and C₂ H₆), corresponding to 611 kinetic reactions. The UV absorption cross-sections and branching yields were adapted from an updated version of the Hébrard et al. (2012) network, which includes 29 photoreactions for the species in our reduced network. The photo-species are H₂O, CO₂, H₂ CO, OH, CO, H₂, CH₄, CH₃, C₂ H₂, C₂ H₃, C₂ H₄, C₂ H₆, N₂, HCN, NH₃, NO, and HCO.

The irradiation field was determined by interpolating between Phoenix stellar atmosphere models (Hauschildt et al. 1999; Husser et al. 2013)¹, assuming the following stellar parameters: T_* = 6441 K, log(g) = 4.02 cm s⁻², and [Fe/H] = 0.15 (Torres et al. 2012). While some UV observations by Hubble Space Telescope’s (HST) the Space Telescope Imaging Spectrograph (STIS) and the Cosmic Origins Spectrograph (COS) instruments are available for similar stars such as HD 27808 (e.g., Bowyer et al. 1994), we used the Phoenix models to take a coherent approach. This means that no data stitching is required in order to construct the host star flux between 1 nm to 300 nm (the wavelengths at which the irradiation becomes relevant to the photolysis of atmospheric constituents). In addition, the dayside of the planet is very hot and therefore likely to be insensitive to the photochemistry (see Sect. 3.3 for details). Therefore, small differences in the UV portion of the constructed stellar spectra using different methods (i.e., data stitching or Phoenix model) are not expected to cause a noticeable change in the results, and the Phoenix model provides an adequately accurate estimate of the HAT-P-7 spectrum. It is also worth noting that small differences in the stellar spectra are inconsequential due to relatively large uncertainties in the data used to model photodissociation. These data, such as cross-sections and branching ratios, are poorly known at temperatures higher than 300 K (Venot et al. 2018b).

Following our hierarchical modeling approach, we used the gas-phase chemical equilibrium results from Paper I as input values. We then computed the thermochemical abundances of species considered in our kinetic network to calculate the initial abundances for the kinetic gas-phase modelling in this paper, as discussed in Sect. 2.2. The results from Paper I include the effect of cloud formation on the element abundances through depletion and/or enrichment (see Sect 2.1. in Paper I). We used the same 1D atmospheric temperature-pressure profiles (Figs. 1 and 3 in Paper I). We present and discuss the results of our kinetic gas-chemical models in Sect. 3.2.

In order to solve the continuity equation (Eq. (2)), we followed Hu et al. (2012; Eq. (2)) and express the vertical transport flux of species i (Φ_i) as (3)

where K_zz is approximated using Eq. (1), n is the total number density, H₀ is the mean scale height of the atmosphere, H is the molecular scale height, T_gas is the gas temperature, f_i is the mixing ratio, D_i is the molecular diffusion coefficient, and α_i is the thermal diffusion factor of species i.

The thermal diffusion factor is usually set to zero, except for light species such as H, H₂, and He, which can be estimated with −0.25 (e.g., Marshall & Plumb 2013). To estimate molecular diffusion coefficients, D_i, we used the Chapman and Enskog equation (Poling et al. 2000), (4)

where M_AB is defined as (M_A and M_B being the molecular weights of species A and B, respectively), Ω_D is the collision integral for diffusion, σ_AB is the characteristic length of the intermolecular force law, f_D is a correction term (usually on the order of unity), n is the number density of molecules in the mixture, k is the Boltzmann constant, and T is the gas temperature of the mixture. A detailed description of how we approximate the molecular diffusion coefficients is provided in Appendix B.

As pointed out, the gas-phase number density profiles of the 38 CHNO-bearing species, as resulting from the kinetic cloud formation presented in Paper I, were used as input values for our chemical kinetic gas-phase calculations. These input molecular abundances already account for elemental depletion or enrichment by cloud formation processes. However, we note here that the comparison with our gas-chemistry results from Paper I bears a conceptional problem: the definition of an interface.

Fig. 2 Calculated eddy diffusion coefficient (Kzz) profiles using Eq. (1). The average eddy diffusion coefficient at pressure higher than 1 bar roughly follows a simple relation as log10(Kzz) = 5.5–1.7log10(pgas). Above this pressure and on the dayside, Kzz has average values about 2 × 107cm2s−1 (vertical reddashed line), and it is weaker by about one order of magnitude on the nightside (vertical blue dashed line).

Defining the interface

A challenge for hierarchical models may be to define an interface that is unambiguous and practical at the same time. The comparison between results from an equilibrium gas-phase code with a kinetic rate network gas-phase code is challenging because the kinetic model only treats a limited number of gaseous species. In this paper, the interface is therefore defined by the number densities, n_x∕n_tot [cm⁻³], of those species that can be considered by the kinetic code (H₂, He, H, CO, H₂O, CH₄, N₂, NH₃, CO₂, HCN, C₂ H₂, O₂, O(3P), O(1D), C, N(4S), N(2D), OH, CH₃, ³CH₂, ¹CH₂, CH, NH₂, NH, H₂ CO, HCO, CN, C₂ H, NO, HNO, NCO, HNCO, HOCN, HCNO, C₂ H₃, C₂ H₄, C₂ H₅, and C₂ H₆). Hence, the input (interface) values for all quenching approaches (Sects. 2.1–2.2) considered here are set by the results from our thermochemical equilibrium calculation in Paper I. The initial abundances in our kinetic model are then computed by calculating the new thermochemical equilibrium abundances of these 38 species, given these input values. Any significant deviation from these initial abundances through time-dependent kinetic calculations is considered as an indication of the fingerprints of disequilibrium processes, as discussed in Sect. 3.2.

This choice of interface, however, results in lower element abundances of C, N, and O if they are derived based on the 38 species alone (rather than from the full list of species used in Paper I). If the element abundances were kept consistent between the two modeling steps, we would have inconsistencies in the initial molecular number densities compared to the values derived in Paper I. See Sect. 4.4 for more details. We also note that our 1D approach neglects horizontal transport, which could alter the exact location of the quenching and abundances (e.g., Agúndez et al. 2014; Drummond et al. 2018a; Mendonça et al. 2018). This limitation is discussed in Sect. 5.

3 Quenching levels

3.1 Results of zeroth-order approximation

The quenching levels by the zeroth-order approximation were determined by following the criteria of at least one of the main species in the Venot et al. (2018a) model (i.e., H₂, H, H₂O, CH₄, CO, CO₂, N₂, NH₃, HCN, CH₃, and OH) that deviates from its thermochemical equilibrium abundance. Therefore, this approximation does not associate the quenching levels with any specific species, unlike the kinetic approach.

The top panels of Fig. 3 reproduce the transection maps of K_zz and the quenched regions at the equator. Despite weaker vertical mixing at the nightside, quenching levels are deeper (10–100 mbar) than at the dayside (1–10 mbar) because the atmospheric environment is hotter by about 2000 K on the dayside. This is more evident in the noon-midnight maps, that is, Fig. 3 middle-panels. The noon-midnight transection map shows a strong increase in the quenching pressure from the dayside (~1 mbar) to the nightside (~100 mbar), with no significant latitudinal quench point variation on either side and a sudden change in the quenching levels over the pole. The transition, however, is more gradual between the evening and morning terminators, see Fig. 3 bottom-panels.

The value of the quenching pressure depends on the mixing efficiency (K_zz) and temperature. Komacek et al. (2019) found that in general, K_zz increases with temperature in the mbar regime. Because both of these quantities decrease at higher zenith angles (i.e., toward the nightside), there should be a region where the quenching pressure is at its minimum values (i.e., the quenching point occurs at its maximum altitude) for the affected species by diffusion. This is evident in Fig. 4 where the quenching levels have a maximum altitude at the evening (hotter) terminator according to our zeroth-order approximation. This shows that the thermochemical equilibrium assumption is the most relevant at the evening terminator. This may partially support the applicability of this assumption for the transmission spectroscopy of UHJs. Further analysis will be performed in a forthcoming paper to investigate this in more detail.

As shown in both Figs. 3 and 4, deeper regions usually remain in chemical equilibrium, but there might be some detached regions at the lower pressures that do not satisfy the Venot et al. (2018a) quenching criteria. For instance, the blue ring in the terminator transection map, Fig. 3 bottom right panel, where the chemistry departs from the quenched profiles and returns to the equilibrium profiles at some pressures. This is mostly due to a combinationof temperature inversion and lower vertical velocity. Nevertheless, if only one quenching level at each 1D trajectory is desired, then the deepest level at which the quenching occurs can be selected as the quenching point.

3.2 Results of chemical kinetic modeling

Temporal variations of the H₂O, CO₂, and CH₄ vertical abundances at the equator and longitude ϕ = 225° in Figs. 5a, d, and g to demonstrate how disequilibrium processes change the composition of the upper atmosphere (see Sect. 2.2 for a description of the chemical kinetic model). The temporal evolution of the molecular abundances of H₂O, CO₂, and CH₄ are color-coded by time, where blue to red shows the progress in time from thermochemical equilibrium to their disequilibrium steady states. All three molecular abundances depart from thermochemical equilibrium values at lower pressures and at some point remain constant up to the top of the atmosphere. CH₄ is quenched at around 0.1 bar for this particular trajectory (Fig. 5g).

The quenching level of H₂O, however, cannot be estimated easily and at first glance might lie somewhere between 0.3 and 30 bar. Its thermochemical equilibrium abundance is almost constant throughout this region, and diffusion therefore does not significantly change its abundance below 0.3 bar. Figure 5d shows yet another complication with the quenching-level determination. In this case, CO₂ departs from its thermochemical equilibrium value at around 0.3 bar, while it does not reach a constant value until around 3 mbar. This usually occurs when the chemical and mixing timescales are comparable in that region. Examination abundance variation of all species at all trajectories reveals that in many cases the abundances do not reach a constant value at lower pressures and accordingly, no quenching level (and consequently, no quenched abundance) can be defined. Now, the question is what shouldbe reported as the quenching level.

We follow Molaverdikhani et al. (2019b) in suggesting the use of the coefficient of variation (CV) to estimate the deviation from chemical equilibrium. This quantity can be calculated as the ratio of the standard deviation of a given species abundance, s, to its temporal mean abundance at any given pressure level. CV is a relative and dimensionless measure of dispersion and hence independent of the absolute value of the mixing ratios. For each 1D trajectory we can therefore calculate a 1D CV profile. Because in our kinetic model the time step for the numerical integration is exponential, we use the geometric coefficient of variation (gCV) of abundances for this purpose, calculated as (Schiff et al. 2014) (5)

where s_ln is the scaled natural log of the standard deviation of abundances and can be estimated as s_ln = sln(10). Figures 5b, e, and h shows the calculated gCV profiles of H₂O, CO₂, and CH₄, respectively, for a latitude θ = 0°/ longitude ϕ = 225° trajectory as an example. A higher gCV value represents a stronger abundance variation at any given pressure level. A choice of gCV = 0.05 seems to reasonably separate the chemical equilibrium and disequilibrium regions. The geometric mean of gCV for all species, , can be also used to determine the quenching level of the atmosphere, (6)

The gCV can be thought of as a measure of normalized abundance deviation from its mean value and within a given time period. The chosen time interval during which the abundances are considered should therefore start at a time at which the atmosphere is at the thermochemical equilibrium and end after a new steady state is reached, such as diffusion equilibrium.

Transection gCV-maps of H₂O, CH₄, and for all species are shown in Fig. 6. Interestingly, the zeroth-order approximation results (red solid lines in Fig. 6 right panels) agree well with our kinetic modeling results, although they generally underestimate the depth of quenching on the dayside. An improved zeroth-order approximation can be foreseen as a fast method to estimate the quenching levels in parameterized atmospheric simulations.

The gCV values can be also used to compare the effectiveness of disequilibrium processes on different atmospheric constituents. For instance, H₂O has been driven away from thermochemical equilibrium less effectively than CH₄ given gCVH₂O < gCVCH₄; Fig. 6, left and middle panels. values also indicate that most of the species have been driven out of thermochemical equilibrium more effectively than these two species on the nightside at the equatorial plane. CH₄ equatorial map shows an asymmetry of gCVCH₄ with a maximum value around longitude of 225°. This is due to lower temperatures on the nighside (Figs. 4 and A.7) and higher vertical mixing (Fig. 3) at this longitude at the equator.

The noon-midnight gCV-maps (Fig. 6 middle panels) support shallower quenching levels on the dayside and a slight latitudinal dependency on the nightside. The anti-substellar point (Lat = 0°, Lon = 180°) appears to behave slightly differently relative to its neighboring trajectories, mostly because of the stronger vertical velocity (w) in this region.

The deepest quenching levels occur at the polar regions, where the temperature is lower and the mixing is relatively strong. The map suggests a morning-evening asymmetry where the morning-side abundances have been driven more significantly out of thermochemical equilibrium and at deeper pressure levels. This asymmetry can potentially shape the transmission spectra of the two limbs differently. Observational consequences of chemical disequilibrium will be discussed in a forthcoming paper.

Fig. 3 Left: calculated Kzz using Eq. (1). Right: calculated quenching levels by following a zeroth-order approach. In the blue regions the atmosphere is in thermochemical equilibrium, and orange highlights the regions that are expected to be affected by disequilibrium processes and at least one major species deviates from its thermochemical equilibrium abundance.

Fig. 4 Variation of quenching pressure with the zenith angle, using the zeroth-order approximation (Sect. 3.1). Temperature profiles of all trajectories are color-coded above the quenching pressures. Minimum quenching pressure occurs at theevening terminator, where the atmosphere contributes significantly in the transmission spectra. The effects of disequilibrium chemistry are therefore expected to be less pronounced in the transmission spectra. Variation of quenching pressure is represented with the dashed line.

3.3 Effect of photochemistry on the quenching levels

Photochemistry is another disequilibrium process. The main source of photolysis on highly irradiated planets is irradiation from the host star. Other sources (such as cosmic rays, GCR, stellar background UV radiation, or solar Lyman-α photons that are scattered from atomic hydrogen in the local interplanetary medium, LIPM) could be also important on weakly irradiated planets, and can be neglected in the case of HAT-P-7b. The photochemistry effect on the quenching levels on a highly irradiated planet remains to be defined.

Figure 7 compares the latitudinal quenching levels, , with and without photochemistry. The bottom panel shows the differences and suggests that the photochemistry plays a negligible role. This is expected because the photochemistry mostly affects the mbar regime and higher altitudes (see, e.g., Molaverdikhani et al. 2019b and references therein) and the HAT-P-7b dayside quenching points are located at deeper levels. The two exceptions are the evening terminator and a portion of the morning terminator, where the quenching pressures can be lower than mbar, as discussed in Sect. 3.1. However, at these regions, the zenith angle is 90°; irradiation was therefore ignored in our 1D kinetic calculations. The nightside does not change either because GCR and LIPM are not included.

Fig. 5 a, d, g: temporal variation of the H2O, CO2, and CH4 vertical abundances at the equator and longitude = 225°, i.e., nightside. b, e, h: gCV of abundances, from initial values to a steady state at ~ 5 × 109 s. A higher value represents a stronger variation of the abundance at any given pressure level. A choice of gCV = 0.05 seems to reasonably separate the chemical equilibrium and disequilibrium regions (dashed gray vertical lines). c, f, i: quenching timescales. An invariant abundance profile at the longer times ensures a new steady state. See Sect. 3.2 for more details. We note that the quenching level cannot be associated with one particular pressure for all species, but usually spans more than one order of magnitude in pressure.

3.4 Effect of mixing strength on the quenching levels

As discussed in Sect. 2.1, we decreased the K_zz values by two orders of magnitude, Eq. (1), to follow the conclusions of Parmentier et al. (2013). The traditional mixing strength approximation, that is, without decreasing K_zz by a factor of 0.01, would then result in different quenching levels. Figure 8 compares the two cases, and as expected, the difference is non-negligible.

The resulted quenching levels from the traditional approach, K_zz = w ⋅ H, are usually deeper than the modified one, K_zz = w ⋅ H × 10⁻². The differences are on the order of 20–40% on the dayside and 40–80% on the nightside. The quenching levels at the terminator trajectories, however, show stronger variations as high as several orders of magnitude.

4 Disequilibrium abundances due to diffusion and photochemistry

4.1 CO and CH₄ abundance in the equatorial and terminator regions

As noted, the initial conditions were set to the equilibrium results from Paper I, where the element depletion due to cloud formation affects the gas-phase abundances in the cloud-forming regions of the atmosphere. These initial abundances, the steady state solution of the kinetic model and their differences are shown in Figs. 9 and A.1 for CO and CH₄.

The carbon monoxide abundance (Fig. 9) on HAT-P-7b is only marginally affected by quenching, and the kinetic modeling results vary by less than a factor of two compared to the number densities in chemical equilibrium for both equatorial and terminator regions. On the other hand, the methane abundance (Fig. A.1) in the kinetic model shows a stronger deviation from the equilibrium results at the ~10 mbar level. The strong vertical mixing causes the number density of methane at P < 1 bar around longitude ϕ = 240° to decrease by a factor of about 1000. An enhancement in methane number density at p_gas > 1 bar by more than ten times is also noticeable. At the morning terminator methane is strongly enhanced between the 1–10 mbar level. We interpret this to mean that the methane enhancement is related to the thermal inversion at the ~10 mbar level (see Figs. 4 and A.7). The vertical mixing homogenizes the relatively methane-rich layer of the cooler region at p_gas > 10 mbar with the hotter region at p_gas < 10 mbar. As a result, vertical mixing leads to a local enhancement in methane number density between the 1–10 mbar level relative to that of the equilibrium result.

4.2 H₂O abundance in the equatorial and terminator regions

We show the effects on the H₂O abundance (Fig. A.2) in the same way because with carbon monoxide, the abundance is relatively unchanged and continues to display an asymmetry between the dayside and nightside. However, some subtle variations are notable. While abundances at pressures higher than 1 bar mostly ramin constant at all longitudes around the equator, abundances at higher altitudes on the nightside increase slightly. On the dayside, the abundance difference between equilibrium and kinetic chemistry is almost invariant at pressures higher than 100 mbar, and at higher altitudes, it decreases by a factor of about 30 on average.

The terminator also shows only small variations, but a prominent change is the shifting of relatively low-abundance regions from around 0.01 bar between ϕ ± 67.5 degrees latitude at the morning side to two vertical columns at the same latitudes on the evening side. The most significant change of a factor of nearly three orders of magnitude decrease occurs just on the dayside of the morning terminator (θ ~ 285 degrees longitude) at 1–10 mbar. This is likely linked to the thermal inversion there, as was the case for methane.

The abundance variation results due to disequilibrium are consistent with the interpretations based on gCV calculations, as discussed in Sect. 3.2 for the two cases of H₂O and CH₄. This means that gCV can be used to determine quenching pressures and as a quantitative measure of the efficiency of disequilibrium processes on the atmospheric constituents.

Fig. 6 Estimation of quenching levels using the chemical kinetic modeling (Sect. 3.2). Top panels: transection maps of the geometric gCV of H2O, CH4, and all species abundances at the equator (latitude = 0°). We choose gCV = 0.05 to separate the chemical equilibrium and disequilibrium regions (dashed white lines). This boundary agrees well with the zeroth-order approximation results (red solid lines in the right panels; Fig. A.7). Middle panels: similar to the top panels, but at the noon-midnight plane (longitude = 0° or 180°). Bottom panels: similar to the top panels, but at the terminator (longitude = 90° and 270°). map suggests a morning-evening asymmetry where the dawn-side abundances have been driven more significantly out of thermochemical equilibrium and at deeper pressure levels.

4.3 CH, C₂ H₂, CN, and HCN abundance in the equatorial and terminator regions

Here we show that the abundances of small carbon molecules CH (methylidyne, Fig. A.3) and CN (cyanide, Fig. A.4) may be affected by vertical transport processes compared to the chemical equilibrium results in Paper I. We note, however, that the atmosphere is globally oxygen rich (bulk carbon-to-oxygen ratio (C/O) = 0.54), therefore oxygen-binding molecules will dominate the chemical composition of the atmosphere. It is worthwhile to consider carbon molecules that can be of relevance for the formation of biomolecules (such as CH₄, C₂ H₂ and HCN) and their precursors, however, in order to identify possible tendencies in the unexplored parameter ranges that exoplanets provide. HCN is furthermore suggested as spectral tracer for lightning events in planetary atmospheres. Knowing itsbackground abundance may help to disentangle such lightning traces (Hodosán et al. 2017; Helling & Rimmer 2019).

CN and CH are precursors for more complex molecules as well, such as C₂ H₂ (Fig. A.5) and HCN (Fig. A.6). CH, CN, C₂ H₂, and HCN are molecules that typically appear in spectra of carbon-rich stars (Eriksson et al. 1984; Helling et al. 1996) and maybe therefore be considered as potential tracers for carbon-enriched planetary environments (e.g., Hron et al. 1998; Kopparapu et al. 2011; Madhusudhan 2012; Moses et al. 2013; Venot et al. 2015; Uttenthaler et al. 2016; Hamren et al. 2016).

First, CH and CN show a clear day-night asymmetry: they are more abundant on the dayside and in the warmer atmospheric regions. This day-night asymmetry is unaffected by quenching because it is driven by the local day-night temperature difference instead of hydrodynamic motions. The CH abundance has decreased on the nightside compared to the thermochemical equilibrium results. Differences are generally small but amount to orders of magnitudes in localized areas in the equatorial plane of as much as a factor of 10¹² (i.e., 12 dex). Differences remain negligible in the low-pressure terminator regions (p_gas < 1 bar).

C₂ H₂ is very lowly abundant, and kinetic chemistry only affects its abundances in localized areas in the terminator regions at p_gas ≈ 10⁻²bar. HCN is somewhat more abundant than C₂ H₂, but disequilibrium effects introduce a day-night asymmetry for HCN. As a result, the nightside HCN abundance has increased by three orders of magnitude compared to the chemical equilibrium results. The absolute number densities, however, remain small compared to CO and H₂O, for example. This asymmetry is also visible in the terminator regions, where the morning terminators show a higher HCN abundance. The terminator regions (right column) are again only locally affected by the kinetic gas chemistry, mainly at the nightside at p_gas ≈ 10⁻² bar. In general, we note that all molecules in the form of C_n H_m are affected by disequilibrium processes on the nightside, where heavier C_n H_m molecules are affected the most. The result in this section emphasizes our finding from the previous sections: quenching of the gas-phase chemistry does take place for some of the molecules, but the actual deviation needs to be considered with care because large deviations might be only local phenomena.

Fig. 7 Top: latitude and longitude dependency of quenching pressures, calculated from the kinetic simulations including only diffusion. The quenching criterion is discussed in Sect. 3.2 Middle: similar to the top panel, but photochemistry is also considered. Bottom: difference between the quenching pressures of models with diffusion alone and models including both diffusion and photochemistry.

Fig. 8 Top: latitude and longitude dependency of quenching pressures, calculated using Kzz = w ⋅ H × 10−2. Middle: similar to the top panel but using Kzz = w ⋅ H, i.e. traditional approach. Bottom: comparing the quenching pressures of models with Kzz = w⋅ H to the quenching pressures of models with Kzz = w ⋅ H × 10−2. The atmosphere is being quenched at deeper levels, usually 20–40% on the dayside and 40–80% on the nightside. Some terminator trajectories, however, show several orders of magnitude change in their quenching levels, deepening from the mbar to 1 bar pressures.

Fig. 9 Carbon monoxide (CO): left column: equatorial slices (θ = 0°), right column: terminator slices. Top: initial condition, . Middle: finalresult from the kinetic simulation, . Bottom: difference, between the initial (equilibrium) values and the kinetic results.

4.4 Atmospheric C/O in HAT-P-7b

The carbon-to-oxygen ratio (C/O) from ≪1 to ≫1 has served as one of the stellar parameters that allowed the link to stellar evolution, and the spectral changes have been studied widely (Gustafsson et al. 2008; Van Eck et al. 2017). Its effects on the spectra of a planetary atmosphere have also been studied extensively (e.g., Seager et al. 2005; Madhusudhan 2012), and it has been shown to change the atmospheric properties significantly (e.g., Mollière et al. 2015; Molaverdikhani et al. 2019a,b). The C/O ratio thus became one of the parameters that was hoped to be useful in characterizing the formation and/or evolution of an exoplanet.

However, C/O cannot be straightforwardly linked to planetary evolutionary states or easily to planet formation when clouds form (Helling et al. 2014). We note that other mineral ratios (Mg/O, Mg/Si, etc.) will be similarly unsuitable as a direct link to planet formation unless cloud formation and related chemical processes are simulated simultaneously (see, e.g., Fig. 4 in Helling 2019).

The main contribution to changes of C/O in our simulations comes from the change in oxygen abundance due to consumption in cloud particles. Our cloud formation simulations for HAT-P-7b start from solar element abundances with C/O = 0.54, that is, a higher oxygen- than carbon-element abundance. The resulting C/O for the cloud-forming atmosphere of HAT-P-7b linked to the equilibrium gas-phase code GGChem (Woitke et al. 2018)² is reproduced in Fig. 10 (top left) for the equatorial slice (θ = 0°): the cloud-forming nightside reaches C/O ≈ 0.78, and the cloud-free parts of the dayside remain at the undepleted solar value, as expected. The oxygen- and carbon-element abundances that result from our kinetic gas simulations are shown in Fig. 10 (bottom left).

We compare in how far the undepleted solar C/O is reflected in the initial conditions for the gas kinetic analysis conducted in this paper. Then, we study how C/O may have been affected by chemical gas kinetic processes. We present our results in the form of slice plotsthrough the equatorial plane in Fig. 10. The corresponding terminator slices can be found in the appendix (Fig. A.8).

The C/O we used as initial condition for the kinetic modeling (top right in Fig. 10) is lower than the C/O derived in Paper I (top left in Fig. 10). The reason is that such rate networks consider fewer molecules than our gas-phase equilibrium calculations with GGChem, for example, and C/O is therefore derived from the C- and O-binding molecules alone considered in the gas kinetic network simulation (see list in Sect. 2.2). The input abundances for these molecules were set to those resulting from our cloud formation calculations in Paper I, as discussed in Sect. 2.2.

The gas kinetic C/Os appear to be relatively homogenized toward a value slightly higher than its solar value. This is expected because the global initial C/O in the gas phase was higher than the solar value due to the depletion of oxygen by cloud formation. The largest differences between (C/O)_eq and (C/O)_kin appear for p_gas < 1 bar on the nightside and to some extent on the morning terminator. All other regions show negligible deviations between (C/O)_eq and (C/O)_kin in the equatorial plane. This homogenization of C/O removes the morning-evening C/O asymmetry due to the cloud formation alone (reported in Paper I). This suggests that the C/O asymmetry may be less pronounced in the UHJs, but an observed global enhancement of C/O could still be an indication of cloud formation on a hot exoplanet.

4.5 Observability of disequilibrium chemistry in the atmosphere of HAT-P-7b

After examining how molecular abundances are modified by accounting for chemical kinetics, we now turn to the prospect of observing the effect of disequilibrium chemistry on the atmospheric spectra. There are many ways in which disequilibrium chemistry can manifest itself in an exoplanet spectrum. The clearest diagnostic would be the appearance of absorption features caused by a molecule that is expected to be undetectable for equilibrium chemical abundances. For example, signatures of NH₃ and HCN can become detectable in hot-Jupiter transmission spectra if disequilibrium mechanisms enhance their abundances to ≳ 1% that of H₂O (MacDonald & Madhusudhan 2017). However, such a dramatic change requires molecular abundances to alter at the order-of-magnitude level. Alternatively, for smaller deviations from chemical equilibrium, we may examine whether the shape of molecule absorption features, or the relative strengths of features due to different molecules, are altered by chemical kinetics.

We quantitatively assesed the effect of disequilibrium chemistry on spectral observations of HAT-P-7b using a similar approach to Paper I. Specifically, we computed the pressure at which the vertically integrated optical depth due to a chemical species, x, reaches unity (i.e., p(τ_x(λ) = 1)). This provides a useful indication of the contributions of different molecules to spectra of HAT-P-7b. We first consider which of the molecules included in our chemical model are likely to be spectrally accessible to observations. We then examine differences due to chemical kinetics, quantified by the percent difference in p(τ_x(λ) = 1) between the initial and final atmospheric abundances.

Figure 11 shows the p(τ_x(λ) = 1) surfaces for prominent molecules at the end of the kinetic calculation. Here we focus on infrared wavelengths because the opacity at visible wavelengths will contain large contributions from atomic and heavy-element molecules that are not included in the present kinetic chemistry model (see Paper I, Fig. 20, for their effect on chemical equilibrium). Three molecules are potentially observable at near-infrared wavelengths: H₂O, CO, and CO₂ (ordered by prominence). The importance of these molecules holds in all four of the plotted atmosphericregions, but their relative influence does vary (e.g., the opacity of CO and CO₂ is more prominent than that of H₂O at the morning terminator than at the evening terminator). At the substellar point, OH also becomes an important opacity source for wavelengths longer than >10 μm (opening the possibility of its detection through mid-infrared dayside emission spectroscopy, such as with the MIRI MRS mode of the James Webb Space Telescope; see Beichman et al. 2014). While H⁻ is reported to be an important opacity source on the highly irradiated exoplanets (see, e.g., Arcangeli et al. 2018), we found that its contribution(at thermochemical equilibrium) is likely to be significant only on the dayside of HAT-P-7b and at wavelengths longer than 15 μm (Helling et al. 2019b). This suggests that the effect of this compound on the transmission spectra of this planet is not significant.

It is important to note that molecular opacity at the antistellar point and morning terminator are likely to be obscured by clouds (see Paper I). This overall picture is essentially the same as in Paper I, indicating that relaxing chemical equilibrium does not alter the prediction for which molecules are expected to be most prominent in spectra of HAT-P-7b.

However, as we showed, chemical kinetics plays an important role in shaping the vertical abundance profiles of various molecules. For example, Fig. 5 demonstrates nightside quenching of the H₂O abundance in the upper atmosphere for pressures ≲10⁻¹ bar. Because H₂O is the most prominent molecular opacity source at infrared wavelengths, any abundance changes at atmospheric pressures such as this that are accessible to transmission and emission spectra are potentially observable. We examine the importance of disequilibrium abundance changes by computing the percent difference between the initial (equilibrium) and final (kinetic) p(τ_x(λ) = 1) surfaces for the molecules that are anticipated to be most prominent in spectra of HAT-P-7b. The results are shown in Fig. 12.

Disequilibrium chemistry tends to manifest itself by enhancing the prominence of molecular absorption features. The differences are most pronounced at the antistellar point and morning terminator, where the p(τ_x(λ) = 1) surfaces occur up to 20% higher in the atmosphere for H₂O (and up to 60% higher for CO₂) in the center of absorption features (e.g., 3 μm) and up to 10% lower inabsorption minima (e.g., 1.7 μm). This disequilibrium-induced contrast in the pressures where the atmosphere becomes optically thick between wavelengths where molecules strongly andweakly absorb would result in stretched spectral absorption features. On the other hand, at the substellar point and evening terminator, the p(τ_x(λ) = 1) surface for H₂O increases in altitude for all wavelengths, which would result in stronger but relatively unstretched absorption features. We thus see that in principle, disequilibrium chemistry can transform the shape of molecular absorption features in different ways for different atmosphericregions. Although the antistellar point tends to be inaccessible to observations (because of its low thermal emission), these differing vertical opacity distributions between the morning terminator (observed in transmission) and the substellar point (observed in emission) could be one avenue to probe disequilibrium chemistry in UHJs. A second avenue is inter-terminator differences observed in transmission spectra that are caused by enhancements in H₂O, CO, and CO₂ absorption in the morning terminator, while the evening terminator sees only minor changes. Ultimately, whether or not these disequilibrium-induced opacity differences will be observable depends crucially on the presence of clouds in the atmosphere of HAT-P-7b. If the morning terminator is dominated by cloud opacity (as predicted in Paper I), it will be difficult to observationally probe disequilibrium chemistry in this region. Dayside emission spectra may detect H₂O features more easily because its abundance is higher, but based on this alone, it would be difficult to uniquely attribute this to disequilibrium chemistry (as opposed to a higher atmospheric O/H ratio but in chemical equilibrium, e.g.). We therefore conclude that it will be challenging to observationally detect signatures of disequilibrium chemistry in the atmosphere of HAT-P-7b. A quantitativeassessment of the detectability of the disequilibrium signature on HAT-P-7b will be presented in a forthcoming paper.

Fig. 10 Carbon-to-oxygen ratio (C/O): top left: C/O after cloud formation as in Paper I according to ɛC ∕ɛO, Top right: C/O for initial condition of gas-kinetic calculations (from CH, CO, CO2, C2 H2, etc. Bottom left: C/O at the end of gas-kinetic calculations (t = 1012 sec) Bottom right: C/O change between initial condition (top right) and final gas-kinetic solution (bottom left).

Fig. 11 Prominent sources of molecular opacity under disequilibrium chemistry. The atmospheric pressure at which the vertical wavelength-dependent optical depth of species x reaches unity, p(τx(λ) = 1), is shown across the infrared. Several molecules included in both the kinetic model and the opacity database that do not reach τx(λ) = 1 by the 10 barlevel (e.g., C2H2) are omitted. See Paper I (Fig. 20) for opacity contributions due to species not included in the kinetic model (heavy-metal oxides, atoms, etc.). All opacities are plotted at a spectral resolution of R = 1000. Four equatorial regions are considered: (i) the antistellar point, (ii) the morning terminator, (iii) the substellar point, and (iv) the evening terminator. For comparison, the τ(λ) = 1 surface due to clouds from Paper I is overlaid.

5 Discussion

Here we discuss the limitations of our work. The caveats that must be taken into account when the results are interpreted are as follow:

Material properties: rate and diffusion coefficients

Gas kinetic rate networks solve diffusion equations, one for each molecular species (e.g., Eq. (22) in Rimmer & Helling 2016). Each equation requires material constants in the form of a rate, an eddy diffusion coefficient, and a molecular diffusion coefficient. In principle, they can all be measured through laboratory experiments, but measuring them for many reactions or reactants is an overwhelming task. We assumed that the rate values are reasonably well known (or have been benchmarked) and did not investigate the propagation of uncertainties in the rate coefficients into our results. The eddy diffusion was assumed to be the same for all species. The molecular diffusion of species was calculated following the Lennard-Jones theory (see Appendix B for a detailed description). This is a critical point for larger molecules in low-pressure regions where molecular diffusion becomes more efficient than turbulent eddy diffusion.

Fig. 12 Differences to the opacity contribution of observable molecular species due to disequilibrium chemical kinetics. The difference metric, , is defined such that a positive difference implies stronger absorption after accounting for disequilibrium chemistry (i.e., p(τx (λ) = 1) occurs higher in the atmosphere). Only H2O, CO, and CO2 are plotted because these molecules offer the greatest detection prospect in the near-infrared (see Fig. 11). All opacities are plotted at a spectral resolution of R = 1000. Four equatorial regions are considered: (i) the antistellar point, (ii) the morning terminator, (iii) the substellar point, and (iv) the evening terminator.

Network size

Gas kinetic simulations are limited by the number of reactions that are included. The more focused a network is, that is, the more reduced a network is, the faster it will be. The ultimate form of such a reduced network manifests itself in the chemical timescale relaxation scheme, where only one rate limiting reaction is considered. The focus of such a simulation is usually on a chemical inversion between only two chemical compounds. This is the most frequently used technique to couple 3D GCM simulations with chemical disequilibrium processes.

As shown by Drummond et al. (2018a), for example, a consistent coupling of dynamics, radiative transfer, and chemistry demands more computational power and may result in uncovered physics by the relaxation scheme. For instance, an outcome of such full coupling leads to wind-driven chemistry that in turn has a significant effect on thermal and dynamical structure of cool atmospheres, such as that of HD 189773b. The effect of full network coupling, however, is reported to be less significant in the case of hotter atmospheres, such as that of HD 209458b. It remains to be decided whether this is also the case for the UHJs based on a full coupling of chemistry and GCM.

Limitations of 1D kinetics models

One important limitation of our approach is that our chemical kinetic model is 1D. Thus, horizontal mixing is neglected. Horizontal mixing can homogenize abundances longitudinally (e.g., Agúndez et al. 2014) and/or latitudinally and lead to horizontal quenching. Recent studies coupling a simplified chemical relaxation scheme to 3D general circulation models of hot Jupiters have found that horizontal mixing dominates vertical mixing (Mendonça et al. 2018). Coupling 3D GCM with reduced chemical networks also confirms this and showed that abundances are determined by a combination of vertical and horizontal mixing (Drummond et al. 2018a,b).

However, all of these studies focused on hot Jupiters with equilibrium temperatures < 1500 K, and it is notclear to what extent these results apply to UHJs. Komacek et al. (2017) demonstrated that with increasing equilibrium temperature, the day-to-night flow becomes more important than the equatorial jet, and vertical velocities increase. As a consequence, the efficiency of vertical mixing is believed to increase by more than an order of magnitude between T_eq = 1500 K and T_eq = 2200 K (Komacek et al. 2019). On the other hand, vertical velocities are set by the horizontal ones through the mass conservation. In principle, stronger vertical velocities should therefore result in stronger horizontal velocities. Future studies applying GCMs including simplified chemical schemes, such as a new reduced network recently proposed by Venot et al. (2019), to a wider range of planets will help clarify this aspect. In the long term, coupling full chemical networks to GCMs will hopefully become possible.

A related but different limitation is that vertical mixing in kinetic models is parameterized as an eddy diffusion coefficient. Zhang &Showman (2018a,b) demonstrated that this approximation is not always valid: the interaction of chemical species with the 3D circulation can produce situations in which they do not behave diffusively. Furthermore, the authors found that species with different chemical lifetimes can have different eddy diffusivities, an effect that commonly has not been taken into account in the kinetic models.

We therefore stress that even though we used a dedicated modeling approach, our results are meant to provide a first exploration of the extent to which disequilibrium chemistry plays a role on HAT-P-7b. Our results provide an estimate of the importance of disequilibrium abundances and in which regions of the atmosphere deviations from equilibrium chemistry can be expected. The derived disequilibrium abundances should not be taken as exact, however.

Photodissociation data

Simulation of photodissociation processes on warm or hot exoplanets require inputs that are relevant to these temperatures. However, measurements of such data, including cross-sections and branching ratios, are mostly available at around 300 K. While neglecting temperature-dependent cross-sections could change the compositions at the upper atmosphere of exoplanets profoundly (Venot et al. 2018b), we showed that on the dayside of UHJs, photodissociation plays a negligible role. Because the data we used to model photodissociation have large uncertainties, small differences in the stellar spectra (whether synthetic or observed) do not strongly affect the atmospheric composition either.

Photochemical hazes

While the formation of photochemical hazes is poorly understood, it is typically assumed (e.g., Rimmer et al. 2013; Gao et al. 2017; Hörst et al. 2018; Edgington et al. 2018; He et al. 2018; West 2018; Kawashima & Ikoma 2018; Helling 2019; Berry et al. 2019; Krasnopolsky 2020) that the formation of hydrocarbon-based hazes starts with the photolysis of haze precursor molecules such as CH₄, C₂ H₂, HCN, and C₆ H₆. Except for C₆ H₆, all of these species are included in our reaction network. The abundances of these species can thus approximate the importance of photochemical hazes. Using a combination of a kinetics model and an aerosol microphysics model, Lavvas & Koskinen (2017) found that hazes can be efficiently formed on some of the colder hot Jupiters, such as HD 189733b, with an equilibrium temperature of ~1100 K. On the hotter HD 209458b, Teq ≈1300 K, photochemical hazes did not form efficiently enough to affect the predicted transit spectra. Comparing the dayside abundances of CH₄, C₂ H₂, and HCN from our kinetic model, we find that the abundances of each of these haze precursor molecules are either comparable to or significantly lower than the abundances in the HD 209458b simulation. We thus conclude that photochemical hazes are even less important for HAT-P-7b. Furthermore, temperatures exceeding 2700 K in the upper atmosphere on large parts of the dayside preclude even hazes with properties similar to soot particles in these regions because they are not stable at these temperatures (Rimmer et al. 2013).

6 Conclusions

We investigated the importance of chemical disequilibrium (vertical mixing and photochemistry) on an ultra-hot Jupiter, HAT-P-7b, as an example of this class of planets. We employed two approaches to investigate the quenching levels: a zeroth-order approximation based on the results of Venot et al. (2018a), and full kinetic modeling. While both approaches usually agree with each other, the zeroth-order approximation underestimates the quenching levels at high temperatures by about one order of magnitude.

By performing 1D chemical kinetic simulations for 38 molecules and at 97 latitude and longitude points on a 3D GCM model, we find that the disequilibrium is a local phenomena and affects the nightside and the morning terminator of the most for specific molecules such as H, H₂O, CH₄, CO₂, HCN, and all C_n H_m molecules; higher mass C_n H_m molecules are more affected by disequilibrium processes. For these molecules, the variation of abundances on the dayside is also noticeable, but the evening terminator is least affected by disequilibrium processes. This may partially explain why the observed transmission spectra of UHJs are consistent with thermochemical equilibrium.

We find that the CO abundance is only marginally affected by vertical mixing, whereas the abundance of other major opacity species such as CH₄, H₂O, and also HCN are influenced by these processes. Photochemistry has almost no effect on the abundance of the species considered in this study.

Based on the locality characteristic of disequilibrium processes, we propose four avenues to search for the effect of these processes on UHJs. All four methods are based on comparative retrieval, where local atmospheric properties are retrieved and compared. We caution, however, that clouds will have to be taken into account for any interpretation. In the first method we suggest to investigate vertical opacity distributions between the morning terminator (in transmission) and the substellar point (in emission). The second method is to determine the differences observed in transmission spectra between the morning and evening terminators. The third method is to measure C/O at the two terminators and search for a lack of asymmetry (C/O may deviate by 0.2 between the two limbs when only cloud formation is considered). A similar C/O at the two terminators might indicate that disequilibrium processes are in action. An enhancement of global observed C/O (by about 0.05 in the case of HAT-P-7b) might also be an indication of cloud formation on this planet. The challenging nature of these methods means that the best chance to observe the effect of clouds and disequilibrium chemistry with the current facilities is probably precise phase curves as the fourth avenue. Through such precision photometry, a longitudinal chemical map of the planet might be constructed, as reported by Tan & Komacek (2019), Arcangeli et al. (2019), Mansfield et al. (2020), and Wong et al. (2019). Finally, it is worth mentioning that the 3D atmospheric structure of HAT-p-7b and possible geometric expansion of the dayside over the nightside, as described by Caldas et al. (2019), could make the retrieval of these parameters even more challenging. We investigate and discuss the observability of the atmosphere of HAT-P-7b and the effects of disequilibrium on it in a forthcoming paper in detail.

Acknowledgements

We acknowledge Cloud Academy I 2018 at the Les Houches School of Physics, France, during which this project was kicked off. We thank Maria Steinrueck for the constructive discussions and contribution. Part of this work was supported by the German Deutsche Forschungsgemeinschaft, DFG project number Ts 17/2–1. We are grateful to the anonymous referee for helpful comments and constructive criticism of the manuscript.

Appendix A Supplementary figures

Additional figures are provided in this section. These figures include initial and final abundances as well as their differences for CH₄, H₂O, CH, CN, C₂ H₂, and HCN, the transection maps of temperature and vertical velocities, and C/O ratio maps in terminator slice representation.

Fig. A.1 Methane (CH4): left column: equatorial slices (θ = 0°), right column: terminator slices. Top: initial condition, . Middle: finalresult from the kinetic simulation, . Bottom: difference, between the initial (equilibrium) values and the kinetic results.

Fig. A.2 Water (H2O): left column: equatorial slices (θ = 0°), right column: terminator slices. Top: initial condition, . Middle: finalresult from the kinetic simulation, . Bottom: difference, between the initial (equilibrium) values and the kinetic results.

Fig. A.3 Methylidyne (carbyne) (CH): left column: equatorial slices (θ = 0°), right column: terminator slices. Top: initial condition, . Middle: finalresult from the kinetic simulation, . Bottom: difference, between the initial (equilibrium) values and the kinetic results.

Fig. A.4 Cyanide (CN): left column: equatorial slices (θ = 0°), right column: terminator slices. Top: initial condition, . Middle: finalresult from the kinetic simulation, . Bottom: difference, between the initial (equilibrium) values and the kinetic results.

Fig. A.5 Acetylene (C2H2): left column: equatorial slices (θ = 0°), right column: terminator slices. Top: initial condition, . Middle: finalresult from the kinetic simulation, . Bottom: difference, between the initial (equilibrium) values and the kinetic results.

Fig. A.6 Hydrogen cyanide (HCN): left column: equatorial slices (θ = 0°), right column: terminator slices. Top: initial condition, . Middle: finalresult from the kinetic simulation, . Bottom: difference, between the initial (equilibrium) values and the kinetic results.

Fig. A.7 Transection maps of temperature (K) and vertical velocities (w) obtained from the GCM solution.

Fig. A.8 Carbon-to-oxygen ratio (C/O) in terminator slices representation: top left: C/O after cloud formation as in Paper I according to ɛC∕ɛO, top right: C/O for initial condition of gas-kinetic calculations (from CH, CO, CO2, C2 H2, etc. Bottom left: C/O at end of gas-kinetic calculation (t = 300 yr), bottom right: C/O change between initial condition (top right) and final gas-kinetic solution (bottom left). See Fig. 10 for terminator slice representation.

Appendix B Approximation of molecular diffusion coefficients

We described our chemical kinetic modeling in Sect. 2.2. Here we expand on how we approximated the molecular diffusion coefficients using the Chapman and Enskog equation, Eq. (4), (Poling et al. 2000). We recall (B.1)

where M_AB is defined as (M_A and M_B being the molecular weights of species A and B, respectively), Ω_D is the collision integral for diffusion, σ_AB is the characteristic length of the intermolecular force law, f_D is a correction term (usually on the order of unity), n is the number density of molecules in the mixture, k is the Boltzmann constant, and T is the gas temperature of the mixture. Of these, Ω_D and σ_AB demand additional steps to be approximated.

Ω_D can be estimated following the relation of Neufeld et al. (1972), (B.2)

where A = 1.06036, B = 0.15610, C = 0.19300, D = 0.47635, E = 1.03587, F = 1.52996, G = 1.76474, and H = 3.89411, and T^* = kT∕ε_AB. ε_AB can be estimated as (B.3)

where ε_A and ε_B are the characteristic Lennard-Jones energy of molecules A and B, respectively, whose values are usually tabulated and known (see, e.g., Poling et al. 2000, Appendix B).

Similarly, σ_AB can be estimated as (B.4)

where σ_A and σ_B are the characteristic Lennard-Jones lengths of molecules A and B, respectively, whose values are also tabulated (see, e.g., Poling et al. 2000, Appendix B). With these approximations, the molecular diffusion coefficients can be estimated in the chemical kinetic models.

References

All Figures

	Fig. 1 Estimation of the eddy coefficient needed for quenching to occur as a function of pressure and temperature, K(p,T). Venot et al. (2018a) reported the quenching pressures for several pairs of temperature-vertical mixing strength (in their Fig. 1), and we interpolate their results to obtain a continuous K function at any chosen pressure and temperature. See Sect. 2.1 for more details.
In the text

Fig. 2 Calculated eddy diffusion coefficient (Kzz) profiles using Eq. (1). The average eddy diffusion coefficient at pressure higher than 1 bar roughly follows a simple relation as log10(Kzz) = 5.5–1.7log10(pgas). Above this pressure and on the dayside, Kzz has average values about 2 × 107cm2s−1 (vertical reddashed line), and it is weaker by about one order of magnitude on the nightside (vertical blue dashed line).

In the text

	Fig. 3 Left: calculated Kzz using Eq. (1). Right: calculated quenching levels by following a zeroth-order approach. In the blue regions the atmosphere is in thermochemical equilibrium, and orange highlights the regions that are expected to be affected by disequilibrium processes and at least one major species deviates from its thermochemical equilibrium abundance.
In the text

Fig. 4 Variation of quenching pressure with the zenith angle, using the zeroth-order approximation (Sect. 3.1). Temperature profiles of all trajectories are color-coded above the quenching pressures. Minimum quenching pressure occurs at theevening terminator, where the atmosphere contributes significantly in the transmission spectra. The effects of disequilibrium chemistry are therefore expected to be less pronounced in the transmission spectra. Variation of quenching pressure is represented with the dashed line.

In the text

Fig. 5 a, d, g: temporal variation of the H2O, CO2, and CH4 vertical abundances at the equator and longitude = 225°, i.e., nightside. b, e, h: gCV of abundances, from initial values to a steady state at ~ 5 × 109 s. A higher value represents a stronger variation of the abundance at any given pressure level. A choice of gCV = 0.05 seems to reasonably separate the chemical equilibrium and disequilibrium regions (dashed gray vertical lines). c, f, i: quenching timescales. An invariant abundance profile at the longer times ensures a new steady state. See Sect. 3.2 for more details. We note that the quenching level cannot be associated with one particular pressure for all species, but usually spans more than one order of magnitude in pressure.

In the text

Fig. 6 Estimation of quenching levels using the chemical kinetic modeling (Sect. 3.2). Top panels: transection maps of the geometric gCV of H2O, CH4, and all species abundances at the equator (latitude = 0°). We choose gCV = 0.05 to separate the chemical equilibrium and disequilibrium regions (dashed white lines). This boundary agrees well with the zeroth-order approximation results (red solid lines in the right panels; Fig. A.7). Middle panels: similar to the top panels, but at the noon-midnight plane (longitude = 0° or 180°). Bottom panels: similar to the top panels, but at the terminator (longitude = 90° and 270°). map suggests a morning-evening asymmetry where the dawn-side abundances have been driven more significantly out of thermochemical equilibrium and at deeper pressure levels.

In the text

	Fig. 7 Top: latitude and longitude dependency of quenching pressures, calculated from the kinetic simulations including only diffusion. The quenching criterion is discussed in Sect. 3.2 Middle: similar to the top panel, but photochemistry is also considered. Bottom: difference between the quenching pressures of models with diffusion alone and models including both diffusion and photochemistry.
In the text

Fig. 8 Top: latitude and longitude dependency of quenching pressures, calculated using Kzz = w ⋅ H × 10−2. Middle: similar to the top panel but using Kzz = w ⋅ H, i.e. traditional approach. Bottom: comparing the quenching pressures of models with Kzz = w⋅ H to the quenching pressures of models with Kzz = w ⋅ H × 10−2. The atmosphere is being quenched at deeper levels, usually 20–40% on the dayside and 40–80% on the nightside. Some terminator trajectories, however, show several orders of magnitude change in their quenching levels, deepening from the mbar to 1 bar pressures.

In the text

	Fig. 9 Carbon monoxide (CO): left column: equatorial slices (θ = 0°), right column: terminator slices. Top: initial condition, . Middle: finalresult from the kinetic simulation, . Bottom: difference, between the initial (equilibrium) values and the kinetic results.
In the text

	Fig. 10 Carbon-to-oxygen ratio (C/O): top left: C/O after cloud formation as in Paper I according to ɛC ∕ɛO, Top right: C/O for initial condition of gas-kinetic calculations (from CH, CO, CO2, C2 H2, etc. Bottom left: C/O at the end of gas-kinetic calculations (t = 1012 sec) Bottom right: C/O change between initial condition (top right) and final gas-kinetic solution (bottom left).
In the text

Fig. 11 Prominent sources of molecular opacity under disequilibrium chemistry. The atmospheric pressure at which the vertical wavelength-dependent optical depth of species x reaches unity, p(τx(λ) = 1), is shown across the infrared. Several molecules included in both the kinetic model and the opacity database that do not reach τx(λ) = 1 by the 10 barlevel (e.g., C2H2) are omitted. See Paper I (Fig. 20) for opacity contributions due to species not included in the kinetic model (heavy-metal oxides, atoms, etc.). All opacities are plotted at a spectral resolution of R = 1000. Four equatorial regions are considered: (i) the antistellar point, (ii) the morning terminator, (iii) the substellar point, and (iv) the evening terminator. For comparison, the τ(λ) = 1 surface due to clouds from Paper I is overlaid.

In the text

Fig. 12 Differences to the opacity contribution of observable molecular species due to disequilibrium chemical kinetics. The difference metric, , is defined such that a positive difference implies stronger absorption after accounting for disequilibrium chemistry (i.e., p(τx (λ) = 1) occurs higher in the atmosphere). Only H2O, CO, and CO2 are plotted because these molecules offer the greatest detection prospect in the near-infrared (see Fig. 11). All opacities are plotted at a spectral resolution of R = 1000. Four equatorial regions are considered: (i) the antistellar point, (ii) the morning terminator, (iii) the substellar point, and (iv) the evening terminator.

In the text

	Fig. A.1 Methane (CH4): left column: equatorial slices (θ = 0°), right column: terminator slices. Top: initial condition, . Middle: finalresult from the kinetic simulation, . Bottom: difference, between the initial (equilibrium) values and the kinetic results.
In the text

	Fig. A.2 Water (H2O): left column: equatorial slices (θ = 0°), right column: terminator slices. Top: initial condition, . Middle: finalresult from the kinetic simulation, . Bottom: difference, between the initial (equilibrium) values and the kinetic results.
In the text

	Fig. A.3 Methylidyne (carbyne) (CH): left column: equatorial slices (θ = 0°), right column: terminator slices. Top: initial condition, . Middle: finalresult from the kinetic simulation, . Bottom: difference, between the initial (equilibrium) values and the kinetic results.
In the text

	Fig. A.4 Cyanide (CN): left column: equatorial slices (θ = 0°), right column: terminator slices. Top: initial condition, . Middle: finalresult from the kinetic simulation, . Bottom: difference, between the initial (equilibrium) values and the kinetic results.
In the text

	Fig. A.5 Acetylene (C2H2): left column: equatorial slices (θ = 0°), right column: terminator slices. Top: initial condition, . Middle: finalresult from the kinetic simulation, . Bottom: difference, between the initial (equilibrium) values and the kinetic results.
In the text

	Fig. A.6 Hydrogen cyanide (HCN): left column: equatorial slices (θ = 0°), right column: terminator slices. Top: initial condition, . Middle: finalresult from the kinetic simulation, . Bottom: difference, between the initial (equilibrium) values and the kinetic results.
In the text

	Fig. A.7 Transection maps of temperature (K) and vertical velocities (w) obtained from the GCM solution.
In the text

Fig. A.8 Carbon-to-oxygen ratio (C/O) in terminator slices representation: top left: C/O after cloud formation as in Paper I according to ɛC∕ɛO, top right: C/O for initial condition of gas-kinetic calculations (from CH, CO, CO2, C2 H2, etc. Bottom left: C/O at end of gas-kinetic calculation (t = 300 yr), bottom right: C/O change between initial condition (top right) and final gas-kinetic solution (bottom left). See Fig. 10 for terminator slice representation.

In the text