A library of self-consistent simulated exoplanet atmospheres

Jayesh M. Goyal, Nathan Mayne, Benjamin Drummond, David K. Sing, Eric
Hébrard, Nikole Lewis, Pascal Tremblin, Mark W. Phillips, Thomas
Mikal-Evans, Hannah R. Wakeford, et al.

► To cite this version:

Jayesh M. Goyal, Nathan Mayne, Benjamin Drummond, David K. Sing, Eric Hébrard, et al. A library of self-consistent simulated exoplanet atmospheres. Monthly Notices of the Royal Astronomical Society, 2020, 498 (4), pp.4680-4704. 10.1093/mnras/staa2300. hal-04488833

HAL Id: hal-04488833

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

Submitted on 22 May 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_1.png" /></figure>

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

# A library of self-consistent simulated exoplanet atmospheres

Jayesh M. Goyal<sup>1,2*</sup> Nathan Mayne<sup>1,2</sup> Benjamin Drummond,<sup>2,3</sup> David K. Sing,<sup>2,4,5</sup> Eric Hébrard,<sup>2</sup> Nikole Lewis,<sup>1</sup> Pascal Tremblin,<sup>6</sup> Mark W. Phillips,<sup>2</sup> Thomas Mikal-Evans<sup>1,7</sup> and Hannah R. Wakeford<sup>1,8</sup>

Accepted 2020 July 29. Received 2020 July 29; in original form 2020 April 7

## ABSTRACT

We present a publicly available library of model atmospheres with radiative-convective equilibrium pressure–temperature (P–T) profiles fully consistent with equilibrium chemical abundances, and the corresponding emission and transmission spectrum with R ~ 5000 at 0.2 µm decreasing to R ~ 35 at 30 µm, for 89 hot Jupiter exoplanets, for four recirculation factors, six metallicities, and six C/O ratios. We find the choice of condensation process (local/rainout) alters the P–T profile and thereby the spectrum substantially, potentially detectable by James Webb Space Telescope. We find H⁻ opacity can contribute to form a strong temperature inversion in ultrahot Jupiters for C/O ratios ≥ 1 and can make transmission spectra features flat in the optical, alongside altering the entire emission spectra. We highlight how adopting different model choices such as thermal ionization, opacities, line-wing profiles and the methodology of varying the C/O ratio, effects the P–T structure, and the spectrum. We show the role of Fe opacity to form primary/secondary inversion in the atmosphere. We use WASP-17b and WASP-121b as test cases to demonstrate the effect of grid parameters across their full range, while highlighting some important findings, concerning the overall atmospheric structure, chemical transition regimes, and their observables. Finally, we apply this library to the current transmission and emission spectra observations of WASP-121b, which shows H₂O and tentative evidence for VO at the limb, and H₂O emission feature indicative of inversion on the dayside, with very low energy redistribution, thereby demonstrating the applicability of library for planning and interpreting observations of transmission and emission spectrum.

**Key words:** techniques: spectroscopic – planets and satellites: atmospheres – planets and satellites: composition – planets and satellites: gaseous planets.

## 1 INTRODUCTION

The thermal or the pressure–temperature (P–T) structure of a planetary atmosphere is a result of constant feedback between radiative, advective, and chemical processes. Determining the P–T structure of a planet is important for understanding underlying thermochemical and dynamical processes. The P–T structure of the planetary atmosphere also governs its spectral signatures, when remotely observed using telescopes, spacecrafts, or satellites. Therefore, it is necessary to constrain the P–T structure of a planet’s atmosphere, to understand the various physical processes occurring within them as shown in Fig. 1. Constraining the 3D P–T structure of the planetary atmosphere is ideally required, but the complexity and the computational resources required for such a model, especially when constructing a library of forward model simulations, motivates 1D P–T profiles.

In Goyal et al. (2018, 2019a,b), we used isothermal P–T profiles and corresponding equilibrium chemical abundances for computing

transmission spectra for a wide range of exoplanet atmospheres. This extensive library of more than a million model simulations is publicly available and has proved to be very useful in interpreting observations of various exoplanets using the Hubble Space Telescope (HST), Spitzer, and Very Large Telescope (see for e.g Wakeford et al. 2018; Alam et al. 2018; Zhang et al. 2019; Carter et al. 2020) and for planning future observations<fnref n="1" /> using various telescope facilities, including the James Webb Space Telescope (JWST). However, the assumption of an isothermal P–T profiles is only accurate for a small region of the atmosphere, specifically the high-altitude, low-pressure regions probed by transmission spectra. In fact, the atmosphere even in this region may rarely be exactly isothermal, but current transmission spectra observations cannot differentiate between such small temperature changes in P–T profiles (Fortney 2005; Heng & Kitzmann 2017; Goyal et al. 2018), which may be revealed by higher precision, higher resolution, and broader wavelength coverage observations. Furthermore, truly isothermal P–T profiles

<footnote n="1"><sup>1</sup>Department of Astronomy and Carl Sagan Institute, Cornell University, 122 Sciences Drive, Ithaca, NY 14853, USA</footnote>

<footnote n="2"><sup>2</sup>Astrophysics Group, School of Physics and Astronomy, University of Exeter, Exeter EX4 4QL, UK</footnote>

<footnote n="3"><sup>3</sup>Met Office, Fitzroy Road, Exeter, EX1 3PB, UK</footnote>

<footnote n="4"><sup>4</sup>Department of Earth and Planetary Sciences, Johns Hopkins University, Baltimore, MD 21210, USA</footnote>

<footnote n="5"><sup>5</sup>Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA</footnote>

<footnote n="6"><sup>6</sup>Centre detude de Saclay, Malson de la Simulation, CEA-CNRS-INRIA-UPS-UVSQ, USR 3441, France</footnote>

<footnote n="7"><sup>7</sup>Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA</footnote>

<footnote n="8"><sup>8</sup>School of Physics, University of Bristol, HH Wills Physics Laboratory, Tyndall Avenue, Bristol BS8 1TL, UK</footnote>

<footnote>* E-mail: jgoyal@astro.cornell.edu</footnote>

<footnote n="1">¹https://exoctk.stsci.edu/generic</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Illustrative figure to show dominant processes in different parts of a hot Jupiter atmosphere with radiative-convective equilibrium *P-T* profiles, with (black) and without inversion (blue).</figcaption></figure>

would give rise to a simple blackbody emission spectrum, devoid of spectral features. The emission spectrum is much more strongly and directly dependent on the temperature than the transmission spectrum (emission flux ∝ *T*⁴). Moreover, emission spectrum reveals the day/nightside of the exoplanet atmosphere, providing altogether different information than the transmission spectrum, which measures the limb. Emission spectrum also probes deeper (higher pressure) regions of the planetary atmosphere in comparison to transmission spectrum. Therefore, to identify features in the emission spectrum as well as to constrain the *P-T* profile using these features, we require computation of a more accurate non-isothermal *P-T* profile.

1D *P-T* profiles of irradiated H₂/He dominated planetary atmospheres with sufficiently high equilibrium temperatures are expected to reach a radiative-convective equilibrium condition (see for e.g. Marley et al. 1996; Burrows et al. 1997; Seager & Sasselov 1998; Sudarsky, Burrows & Hubeny 2003; Marley & Robinson 2015; Mollière et al. 2015; Gandhi & Madhusudhan 2017; Malik et al. 2019). This is basically a result of two factors, first the exposure to strong irradiation from the host star, rapidly forcing perturbations in the *P-T* structure back to a radiative equilibrium in the lower pressure regions (<~100 bar), and second the dominance of convection in the deep atmosphere forcing perturbations in the *P-T* structure to convective equilibrium in the higher pressure regions (>~100 bar). Therefore, combining the lower pressure regions, and higher pressure regions, these atmospheres are likely to exist close to a radiative-convective equilibrium, termed RCE hereafter. Additionally, assuming chemical equilibrium to model hot Jupiter atmospheres is likely a reasonable assumption, due to the high temperature of these planets, especially where chemical timescales are likely to be short, particularly for temperatures above ~2000 K. The assumption of chemical equilibrium is also a good starting point to constrain the main atmospheric constituents for lower temperature planets (see for e.g. Burrows & Sharp 1999; Lodders & Fegley 2002; Visscher, Lodders & Fegley 2006; Marley & Robinson 2015; Madhusudhan et al. 2016).

RCE *P-T* profiles for different planets can vary depending on their gravity, host star distance and spectral type, circulation in the planet's atmosphere, and the chemical composition of the atmosphere. This can then lead to a wide range of possible spectra for a given planet, governed by its *P-T* profile. Therefore, a library of RCE *P-T* profiles and the corresponding chemical abundances and simulated spectra is required to interpret the observations of exoplanet atmospheres and constrain the important physical processes occurring within them. In this work, we compute *P-T* profiles in RCE for 89 observationally

significant exoplanets, along with their corresponding equilibrium chemical abundances, simulated transmission and emission spectra with resolution of *R* ~ 5000 at 0.2 µm decreasing to *R* ~ 35 at 30 µm, and contribution functions (CF), which are publicly available here.<fnref n="2" />,<fnref n="3" />

We start by describing the 1D–2D atmosphere model ATMO (Amundsen et al. 2014; Tremblin et al. 2015, 2016, 2017; Drummond et al. 2016; Goyal et al. 2018) and the procedure to compute RCE *P-T* profiles, followed by recent additions to ATMO in Section 2. We detail the numerical setup of this grid with RCE *P-T* profiles in Section 3, including the chemistry and opacity setup. In Section 4 we describe the parameter space of the grid. We discuss the results in Section 5, where we show the sensitivity of model simulations to different model choices by comparing RCE *P-T* profiles, equilibrium chemical abundances, and the transmission/emission spectra in Section 5.1. In Section 5.2, we show the effects caused by high levels of irradiation by using extremely irradiated hot Jupiter WASP-121b as the test case. Sensitivity to grid parameters, namely; recirculation factor, metallicity and C/O ratio is discussed is Section 5.3 using WASP-17b and WASP-121b as test cases. In Section 6, transmission and emission spectrum observations of WASP-121b are interpreted using the grid of models presented in this work. Finally, we summarize and conclude in Section 7. The implementation and validation of thermal ionization, H⁻ and Fe opacity used in ATMO is detailed in Appendix A.

## 2 MODEL DESCRIPTION

As in Goyal et al. (2018, 2019b), in this work we use ATMO, a 1D–2D RCE model for planetary atmospheres. A nice comparison between ATMO and other models, before some of the additions described in this paper, is shown in Baudino et al. (2017) and Malik et al. (2019). In this section, we describe some of the details of ATMO along with recent additions to the model.

### 2.1 Computing *P-T* profiles in radiative-convective equilibrium

We here discuss the basic methodology used in ATMO to compute *P-T* profiles in RCE, consistently with equilibrium chemistry. Radiative transfer and absorption cross-section calculations along with equilibrium chemistry computations performed iteratively, all contribute to the computation of consistent RCE *P-T* profiles, discussed in detail in Drummond et al. (2016) and Goyal (2019). We initialize the model with an estimated (initial guess) *P-T* profile (the starting profile can also be isothermal although this will likely take longer computation time to reach convergence). We note that the final RCE *P-T* in independent of the initial starting *P-T* profile. Chemical equilibrium abundances are then calculated for this *P-T* profile. Using these abundances along with the absorption cross-sections of different chemical species, in the form of *k*-coefficients (Amundsen et al. 2014), the total opacity is calculated for each layer of the atmosphere. Following this, the radiative transfer equation is solved to compute incoming and outgoing radiative fluxes for each layer of the atmosphere. The same approach is followed for convective flux. This is then checked for energy flux balance in each layer of the atmosphere as well as the atmosphere as a whole, using the energy

<footnote n="2">²https://drive.google.com/drive/folders/1zCCe6HICuK2nLgnYJFal7W4lyu
njU4JE</footnote>

<footnote n="3">³https://noctis.erc-atmo.eu:5001/fsdownload/hq0z4udQJ/goyal2020</footnote>
conservation equation given by

$$F_{\text{rad}} + F_{\text{conv}} - \sigma T_{\text{int}}^4 = 0, \quad (1)$$

where $F_{\text{rad}}$ the radiative flux and $F_{\text{conv}}$ is the convective flux, respectively, as derived in Drummond et al. (2016) and Goyal (2019). $T_{\text{int}}$ represents the internal temperature of planet or the temperature at which the planet is cooling. Isotropic scattering is also included in the radiative flux computation. The P-T profile is also constrained by hydrostatic equilibrium, which defines the pressure structure as a function of altitude, implemented using

$$\frac{d}{dz}(P_{\text{gas}} + P_{\text{turb}}) - \rho g = 0, \quad (2)$$

where $P_{\text{gas}}$ is pressure due to gaseous species, $P_{\text{turb}}$ is turbulent pressure, $z$ is altitude, $\rho$ is density, and $g$ is gravity. If the conditions for energy conservation and hydrostatic equilibrium are not satisfied within the required numerical accuracy, corresponding to an error in flux balance (typical value of ~$10^{-3}$), the P-T profile is perturbed within minimum and maximum step sizes (~0.1–0.9) for the next iteration. This step size is basically the multiplying factor to the temperature perturbation while the model is iterating to obtain a converged solution. All the steps are repeated for each iteration, until a P-T profile that satisfies hydrostatic equilibrium and conservation of energy is obtained, consistent with equilibrium chemistry and corresponding opacities, for a given set of planetary characteristics.

