Movatterモバイル変換

• Top
• Abstract
• 1 Introduction
• 2 Method
• 3 Results
• 4 Discussion
• 5 Conclusions
• Acknowledgements
• Appendix A
• References
• List of tables
• List of figures

© ESO 2021

1 Introduction

Since the first discovery of an exoplanet in 1995 (Mayor & Queloz 1995), atmospheric spectra have been observed for dozens of exoplanets via transmission, emission, and direct imaging spectroscopy by both space- and ground-based telescopes. From these spectra, the abundances of several chemical species have been determined for some of the observed planets, most of which are expected to be abundant in hydrogen and helium-dominated atmospheres, such as H₂O, CH₄, CO, and CO₂ (e.g.,Line et al. 2014).

Recently, the use of spectral retrieval models has become standard when deriving the atmospheric properties from observed spectra (e.g.,Irwin et al. 2008;Benneke & Seager 2012;Line et al. 2013;Waldmann et al. 2015;Fisher & Heng 2018;Mollière et al. 2019;Min et al. 2020). Since the computational cost of the spectral retrieval codes is quite high, some assumptions and/or simplifications have to be made for atmospheric chemistry. For simplicity, most of the current spectral retrieval models assume constant abundances of chemical species throughout the atmosphere or their chemical equilibrium abundances, ignoring any disequilibrium chemistry.

In reality, however, the abundances of chemical species in the atmospheres cannot always be determined from the equilibrium chemistry, as has been studied by various 1D (photo-)chemical models (e.g.,Kasting et al. 1985;Moses et al. 2011;Venot et al. 2012;Hu et al. 2012;Grassi et al. 2014;Rimmer & Helling 2016;Tsai et al. 2017;Kawashima & Ikoma 2018) and 2D or 3D atmospheric circulation models with chemistry (e.g.,Cooper & Showman 2006;Agúndez et al. 2012,2014;Drummond et al. 2018b,a,2020;Mendonça et al. 2018). Among the chemical disequilibrium processes, such as photochemistry, the vertical quenching effect due to eddy diffusion transport is regarded as one of the most important because it significantly affects the atmospheric region where we can probe via observation. Indeed,Baxter et al. (2021) recently analyzed the transit depths of 49 gas giants measured at the 3.6 and 4.5 μm bands ofSpitzer Space Telescope and found evidence of disequilibrium chemistry, which can be explained by the quenching effect of CH₄. Quenching happens at altitudes where the thermochemical reaction timescale,τ_chem, is equal to that of the eddy diffusion,τ_diff. Below this altitude, whereτ_chem <τ_diff, the volume mixing ratio of a chemical species is consistent with that in thermochemical equilibrium due to high temperature and number density. On the other hand, above that altitude, whereτ_diff <τ_chem, the abundance is “frozen” to the value at the quenching altitude since vertical transport due to eddy diffusion tends to smooth out the gradient of the volume mixing ratio. Thus, deriving the metallicity or elemental abundance ratios of the atmosphere from those frozen molecular abundances without considering the quenching effect can lead to over- or under-estimates of the metallicity or elemental abundance ratios. Because of the longer chemical timescale at lower temperatures, this quenching effect is especially important for cooler (≲1000 K) atmospheres, which are the primary targets for upcoming characterizations of exoplanet atmospheres.

Morley et al. (2017) introduced a quenching pressure as a parameter for spectral retrieval, assuming the same quenching pressures for all species. This treatment of the quenching effect was also adopted inMollière et al. (2020). However, the quenching pressure should be different for each species because the chemical timescale is different for each species, especially between CH₄/CO and NH₃/N₂ (e.g.,Moses 2014). Also, as a more general approach for capturing discontinuities in the abundance profiles, caused not only by vertical mixing but also by photodissociation and the formation of clouds and haze,Changeat et al. (2019) recently proposed a “two-layer” retrieval approach that allows for different constant abundances in the upper and lower atmospheres.

Predicting the quenching pressure for each species with a physical basis involves calculating chemical reactions under the effect of eddy diffusion transport. It is, however, unrealistic to couple full kinetic chemistry to spectral retrieval codes, even for 1D modeling, due to its high computational cost. As such, we can refer to previous works that attempted to couple chemistry with computationally expensive 2D or 3D atmospheric circulation models to study the disequilibrium effect due to atmospheric circulation.Cooper & Showman (2006) adopted the chemical relaxation method, which has since been used in various studies (e.g.,Drummond et al. 2018b,a;Mendonça et al. 2018). The chemical relaxation scheme replaces numerous chemical production and loss terms in the continuity equation of a species with a single term given by the deviation of the abundance from the equilibrium value divided by its chemical timescale (Smith 1998;Cooper & Showman 2006). Other ways to reduce the computational cost have also been adopted, such as simplifying the atmospheric dynamics model (Agúndez et al. 2012,2014) or using the reduced chemical network (Drummond et al. 2020).

When using the chemical relaxation method, adopting an appropriate chemical timescale for each species is important. The chemical timescale can be approximated by that of the rate-limiting or the slowest reaction along the conversion pathway from one species to another. Previous works have investigated the rate-limiting reaction for the conversion, such as that from CH₄ to CO and/or that from NH₃ to N₂ (e.g.,Prinn & Barshay 1977;Yung et al. 1988;Visscher et al. 2010;Moses et al. 2011;Zahnle & Marley 2014). To revisit the chemical relaxation method,Tsai et al. (2018) derived the chemical timescales of the major species in hydrogen and helium-dominated atmospheres valid for the wide pressure and temperature ranges of currently observable exoplanet atmospheres, namely temperatures from 500 K to 3000 K and pressures from 0.1 mbar to 1000 bar.

In this paper, we adopt the chemical relaxation method to incorporate the effect of vertical mixing or quenching for each major species in hydrogen and helium-dominated atmospheres into the spectral retrieval code ARCiS (ARtful modelling Code for exoplanet Science;Min et al. 2020) with a physical basis. This enables us to directly retrieve the eddy diffusion coefficient of exoplanet atmospheres, unlike for the various previous spectral retrieval codes. We note that we consider the quenching effect only as a disequilibrium process since this is what most affects the composition of the atmospheric region that we can probe by transmission spectroscopy in the optical and infrared. We neglect photochemistry, which is nevertheless important in the upper atmosphere. We use the chemical timescales derived byTsai et al. (2018) as they were validated for the broad pressure and temperature ranges relevant to exoplanet atmospheres. Then, using the updated ARCiS code with the consideration of disequilibrium chemistry, we perform retrievals of the observed transmission spectra of exoplanets with sizes ranging from that of Jupiter to mini-Neptunes and investigate the trend of disequilibrium chemistry in exoplanet atmospheres.

The rest of this paper is organized as follows. In Sect.2, we describe our method and simulation setup. In Sect.3, we present the results. In Sect.4, we discuss several points that should be addressed in future works. Finally, we conclude this paper in Sect.5.

2 Method

2.1 Model description

The 1D continuity-transport equation for the number density of speciesi,n_i, is written as$$\frac{\partial n_i}{\partial t}=P_i-L_i-\frac{\partial {\Phi }_i}{\partial z},$$(1)

wheret andz are the time and altitude, respectively,P_i andL_i are the production and loss rates ofn_i due to chemical reactions, respectively, and Φ_i is the vertical transport flux. Assuming that eddy diffusion is the dominant transport mechanism, which is usually the case for the atmospheric region of exoplanets that we can observe in the optical and infrared, Φ_i is given as$${\Phi }_i=-K_{\text{zz}}N\frac{\partial f_i}{\partial z}.$$(2)

Here,K_zz is the eddy diffusion coefficient, which we assume to be constant throughout the atmosphere for simplicity, andN andf_i are the total number density of the gaseous species and the volume mixing ratio of speciesi, namelyf_i =n_i∕N, respectively.

