# Toward a multidimensional analysis of transmission spectroscopy

Tiziano Zingales, Aurélien Falco, William Pluriel, Jérémy Leconte

► To cite this version:

Tiziano Zingales, Aurélien Falco, William Pluriel, Jérémy Leconte. Toward a multidimensional analysis of transmission spectroscopy. Astronomy & Astrophysics - A&A, 2022, 667, pp.A13. 10.1051/0004-6361/202243492. hal-04308955

HAL Id: hal-04308955

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

Submitted on 2 Dec 2023

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Toward a multidimensional analysis of transmission spectroscopy

## III. Modeling 2D effects in retrievals with TauREx*

Tiziano Zingales<fnref n="1" />,<fnref n="2" />✉, Aurélien Falco<fnref n="1" />, William Pluriel<fnref n="1" />,<fnref n="3" />, and Jérémy Leconte<fnref n="1" />

Received 7 March 2022 / Accepted 26 July 2022

## ABSTRACT

New-generation spectrographs dedicated to the study of exoplanetary atmospheres require a high accuracy in the atmospheric models to better interpret the input spectra. Thanks to space missions such as James Webb Space Telescope (JWST), ARIEL, and Twinkle, the observed spectra will indeed cover a large wavelength range from visible to mid-infrared with an higher precision compared to the old-generation instrumentation, revealing complex features coming from different regions of the atmosphere. For hot and ultra hot Jupiters (HJs and UHJs), the main source of complexity in the spectra comes from thermal and chemical differences between the day and the night sides. In this context, 1D plane parallel retrieval models of atmospheres may not be suitable to extract the complexity of such spectra. In addition, Bayesian frameworks are computationally intensive and prevent us from using complete 3D self-consistent models to retrieve exoplanetary atmospheres, and they require us to use simplified models to converge at a set of atmospheric parameters. We thus propose the TauREx 2D retrieval code, which uses 2D atmospheric models as a good compromise between computational cost and model accuracy to better infer exoplanetary atmospheric characteristics for the hottest planets. TauREx 2D uses a 2D parametrization across the limb which computes the transmission spectrum from an exoplanetary atmosphere assuming azimuthal symmetry. It also includes a thermal dissociation model of various species. We demonstrate that, given an input observation, TauREx 2D mitigates the biases between the retrieved atmospheric parameters and the real atmospheric parameters. We also show that having prior knowledge of the link between local temperature and composition is instrumental in inferring the temperature structure of the atmosphere. Finally, we apply such a model on a synthetic spectrum computed from a global climate model (GCM) simulation of WASP-121b and show how parameter biases can be removed when using 2D forward models across the limb.

**Key words.** atmospheric effects – instrumentation: spectrographs – methods: data analysis – radiative transfer – methods: statistical – occultations

## 1. Introduction

Exoplanets have diverse and complex atmospheres, which often require very sophisticated methods to be interpreted correctly. In particular, hot Jupiters (HJs) and ultra-hot Jupiters (UHJs) appear to have a strong dichotomy between the inflated hot day side and the cold and shrunken night side (Showman & Guillot 2002; Showman et al. 2008; Menou & Rauscher 2009; Wordsworth et al. 2011; Heng et al. 2011; Charnay et al. 2015; Kataria et al. 2016; Drummond et al. 2016; Tan & Komacek 2019; Pluriel et al. 2020, 2022). For such planets, it is necessary to use more-than-1D models to, at least, reproduce the different chemistry and physical properties realistically.

Transit spectroscopy allows the observer to infer atmospheric features by looking at a limited volume around the terminator line (Brown 2001; Kreidberg 2018). The transmitted light carries the information of both the day and the night side of the planet of a small region around the terminator line (Fortney et al. 2010; Caldas et al. 2019; Wardenier et al. 2022). For cold planets, a 1D plane parallel geometry can be enough to

interpret low-resolution data since the atmosphere is more homogeneous (MacDonald et al. 2020). Nevertheless, in the case of HJs and UHJs, the use of a higher-dimensional model becomes crucial to mitigate the biases between the retrieved model and the real 3D observed atmosphere. The impact of 3D atmospheric structures cannot be neglected and will have an impact on the spectra, and thus on the retrieval analysis because of chemistry heterogeneities (Changeat et al. 2019; Baeyens et al. 2021), or due to morning-evening asymmetry (Line & Parmentier 2016; Parmentier et al. 2016; MacDonald et al. 2020; Espinoza & Jones 2021).

Atmospheric retrievals encompass a series of algorithms which try to associate an input exoplanetary spectrum to an atmospheric model (Rodgers 2000; Irwin et al. 2008; Madhusudhan & Seager 2009; Line et al. 2013; Waldmann et al. 2015b,a; Gandhi & Madhusudhan 2018; Al-Refaie et al. 2021; Mollière et al. 2019; Zhang et al. 2019; Min et al. 2020). Bayesian frameworks are very useful to find the best atmospheric models which will reproduce the input spectra. They are often the main analysis tool to infer atmospheric properties of exoplanets. The use of such tools poses two main challenges: computational cost and model accuracy. First, Bayesian analysis is typically computationally intensive, especially when one tries to fit several model

<footnote n="1">¹ Laboratoire d'Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France
e-mail: tiziano.zingales@unipd.it</footnote>

<footnote n="2">² Università di Padova, Dipartimento di Astronomia, vicolo dell'Osservatorio 3, 35122 Padova, Italy</footnote>

<footnote n="3">³ Observatoire astronomique de l'Université de Genève, chemin Pegasi 51, 1290 Versoix, Switzerland</footnote>

<footnote>* The TauREx 2D plug-in is available at https://forge.oasu.u-bordeaux.fr/falco/taurex_2d</footnote>
<table><caption>**Table 1.** Input parameters for the three atmosphere configurations.</caption><thead><tr><th>Test</th><th><i>T</i><sub>deep</sub></th><th><i>T</i><sub>up, day side</sub></th><th><i>T</i><sub>up, night side</sub></th><th><i>P</i><sub>iso</sub></th><th>log<sub>10</sub>[H<sub>2</sub>O]<sup>(*)</sup></th><th>log<sub>10</sub>[CO]</th><th>β</th></tr></thead><tbody><tr><td>1D fully isothermal</td><td>2500 K</td><td>2500 K</td><td>2500 K</td><td>10<sup>4</sup> Pa</td><td>-3.30</td><td>-3.36</td><td>—</td></tr><tr><td>Two layers</td><td>1800 K</td><td>2600 K</td><td>2600 K</td><td>10<sup>4</sup> Pa</td><td>-3.30</td><td>-3.36</td><td>—</td></tr><tr><td>Day/night</td><td>2500 K</td><td>3300 K</td><td>500 K</td><td>10<sup>4</sup> Pa</td><td>-3.30</td><td>-3.36</td><td>20°</td></tr></tbody></table>