Fig. 1 shows RCE P-T profiles with and without a temperature inversion. Unlike a static table of chemical abundances used in various models, these P-T profiles (and all models in this grid) are consistent with equilibrium chemical abundances, as explained above. In this schematic, we also show the dominant physical processes in each region. Convection plays an important role only in the deepest parts of the irradiated hot Jupiter atmospheres (>100 bar), which current observations cannot probe. Radiation governs the P-T structure across a wide range of pressures. Radiative diffusion leads to an isothermal structure in the deep atmosphere as shown in Fig. 1, primarily because the high opacity in this region decreases the mean-free path for the photons as in stellar atmospheres (Rybicki & Lightman 1986). The interplay between absorption of stellar radiation due to optical opacities and planetary emission due to infrared opacities, governs the P-T profile in the radiative region between ~1 bar to 0.1 mbar, which is the only region probed by current observations, either using transmission or emission spectra. Sufficiently high optical opacity can lead to a temperature inversion in this observed radiative region. At the top of the atmosphere (ToA) where the pressure becomes less than ~0.1 mbar the atmosphere becomes optically thin and most of the radiation penetrates this part of the atmosphere.

## 2.2 Recirculation factor

In a planetary atmosphere advection due to winds is one of the major processes responsible for transporting energy. To incorporate the 3D effect of advection in 1D models as adopted by Fortney & Marley (2007), we simply reduce the incoming flux in the 1D column of the atmosphere by a factor called the ‘recirculation factor’, hereafter termed $f_c$. It parametrizes the redistribution of input stellar energy in the planetary atmosphere, by the dynamics, where a value of 1 equates to no redistribution, while 0.5 represents efficient redistribution. The value of 0.5 $f_c$ indicates 50 per cent of the total incoming stellar energy is advected to the nightside (the side of the planet facing away from the star), while 0.25 $f_c$ indicates 75 per cent of the total incoming stellar energy is advected to the nightside. It

must also be noted that an additional factor, the incidence angle $\theta_o$ also contributes to the reduction in this total incoming stellar energy. $\theta_o = 60^\circ$ equating to the dayside (the side of the planet facing towards the star) average is the most commonly adopted value of incidence angle, contributing to 50 per cent reduction in the total incoming stellar energy, since $\cos 60^\circ = 0.5$.

## 2.3 Contribution function

The emission spectrum represents the ToA flux at different wavelengths for a given planet. However, this is a combination of flux from the different layers of the atmosphere. To identify the levels of the atmosphere contributing the most to this ToA emission, the CF (Chamberlain & Hunten 1987; Knutson et al. 2009; Drummond et al. 2018) is calculated given by,

$$CF = B(v, T) \frac{d(e^{-\tau_v})}{d(\ln(P))}, \quad (3)$$

where $v$ is the frequency, $T$ is the temperature, $B(v, T)$ is the Planck emission, $\tau_v$ is the optical depth, and $P$ is the pressure at each level of the atmosphere. The vertical P-T profile and the wavelength-dependent optical depth, are the primary quantities required to calculate the CF. Optical depth is a function of transmittance which decreases as we go deeper in the planet’s atmosphere. Therefore, the CF is higher in the region where there is a larger change in the optical depth or transmittance for a unit change in the pressure (altitude) over the same region. In simpler terms, the CF peaks in the region where the wavelength-dependent optical depth is one, when the planet is being probed from the ToA. Although we define CF here, for ease of understanding and plotting we compute the normalized CF (NCF), by normalizing using the largest value of the CF along the P-T profile.

## 2.4 High temperature additions

We include thermal ionization chemistry in ATMO by including ion species while computing equilibrium chemical abundances. We also included iron (Fe) and $H^-$ opacity. We detail their implementation and validation in Appendix A.

## 3 NUMERICAL SETUP FOR THE GRID

Computing RCE P-T profiles consistent with equilibrium chemistry is not a trivial task. The P-T profile and the chemistry are intricately linked as they depend on each other. The chemical composition is largely dependent on the temperature, and the temperature is largely dependent on composition (via opacities). Moreover, the P-T profile as well as the chemical abundance profile continuously change as the simulation progresses towards the solution. In such a scenario a large number of temperature and pressure points are encountered. Therefore, it is extremely difficult to obtain converged (satisfying all constraints) RCE P-T profiles for values across a large parameter space. There are always the regions of the parameter space where the simulations tend not to reach a converged solution due to many factors such as the boundary conditions, numerical instabilities, non-convergence of equilibrium chemistry especially with condensation and many more. We deal with such problems by incrementally adjusting the numerical setup for some of the failed model simulations as described in the next section.

We use 50 vertical model levels with a maximum optical depth of $5 \times 10^5$ at 1 µm. Since ATMO calculates quantities on an optical
depth grid, the minimum and maximum optical depths govern the pressure domain (extent of the atmosphere). An increase in the maximum optical depth leads to an increase in the pressure domain of the *P-T* profile, for a given set of parameters. The model stability when solving for RCE, consistently with equilibrium chemistry, is very sensitive to the selected ToA (minimum) optical depth boundary condition, as the atmosphere can become very sparse (less dense) in this region. Therefore, the ToA optical depth is varied to achieve convergence. The typical values used for ToA optical depth are 10⁻², 10⁻<fnref n="5" />, 10⁻¹⁰, and 10⁻¹³ at 1 µm. Although this is a wide range for ToA optical depth, these extremely low optical depth regions are outside the domains of the region probed by either transmission or emission spectra. Moreover, we find that a value of 10⁻<fnref n="5" /> is sufficient for most of the model simulations. The ToA pressure is restricted to 10⁻<fnref n="6" /> bar, which corresponds to the ToA minimum optical depth. Even though we vary the minimum optical depth to achieve convergence, the pressure is always set at 10⁻<fnref n="6" /> bar for this minimum optical depth, serving as a reference for the atmospheric *P-T* profile.

We use 32 band correlated-*k* cross-sections (Amundsen et al. 2014; Goyal et al. 2018) for generating consistent RCE *P-T* profiles and 5000 bands to generate transmission spectra, emission spectra, and CF. The internal temperature of the planet (*T*<sub>int</sub>) defined in Section 2.1 is set at 100 K following (Guillot & Showman 2002; Fortney, Marley & Barnes 2007). However, we note that following some recent studies (Thorngren, Gao & Fortney 2019; Sing et al. 2019), this assumption might be debatable. We adopt a mixing length constant α = 1.5 for calculating convective flux (Baraffe et al. 2015), as used in previous ATMO simulations (e.g Drummond et al. 2016). To standardize the comparison of transmission spectra for a range of variables, we set the pressure at which the radius of the planet is defined at 1 mbar (Lecavelier Des Etangs et al. 2008), only while computing the transmission spectra. We note that there exists a degeneracy between reference transit radius and associated reference pressure as highlighted by Lecavelier Des Etangs et al. (2008).

The target planet selection technique for this library of models is the same as that adopted in Goyal et al. (2018) with the additional constraint of the planetary equilibrium temperature. In this library only the planets with equilibrium temperatures greater than 1200 K, as computed in the TEPCAT data base (Southworth 2011) are selected, all shown in Table S3 of the online Supporting Information. The choice of 1200 K is based on the capability of ATMO to obtain converged solution (RCE *P-T* profile consistent with equilibrium chemistry) for a range of grid parameter space values. It was found that for planets with equilibrium temperatures less than 1200 K a large number of model simulations across the grid parameter space fail to converge, most likely due to numerical instabilities arising due to inclusion of various condensates in our equilibrium chemistry computation.

The input stellar spectra for each planetary model grid are taken from the BT-Settl<sup>4</sup> models (Allard, Homeier & Freytag 2012; Rajpurohit et al. 2013). These stellar spectra are selected according to the closest obtained host star temperature, gravity and metallicity from the TEPCAT data base (Southworth 2011). All the parameters required for model initialization like stellar radius, planetary radius, planetary equilibrium temperature, surface gravity, and semimajor axis are also adopted from the TEPCAT<sup>5</sup> data base, shown in the table 1 of the online Supporting Information.

## 3.1 Chemistry setup

Similar to Goyal et al. (2018, 2019b), we restrict our calculations to equilibrium chemistry in this work, but with self-consistent RCE *P-T* profiles. However, in this work we also include thermal ionization in the equilibrium chemistry computation, detailed, and validated in Appendix A1. Therefore, in addition to the 258 species as used in Goyal et al. (2018, 2019b) which are listed here,<sup>6</sup> the list of species for equilibrium computation now also includes H<sup>+</sup>, H<sup>-</sup>, Na<sup>+</sup>, K<sup>+</sup>, e<sup>-</sup>, C<sup>+</sup>, He<sup>+</sup>, Ca<sup>+</sup>, and Si<sup>+</sup> ions along with additional gaseous species NaF, KF, SiO, SiS, CaH, CaOH, and condensate species Na<sub>3</sub>AlF<sub>6</sub> in three different forms and LiF in the crystalline form. Thus the total number of species used in the equilibrium chemistry calculation adds up to 277 for the simulations presented here.

## 3.2 Opacity setup

The potential of a particular species to absorb/emit radiation at a particular spectral interval is governed by its absorption cross-sections. Compared to the simulations presented in Goyal et al. (2018, 2019b), the opacities due to H⁻ and Fe have been newly included in this work, in addition to H<sub>2</sub>-H<sub>2</sub> and H<sub>2</sub>-He collision-induced absorption (CIA) opacities, and opacities due to H<sub>2</sub>O, CO<sub>2</sub>, CO, CH<sub>4</sub>, NH<sub>3</sub>, Na, K, Li, Rb, Cs, TiO, VO, FeH, CrH, PH<sub>3</sub>, HCN, C<sub>2</sub>H<sub>2</sub>, H<sub>2</sub>S, and SO<sub>2</sub>. The implementation of H⁻ opacity is described in Appendix A2 and that of Fe in Appendix A3. The line-list sources for all opacities used in this library and the type/sources of pressure broadening, are shown in Tables S1 and S2 of the Supporting Information, respectively. The effects of these newly added opacities, H⁻ and Fe on the *P-T* profiles and thereby the spectra are discussed in Sections 5.2.2 and 5.2.3, respectively.

Absorption cross-sections for all the species included in ATMO in each of the 5000 correlated-*k* bands at 1 bar and 2000 K are shown in Fig. 2. These plots aid identifying the major absorbing species in various parts of the spectrum. However, it must be noted that the final opacities are the product of absorption cross-sections and chemical abundances. Therefore, although species such as TiO, VO, FeH, and Fe have strong absorption cross-sections in the optical their contribution to total absorption will be zero if they condense or do not form at any given temperature.

## 4 GRID PARAMETER SPACE

The library of models with RCE *P-T* profiles and corresponding equilibrium chemical abundances, transmission spectra, emission spectra, and CF for various planets are computed at four different recirculation factors (*f*<sub>c</sub> = 0.25, 0.5, 0.75, 1.0), six metallicities (0.1, 1, 10, 50, 100, 200; all in × solar), and six C/O ratios (0.35, 0.55, 0.7, 0.75 1.0, 1.5), giving a total of 144 model RCE *P-T* profiles and spectra per planet. However, as mentioned earlier, not all models in the parameter space achieve convergence, thus leading to absence of some models in the library (~5–10 per cent) which varies for different planets. The RCE *P-T* profiles and corresponding equilibrium chemical abundances used to generate our transmission and emission spectra are the same, that is they span a similar parameter space. Since we are presenting here a library of model atmospheres and corresponding spectra, they are intended to cover a wide range of possibilities and all the extremes.

<footnote n="4"><sup>4</sup>https://phoenix.ens-lyon.fr/Grids/BT-Settl/AGSS2009/SPECTRA/</footnote>

<footnote n="5"><sup>5</sup>http://www.astro.keele.ac.uk/jkt/tepcat/allplanets-ascii.txt</footnote>

<footnote n="6"><sup>6</sup>https://drive.google.com/drive/folders/1g7Bc6pbwvLUDf-QFJOCEOFKX6glKLuqF</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Absorption cross-sections ($k_{band}$) of all species used in ATMO in each of the 5000 correlated-$k$ bands at 1 bar and 2000 K. These are absolute cross-sections without any dependance on chemical abundances. (a) Cross-sections for H$_2$–H$_2$ (blue), H$_2$–He (green), H$_2$O (cyan), CO$_2$ (red), CO (deep sky blue), CH$_4$ (lawn green), NH$_3$ (saddle brown), Na (purple), K (magenta), and Li (orange). (b) Cross-sections for Rb (blue), Cs (khaki), TiO (cyan), VO (red), FeH (deep sky blue), PH$_3$ (magenta), H$_2$S (saddle brown), HCN (purple), C$_2$H$_2$ (black), SO$_2$ (orange), and Fe (grey).</figcaption></figure>

The choice of recirculation factor ($f_c$) covers all possible scenarios from no recirculation (1.0) to extremely fast winds (0.25). We do not extend the grid to metallicities greater than 200× solar, because above this metallicity the atmosphere becomes abundant in species other than H$_2$ and He, such as CO$_2$, H$_2$O, CO, etc. This would require the inclusion of pressure broadening effects due to these species, and no existing studies have solved this problem due to lack of lab-based observational data (Fortney et al. 2016). The choice of C/O ratios in the grid is guided by important transition regimes as found by previous studies (Madhusudhan 2012; Mollière et al. 2015; Goyal et al. 2018; Molaverdikhani, Henning & Mollière 2019).

As described in Goyal et al. (2018, 2019a,b), there are two approaches to treat condensation in our library of models, rainout, and local condensation. While computing equilibrium chemical abundances to obtain RCE *P*-*T* profiles with rainout condensation, each layer is dependent on other layers, specifically only on layers that lie at higher pressures, in contrast to the local condensation approach where each layer is independent. This makes the assumption of rainout with RCE *P*-*T* profiles more realistic for a planetary atmosphere as compared to just local condensation. Therefore, we generate RCE *P*-*T* profiles with rainout condensation for the library of models presented in this work. However, we show the differences in the *P*-*T* profiles and the spectra due to the different approaches in the next section.