The chemical relaxation method replaces the chemical source and sink terms ofP_i −L_i by the deviation from the equilibrium number density divided by the chemical timescale (Smith 1998;Cooper & Showman 2006). In that case, Eq. (1) can be rewritten as$$\frac{\partial n_i}{\partial t}=-\frac{n_i-n_{i,\text{eq}}}{{\tau }_{i,\text{chem}}}-\frac{\partial {\Phi }_i}{\partial z}.$$(3)

Here,n_i,eq andτ_i,chem denote the equilibrium number density and chemical timescale of speciesi, respectively.

If we assume a steady-state condition, Eq. (3) can be further transformed as$$n_i+{\tau }_{i,\text{chem}}\frac{\partial {\Phi }_i}{\partial z}=n_{i,\text{eq}}.$$(4)

Next, we discretize Eq. (4) using the subscriptj for thejth altitude layer, which we assume to have the same thickness of Δz as the other layers, as$$n_{i,j}+{\tau }_{i,j,\text{chem}}\frac{{\Phi }_{i,j+\frac{1}{2}}-{\Phi }_{i,j-\frac{1}{2}}}{\Delta z}=n_{i,j,\text{eq}}.$$(5)

Here, we approximate the diffusion fluxes at the boundary between thej − 1th andjth layers and that between thejth andj + 1th layers,${\Phi }_{i,j-\frac{1}{2}}$ and${\Phi }_{i,j+\frac{1}{2}}$, as$${\Phi }_{i,j-\frac{1}{2}}=-K_{\text{zz}}{N}_{j-\frac{1}{2}}\frac{f_{i,j}-f_{i,j-1}}{\Delta z}$$(6)

and$${\Phi }_{i,j+\frac{1}{2}}=-K_{\text{zz}}{N}_{j+\frac{1}{2}}\frac{f_{i,j+1}-f_{i,j}}{\Delta z},$$(7)

where$N_{j-\frac{1}{2}}$ and$N_{j+\frac{1}{2}}$ are the total number densities at those boundaries.

Finally, we can convert Eq. (5) into the matrix form as$$(E+T_i{M}_i)n_i=n_{i,\text{eq}}.$$(8)

Here,E is the identity matrix, and the other matrices are given as follows with the number of the altitude layers$\mathcal{N}$. “diag” indicates a diagonal matrix.$$n_i={\left(n_{i,1},n_{i,2},n_{i,3},\cdots ,n_{i,j},\cdots ,n_{i,\mathcal{N}}\right)}^{\text{T}}$$(9)$$n_{i,\text{eq}}={\left(n_{i,1,\text{eq}},n_{i,2,\text{eq}},n_{i,3,\text{eq}},\cdots ,n_{i,j,\text{eq}},\cdots ,n_{i,\mathcal{N},\text{eq}}\right)}^{\text{T}}$$(10)$$T_i=\text{diag}\left({\tau }_{i,1,\text{chem}},{\tau }_{i,2,\text{chem}},{\tau }_{i,3,\text{chem}},\cdots ,{\tau }_{i,j,\text{chem}},\cdots ,{\tau }_{i,\mathcal{N},\text{chem}}\right)$$(11)$$M_i=-K_{\text{zz}}\left(\begin{array}{ccccccccc}0& 0& 0& \cdots & \cdots & \cdots & \cdots & \cdots & 0\\ \frac{1}{\Delta z_2{}^2}\frac{N_{\frac{3}{2}}}{N_1}& -\frac{1}{\Delta z_2{}^2}\frac{N_{\frac{3}{2}}+N_{\frac{5}{2}}}{N_2}& \frac{1}{\Delta z_2{}^2}\frac{N_{\frac{5}{2}}}{N_3}& 0& \cdots & \cdots & \cdots & \cdots & 0\\ 0& \frac{1}{\Delta z_3{}^2}\frac{N_{\frac{5}{2}}}{N_2}& -\frac{1}{\Delta z_3{}^2}\frac{N_{\frac{5}{2}}+N_{\frac{7}{2}}}{N_3}& \frac{1}{\Delta z_3{}^2}\frac{N_{\frac{7}{2}}}{N_4}& 0& \cdots & \cdots & \cdots & 0\\ ⋮& & & \ddots & & & & & ⋮\\ 0& \cdots & 0& \frac{1}{\Delta z_j{}^2}\frac{N_{j-\frac{1}{2}}}{N_{j-1}}& -\frac{1}{\Delta z_j{}^2}\frac{N_{j-\frac{1}{2}}+N_{j+\frac{1}{2}}}{N_j}& \frac{1}{\Delta z_j{}^2}\frac{N_{j+\frac{1}{2}}}{N_{j+1}}& 0& \cdots & 0\\ ⋮& & & & & \ddots & & & ⋮\\ ⋮& & & & & & \ddots & & ⋮\\ ⋮& & & & & & & \ddots & ⋮\\ 0& \cdots & \cdots & \cdots & \cdots & \cdots & 0& \frac{1}{\Delta z_{\mathcal{N}}{}^2}\frac{N_{\mathcal{N}-\frac{1}{2}}}{N_{\mathcal{N}-1}}& -\frac{1}{\Delta z_{\mathcal{N}}{}^2}\frac{N_{\mathcal{N}-\frac{1}{2}}}{N_{\mathcal{N}}}\end{array}\right)$$(12)

In the above, we adopted the lower boundary condition of the chemical equilibriumn_i,1 =n_i,1,eq because of the short chemical timescale in the deeper atmosphere. For the upper boundary condition, we adopted zero flux, namely without atmospheric escape. By multiplying both sides of Eq. (8) by the inverse matrix of (E +T_iM_i), we finally get the solution as$$n_i={(E+T_i{M}_i)}^{-1}{n}_{i,\text{eq}}.$$(13)

We note that this equation leads ton_i ~n_i,eq, namely in equilibrium, whenK_zz has a negligible value.

2.2 Implementation to ARCiS

We implemented this module to the spectral retrieval code ARCiS (Min et al. 2020) using the chemical timescales of CH₄, CO, H₂O, NH₃, and N₂ fromTsai et al. (2018). We note that this disequilibrium calculation part takes almost negligible time (~ 0.01% of the total computation time for one forward model). We also note that for the chemical timescales of CH₄, CO, and H₂O, their different expressions above and below the C/O ratio of unity are used. CO₂ tends to remain in pseudo-equilibrium even after its related molecules CO and H₂O are quenched (Moses et al. 2011;Tsai et al. 2018). Thus, for the calculation of CO₂ abundance, instead of using Eq. (12), we adopt the following pseudo-equilibrium abundance formula ofTsai et al. (2018),$$n_{{\text{CO}}_2}=\frac{n_{\text{CO}}{n}_{{\text{H}}_2\text{O}}}{n_{\text{CO},\text{eq}}{n}_{{\text{H}}_2\text{O},\text{eq}}}{n}_{{\text{CO}}_2,\text{eq}},$$(14)

which modifies the CO₂ number density based on the quenched abundances of CO and H₂O when they are quenched. All the species except for the above are assumed to have their equilibrium abundances.

ARCiS uses GGchem (Woitke et al. 2018) for the calculation of thermochemical equilibrium abundances. The original version of GGchem did not include CH₂OH, CH₃OH, or N₂H₃, the abundances of which are needed to calculate the chemical timescale of the species mentioned above. Taking their thermodynamic data fromBurcat & Ruscic (2005)¹, we added these species to GGchem, available on the GitHub page of GGchem².

2.3 Application to the observed transmission spectra of 16 planets

With the updated ARCiS, we performed spectral retrievals where disequilibrium chemistry is allowed to affect the profiles of the molecular abundances (hereafter, “disequilibrium retrieval”). To investigate the effect of the inclusion of disequilibrium chemistry to the retrieval code, we also performed retrievals imposing equilibrium chemistry (hereafter, “equilibrium retrieval”) and compared the results of the two retrievals.

