<figure><img src="image_1.png" /></figure>

HAL
open science

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

Jean-Loup Baudino, Paul Mollière, Olivia Venot, Pascal Tremblin, Bruno

Bézard, Pierre-Olivier Lagage

► To cite this version:

Jean-Loup Baudino, Paul Mollière, Olivia Venot, Pascal Tremblin, Bruno Bézard, et al. Toward the Analysis of JWST Exoplanet Spectra: Identifying Troublesome Model Parameters. The Astrophysical Journal, 2017, 850 (2), 10.3847/1538-4357/aa95be. [hal-01834166](https://hal.science/hal-01834166v1)

HAL Id: hal-01834166

https://hal.science/hal-01834166v1

Submitted on 27 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

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

Jean-Loup Baudino¹,²,³✉, Paul Mollière⁴, Olivia Venot⁵,⁶, Pascal Tremblin⁷,⁸✉, Bruno Bézard⁹✉, and Pierre-Olivier Lagage¹,²

¹Irfu, CEA, Université Paris-Saclay, F-91191 Gif-sur Yvette, France; jean-loup.baudino@physics.ox.ac.uk, jean-loup.baudino@obspm.fr

²AIM, Université Paris Diderot, F-91191 Gif-sur-Yvette, France

³Department of Physics, University of Oxford, Oxford, UK

⁴Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany

⁵Laboratoire Interuniversitaire des Systèmes Atmosphériques (LISA), CNRS UMR 7583, Université Paris-Est Créteil, Université Paris Diderot, 61 avenue du Général de Gaulle, F-94010 Créteil Cedex, France

⁶Institut voor Sterrenkunde, Katholieke Universiteit Leuven, Leuven, Belgium

⁷Astrophysics Group, University of Exeter, Exeter, EX4 4QL, UK

⁸Maison de la Simulation, CEA, CNRS, Univ. Paris-Sud, UVSQ, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

⁹LESLIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, Univ. Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, F-92195 Meudon, France

Received 2017 May 12; revised 2017 October 5; accepted 2017 October 21; published 2017 November 28

## Abstract

Given the forthcoming launch of the *James Webb Space Telescope* (*JWST*), which will allow observing exoplanet atmospheres with unprecedented signal-to-noise ratio, spectral coverage, and spatial resolution, the uncertainties in the atmosphere modeling used to interpret the data need to be assessed. As the first step, we compare three independent 1D radiative-convective models: *ATMO*, *Exo-REM*, and *petitCODE*. We identify differences in physical and chemical processes that are taken into account thanks to a benchmark protocol we have developed. We study the impact of these differences on the analysis of observable spectra. We show the importance of selecting carefully relevant molecular linelists to compute the atmospheric opacity. Indeed, differences between spectra calculated with Hitran and ExoMol exceed the expected uncertainties of future *JWST* observations. We also show the limits of the precision of the models due to uncertainties on alkali and molecule lineshape, which induce spectral effects that are also larger than the expected *JWST* uncertainties. We compare two chemical models, *Exo-REM* and *Venot Chemical Code*, which do not lead to significant differences in the emission or transmission spectra. We discuss the observational consequences of using equilibrium or out-of-equilibrium chemistry and the major impact of phosphine, detectable with the *JWST*. Each of the models has benefited from the benchmarking activity and has been updated. The protocol developed in this paper and the online results can constitute a test case for other models.

**Key words:** planets and satellites: atmospheres – planets and satellites: gaseous planets – radiative transfer

**Supporting material:** tar.gz file

## 1. Introduction

Since the publication of the detection of an extrasolar planet orbiting the 51 Peg star by Mayor & Queloz (1995), the field of exoplanet detection has developed very rapidly; more than 3700 exoplanets have been detected so far (e.g., exoplanet.eu). The field is now shifting more and more from the detection of exoplanets to the characterization of known exoplanets, especially their atmosphere, thanks to spectroscopic observations in the visible and infrared. One of the prime facilities used so far to observe the atmosphere of transiting exoplanets has been the *Hubble Space Telescope* (see the review by Deming & Seager 2017, and references therein). Thanks to the use of the so-called scanning observing techniques with the WFC3 instrument, as well as sophisticated data reduction methods, variation of the star flux due to the transit of an exoplanet has been detected down to a few tens of ppm (Kreidberg et al. 2014; Tsiaras et al. 2016). In the near future, the *James Webb Space Telescope* (*JWST*, see jwst.stsci.edu), thanks to a large collecting area and a suite of state-of-the-art instruments covering a large wavelength range (0.6–28 µm), will be a key tool to study the atmosphere of exoplanets (e.g., Greene et al. 2016; Mollière et al. 2017). The *JWST* will be used to study not only transiting/eclipsing exoplanets, but also

exoplanets observed by direct imaging, thanks to its high angular resolution. For those exoplanets far enough from stars to allow spectroscopic observations, spectra of exoplanet emission with a signal-to-noise ratio greater than 100 are expected.

Atmospheric models are needed to quantitatively interpret the observations. Several models have been developed over the years (for example, see the reviews by Helling et al. 2008; Marley & Robinson 2015; Hubeny 2017, and the references therein). Given the high observational precision achieved nowadays with current facilities, and in the near future with the *JWST*, the question of comparing the uncertainties in the model predictions with the precision of the observations has to be considered. To tackle this question, we have, as a first step, only considered one-dimensional models. To investigate the uncertainties in the model’s predictions, our approach has been to start with a comparison of the results from three models developed independently: the *petitCODE* model (Mollière et al. 2015, 2017), the *ATMO* model (Tremblin et al. 2015), and the *Exo-REM* model (Baudino et al. 2015). To identify the differences between the models and their impact on the predicted spectra, we have adopted a benchmark approach to give a common basis to the models, allowing us to disentangle the impact of such differences.
<table><caption>**Table 1**
References of the Different Sources of Molecular Opacities for the Three Models Considered (ATMO, Exo-REM, petitCODE)</caption><thead><tr><th>Opacities</th><th colspan="3">Common references</th></tr></thead><tbody><tr><td>CH<sub>4</sub></td><td colspan="3">Yurchenko & Tennyson (2014)</td></tr><tr><td>TiO and VO</td><td colspan="3">Plez (1998; with update from private communication)</td></tr><tr><td>H<sub>2</sub>-H<sub>2</sub></td><td colspan="3">HITRAN (Richard et al. (2012))</td></tr><tr><td>H<sub>2</sub>-He</td><td colspan="3">Borysow et al. (2001) and Borysow (2002)</td></tr><tr><th>Opacities</th><th>ATMO</th><th>Exo-REM</th><th>petitCODE</th></tr><tr><td>H<sub>2</sub>O</td><td>Barber et al. (2006)</td><td>Rothman et al. (2010)</td><td>Rothman et al. (2010)</td></tr><tr><td>CO</td><td>Rothman et al. (2010)</td><td>Rothman et al. (2010)</td><td>Rothman et al. (2010); Kurucz (1993)</td></tr><tr><td>CO<sub>2</sub></td><td>Tashkun & Perevalov (2011)</td><td>Rothman et al. (2013)</td><td>Rothman et al. (2010)</td></tr><tr><td>NH<sub>3</sub></td><td>Yurchenko et al. (2011)</td><td>Yurchenko et al. (2011)</td><td>Rothman et al. (2013)</td></tr><tr><td>PH<sub>3</sub></td><td>Sousa-Silva et al. (2015)</td><td>Sousa-Silva et al. (2015)</td><td>Rothman et al. (2013)</td></tr></tbody></table>

In the Section 2 of the paper, the current versions of the codes are briefly discussed and their differences outlined. The next section describes the protocol developed to benchmark the models. The results of the benchmark are discussed in Section 4. In Section 5, we focus on evaluating uncertainties that are due to the ways the shape of resonant-line far wings of alkalies are treated as well as to the molecule far wing lineshape in use. In Section 6, we consider out-of-equilibrium chemistry (see Moses et al. 2011; Venot et al. 2012; Line & Yung 2013; Zahnle & Marley 2014; Hu et al. 2015) using two approaches coming from Exo-REM and Venot et al. (2012) and the observational consequences of this phenomenon. Then, in Section 7, we compare the deviation due to the different modeling approaches to the expected JWST uncertainties. Our conclusions and perspectives are drawn in Section 8.

## 2. Description of the Radiative-convective Equilibrium Models

Before focusing on the benchmark itself, we describe here the current state of the three radiative-convective equilibrium models that have been benchmarked. This is the state of the models, taking into account the evolution resulting from the benchmark.

These models compute the 1D vertical structure of the atmosphere of giant planets and generate emission and/or transmission spectra given a set of input parameters: the effective temperature $T_{eff}$, the surface gravity $log_{10}(g)$, the radius, and the elemental composition.

The atmosphere is discretized in a number of layers (50, 64, and 120 for ATMO, Exo-REM, and petitCODE, respectively). The models calculate the net energy flux as a function of pressure level using the radiative transfer equation, and solve it iteratively for radiative-convective equilibrium, i.e., conservation of the flux. The codes also incorporate a thermochemical model that calculates layer-by-layer the molecular mole fractions, given the elemental abundances and temperature profile. Heating sources, which can be internal (e.g., coming from the formation and evolution process) or external (e.g., coming from the star), are taken into account.

All the models take the internal heating due to the contraction of the planet into account:

1. ATMO and petitCODE consider additional heating by the star radiation; they are thus adapted to model close-in exoplanets and to calculate transit spectra.

2. Exo-REM is only suited to young, giant exoplanets for which stellar heating can be neglected with respect to the internal heat flux. This model is not able to compute transmission spectra.

One of the key drivers to determine the structure of an atmosphere is its opacity, which depends on its atomic and molecular composition, which, in turn, depends on the chemistry at work in the atmosphere.

### 2.1. Opacities

All the models consider continuum collisional induced absorption (CIA) for H<sub>2</sub>-H<sub>2</sub> and H<sub>2</sub>-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<sub>2</sub>O, CH<sub>4</sub>, CO, CO<sub>2</sub>, NH<sub>3</sub>, PH<sub>3</sub>, 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<sub>2</sub>, 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<sub>2</sub> and He for ATMO and Exo-REM.

**Note** PH<sub>3</sub> and CO<sub>2</sub> 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).