**Notes.** (*)These values change as a function of the pressure and temperature when water dissociation is assumed.

of an isothermal annulus ($P_{iso}$). This parametrization, although simplistic, is based on GCM simulations (Showman et al. 2015; Parmentier et al. 2018; Tan & Komacek 2019) to capture the most salient features of the thermal structure of HJs.

The temperature is shown by the colors (yellow to purple) in Fig. 1. The data (pressure, temperature, and abundances) in the input grid ($\alpha^*, P$) are then interpolated in an altitude-based coordinate system ($\alpha^*, z$), shown by the red grid in Fig. 1. Using the new coordinate system ($\alpha^*, z$), rays of light crossing the atmosphere may be subdivided into a list of segments, as illustrated in Falco et al. (2022), using simple geometric formula. There is one ray per layer in the model. The length of a ray $r$, which crosses the $i$th atmospheric layer, is identified with $\Delta\ell_{r,i}$. The optical depth can then be computed for each ray $r$ with the following:

$$ \tau_r^\lambda = \sum_i \frac{P_{r,i}}{k_B T_{r,i}} \left( \sum_{m=1}^{N_{\text{gas}}} \chi_{m,r,i} \sigma_{m,\lambda} + \sum_{j=1}^{N_{\text{con}}} k_{\text{mie},j} \right) \Delta\ell_{r,i}, \quad (2) $$

where $P_{r,i}$ and $T_{r,i}$ are the pressure and temperatures of the cell corresponding to the segment $i$ of the ray $r$, $k_B$ is the Boltzmann constant, $\chi_{m,r,i}$ is the volume-mixing ratio of the $m$-th molecule, $\sigma_{m,\lambda}$ is the total cross section of Rayleigh scattering and molecular and continuum absorptions, and $k_{\text{mie},j}$ is the cross section associated with the Mie scattering for the $j$th aerosol. In addition, $N_{\text{gas}}$ is the number of chemical species and $N_{\text{con}}$ is the number of aerosols. We take collision-induced absorption for H$_2$–H$_2$ and H$_2$–He pairs from HITRAN into account (Richard et al. 2012). Conventionally, the contribution of collision-induced absorption (CIA) opacities is normalized to the number of H$_2$ molecules and included in $\sigma_{H_2,\lambda}$.

The volume-mixing ratios for all the species can also be differentiated with a day and a night value array. In this work, we have set the number of layers to 100 and the number of angular points (or slices) to 24. These two numbers are a good compromise between precision of the model and computational speed (Falco et al. 2022). Two angular points correspond to the full day and night sides, while the rest equally subdivide the space between $-\beta/2$ and $\beta/2$, accounting for the linear transition in between.

We tested TauREx 2D retrievals using three different configurations described in Table 1. In all three configurations, the thermal structure of the atmosphere is described using the following free parameters:

1. $\beta$ angle, the angle across the terminator line within which we computed the radiative transfer equation;

2. $P_{iso}$, the pressure of the boundary between the lower and the upper level;

3. $R_p$, the 10 bar radius of the planet;

4. $T_{\text{day/night/deep}}$, the temperature in the day and night side of the planet and in the deeper isothermal region, respectively; and

5. $P_{\text{clouds}}$, the pressure of the upper level of the cloud deck.

The TauREx 2D retrieval code is configured to have a deep atmospheric region where the molecular abundances and the temperature remain constant: at pressures higher than $P_{iso}$, the temperature is equal $T_{\text{deep}}$ and the volume-mixing ratio of each molecule has a constant value. In the upper part of the atmosphere, where $P < P_{iso}$, the temperature is different between the day side ($T_{\text{day}}$) and the night side ($T_{\text{night}}$).

For the chemistry, we can use two different schemes:

1. Constant (or “free”) chemistry mode, where volume-mixing ratios are constant within each column ($\alpha$) for $P < P_{iso}$. In this mode, TauREx 2D can fit the abundances on each hemisphere independently; this implies three free parameters per molecule ([mol]$_{side1}$, [mol]$_{side2}$, and [mol]$_{deep}$);

2. Thermal dissociation mode, where we assume local equilibrium chemistry. In this mode, the free parameters are the deep abundances of each molecule ([mol]), which are considered uniform over the planet, but we recomputed the local abundances as a function of the local pressure and temperature following Parmentier et al. (2018). In the regime we consider, these changes in abundances are mainly driven by thermal dissociation.

The thermal dissociation module we use in this work follows the parametrization of Parmentier et al. (2018), which uses a modified version of the NASA CEA Gibbs minimization code (Gordon et al. 1984) in a grid of models to investigate gas and condensate equilibrium chemistry in different atmospheric conditions in substellar objects (Moses et al. 2013; Kataria et al. 2015; Skemer et al. 2016; Wakeford et al. 2017; Burningham et al. 2017; Marley et al. 2017). While the constant chemistry mode could seem better as it offers more flexibility, we show that its higher number of parameters creates more degeneracies than the other parametrization (see Sect. 4) so that the thermal dissociation mode is our default except when stated otherwise. In the tests we show in the next sections, we assumed constant 30 ppm error bars in the considered wavelength range 0.5–10 µm. The assumed error was kept throughout the whole spectral range with a normal distribution (Greene et al. 2016). This value was used as a noise floor in Pluriel et al. (2022, 2020) and it is a good representation of the average noise computed using the estimated mean number of photoelectrons in Cowan et al. (2015) for a WASP-121b-like star at 270 pc during a single transit. This wavelength range will be covered using MIRI and NIRSpec instruments on board of the JWST space telescope (Gardner et al. 2006).

## 3. Test and validation