For the planet samples, we selected ten hot Jupiters and six Neptunes compiled inSing et al. (2016) andCrossfield & Kreidberg (2017), respectively. High-precision transmission spectra were observed by theHubble Space Telescope (HST) andSpitzer for those planets. The samples are listed inTable 1, along with the references for their observed transmission spectrum data used in the retrieval. These planet samples range from clear to cloudy atmospheres, and we note that the hot Jupiter samples and their spectral data are the same as those used inMin et al. (2020) except for HAT-P-12b, for which we adopt the observation data recently analyzed byWong et al. (2020).

We adopt the same settings as the “constrained retrieval” ofMin et al. (2020), which considered the thermal and cloud structures with the models ofGuillot (2010) andOrmel & Min (2019), respectively, and adopted the precomputed opacities presented inChubb et al. (2021). The retrieval parameters and the employed ranges and priors of those parameters are presented inTable 2. They are also the same as inMin et al. (2020), except for the eddy diffusion coefficient used in the disequilibrium chemistry module,$K_{zz}^{\text{Ch}}$, which was not considered in the previous retrieval. In this study, as a first step to consider disequilibrium chemistry in spectral retrieval, we define a “standard” disequilibrium retrieval as a retrieval treating the eddy diffusion coefficient used in the chemistry calculation,$K_{zz}^{\text{Ch}}$, and that used in the cloud simulation,$K_{zz}^{\text{Cl}}$, independently. This is because the atmospheric region where clouds exist and the region where chemical species affecting the spectra are located are usually different, and hence the two eddy diffusion coefficients reflect the coefficients at different parts of the atmosphere. Moreover, we also examine the effect of imposing the condition that the two eddy diffusion coefficients are the same, namely$K_{zz}^{\text{Ch}}=K_{zz}^{\text{Cl}}$ (hereafter, “same-K_zz disequilibrium retrieval”). Retrieval with a nonconstant coefficient along the altitude is left as a future work. We note that for the retrieval of WASP-12b, in the same way asMin et al. (2020), we adopt an additional parameter, which allows the scaling of the data fromKreidberg et al. (2015) to match the remaining data points due to the existing slight offset of their data when compared to the other data (Sing et al. 2013,2016). (SeeMin et al. 2020 for a detailed treatment of this additional parameter.)

Our disequilibrium retrieval can also treat the equilibrium case since a negligible value ofK_zz leads to the molecular abundances in equilibrium, as indicated in Eq. (12). We confirmed that the lower bound of$K_{zz}^{\text{Ch}}$ we adopted, namely 10⁰ cm² s⁻¹, is small enough to give molecular abundances consistent with the equilibrium values for the hotter planets of our samples. For some cooler planets, however, a few species deviate from their equilibrium abundances, even for$K_{zz}^{\text{Ch}}={10}^0$ cm² s⁻¹. This is because the chemical timescale becomes exceedingly long at lower temperatures (see Fig. 4 ofTsai et al. 2018). Considering a ~Gyr planetary age timescale,$K_{zz}^{\text{Ch}}={10}^0$ cm² s⁻¹ is approximately the limit for the vertical transport by eddy diffusion to affect molecular abundances within the planetary age. Here, we estimate the diffusion timescale asτ_diff ~H²∕K_zz using the atmospheric scale heightH of each of our samples. Thus, we set the lower bound of$K_{zz}^{\text{Ch}}$ to be 10⁰ cm² s⁻¹ even though the exploration range of$K_{zz}^{\text{Ch}}$ does not include the “true” equilibrium regime for a few minor species in relatively cool planets, such as GJ 436b, GJ 1214b, GJ 3470b, HAT-P-1b, HAT-P-11b, HAT-P-12b, HAT-P-26b, HD 97658b, HD 209458b, WASP-31b, and WASP-39b. We estimated this using the best-fit parameter set that yields the minimumχ² compared to the observed data from our standard disequilibrium retrieval calculations for each of our samples.

Finally, in addition to the aforementioned two disequilibrium retrievals, we also conduct a disequilibrium retrieval fixing one of the highly uncertain parameters, namely the intrinsic temperatureT_int, to an estimated value based on current knowledge. This is done based on the expectation that quenching often happens in the deeper atmosphere, especially for the nitrogen species, which means that the retrieval of$K_{zz}^{\text{Ch}}$ can degenerate withT_int. We assume the value ofT_int for each of our planet samples using its relation to the equilibrium temperatureT_eq at a thermal equilibrium inferred byThorngren et al. (2019), who derived the fractionϵ of the incident flux on the planetF that heats the interior sufficiently to reproduce the observed radii of their hot Jupiter samples as a function ofF. When assuming thermal equilibrium for the planet interiors, which will be reached in as little as tens of megayears, that fraction directly relatesT_eq andT_int as (Thorngren & Fortney 2018)$$\begin{array}{ll}{T}_{\text{int}}\hfill & =ϵ{\left(F\right)}^{1/4}{T}_{\text{eq}}\hfill \\ \hfill & ~0.39T_{\text{eq}}\exp \left(-\frac{{\left({\log }_{10}F-0.14\right)}^2}{1.095}\right),\hfill \end{array}$$(15)

where$F=4\sigma T_{\text{eq}}^4$ in units of Gerg cm⁻²s⁻¹ with the Stefan–Boltzmann constantσ. We adopt this assumption for this retrieval. Since their formula was inferred from planet samples with masses larger than half Jupiter mass (Thorngren & Fortney 2018), we apply this type of disequilibrium retrieval only to the planets that satisfy this criterion, namely HAT-P-1b, HD 189733b, HD 209458b, WASP-12b, and WASP-19b. For those planets, the adopted values ofT_int are 476, 376, 561, 517, and 653 K, respectively.

3 Results

3.1 Standard disequilibrium retrieval in comparison with equilibrium retrieval

In this section, we present the results of the standard disequilibrium retrieval and their differences from the equilibrium retrieval. As mentioned in Sect.2.3, for hotter planet samples, our disequilibrium retrieval can also treat the equilibrium cases. InFig. 1, we show the results of the spectra from the standard disequilibrium and equilibrium retrievals along with those of the disequilibrium retrieval with the sameK_zz and that with fixedT_int (only shownfor the planets with masses larger than half Jupiter mass), which we discuss in Sects.3.2 and3.3, respectively. Also, the retrieved pressure–temperature structure and the abundance profiles for four representative planets are shown in the left and right panels ofFig. 2, respectively. The profiles for the other planets are presented inFig. A.1 In addition, plots of the posterior distributions of the parameters from the sampling are presented inFig. B.1.

The second column ofTable 3 presents the natural logs of the Bayes factors between the standard disequilibrium and equilibrium retrievals. It can be seen that except for HD 209458b, and tentatively GJ 3470b and WASP-39b, the differences are almost negligible (cf.Trotta 2008), which means that the retrieved parameters from either of the two retrievals are barely favored over the parameters from the other retrieval. For HD 209458b and WASP-39b, the natural logs of the Bayes factors indicate that the disequilibrium scenario is favored over the equilibrium scenario by ≥ 4.1σ and ≥ 2.7σ, respectively, calculated with Eq. (27) ofTrotta (2008). For GJ 3470b, the value of the Bayes factor implies the preferability of the equilibrium scenario with ≥ 3.0σ. However, as mentioned in Sect.2.3, for this planet, a few species still deviate from their equilibrium abundances, even for the lower bound of the adopted$K_{zz}^{\text{Ch}}$ range, namely 10⁰ cm² s⁻¹. Thus, it is uncertain whether the “true” equilibrium scenario is indeed favored. Retrieval simulation extending the lower bound of$K_{zz}^{\text{Ch}}$ is needed to confirm this. For the retrieved spectra, thermal structures, abundance profiles, and parameters (seeFigs. 1,2,A.1, andB.1), HD 209458b and WASP-39b exhibit relatively large differences between the disequilibrium and equilibrium retrievals while they are largely similar within the uncertainty range for the other planets. This is consistent with the relatively large eddy diffusion coefficients retrieved for those two planets, such as${\log }_{10}\left(K_{zz}^{\text{Ch}} ({\text{cm}}^2 {\text{s}}^{-1})\right)={13.12}_{-4.61}^{+1.39}$ for HD 209458b and${8.61}_{-4.22}^{+3.68}$ for WASP-39b, inferring that disequilibrium chemistry is playing an essential role in their atmospheres.

