Toward the Analysis of JWST Exoplanet Spectra: Identifying Troublesome Model Parameters

All the models consider continuum collisional induced absorption (CIA) for H₂–H₂ and H₂–He from the same references, atomic absorption coming from Na and K with various approaches for the resonant lines (see Section 5.1), and molecular absorption coming from H₂O, CH₄, CO, CO₂, NH₃, PH₃, TiO, and VO. A lot of work has been recently devoted to update molecular linelists, especially by the Exomol team (see exomol.com). The references for the molecular linelists used in the various models are given in Table 1. It is interesting to note that, except for CO₂, at least two of the three models use the same linelist for a given molecule, which makes it possible to infer the influence of using different linelists. Absorption coefficients are computed for a broadening due to:

• 1.  
  terrestrial atmosphere (air) for petitCODE
• 2.  
  a mix of H₂ and He for ATMOand Exo-REM.

Note PH₃ and CO₂ have been added in Exo-REM since the first version described in Baudino et al. (2015).

After considering the linelist sources, we focus on how each model computes the shape of the wings of each line (except for alkalies fully described in Section 5.1, because the different approaches for the wings treatment of the atomic lines used in the models induce significant differences in the model predictions). Our three codes use the Voigt profile, but with various "cut-off" (or sub-Lorentzian lineshape) implementations:

• 1.  
  For ATMO, the line cut-off is described in detail in Amundsen et al. (2014). The cut-off distance is calculated on-the-fly by estimating when the line mass absorption coefficient has reached a critical value.
• 2.  
  For petitCODE and Exo-REM, the applied sub-Lorentzian lineshape comes from Hartmann et al. (2002) for all molecules except CO₂, for which a profile from Burch et al. (1969) is used (see Mollière et al. 2017).

For fast computation, the opacity is treated by using opacity distribution functions and the correlated-k assumption. This technique is very useful, especially to manage the large number of lines in modern molecular linelists. In our case, all models use correlated-k coefficients computed molecule-by-molecule.

Our three models combine the molecules, assuming no correlation between species in any spectral interval. This is done using the method named "reblocking of the joint k-distribution" by Lacis & Oinas (1991) and more extensively described in Mollière et al. (2015) ("R1000" method in their Appendix B.2.1), with a few differences for ATMO as described in Amundsen et al. (2014, 2016).

The three models compute the equilibrium chemistry using various methods.

ATMO and petitCODE use a Gibbs energy-minimization scheme following the method of Gordon & McBride (1994).

ATMO uses the same thermochemical data as Venot et al. (2012), in the form of NASA polynomial coefficients (see McBride et al. 1993). The Gibbs minimizing method allows for depletion of gas phase species due to condensation and is included in the chemical solver.

petitCODE uses the same NASA polynomials, but extends them to low temperatures (<200 K) using the JANAF database (see also Appendix A2 of Mollière et al. 2017).

Exo-REM solves the thermochemical equilibrium layer-by-layer using on-line JANAF data from Chase (1998). Exo-REM does not use a general Gibbs energy minimization scheme, but rather a simplified approach. The molecules are grouped by independent subsets (see Table 2 in Baudino et al. 2015). In a given subset, the molecular mole fractions are found by solving a system of equations consisting of equilibrium constants and element conservation. When condensation occurs in a layer, the elements trapped in the condensed species are considered unavailable to the chemistry occurring at higher layers. In doing so, ExoREM includes the so-called cold trap phenomenon, in contrast to ATMO and petitCODE.

The starting point of a benchmark is the definition of a common set of parameters and of physical, chemical processes to be taken into account. Given that one of the key parameters determining the vertical structure of an atmosphere is its atmosphere composition, we must first define a common set of absorbers.

The common model considers the opacities from the most important molecules and atoms: NH₃, CH₄, CO, CO₂, H₂O, PH₃, Na, and K, as well as the CIA for H₂–H₂ and H₂–He. For the alkalies, we use the following H₂ broadening coefficients:

• 1.  
  Na: $0.37\times {(T/296)}^{-0.65}$ cm⁻¹ atm⁻¹ fitted to Allard et al. (2012),
• 2.  
  K: $0.40\times {(T/296)}^{-0.64}$ cm⁻¹ atm⁻¹ fitted to Allard et al. (2016).

For the far wings, we use a Voigt profile up to 4500 cm⁻¹ from line center for all Na and K lines, and zero absorption beyond. Other line shapes used in the literature are explored in Section 5.1. TiO and VO are not taken into account, in order to keep a minimal "common" model; we consider them only in the case of an extremely hot atmosphere. We have not considered the presence of clouds.