The structure of the pressure broadened line wings of Na and K can have a substantial effect on the *P*-*T* profiles and thereby the emission spectrum of brown dwarfs and hot Jupiter exoplanet atmospheres (Burrows, Marley & Sharp 2000; Allard et al. 2003). Even with their high-resolution measurements, the shape of the pressure broadened wings are still a matter of debate for brown dwarfs (Burrows et al. 2002; Burgasser et al. 2003; Allard et al. 2003). For hot Jupiters, we have very recently started to observationally probe the line wings of Na/K (Nikolov et al. 2018). In Section 2.1 of the online Supporting Information of this paper we show the differences in the *P*-*T* profiles and the spectra due to these two line-wing profiles for Na and K. The differences are negligible and unlikely to be detectable by observations. We adopt the Na and K line-wing profiles from Allard et al. (2003) for the library of models presented in this work, which includes detailed quantum mechanical calculations while computing these profiles (Amundsen 2015).

Although convection plays an important role in determining the *P*-*T* profile of brown dwarfs, we see for hot Jupiter exoplanet atmospheres the effect of convection on the *P*-*T* structure of the observable atmosphere is negligible. This is because of the strong irradiation on these planets from their host stars, which reduces the radiative time-scale, thus making the atmospheric *P*-*T* profile almost entirely dependent on the ToA irradiation (along with the atmospheric pressure, temperature, and heat capacity), at-least in the region where observations can probe the atmosphere. Therefore, we conclude that it is not necessary to include convection while computing RCE *P*-*T* profiles for irradiated hot Jupiter exoplanet atmospheres. However, in this work all the *P*-*T* profiles include parametrized convection for completeness, as it is computationally inexpensive.

As described in Drummond et al. (2019), there are three different methodologies to vary C/O ratio relative to the solar C/O ratio. In the simulations presented in Goyal et al. (2018, 2019b), we varied C/O ratio by varying O/H. However, varying C/O by varying C/H can lead to differences in *P*-*T* profiles and equilibrium chemical abundances as shown by Drummond et al. (2019). Therefore, we investigated the effect of varying C/O ratio, by varying O/H and those by varying C/H using WASP-17b as the test case, with the results presented in the online Supporting Information Section 2.3. Although there are some differences in the results obtained using these two different methodologies, in the parameter space we consider, they are smaller compared to the effects of other parameters in the grid and other model choices (e.g. local or rainout condensation). Therefore, to enable fair comparison between the model spectra generated in our previous works (Goyal et al. 2018, 2019b) and keep the library of models consistent with them, we again in this work adopt the methodology of varying C/O ratio by varying O/H. Ideally, we could use O/H and C/H as separate parameters in the grid, however, that increases the size of the grid substantially and makes it computationally expensive for a large number of exoplanets. Therefore, we select one methodology over the other.

## 5 RESULTS

In this section, we initially show the effects of some of the important model choices on the RCE *P*-*T* profiles and thereby the equilibrium
chemical abundances, transmission and emission spectra, which
could effect interpretation of observations. Additional tests showing
the sensitivity of some of the other model choices such as the Na and
K line-wing profiles, VO line-list sources, and the methodology to
vary C/O ratio either by varying O/H or C/H, are all detailed in the
online Supporting Information of this paper. Decoupled emission
spectrum are also shown and described in the online Supporting
Information, to aid identification of different spectral features in the
emission spectra. In the following sections, we show the sensitivity
to the choice of condensation process (rainout/local), high levels of
irradiation due to thermal ionization, addition of H⁻ and Fe opacity,
and formation of temperature inversion in the atmosphere. Finally,
we show the sensitivity of the model simulations, that is P-T profiles,
chemical abundances, spectra and the CF, to all the grid parameters
across their full range, using WASP-121b and WASP-17b as test
cases.

**5.1 Sensitivity to model choices**

5.1.1 Comparing model simulations with rainout and local condensation

Goyal et al. (2019b) have demonstrated that the adoption of either a
rainout or local condensation approach results in differences in the
transmission spectra when adopting an isothermal *P*-*T* profile. Here,
we investigate the effects of these two condensation approaches on
the RCE *P*-*T* profiles and thereby the equilibrium chemical abun-
dances, and transmission/emission spectra. Figs 3(a)–(e), show the
RCE *P*-*T* profiles, equilibrium chemical abundances, transmission
spectra, emission spectra, and NCF, respectively, for WASP-17b at
0.25 *f*<sub>c</sub>, solar metallicity and solar C/O ratio, adopting both rainout
and local condensation. A difference of ~100–300 K can be seen
between the 1 and 100 mbar pressure levels in Fig. 3(a), with the
*P*-*T* profile adopting the local condensation approach having higher
temperatures. However, in the deeper atmosphere (~0.1–1000 bar)
the *P*-*T* profile with rainout condensation is hotter by ~200 K. This
difference in temperature leads to lower abundance of Na and K
in the upper atmosphere and a higher abundance of TiO/VO in the
local condensation case, as compared to rainout condensation case,
shown in Fig. 3(b). This higher TiO/VO abundance can also be seen
in the transmission spectra shown in Fig. 3c where the spectra with
local condensation show TiO/VO features in the optical, missing
in the rainout condensation case. This also strengthens the findings
of Goyal et al. (2019b), that TiO/VO features can reveal dominant
physical process (rainout or local condensation) in the planet’s
atmosphere.

The differences in the emission spectra shown in Fig. 3(d) are
substantial, primarily due to the differences in the *P*-*T* structure. At
1.1 µm deeper parts of the atmosphere (~1 bar) are being probed
as seen in the NCF shown in Fig. 3(e), mainly because none of the
species have a strong opacity to absorb in this region. Therefore, it
can be seen in the emission spectra that the flux is higher for the
rainout condensation as compared to local condensation, since the
temperature at ~1 bar is higher for the simulation adopting rainout
condensation. At 1.4 µm, the wavelength of one of the strongest
water opacity bands, we see the NCF moves to the upper atmosphere
and the difference between the pressure levels being probed using
emission spectrum in the case of the rainout (0.1 bar) and local
condensation (0.01bar) simulations are also substantially different.
This is mainly due to the large difference in the temperatures between
the rainout and local condensation simulations between 0.1 and
0.01 bar. The difference in the peak pressure level of the NCF as

observed between 1.1 and 1.4 µm, the wing and core of water band,
respectively, can be observed more strongly at 2.25 and 2.8 µm, since
this is the region of peak emission for a body with temperature similar
to the equilibrium temperature of WASP-17b. Since the differences
in the spectra are substantial, it might be possible to distinguish
between emission spectra due to rainout and local condensation, and
therefore constrain the *P*-*T* profiles and thereby the condensation
processes using JWST.

**5.2 Effects caused by high levels of irradiation**

5.2.1 Effect of thermal ionization

Thermal ionization of certain species can have a substantial effect on the chemistry of the planetary atmospheres, depending on the atmospheric temperature. Fig. 4(a) shows the equilibrium chemical abundances of certain important species for an extremely irradiated hot Jupiter WASP-121b, with (dashed line) and without (solid line) thermal ionic species included in the equilibrium chemistry computation. H⁻ opacity has been included in the model simulation with thermal ionic species and Fe opacity in both. The abundance of Na decreases by about 3 orders of magnitude in the upper atmosphere (transmission spectra probed region) when thermal ionic species are included, as Na is ionized to form Na⁺. Fig. 4(a) shows that the abundance of Na⁺ becomes almost equal to Na without thermal ionic species model simulation, in the upper atmosphere (P < 10⁻³ bar). Similar behaviour can be observed for K and K⁺ ions. The abundance of H₂O also drastically decreases for P < 10⁻³ bar, however, is still substantial in the deeper atmosphere probed by emission spectrum. Inclusion of thermal ionic species also has effects on the transmission spectra, where the narrow Na and K features (cores) seen in the model spectra without thermal ionic species, disappear in the model spectra with thermal ionic species as shown in Fig. 4(b). We note that the equilibrium abundance of H₂O at very high temperatures and low pressures as for WASP-121b is almost similar (see Fig. 4a), with and without thermal ionic species, because, thermal decomposition of H₂O is taken care of in the equilibrium chemistry computation (by neutral products H₂ and O₂), even if ions are not included in the computation.

5.2.2 Effect of H⁻ opacity

H⁻ opacity contributes to the absorption of radiation in hot Jupiters
via bound-free and free-free cross-sections as explained in Ap-
pendix A2. There is also a strong observational evidence for H⁻
opacity in Mikal-Evans et al. (2019). To understand the effect of
H⁻ opacity we compute *P*-*T* profiles with and without H⁻ opacities
included in the model. We also do not include Fe opacity specifically
for models in this tests, to isolate the effects of H⁻ opacity. Fig. 5(a)
shows these *P*-*T* profiles at 0.5 *f*<sub>c</sub> with solar metallicity and C/O
ratio for WASP-121b. It can be noticed from the *P*-*T* profiles that
H⁻ opacity tends to cool the deeper atmosphere (>1 bar) by about
100 K which increases to 400 K with an *f*<sub>c</sub> value of 1 (not shown
here). This can be attributed to an increase in H⁻ abundance in the
lower atmosphere (>1 bar), with a higher *f*<sub>c</sub> value, as the temperature
increases. In Fig. 5(b), the *P*-*T* profile obtained using an *f*<sub>c</sub> value of
1.0, solar metallicity and a C/O ratio of 1.5 is shown. This has been
particularly chosen to show the extreme effects due to H⁻ opacity.
Fig. 5(b) shows that there is a substantial difference in *P*-*T* profiles
with and without H⁻ opacity.

Without H⁻ opacity there is a weak temperature inversion as
compared to that with H⁻ opacity. At high C/O ratios (> 1) as shown
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) Figure showing RCE *P*-*T* profiles for WASP-17b at 0.25 *f*<sub>c</sub> and solar metallicity and solar C/O ratio (0.55), with rainout (blue) and local condensation (red). (b) Figure showing equilibrium chemical abundances of important species using *P*-*T* profiles shown in (a), with rainout (solid) and local condensation (dashed). (c) Figure showing transmission spectra using *P*-*T* profiles shown in (a) and chemical abundances shown in (b), with rainout (blue) and local condensation (red). (d) Figure showing emission spectra using *P*-*T* profiles shown in (a) and chemical abundances shown in (b), with rainout (blue) and local condensation (red). (e) Figure showing CF for emission spectra shown in (d), with rainout (blue) and local condensation (red) at 1.1 µm (solid), 1.4 µm (dotted), 2.25 µm (dashed), and 2.8 µm (dotted-dashed).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Figure showing equilibrium chemical abundances of various species without (solid) and with (dashed) thermal ionic species (ionic species are shown using marked lines) included in the model for WASP-121b at solar metallicity, solar C/O ratio (0.55), and $0.5 f_c$. (b) Figure showing transmission spectra models, with and without thermal ionic species at solar metallicity, solar C/O ratio, and $0.5 f_c$.</figcaption></figure>

in Fig. 14(b), the abundance of TiO/VO decreases dramatically, the major absorbers likely to cause an inversion in extremely irradiated hot Jupiters like WASP-121b. Therefore, other species start contributing to form a temperature inversion. Without $H^-$ opacity the strong inversion is not sustained but there is a weak inversion due to Na and K (as shown in Mollière et al. 2015) discussed in detail in (ii) under Section 5.3.3. With $H^-$ opacity the strength of inversion increases dramatically at $f_c = 1.0$ and C/O ratio of 1.5 as shown in Fig. 5(b), because the $H^-$ abundance between pressure levels $\sim 1$ to $10^{-2}$ bar increases (see Fig. 14b for trends even though at $f_c = 0.5$), such that its opacity can create an inversion similar to TiO/VO. However, thermal inversion due to $H^-$ opacity lies deeper (higher pressure) in the atmosphere ($\sim 1$ to $10^{-2}$ bar), in comparison to the inversion due to TiO/VO ($\sim 10^{-1}$ to $10^{-3}$ bar). To investigate the potential of $H^-$ opacity by its own to form temperature inversions at different C/O ratios, we performed some model simulations by removing dominant optical absorbers such as TiO, VO, Na, K, Li, Rb, Cs, Fe, and FeH. We find $H^-$ opacity is able to create a strong thermal inversion for C/O ratios $\ge 1$, we also see a very weak temperature inversion ($\sim 100$ K) at C/O ratio of 0.75 and almost negligible at solar C/O ratio (0.55), highlighting the C/O ratio dependency of $H^-$ opacity to create thermal inversions.

The effect of $H^-$ opacity on the transmission spectrum is also substantial as shown in the Fig. 5(c) where it tends to obstruct the deeper atmosphere due to its high opacity as shown for a high C/O ratio case for WASP-121b. In this case it tends to mute the wings of Na and K mimicking the effect of cloud at optical wavelengths. Fig. 5(d) shows the emission spectrum for WASP-121b with and without $H^-$ opacity using the *P-T* profile shown in Fig. 5(b). Fig. 5(d) shows that there is a substantial difference in the emission spectrum with and without $H^-$ opacity, due to difference in *P-T* profiles.

### 5.2.3 Effect of Fe opacity

Gaseous iron (Fe) opacity dramatically affects the *P-T* profiles of extremely irradiated hot Jupiters (Lothringer, Barman & Koskinen 2018). It tends to heat the upper, lower pressure, atmosphere as shown in the Fig. 6(a) and cool the lower, high pressure, atmosphere.