In all of the test cases, we assume the same planetary and stellar parameters as of the system WASP-121, whose values are specified in Table 2. TauREx 2D works within the same Bayesian framework as the 1D TauREx retrieval code (Al-Refaie et al. 2021). The atmospheric retrieval, which uses a 2D forward
<table><caption>**Table 2.** Input planetary and stellar parameters.</caption><thead><tr><th colspan="2">Planetary parameters</th></tr></thead><tbody><tr><td>R<sub>p</sub></td><td>1.807R<sub>Jupiter</sub></td></tr><tr><td>M<sub>p</sub></td><td>1.183M<sub>Jupiter</sub></td></tr><tr><td>g<sub>s</sub></td><td>9.39ms<sup>-2</sup></td></tr><tr><td>H<sub>e</sub>/H<sub>2</sub></td><td>0.25884514</td></tr><tr><td>[log P<sub>min</sub>, log P<sub>max</sub>]</td><td>[-4, 6](Pa)</td></tr></tbody><thead><tr><th colspan="2">Stellar parameters</th></tr></thead><tbody><tr><td>R<sub>*</sub></td><td>1.458 R<sub>☉</sub></td></tr><tr><td>T<sub>eff</sub></td><td>6460 K</td></tr></tbody><thead><tr><th colspan="2">Retrieval bounds</th></tr></thead><tbody><tr><td>β (degrees)</td><td>0–30</td></tr><tr><td>T (K)</td><td>100–3600</td></tr><tr><td>log P<sub>iso</sub> (Pa)</td><td>-4 to 6</td></tr><tr><td>R<sub>p</sub> (R<sub>jupiter</sub>)</td><td>0.9–2.7</td></tr><tr><td>log mixing ratios</td><td>-10 to -1</td></tr><tr><td>log P<sub>clouds</sub> (Pa)</td><td>-4 to 6</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Fig. 3.** Best-fit spectrum for the day and night convergence test.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Pressure temperature (**left**) and log [H₂O] mixing ratio (**right**) profiles for the three different test cases: FI, fully isothermal; TL, two layers; and DN, day-night cases. The DN case has been split into three curves showing the profiles at the terminator, the day side, and the night side.</figcaption></figure>

model, is more computationally intensive than the one calculated with a 1D forward model due to its additional number of parameters and geometrical dimensions.

Before focusing on a realistic test case computed using SPARC/MITgcm (Showman et al. 2009), it is fundamental to test the reliability of the TauREx 2D code on some synthetic spectra generated by itself. We decided to test the code using the three cases described in Table 1. For these tests, the input spectra have been generated by the TauREx 2D forward model itself. The three cases represent three different scenarios with different pressure-temperature (PT) profiles and then H₂O mixing ratio profiles as shown in Fig. 2. All the simulated spectra have been generated assuming a resolution R = λ₀/Δλ = 100, where λ₀ = 5.25µm and Δλ = [0.5µm–10µm].

The retrieval tests were performed using the MultiNest algorithm (Feroz et al. 2009) with 1000 live points and uniform prior distribution between the bounds shown in Table 2. In the next three subsections, we discuss the results for the three test configurations. For the sake of conciseness, the fitted spectrum and

posterior distributions are only shown for the third, most complex case. To summarize, the three tests show that TauREx 2D can perfectly retrieve itself when there is enough information in the input spectrum.

### 3.1. Fully isothermal atmosphere

A fully isothermal atmosphere represents the simplest assumption in our test set. We assumed a constant vertical as well as day and night composition with H₂O and CO defined in Table 1.

The solution found is consistent with the input model. TauREx 2D found a solution which is consistent with the input parameters. TauREx found a volume-mixing ratio for log [H₂O] = -3.30<sup>+0.01</sup><sub>-0.01</sub> for water and log [CO] = -3.36<sup>+0.05</sup><sub>-0.06</sub> for CO, centered around the ground truth for this case. Since we used a forward model with two vertical layers, divided at pressure P<sub>iso</sub>, TauREx found the best model where the upper layer is suppressed at the top of the atmosphere around 10<sup>-3.07</sup> Pa and treated the rest of the atmosphere as isothermal with T = 2502<sup>+10</sup><sub>-11</sub>, which is consistent with the input model.

The 2D model of TauREx is equivalent to the 1D version of TauREx in the case of isothermal atmospheres, as demonstrated in Falco et al. (2022). Therefore, this retrieval test is a confirmation of the consistency of the implementation with its 1D predecessor.

### 3.2. Two vertical layers atmosphere

The test case with two layers is an atmosphere where the lower layer has a temperature T<sub>deep</sub> = 1800 K and the upper temperature has T<sub>upper, day/night</sub> = 2600 K, but we do not have any day or night difference. We used the deep abundances defined in Table 1 and let molecules dissociate as a function of the pressure and temperature. This test is important to verify whether our model can find the input vertical structure of the atmosphere, without including day or night effects. When T<sub>upper</sub> ≠ T<sub>deep</sub>, the code should recover the input P<sub>iso</sub> and two different temperatures. When T<sub>upper</sub> = T<sub>deep</sub>, it should converge at the isothermal configuration.

We chose these parameters for the two-layer configuration to demonstrate how TauREx 2D can converge at the correct temperature structure, given that the composition remains the same in both hemispheres and with a temperature inversion between the lower and the upper part of the atmosphere. For such a test, we put the pressure threshold P<sub>iso</sub> between the two layers at 1 mbar, which corresponds to the range of pressures probed in such hot atmospheres, thus ensuring that the input spectrum should encode some information about the temperature change.
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Posterior distribution for the day and night case. The input spectrum refers to a hot atmosphere with three different equilibrium temperatures: lower temperatures, upper day, and upper night. As in the previous models, H₂O and H₂ can dissociate as a function of P (Pa) and T (K). The red lines represent the input values to generate the input spectrum.</figcaption></figure>

TauREx 2D found a $T_{side1} = 2500^{+700}_{-1400}$, $T_{side2} = 2600^{+600}_{-1400}$, $T_{deep} = 2594^{+8}_{-12}$, and $\log P_{iso} = -2.89^{+0.82}_{-0.51}$. TauREx favors a thin upper atmosphere and fits the lower layer with an isothermal at ~2600 K, which is the input upper temperature. This means that we managed to probe only the upper layer of the atmosphere. The retrieved molecular abundances found are compatible with the input values $\log [\text{H}_2\text{O}] = -3.33^{+0.01}_{-0.01}$ for water and $\log [\text{CO}] = -3.34^{+0.05}_{-0.04}$.

### 3.3. Day and night atmosphere

With this complex test, we explored the most complex atmospheric configuration in our test cases. Indeed, the input