For HD 209458b, while the retrieved spectra from both retrievals match relatively well in the wavelength range below 1.6 μm (seeFig. 1a), where precise observational data exists, differences are observed around 2.5–4.0, 6.5–9, and 9–13 μm. The discrepancies at the first two wavelength regions mainly come from the difference in the CH₄ abundance. A comparison of the pink thick solid line and pink thin solid line inFig. 2a shows that the CH₄ abundance is quenched in the deeper atmosphere with a pressure around 1 bar, resulting in a smaller abundance compared to the equilibrium retrieval case (pink thick dashed line). While quenching reduces the CH₄ abundance in the observable pressure region, it works to maintain a high NH₃ abundance (compare yellow thick solid line and yellow thin solid line inFig. 2a). This results in the prominent NH₃ absorption feature at 10.5 μm for the disequilibrium retrieval case (blue line inFig. 1a), while that feature hardly exists in the equilibrium retrieval case (orange line).

To understand which observed features require the condition of disequilibrium chemistry for HD 209458b, we additionally performed standard disequilibrium and equilibrium retrievals, systematically excluding each observed data set. Natural logs of the Bayes factors for those additional retrievals are presented inTable 4. It can be seen that the data of HST/WFC3/G141 (1.1–1.6 μm) require the disequilibrium chemistry condition most for the case of HD 209458b while it is also important for the HST/STIS data (0.3–0.9 μm). We confirmed that the best-fit parameter set, which yields the minimumχ² when compared to the observed data, from the equilibrium retrieval explains the observed 1.3–1.6 μm feature by H₂O alone. On the other hand, in the disequilibrium retrieval case, NH₃ also partially contributes to reproducing the feature, yielding a better match to the observed data and thus demonstrating the preferability of the disequilibrium retrieval. Here we mention that this weak evidence regarding the NH₃ feature was already noted byMacDonald & Madhusudhan (2017), who used the same data analyzed bySing et al. (2016) for their retrieval simulations and raised the possibility of disequilibrium chemistry playing a role in the atmosphere of this planet.

Recently,Giacobbe et al. (2021) performed high-resolution transmission spectroscopy of HD 209458b and reported the detection of six species, H₂O, CO, HCN, CH₄, NH₃, and C₂H₂. In their analysis, the disequilibrium scenario was strongly disfavored, though they did not deny the possibility that disequilibrium processes were in effect to some extent. This discrepancy may partly be because they assumed specific thermal andabundance profiles for their disequilibrium scenario while we allow those profiles to vary within our retrieval.

Next, regarding WASP-39b, noticeable differences are found around 0.3–0.5, 2.7, 4.3, and 15 μm (seeFig. 1b). In the case of the equilibrium retrieval (orange line), the favored cloud distributions contribute to flattening the Rayleigh scattering slope and fail to reproduce the observed steepness of the optical slope. The small difference at 2.7 μm arises from the smaller H₂O abundance in the disequilibrium retrieval case (compare the purple thick solid line and purple thick dashed line inFig. 2b). Due to the absence of CO₂ in the upper atmosphere (gray thick solid line in the same figure), its strong features at 4.3 and 15 μm hardly exist in the disequilibrium retrieval case (blue line inFig. 1b).

For the retrieved thermal profiles, although the number of samples is limited, a higher temperature in the lower atmosphere is favored for the disequilibrium retrieval, while the equilibrium retrieval prefers an isothermal-like profile for several planet samples, such as HD 209458b, WASP-39b, GJ 1214b, GJ 3470b, HAT-P-26b, and HD 97658b (seeFigs. 2 andA.1). We speculate that this is due to the absence of NH₃ and/or CH₄ absorption features for some of the above planets. This is because if moderate vertical diffusion is imposed for the cool lower atmosphere,such as the equilibrium retrieval cases of the above planets, their abundances would be quenched to the larger values because of their stability at low temperatures and high pressures. Thus, the temperature in the lower atmosphere needs to be increased to reduce their abundances.

The planets for which at least one of the retrieved parameters differs by more than 1σ between the standard disequilibrium and equilibrium retrievals are GJ 436b, GJ 1214b, GJ 3470b, HAT-P-11b, HAT-P-12b, HAT-P-26b, HD 209458b, WASP-6b, and WASP-39b, and their differences are summarized inTable 5. Aside from HD 209458b and WASP-39b, some differences are also found for HAT-P-11b and GJ 436b. For HAT-P-11b, the median values of C/O ratio and metallicity differ by 2.9 and 1.4σ, respectively. While the retrieved C/O ratio from the disequilibrium retrieval is${0.82}_{-0.42}^{+0.26}$, that from the equilibrium retrieval is${0.27}_{-0.12}^{+0.19}$. For the metallicity, the retrieved values from the disequilibrium and equilibrium retrievals are${\log }_{10}\left(Z (\text{dex})\right)={2.39}_{-0.34}^{+0.27}$ and${0.67}_{-0.68}^{+1.24}$, respectively. For this planet, even though the retrieved eddy diffusion coefficient is relatively small (${\log }_{10}\left(K_{zz}^{\text{Ch}} ({\text{cm}}^2 {\text{s}}^{-1})\right)={3.74}_{-2.29}^{+2.67}$), because of the somewhat low temperature, disequilibrium chemistry plays an important role (compare the thick and thin solid lines inFig. 2c). However, the quite small value of the natural log of the Bayes factor inTable 3 indicates that the retrieved parameters from the disequilibrium retrieval are never favored over those from the equilibrium case and vice versa. A further observational constraint is needed to determine the atmospheric properties of this planet. For this purpose, the search for the CO₂ feature at 4.3 μm is promising because of its absence in the disequilibrium retrieval case (Fig. 1c). For GJ 436b, while the irradiation parameter from the equilibrium retrieval is${0.01}_{-0.01}^{+0.01}$, that from the disequilibrium retrieval is${0.06}_{-0.04}^{+0.06}$, resulting in a 4.5σ difference of this parameter and different retrieved thermal structure (Fig. A.1e). However, as for the case of HAT-P-11b, the natural log of the Bayes factor is too small to draw a firm conclusion on the parameter. High-precision observations in the optical wavelength range, where the spectra from the two retrievals show a difference, could allow this to be distinguished. These examples raise the possibility that ignoring the effect of disequilibrium chemistry can lead to an incorrect constraint on the atmospheric properties.

Fig. 1 Transmission spectra retrieved from the standard disequilibrium retrieval (blue), equilibrium retrieval (orange), and disequilibrium retrievals with the sameKzz (green) and with fixedTint (red; only shown for planets with masses larger than half Jupiter mass). The solid lines show the median models while the shaded regions indicate the 1, 2, and 3σ confidence intervals. The black points denote the observed data used.

Fig. 1 continued. We note that for (m) WASP-12b, even thoughthe observed data ofKreidberg et al. (2015) are allowed to vertically shift in each retrieval simulation (see Sect.2.3 for details), their original values are plotted.