## 2.2. Equilibrium Chemistry

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*.

# 3. Benchmark Protocol

The two following sections (Sections 3 and 4) are focused on the benchmark itself. First, we define a “common” model to compare the results of our models (Section 3.1) and the test protocol used to obtain the results (Sections 3.2 and 3.3). We then compare the results of this benchmark in Section 4.

## 3.1. Definition of a “Common” Model

The starting point of a benchmark is the definition of a common set of parameters and of physical, chemical processes

<figure><img src="image_4.png" /><figcaption>Figure 1. Input pressure-temperature profiles used in the benchmark.</figcaption></figure>

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} \text{ cm}^{-1} \text{ atm}^{-1}$ fitted to Allard et al. (2012),

2. K: $0.40 \times (T/296)^{-0.64} \text{ cm}^{-1} \text{ atm}^{-1}$ fitted to Allard et al. (2016).

For the far wings, we use a Voigt profile up to $4500 \text{ cm}^{-1}$ 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<sub>jup</sub> (where R<sub>jup</sub> is the radius of Jupiter; we use 69,911 km for 1 R<sub>jup</sub>) 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.

## 3.2. Comparison with Fixed Temperature Profiles

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 $log_g(g[cs]) = 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_{int} = 200$ K. For $T_{eff} = 1000$ K, we use two additional values of the metallicity: 3 and 30 times the solar value.

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