spectrum was generated assuming an isothermal inner annulus at $T_{deep} = 2500$ K, a day-side temperature $T_{day} = 3300$ K, and a night side of $T_{night} = 500$ K with a $\beta$ angle of 20° (see Table 1 for the exact setting). The day and night thermal difference is realistic according to GCM models of this type of planet (e.g., WASP-121b) and, within this scenario, the $\beta$ parameter can give us a hint as to the 3D structure of the planet, since within this angle we probe both the day side and the night side around the terminator line. Based on these GCM equatorial temperature maps, the atmosphere can be approximated to an isothermal day and night sides, with a sharp transition region between the two sides (represented by the $\beta$ parameter), and an isothermal ring under a certain pressure. Within this scenario, the parameter can
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Posterior distribution for a particular case of the fully isothermal atmosphere with different water abundances in the day (log[H₂O]day = -5.3) and the night side (log[H₂O]night = -4.3).</figcaption></figure>

also give us a hint as to the 3D structure of the planet, since within this angle we probe both the day side and the night side around the terminator line. In Fig. 3, we show the best-fit model for this last test case.

In Fig. 4, we can see how TauREx finds a unique solution around the ground truth represented by the vertical red lines. With this final convergence test, we can see the reliability of TauREx 2D with handling spectra from bidimensional models. We let the P<sub>iso</sub> and β as free parameters and, we gave TauREx the possibility to look for the best spectrum.

The β parameter retrieved with TauREx 2D, in particular, would give us valuable information to model 3D effects of an

exoplanetary atmospheres. We see in Fig. 4 that TauREx found the parameters of the input spectrum.

**4. Intrinsic degeneracies in transmission spectra: how to recover independently the compositions of the day and night sides**

By construction, the transmission spectrum of a planet only depends on the total optical depth of the atmosphere along the various rays of light passing through the limbs (see Eq. (2)). Therefore, it can be seen right away that the transmission
<figure><img src="image_6.png" /><figcaption>Fig. 6. Posterior distribution for a retrieval test made using the same “day and night” input spectrum as in Fig. 4, but retrieved using a forward atmospheric model with independent compositions for the day and night sides. One can see the strong cross correlations between the abundances of a given molecule on the day and night sides.</figcaption></figure>

spectrum of a planet would remain unaffected if we took a mirror image of the planet, moving the day side from the star.

This intrinsic symmetry of the problem tells us that, technically, a retrieval algorithm cannot determine the conditions on the day and night sides. At best, it can infer the presence of different conditions on both hemispheres and we have to use our physical insight to decipher what conditions correspond to what side, for example saying that the hottest hemisphere is probably the dayside one. This is why we used $T_{side1}$ and $T_{side2}$ as independent parameters in our retrievals.

However, the fact that we measured only the total optical depth creates another degeneracy: to first order, only the total

column density of a species along the ray matters. To take an extreme example, if we had a completely isothermal planet with two hemispheres with different chemical compositions, the effect of a trace molecule (mol) on the spectrum would be unaffected as long as

$$[\text{mol}]_{\text{day}} + [\text{mol}]_{\text{night}} = \text{cst.} \quad (3)$$

A homogeneous atmosphere would have the same signature as an atmosphere where all the molecules would be in a single hemisphere.

To check the performance of TauREx 2D with a fully isothermal atmosphere and a compositional difference between the day
<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Spectral comparison between the two day and night retrieval cases. We can see that both the coupled model (the model that considers thermal dissociation in the atmosphere) and the “free” model (which is the one that considers a constant chemistry, but independent on both sides of the planets) lead to a consistent spectrum.</figcaption></figure>

side and the night side, we simulated an input synthetic spectrum as in the fully isothermal test case, with the difference that the water abundance changes from the day (log[H₂O]day = -4.3) to the night side (log[H₂O]night = -3.3). As shown in the posterior distributions (Fig. 5), TauREx 2D favors the solution in which the upper layer is flattened toward the upper edge of the planetary atmosphere, in a small annulus with a pressure lower than log[P](Pa) = -2.81+1.17-0.80, and it fits an isothermal atmosphere using only the lower layer. The input temperature is found within two sigma and the lower water abundance log[H₂O]deep = -3.53+0.04-0.05 is within one sigma from the arithmetic average between the input abundances log[H₂O]deep input = -3.56.

Of course, having different temperatures on both sides formally breaks this degeneracy because of the scale height effect as well as the temperature dependence of the molecular opacities (MacDonald & Lewis 2022). But the question remains whether there is enough information in a transmission spectrum to retrieve the composition of both sides independently.

To test this, we performed another retrieval on the “day and night” input spectrum used in Sect. 3.3, using TauREx 2D’s ability to model the planet with independent compositions for the day and night sides as described in Sect. 2 (“free” mode). The posterior distributions for this retrieval are shown in Fig. 6, and the best-fit spectra found by the two models (along with the residuals) are shown in Fig. 7.

As can be seen in the bottom panel of Fig. 7, the residuals are well below the 30 ppm noise floor used as uncertainty in these retrievals. The two approaches thus give an acceptable fit. However, it can be seen in Fig. 6 that the posterior distributions of the “free” retrieval show a lot of cross correlations between several parameters. In particular, as can be expected from the argument above, when the CO or H₂O abundances are high on one side, they can be very low on the other. This model also seems to miss the temperature dichotomy between the two hemispheres of the planet.

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Best-fit model for an atmosphere without optical absorbers. The TauREx 2D retrieval code found two possible solutions for this case. The comparison with the retrieval made with a 1D isothermal forward model shows that TauREx 2D leads to a better interpretation of the atmosphere.</figcaption></figure>

To see whether these two forward models can be statistically differentiated, we introduced the (logarithmic) Bayes factor:

$$ \mathcal{B} = \log \frac{E_d}{E_f} = \log E_d - \log E_f, \qquad (4) $$

where E_d and E_f are Bayesian evidence for the retrievals with dissociation (d) and free (f) models.

Given that Bayesian evidence of the retrieval that uses thermal dissociation and coupled day and night sides is E_d = 2914, while the one related to the model that uses independent day and night sides is E_f = 2908, it follows that the Bayes factor is ~6, giving a slight yet significant advantage to the coupled model. In this case, it seems that the higher number of free parameters penalizes the model with “free” composition.

This seems to show that there is not enough information in the input spectrum (at least with the level of uncertainty expected for a JWST-like observation) to retrieve the composition of both sides of the planet simultaneously and independently. Going further, the results even imply that the information on the temperature dichotomy in a transmission spectrum mainly comes from the link between temperature and composition in the observed atmosphere. Assuming chemical equilibrium brings a considerable amount of prior knowledge on the atmosphere. If this link seems robust in UHJs where molecular dissociation is fast enough to happen close to chemical equilibrium, it should be revisited for cooler atmospheres. For this reason, the ability to retrieve the temperature dichotomy on cooler planets might be greatly reduced compared to what is shown in Lacy & Burrows (2020), where they assume chemical equilibrium to hold down to day-side temperatures near 1200 K. This issue should be thoroughly investigated.