The cross-section of Fe is quite large in the optical leading to strong absorption and therefore the heating in the upper atmosphere. This blocking of radiation higher up in the atmosphere is one of the causes of the cooling in the deeper atmosphere similar to $H^-$. The transmission spectra in Fig. 6(b), shows that Fe opacity is important in the optical part of the spectrum, generally dominated by Rayleigh scattering slope, with very sharp spectral features. The emission spectrum in Fig. 6(c) with and without Fe opacity, shows that the strength of the thermal inversion in the emission increases after adding Fe opacity. Thus, Fe opacity can play an important role in creating secondary inversions, but in the lower pressure regions ($\sim 0.1$ mbar) as compared to that created due to $H^-$ or TiO/VO opacity (see (ii) under Section 5.2.4 for more details). The absence of $H_2O$ due to thermal dissociation also contributes to this secondary inversion at extremely low pressures, as there is no strong infrared emitting species available to re-emit the energy absorbed by Fe. The difference between the emission spectrum shown in Fig. 6(c) is mainly due to the increase in the size of emission features, due to a larger temperature difference between the lower (high pressure) and the upper (low pressure) atmosphere in simulations including Fe opacity, as compared to those without Fe opacity.

### 5.2.4 Inversions in hot Jupiters without TiO/VO

The formation of temperature inversion in irradiated planetary atmosphere is governed by the interplay between absorbed stellar radiation and planetary emission. The presence of temperature inversions in hot Jupiter exoplanet atmospheres have been predicted for a long time (Hubeny, Burrows & Sudarsky 2003; Fortney et al. 2008), primarily due to TiO/VO opacities, as shown in Fig. 8(a). The presence of temperature inversion due to TiO/VO also leads to substantial changes in the emission spectrum of the planet as shown in Fig. 8(b). However, it is only very recently that a definitive evidence of an inversion layer has been seen in the atmosphere of an ultrahot Jupiter WASP-121b (Evans et al. 2017). Although, there is a tentative evidence of VO in the transmission spectrum of WASP-121b (Evans et al. 2018), it is still unclear what opacity source is causing this additional absorption creating a temperature inversion in WASP-121b, Therefore, in this section, the opacity required to
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) Figure showing RCE *P-T* profiles without (blue) and with (red) *H*⁻ opacity included in the model for WASP-121b at solar metallicity, solar C/O ratio, and 0.5 *f*<sub>c</sub>. Fe opacity is not included in both these models (see Section 5.2.2 for details). (b) Same as (a) but at C/O ratio of 1.5 and 1.0 *f*<sub>c</sub>. (c) Figure showing transmission spectra using *P-T* profiles shown in (b) and corresponding chemical abundances with and without *H*⁻ opacity. (d) Figure showing emission spectra using *P-T* profiles shown in (b) and corresponding chemical abundances with and without *H*⁻ opacity.</figcaption></figure>

produce an inversion, its impact on the *P-T* profile and thereby the emission spectra is investigated.

(i) *Temperature inversions with grey and H⁻ opacity:* following Burrows, Budaj & Hubeny (2008) and Spiegel, Silverio & Burrows (2009), we add an arbitrary grey absorbing opacity across the optical wavelengths (0.44–1 µm) throughout the atmosphere (all model layers) of WASP-121b with 0.8 *f*<sub>c</sub>, solar metallicity, and C/O ratio. Varying the magnitude of this opacity then allows us to explore the evolution of the *P-T* profile from being non-inverted to being inverted as a function of opacity, along with the evolution of the emission spectrum for these different atmospheric structures. The value of 0.8 *f*<sub>c</sub> is chosen as the best-fitting value to observations in Evans et al. (2017). For this particular test along with a grey absorbing opacity in the optical, only the H<sub>2</sub>–H<sub>2</sub> and H<sub>2</sub>–He CIA, and H<sub>2</sub>O, CO<sub>2</sub>, CO, CH<sub>4</sub>, NH<sub>3</sub>, Na, K, Li, Rb, Cs opacities were used for simplicity. Figs 7(a) and (b) show the *P-T* profile and emission spectra, respectively, for varying levels of grey opacity for WASP-121b. We note that we have omitted TiO/VO opacity in this model simulation. It can be seen that as the grey opacity increases

from 0.002 to 0.02 cm² g⁻¹ the *P-T* profile changes from being non-inverted to inverted. We also compute the ratio

$$ \gamma = \frac{\kappa_{vis}}{\kappa_{IR}}, \qquad (4) $$

where *κ*<sub>vis</sub> and *κ*<sub>IR</sub> are the mean opacities in the visible and infrared wavelengths, as used in previous works (Mollière et al. 2015; Gandhi & Madhusudhan 2019), to quantify the strength of the inversion. We compute *κ*<sub>IR</sub> using the Planck mean opacity for the temperature of each atmospheric level, whereas *κ*<sub>vis</sub> is computed using the Planck mean opacity at the stellar effective temperature of WASP-121, which is ~6500 K. At the photospheric level (~0.1 bar) the *γ* values are 4.78, 5.5, 6.06, and 8.3 for grey opacity values of 0.002, 0.004, 0.006, and 0.02 cm² g⁻¹, respectively. The strength of the inversion increases as the value of *γ* increases. However, we note that this method to quantify strength of inversion can be too simplistic, especially when comparing strength of inversion caused by different species, mainly due to two reasons. First, the altitude/pressure of the inversion can vary depending on the RCE *P-T* profile and the corresponding chemical equilibrium abundances, therefore using the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Figure showing RCE *P-T* profiles without (blue) and with (red) Fe opacity included in the model for WASP-121b at solar metallicity, solar C/O ratio, and 0.5 *f*<sub>c</sub>. (b) Figure showing transmission spectra using *P-T* profiles shown in (a) and corresponding chemical abundances with and without Fe opacity. (c) Figure showing emission spectra using *P-T* profiles shown in (a) and corresponding chemical abundances with and without Fe opacity.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) Figure showing *P-T* profiles for WASP-121b with 0.8 *f*<sub>c</sub>, solar metallicity and solar C/O ratio, with different amount of optical grey opacity added throughout the atmosphere. (b) Figure showing emission spectra with different amount of optical grey opacity using *P-T* profiles shown in (a).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** (a) Figure showing RCE *P-T* profiles with all opacities (blue) including Fe and without TiO/VO opacities (red), at solar metallicity and solar C/O ratio with 0.5 *f*<sub>c</sub> for WASP-121b. (b) Figure showing emission spectra using *P-T* profiles shown in (a) with all opacities (blue) including Fe and without TiO/VO opacities (red) for WASP-121b.</figcaption></figure>

$\gamma$ value at one particular pressure level may not be an accurate way to compare strength of inversions. Secondly, for hot planets like WASP-121b there will be some overlap between the visible (stellar spectrum) and infrared (planet spectrum), which again means this simplistic approach is not likely to provide accurate measure of the strength of the inversion.

The change in the emission spectrum is more interesting where the H₂O absorption feature at 1.4 µm gradually changes into an emission feature, as the amount of grey opacity is increased. Moreover, with 0.006 cm² g⁻¹ grey opacity the spectrum almost resembles a blackbody spectrum indicating an isothermal atmosphere, which can be seen in Fig. 7(a).

As discussed in detail in Section 5.2.2 and shown in Fig. 5(b), H⁻ opacity which almost acts like a grey opacity in the optical can produce a strong inversion in the absence of TiO/VO opacities, which only happens at high C/O ratios (>1), under chemical equilibrium conditions. Therefore, H⁻ opacity could possibly be the opacity source that leads to inversion in the atmosphere of WASP-121b and many other ultra hot Jupiter planets.

(ii) *Can Fe without TiO/VO form an inversion?* Iron (Fe) has very strong opacity in the UV and the optical part of the spectrum. Therefore, we investigated whether Fe can lead to inversions without TiO/VO in extremely irradiated hot Jupiter such as WASP-121b. For this we removed TiO/VO opacities from the model while computing RCE *P-T* profiles. We found that the Fe opacity we include is unable to produce a strong inversion like TiO/VO at solar metallicity as shown in Fig. 8. Even if we increase the metallicity to 200 times solar, leading to an increase in Fe abundance to about ~30 × than that at solar metallicity, it does not lead to an inversion. Fe opacity however, leads to a sharp increase in the temperature at lower pressures (~10⁻⁵ bar like a thermospheric secondary inversion). The primary reason behind this secondary inversion at low pressures is the low abundance of species that emit strongly in the infrared (for e.g. thermal decomposition of H₂O), which reduces atmospheric cooling (see Fig. 14b). We also observe this secondary inversion due to various other strong optical absorbers such as Na and K, when the abundance of strong infrared emitters such as H₂O or CH₄ decreases. This overall result due to Fe opacity might seem to contradict that of Lothringer et al. (2018), where they see an inversion due to Fe

opacity at around 10⁻² bar (10 mbar). However, Lothringer et al. (2018) used a fiducial model with a planet at 0.025 au from a star with effective temperature of 7200 K, about ~750 K hotter than WASP-121. Moreover, they include bound-free Fe opacity in their model simulations which becomes important for wavelengths shorter than 0.2 µm (Sharp & Burrows 2007). The combination of all these factors may have led to differences between our results, alongside other factors such as the equilibrium chemical abundances of Fe as well as the treatment of condensation. However, in KELT-9b the planet with highest equilibrium temperature in our library of models, Fe opacity leads to formation of a strong temperature inversion starting at ~0.1 bar, even at solar metallicity and without requiring TiO/VO/H⁻. Thus, the capability of Fe opacity to form a temperature inversion is system specific showing a strong dependance on the host star effective temperature.

We note that there are many other potential species such as AlO, CaO, SiO, NaH, AlH, SiH, MgH, etc. shown in recent studies (Lothringer et al. 2018; Malik et al. 2019; Gandhi & Madhusudhan 2019) with strong optical opacities that could lead to the formation of an inversion in ultrahot Jupiters. However, these opacities have not been included in the current version of our library. Among these species there have been indications of AlO in some of the planets (von Essen et al. 2019; Chubb et al. 2020). Therefore, some of these opacities will be considered in the future version of the library.

## 5.3 Sensitivity to grid parameters

In this section, we show the sensitivity of the model simulation, that is *P-T* profiles and the spectra, to all the grid parameters, namely recirculation factor, metallicity, and C/O ratio across their full range used in this library, using WASP-17b and WASP-121b as test cases.

### 5.3.1 Effect of recirculation factor

The recirculation factor ($f_c$) described in Section 2.2 governs the efficiency of re-distribution of energy (by winds) received from the host star in a column. The value of 1 corresponds to no-redistribution, with an increase in redistribution as this factor decreases. Here,
we show the effect of varying the recirculation factor on the P-T profiles and thereby the chemical abundances and spectra, using the extremely irradiated hot Jupiter WASP-121b as an example. As can in seen in the P-T profiles in Fig. 9(a) for WASP-121b, the strength of the atmospheric temperature inversion increases with $f_c$ as expected, since more energy is available to create an inversion at higher values of $f_c$. At 0.25 $f_c$, the inversion is absent in the P-T structure due to reduced irradiation (energy). Absorption due to TiO/VO is the primary reason for the inversion, but surprisingly the abundance of TiO/VO starts decreasing, as the inversion is formed and increases in strength, as can be seen in Fig. 9(b) for increasing $f_c$. However, the abundance of H$^-$ and Fe increases with increasing $f_c$. This increase in H$^-$ maintains the temperature inversion even though the abundance of TiO/VO decreases, as described in Sections 5.2.2 and (ii) under Section 5.3.3.

The transmission spectrum of WASP-121b in Fig. 9(c), shows that the strength of the H₂O features decreases with increasing $f_c$ as H₂O abundance decreases in the region where transmission spectra probes (~1 mbar), since H₂O starts becoming thermally unstable with increasing temperatures. However, CO features start appearing near ~2.5 µm and broad-band CO features between 4 and 6 µm. The strength of these CO features increases with increasing $f_c$. As expected, the flux in the planetary emission spectrum shown in Fig. 9(d) increases with increasing $f_c$ as the temperature of emission increases. The P-T profile at 0.25 $f_c$ is very close to isothermal, therefore its emission spectrum also resembles a blackbody curve, with small dips in the strong water vapour absorption bands at 1.4, 2, and 3 µm. In contrast, P-T profiles at other $f_c$ values have temperature inversions, thus leading to a bump instead of dip in the strong H₂O and CO absorption bands.

The NCF at 1.7 and 1.4 µm are shown in Figs 9(e) and (f), respectively, for WASP-121b. At 1.7 µm, the NCF peaks deeper in the atmosphere as compared to that at 1.4 µm, indicating emission at 1.7 µm is from comparatively deeper parts of the atmosphere, since 1.7 µm is at the edge of strong water absorption band centred at 1.4 µm. For the profiles with a temperature inversion ($f_c = 0.5, 0.75, \text{and } 1.0$), the deeper and cooler isothermal part (not the inversion) is primarily probed at 1.7 µm as compared to 1.4 µm which probes the inversion layer, thus leading to a emission feature (bump) instead of a absorption feature (dip) in the emission spectrum from 1.2 to 1.7 µm. This also happens for other strong water absorption bands as shown in the emission spectra. The 1.4 µm emission feature has been detected in the atmosphere of WASP-121b (Evans et al. 2017), but other such potential emission features indicative of inversion for wavelengths > 2 µm are still to be detected for WASP-121b and can only be possible with JWST.

At 0.25 $f_c$, the inversion is absent as we reduce the irradiation received from the host star to 25 per cent of its original value, mimicking the transport of energy by advection (strong winds). This also motivates accurate 3D modelling of extremely irradiated hot Jupiter exoplanets to predict inversions as well as to infer wind velocities based on the presence or absence of inversions.

### 5.3.2 Effect of metallicity

Metallicity fundamentally effects the chemistry of the atmosphere and thereby its P-T structure and observed spectrum. The effect of metallicity on the P-T structure, chemistry, transmission, and emission spectra for two different planets is discussed in this section.

(i) WASP-17b: WASP-17b is a hot Jupiter planet with an equilibrium temperature of 1755 K. Adopting a value of 0.5 for the $f_c$