Fig. 2 Atmospheric profiles for four representative planets, while those for the other planets are presented inFig. A.1.Left: Pressure–temperature profiles retrieved from the standard disequilibrium retrieval (blue), equilibrium retrieval (orange), and disequilibrium retrievals with the sameKzz (green) and with fixedTint (red; only shown for planets with masses larger than half Jupiter mass). The solid lines show the median models while the shaded regions indicate the 1, 2, and 3σ confidence intervals.Right: Median abundance profiles of H2O, CO, and CH4 (left) and those of CO2, NH3, and N2 (right) from the standard disequilibrium retrieval (thick solid lines), equilibrium retrieval (thick dashed lines), and disequilibrium retrievals with the sameKzz (thick dotted lines) and with fixedTint (thick dash-dotted lines; only shown for planets with masses larger than half Jupiter mass). For reference, the median abundance profiles for the case where the equilibrium chemistry is imposed are shown with the corresponding thin lines.

3.2 Comparison of disequilibrium retrievals with and without imposing the sameK_zz for chemistry and cloud calculations

In this subsection, we compare the results from the standard disequilibrium retrieval and a retrieval where the condition$K_{zz}^{\text{Ch}}=K_{zz}^{\text{Cl}}$ is imposed to examine how this affects the retrieval of parameters. As shown in the third column ofTable 3, the absolute values of the natural logs of the Bayes factors between the two retrievals are almost negligible for all of our samples (cf.Trotta 2008), indicating that the retrieved parameters from either of the two retrievals are barely favored over those from the other.

The planets for which at least one of the retrieved parameters differs by more than 1σ between the same-K_zz and standard disequilibrium retrievals, and their differences are summarized inTable 6. As expected, employing the condition$K_{zz}^{\text{Ch}}=K_{zz}^{\text{Cl}}$ strongly affects the retrieval of these two parameters and the cloud nucleation rate. Large (≳ 3σ) differences are found for GJ 436b, GJ 1214b, HAT-P-12b, HD 189733b, HD 209458b, WASP-17b, WASP-19b, and WASP-39b. For all the planets, those differences are for$K_{zz}^{\text{Cl}}$ and$K_{zz}^{\text{Ch}}$ as well as$\dot{\Sigma }$ for HD 189733b.

For HD 209458b, for which we found an indication of disequilibrium chemistry from the comparison between the standard disequilibrium and equilibrium retrievals in Sect.3.1, the retrieved value of${\log }_{10}\left(K_{zz}^{\text{Ch}} ({\text{cm}}^2 {\text{s}}^{-1})\right)$ decreases from${13.12}_{-4.61}^{+1.39}$ to${4.17}_{-1.05}^{+1.15}$ with a 1.9σ difference. Despite this decrease, the spectrum from the same-K_zz retrieval (green line inFig. 1a) still exhibits certain disequilibrium features, as discussed in Sect.3.1, such as those at 2.5–4.0, 6.5–9, and 9–13 μm.

Also, for WASP-39b, for which we also found an indication of disequilibrium chemistry in Sect.3.1, the retrieved spectra, thermal structures, abundance profiles, and parameters (seeFigs. 1b,2b, andB.1b) are quite similar to the standard disequilibrium retrieval case. Thus, our findings of indications of disequilibrium chemistry for HD 209458b and WAS-39b still hold.

3.3 Comparison of disequilibrium retrievals with and without a fixed intrinsic temperature

In this subsection, we compare the results from the standard disequilibrium retrieval and that with a fixed intrinsic temperature. The motivation to fix the intrinsic temperature is based on the expectation that the retrieval of the chemical eddy diffusion coefficient might degenerate with the intrinsic temperature since the quenching often happens in the deeper atmosphere where the intrinsic temperaturehas a large influence on the thermal structure. Despite our expectation, however, such a correlation is not found (seeFig. B.1), which may be due to the limits of the current observational precision. The correlation coefficients between${\log }_{10}\left(T_{\text{int}} (K)\right)$ and${\log }_{10}\left(K_{zz}^{\text{Cl}} ({\text{cm}}^2 {\text{s}}^{-1})\right)$ are 0.12, 0.071, 0.17, 0.0037, and 0.10 for HAT-P-1b, HD 189733b, HD 209458b, WASP-12b, and WASP-19b, respectively. We note that for all the five targets to which we appliedT_int-fixed retrieval, the standard retrieval favors a lower value ofT_int than that inferred from the formula derived byThorngren et al. (2019) (seeTable 7). As shown in the fourth column ofTable 3, the absolute values of the natural logs of the Bayes factors between the two retrievals are quite small (cf.Trotta 2008) except for HD 189733b.

HD 189733b exhibits significantly weaker Bayesian evidence of disequilibrium retrieval with fixedT_int, indicating that the retrieved parameters from that retrieval are strongly (≳ 7.8σ) disfavored compared to those from the standard disequilibrium retrieval. The significant difference in the Bayesian evidence is partly due to the worse fit to the observed data points in the bluest part (0.3–0.5 μm) for theT_int-fixed retrieval (red line inFig. 1d). This discrepancy is due to the different retrieved parameters for the clouds. We note, however, that there is room for discussion regarding the quite steep optical slope of this planet (Pont et al. 2013). Several origins aside from the planetary atmosphere, such as an unknown systematic instrument offset between the optical and near-infrared or the presence of starspots on the host star (e.g.,Oshagh et al. 2014,2020) have been proposed as being capable of artificially producing such a steep slope in the planetary spectrum.

Comparing the retrieved spectra of HD 189733b from the two retrievals, differences can also be found above 8.0 μm. In the standard disequilibrium retrieval case (blue line), cloud features exist around 9.4 and 20 μm while the C₂H₂ absorption feature is visible at 13.6 μm in theT_int-fixed retrieval case (red line). A striking difference is also found in the retrieved thermal structure (seeFig. 2d). A strong thermal inversion is retrieved from theT_int-fixed retrieval (red line), whereas no such inversion is derived from the standard disequilibrium retrieval (blue line). We note that in our retrieval, we retrieved the parameters regarding the thermal structure independently from the atmospheric constituents. A retrieval employing a thermal structure consistent with the opacity of the atmospheric composition will be the subject of future work.

The retrieved parameters with more than 1σ difference between theT_int-fixed and standard disequilibrium retrievals are summarized inTable 7. Indeed, quite significant differences are found for the parameters of HD 189733b.

3.4 Trend of the retrieved parameters

In this subsection, we explore the trend of the atmospheric parameters derived from our retrieval calculations.

3.4.1 Eddy diffusion coefficient

Figure 3 shows the retrieved values of the eddy diffusion coefficients as a function of the equilibrium temperature for our samples (top panel for the chemical eddy diffusion coefficient,$K_{zz}^{\text{Ch}}$, and bottom panel for the cloud eddy diffusion coefficient,$K_{zz}^{\text{Cl}}$). First, neither$K_{zz}^{\text{Ch}}$ (blue circle points in thetop panel) nor$K_{zz}^{\text{Cl}}$ (purple triangle points in the bottom panel) from the standard disequilibrium retrieval exhibits any clear trend over the equilibrium temperature, at least for the current observed data. On the other hand, a larger diffusion coefficient for higher temperature has been predicted from both theory and numerical GCM simulations because of the increasing speed of vertical winds (Komacek et al. 2019). This trend is also tentatively indicated by the recent analysis of theSpitzer transit depths of about fifty gas giants byBaxter et al. (2021). The higher observational precision achievable with theJames Webb Space Telescope (JWST;Gardner et al. 2006) and the Atmospheric Remote-sensing Infrared Exoplanet Large-survey (Ariel;Tinetti et al. 2018) will enable us to further explore such predictions.

Compared to$K_{zz}^{\text{Ch}}$ (blue circle points),$K_{zz}^{\text{Cl}}$ (purple triangle points) is retrieved with much smaller uncertainties due to the strong effect of the diffusion coefficient on the cloud structure (Gao & Benneke 2018;Ormel & Min 2019). Also, the comparison of$K_{zz}^{\text{Cl}}$ between the standard disequilibrium and equilibrium retrievals (purple and orange triangle points in the bottom panel) shows that they are quite similar.