For all of these reasons, in the following, we continue to use the model assuming thermal dissociation at chemical equilibrium as it possesses fewer free parameters and results in less correlated posterior distributions between the parameters. We should, however, keep in mind that this assumption should be revisited when modeling cooler atmospheres.
<figure><img src="image_9.png" /><figcaption>Fig. 9. Posterior distribution for the case where we used a WASP-121b GCM’s synthetic atmosphere as input and employing only H₂O and CO as trace gases. The red lines correspond to the input deep abundances from the GCM model atmosphere.</figcaption></figure>

**5. A GCM retrieval**

After validating the TauREx 2D retrieval code on some self-
generated spectra, we tested its abilities on a more realistic atmo-
sphere, generated by the SPARC/MIT GCM and described in
detail in Parmentier et al. (2018). We simulated the transmission
spectra using 30ppm photon noise as described in Sect. 2. We
used the temperature map obtained by the GCM simulation
and computed the radiative transfer equation through an HJ
atmosphere with Pytmosph3R 2.0 (Falco et al. 2022). The rep-
resentative planet for this test is a WASP-121 b model, where
molecular species can thermally dissociate as a function of the
pressure and temperature. In this model, we included TiO and
VO in the atmosphere as they have been detected on WASP-121 b

by Evans et al. (2018). These molecules are optical absorbers
that heat the day side of the planet and create a strong thermal
inversion (Fortney et al. 2008; Parmentier et al. 2015, 2018).
Other molecules, ions, and metals (such as Iron detected in
the UHJ WASP-76 b by Ehrenreich et al. (2020)) also sig-
nificantly contribute to thermal inversion in UHJs, as pointed
out by Lothringer et al. (2018). Based on this GCM model
of WASP-121 b, we simulated two spectra from two different
atmospheric configurations to test its impact on the retrieval
analysis:

- A solar metallicity atmosphere with H₂O and CO as trace gases; and

- A solar metallicity atmosphere with H₂O, CO, TiO, and VO as trace gases.
<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Mixing ratio for H₂O and CO for the GCM case without optical absorbers. The dashed lines are the retrieved values, while the continuous lines represent the input value. The results from the 1D isothermal retrieval are in yellow. For this retrieval, the terminator and the night mixing ratios converged at the same values.</figcaption></figure>

The presence of TiO and VO does indeed cover some fundamental wavelengths where Rayleigh scattering otherwise dominates. In the absence of clouds, this wavelength range allows the model to calibrate the pressure in the atmosphere and thus to have better absolute abundances. We thus want to test whether hiding these regions leads to additional biases on the retrieved parameters.

A crucial parameter related to the three-dimensionality of the planet is the $\beta$ parameter. The $\beta$ angle is a 2D parameter and therefore its equivalence for a 3D GCM simulation is hard to define. For the following input spectra, the best approximation we can make is $5^{+3}_{-3}$ degrees based on the GCM equatorial temperature map, indicating a very sharp transition between the day and the night sides.

## 5.1. Atmosphere without optical absorbers

We started analyzing an atmosphere without the presence of optical absorbers (e.g., molecules such as TiO and VO). The retrieval performed by TauREx 2D converged to the best-fit solution shown in Fig. 8. In this test we can also see a degenerate solution converging at two different $\beta$ angles and temperatures for the two sides of the planet.

The two solutions found by TauREx have different physical properties and the best-fit parameters are shown in Fig. 9. TauREx 2D managed to converge around a realistic solution which is compatible with the ground truth, that is the parameters used to generate the input spectrum. Indeed, as shown by the red line in Fig. 9, the deep atmosphere trace gas parameter volume-mixing ratios are compatible with the ground truth of the input model within $1\sigma$ error bars for both solutions. In addition, the temperatures obtained in both solutions realistically reflect the ones of WASP-121b. This 2D model is more realistic and accurate than a simpler 1D model since it is able to catch day-night differences that, as shown in Pluriel et al. (2020), are crucial in HJs such as WASP-121b.

<figure><img src="image_11.png" /><figcaption>**Fig. 11.** PT profile retrieved compared to the input PT profile. In both scenarios, with and without TiO and VO, the PT profile with temperature inversion has been used. The night, terminator, and day-side PT profile of both the input (continuous line) and the retrieved models (dashed lines) are in blue, purple, and red, respectively. In the upper and lower part, we show the first and the second solution from the 2D retrievals, respectively. The solutions from a 1D retrieval using the four PT point profile (only for the GCM case with optical absorbers) and the isothermal PT profile are in green and yellow, respectively.</figcaption></figure>

We find $\beta_0 = 1.2^{+1.2}_{-0.9}$ and $\beta_1 = 11.4^{+2.5}_{-2.8}$ for the $\beta$ angle for the first and second solution, respectively, whose measures are correlated with the temperature values found in the two sides, that is $T_{\text{night}} = 190^{+150}_{-60}$ and $T_{\text{day}} = 2000^{+100}_{-100}$ for the first solution and $T_{\text{night}} = 980^{+80}_{-80}$ and $T_{\text{day}} = 2320^{+30}_{-30}$ for the second solution. Interestingly, both solutions give the same pressure level and the same temperature for the deep layer ($P_{\text{iso}} \sim 10$ Pa and $T_{\text{deep}} \sim 1780$ K). This means that the retrieval can probe approximately up to that atmospheric region, after which we approximated the atmosphere with an isothermal annulus of $T_{\text{deep}} \sim 1780$ K. The $\beta$ parameter gives us a hint as to the 3D structure of the planet since it spans above the day side and the night side, and a dramatic change in this parameter can lead to a dramatic change in the retrieved temperature on both sides of the planet. In both solutions, the transition between the day and the night side retrieved is sharp since it is less than $12^\circ$.