### 3.3. Comparison of Self-consistent Calculations

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, $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, $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.

## 4. Results

### 4.1. Results Comparison

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 \text{ cm}^{-1}$). 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).

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, it comes mainly from linelist differences, as present in the next section.

#### 4.1.1. Differences Due to Opacities

A large part of the differences originates from the absorption of $\text{PH}_3$, which is more observable for low-temperature cases (4–5 µm; see Figure 4 and Appendix A Figures 21(a) and (b), 25(a) and (b)) and treated differently, depending on the code: *petitCODE* uses the linelist from HITRAN (Rothman & Gordon 2013) while *ATMO* and *Exo-REM* use data from Exomol (Yurchenko & Tennyson 2014).

There are also some differences between spectra (Figure 2) at ~6 µm and 10 µm (and, in Appendix A, Figures 25(a) and (b)), which are attributed to the use of different $\text{NH}_3$ linelists; *petitCODE* uses HITRAN, whereas the other codes use Exomol (Yurchenko et al. 2011). The use of only the main CO isotopologue by *petitCODE* is at the root of the differences observed between 4 and 5 µm in the high-temperature cases (Figure 4 and Appendix A, Figures 23(b) and (c), 25(c)).

The use of different linelists also affects the integral flux of the spectrum, i.e., the effective temperature (e.g., in Appendix A, the temperature profiles of the two self-luminous exoplanets, Figure 27, are a bit shifted toward lower temperatures for *petitCODE*, to compensate for the lack of absorption).

These results point out the importance of updating line lists in the models in the *JWST* era.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Emission (a) and transmission (b) spectra for models at 1000 K with solar metallicity (solid), and 1000 K with 30 × solar metallicity (dashed).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Emission spectra (bottom to top) for effective temperatures of 1500, 2000, and 2500 K.</figcaption></figure>

At 2000 K, and much more so at 2500 K (Figure 3 and Appendix A, Figure 23(c)), *Exo-REM* has more flux than the other models. The reason for this is the range of validity of this model. Not designed for irradiated objects, correlated-k coefficients of *Exo-REM* are not computed for high temperature/low pressure. A temperature higher than 1800 K, at high altitude (for a pressure < $10^{-5}$ bar), is out of the range of validity of *Exo-REM*.

Note that we do not find differences stemming from the fact that *petitCODE* uses air broadening, or from the correlated-k approaches of *ATMO*.

### 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_3$ 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).