the transmission spectra for this planet shows TiO/VO features due to hotter P-T structure as compared to that using 0.25 $f_c$. However, observations from Sing et al. (2016) shows the absence of TiO/VO features, therefore we restrict the following analysis to simulations adopting an $f_c$ of 0.25 (although all values are available in the model grid). The C/O ratio is fixed to the solar value and the metallicity varied from across the grid range to investigate the effect of metallicity. For this planet, an increasing metallicity leads to an increase in the temperature throughout the atmosphere as shown in Fig. 10(a). This is a result of increased absorption of radiation at lower pressures, due to increased opacity, as the mean molecular weight of the atmosphere increases, driven by an increase in the abundances of species heavier than H₂ and He, namely H₂O, CO₂, and Na as shown in Fig. 10(b). A sharp increase in temperature for subsolar metallicity can be seen at pressures less than ~$10^{-5}$ bar, but it is observationally insignificant either in transmission or emission spectra, due to the low atmospheric density at these pressures.

The effect of varying the metallicity on the transmission of WASP-17b is shown in Fig. 10(c). The transmission spectra are primarily dominated by Na, K, and narrow Li features in the optical for all metallicities, with weaker TiO/VO features as metallicity increases. The infrared part of the spectrum is primarily dominated by H₂O features for all metallicities. The size of the H₂O features initially increases with an increase in metallicity as H₂O abundance increases shown in Fig. 10(b), but then again decreases with increasing metallicity. This is caused by the increase in the mean molecular weight of the atmosphere, leading to a decrease of the atmospheric scale height, which, in turn, shrinks the spectral features in transmission. Pressure broadened wings of Na and K are also effected by change in metallicity. Due to decreasing scale height associated with increasing metallicity, transmission spectra probes high pressure levels of the atmosphere, resulting in enhanced broadening of Na and K line wings with increasing metallicity. The CO₂ feature near 4.5 and 2.5–3 µm increases in amplitude, which can also be seen in the emission spectrum, primarily due to a rapid increase in the CO₂ abundances. This shows that even under chemical equilibrium conditions the atmosphere rapidly tends to migrate towards a CO₂ abundant atmosphere with increasing metallicity (Moses et al. 2013), offering potential reasons for the CO₂ dominated compositions of Mars and Venus in our Solar system. Even Earth in the past may have had a CO₂ dominated atmosphere, currently captured in the oceans and rocks by various geological processes (Zahnle, Schaefer & Fegley 2010). The transmission spectra can also be compared with Fig. 12(a) in Goyal et al. (2019a) which uses isothermal P-T profiles. The comparison might not be completely accurate as the RCE P-T profiles shown here are overall cooler than the equilibrium temperature of WASP-17b (1755 K) in the transmission spectra probed region (~1 mbar). The ‘best-fitting’ model transmission spectra of WASP-17b using the grid of models presented in this paper and observations from Sing et al. (2016) is shown in Fig. S6 of the online Supporting Information.

The emission spectra of WASP-17b with varying metallicity is shown in Fig. 10(d). The blackbody emission flux increases with increasing metallicity, since the overall temperature of P-T profiles shown in Fig. 10(a) increases with increasing metallicity. This increase in temperature combined with the increase in abundances of species such as H₂O and CO₂, with increasing metallicity, leads to deeper absorption features. Similar to the transmission spectrum, the emission spectrum is also dominated by water absorption features for WASP-17b in the infrared with CO₂ features around 4.5 and 2.5–3 µm region for metallicities greater than 10 times solar value. The atmospheric level which contributes most to the emission at different wavelengths can
<figure><img src="image_11.png" /><figcaption>**Figure 9.** (a) Figure showing P-T profiles for a range of *f*<sub>c</sub> values (0.25, 0.5, 0.75, and 1) at solar metallicity and solar C/O ratio (0.55) for WASP-121b (b) Figure showing equilibrium chemical abundances for some important species for various *f*<sub>c</sub> values obtained using P-T profiles shown in (a). (c) Figure showing transmission spectra for various *f*<sub>c</sub> values obtained using P-T profiles shown in (a) and corresponding equilibrium chemical abundances shown in (b). (d) Figure showing planetary emission spectra for various *f*<sub>c</sub> values obtained using P-T profiles shown in (a) and corresponding equilibrium chemical abundances shown in (b). (e) Figure showing NCF at 1.7 µm for a range of *f*<sub>c</sub> values for emission spectra as shown in (d). (f) Same as (e), but at 1.4 µm.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** (a) Figure showing *P-T* profiles for a range of metallicities at 0.25 *f*<sub>c</sub> and a solar C/O ratio (0.55) for WASP-17b. (b) Figure showing equilibrium chemical abundances for some important species for various metallicity values obtained using *P-T* profiles shown in (b). (c) Figure showing transmission spectra for WASP-17b for different values of metallicity obtained using *P-T* profiles shown in (a) and corresponding equilibrium chemical abundances shown in (b). (d) Figure showing emission spectra for WASP-17b for different values of metallicity obtained using *P-T* profiles shown in (a) and corresponding equilibrium chemical abundances shown in (b). (e)–(h) Figures are showing NCF at 2.25, 2.8, 3.8, and 4.5 for a range of metallicity values for emission spectra shown in (d).</figcaption></figure>

be found using the NCF. The increase in the overall temperature of
the *P-T* profile causes the NCF to consistently shifts towards lower
pressure with increase in metallicity as shown in Figs 10(e)–(h). The
NCF at the core of the CO₂ absorption band at 4.5 µm shown in the

Fig. 10(h), peaks at lower pressures (<~1 mbar) in the atmosphere,
for metallicities greater than 50 times solar, than the wing of the
band at 3.8 µm as shown in the Fig. 10(g). The dramatic drop in
the emission flux between these two wavelengths is also shown in
the emission spectra for metallicities greater than 50 times solar. The NCF for the 2.8 µm H₂O absorption band also shows similar effect as 4.5 µm CO₂ absorption band, shown in Figs 10(e) and (f), for 2.25 and 2.8 µm, respectively. The atmospheric level probed by emission spectrum shifts to lower pressure levels with increasing metallicity, as shown by CF in the H₂O and CO₂ bands. The wings and cores of absorption/emission bands also probe different atmospheric layers for a given metallicity. Therefore, the absorption/emission features of various species can be used to constrain metallicity, as well as the *P*-T temperature structure of the atmosphere. For some cases (for e.g. Fig. 10h) particularly at high metallicity, NCF shows that the emission spectrum is probing quite high in the atmosphere (~10⁻⁵ bar). In such cases, photochemistry can alter the emission spectrum via the formation of photochemical products and change in chemical abundances (and thereby the *P*-T profile) at high altitudes (*P* < ~10⁻⁴ bar) (Madhusudhan et al. 2016; Hobbs et al. 2019). However, since photochemistry is a dis-equilibrium effect, it is not considered in this work and is part of our future work as discussed in the conclusions.

(ii) **WASP-121b**: WASP-121b is an extremely irradiated hot Jupiter planet with an equilibrium temperature (*T*<sub>eq</sub>) of ~2400 K (Southworth 2011). Figs 11(a)–(d) show the *P*-T structure, chemical abundances, transmission, and emission spectra, respectively for WASP-121b at 0.5 *f*<sub>c</sub> and solar C/O ratio. Figs 11(e)–(h) show the CF at 2.25, 2.8, 3.8, and 4.5 µm, respectively.

For WASP-121b we see a temperature inversion in the *P*-T profile, which moves towards higher pressure levels with increasing metallicit-y, primarily driven by TiO/VO absorption. At subsolar metallicity the inversion is very weak, due to the low abundance of TiO/VO, seen in Fig. 11(b). The Fe abundance is also substantial, mainly contributing to upper atmosphere heating, as the Fe absorption cross-sections are largest in the UV-optical spectrum. This upper atmosphere heating leads to the formation of a second inversion layer but with weak observational signatures due to the very low density of the atmosphere in this region.

The transmission spectra shown in Fig. 11(c) is dominated by TiO/VO features in the optical, the size of which decreases with increasing metallicity due to a reduction in the atmospheric scale height. Sharp Fe features dominate the optical spectrum shortward of ~0.4 µm with H₂O features in the infrared. The CO features are also seen in the transmission spectrum particularly around 2.5 and 4–6 µm, which can also be seen in the emission spectrum in Fig. 11(d). In the emission spectra shown in Fig. 11(d) for WASP-121b, due to an inversion layer in the *P*-T profile, most of the molecular features are seen as emission features as opposed to the absorption features seen for WASP-17b. The amplitude of the CO features increases with increasing metallicity. The TiO/VO features can also be seen as emission features in the emission spectra in the optical. The H₂O emission features dominate the infrared, where the 1.4 µm feature has led to the detection of an inversion layer for the first time in an exoplanet atmosphere (Evans et al. 2017). It can also be noticed from the NCF that the wings of strong absorption bands at 2.25 (H₂O) and 3.8 µm (CO) shown in Figs 11(e) and (g), respectively, mainly probe the region below the inversion layer, while the cores of absorption bands at 2.8 and 4.5 µm shown in Figs 11(f) and (h) probe the inversion layer.

### 5.3.3 Effect of C/O ratio

(i) **WASP-17b**: the *P*-T profiles for WASP-17b for a range of C/O ratios are shown in Fig. 12(a) at 0.5 *f*<sub>c</sub> and solar metallicity. We choose to show simulations with 0.5 *f*<sub>c</sub> instead of 0.25 *f*<sub>c</sub> as one of

the models with 0.25 *f*<sub>c</sub> and a C/O ratio of 1.0 failed to converge in the grid. It can be seen that with increasing C/O ratio the *P*-T structure consistently cools for *P* ≤ 10⁻¹ bar. However, for *P* ≥ 10⁻¹ bar, the *P*-T structure first cools upto C/O ratio of 0.75 and then the temperature increases for a C/O ratio of 1 and 1.5. The sharp heating at around 10⁻⁶ bar is due to Fe opacity as explained earlier in (ii) under Section 5.2.4.

The change in the equilibrium chemical abundances due to the change in the C/O ratio is drastic, as it effects all the major carbon and oxygen bearing molecules. As expected the abundances of H₂O drop with increasing C/O ratio. Although CO₂ bears a carbon atom it needs two oxygen atoms per carbon atom, therefore the equilibrium abundance of CO₂ also drops with increasing C/O ratio, but in smaller increments as compared to H₂O. The abundance of carbon bearing species such as CH₄, HCN, C₂H₂ increases with increasing C/O ratio, while the abundance of CO is almost constant, since it has one atom of carbon and oxygen each. This transition from H₂O dominated spectra, to spectra dominated by various carbon bearing species occurs between C/O ratios of 0.75 and 1, slightly higher than that found with isothermal *P*-T profiles where it was between 0.7 and 0.75 as shown in Goyal et al. (2018). However, this value might change with change in the *f*<sub>c</sub> value as the C/O transition point is a strong function of temperature (Mollière et al. 2016; Goyal et al. 2018; Molaverdikhani et al. 2019).

This C/O transition is also seen in the transmission spectrum shown in Fig. 12(c), where the spectrum transitions from being H₂O dominated to being dominated by CH₄, HCN, and C₂H₂ between a C/O ratio of 0.75 and 1. Figs 13(a) and (b) show this transmission spectrum at a C/O ratio of 1.5 decoupled into various molecules. It can be seen that at a C/O ratio of 1.5, the transmission spectrum is dominated by CH₄ features in the infrared, with contributions from CO, HCN, and C₂H₂. There is a strong HCN and C₂H₂ feature at ~1.6 µm and the most common CO feature at 4.5 µm. The emission spectrum shown in Fig. 12d also shows this C/O transition between 0.75 and 1, from deep H₂O absorption features to more deeper CH₄ absorption features, in the peak region of emission around ~3 µm

(ii) **WASP-121b**: the *P*-T profiles for WASP-121b for a range of C/O ratios are shown in Fig. 14(a). With increasing C/O ratio the major temperature inversion shifts to higher pressures. The major temperature inversion refers to the inversion that has a potential observational signature unlike the inversion due to Fe opacity at extremely low pressures (~10⁻⁵ bar) as described in Sections 5.2.3 and (ii) under Section 5.2.4 and also seen in Fig. 14(a).

The abundance of TiO/VO species which are the primary absorbers for forming the temperature inversion layer, decrease with increasing C/O ratio. Their abundance at a C/O ratio of 1.0 and 1.5 is low, still the inversion layer is maintained, albeit not as hot as at other low C/O ratios. We investigated this phenomenon by using two tests, first removing TiO/VO opacities and second removing Na, K, TiO, and VO opacities. We found that this inversion at high C/O ratio can be maintained due to Na and K opacities in the absence of TiO/VO opacities or their low abundance, as found by Mollière et al. (2015). Furthermore, H⁻ opacity also contributes to this inversion at high C/O ratio and high value of *f*<sub>c</sub> (1.0) more than Na and K opacities, as discussed in detail in Section 5.2.2. The increase in the H⁻ abundance at C/O ratio of 1.5 around the photosphere region (~0.1 to 1 bar) is the primary reason for this inversion due to H⁻ opacity. H⁻ features also can be clearly seen in the transmission spectrum for C/O ratio of 1.5.

At high C/O ratios it is interesting to see that the HCN abundance is substantial even in the low pressure regions (*P* > ~1 mbar). This can also be noticed in the transmission spectrum shown in Fig. 14(c) and
<figure><img src="image_13.png" /><figcaption>**Figure 11.** (a) Figure showing *P-T* profiles for a range of metallicity at 0.5 *f*<sub>c</sub> and a solar C/O ratio (0.55) for WASP-121b. (b) Figure showing equilibrium chemical abundances for some important species for various metallicity values obtained using *P-T* profiles shown in (a). (c) Figure showing transmission spectra for WASP-121b for different values of metallicity obtained using *P-T* profiles shown in (a) and corresponding equilibrium chemical abundances shown in (b). (d) Figure showing emission spectra for WASP-121b for different values of metallicity obtained using *P-T* profiles shown in (a) and corresponding equilibrium chemical abundances shown in (b). (e)–(h) Figures showing NCF at 2.25, 2.8, 3.8, and 4.5 µm for a range of metallicity values for emission spectra shown in (d).</figcaption></figure>