The mixing ratio retrieved by the 2D model and the PT profile approximates the input profiles better than the 1D model, as shown in Figs. 10 and 11. Both solutions from the 2D retrieval are closer to the input mixing ratio profiles than the solution of the 1D isothermal model. In particular, TauREx 2D found the correct abundance of water in the night side and the correct value of CO within a $2\sigma$ error.
<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Best-fit model for a GCM's synthetic atmosphere with TiO and VO as optical absorbers. TauREx 2D retrieval code found two possible solutions for this case (the red and the blue one for the solution with $\beta = 4.3$ and $\beta = 10.0$ degrees, respectively). The solutions from the 1D retrieval using an isothermal profile and a four PT point profile are shown in yellow and green, respectively.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 13.** Top and middle: temperature model for the two different solutions found in the GCM retrieval with optical absorbers. For visual reasons, the size of the atmosphere has been doubled. The temperature structure in the top represents the one of the first solution, the green one in Fig. 9. In the middle, there is the thermal structure corresponding to the second, violet solution of Fig. 14. Bottom: temperature map (equatorial cut) of WASP-121 b GCM simulation.</figcaption></figure>

To evaluate the contribution of a 2D retrieval model, we compared the solution given by a 2D model with the one provided with a 1D simple model, where we also took molecular dissociation along the vertical dimension of the atmosphere into account. From Pluriel et al. (2022), we saw how 1D models cannot interpret the input spectrum of 2D models. To quantify how compares the 2D model is to the 1D one, we can calculate the logarithmic Bayes factor expressed in Eq. (4).

The resulting logarithmic Bayes factor for the 2D model compared to the 1D one is $\mathcal{B} = +89.2$, which favors the solution of the 2D retrieval over the 1D one with a $\sigma$ significance higher than $5\sigma$ (Trotta 2008). For this test case, without molecules hiding the spectral part with a significant contribution from the Rayleigh scattering, the 2D retrieval converges at the input spectrum compatibly with the input parameters of the synthetic spectrum. An important goal of this result is to have managed to unravel the biases observed in Pluriel et al. (2020) in the molecular abundances of the atmospheric gases. In addition, the parametrization of our 2D retrieval model develops a thermal structure of the atmosphere which is compatible with 3D GCM models.

## 5.2. Atmosphere with optical absorbers

After studying the case of an atmosphere without TiO and VO as optical absorbers, we wanted to measure how the solution of a spectral retrieval in the presence of optical absorbers is impacted. Given the degeneracy of the TiO and VO and the Rayleigh scattering at short wavelengths and, depending on the spectral resolution and wavelength range of the data set, it could be common to not be able to disentangle the contribution of the Rayleigh scattering from the one of the optical absorbers correctly. The suggested Bayesian model has 11 free parameters and runs in ~ 20 days on 20 CPU cores using 1000 live points. The best-fit models found by TauREx 2D are shown in Fig. 12.

The results are similar to the previous retrieval: (i) TauREx 2D finds two different solutions that provide the same atmospheric spectrum (ii) converging around two beta angles, $\beta_0 = 4.3_{-0.9}^{+0.9}$ and $\beta_1 = 10.0_{-0.9}^{+0.8}$, (iii) with the expected deep abundance of each molecule getting considered. As we increase the degrees of freedom of the model (in our case, we have an additional spatial dimension compared to the 1D atmospheric model), the degeneracy of the solution could be a normal consequence that we need to handle. The $\beta_0$ angle found in this simulation is ~3 degrees higher than the $\beta_0$ solution found in the retrieval without optical absorbers. In the presence of TiO and VO, we can probe a more extensive part of the day side in the high atmosphere and, then, TauREx finds a greater retrieved $\beta$ angle for this latter case.

From Fig. 12 it is difficult to understand the contribution of the Rayleigh scattering under the TiO and VO absorption bands. In that spectral region, the H$_2$–H$_2$ CIA also provides a continuum contribution, with a peak at around 2.2 µm. Compared to the case without optical absorption, we have a precision decrease from 0.02 dex to 0.10 dex, which is related to the loss of a clear Rayleigh slope signal.
<figure><img src="image_14.png" /><figcaption>**Fig. 14.** Best-fit model for an atmosphere with TiO and VO as optical absorbers. The TauREx 2D retrieval code found two possible solutions for this case. The red lines correspond to the known deep input abundances in the GCM model.</figcaption></figure>

In Fig. 13, we show the retrieved thermal structure for this simulation, together with the equatorial cut view of the model. In Fig. 14, we can see that the 2D retrieval managed to converge around the input mixing ratios for the four molecules. The temperatures retrieved are $T_{\text{night}} = 780_{-80}^{+80}$ K and $T_{\text{day}} = 2730_{-40}^{+400}$ K for the first solution, $T_{\text{night}} = 690_{-40}^{+40}$ K and $T_{\text{day}} = 3470_{-120}^{+80}$ K for the second solution, and $T_{\text{deep}} \sim 1800$ K for both solutions.

Also, in this test, the volume-mixing ratio retrieved by the 2D model and the PT profile approximates the input profiles better than the 1D retrieval, as shown in Figs. 11 and 15. TauREx 2D found two solutions, approximating the input day

and night temperatures and estimating the deep atmosphere around its average temperature (Fig. 11). Moreover, TauREx 2D found a chemical distribution of elements closer to the input distribution than the one found with the 1D retrievals (Fig. 15).

The crucial result of this test is that TauREx 2D manages to retrieve the expected deep abundance of each molecule with a very good agreement to fit the input data. TauREx 2D correctly separates the contribution of the Rayleigh scattering from the one of the optical absorbers, also unraveling the biases explained in Pluriel et al. (2020). Also, in this case, a comparison with an equivalent 1D retrieval leads to a Bayes factor $\mathcal{B} = +47.9$ in favor of the 2D retrieval.
<figure><img src="image_15.png" /><figcaption>**Fig. 15.** Mixing ratio for H₂O, CO, TiO, and VO in the GCM test retrieval. The dashed lines are the values retrieved by TauREx, while the continuous lines represent the input value. The results from the 1D isothermal retrieval are in yellow and we show the solution found using a four-point PT profile forward model in green. For this retrieval, the terminator and the night-mixing ratios converged at the same values.</figcaption></figure>

<table><caption>**Table 3.** Comparison of models.</caption><thead><tr><th>Retrieval type</th><th>Number of parameters</th><th>Log Bayesian evidence</th><th>χ̃²</th><th>β with 1D iso</th><th>β with 1D four-point PT</th></tr></thead><tbody><tr><td>2D sol<sub>0</sub></td><td>11</td><td>2795.0</td><td>0.07</td><td>47.9</td><td>14.5</td></tr><tr><td>2D sol<sub>1</sub></td><td>11</td><td>2795.0</td><td>0.07</td><td>47.9</td><td>14.5</td></tr><tr><td>1D 4 PT</td><td>12</td><td>2780.5</td><td>0.18</td><td>33.4</td><td>—</td></tr><tr><td>1D iso</td><td>7</td><td>2747.1</td><td>0.45</td><td>—</td><td>-33.4</td></tr></tbody></table>