The abundance of the various absorbers depends on the chemical reactions at work in the atmosphere; we have considered a chemistry including the following reactant species: H₂, H, He, H₂O, H₂O(s), CO, CO₂, CH₄, NH₃, N₂, Na₂S(s), H₂S, Na, HCl, NaCl, K, KCl, KCl(s), NH₄Cl(s), SiO(s), Mg, Mg₂SiO₄(s), MgSiO₃(s), SiO₂(s), Fe, Fe(l), PH₃, H₃PO₄(l), P, P₂, PH₂, CH₃, SiH₄, and PO.

Unless specified otherwise, we consider an atmosphere with a solar metallicity. We use solar elemental abundances from Asplund et al. (2009).

By default, the radius, at 10 bars, of the exoplanet is fixed at 1.25 R_jup (where R_jup is the radius of Jupiter; we use 69,911 km for 1 R_jup) and the distance between the exoplanet and the observer is fixed at 10 pc.

After adapting (upgrading or downgrading) our models to this common model, we generated spectra and atmospheric structures (composition and temperature) in the following conditions.

First, we compute the abundance profiles at chemical equilibrium and spectra for five given temperature profiles. Indeed, our models are such that it is possible to impose a pressure–temperature profile (PT profile), and then run them without reaching the radiative-convective equilibrium, by just computing the chemistry and only doing the radiative transfer. We compare the result without iteration of the models in exactly the same conditions (in Section 3.3, we apply the iterations needed for the full convergence of a model to radiative-convective equilibrium).

We use the thermal profiles of Guillot (2010) (Figure 1), with surface gravity ${\mathrm{log}}_{10}(g[\mathrm{cgs}])=3.7$, solar metallicites, and effective temperatures (corresponding to the model of Guillot, these T_eff are not necessarily consistent with models used here) of: 500, 1000, 1500, 2000, and 2500 K. More specifically, we use Equation (29) of Guillot (2010), with f = 1/4, $\gamma =0.4$, $\kappa =0.01$ and ${T}_{\mathrm{int}}=200$ K. For T_eff = 1000 K, we use two additional values of the metallicity: 3 and 30 times the solar value.

Zoom In Zoom Out Reset image size

This allows us to explore how our results compare for a broad range of temperatures, as well as different enrichment values.