decoupled transmission spectrum shown in Figs 15(a) and (b), which is dominated by HCN features in the infrared at a C/O ratio of 1 and 1.5, as opposed to CH<sub>4</sub> for WASP-17b. This shows an important result that at high temperatures HCN dominates over CH<sub>4</sub> in the

atmosphere at high C/O ratios. Therefore, HCN features provide a very strong signature to constrain high C/O ratios in exoplanet atmospheres. The transition from an H<sub>2</sub>O dominated spectrum to that dominated by HCN happens between C/O ratio of 0.75–1.0.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** (a) Figure showing *P-T* profiles for a range of C/O ratios at 0.5 *f*<sub>c</sub> and solar metallicity for WASP-17b. (b) Figure showing equilibrium chemical abundances for some important species for various C/O values obtained using *P-T* profiles shown in (a). (c) Figure showing transmission spectra for WASP-17b for different values of C/O ratios obtained using *P-T* profiles shown in (a) and corresponding equilibrium chemical abundances shown in (b). (d) Figure showing emission spectra for WASP-17b for different values of C/O ratios obtained using *P-T* profiles shown in (a) and corresponding equilibrium chemical abundances shown in (b).</figcaption></figure>

Therefore, the temperature dependence of the C/O transition as seen in Mollière et al. (2015) and Goyal et al. (2018) does not seem to be holding in this case, as the transition for WASP-17b also happens in this C/O regime (0.75–1.0). The transmission spectrum also shows FeH features between 0.8 and 1.2 µm at a C/O ratio of 1.0 and 1.5. The absence of TiO/VO features makes it possible for FeH to appear, without being masked.

In the emission spectrum shown in Fig. 14(d) most of the molecular features are seen in emission due to the presence of an inversion layer, as explained before. However, between a C/O ratio of 1 and 1.5 surprising differences can be seen, especially at 3.1 µm which is the wavelength of a strong HCN absorption band. At a C/O ratio of 1 it is an emission feature, however it transforms to an absorption feature at a C/O ratio of 1.5. This is because at a C/O ratio of 1.5 slightly cooler upper atmosphere is being probed as can be seen in the NCF at 3.1 µm in Fig. 15(d) as compared to a C/O ratio of 1. Therefore, at a C/O ratio of 1.5, HCN absorbs the radiation from the lower (high pressure atmosphere) in comparison to for a C/O ratio of 1.0, which probes deeper atmosphere where it has

an inversion and leads to HCN emission. Interestingly, at 1.6 µm HCN leads to an emission feature both at a C/O ratio of 1 and 1.5, since at 1.6 µm, the inversion layer is being probed at both C/O ratios as seen in the NCF in Fig. 15(c) and *P-T* profile in Fig. 14(a). Therefore, HCN should be observed in emission as well as absorption at high C/O ratios (1.5) for extremely irradiated planets such as WASP-121b.

## 6 INTERPRETING THE OBSERVATIONS OF WASP-121B

WASP-121b is not expected to be cloudy due to its extremely high predicted dayside temperature (*T*<sub>eq</sub> = 2400 K). Moreover, Evans et al. (2018) showed that the extremely steep slope in the optical is not due to enhanced Rayleigh scattering from a haze. In Evans et al. (2018) observations were interpreted using isothermal *P-T* profiles (see figs 11 and 12 in Evans et al. 2018) resulting in detection of H<sub>2</sub>O, evidence of VO and the possibility of SH causing significant
<figure><img src="image_15.png" /><figcaption>**Figure 13.** (a) Figure showing transmission spectra features of each individual molecule used in ATMO (1–10) for WASP-17b transmission spectra shown in Fig. 12(c) at C/O ratio of 1.5. H<sub>2</sub>–H<sub>2</sub> (blue), H<sub>2</sub>–He (green), H<sub>2</sub>O (red), CO<sub>2</sub> (cyan), CO (magenta), CH<sub>4</sub> (yellow), NH<sub>3</sub> (light blue), Na (purple), K (brown), Li (light green), and all 21 opacities (black). The lower continuum boundary in the spectrum sharply dropping around 1.6 µm is due to H<sup>−</sup> opacity. No R<sub>p</sub>/R<sub>⋆</sub> offset was applied while plotting. Individual simulations are divided into blocks of 10 while plotting for clarity. (b) Same as (a) but for Rb (blue), Cs (green), TiO (red), VO (cyan), FeH (magenta), PH<sub>3</sub> (yellow), H<sub>2</sub>S (light blue), HCN (purple), C<sub>2</sub>H<sub>2</sub> (brown), SO<sub>2</sub> (light green), Fe (violet), and all 21 opacities (black).</figcaption></figure>

absorption in the UV and optical. Therefore, here we use our grid of model transmission spectra with RCE *P*-*T* profiles and additional H⁻ and Fe opacity, without grey clouds and enhanced scattering, to fit the WASP-121b data from Evans et al. (2018). We have excluded the data shortward of 0.47 µm while performing the fitting, since the steep slope in this region is not explained by Rayleigh scattering and is likely in part due to absorption of atomic species like Fe in the thermospheric and higher atmospheric layers (Sing et al. 2019; Gibson et al. 2020) which are not modelled here. Fig. 16(a) shows the best-fitting model transmission spectra for observations from Evans et al. (2018). It shows an H₂O feature at 1.4 µm and features in the optical resembling TiO/VO features. The decoupled spectra in Figs 16(c) and (d) show that the optical spectra is dominated by VO features. This fitting result corroborates the result using isothermal *P*-*T* profiles discussed in Evans et al. (2018) and consistent with VO being present in the atmosphere of WASP-121b, since the limb *P*-*T* structure of this planet is in that narrow temperature regime where VO abundance dominates the TiO abundance, as discussed in Goyal et al. (2019b).

The best-fitting transmission spectra model with RCE *P*-*T* profiles gives a reduced χ² value of 1.37 with 74 dof (75 data points minus 1). It gives supersolar metallicity and supersolar C/O ratio as shown in Fig. 16(a) (both when varying O/H and C/H). The ‘best-fitting’ *P*-*T* profile shown in the subplot of Fig. 16(a) does not show a temperature inversion, since the limb *P*-*T* profile (in the observed region) is expected to be cooler than the dayside *P*-*T* profile shown in subplot of Fig. 16(b). Although the 100 times solar metallicity case is the ‘best fit’ (lowest chi-squared), the 50 times solar metallicity case also fits within 3σ and a solar metallicity case is consistent with the data at the 6σ level, as shown in Fig. B1 of Appendix B, which shows χ² maps for all the model simulations in the grid for WASP-121b with confidence intervals, similar to those described in details in Goyal et al. (2018). There are two reasons why high metallicities are preferred, one is the increase in the VO abundances as metallicity increases (more than that of TiO) and second the decrease in the scale height compared to solar metallicity case, both required to better fit the observations, driven by the low R<sub>p</sub>/R<sub>s</sub> points in the

optical. We also performed tests where we excluded the optical data (< 1 µm) and just used the near-infrared data while fitting and in this case the subsolar metallicity model simulation (0.1 times solar) was the ‘best fit’, thus again concluding that the optical data lead to a higher metallicity being preferred. We note that we have not included relative levels between STIS and WFC3 as a free parameter while fitting, which can potentially alter the results.

These results demonstrate that a reasonably good match to the observations can be achieved with self-consistent 1D models assuming chemical equilibrium for a range of plausible metallicities and C/O values. The reality is undoubtedly significantly more complex. In particular, we have assumed a single RCE *P*-*T* profile for the entire limb of the planet, but there should be a strong temperature contrast between the east and west limbs, in turn affecting which species are in the gas phase or condensed. However, limb to limb variations can be incorporated while fitting, simply by the linear combination of transmission spectra from our library for a specific planet (MacDonald, Goyal & Lewis 2020). The ‘best-fitting’ limb *P*-*T* profile of WASP-121b is cold enough for Fe to rainout from the atmosphere, at least in the region of the atmosphere where local thermodynamic equilibrium (LTE) can be assumed. However, Sing et al. (2019) and Gibson et al. (2020) detected the presence of Fe II and Fe I in UV and optical high-resolution observations of WASP-121b, which probes much higher altitudes than the optical transmission spectrum and may well be in the non-LTE region of the atmosphere. The deep interior temperature may also be hotter (Thorngren et al. 2019; Sing et al. 2019) than assumed in our work (T<sub>int</sub>=100 K), which can prevent Fe from raining out of the atmosphere. The other possibility is that the true planet limb *P*-*T* profile is much hotter than our ‘best-fitting’ limb *P*-*T* profile, since we have not modelled the transmission spectrum of WASP-121b shortward of 0.47 µm for this fitting and not considered interterminator variations, as detailed earlier.

Fig. 16(b) shows the best-fitting model emission spectra for observations from Mikal-Evans et al. (2020), with a χ² value of 88.63 and reduced χ² value of 1.81 with 49 dof. The best-fitting model has solar metallicity and C/O ratio of 0.75, with *f*<sub>c</sub> value of 1.0, indicative
<figure><img src="image_16.png" /><figcaption>**Figure 14.** (a) Figure showing *P-T* profiles for a range of C/O ratios at 0.5 *f*<sub>c</sub> and solar metallicity for WASP-121b. (b) Figure showing equilibrium chemical abundances for some important species for various C/O values obtained using *P-T* profiles shown in (a). (c) Figure showing transmission spectra for WASP-121b for different values of C/O ratios obtained using *P-T* profiles shown in (a) and corresponding equilibrium chemical abundances shown in (b). (d) Figure showing emission spectra for WASP-121b for different values of C/O ratios obtained using *P-T* profiles shown in (a) and corresponding equilibrium chemical abundances shown in (b).</figcaption></figure>

of very low energy redistribution in the atmosphere. The dayside best-fitting *P-T* profile in the subplot of Fig. 16(b) shows a temperature inversion, while the limb *P-T* profile obtained using the transmission spectra fit (see Fig. 16a) is cooler without a temperature inversion, again supporting very low energy redistribution in the atmosphere of WASP-121b. It shows a 1.4 µm H<sub>2</sub>O emission feature, indicative of inversion layer, as explained in previous sections.

## 7 CONCLUSIONS

The approximation of an isothermal *P-T* profile has many limitations for a real planetary atmosphere and is a reasonable approximation only to generate model transmission spectra, which probes a very small part of the overall atmosphere. Therefore, in this work we create a library of RCE *P-T* profiles and corresponding equilibrium chemical abundances, transmission spectra, emission spectra, and CF, for a range of observationally significant exoplanets. Unlike a static table of chemical abundances used in various models, these *P-T* profiles (and all models in this grid) are consistent with equilibrium

chemical abundances, meaning *P-T* profile, and chemical abundances are tied together during each iteration, making it a state of the art model grid. The library of models presented in this work also included H<sup>-</sup> and Fe opacity, in addition to opacities included in previous library of models presented in Goyal et al. (2018, 2019b). The major conclusions from this work are divided into three parts, first the sensitivity to model choices, second the sensitivity to grid parameters, and third the general conclusions.

### 7.1 Sensitivity to model choices

We have to make different model choices while interpreting the observations of exoplanet atmospheres, such as the type of condensation, opacities to be included, source of line-list, line-broadening profile, etc. Here, we conclude on the effect of these different model choices.

(i) The investigation of differences in RCE model simulations with local and rainout condensation revealed that adopting different condensation approaches can result in substantial differences in the
<figure><img src="image_17.png" /><figcaption>**Figure 15.** (a) Same as Fig. 13(a), but for WASP-121b at 0.5 *f*<sub>c</sub>, solar metallicity, and C/O ratio of 1.5. H<sup>-</sup> opacity is quite dominant especially between 1 and 1.4 µm. (b) Same as Fig. 13(b), but for WASP-121b at 0.5 *f*<sub>c</sub>, solar metallicity and C/O ratio of 1.5. (c) and (d) Figures showing NCF at 1.6 and 3.1 µm for a range of C/O ratios for emission spectra as shown in Fig. 14(d).</figcaption></figure>

resulting *P*-*T* profiles and thereby the model spectra, for a typical hot Jupiter planet like WASP-17b.

(ii) Inclusion of thermal ionic species reduces the abundance of Na and K at atmospheric temperatures greater than ~2000 K and low pressures (~10<sup>-1</sup> bar), as they are ionized to form Na<sup>+</sup> and K<sup>+</sup>. Therefore, Na and K features are absent in the transmission spectrum at these temperatures. It also allows formation of other ions with significant opacities such as H<sup>-</sup>. Thermal decomposition leads to reduction in the H<sub>2</sub>O abundance for planets such as WASP-121b, but only in the region of the atmosphere probed by transmission spectrum and not the emission spectrum.

(iii) Iron (Fe) opacity leads to very sharp inversion at the top of the model atmosphere (~10<sup>-5</sup> to 10<sup>-6</sup> bar). However, it is unable to produce a temperature inversion in the region of the atmosphere being probed by emission spectrum observations even at 200× solar metallicity, as produced by TiO/VO or H<sup>-</sup> opacities for WASP-121b. Moreover, Fe opacity tends to cool the *P*-*T* profile for *P* > 0.1 bar, due to absorption of radiation in the upper atmosphere. It also imparts narrow spectral features in the UV-optical part of the transmission spectrum shortward of ~0.4 µm. However, Fe opacity

is sufficient to create a strong thermal inversion in KELT-9b, the planet with the highest equilibrium temperature in this library of models, even at solar metallicity and without TiO/VO/H<sup>-</sup> opacities. Thus, the capability of Fe opacity to form a temperature inversion is system specific showing a strong dependance on the host star effective temperature.