Next,K_zz from the disequilibrium retrieval employing the condition$K_{zz}^{\text{Ch}}=K_{zz}^{\text{Cl}}$ (green square points in both panels) does not exhibit any clear trend over the equilibrium temperature either. In some cases, the 1σ uncertainty range ofK_zz from the same-K_zz retrieval lies outside the median values of$K_{zz}^{\text{Ch}}$ and$K_{zz}^{\text{Cl}}$ from the standard disequilibrium retrieval. This is probably because we imposed the strong condition$K_{zz}^{\text{Ch}}=K_{zz}^{\text{Cl}}$ with a constant value throughout the atmosphere even though the location of the quenching and clouds affecting the spectrum are usually different and several orders of magnitude difference ofK_zz within the atmosphere is expected from GCM simulations (Parmentier et al. 2013;Zhang & Showman 2018;Komacek et al. 2019).

Among our samples, HD 209458b has a markedly large retrieved$K_{zz}^{\text{Ch}}$ (blue circle point) from the standard disequilibrium retrieval. When the condition$K_{zz}^{\text{Ch}}=K_{zz}^{\text{Cl}}$ is imposed, the retrieved value ofK_zz (green square point) decreases, especially when compared to$K_{zz}^{\text{Ch}}$. However, as we mentioned in Sect.3.2, even with this smallK_zz, the retrieved spectrum still exhibits disequilibrium features. Thus, we propose that this planet is an ideal target for studying disequilibrium chemistry in exoplanet atmospheres.

$K_{zz}^{\text{Ch}}$ retrieved from theT_int-fixed disequilibrium retrieval (red circle points in the top panel) is tightly constrained compared to the standard disequilibrium retrieval for HD 189733b and HD 209458b, while this is not the case for HAT-P-1b, WASP-12b, and WASP-19b. There are several possible reasons for this. First, for HD 189733b, the value ofT_int we fixed with the formula ofThorngren et al. (2019) is different from the retrieved value from the standard equilibrium retrieval by more than 3σ (seeTable 7), and theT_int-fixed disequilibrium retrieval results in a different thermal profile with high temperatures in the upper atmosphere (seeFig. 2d). This requires a larger$K_{zz}^{\text{Ch}}$ to reproduce the observed amplitude of the 1.4 μm H₂O absorption feature by overcoming its thermal dissociation due to high temperatures (compare the purple thick dash-dotted line and purple thin dash-dotted line inFig. 2d). On the other hand, in the case of the standard disequilibrium retrieval, the retrieved small value of$K_{zz}^{\text{Ch}}$ indicates that disequilibrium chemistry is not significantly preferred over equilibrium chemistry. This is also seen in the almost consistent median abundance profiles between the standard disequilibrium and equilibrium retrievals (compare the thick solid lines and dashed lines in the right panel ofFig. 2d). Thus,$K_{zz}^{\text{Ch}}$ remains unconstrained, though it needs to be sufficiently small. Next, for HD 209458b, disequilibrium chemistry plays an important role in its atmosphere, unlike the cases of WASP-12b and WASP-19b (seeFigs. 2 andA.1). Since the quenching happens in the deep atmosphere where the intrinsic temperature has a big influence on the thermal structure, the value of$K_{zz}^{\text{Ch}}$ could be tightly constrained whenT_int is fixed to a certain value for the atmospheres with disequilibrium chemistry. Given the somewhat stronger constraint on$K_{zz}^{\text{Ch}}$ for HD 209458b, we ideally expect a correlation betweenT_int and$K_{zz}^{\text{Ch}}$ in the results of the standard disequilibrium retrieval for this planet, which we do not find, as mentioned in Sect.3.3. We considerthat the complexity of the retrieval we performed, namely retrieval with more than ten retrieval parameters, could cause this apparent absence of the correlation. Moreover, we speculate that the precision of the observational data is also related because if the observational uncertainty is large, it would not have an impact on the constraint on the parameters even if we fix any of the parameters. The observational data for HD 189733b and HD 209458b are relatively precise when compared to those of the other three planets. Indeed, we found that when all the observational errors are artificially reduced to half their values,$K_{zz}^{\text{Ch}}$ becomes somewhat tightly constrained for HAT-P-1b while not for WASP-12b and WASP-19b whenT_int is fixed. In theatmosphere of HAT-P-1b, disequilibrium chemistry has a modest impact (seeFig. A.1h). Future investigations are needed to draw firm conclusions about the points raised.

Finally,$K_{zz}^{\text{Cl}}$ retrieved from theT_int-fixed disequilibrium retrieval (brown triangle points in the bottom panel) are largely similar to those of the standard disequilibrium retrieval except for HD 189733b, for which the retrieved thermal structure differs significantly (seeFig. 2d).

Fig. 3 Top: Eddy diffusion coefficients for the chemical structure simulations,$K_{zz}^{\text{Ch}}$, retrieved from the standard disequilibrium retrieval (blue circle points) and that with fixedTint (red circle points; only shown for the planets with masses larger than half Jupiter mass). Also shown are the eddy diffusion coefficients used for both chemistry and cloud simulations in the same-Kzzretrieval (green square points).Bottom: Eddy diffusion coefficients for the cloud structure simulations,$K_{zz}^{\text{Cl}}$, retrieved from the standard disequilibrium retrieval (purple triangle points), equilibrium retrieval (orange triangle points), and disequilibrium retrieval with fixedTint (brown triangle points; again, only shown for the planets with masses larger than half Jupiter mass). The eddy diffusion coefficients used for both chemistry and cloud simulations in the same-Kzz retrieval are also plotted (green square points, as in the top panel). The unfilled points indicate our hot Jupiter samples, while the filled points are the Neptune samples, which we define as masses smaller than 0.1 times the Jupiter mass. The error bars indicate 1σ uncertainty from the retrievals.

Fig. 4 Atmospheric metallicity (a) and C/O ratio (b) derived from the standard disequilibrium retrieval (blue circle points), equilibrium retrieval (orange triangle points), disequilibrium retrieval with sameKzz (green square points), andTint-fixed disequilibrium retrieval (red diamond points; only shown for the planets with masses larger than half Jupiter mass). The solar metallicity and C/O ratio from the solar photospheric elemental abundance ratios ofAsplund et al. (2021) are marked with dotted lines for reference.

3.4.2 Metallicity and C/O ratio

Figure 4 shows the derived values of (a) metallicity and (b) C/O ratio from our four different retrievals as a function of the planet mass. For most planets, the values of the two parameters derived from the different retrievals are consistent within the 1σ error. As for the metallicity trend, our results generally agree with the previous studies that found a tentative trend of higher metallicity with decreasing planet mass, such as seen for the solar system planets (Wakeford et al. 2017;Nikolov et al. 2018;Welbanks et al. 2019). On the other hand, we see no clear trend for the C/O ratio, at least for the current observational precision, which is also consistent with previous work (Min et al. 2020). Since both metallicity and C/O ratio are regarded as important indicators of planet formation (Öberg et al. 2011;Eistrup et al. 2016,2018;Mordasini et al. 2016;Notsu et al. 2020), future higher-precision observations are essential to further explore the trends of these parameters toward understanding the origin and formation mechanism of planets.

4 Discussion

4.1 Eddy diffusion coefficient

4.1.1 Range used in retrieval