## 6. Comparison with a more complex 1D model retrieval

The 2D retrieval, in this case, helps us explore more on the 3D structure of the planet. We see that the two different solutions differ mainly from two possible β angles which lead to two possible values of $T_{day}$. We remark that a higher β angle results in a higher day-side temperature.

To quantify the goodness of the fit and see the advantages of a 2D retrieval over the “classic” 1D one, we fit the same input spectrum with a 1D model and calculated the Bayes factor between the two best-fit models. We tested the 1D model twice, assuming two different pressure-temperature profiles: one isothermal atmosphere and **one** four-point PT profile atmosphere.

In Table 3, we can find a resume of all the comparisons made between the 2D retrieval, including TiO and VO and the 1D retrievals. Together with the Bayesian evidence, we also provide the reduced $\tilde{\chi}^2$ of each best model which is defined as follows:

$$ \tilde{\chi}^2 = \frac{1}{N-p} \sum_{i=1}^{N} \left( \frac{O_i - C_i}{\sigma_i} \right)^2, \qquad (5) $$

where O and C are the observed and the computed spectral points, respectively, σ is the uncertainty, p is the number of

parameters in the model, and N is the total number of measurements. As described in Table 3, the 2D models approximate the input spectrum better, returning the lowest $\tilde{\chi}^2$ with respect with the 1D retrievals.

In Fig. 12, we can see that the four PT points retrieval leads to a better retrieval than the one done using an isothermal atmosphere. The logarithmic Bayes factor $\mathcal{B} = \log E_{4pt} - \log E_{iso} = 33$ confirms that the four-points retrieval is statistically better than the isothermal retrieval with a σ significance higher than $5\sigma$.

A comparison between the 2D and the four-points retrieval leads to a logarithmic Bayes factor of $\mathcal{B} = \log E_{2D} - \log E_{4pt} = 14$ and, thus, the 2D retrieval is also better than the four-points retrieval with a σ significance higher than $5\sigma$. The reduced $\tilde{\chi}^2$ of the four-points retrieval tells us that the atmospheric model found is compatible with the input spectrum.

All the molecular abundances appear to be overestimated and, thus, biased compared to the input parameters used to generate the input spectrum.

The calculation of the Bayes factor allowed us to quantify the goodness of the 2D retrieval compared to the 1D one. Indeed, if we consider only the reduced $\tilde{\chi}^2$ of the fit, we see that the single models are compatible with the input spectra with a high significance. However, it is impossible to get rid of the parameter biases when using the 1D model. The key point of the calculation of the
<figure><img src="image_16.png" /><figcaption>**Fig. 16.** Posterior distributions related to the four-points retrieval. Confirming the parameter biases found by Pluriel et al. (2022).</figcaption></figure>

logarithmic Bayes factor is to highlight the fundamental differences of the two model and the limits given by their geometry. With 30ppm error bars, it is indeed possible to reveal and explain some spectral features by introducing an additional geometric dimension on the atmospheric model.

Indeed, even though the 1D retrieval converges toward a statistically significant atmospheric model, the parameters found can differ from the ground truth by several orders of magnitude. In Fig. 16, as an example, we show the posterior distributions for a retrieval of the GCM input spectrum with optical absorbers, using a 1D atmospheric model and four

temperature-pressure points. The best model fits the input spectrum very well, but the retrieved amount of atmospheric trace gases differs from the ground truth by up to two orders of magnitude.

We see that the logarithmic Bayes factor defined in Eq. (4), which expresses the difference of the 2D model compared with the 1D one, is $\mathcal{B} = +47.9$. This favors the solution of the bi-dimensional model with a $\sigma$ significance higher than $5\sigma$. Also, in this case, the 2D model leads to a statistically better fit than the 1D model. The statistical significance of the 2D model is higher than that of the 1D model. In addition, it is less affected by the
parameter biases since it approximates the 3D geometry of the input model better.

The comparison with the 1D retrieval emphasizes the importance of running more-than-one dimensions retrieval models when using 30 ppm noise input spectra. In the era of the JWST, ARIEL, and Twinkle space missions it will be of fundamental importance to use better models which rely on more complex geometrical assumptions than the plane parallel atmospheres.

## 7. Conclusion

More-than-1D retrievals will be crucial to interpret the input spectral observation of new generation instruments better. Even though the approximation of a plane-parallel atmosphere can lead to a statistically acceptable fit, it inevitably leads to parameter biases that can be very far away from the reality (Pluriel et al. 2020, 2022). This is especially true for HJs and UHJs, where the great dichotomy between the day and the night side of the planet may lead to a misinterpretation of the input observations, without an accurate model.

We developed a 2D version of TauREx, which demonstrates that a 2D forward model is a good compromise between model accuracy and computational requirements. This model can solve the issues of the parameter biases induced by a 1D model which is essential in the context of future mission such as JWST, Twinkle, or Ariel. Thanks to TauREx 2D, we are able to determine the abundances of the species with much higher confidence than with 1D retrieval models.

We have shown how the model choice could be crucial to retrieve the thermal structure of an exoplanetary atmosphere properly (in Sect. 4), due to the intrinsic degeneracies of a high dimensional transmission model. Finally, in Sect. 6 we demonstrate how a 2D retrieval with an isothermal atmosphere returns a statistically more significant solution than a 1D retrieval using a more complex four-point pressure-temperature profile.

**Acknowledgements.** This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 679030/WHIPLASH). T.Zi acknowledges the funding support from Italian Space Agency (ASI) regulated by 'Accordo ASI-INAFF no. 2013-016-R.0 del 9 luglio 2013 e integrazione del 9 luglio 2015 CHEOPS Fasi A/B/C'.

## References

Al-Refaie, A. F., Changeat, Q., Waldmann, I. P., & Tinetti, G. 2021, ApJ, 917, 37

Baeyens, R., Decin, L., Carone, L., et al. 2021, MNRAS, 505, 5603

Barstow, J. K., Changeat, Q., Garland, R., et al. 2020, MNRAS, 493, 4884

Beichman, C., Benneke, B., Knutson, H., et al. 2014, PASP, 126, 1134

Brown, T. M. 2001, ApJ, 553, 1006

Burningham, B., Marley, M. S., Line, M. R., et al. 2017, MNRAS, 470, 1177

Caldas, A., Leconte, J., Selsis, F., et al. 2019, A&A, 623, A161