(iv) H<sup>-</sup> opacity plays an important role in shaping the *P*-*T* profile, transmission, and emission spectrum, at temperatures greater than ~2200 K. In the presence of TiO/VO in the atmosphere H<sup>-</sup> opacity leads to cooling in the deeper atmosphere (*P* > 0.1 bar). At C/O ratios ≥ 1, H<sup>-</sup> opacity leads to formation of a strong thermal inversion layer, but not at solar C/O ratio even if TiO/VO is absent. Thermal inversion due to H<sup>-</sup> opacity lies deeper (higher pressure) in the atmosphere (~1 to 10<sup>-2</sup> bar), in comparison to the inversion due to TiO/VO opacities (~10<sup>-1</sup> to 10<sup>-3</sup> bar). In the transmission spectrum, H<sup>-</sup> opacity tends to mute the line wings of Na, K, and all features up-to 1.6 µm (limit of H<sup>-</sup> bound free absorption) mimicking the effect of grey cloud deck opacity. The effect of H<sup>-</sup> opacity on *P*-*T* profiles also substantially alters the emission spectrum, with increase in peak planet flux due to thermal inversion.
<figure><img src="image_18.png" /><figcaption>**Figure 16.** (a) Figure showing best-fitting model transmission spectra using the grid of model transmission spectra for WASP-121b presented in this paper and observations from Evans et al. (2018) giving $\chi^2$ value of 101.69 with $f_c = 0.25$, 100 times solar metallicity and C/O ratio of 0.75. Best-fitting *P-T* profile is shown in the subplot. (b) Figure showing best-fitting model emission spectra using the grid of model emission spectra for WASP-121b presented in this paper and observations from Mikal-Evans et al. (2020) with $\chi^2$ value of 88.63 and $f_c = 1.0$, solar metallicity, and C/O ratio of 0.75. Best-fitting *P-T* profile is shown in the subplot. (c) Best-fitting decoupled transmission spectra similar to Fig. 13(a). The lower continuum boundary in the spectrum sharply dropping around 1.6 µm is due to H$^{-}$ opacity. (d) Best-fitting decoupled transmission spectra similar to Fig. 13(b).</figcaption></figure>

Adopting different line-wing profiles of Na and K, led to differences of ~50–100 K in the deeper atmosphere temperature as shown in the online Supporting Information, thus leading to some differences in the emission spectra, which may be detectable with future observations. The differences in the derived *P-T* profiles and the spectra between the simulations, varying C/O ratio by varying O/H and those by varying C/H, shown in the online Supporting Information, revealed some differences within the range of C/O ratio adopted in the library, but to a lesser extent compared to other parameters in the library and other model choices.

## 7.2 Sensitivity to grid parameters

The library of simulations presented in this work was developed for four different recirculation factors, six metallicities, and six C/O ratios for 89 hot Jupiter and warm Neptune exoplanets with equilibrium temperatures greater than 1200 K. The effect of varying the recirculation factor, metallicity, and C/O ratio were investigated and shown for WASP-17b and WASP-121b in Sections 5.3, showing the variation of the *P-T* profile, equilibrium chemical abundances,

and thereby the transmission and the emission spectrum. The major conclusion from these sensitivity tests for different grid parameters are:

(i) The increase in recirculation factor (lower redistribution of energy) leads to transitioning of the *P-T* profile from without a temperature inversion to isothermal to a profile with temperature inversion for WASP-121b, which is also apparent in the emission and transmission spectrum. Therefore, emission spectrum and to a lesser extent transmission spectrum (due to other degeneracies) can be used to constrain the energy redistribution for all the planets presented in this grid.

(ii) The increase in metallicity leads to the *P-T* profile becoming hotter for both WASP-17b and WASP-121, due to the increase in the total opacity, driven by the increase in the mean molecular weight of the atmosphere. The increase in metallicity leads to an increase in the size of the transmission spectral features at first (up to 10× solar metallicity) due to higher abundances and higher temperatures, but then they decrease again as the scale height of the atmosphere decreases (due to increased mean molecular weight).
with increased metallicity. However, this effect is not seen in the emission spectrum where absorption/emission features increase with increasing metallicity, since emission spectral features do not have scale height dependance as transmission spectral features. The shift in the atmospheric level probed by emission spectrum to lower pressure levels with increasing metallicity, as shown by CF in the H₂O and CO₂ bands, as well as the difference in the level probed by the wings and cores of absorption/emission bands for a given metallicity, both can be used to constrain metallicity and the *P-T* profile of the atmosphere.

(iii) Increasing the C/O ratio leads to a cooling of *P-T* profile for WASP-17b, with a transition seen from a H₂O dominated spectrum to a spectrum dominated by carbon bearing species (primarily CH₄) between C/O ratios of 0.75 and 1. For WASP-121b, the strength of the temperature inversion decreases as the abundance of TiO/VO decreases with increasing C/O ratio. At a C/O ratio greater than 1 even though TiO/VO abundances are low, a weak temperature inversion is maintained in the atmosphere due to Na and K at *f*<sub>c</sub> value of 0.5 as shown by previous works. Whereas, for an *f*<sub>c</sub> value of 1 (hotter *P-T* profile) H⁻ opacity contributes substantially to the development of an inversion. Interestingly, extremely irradiated exoplanet atmospheres like WASP-121b show a spectrum dominated by HCN features at infrared wavelengths for high C/O ratios (>1) in comparison to comparatively cooler planets like WASP-17b which are dominated by CH₄ features. HCN and H⁻ features can be seen as an emission as well as absorption feature at different wavelengths in their emission spectrum. Therefore, HCN and H⁻ features for these planets can be used to constrain their C/O ratios as well as their thermal structure.

## 7.3 General conclusions

The interpretation of observations using the library of models with RCE *P-T* profiles presented in this work was demonstrated using WASP-121b observations. The transmission spectra probing the planetary atmospheric limb reveal H₂O features in the infrared with evidence of VO in the optical. The best-fitting model indicates supersolar metallicity driven by optical observations and VO opacity, with non-inverted limb *P-T* profile. The emission spectra probing the dayside of the planetary atmosphere reveals an H₂O feature in emission giving evidence of a potential temperature inversion and very low energy redistribution in the atmosphere.

One of the major limitations of the model simulations presented in this paper is that they are 1D. Ideally 2D or 3D approaches would more realistically represent a planetary atmosphere. However, 1D models allow us to explore a very large parameter space in chemistry and radiative transfer for a large number of planetary systems, not possible currently using 2D or 3D approaches due to computational limitations. Moreover, the combination of different *P-T* profiles, transmission and emission spectra for a given planet can be used to represent different parts of the atmosphere (see for e.g. MacDonald et al. 2020), thus paving the way to interpret 2D or 3D observations, like for example using spectral phase curves. We have also not modelled disequilibrium processes such as photochemistry and vertical mixing in this library of simulated atmospheres, which is planned for the future development. The current model spectra in the library has a resolution of *R* ~ 5000 at 0.2 µm decreasing to *R* ~ 35 at 30 µm, which cannot be used to interpret high-resolution observations. However, the model atmospheric RCE *P-T* profiles and the corresponding chemical abundances from our library can be used to generate synthetic high-resolution spectra (transmission or emission), which can then be used to interpret high-resolution

observations. We also plan to update the spectra in our library online, to achieve the highest possible *JWST* resolution over some specific wavelength ranges.

The library of models presented in this work can be used to plan and interpret a broad range of transiting exoplanet observations with current and future facilities such as *HST* and *JWST*, since it covers the wavelength range from 0.2 to 30 µm. It represents a valuable resource to guide choices of observational targets and instrument modes. The library of models can also be used to initialize retrieval models as well as to constrain priors in them. Our development of an expansive grid of thermochemically consistent 1D models for exoplanet atmospheres provides the community with another tool with which to tackle our theoretical understanding of hot Jupiter atmospheres and interpret current and future observations of these planets.

## ACKNOWLEDGEMENTS

We would like to thank the anonymous referee for their detailed and constructive comments, that improved the clarity and quality of this manuscript. JMG and NM acknowledge funding by a Leverhulme Trust Research Project Grant and University of Exeter College of Engineering, Mathematics and Physical Sciences PhD studentship. MWP acknowledges support through a UK Research and Innovation-Science and Technology Facilities Council (UKRISTFC) studentship. This work used the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility. This work also used the University of Exeter Supercomputer, a DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS, and the University of Exeter.

## DATA AVAILABILITY

The data underlying this article are currently available publicly here https://drive.google.com/drive/folders/1zCCe6HICuK2nLgnYJFa17W4lyunjU4JE and here https://noctis.erc-atmo.eu:5001/fsdow nload/hq0z4udQJ/goyal2020.

## REFERENCES

* Alam M. K. et al., 2018, AJ, 156, 298

* Allard N. F., Allard F., Hauschildt P. H., Kielkopf J. F., Machin L., 2003, A&A, 411, L473

* Allard F., Homeier D., Freytag B., 2012, Philos. Trans. R. Soc. Lond. Ser. A, 370, 2765

* Allard N.-F., Spiegelman F., Kielkopf J.-F., 2007, A&A, 465, 1085

* Amundsen D. S., 2015, PhD thesis, University of Exeter

* Amundsen D. S., Baraffe I., Tremblin P., Manners J., Hayek W., Mayne N. J., Acreman D. M., 2014, A&A, 564, A59

* Arcangeli J. et al., 2018, ApJ, 855, L30

* Baraffe I., Homeier D., Allard F., Chabrier G., 2015, A&A, 577, A42

* Baudino J.-L., Mollière P., Venot O., Tremblin P., Bézard B., Lagage P.-O., 2017, ApJ, 850, 150

* Bell K. L., Berrington K. A., 1987, J. Phys. B At. Mol. Phys., 20, 801

* Burgasser A. J., Kirkpatrick J. D., Liebert J., Burrows A., 2003, ApJ, 594, 510

* Burrows A. et al., 1997, ApJ, 491, 856

* Burrows A., Sharp C. M., 1999, ApJ, 512, 843

* Burrows A., Marley M. S., Sharp C. M., 2000, ApJ, 531, 438

* Burrows A., Ram R. S., Bernath P., Sharp C. M., Milsom J. A., 2002, ApJ, 577, 986

* Burrows A., Budaj J., Hubeny I., 2008, ApJ, 678, 1436
Carter A. L. et al., 2020, MNRAS, 494, 5449

Chamberlain J. W., Hunten D. M., 1987, Theory of Planetary Atmospheres. An Introduction to their Physics and Chemistry, Vol. 36, Academic Press Inc., Orlando, FL, USA

Chubb K. L., Min M., Kawashima Y., Helling C., Waldmann I., 2020, A&A, 639, A3

Drummond B., Tremblin P., Baraffe I., Amundsen D. S., Mayne N. J., Venot O., Goyal J., 2016, A&A, 594, A69

Drummond B., Mayne N. J., Manners J., Baraffe I., Goyal J., Tremblin P., Sing D. K., Kohary K., 2018, ApJ, 869, 28

Drummond B., Carter A. L., Hébrard E., Mayne N. J., Sing D. K., Evans T. M., Goyal J., 2019, MNRAS, 486, 1123

Evans T. M. et al., 2017, Nature, 548, 58

Evans T. M. et al., 2018, AJ, 156, 283

Fortney J. J. et al., 2016, preprint (arXiv:1602.06305)

Fortney J. J., 2005, MNRAS, 364, 649

Fortney J. J., Marley M. S., 2007, ApJ, 666, L45

Fortney J. J., Marley M. S., Barnes J. W., 2007, ApJ, 659, 1661

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

Gandhi S., Madhusudhan N., 2017, MNRAS, 472, 2334

Gandhi S., Madhusudhan N., 2019, MNRAS, 485, 5817

Gibson N. P. et al., 2020, MNRAS, 493, 2215

Goyal J. M., 2019, PhD thesis, University of Exeter

Goyal J. M. et al., 2018, MNRAS, 474, 5158

Goyal J. M. et al., 2019a, MNRAS, 486, 783

Goyal J. M., Wakeford H. R., Mayne N. J., Lewis N. K., Drummond B., Sing D. K., 2019b, MNRAS, 482, 4503

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

Heiter U. et al., 2008, J. Phys. Conf. Ser, 130, 012011

Heiter U. et al., 2015, Phys. Scr., 90, 054010

Heng K., Kitzmann D., 2017, MNRAS, 470, 2972

Hobbs R., Shorttle O., Madhusudhan N., Rimmer P., 2019, MNRAS, 487, 2242

Hubeny I., Burrows A., Sudarsky D., 2003, ApJ, 594, 1011

John T. L., 1988, A&A, 193, 189

Knutson H. A. et al., 2009, ApJ, 690, 822

Kreidberg L. et al., 2018, AJ, 156, 17

Lecavelier Des Etangs A., Pont F., Vidal-Madjar A., Sing D., 2008, A&A, 481, L83

Lodders K., Fegley B., 2002, Icarus, 155, 393

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

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

Madhusudhan N., 2012, ApJ, 758, 36

Madhusudhan N., Agúndez M., Moses J. I., Hu Y., 2016, Space Sci. Rev., 205, 285

Malik M., Kitzmann D., Mendonça J. M., Grimm S. L., Marleau G.-D., Linder E. F., Tsai S.-M., Heng K., 2019, AJ, 157, 170

Mansfield M. et al., 2018, AJ, 156, 10

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

Marley M. S., Saumon D., Guillot T., Freedman R. S., Hubbard W. B., Burrows A., Lunine J. I., 1996, Science, 272, 1919

McKemmish L.-K., Yurchenko S.-N., Tennyson J., 2016, MNRAS, 463, 771

Mikal-Evans T. et al., 2019, MNRAS, 488, 2222