As mentioned in Sect.2.3, for some of the cooler planets of our samples, the lower bound of the$K_{zz}^{\text{Ch}}$ range adopted in the retrievals, namely 10⁰ cm² s⁻¹, is not sufficiently small for all the species to be in chemical equilibrium. To confirm our findings, we performed additional standard disequilibrium retrievals extending the lower bound of$K_{zz}^{\text{Ch}}$ to an extremely small value of 10⁻⁴⁰ cm² s⁻¹. We chose this value so that the atmosphere lower than the 0.1 mbar level, which is the lower bound of the valid pressure range for the chemical timescale formula ofTsai et al. (2018), is in chemical equilibrium when calculated with the best-fit parameter set of the standard disequilibrium retrieval for all of our samples except WASP-39b. We note that while we found that the best-fit parameter set for WASP-39b requires$K_{zz}^{\text{Ch}}\lesssim {10}^{-60}$ cm² s⁻¹, the resultant temperature in its upper atmosphere is lower than the lower bound of the valid temperature range for the chemical timescale ofTsai et al. (2018), namely 500 K. Thus it is uncertain whether such a small value of$K_{zz}^{\text{Ch}}$ is needed toachieve the equilibrium condition in its atmosphere. Regarding our finding of the indication of disequilibrium chemistry for WASP-39b, we confirmed that the temperature in the relatively low atmosphere where the quenching happens is within the validated range, so this does not affect our conclusion.

Figure 5 shows the posterior probability distribution of$K_{zz}^{\text{Ch}}$ for (a) HD 209458b, (b) WASP-39b, and (c) HD 189733b. For HD 209458b and WASP-39b, for which we found indicative evidence of disequilibrium chemistry, the peak of the posterior distribution remains at a similar value regardless of the adopted wider range, ensuring the relatively large values of$K_{zz}^{\text{Ch}}$ we retrieved for these two planets. On the other hand, for all the other planets, as for the case of (c) HD 189733b, we find that adopting a smaller lower bound for the$K_{zz}^{\text{Ch}}$ range, the posterior distribution flattens and the median and 1σ confidence interval shift toward smaller values, inferring that the retrieval of$K_{zz}^{\text{Ch}}$ is affected by the choice of the$K_{zz}^{\text{Ch}}$ range and thus its retrieved value is not reliable. This is consistent with our results that the disequilibrium scenarios are never favored over the equilibrium ones for all of our samples except HD 209458b and WASP-39b in Sect.3.1.

Fig. 5 Posterior probability distribution of${\log }_{10}\left(K_{zz}^{\text{Ch}} ({\text{cm}}^2 {\text{s}}^{-1})\right)$from the standard disequilibrium retrievals with different ranges adopted, [0, 15] (blue), [−20, 15] (orange), and [−40, 15] (green). The vertical solid and dashed lines indicate the median value and the 1σ confidence interval, respectively.

4.1.2 Profile

While we used a constant eddy diffusion coefficient throughout the atmosphere in this study,K_zz is expected to vary over several orders of magnitude within the atmosphere (Parmentier et al. 2013;Zhang & Showman 2018;Komacek et al. 2019), which is also observed in the atmospheres of the solar system planets (e.g.,Allen et al. 1981;Moses et al. 2005). Employing an eddy diffusion coefficient profile, such as an increasing coefficient with decreasing pressure, will be the subject of future work.

4.2 Effect of photochemistry

As mentioned in the introduction, we ignored photochemistry. Its effect is more significant for cooler atmospheres because of the longer timescale of thermochemical reactions, and the atmospheres of the planets close to the host stars. The UV irradiation from the host star dissociates the molecules, such as H₂O, in the upper atmosphere. Thus, if we think naively, neglecting the photochemistry effect could lead to an underestimate of the chemical eddy diffusion coefficient needed to reproduce the observed absorption features of molecules, such as H₂O³. On the other hand, in hot (≳1000 K) hydrogen-dominated atmospheres, photodissociation of CO contributes to the formation of additional CO₂ in the upper atmosphere (e.g.,Moses et al. 2011;Venot et al. 2012) while in cool (≲ 1000 K) atmospheres, photodissociation of CH₄ forms both CO₂ and CO at high altitudes (e.g.,Miller-Ricci Kempton et al. 2012;Kawashima & Ikoma 2018). Thus, for the spectra with prominent CO₂ and/or CO features, neglecting the photochemistry effect could instead lead to overestimation of the chemical eddy diffusion coefficient. Among theHubble andSpitzer data we used in this study, CO₂ and CO can exhibit their features inSpitzer’s IRAC/Ch2 (4.5 μm) band.

To examine the effect of neglecting photochemistry on the observed indications of disequilibrium chemistry, we performed additional standard disequilibrium and equilibrium retrievals, excluding the data ofSpitzer’s IRAC/Ch2 (4.5 μm) band for WASP-39b, while the results of these retrievals for HD 209458b have already been presented inTable 4. We find that disequilibrium chemistry remains favored for both planets HD 209458b and WASP-39b, with natural logs of the Bayes factors of 6.9 and 2.3, respectively. These values are almost the same as those inTable 3, indicating that the features of CO₂ and CO hardly affect the need for disequilibrium chemistry in the atmospheres of those planets. The retrieved values of the chemical eddy diffusion coefficients are${\log }_{10}\left(K_{zz}^{\text{Ch}} ({\text{cm}}^2 {\text{s}}^{-1})\right)={13.05}_{-2.45}^{+1.37}$ and${8.65}_{-4.34}^{+3.76}$ for HD 209458b and WASP-39b, respectively, which are also similar to the retrieved values for the data including IRAC/Ch2 of${13.12}_{-4.61}^{+1.39}$ and${8.61}_{-4.22}^{+3.68}$. Thus, the retrieved large value of$K_{zz}^{\text{Ch}}$ for HD 209458b from the standard disequilibrium retrieval still holds and is not greatly affected by the effect of neglecting photochemistry.

4.3 Additional chemical species

In this study, we considered the quenching process for the molecules containing the most abundant elements next to hydrogen and helium, namely carbon, oxygen, and nitrogen. Currently, the reaction rate coefficients are less constrained for the molecules consisting of the other elements, while they are also expected to be subject to the quenching effect. Derivation of the chemical timescale for molecules, such as H₂S and PH₃, is urgently needed given their relatively good observability and large abundances.

4.4 Effect of photochemical haze

While we considered the effect of clouds in this study using the model ofOrmel & Min (2019), we ignored the effect of photochemical haze that is also expected to form in the atmospheres of exoplanets (Zahnle et al. 2009,2016;Miller-Ricci Kempton et al. 2012;Hu et al. 2013;Morley et al. 2013,2015;Gao et al. 2017,2020;Lavvas & Koskinen 2017;Kawashima & Ikoma 2018). Clouds and haze are expected to have opposite dependences on the eddy diffusivity. A largeK_zz can raise clouds,which are usually formed in the relatively low atmosphere, thus resulting in an optically thick atmosphere (Gao & Benneke 2018;Ormel & Min 2019). On the other hand, vigorous mixing efficiently removes photochemical haze from the upper atmosphere,making the atmosphere optically thinner (Kawashima & Ikoma 2019;Ohno & Kawashima 2020). Thus, the inclusion of haze particles, which we will leave for future work, might lead to different results, especially forK_zz.

4.5 Toward 2D and 3D retrieval modeling with disequilibrium chemistry

In this study, we considered 1D atmospheres for the sake of simplicity and reducing the computational cost, consistent with most current spectral retrieval models. Recently,Irwin et al. (2020) andFeng et al. (2020) explored an extension of the retrieval method beyond 1D modeling since close-in exoplanets are subject to tidal locking and thus significant horizontal and latitudinal variations in the atmospheric properties, which can be explored by phase curve observations, are expected. The approach we proposed in this study of including disequilibrium chemistry in spectral retrievals with a physical basis is also applicable to future 2D and 3D retrieval modeling, by adding horizontal or latitudinal transport term to the continuity-transport equation, Eq. (1) (e.g.,Cooper & Showman 2006;Mendonça et al. 2018).

5 Conclusions

In this study, we have implemented the disequilibrium effect of vertical mixing or quenching to the spectral retrieval code ARCiS (Min et al. 2020) with a physical basis. Adopting a chemical relaxation method with a chemical timescale derived byTsai et al. (2018), we have developed a module to compute the profiles of molecular abundances taking the disequilibrium effect into account for the major species in hydrogen and helium-dominated atmospheres, namely CH₄, CO, H₂O, NH₃, N₂, and CO₂. Then using ARCiS updated with this module, we have performed retrievals of the observed transmission spectra of 16 exoplanets with sizes ranging from Jupiter to mini-Neptune.