As a next step, we compare self-consistent atmospheric solutions of our codes, calculating radiative-convective equilibrium structures for known exoplanets. In the previous section, we fixed the same temperature-profile in each model, while in this section, we impose the T_eff itself. This means that we have to iterate on the P–T profile to get a self-consistent solution to the transfer equation. We do this for two direct imaging and two transiting exoplanets. We run our models considering published physical parameters for the exoplanets (based on http://exoplanet.eu).

For the two self-luminous planets, we stay in the benchmark condition, except for the metallicity. GJ 504 b (Kuzuhara et al. 2013) is a cold planet with T_eff = 510 K, $\mathrm{log}(g)=3.9$ and a metallicity z = 0.28 dex. VHS 1256–1257 b (Gauza et al. 2015) is an object hotter than GJ 504 b, with T_eff = 880 K, $\mathrm{log}(g)=4.24$ and a metallicity z = 0.21 dex. In the absence of strong constraints on the radius, we keep the benchmark value of 1.25 R_Jup.

For the two irradiated transiting exoplanets, we keep the solar metallicity but adapt the radius to the published values. For GJ 436 b (Butler et al. 2004), the host star parameters are an effective temperature T_eff of 3684 K and a radius of 0.464 solar. The planet has a T_eff = 712 K, a radius of 0.38 R_Jup, and a mass of 0.07 M_Jup. The incidence angle of the irradiation $\cos (\theta )=0.5$. The irradiation of the second transiting exoplanet considered here, WASP 12 b (Hebb et al. 2009), is more extreme and could possibly lead to a temperature inversion at high altitude. Thus, for this planet, we decided to add the chemistry and opacities for TiO and VO. The host star parameters in this case are an effective temperature T_eff of 6300 K and a radius of 1.599 solar. The planet has a high T_eff = 2536 K, a radius of 1.736 R_Jup, and a mass of 1.04 M_Jup. We assume an incidence angle for the irradiation of $\cos (\theta )=0.5$. For the Ti and V chemistry, we considered Ti, TiO, V, and VO.

In this section, we present the results of the comparison between our models following the benchmark protocol. We show here a selection of plots to illustrate these results. In Appendix A, we attached the complete set of spectra and profiles for all the various test cases. All the data are available online (as supplementary material).

For each model, we compare the spectra, the abundance profiles of the absorbers, and when relevant, i.e., for the four targets modeled in a self-consistent way, the temperature profiles. All spectra are plotted at the same spectral resolution (corresponding to that of Exo-REM, i.e., with a step and resolution of 20 cm⁻¹). The spectra are shown from near to mid-infrared, the wavelength range of the JWST.

At the end of our convergence process (fully described in Section 4.2), we reached a good agreement between all models in the common model conditions. This is observable in Figure 2, which shows the spectra of the cases with a T_eff = 1000 K at solar and super-solar metallicities. Figure 3 shows the spectra for the three hottest cases (T_eff = 1500 K, 2000 K, 2500 K), also with good agreement between ATMO, petitCODE, and Exo-REM. If we consider the general trend, the emission, transmission spectra, and molecular abundances calculated by the three radiative-convective equilibrium models are very similar. It is especially interesting to see how similar spectra of exoplanet targets are without a pre-defined temperature profile in the models (Figure 4, and in Appendix A Figure 25).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Nonetheless, we still observe some differences and will be focusing on these in the following paragraphs.

To highlight the impact of the remaining differences, we investigate the predicted error bars from JWST observations for VHS 1256–1257 b in Figure 4.

To obtain these error bars, we use the JWST Exposure Time Calculator (https://jwst.etc.stsci.edu/) to simulate a half-hour observation with the NIRSpec JWST instrument using the prism mode, and with the MIRI JWST instruments using the low-resolution spectroscopic mode.

While ATMO and Exo-REM codes provide similar results within the error bars, petitCODE provides significantly different results. In the case of Figure 4, its comes mainly from linelist differences, as present in the next section.

4.1.1. Differences Due to Opacities

4.1.2. Differences Due to Chemistry

While there is no large difference in the spectra, large differences in some molecule and atom abundances can be observed (See Figure 5), especially for PH₃ up to 0.1 bar for the case at 500 K (Figure 5(a), and in Appendix A Figures 26(a) and (b)), and for the alkalies in general at 1000 K (Figure 5(b) for 30 solar metallicity and Appendix A, Figure 22(b)–(d)). These differences occur each time the condensation curve of a compound is crossed, and Exo-REM gives systematically lower abundances values. This effect is explained by the implementation of the cold trap in Exo-REM. Indeed, Exo-REM makes the following assumption: if a molecule condensates at a given level in the atmosphere, the material is unavailable to the upper layers. For instance, as alkali condensation occurs deep in the atmosphere, Exo-REM predicts a lower amount of Na and K than do ATMO and petitCODE. These differences do not have a significant effect on spectra (Figure 2).

Zoom In Zoom Out Reset image size

To arrive at the level of similarity presented in the previous section, we had to make several changes in our models. In this section, we describe how we updated our models and the benchmark protocol.

Before the beginning of the benchmark, Exo-REM was updated by adding CO₂ to the absorbing molecules list and by updating the H₂–He CIA sources. Exo-REM also implemented a new way to combine correlated-k coefficients (described in Appendix B).

When we began this benchmark, we identified a major difference between the models concerning how the alkali far wings were accounted for (see Section 5.1). Without any good answer to this question (see Section 5.1), we decided to consider the same simple alkali treatment (same Voigt profiles with cut-off). The other major difference was identified as coming from the molecular far-wing lineshape. In fact, petitCODE did not initially include any cut-off in the line profile (unlike ATMO and Exo-REM). Across the full considered wavelength range, all lines had a Voigt profile (see Mollière et al. 2015, Appendix A, for the scheme). The first consequence of this comparison is that petitCODE now also includes a line cut-off as its default line opacity treatment, implemented in the same way as in Exo-REM. Moreover, we have studied the effect of considering a cut-off, when compared to the case without any cut-off application, in Section 5.2.

These modifications led to a significant improvement of the convergence. We then spotted a strong absorption in petitCODE results (between 4 and 5 μm), not visible in the results from other models. The corresponding absorber was PH₃, not added at this stage in ATMO or Exo-REM. This implied an important modification of ATMO and Exo-REM to add the chemistry and absorption of phosphine. This molecule has now been included in these models.

At this step, spectra obtained with the three models began to show good agreements, but we still observed differences that were identified as due to methane features and H₂–He CIA. These differences came from the use of different CH₄ line lists and continuum CIA coefficients. Indeed, at that time, petitCODE did not use the latest version of CH₄ linelist, but rather that of Rothman et al. (2013); it also was not using the latest CIA values, but rather those of Borysow et al. (2001) and Borysow (2002). ATMO used an updated version of H₂–He CIA (Richard et al. 2012), but not a complete one compared to Exo-REM (see Appendix B). Once each model used the CH₄ lines list, and the same H₂–He CIA as Exo-REM, the differences disappeared.

Three other effects have been identified as being important to converge on the results.

In the extreme case with a T_eff = 2500 K, petitCODE showed a strong effect due to ionization, a process not taken into account by the other models; hence, we decided to deactivate it for the benchmark protocol.

The next two effects have an impact on transmission spectra and were observable only at the end of the whole convergence process. They are taken into account to achieve the remarkably similar spectra that we obtained for the self-consistent transiting targets WASP 12 b and GJ 436 b (see Appendix A, Figure 25(c)).

First, the mean molecular weight (MMW) is known to have a strong effect on transmission spectra (Miller-Ricci et al. 2009; Croll et al. 2011). In the beginning of the benchmark study, the chemistry code written for and used for the first time in petitCODE (Mollière et al. 2017), erroneously included condensate species for calculating the MMW, for which the condensate mass fraction and condensate monomer mass was used. This mistake was spotted quickly during our benchmark process, and corrected before the publication of Mollière et al. (2017), and therefore never affected any published results or models of petitCODE.

The reason for the difference is that the monomers can be significantly more massive than the gas species molecules from which they are constructed, e.g., for silicates such as Mg₂SiO₄, leading to an overestimation of the MMW, which resulted in noticeable differences in the output transmission spectra between ATMO and petitCODE. In reality, condensible species only affect the MMW by depleting heavier molecules and atoms from the gas phase, hence leading to a decrease in MMW, rather than an increase. As the condensate particles in the atmosphere, which have masses much larger than the monomer mass, will usually not be collisionally coupled sufficiently well to the atmospheric gas, they should not be taken into account in the calculation of the MMW.

Second, transit spectroscopy is also sensitive to the pressure-radius profile used. It appears that we need to have exactly the same radius at a pressure of 10 bars to obtain similar transmission spectra without an offset.

At temperatures below the condensation of TiO/VO (<2000 K), the opacity in the visible part of the spectrum is dominated by the wings of the Na-D (5890 Å) and KI (7700 Å) resonance lines. While an accurate treatment of the line shape of usual atomic lines is only needed at ∼25 Å from the line center, the alkali resonant lines far-wings profile at thousands of Ångström from the line center still induce an important absorption effect. In the past literature from brown dwarfs atmosphere modeling, two groups have tried to address the issue of the shape of the alkali-resonant-line far wings:

• 1.  
  Burrows & Volobuyev (2003) have used standard quantum chemical codes and the unified Franck–Condon model in the quasi-static limit to calculate the interaction potentials and wing shapes for the alkali resonant lines in H₂/He atmospheres.
• 2.  
  Allard et al. (2007) have used a unified line shape semi-classical theory (Allard et al. 1999) using molecular potentials computed using a valence pseudo-potential (see Allard et al. 2007).

To our knowledge, these two works have never been benchmarked in the astrophysics literature. In this section, we illustrate the effects of the different alkali treatments by using the ATMO code on the Guillot–1500K equilibrium PT profile, and on the radiative/convective solution for GJ 504 b. We have used three treatments:

• 1.  
  The Voigt profiles, as defined in Section 3.1.
• 2.  
  The "Burrows" profiles, as implemented by Baudino et al. (2015).
• 3.  
  The "Allard" profiles, as used in Tremblin et al. (2015).

5.1.1. Guillot Profile T_eff = 1500 K

We first perform a benchmark study on a fixed PT profile, the Guillot T_eff = 1500 K profile used in Section 2.1. Figures 6 and 7 show the emission and transmission spectra. Both spectra highlight a low absorption in the wings of the Voigt profile; both Burrows and Allard profiles increase the absorption in the far wings. However, this increase is significantly bigger with the Allard profile, and the use of one or the other is therefore not a matter of indifference. For a Jupiter-like exoplanet transiting a solar-type star, the difference between the transmission spectrum using Allard profiles or Burrow profiles is found to be in the 20 ppm range, which is accessible to observations with the JWST. However, these differences should be attenuated when using models with clouds, so it might be difficult to disentangle which profile fits better the observations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.1.2. GJ 504b

We have used the different alkali treatments in the radiative/convective equilibrium model for GJ 504 b. Figure 8 shows the effects on the converged PT profile: the Allard profile leads to a slightly warmer deep atmosphere than the Voigt profile, and the Burrows profile to a significantly cooler deep atmosphere, with a difference of 100–200 K at 100 bars. These effects are directly linked to the differences in absorption caused by the line profile. Figure 9 shows the effects of the different alkali treatments on the emission spectrum of GJ 504 b.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Allard profile leads to a higher absorption both below 0.9 μm and above 1 μm, compared to a simple Voigt profile. The Burrows profile leads to a higher absorption, compared to Voigt profile, below 0.9 μm—but a lower one above 1 μm. In fact, with the Burrows profile, there is no absorption caused by the potassium far wing in the Y band, while a strong absorption appears with the Allard profile.

Comparisons with observations, as done in Tremblin et al. (2015) for T dwarfs, would favor the Burrows profile for this kind of object. The authors have used the Allard profile, but they had to assume a condensation of potassium in KAlSi₃O₈ to remove all effects of the far wing on the Y band. We emphasize again that this conclusion depends on the model used to produce the NIR reddening (temperature gradient reduction in the case of Tremblin et al. 2015). In the other way, a comparison of the profile effects on the spectra of L dwarfs seems to favor the Allard profile (private communication B. Burmingham). The limits we have faced in this alkali-treatment benchmark strongly advocate for a community effort to gain a better understanding on the far-wing behavior as a function of pressure and temperature. It is also of great importance to place better constraints on the different cloud models or alternatives.

In this section, we estimate the effect of applying a line "cut-off" (a sub-Lorentzian far wing profile) in the Voigt profile of the molecular lines, as described in Section 3. For that purpose, we hereafter compare to a Voigt profile without cut-off. Doing so, we face a difficulty with regard to CH₄. Indeed, the huge number of CH₄ lines made the calculation of CH₄ opacities without line cut-off infeasible (for this paper). To deal with that difficulty, we have estimated the effect of having no cut-off in two ways: (1) by considering the effect only on a limited wavelength range and with a fixed P–T profile (Section 5.2.1); and (2) by taking into account having no cut-off on all molecules except CH₄, whose opacity is not considered, and doing a self-consistent calculation.

5.2.1. Comparison at Fixed Temperature Structure

As a first step, we compare the results from a self-consistent atmospheric calculation for GJ 504 b, as obtained by petitCODE, including the sub-Lorentzian lineshape (cut-off), to those with full Voigt profiles, using the same temperature structure. We use the same set of absorbing species as in the baseline case. Due to the extensive line list of CH₄ (Yurchenko & Tennyson 2014), we calculated the opacity for this molecule in the full Voigt case only within 1.1–1.4 μm, and neglected all lines outside of this spectral region. Therefore, we will concentrate on this spectral region for our comparison. This spectral interval was chosen because both water and methane have an opacity minimum ranging from 1.2 to 1.3 μm, and a change in the line-wing continuum should be most noticeable in such regions.

The resulting emission spectra can be seen in Figure 10. The difference between the two cases is strongest in the peak of emission, i.e., the location where the total line opacity was the smallest and the contribution of the line-wing continuum the largest. We find maximum differences in the flux on the order of ∼20%. If we had included methane lines outside of the 1.1 to 1.4 μm region, the decrease in emission may have been even stronger, with lines further away potentially still non-negligibly contributing to the line-wing continuum.

Zoom In Zoom Out Reset image size

5.2.2. Comparison for Self-consistent Structures

The spectral test we carried out in the previous section was done only within a small spectral window, due to the extensive CH₄ line list, which makes the calculation of opacities applying a Lorentzian lineshape over the full spectral range infeasible and did not allow for comparing self-consistent structures and spectra arising from the Lorentzian versus sub-Lorentzian lineshape calculations. In this section, we therefore vary our baseline model and calculate hot Jupiter atmospheric structures without considering CH₄ opacities, because for all other species an opacity calculation with full Lorentzian lineshape is possible. We do not claim that an exclusion of CH₄ is realistic, although CH₄ should not be dominant at the elevated temperatures studied here; we merely carry out this test to study the effect of a line wing cut-off. The hot Jupiter is defined as an irradiated planet around a Sun-like star, with ${T}_{\mathrm{eff}}=1500$ K, ${M}_{\mathrm{Pl}}=1.2\ {M}_{\mathrm{Jup}}$, ${R}_{\mathrm{Pl}}=1\ {R}_{\mathrm{Jup}}$ (at 10 bar), ${T}_{\mathrm{int}}=100$ K, and an atmospheric enrichment of three times solar.

The P-T structures, emission, and transmission spectra of the self-consistent hot Jupiter cases with sub-Lorentzian and Lorentzian far wings, neglecting the methane opacity, are shown in Figures 11–13.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It can immediately be seen that the deep isothermal region of the full Voigt atmosphere is cooler, although the total opacity there is higher: the opacity increase arising from pure Voigt profiles causes the stellar light to be absorbed already at higher layers, leading to a small temperature increase between 0.1 and 3 bar in this atmosphere. Consequently, less flux has to be emitted from the optically thick, deeper regions of the planet, leading to a cooler temperature structure for pressures larger than 3 bar. Thus, in the presence of radiation, the calculated temperature structure of an atmosphere can change significantly, depending on whether or not a Lorentzian lineshape is applied.

Because the higher regions of the atmosphere are not so strongly affected by the line far wings, the flux changes in the absorption features of the emission spectra are small. The emission peaks with radiation originating in the deeper and hotter regions of the atmosphere are affected more strongly, with maximum differences of 7%, e.g., at 2.25 μm.

Over a wide wavelength range, the differences between the two cases in the transmission spectrum are small, due to the fact that the transmission probes high layers of the atmosphere where (i) the temperature of the two cases is very similar and (ii) pressure broadening is less important. However, the largest differences in the optical and NIR, again at locations of minimum opacity, can be as large as ∼10% of the transmission signal amplitude (e.g., Y-band in Figure 13). Such differences are greater than the expected precision of JWST observations of hot Jupiter around bright stars.

Known since the 1970s in the atmospheres of the giant planets of our solar system (Prinn & Barshay 1977), non-equilibrium chemistry was discovered in a cool brown dwarf, Gl 229 b, in 1997 by Noll et al. (1997). More recently, it started to be considered for transiting exoplanets (e.g., Moses et al. 2011; Line & Yung 2013), and evidence for non-equilibrium chemistry was also found in extrasolar directly imaged planets of the HR 8799 system (e.g., Barman et al. 2011, 2015; Konopacky et al. 2013) through the concomitant detection of CO absorption features and the unexpectedly weak absorption by CH₄.

This important phenomenon, which can strongly affect the atmospheric chemical composition, is due to a competition between vertical mixing and chemical kinetics. When the dynamical timescale is shorter than the chemical timescale, we observe a freezing of the abundance of some species. The "quenched" species keep the same abundances as the deep hot layers, where chemical equilibrium is achieved. This phenomenon, called transport-induced quenching, affects particularly the CO/CH₄ and N₂/NH₃ ratios, which can be driven away from what is expected from equilibrium chemistry. For the atmospheres highly affected by quenching, interpretation of observational spectra can thus be quite difficult.

To model the non-equilibrium chemistry that takes place in warm exoplanet atmospheres, two different approaches can be used: a complete modeling of kinetics with a model accounting for atmospheric mixing (like in the Venot chemical model, Venot et al. 2012, hereafter V12), or a chemical equilibrium model coupled to a parametrization of the quenching level for the main species (like in Exo-REM).

Venot chemical model. The thermo-photochemical model developed by V12 is a full 1D time-dependent model. It includes kinetics, photodissociation, and vertical mixing (eddy diffusion and molecular diffusion). To determine the atmospheric composition of a planet, the thermal profile is divided in discrete layers (∼100) of thickness equal to a fixed fraction of the pressure scale height. The code solves the continuity equation (Equation (1)) for each species and for each atmospheric layer, until a steady state is reached

$$\begin{eqnarray}&&\displaystyle \frac{\partial {n}_{i}}{\partial t}={P}_{i}-{n}_{i}{L}_{i}-\mathrm{div}({{\rm{\Phi }}}_{i}{{\boldsymbol{e}}}_{z})\end{eqnarray} \tag{ 1 }$$

where n_i the number density of the species i (${\mathrm{cm}}^{-3}$), P_i its production rate (${\mathrm{cm}}^{-3}\,{{\rm{s}}}^{-1}$), L_i its loss rate (${{\rm{s}}}^{-1}$), and ${{\rm{\Phi }}}_{i}$ its vertical flux (${\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$) that follows the diffusion equation,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{i}=-{n}_{i}{D}_{i}\left[\displaystyle \frac{1}{{n}_{i}}\displaystyle \frac{\partial {n}_{i}}{\partial z}+\displaystyle \frac{1}{{H}_{i}}+\displaystyle \frac{1}{T}\displaystyle \frac{{dT}}{{dz}}\right]-{n}_{i}{K}_{{zz}}\left[\displaystyle \frac{1}{{y}_{i}}\displaystyle \frac{\partial {y}_{i}}{\partial z}\right],\end{eqnarray} \tag{ 2 }$$

where y_i is the mixing ratio, K_zz is the eddy diffusion coefficient (${\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$), D_i is the molecular diffusion coefficient (${\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$), and H_i the scale height of the species i.

At both upper and lower boundaries, we impose a zero flux for each species.

The strength of this model relies on the chemical scheme it uses: the C₀–C₂ scheme. It describes the kinetics of 105 species made of H, C, O, and N, which are linked by ∼2000 reactions. This scheme has been implemented in close collaboration with specialists of combustion and validated experimentally as a whole (not only each reaction individually). The completeness of this scheme (both the forward and the reverse directions of each reaction are considered) implies that thermochemical equilibrium is achieved kinetically—or, conversely, that out-of-equilibrium processes are considered. The temperatures found in warm exoplanet atmospheres (particularly where quenching happens) are within the very large range of validation of the scheme ([300–2500] K and [0.01–100] bar), leading to a high level of confidence in the results obtained. Because the model is time-dependent, it naturally covers the phenomenon of quenching.

It is important to note that, contrary to Exo-REM, the chemical network of V12 does not account for P-bearing compounds, alkalies, and silicates.

Exo-REM model. The way Exo-REM accounts for non-equilibrium chemistry of the CO–CO₂–CH₄ and N₂–NH₃ networks is based on the analysis performed by Zahnle & Marley (2014) (hereafter ZM14). These authors applied a 1D full kinetics model to a set of atmospheric models for self-luminous objects with T_eff between 500 and 1100 K, along with a range of g, metallicity and vertical diffusivity. The objective was to determine the quench conditions that yield the asymptotic abundances of various non-equilibrium species in the upper atmosphere and describe those with simple equations. These equations provide a so-called chemical time (t_chem) expressed as a function of temperature, pressure, and in some cases, metallicity (Equations (12)–(14) for the CO–CH₄ system, Equation (44) for CO–CO₂, Equation (32) for the N₂–NH₃ system in ZM14). Equating t_chem to a dynamical time constant t_mix defined as H²/K_zz, where H is the atmospheric scale height and K_zz the eddy mixing coefficient, determines the quench level for the considered chemical system. The asymptotic abundances at high atmospheric levels are then given by those at this quench level. In Exo-REM, the mole fractions are simply set to those at thermochemical equilibrium up to the quench level and kept constant above this level. Doing so, we do not accurately represent the transition region located around the quench level, in which the composition gradually freezes, but reproduce the full kinetics model calculations of ZM14 below and above the quench level. Note that, as the CO–CO₂ system is quenched at lower temperatures, and thus higher levels, than CO–CH₄, the sum of the CO and CO₂ mole fractions is kept constant between these two quench levels while the CO/CO₂ ratio is still determined by thermochemical equilibrium.

In addition to the chemical networks investigated by ZM14, Exo-REM now includes that of P-bearing compounds. In Jupiter and Saturn, phosphine (PH₃) is observed with mole fractions of a few ppmv (Ridgway et al. 1976; Larson et al. 1980), orders of magnitude larger than the chemical equilibrium values (Fegley & Lodders 1994). Investigating a reaction network involving H/P/O compounds, Wang et al. (2016) concluded that, at equilibrium, H₃PO₄ is the major P-bearing species at temperatures below 700 K in Jupiter and Saturn and identified the limiting step for PH₃ to H₃PO₄ conversion. From their kinetic model, they concluded that, for any reasonable value of K_zz, the quench level for this reaction is at $T\approx$ 900 K, well below that where condensation of H₃PO₄ is expected (≈700 K), such that PH₃ is still the dominant P-bearing species at observable levels. We checked that this conclusion also applies to young giant exoplanets, which have T_eff much higher than Jupiter and Saturn. For example, in GJ 504 b, the quench level is around 1000 K while the H₃PO₄ condensation level is ≈520 K. Thus, the non-equilibrium scheme of Exo-REM simply discards formation of liquid H₃PO₄ and removes this compound from the list of P-bearing species used to calculate the phosphine vertical profile.

We first compare the two approaches previously introduced and then (Section 6.4) show the general effects (independent of the model) of non-equilibrium chemistry.

In this part, we use the same elemental abundances defined in the benchmark protocol. To compare the results obtained with the two chemical models, we apply the following process. First, we generate a temperature profile with the radiative-convective equilibrium model Exo-REM. This input profile is then used both by Exo-REM and V12 to compute the abundance profiles. Finally, Exo-REM performs the radiative transfer and generates the corresponding synthetic spectra. Because V12 does not consider Na, K, and PH₃, which are nevertheless important in radiative transfer, we used for these three species the abundances computed by Exo-REM in all spectra. We do this work for the two directly imaged planets and for the irradiated planet with T_eff = 1000 K.

We begin this comparison with the thermochemical equilibrium composition. In all simulated cases at equilibrium, the abundances profiles obtained with the two codes are very similar. However, we observe some differences for H₂O, CO, CH₄ (about a factor of 2 for CO and CH₄) in the deep atmosphere due to the silicate condensation taken into account in Exo-REM but not in the V12 model. To overcome this difference, we decided to apply a modification of the initial elemental abundances used in V12, in all the cases of this section, corresponding to a sequestration of ∼20% of the amount of oxygen in silicates (see in Appendix A, Figures 5(a), 29(a), 30 are after the correction).

For the non-equilibrium cases, the abundance profiles show us that for the two self-luminous targets we have a good agreement. However, this is not the case for the atmospheric model corresponding to T_eff = 1000 K (Figure 14 and Appendix A, Figures 28(b)–(d)). In fact we observe a strong difference in CO and CO₂ treatment (we find 1000× more CO with V12), stemming from the fact that ZM14, designed for self-luminous planets, does not address the kinetic inhibition against oxidizing CH₄ (we find 18% less CH₄ with V12) to CO contrary to V12.

Other specific and limited differences arise from different physical assumptions in the two models. For GJ 504 b, up to an altitude corresponding to 8 mbar, there is a clear effect of water cold trap (see description of Exo-REM in Section 2) with Exo-REM, not seen in V12 where it is not used (Figure 16 and, in Appendix A, Figure 29). On the other side, at the top, all atmosphere at K_zz = 10⁷ cm² s⁻¹. The difference comes from molecular diffusion, not implemented in Exo-REM (this difference appears at pressures lower than shown in the figures). Condensation of NH₃ is also a difference; Exo-REM takes into account the possibility of NH₄Cl and NH₄SH formation, explaining the difference at 8 mbar for VHS 1256–1257 b and 40 mbar for GJ 504 b (V12 does not consider these species; Figures 14, 17 and Appendix A, Figures 29(b)–(d), 30(b)–(d)).

Zoom In Zoom Out Reset image size

While not obvious for GJ 504 b (Figure 16), for VHS1256–1257 b (Figure 15) we see a difference coming from the fact that the V12 approach exactly computes the quenching process (giving a realistic transition across the region of chemical quenching) when compared to Exo-REM where the quenching level is imposed at one precise location. We observe more water (10%) and less methane (50%) up to 1 bar in Exo-REM.

To conclude, for directly imaged planets, the two chemical approaches give similar results, but this comparison confirms that Exo-REM, which uses the approach of ZM14, is not adapted to study the chemical composition of irradiated planets.

Effect on emission spectra. In this part, we detail the effect of non-equilibrium chemistry on emission and transmission spectra.

Figures 15 and 18 present spectra obtained with the abundance profiles shown in Figures 14 and 17. We can observe that the two chemical approaches (V12 and Exo-REM) yield very similar spectra. It is difficult to differentiate them, and the differences are within the error bars expected for VHS 1256–1257 b introduced in Figure 4. High-altitude departures observed in the abundance profiles (see previous subsection) have no visible effect on the spectra because there is not enough material to absorb efficiently at pressures lower than 1 mbar.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We use Exo-REM in the full benchmark condition, except for the non-equilibrium chemistry, to compute GJ 504 b emission spectra, and the full petitCODE (no constrained to benchmark conditions), combined with V12 non-equilibrium chemistry (and also at equilibrium), for the "Guillot" 1000 K case in emission and transmission (Figure 19).

Important differences are observed between equilibrium and non-equilibrium chemistry at two locations of the spectrum: between 4 and 5 μm, and between 8 and 20 μm. The first location is dominated by absorption of PH₃, CO, and CO₂, in out-of-equilibrium chemistry conditions. In the mid-infrared, it seems easier to discriminate between out-of-equilibrium cases regarding K_zz. CH₄ and NH₃ have the greatest impact at this location, and the depth of the line of the latter near 10 μm offers a possible way of constraining K_zz (Figures 15, 18). PH₃ also impacts the wavelength range just beyond 10 μm, and at 15 μm, the signature of CO₂ appears with non-equilibrium chemistry (Figure 17).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The coronagraphic mode of JWST's MIRI instrument is designed to be sensitive to the depth of the NH₃ feature, especially between the bottom (10.65 μm) and the edge of the band at (11.4 μm) (Figures 15, 18). The depth of this feature seems to be a good way to investigate non-equilibrium chemistry.

Keeping in mind the fact that temperature profiles of GJ 504 b (Figure 17) and VHS 1256-1257 b (Figure 18), for each K_zz, are computed self-consistently, we see some differences in the spectral energy distribution. This occurs because the model computes the spectrum by keeping the same integrated spectrum (i.e., T_eff). The difference between non-equilibrium and equilibrium chemistry is more important for GJ 504 b; around 10 μm, we observe a difference of 50% between equilibrium and non-equilibrium cases, but the normalized spectra of the non-equilibrium cases do not show a lot of variation (Figure 17). It seems to be easier to determine the K_zz with the VHS 1256–1257 b spectrum (spectra around 10 μm exhibit more variation in flux and the depth of NH₃ line is more dependent on the considered K_zz, Figure 18).

In a low T_eff case (Figure 17), non-equilibrium chemistry has an important global effect weakly dependent on the eddy coefficient (the non-equilibrium effect is something like saturated). At higher T_eff (Figure 18), non-equilibrium chemistry has less effect but deep and hot molecules are effectively transported, proportionally to the K_zz, to the top of the atmosphere. For GJ 504 b (Figure 17), around 10 μm the three spectra at non-equilibrium just seem to be shifted without any strong shape modification. In VHS 1256–1257 b (Figure 18), the shape of the NH₃ features varies significantly with the K_zz, with a significant effect of PH₃ on the low-wavelength wing.

PH₃ is the only molecule included in Exo-REM and not in V12 with a profile impacted by non-equilibrium chemistry (it is not the case for the alkalis). In the following paragraph, we study the effect of PH₃ on the spectra.

We compute first the full model of GJ 504 b with Exo-REM, at equilibrium, with PH₃. We then use the obtained PT profile in this first case in Exo-REM to generate spectra at equilibrium or with K_zz = 10¹¹ cm² s⁻¹, with or without PH₃. Figure 20 shows all the corresponding spectra. We see the effect of phosphine around 4.5 and 9 μm. At thermochemical equilibrium, the right wing of the PH₃ feature is visible at 4.5 μm, but the signature of PH₃ at 9 μm is completely hidden by the ammonia absorption. In the case of non-equilibrium chemistry, both PH₃ features are visible, especially the left wing around 4.5 μm, but the red wing is hidden by CO and CO₂.

Zoom In Zoom Out Reset image size