Mikal-Evans T. et al., 2020, MNRAS, 496, 1638

Molaverdikhani K., Henning T., Mollière P., 2019, ApJ, 873, 32

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

Mollière P., van Boekel R., Bouwman J., Henning T., Lagage P.-O., Min M., 2016, A&A, 600, A10

Moses J. I., Madhusudhan N., Visscher C., Freedman R. S., 2013, ApJ, 763, 25

Nikolov N. et al., 2018, Nature, 557, 526

Parmentier V. et al., 2018, A&A, 617, A110

Plez B., 1999, in Le Bertre T., Lebre A., Waelkens C., eds, Asymptotic Giant Branch Stars, IAU Symposium, p. 75

Rajpurohit A. S., Reylé C., Allard F., Homeier D., Schultheis M., Bessell M. S., Robin A. C., 2013, A&A, 556, A15

Rau A. R. P., 1996, J. Astrophys. Astron., 17, 113

Rybicki G. B., Lightman A. P., 1986, Radiative Processes in Astrophysics, Wiley-VCH, New York, p. 400

Sauval A. J., Tatum J. B., 1984, ApJS, 56, 193

Seager S., Sasselov D. D., 1998, ApJ, 502, L157

Sharp C. M., Burrows A., 2007, ApJS, 168, 140

Sing D. K. et al., 2016, Nature, 529, 59

Sing D. K. et al., 2019, AJ, 158, 91

Southworth J., 2011, MNRAS, 417, 2166

Spiegel D. S., Silverio K., Burrows A., 2009, ApJ, 699, 1487

Sudarsky D., Burrows A., Hubeny I., 2003, ApJ, 588, 1121

Thorngren D., Gao P., Fortney J. J., 2019, ApJ, 884, L6

Tremblin P. et al., 2017, ApJ, 841, 30

Tremblin P., Amundsen D. S., Mourier P., Baraffe I., Chabrier G., Drummond B., Homeier D., Venot O., 2015, ApJ, 804, L17

Tremblin P., Amundsen D. S., Chabrier G., Baraffe I., Drummond B., Hinkley S., Mourier P., Venot O., 2016, ApJ, 817, L19

Visscher C., Lodders K., Fegley B., Jr, 2006, ApJ, 648, 1181

von Essen C., Mallonn M., Welbanks L., Madhusudhan N., Pinhas A., Bouy H., Weis Hansen P., 2019, A&A, 622, A71

Wakeford H. R. et al., 2018, AJ, 155, 29

Woitke P., Helling C., Hunter G. H., Millard J. D., Turner G. E., Worters M., Blecic J., Stock J. W., 2018, A&A, 614, A1

Zahnle K., Schaefer L., Fegley B., 2010, Cold Spring Harb. Perspect. Biol., 2, a004895

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

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

**Figure S1.** Figure showing WASP-17b emission spectrum with 0.25 $f_c$, solar metallicity, and solar C/O ratio, when removing opacity contributions of single species at a time from the calculation, shown by their respective colours.

**Figure S2.** (a) Same as Fig. 1(a), but for WASP-121b, at solar metallicity and solar C/O ratio with 0.5 $f_c$.

**Figure S3.** (a) Figure showing absorption cross-sections ($cm^2$ molecule$^{-1}$) of Na and K using two different pressure broadening treatments from Allard et al. (2007) and Burrows et al. (2000). (b) Figure showing RCE P-T profiles using two different pressure broadening profiles (Allard and Burrows) for Na and K, shown in Fig. 3(a), for WASP-17b at solar metallicity and solar C/O ratio with 0.25 $f_c$.

**Figure S4.** (a) Figure showing absorption cross-sections ($cm^2$ molecule$^{-1}$) of VO using different two different line-lists VOMYT (McKemmish et al. 2016) and PLEZ (Plez 1999) (b) Figure showing RCE P-T profiles using different VO line-list named PLEZ and VOMYT for WASP-121b at 200 × solar metallicity and solar C/O ratio with 0.5 $f_c$.

**Figure S5.** (a) Figure showing P-T profiles with C/O ratio of 0.35 by varying C/H (red) and O/H (blue) relative to solar C/O ratio (0.55) at 0.25 $f_c$ and solar metallicity (0.56) for WASP-17b.

**Figure S6.** Figure showing best-fitting model transmission spectra using the grid of model transmission spectra for WASP-17b presented in this paper and observations from Sing et al. (2016) giving $\chi^2$ value of 64.64 (reduced $\chi^2$ = 1.7) with $f_c$ = 0.25, solar metallicity, and C/O ratio of 0.35 (subsolar).

**Table S1.** Line-list (opacity) sources of different species used in this library. Pressure broadening sources are shown in Table S2.

**Table S2.** Type and source of pressure broadening for all opacities used in this library.
**Table S3.** All the stellar and planetary parameters adopted from TEPCAT (Southworth 2011a) data base, for the model simulations of 89 exoplanets in the grid are listed here.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

# APPENDIX A: HIGH TEMPERATURE ADDITIONS

## A1 Thermal ionization

Thermal ionic species were recently included in ATMO to compute the chemical equilibrium abundance of H⁺, H⁻, Na⁺, K⁺, Ca⁺, Si⁺, and e⁻. This was simply done by including thermodynamic data of these species in our database and using Gibbs energy minimization to compute the equilibrium abundances in each layer of the atmosphere, including the ionic species. The equilibrium chemical abundances of ionic species from ATMO were validated by comparing with those from the GGchem model (Woitke et al. 2018) at 2000, 3000, and 4000 K as shown in Fig. A1 for 3000 K. We used the same input P-T profile for both ATMO and GGchem for this validation. It can be seen that there is a very good agreement between the abundances from both the models. Only small differences exist in the deeper parts of the atmosphere, currently not probed by either emission or transmission spectra.

## A2 Implementation and validation of H⁻ opacity

Many publications highlight the importance of H⁻ opacity in the atmosphere of extremely irradiated hot Jupiters (e.g Arcangeli et al. 2018; Mansfield et al. 2018; Kreidberg et al. 2018; Parmentier et al. 2018). We also included H⁻ opacities in ATMO and investigated its effects on the P-T profiles and thereby the spectra. The formation of species H⁻ (hydrogen anion) is basically a result of electron attachment, driven by the presence of abundant hydrogen (H) and low energy electrons in the ionized atmospheres of stars or hot Jupiters (Rau 1996). The absorption of electromagnetic radiation by H⁻ is driven by photo-detachment (bound-free) and free-free transition processes. These are computed using the analytical equations of John

<figure><img src="image_19.png" /><figcaption>**Figure A1.** Figure showing validation of ATMO chemical equilibrium abundances of ionic species with GGChem (Woitke et al. 2018) at 3000 K.</figcaption></figure>

(1988) derived from the original derivation of free-free transition in Bell & Berrington (1987).

The photo-detachment process of H⁻ absorption is given by

$$h\nu + H^{-} \rightarrow H + e^{-}, \qquad (A1)$$

where h is the Planck's constant and v is the frequency of radiation. In this process, H⁻ ions absorb radiation of frequency v, to form atomic hydrogen and an electron. As shown in John (1988) for wavelengths (λ) less than the ionization threshold λ of H⁻ (λ < 1.6419 µm) this is computed using,

$$k_{\text{bf}}(\lambda) = 10^{-18} \lambda^3 \left(\frac{1}{\lambda} - \frac{1}{\lambda_0}\right)^{3/2} f(\lambda),$$

where,

$$f(\lambda) = \sum_{n=1}^{6} C_n \left[ \frac{1}{\lambda} - \frac{1}{\lambda_0} \right]^{(n-1)/2}, \qquad (A2)$$

where $k_{\text{bf}}(\lambda)$ is the bound-free cross-section of H⁻ in units of cm², λ is the wavelength, $\lambda_0 = 1.6419$ µm is the threshold wavelength, and $C_n$ are the coefficients for n different values given in table 2 of John (1988). The bound-free absorption cross-section above the threshold λ of 1.6419 µm is zero. The total opacity due to bound-free absorption is then computed using,

$$\kappa_{\text{bf}} = k_{\text{bf}}(\lambda, T) A[H^{-}]^{n_d} \rho, \qquad (A3)$$

where $\kappa_{\text{bf}}$ is the total bound-free opacity, $A[H^{-}]$ is the abundance of H⁻ (mixing ratio), $n_d$ is the atmospheric number density (cm⁻³), and $\rho$ is the atmospheric mass density (g cm⁻³).

The free-free transition process of H⁻ is given by

$$h\nu + e^{-} + H^{-} \rightarrow H + e^{-}. \qquad (A4)$$

In this reaction photons can be absorbed by electrons interacting with neutral hydrogen atom across the whole spectral range (0 < λ < ∞). This process is solely responsible for H⁻ opacity beyond the ionization threshold wavelength (1.6419 µm). This is computed using

$$k_{\text{ff}}(\lambda, T) = 10^{-29} \sum_{n=1}^{6} \left(\frac{5040}{T}\right)^{(n+1)/2} \times (\lambda^2 A_n + B_n + C_n/\lambda + D_n/\lambda^2 + E_n/\lambda^3 + F_n/\lambda^4), \qquad (A5)$$

where $k_{\text{ff}}(\lambda, T)$ is the free-free cross-section of H⁻ in units of cm⁴ dyne⁻¹ and T is the temperature. $A_n, B_n, C_n, D_n, E_n$, and $F_n$ are coefficients as given in tables 3a and b of John (1988) for λ > 0.3645 µm and 0.1823 < λ < 0.3645 µm, respectively. By multiplying equation (A5) by the Boltzmann's constant (1.38 × 10⁻¹⁶ erg s⁻¹) and the temperature, $k_{\text{ff}}(\lambda, T)$ is obtained in the units of cm⁵. The total opacity due to free-free absorption is then computed using,

$$\kappa_{\text{ff}} = k_{\text{ff}}(\lambda, T) A[H^{-}] \frac{n_d}{\rho} A[e^{-}] n_d, \qquad (A6)$$

where $\kappa_{\text{ff}}$ is the total free-free opacity, $A[H^{-}]$ is the abundance of neutral hydrogen (mixing ratio), and $A[e^{-}]$ is the abundance of electron (mixing ratio). The total opacity of the H⁻ ion ($\kappa_{\text{tot}}$) due to bound-free and free-free transitions is then given by $\kappa_{\text{tot}} = \kappa_{\text{bf}} + \kappa_{\text{ff}}$.

We validate the abundance weighted H⁻ opacity by comparing with the results from Parmentier et al. (2018) and Mansfield et al. (2018). Fig. 4 of Parmentier et al. (2018) shows the abundance
<figure><img src="image_20.png" /><figcaption>**Figure A2.** (a) Figure showing abundance weighted cross-section of H⁻ at 0.042 bar and 3100 K from fig. 4 of Parmentier et al. (2018, green) and from ATMO (blue). (b) Figure showing abundance weighted cross-section of H⁻ at 0.084 bar and 2756 K from fig. 7 of Mansfield et al. (2018, green) and from ATMO (blue).</figcaption></figure>

weighted absorption cross-section of H⁻ (cm² molecule⁻¹) as a function of wavelength at the P-T point of 0.042 bar and 3100 K. We also compute the abundance weighted absorption cross-section of H⁻ at this P-T point within ATMO using $κ_{tot} × \frac{\mu_{mean}}{N_A}$, where $κ_{tot}$ is the total opacity of the H⁻ ion computed as shown in the previous section, $\mu_{mean}$ is the mean molecular weight of the atmosphere, and $N_A$ is the Avogadro's constant. The comparison is shown in Fig. A2(a), demonstrating that the agreement is good and even the equilibrium chemical abundances at this P-T point are similar in ATMO and Parmentier et al. (2018, from fig. 3 in their paper), thus validating the implementation of H⁻ opacity in ATMO.

Fig. 7 of Mansfield et al. (2018) shows the abundance weighted absorption cross-section of H⁻ at a P-T point of 0.084 bar and 2756 K. When compared with this, there is a substantial difference between the abundance weighted absorption cross-section of H⁻ as shown in Fig. A2(a). The primary reason being the difference in equilibrium chemical abundances at this P-T point, which is ~3.4 times larger in ATMO as compared to Mansfield et al. (2018). However, when this factor of 3.4 is taken into consideration while comparing the abundance weighted absorption cross-section of H⁻ from both models, there is a good agreement (not shown here). The reason for differences in equilibrium chemical abundances is still unclear and can be due to many factors, such as the differences in input elemental abundances, polynomial coefficients, etc., as shown in Goyal et al. (2019a). We note that the equilibrium chemistry scheme used in ATMO has been validated by comparing to various numerical and analytical equilibrium chemistry models, with local and rainout condensation (Drummond et al. 2016; Goyal et al. 2019a).

## A3 Fe absorption cross-sections

Absorption cross-sections of gaseous iron (Fe) have been included in ATMO using the Fe line list from the Vienna Atomic Line Database (VALD) data base (Heiter et al. 2008, 2015) and the

partition function from Sauval & Tatum (1984). The H₂ and He pressure broadening line-widths for Fe have been computed using the van der Waals coefficient contained within the VALD data base and using equation (23) in Sharp & Burrows (2007). The absorption cross-sections of Fe included in ATMO have also been validated by comparing with the absorption cross-sections of Sharp & Burrows (2007).

## APPENDIX B: WASP-121B χ² MAP

<figure><img src="image_21.png" /><figcaption>**Figure B1.** Figure showing χ² map for WASP-121b for the transmission spectra fitting described in Section 6. Contours of χ² are shown for all the combinations of grid parameters. Metallicity is also log-scaled, 0 being solar metallicity, and 2 being 100 times solar metallicity. Colours indicate confidence intervals as shown in colour map to the right. See Goyal et al. (2018) for more details.</figcaption></figure>

This paper has been typeset from a TEX/LATEX file prepared by the author.