We have found indicative evidence of disequilibrium chemistry for HD 209458b and WASP-39b. For HD 209458b, the retrieved value of the eddy diffusion coefficient, which is used in the chemistry calculation, is as large as${\log }_{10}\left(K_{zz}^{\text{Ch}} ({\text{cm}}^2 {\text{s}}^{-1})\right)={13.12}_{-4.61}^{+1.39}$ from our standard disequilibrium retrieval, indicating that disequilibrium chemistry plays a significant role in determining the molecular abundance profiles in its atmosphere. Owing to the enhanced abundance of NH₃ due to the quenching effect, its retrieved spectrum exhibits a strong NH₃ absorption feature at 10.5 μm, which is absent in the retrieved spectrum from the equilibrium retrieval. This feature is accessible by JWST/MIRI. Thus, HD 209458b offers a unique opportunity to study disequilibrium chemistry in exoplanet atmospheres. Moreover, for HAT-P-11b and GJ 436b, we obtained relatively different results between the disequilibrium and equilibrium retrievals, such as a 2.9σ difference for the C/O ratio. This demonstrates the importance of taking disequilibrium chemistry into account for spectral retrieval, where we might otherwise misinterpret results. We have also examined the trend of the retrieved eddy diffusion coefficients over the equilibrium temperature, though no trend was found, possibly due to the limits of the current observational precision. This study makes clear that including a consideration of disequilibrium chemistry in spectral retrieval is essential in the coming era of JWST and Ariel.

Acknowledgements

We wish to thank S.-M. Tsai for kindly providing his calculation data for the code validation. The insightful advice and valuable comments on the manuscript received from K. Ohno are greatly appreciated. Also, we are grateful to P. Woitke, C. Visscher, and J. Moses for their helpful comments. Finally, we wish to thank the anonymous referee for his/her careful reading and constructive comments and the editor E. Lellouch for his helpful comments, both of which significantly helped improve this paper. The numerical computations were carried out on the PC cluster at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. This work was supported by JSPS KAKENHI Grant Numbers JP21K13984 and JP21J04998. Y.K. acknowledges support from the European Union’s Horizon 2020 Research and Innovation Programme under Grant Agreement 776403 and Special Postdoctoral Researcher Program at RIKEN. This work has made use of Numpy (Harris et al. 2020), Matplotlib (Hunter 2007), and corner.py (Foreman-Mackey 2016), and we are grateful to the developers of those packages.

Appendix A Atmospheric profiles

Fig. A.1 As forFig. 2 for the remaining planets.

Appendix B Corner plots

Fig. B.1 Corner plots for the standard disequilibrium retrieval (blue). Also shown in the posterior probability distributionplots are the results of the equilibrium retrieval (orange) and disequilibrium retrievals with the sameKzz (green) and with fixedTint (red; only shown for planets with masses larger than half Jupiter mass). The vertical solid and dashed lines in the posterior plots indicate the median value and the 1σ confidence interval for each parameter, respectively.Kzz retrieved from the same-Kzz retrieval is shown in the place of$K_{zz}^{\text{Ch}}$, and the vertical red dotted line in the posterior plot of log10Tint is the value assumed in theTint-fixed retrieval. The unit for each parameter is the same as that used inTable 2.

References

All Tables

All Figures

Fig. 1 Transmission spectra retrieved from the standard disequilibrium retrieval (blue), equilibrium retrieval (orange), and disequilibrium retrievals with the sameKzz (green) and with fixedTint (red; only shown for planets with masses larger than half Jupiter mass). The solid lines show the median models while the shaded regions indicate the 1, 2, and 3σ confidence intervals. The black points denote the observed data used.

In the text

	Fig. 1 continued. We note that for (m) WASP-12b, even thoughthe observed data ofKreidberg et al. (2015) are allowed to vertically shift in each retrieval simulation (see Sect.2.3 for details), their original values are plotted.
In the text

Fig. 2 Atmospheric profiles for four representative planets, while those for the other planets are presented inFig. A.1.Left: Pressure–temperature profiles retrieved from the standard disequilibrium retrieval (blue), equilibrium retrieval (orange), and disequilibrium retrievals with the sameKzz (green) and with fixedTint (red; only shown for planets with masses larger than half Jupiter mass). The solid lines show the median models while the shaded regions indicate the 1, 2, and 3σ confidence intervals.Right: Median abundance profiles of H2O, CO, and CH4 (left) and those of CO2, NH3, and N2 (right) from the standard disequilibrium retrieval (thick solid lines), equilibrium retrieval (thick dashed lines), and disequilibrium retrievals with the sameKzz (thick dotted lines) and with fixedTint (thick dash-dotted lines; only shown for planets with masses larger than half Jupiter mass). For reference, the median abundance profiles for the case where the equilibrium chemistry is imposed are shown with the corresponding thin lines.

In the text

Fig. 3 Top: Eddy diffusion coefficients for the chemical structure simulations,$K_{zz}^{\text{Ch}}$, retrieved from the standard disequilibrium retrieval (blue circle points) and that with fixedTint (red circle points; only shown for the planets with masses larger than half Jupiter mass). Also shown are the eddy diffusion coefficients used for both chemistry and cloud simulations in the same-Kzzretrieval (green square points).Bottom: Eddy diffusion coefficients for the cloud structure simulations,$K_{zz}^{\text{Cl}}$, retrieved from the standard disequilibrium retrieval (purple triangle points), equilibrium retrieval (orange triangle points), and disequilibrium retrieval with fixedTint (brown triangle points; again, only shown for the planets with masses larger than half Jupiter mass). The eddy diffusion coefficients used for both chemistry and cloud simulations in the same-Kzz retrieval are also plotted (green square points, as in the top panel). The unfilled points indicate our hot Jupiter samples, while the filled points are the Neptune samples, which we define as masses smaller than 0.1 times the Jupiter mass. The error bars indicate 1σ uncertainty from the retrievals.

In the text

Fig. 4 Atmospheric metallicity (a) and C/O ratio (b) derived from the standard disequilibrium retrieval (blue circle points), equilibrium retrieval (orange triangle points), disequilibrium retrieval with sameKzz (green square points), andTint-fixed disequilibrium retrieval (red diamond points; only shown for the planets with masses larger than half Jupiter mass). The solar metallicity and C/O ratio from the solar photospheric elemental abundance ratios ofAsplund et al. (2021) are marked with dotted lines for reference.

In the text

	Fig. 5 Posterior probability distribution of${\log }_{10}\left(K_{zz}^{\text{Ch}} ({\text{cm}}^2 {\text{s}}^{-1})\right)$from the standard disequilibrium retrievals with different ranges adopted, [0, 15] (blue), [−20, 15] (orange), and [−40, 15] (green). The vertical solid and dashed lines indicate the median value and the 1σ confidence interval, respectively.
In the text

	Fig. A.1 As forFig. 2 for the remaining planets.
In the text

Fig. B.1 Corner plots for the standard disequilibrium retrieval (blue). Also shown in the posterior probability distributionplots are the results of the equilibrium retrieval (orange) and disequilibrium retrievals with the sameKzz (green) and with fixedTint (red; only shown for planets with masses larger than half Jupiter mass). The vertical solid and dashed lines in the posterior plots indicate the median value and the 1σ confidence interval for each parameter, respectively.Kzz retrieved from the same-Kzz retrieval is shown in the place of$K_{zz}^{\text{Ch}}$, and the vertical red dotted line in the posterior plot of log10Tint is the value assumed in theTint-fixed retrieval. The unit for each parameter is the same as that used inTable 2.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.