## 4.2. Convergence Process: Step by Step

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_2$ to the absorbing molecules list and by updating the $H_2$–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,
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Emission spectra of VHS 1256–1257 b (top) and GJ 504 b (bottom). The black uncertainties correspond to the combination of simulated NIRSpec/Prism and MIRI/LRS noise level for VHS 1256–1257 b for 0.5 hr of integration.</figcaption></figure>

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_{\text{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.

## 5. Focus on Key Parameters Effects

During the benchmark, we realized that the uncertainties in two inputs of the models, i.e., the alkali line profiles, and the molecules far line-wing shape, lead to significant uncertainties in the results. In this section, we detail the effects of these uncertainties.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Abundance profiles of the defined molecules for the case with a $T_{\text{eff}} = 500$ K at solar metallicity (a) and 1000 K at 30× solar metallicity (b). The curves for the models are often superposed, except for alkalies, where Exo-REM is not.</figcaption></figure>

### 5.1. Alkali Line Profiles

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_{\text{eff}} = 1500$ K

We first perform a benchmark study on a fixed PT profile, the Guillot $T_{\text{eff}} = 1500$ K profile used in Section 2.1. Figures 6 and 7 show the emission and transmission spectra. Both spectra

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Effects of the different alkali treatments on the emission spectra for the Guillot profile with $T_{\text{eff}} = 1500$ K.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Effects of the different alkali treatments on the transmission spectra for the Guillot profile with $T_{\text{eff}} = 1500$ K.</figcaption></figure>

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.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Effects of the different alkali treatments on the PT profiles computed in radiative/convective equilibrium for GJ 504b.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Effects of the different alkali treatments on the emission spectrum of GJ 504b.</figcaption></figure>

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.

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_3O_8$ 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

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Emission spectra for GJ 504 b, calculated using the nominal temperature structure of the case including sub-Lorentzian far wings. We show the resulting spectrum of the nominal case (blue solid line), as well as the case with full Voigt profiles (red dashed line).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Pressure–temperature structures of the self-consistent hot Jupiter cases’ sub-Lorentzian and Lorentzian far wings, neglecting the methane opacity.</figcaption></figure>

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.

5.2. Molecular Far-wing Lineshape

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_4$. Indeed, the huge number of $CH_4$ lines made the calculation of $CH_4$ 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);
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Emission spectra of the self-consistent hot Jupiter cases' sub-Lorentzian and Lorentzian far wings, neglecting the methane opacity.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Transmission spectra of the self-consistent hot Jupiter cases' sub-Lorentzian and Lorentzian far wings, neglecting the methane opacity.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Abundance profiles for the case at 1000 K and $K_{zz} = 10^7$ cm$^2$ s$^{-1}$. Solid lines correspond to the simplified approach of Exo-REM while dotted lines correspond to the chemical kinetic model of V12.</figcaption></figure>

and (2) by taking into account having no cut-off on all molecules except CH$_4$, 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 petit-CODE, 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$_4$ (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

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Abundance profiles of the defined molecules for the case VHS 1256–1257 b and $K_{zz} = 10^{11}$ cm$^2$ s$^{-1}$. Solid lines correspond to the simplified approach of Exo-REM while dotted lines correspond to the chemical kinetic model of V12.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Abundance profiles of the defined molecules for the case GJ 504 b and $K_{zz} = 10^{11}$ cm$^2$ s$^{-1}$. Solid lines correspond to the simplified approach of Exo-REM while dotted lines correspond to the chemical kinetic model of V12.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Emission spectra of GJ 504 b considering an atmosphere at chemical equilibrium (dark blue) and at disequilibrium ($K_{zz} = 10^7$, $10^9$ and $10^{11}$ cm$^2$ s$^{-1}$, respectively cyan, green and red). Spectra have been calculated using the chemical composition determined by the kinetic model of V12 (dotted lines) and Exo-REM (solid lines).</figcaption></figure>

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.
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Emission spectra of VHS 1256–1257 b at chemical equilibrium (dark blue) or in non-equilibrium chemistry ($K_{zz} = 10^7, 10^9$, and $10^{11} \text{ cm}^2 \text{ s}^{-1}$, shown respectively in cyan, green, and red) with our two chemical models. Location of molecular features, other than H$_2$O, are indicated. Solid lines correspond to the simplified approach of Exo-REM while dotted lines correspond to the chemical kinetic model of V12.</figcaption></figure>

### 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$_4$ 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$_4$ opacities, because for all other species an opacity calculation with full Lorentzian lineshape is possible. We do not claim that an exclusion of CH$_4$ is realistic, although CH$_4$ 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_{\text{eff}} = 1500$ K, $M_{\text{Pl}} = 1.2 M_{\text{Jup}}$, $R_{\text{Pl}} = 1 R_{\text{Jup}}$ (at 10 bar), $T_{\text{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.

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.

## 6. Non-equilibrium Chemistry

### 6.1. Impact of Non-equilibrium Chemistry

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$_4$.

This important phenomenon, which can strongly affect the atmospheric chemical composition, is due to a competition between vertical mixing and chemical kinetics. When the
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Emission (up) and transmission (bottom) spectra of the test exoplanet with a T<sub>eff</sub> = 1000 K, computed by petitCODE using equilibrium (blue) or non-equilibrium chemistry with K<sub>zz</sub> = 10<sup>11</sup> cm<sup>2</sup> s<sup>-1</sup> (red) using V12.</figcaption></figure>

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<sub>4</sub> and N<sub>2</sub>/NH<sub>3</sub> 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.

## 6.2. Description of the Chemicals Models

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

$$ \frac{\partial n_i}{\partial t} = P_i - n_i L_i - \text{div}(\Phi_i e_z) \quad (1) $$

where n<sub>i</sub> the number density of the species i (cm<sup>-3</sup>), P<sub>i</sub> its production rate (cm<sup>-3</sup> s<sup>-1</sup>), L<sub>i</sub> its loss rate (s<sup>-1</sup>), and Φ<sub>i</sub> its vertical flux (cm<sup>-2</sup> s<sup>-1</sup>) that follows the diffusion equation,

$$ \Phi_i = -n_i D_i \left[ \frac{1}{n_i} \frac{\partial n_i}{\partial z} + \frac{1}{H_i} + \frac{1}{T} \frac{dT}{dz} \right] - n_i K_{zz} \left[ \frac{1}{y_i} \frac{\partial y_i}{\partial z} \right], \quad (2) $$

<figure><img src="image_23.png" /><figcaption>**Figure 20.** Spectra of GJ 504 b, computed by Exo-REM with and without PH<sub>3</sub>, assuming thermochemical equilibrium or in non-equilibrium conditions (K<sub>zz</sub> = 10<sup>11</sup> cm<sup>2</sup> s<sup>-1</sup>) with the same PT profile (case with PH<sub>3</sub> at equilibrium).</figcaption></figure>

where y<sub>i</sub> is the mixing ratio, K<sub>zz</sub> is the eddy diffusion coefficient (cm<sup>2</sup> s<sup>-1</sup>), D<sub>i</sub> is the molecular diffusion coefficient (cm<sup>2</sup> s<sup>-1</sup>), and H<sub>i</sub> 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<sub>0</sub>–C<sub>2</sub> 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
<figure><img src="image_24.png" /><figcaption>**Figure 21.** Emission and transmission spectra for the case $T_{\text{eff}} = 500$ K (respectively (a) and (b)) and the cases $T_{\text{eff}} = 1000$ K (respectively (c) and (d)) with variation of metallicity (solid = solar, small dash = 3× solar, large dash = 30× solar). Exo-REM is in blue, petitCODE in green, ATMO in red.</figcaption></figure>

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$_2$–CH$_4$ and N$_2$–NH$_3$ 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_{\text{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_{\text{chem}}$) expressed as a function of temperature, pressure, and in some cases, metallicity (Equations (12)–(14) for the CO–CH$_4$ system, Equation (44) for CO–CO$_2$, Equation (32) for the N$_2$–NH$_3$ system in ZM14). Equating $t_{\text{chem}}$ to a dynamical time constant $t_{\text{mix}}$ defined as $H^2/K_{\text{zz}}$, where H is the atmospheric scale height and $K_{\text{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
<figure><img src="image_25.png" /><figcaption>**Figure 22.** Abundance profiles for the case $T_{\text{eff}} = 500$ K (a) and the cases $T_{\text{eff}} = 1000$ K with variation of metallicity (b) solar, (c) 3× solar, (d) 30× solar.</figcaption></figure>

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<sub>zz</sub>, the quench level for this reaction is at T ≈ 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<sub>eff</sub> 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
<figure><img src="image_26.png" /><figcaption>**Figure 23.** Emission and transmission spectra for the case $T_{\text{eff}} = 1500$ K (respectively (a) and (b)) and the cases $T_{\text{eff}} = 2000$ K (plain) and 2500 K (dash) (respectively (c) and (d)). Exo-REM is in blue, petitCODE in green, ATMO in red.</figcaption></figure>

$H_3PO_4$ and removes this compound from the list of P-bearing species used to calculate the phosphine vertical profile.

## 6.3. *Comparison of Results from Different Models*

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_3$, 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_{\text{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_2O$, CO, $CH_4$ (about a factor of 2 for CO and $CH_4$) 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).
<figure><img src="image_27.png" /><figcaption>**Figure 24.** Abundance profiles for the cases $T_{\text{eff}} = 1500$ K (a), 2000 K (b), and 2500 K (c).</figcaption></figure>

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_{\text{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^7 \text{ cm}^2 \text{ s}^{-1}$. 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)).

While not obvious for GJ 504 b (Figure 16), for VHS 1256–1257 b (Figure 15) we see a difference coming from the fact that the V12 approach exactly computes the quenching
<figure><img src="image_28.png" /><figcaption>**Figure 25.** Emission spectra of the two directly imaged planets (a). Emission (b) and transmission (c) spectra of the two transiting planets. For GJ 436 b, the flux is multiplied by 900 and the radius by 4, to put the two cases in the same plot.</figcaption></figure>

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.

## 6.4. Intrinsic Effect of Non-equilibrium Chemistry

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<sub>zz</sub>. 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<sub>zz</sub> (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).

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<sub>zz</sub>, are computed self-consistently, we see some differences in the spectral energy distribution. This occurs
<figure><img src="image_29.png" /><figcaption>**Figure 26.** Abundance profiles of the potential targets: (a) GJ 504 b, (b) VHS 1256–1257 b, (c) GJ 436 b, and (d) WASP 12 b.</figcaption></figure>

because the model computes the spectrum by keeping the same integrated spectrum (i.e., $T_{\text{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_{\text{zz}}$ with the VHS 1256–1257 b spectrum (spectra around 10 µm exhibit more variation in flux and the depth of NH$_3$ line is more dependent on the considered $K_{\text{zz}}$, Figure 18).

In a low $T_{\text{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_{\text{eff}}$ (Figure 18), non-equilibrium chemistry has less effect but deep and hot molecules are effectively transported, proportionally to the $K_{\text{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$_3$ features varies significantly with the $K_{\text{zz}}$, with a significant effect of PH$_3$ on the low-wavelength wing.

PH$_3$ 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$_3$ on the spectra.

We compute first the full model of GJ 504 b with Exo-REM, at equilibrium, with PH$_3$. We then use the
<figure><img src="image_30.png" /><figcaption>**Figure 27.** Temperature profiles computed by our three models (Exo-REM is in blue, petitCODE in green, ATMO in red) of the potential targets.</figcaption></figure>

obtained PT profile in this first case in Exo-REM to generate spectra at equilibrium or with $K_{zz} = 10^{11} \text{ cm}^2 \text{ s}^{-1}$, with or without PH$_3$. 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$_3$ feature is visible at 4.5 µm, but the signature of PH$_3$ at 9 µm is completely hidden by the ammonia absorption. In the case of non-equilibrium chemistry, both PH$_3$ features are visible, especially the left wing around 4.5 µm, but the red wing is hidden by CO and CO$_2$.

## 7. Influence of the Differing Model Assumptions on the Analysis of JWST Observations

In this section, we quantify differences coming from the previously identified model differences and compare them to the expected direct imaging uncertainties of JWST observations. For this analysis, we use the signal-to-noise ratio (S/N) that we previously estimated for the NIRSpec and MIRI instrument modes, described by Figure 4 in Section 4.1 for VHS 1256–1257 b.

We compute the $\chi^2$ in the wavelength range where we observe an effect between the average spectrum (between all the approaches) and the other spectra taken separately. We perform, in parallel, a Kolmogorov–Smirnov test (KStest) as done in Mollière et al. (2017). If the $\chi^2$ or the $p_{KStest}$ is around 1, then the differences between the approaches correspond to the uncertainties of the JWST. If, on the other hand, $\chi^2 \gg 1$ or $p_{KStest} \ll 1$, then the differences between the approaches are larger than the JWST uncertainties.

First, for the different alkali line profiles, we use the spectra plotted in Figures 6 and 9, restrained to the wavelength range between 0.6 and 0.9 µm (i.e., in the part of the NIRSpec range that is mostly impacted by alkalies).

We find that, in the case of a directly imaged exoplanet, the average $\chi^2$, as defined as above, is always larger than 20 and ($p_{KStest}$ is lower than $10^{-30}$).

Therefore, the precision expected with JWST allows us to hope to be able to learn more about the profile of the alkalies, even if this wavelength range is also impacted by the clouds (we do not study that effect in this paper).

We focus then on the effect of the selection of a given molecular far-wing lineshape, using the spectra plotted in Figure 12. We combine here wavelength ranges from NIRSpec

and MIRI: 0.6–13.9 µm (the same as for VHS 1256–1257b previously used in this paper).

For a directly imaged planet, we find that using or not using a far wing cut-off yields a significant difference (average $\chi^2 > 20$ and $p_{KStest} < 10^{-9}$).

The effect is more important in spectral windows outside the main molecular features, acting like a continuum, and may induce uncertainties in the retrievals of some atmospheric parameters, such as those of clouds, larger than those due to the uncertainties on the data themselves.

For directly imaged planets, differences between equilibrium and out of equilibrium compositions ($K_{zz} = 10^{11} \text{ cm}^2 \text{ s}^{-1}$) are largely higher than the observational uncertainties (the average $\chi^2$ is over 15 and $p_{KStest} < 10^{-4}$). It should be possible to determine what is the chemical regime of the observed atmosphere.

Finally, we focus on PH$_3$ in two wavelength ranges, 4–6 (NIRSpec) and 9–10 (MIRI) µm, using Figure 20. This molecule has a strong effect on the observable spectra and it should be easy to discriminate whether the atmosphere contains PH$_3$ or not (with a minimal average $\chi^2$ of 40 and $p_{KStest} < 10^{-10}$).

## 8. Conclusions and Perspectives

The forthcoming JWST will open a new era in the characterization of exoplanet atmospheres. The observed spectra with high S/N (> 100 for exoplanets detected in direct imaging, such as VHS 1256–1257b) expected with the JWST instruments will require precise and well-tested models to interpret the results. We need to be able to understand differences between models. To do that, we define in this paper a benchmark protocol.

We compared and updated three radiative-convective equilibrium models (ATMO, Exo-REM, and petitCODE) and compared two chemical models (V12 and Exo-REM). We showed that the analysis of JWST observations will be sensitive to differences in the modeling of some key parameters introduced in the various models. This fact should be kept in mind when drawing conclusions.

Many iterations were needed to achieve similar results in our benchmark conditions. We have spotted the importance of PH$_3$ in atmospheres, and also the crucial need for complete molecular linelists. Exomol is often considered as the most up-to-date database currently available, and is indeed better adapted to hot atmospheres than Hitran (built for Earth-like temperatures). However, the current version of Exomol is not the ultimate linelist database. Improvements are still going on; for example, those of Rey et al. (2017) for CH$_4$. During the definition of the benchmark protocol, we noticed that the chemistry model can also imply important differences, if one takes into account ionic species at hot temperature, isotopes, etc.

The first major difference spotted early in the process was the different treatment applied for the far-wing lineshape of alkali lines. The definition of the alkali treatment is crucial. Using one approach or the other can shift the temperature by more than 100 K and change the flux in spectral region where we also expect cloud effects (reddening of the spectrum observed, for example, in brown dwarf observations by Tsiuji et al. (1996)). However, comparison to observations does not allow us to select the best approach.

We also investigated a second effect. We found that applying a sub-Lorentzian lineshape in the far wings of molecular species can change the resulting spectra both in transmission and emission by up to ~20%, but usually less. This is comparable to what has been
<figure><img src="image_31.png" /><figcaption>**Figure 28.** Abundance profiles for T<sub>eff</sub> = 1000 K at equilibrium chemistry (a) and for K<sub>zz</sub> = 10<sup>7</sup> (b), 10<sup>9</sup> (c), and 10<sup>11</sup> (d) cm<sup>2</sup> s<sup>-1</sup>.</figcaption></figure>

found by Grimm & Heng (2015). In the self-consistent temperature structures obtained for irradiated planets, an increased opacity (due to neglecting the cut-off) leads to an inverse greenhouse effect in the atmosphere, because less stellar radiation is absorbed in the deep atmosphere if full Voigt profiles are used. Uncertainties in the far-wing lineshape are also a limitation in analyzing thermal emission spectra from solar system planetary atmospheres. This is the case for relatively transparent spectral windows, such as the methane windows in Titan between 1.3 and 5 µm (e.g., Hirtzig et al. 2013) or near-infrared windows observed on Venus' night side between 1 and 2.3 µm (e.g., Fedorova et al. 2015). Sub-Lorentzian lineshapes beyond a few wavenumbers from line core are clearly needed to reproduce the observed spectra, but precise

constraints on the profiles cannot be retrieved independently of other atmosphere or surface parameters. For giant exoplanets, while sub-Lorentzian far-wing profiles are also likely to occur, no prescription can be made today due to the lack of laboratory measurements of H<sub>2</sub>-broadened line far wings at high temperatures.

To compare models carefully, we need to use the same alkali treatment and far-wing lineshape for other absorbers.

We have also compared various chemical approaches. While we see clear differences between non-equilibrium and equilibrium abundances, we did not spot significant differences in the spectra based on the different implementations used for the two processes. The major differences occur at high altitudes: cold-trap and molecular scattering, without strong impact on the spectrum.
<figure><img src="image_32.png" /><figcaption>Figure 29. Abundance profiles for GJ 504 b at equilibrium chemistry (a) and for $K_{zz} = 10^7$ (b), $10^9$ (c), and $10^{11}$ (d) cm$^2$ s$^{-1}$.</figcaption></figure>

Among all the models parameters discussed and analyzed in the paper, those with a major impact on the observable are:

1. Alkali wing lineshape

2. PH$_3$ opacity

3. sub-Lorentzian lineshape

4. non-equilibrium chemistry

We showed that these four modeling parameters will have an impact on spectra with differences larger than the JWST uncertainties, especially for directly imaged exoplanets. However, in the case of alkali lineshape or molecular-line far wing, the differences in the modeled spectra could be translated as an artificial offset, without an easy way to identify it.

To improve the reliability of atmospheric models, particular efforts need to be undertaken to enhance experimental and theoretical work at high-pressure and high-temperature conditions in terms of:

1. wing lineshape, especially for alkalies

2. linelist completeness

This work is urgently needed if we want to be able to interpret correctly the observations of the JWST.

We also identified some other effects coming from:

1. adding isotopes

2. adding ions

3. using up-to-date linelists
<figure><img src="image_33.png" /><figcaption>**Figure 30.** Abundance profiles for VHS 1256–1257 b at equilibrium chemistry (a) and for $K_{zz} = 10^7$ (b), $10^9$ (c), and $10^{11}$ (d) cm$^2$ s$^{-1}$.</figcaption></figure>

4. definition of the reference radius for transmission spectrum

5. computation of the MMW.

Finally, the differences spotted in the abundance profiles
without any significant spectral effect found in the analysis are:

1. cold trap

2. molecular diffusion

Spectra and profiles generated for the benchmark (i.e., the
data behind the Figures 21–27 of the Appendix A) are available
as supplementary material. We encourage the community to
compare them to their models in these benchmark conditions

and iterate with us to continue to improve our respective
models and identify the existing differences.

P.O.L. acknowledges support from the LabEx P2IO,
the French ANR contract 05-BLAN-NT09-573739. O.V.
acknowledges support from the KU Leuven IDO project
IDO/10/2013, the FWO Postdoctoral Fellowship program,
and the Centre National d'Etudes Spatiales (CNES). J.L.B.
acknowledges the support of the UK Science and Technology
Facilities Council. We thank the anonymous referee for the
fast report and for the useful comments.
# Appendix A
Complete Set of Figures

The following Figures 21–30 show all the radiative-convective equilibrium models in the benchmark condition (see Sections 3 and 4).

These figures complement those shown in Section 6.

# Appendix B
Updates of Exo-REM

Exo-REM has been updated in several aspects since the version described in Baudino et al. (2015):

1. Correlated k-coefficients of different molecules are now combined 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)

2. CO₂ and PH₃ have been added to the list of molecular absorbers, using linelists from Rothman et al. (2013) and Sousa-Silva et al. (2015) respectively. For CO₂ pressure-broadened halfwidths (γ), we used the air-broadened values multiplied by 1.34 (Burch et al. 1969) to obtain the broadening by H₂ and He. For PH₃, we used: for J ≤ 20, γ = 0.1009–0.00137 × J if J ≠ K and 0.9 time this result if J = K with a temperature exponent n = 0.643 – 0.00912 × J; for J > 20, we fixed γ = 0.073 and n = 0.46. Relevant references are Levy et al. (1993), Salem et al. (2004), and Levy et al. (1994).

3. CIA of H₂–H₂ and H₂–He is now modeled following HITRAN (Richard et al. 2012), complemented for H₂–H₂ beyond 10,000 cm⁻¹ by Borysow et al. (2001) and Borysow (2002), with a small rescaling to ensure continuity at 10,000 cm⁻¹.

# ORCID iDs

Jean-Loup Baudino https://orcid.org/0000-0003-4061-2514

Pascal Tremblin https://orcid.org/0000-0001-6172-3403

Bruno Bézard https://orcid.org/0000-0002-5433-5661

# References

Allard, N. F., Kielkopf, J. F., & Allard, F. 2007, EPJD, 44, 507

Allard, N. F., Kielkopf, J. F., Spiegelman, F., Tinetti, G., & Beaulieu, J. P. 2012, A&A, 543, A159

Allard, N. F., Royer, A., Kielkopf, J. F., & Feautrier, N. 1999, PhRvA, 60, 1021

Allard, N. F., Spiegelman, F., & Kielkopf, J. F. 2016, A&A, 589, A21

Amundsen, D., Tremblin, P., Manners, J., Baraffe, I., & Mayne, N. 2016, A&A, 598, A97

Amundsen, D. S., Baraffe, I., Tremblin, P., et al. 2014, A&A, 564, A59

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Barber, R. J., Tennyson, J., Harris, G. J., & Tolchenov, R. N. 2006, MNRAS, 368, 1087

Barman, T. S., Konopacky, Q. M., Macintosh, B., & Marois, C. 2015, ApJ, 804, 61

Barman, T. S., Macintosh, B., Konopacky, Q. M., & Marois, C. 2011, ApJ, 733, 65

Baudino, J.-L., Bézard, B., Boccaletti, A., et al. 2015, A&A, 582, A83

Borysow, A. 2002, A&A, 390, 779

Borysow, A., Jørgensen, U. G., & Fu, Y. 2001, JQSRT, 68, 235

Burch, D. E., Gryvnak, D. A., Patty, R. R., & Bartky, C. E. 1969, JOSA, 59, 267

Burch, D. E., Gryvnak, D. A., Patty, R. R., & Bartky, C. E. 1969, JOSA, 59, 267

Burrows, A., & Volobuyev, M. 2003, ApJ, 583, 985

Butler, R. P., Vogt, S. S., Marcy, G. W., et al. 2004, ApJ, 617, 580

Chase, M. W., Jr. 1998, NIST-JANAF Thermochemical Tables, Monograph No. 9 (4th ed.; Woodbury, NY: American Institute of Physics)

Croll, B., Albert, L., Jayawardhana, R., et al. 2011, ApJ, 736, 78

Deming, L. D., & Seager, S. 2017, JGRE, 122, 53

Fedorova, A., Bézard, B., Bertaux, J.-L., Korablev, O., & Wilson, C. 2015, P&SS, 113, 66

Fegley, B., Jr., & Lodders, K. 1994, Icar, 110, 117

Gauza, B., Béjar, V. J. S., Pérez-Garrido, A., et al. 2015, ApJ, 804, 96

Gordon, S., & McBride, B. J. 1994, NASA Tech. Rep., Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications II. Users Manual and Program Description. 2: Users Manual and Program Description, NASA-RP-1311

Greene, T. P., Line, M. R., Montero, C., et al. 2016, ApJ, 817, 17

Grimm, S. L., & Heng, K. 2015, ApJ, 808, 182

Guillot, T. 2010, A&A, 520, A27

Hartmann, J.-M., Boulet, C., Brodbeck, C., et al. 2002, JQSRT, 72, 117

Hebb, L., Collier-Cameron, A., Loeillet, B., et al. 2009, ApJ, 693, 1920

Helling, C., Ackerman, A., Allard, F., et al. 2008, MNRAS, 391, 1854

Hirtzig, M., Bézard, B., Lellouch, E., et al. 2013, Icar, 226, 470

Hu, R., Seager, S., & Yung, Y. L. 2015, ApJ, 807, 8

Hubeny, I. 2017, MNRAS, 469, 841

Konopacky, Q. M., Barman, T. S., Macintosh, B. A., & Marois, C. 2013, Sci, 339, 1398

Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, Natur, 505, 69

Kurucz, R. L. 1993, SYNTHETHE Spectrum Synthesis Programs and Line Data (Cambridge, MA: Smithsonian Astrophysical Observatory)

Kuzuhara, M., Tamura, M., Kudo, T., et al. 2013, ApJ, 774, 11

Lacis, A. A., & Oinas, V. 1991, JGR, 96, 9027

Larson, H. P., Fink, U., Smith, H. A., & Davis, D. S. 1980, ApJ, 240, 327

Levy, A., Lacome, N., & Tarrago, G. 1993, JMoSp, 157, 172

Levy, A., Lacome, N., & Tarrago, G. 1994, JMoSp, 166, 20

Line, M. R., & Yung, Y. L. 2013, ApJ, 779, 3

Marley, M. S., & Robinson, T. D. 2015, ARA&A, 53, 279

Mayor, M., & Queloz, D. 1995, Natur, 378, 355

McBride, B. J., Gordon, S., & Reno, M. A. 1993, NASA Tech. Rep., Thermodynamic Data for Fifty Reference Elements, NASA/TP-3287/ REV1

Miller-Ricci, E., Meyer, M. R., Seager, S., & Elkins-Tanton, L. 2009, ApJ, 704, 770

Mollière, P., van Boekel, R., Bouwman, J., et al. 2017, A&A, 600, A10

Mollière, P., van Boekel, R., Dullemond, C., Henning, T., & Mordasini, C. 2015, ApJ, 813, 47

Moses, J. I., Visscher, C., Fortney, J. J., et al. 2011, ApJ, 737, 15

Noll, K. S., Geballe, T. R., & Marley, M. S. 1997, ApJL, 489, L87

Plez, B. 1998, A&A, 337, 495

Prinn, R. G., & Barshay, S. 1977, Sci, 198, 1031

Rey, M., Nikitin, V. A., Rannou, P., et al. 2017, Icar, submitted

Richard, C., Gordon, I., Rothman, L., et al. 2012, JQSRT, 113, 1276

Ridgway, S. T., Wallace, L., & Smith, G. R. 1976, ApJ, 207, 1002

Rothman, L., Gordon, I., Barber, R., et al. 2010, JQSRT, 111, 2139

Rothman, L., Gordon, I., Babikov, Y., et al. 2013, JQSRT, 130, 4

Rothman, L. S., & Gordon, I. E. 2013, in AIP Conf. Proc. 1545, Eighth International Conference on Atomic and Molecular Data and their Applications: ICAMDATA-2012, ed. J. D. Gillaspy, W. L. Wiese, & Y. A. Podpaly (Melville, NY: AIP), 223

Salem, J., Bouanich, J.-P., Walrand, J., Aroui, H., & Blanquet, G. 2004, JMoSp, 228, 23

Sousa-Silva, C., Al-Refaie, A. F., Tennyson, J., & Yurchenko, S. N. 2015, MNRAS, 446, 2337

Tashkun, S., & Perevalov, V. 2011, JQSRT, 112, 1403

Tremblin, P., Amundsen, D. S., Mourier, P., et al. 2015, ApJL, 804, 17

Tsiaras, A., Rocchetto, M., Waldmann, I. P., et al. 2016, ApJ, 820, 99

Tsuji, T., Ohnaka, K., & Aoki, W. 1996, A&A, 305, L1

Venot, O., Hébrard, E., Agúndez, M., et al. 2012, A&A, 546, A43

Wang, D., Lunine, J. I., & Mousis, O. 2016, Icar, 276, 21

Yurchenko, S. N., Barber, R. J., & Tennyson, J. 2011, MNRAS, 413, 1828

Yurchenko, S. N., & Tennyson, J. 2014, MNRAS, 440, 1649

Zahnle, K. J., & Marley, M. S. 2014, ApJ, 797, 41