Changeat, Q., & Al-Refaie, A. 2020, ApJ, 898, 155

Changeat, Q., & Edwards, B. 2021, ApJ, 907, L22

Changeat, Q., Edwards, B., Waldmann, I. P., & Tinetti, G. 2019, ApJ, 886, 39

Charnay, B., Meadows, V., Misra, A., Leconte, J., & Arney, G. 2015, ApJ, 813, L1

Cowan, N. B., Greene, T., Angerhausen, D., et al. 2015, PASP, 127, 311

Deming, D. 2010, Nature, 468, 636

Drummond, B., Tremblin, P., Baraffe, I., et al. 2016, A&A, 594, A69

Ehrenreich, D., Lovis, C., Allart, R., et al. 2020, Nature, 580, 597

Espinoza, N., & Jones, K. 2021, AJ, 162, 165

Evans, T. M., Sing, D. K., Goyal, J. M., et al. 2018, AJ, 156, 283

Falco, A., Zingales, T., Pluriel, W., & Leconte, J. 2022, A&A 658 A41

Feng, Y. K., Line, M. R., Fortney, J. J., et al. 2016, ApJ, 829, 52

Feng, Y. K., Line, M. R., & Fortney, J. J. 2020, AJ, 160, 137

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601

Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419

Fortney, J. J., Shabram, M., Showman, A. P., et al. 2010, ApJ, 709, 1396

Gandhi, S., & Madhusudhan, N. 2018, MNRAS, 474, 271

Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485

Gordon, S., McBride, B., & Zeleznik, F. J. 1984, Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications. Supplement 1: Transport Properties

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

Heng, K., Frierson, D. M. W., & Phillipps, P. J. 2011, MNRAS, 418, 2669

Irwin, P. G. J., Teanby, N. A., de Kok, R., et al. 2008, J. Quant. Spec. Radiat. Transf., 109, 1136

Kataria, T., Showman, A. P., Fortney, J. J., et al. 2015, ApJ, 801, 86

Kataria, T., Sing, D. K., Lewis, N. K., et al. 2016, ApJ, 821, 9

Kreidberg, L. 2018, in Handbook of Exoplanets, eds. H. J. Deeg, & J. A. Belmonte, 100

Lacy, B. I., & Burrows, A. 2020, ApJ, 905, 131

Line, M. R., & Parmentier, V. 2016, ApJ, 820, 78

Line, M. R., Wolf, A. S., Zhang, X., et al. 2013, ApJ, 775, 137

Lothringer, J. D., Barman, T., & Koskinen, T. 2018, ApJ, 866, 27

MacDonald, R. J., & Lewis, N. K. 2022, ApJ, 929, 20

MacDonald, R. J., & Madhusudhan, N. 2017, ApJ, 850, L15

MacDonald, R. J., Goyal, J. M., & Lewis, N. K. 2020, ApJ, 893, L43

Madhusudhan, N., & Seager, S. 2009, ApJ, 707, 24

Marley, M. S., Saumon, D., Fortney, J. J., et al. 2017, in American Astronomical Society Meeting Abstracts, 230, 315.07

Menou, K., & Rauscher, E. 2009, ApJ, 700, 887

Mikal-Evans, T., Sing, D. K., Barstow, J. K., et al. 2022, Nat. Astron., 6, 471

Min, M., Ormel, C. W., Chubb, K., Helling, C., & Kawashima, Y. 2020, A&A, 642, A28

Molaverdikhani, K., Henning, T., & Mollière, P. 2019, ApJ, 883, 194

Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67

Moses, J. I., Line, M. R., Visscher, C., et al. 2013, ApJ, 777, 34

Nixon, M. C., & Madhusudhan, N. 2022, ApJ 935 73

Parmentier, V., Guillot, T., Fortney, J. J., & Marley, M. S. 2015, A&A, 574, A35

Parmentier, V., Fortney, J. J., Showman, A. P., Morley, C., & Marley, M. S. 2016, ApJ, 828, 22

Parmentier, V., Line, M. R., Bean, J. L., et al. 2018, A&A, 617, A110

Pluriel, W., Zingales, T., Leconte, J., & Parmentier, V. 2020, A&A, 636, A66

Pluriel, W., Leconte, J., Parmentier, V., et al. 2022, A&A, 658, A42

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

Rodgers, C. D. 2000, Inverse Methods for Atmospheric Sounding: Theory and Practice

Showman, A. P., & Guillot, T. 2002, A&A, 385, 166

Showman, A. P., Cooper, C. S., Fortney, J. J., & Marley, M. S. 2008, ApJ, 682, 559

Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564

Showman, A. P., Lewis, N. K., & Fortney, J. J. 2015, ApJ, 801, 95

Skemer, A. J., Morley, C. V., Zimmerman, N. T., et al. 2016, ApJ, 817, 166

Tan, X., & Komacek, T. D. 2019, ApJ, 886, 26

Taylor, J., Parmentier, V., Irwin, P. G. J., et al. 2020, MNRAS, 493, 4342

Tessenyi, M., Savini, G., Tinetti, G., et al. 2016, in AAS/Division for Planetary Sciences Meeting Abstracts, 48, 123.33

Tinetti, G., Drossart, P., Eccleston, P., et al. 2016, in Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave, SPIE, 9904, 658

Tinetti, G., Drossart, P., Eccleston, P., et al. 2018, Exp. Astron., 46, 135

Trotta, R. 2008, Contemp. Phys., 49, 71

Tsiaras, A., Waldmann, I. P., Zingales, T., et al. 2018, AJ, 155, 156

Wakeford, H. R., Visscher, C., Lewis, N. K., et al. 2017, MNRAS, 464, 4247

Wakeford, H. R., Sing, D. K., Deming, D., et al. 2018, AJ, 155, 29

Waldmann, I. P., Rocchetto, M., Tinetti, G., et al. 2015a, ApJ, 813, 13

Waldmann, I. P., Tinetti, G., Rocchetto, M., et al. 2015b, ApJ, 802, 107

Wardenier, J. P., Parmentier, V., & Lee, E. K. H. 2022, MNRAS, 510, 620

Welbanks, L., & Madhusudhan, N. 2022, ApJ, 933, 79

Welbanks, L., Madhusudhan, N., Allard, N. F., et al. 2019, ApJ, 887, L20

Wordsworth, R. D., Forget, F., Selsis, F., et al. 2011, ApJ, 733, L48

Zhang, M., Chachan, Y., Kempton, E. M. R., & Knutson, H. A. 2019, PASP, 131, 034501