# A spectroscopic thermometer: individual vibrational band spectroscopy with
the example of OH in the atmosphere of WASP-33b

Sam O.M. Wright Department of Physics and Astronomy, University College
London, Gower Street, WC1E 6BT London, United Kingdom<EMAIL_ADDRESS>Stevanus K. Nugroho Astrobiology Center, NINS, 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan  National Astronomical Observatory of Japan, NINS, 2-21-1
Osawa, Mitaka, Tokyo 181-8588, Japan  Matteo Brogi Dipartimento di Fisica,
Università degli Studi di Torino, via Pietro Giuria 1, I-10125, Torino, Italy
Department of Physics, University of Warwick, Coventry CV4 7AL, UK  INAF-
Osservatorio Astrofisico di Torino, Via Osservatorio 20, I-10025 Pino
Torinese, Italy  Neale P. Gibson School of Physics, Trinity College Dublin,
The University of Dublin, Dublin 2, Ireland  Ernst J. W. de Mooij Astrophysics
Research Centre, Queen’s University Belfast, Belfast BT7 1NN, UK  Ingo
Waldmann Department of Physics and Astronomy, University College London, Gower
Street, WC1E 6BT London, United Kingdom  Jonathan Tennyson Department of
Physics and Astronomy, University College London, Gower Street, WC1E 6BT
London, United Kingdom  Hajime Kawahara Department of of Space Astronomy and
Astrophysics, Institute of Space and Astronautical Science, Japan Aerospace
Exploration Agency, 3-1-1, Yoshinodai, Chuo-ku, Sagamihara, Kanagawa,
252-5210, Japan  Masayuki Kuzuhara Astrobiology Center, NINS, 2-21-1 Osawa,
Mitaka, Tokyo 181-8588, Japan  National Astronomical Observatory of Japan,
NINS, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  Teruyuki Hirano
Astrobiology Center, NINS, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
National Astronomical Observatory of Japan, NINS, 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan  Department of Astronomy, School of Science, The Graduate
University for Advanced Studies (SOKENDAI), 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan  Takayuki Kotani Astrobiology Center, NINS, 2-21-1 Osawa,
Mitaka, Tokyo 181-8588, Japan  National Astronomical Observatory of Japan,
NINS, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  Department of Astronomy,
School of Science, The Graduate University for Advanced Studies (SOKENDAI),
2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  Yui Kawashima Cluster for
Pioneering Research, RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan  Kento
Masuda Department of Earth and Space Science, Osaka University, Osaka
560-0043, Japan  Jayne L. Birkby Astrophysics, Department of Physics,
University of Oxford, Keble Road, Oxford OX1 3RH, UK  Chris A. Watson School
of Mathematics and Physics, Queen’s University Belfast, University Road,
Belfast, BT7 1NN, United Kingdom  Motohide Tamura Department of Astronomy,
Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku,
Tokyo 113-0033, Japan  Astrobiology Center, NINS, 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan  National Astronomical Observatory of Japan, NINS, 2-21-1
Osawa, Mitaka, Tokyo 181-8588, Japan  Konstanze Zwintz Institute for Astro-
and Particle Physics, University of Innsbruck, Technikerstrasse 25/8, A-6020
Innsbruck, Austria  Hiroki Harakawa Subaru Telescope, 650 N. Aohoku Place,
Hilo, HI 96720, USA  Tomoyuki Kudo Subaru Telescope, 650 N. Aohoku Place,
Hilo, HI 96720, USA  Klaus Hodapp  University of Hawaii, Institute for
Astronomy, 640 N. Aohoku Place, Hilo, HI 96720, USA  Shane Jacobson
University of Hawaii, Institute for Astronomy, 640 N. Aohoku Place, Hilo, HI
96720, USA  Mihoko Konishi  Faculty of Science and Technology, Oita
University, 700 Dannoharu, Oita 870-1192, Japan  Takashi Kurokawa
Astrobiology Center, NINS, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
Institute of Engineering, Tokyo University of Agriculture and Technology,
2-24-16, Nakacho, Koganei, Tokyo, 184-8588, Japan  Jun Nishikawa National
Astronomical Observatory of Japan, NINS, 2-21-1 Osawa, Mitaka, Tokyo 181-8588,
Japan  Astrobiology Center, NINS, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
Department of Astronomy, School of Science, The Graduate University for
Advanced Studies (SOKENDAI), 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
Masashi Omiya Astrobiology Center, NINS, 2-21-1 Osawa, Mitaka, Tokyo 181-8588,
Japan  National Astronomical Observatory of Japan, NINS, 2-21-1 Osawa, Mitaka,
Tokyo 181-8588, Japan  Takuma Serizawa  Institute of Engineering, Tokyo
University of Agriculture and Technology, 2-24-16, Nakacho, Koganei, Tokyo,
184-8588, Japan  National Astronomical Observatory of Japan, NINS, 2-21-1
Osawa, Mitaka, Tokyo 181-8588, Japan  Akitoshi Ueda  Astrobiology Center,
NINS, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  National Astronomical
Observatory of Japan, NINS, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
Department of Astronomy, School of Science, The Graduate University for
Advanced Studies (SOKENDAI), 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
Sébastien Vievard Astrobiology Center, NINS, 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan  Subaru Telescope, 650 N. Aohoku Place, Hilo, HI 96720, USA
Sergei N. Yurchenko Department of Physics and Astronomy, University College
London, Gower Street, WC1E 6BT London, United Kingdom

###### Abstract

Individual vibrational band spectroscopy presents an opportunity to examine
exoplanet atmospheres in detail by distinguishing where the vibrational state
populations of molecules differ from the current assumption of a Boltzmann
distribution. Here, retrieving vibrational bands of OH in exoplanet
atmospheres is explored using the hot Jupiter WASP-33b as an example. We
simulate low-resolution spectroscopic data for observations with the JWST’s
NIRSpec instrument and use high resolution observational data obtained from
the Subaru InfraRed Doppler instrument (IRD). Vibrational band-specific OH
cross section sets are constructed and used in retrievals on the (simulated)
low and (real) high resolution data. Low resolution observations are simulated
for two WASP-33b emission scenarios: under the assumption of local thermal
equilibrium (LTE) and a toy non-LTE model for vibrational excitation of
selected bands. We show that mixing ratios for individual bands can be
retrieved with sufficient precision to allow the vibrational population
distributions of the forward models to be reconstructed. A simple fit for the
Boltzmann distribution in the LTE case shows that the vibrational temperature
is recoverable in this manner. For high resolution, cross-correlation
applications, we apply the individual vibrational band analysis to an IRD
spectrum of WASP-33b, applying an ‘un-peeling’ technique. Individual detection
significances for the two strongest bands are shown to be in line with
Boltzmann distributed vibrational state populations consistent with the
effective temperature of the WASP-33b atmosphere reported previously. We show
the viability of this approach for analysing the individual vibrational state
populations behind observed and simulated spectra including reconstructing
state population distributions.

Astronomy data modeling (1859); Exoplanet atmospheres (487); Exoplanet
atmospheric composition (2021); Hot Jupiters (753); High resolution
spectroscopy (2096); Near infrared astronomical observations (1093)

##  1 Introduction

The study of exoplanet atmospheres is motivated by a desire to understand the
physical properties of exoplanets and their origins. These properties reflect
a range of planetary environments which can differ drastically from those in
our Solar System and are yet to be fully understood; to gain a better
understanding we must determine their atmospheric composition and thermal
structure. We can also extract valuable information about the origins of these
planets and their formation histories, including information about planetary
migration and providing constraints on the modelling of protoplanetary disks
(Öberg et al., 2011). Beyond enhancing our current understanding,
characterising exoplanets can inform the direction of future research, targets
for follow-up study and the priorities for designing future space missions and
ground-based facilities.

We deduce these characteristics via the process of retrieval; it allows us to
infer parameters related to the atmosphere’s composition and structure from
spectroscopic data acquired by observing the planet. Spectra are produced
based on atmospheric models of physical processes taking place and the
radiative transfer of photons through the planet’s atmosphere. Statistical
techniques (typically Bayesian inference methods) are used to match the
modelled spectra with the observed spectral data given models produced with
the chosen parameterisation. Optimal values for the parameters are found via
these methods, where optimal is quantified as the maximum a posteriori
probability.

Typically the temperature profiles in exoplanetary atmospheric retrievals are
obtained via the fitting of a complex atmospheric model to an observed
absorption (transit) or emission spectrum, where the temperature is one of the
factors shaping the spectral profile. Most notably, temperature affects the
magnitude (intensity) of different absorption/emission bands as well as their
shapes (e.g. Tessenyi et al., 2013). Here we explore an alternative,
‘spectroscopic thermometer’ approach based on the relative intensities of
individual molecular bands as an estimate of their vibrational populations and
therefore the corresponding vibrational temperature. Knowledge of vibrational
state populations obtained via band-by-band retrievals provides an opportunity
to infer the molecular vibrational temperature directly from their spectral
features. Additionally, any outlying band population would indicate the
molecule’s departure from the local thermal equilibrium (LTE) within the
planet’s atmosphere. Especially favourable for this purpose are the (hot)
bands of some molecules which are characterised by substantial displacement
from each other within a relatively small spectroscopic region (Wright et al.,
2022); typical examples are diatomics and linear polyatomics. The ability to
analyse molecular spectra with such granularity will play an increasingly
important role in the characterization of exoplanet atmospheres as accounting
for divergences from local thermodynamic equilibrium and could greatly impact
our understanding of atmospheric conditions, such as temperature structure
(Fossati et al., 2021; Fisher & Heng, 2019) and wind characteristics (Borsa et
al., 2021).

While the approach outlined here naturally enables non-LTE analyses, it also
provides advantages for LTE atmospheric analysis with high and low resolution
observations. For high resolution data it allows temperature measurement
without needing to run computationally expensive T-P profile retrievals. When
analysing low resolution data, this approach can serve as a separate
temperature diagnostic tool with the potential to avoid degeneracies in
temperature retrieval results.

Here we consider the hydroxyl radical (OH) in the atmosphere of WASP-33b both
for low resolution and high resolution applications. OH has recently been
detected in the atmosphere of WASP-33b using the cross correlation technique
(Nugroho et al., 2021) with data from the Subaru telescope’s InfraRed Doppler
instrument (IRD; R $\approx$ 70,000, Tamura et al., 2012; Kotani et al., 2018)
in the spectral region of 0.97–1.75 $\mu$m which was later confirmed by Cont
et al. (2022) using CARMENES data. This near infrared region contains the OH
emission bands formed by the $\Delta v=2$ rotation-vibrational transitions
from vibrational states populated at the relevant temperature of the
atmosphere. The OH spectrum has good wavelength separation between hot
vibrational bands. In this work we aim to use this feature to disentangle
individual band intensities in the atmospheric retrievals in order to target
their individual vibrational populations, i.e. to use information beyond the
aggregated band level information found in the molecular data typically
employed for exoplanet atmospheric analyses. We will show that the separation
between the OH bands in the NIR is sufficient for the bands to be individually
retrievable even with the current observation technology, despite the usual
degeneracy of the atmospheric factors in a typical retrieval of the
atmospheres of hot Jupiters. While spectral signatures have been shown to
differ in non-LTE for other molecules (Wright et al., 2022), it is also worth
investigating OH in particular when considering non-LTE in exoplanet
atmospheres since it has a well documented airglow effect in Earth’s
atmosphere (Noll et al., 2020; Meinel, 1950) and the Martian atmosphere
(Clancy et al., 2013). In non-LTE circumstances, vibrational state populations
diverge substantially from a Boltzmann distribution (as shown in Chen et al.,
2021) which makes the selection of OH particularly relevant for this analysis.

##  2 NIR spectrum of OH and cross section generation

Our goal is to investigate the possibility of retrieving individual
vibrational populations of OH in the NIR emission spectrum of a hot-Jupiter
like atmosphere, using WASP-33b as an example. The MoLLIST line list (Brooke
et al., 2016) for OH is used, along with the ExoCross software (Yurchenko et
al., 2018) for molecular cross section generation. The software provides two
capabilities crucial to this investigation; the means to filter individual
vibrational bands (i.e select isolated bands while still assuming a Boltzmann
distribution under the LTE condition) and also to supply custom state
populations from which cross sections can be calculated.

The NIR spectrum of OH is mainly composed of the $\Delta v=2$ vibrational
bands, including the overtone band $(2,0)$, as well as the weaker $\Delta v=3$
bands, see Fig. 2 where an LTE spectrum of OH for the equilibrium temperature
of $T=3050$ K is shown, with the individual bands ($\Delta v=2$ and $\Delta
v=3$) depicted in different colours. Here, assuming the LTE conditions, cross
sections of 7 (in the first instance) individual vibrational bands ($0\leq
v^{\prime\prime}\leq 6$) of OH were generated using a grid spacing of 0.03
cm-1; this grid spacing corresponds to a minimum resolving power of $R\approx
190,000$ (at 1.75 $\mu$m) which allowed us to retain flexibility at the
modelling stage and the potential for applications of the cross section set to
higher resolution data in the future. The constant grid spacing allows us to
guarantee that we did not introduce artefacts from the grid subdivision
performed by ExoCross to achieve a fixed resolution across the wavelength
range. The vibrational bands were individually isolated via the filtration
method where ExoCross selects only those spectroscopic lines originating from
the desired vibrational state. In line with the LTE approximation, their
magnitudes drop due to the Boltzmann factor as well as the corresponding
vibrational band strength. For reference, we also retained the standard set of
cumulative cross sections of OH for modelling the contribution of OH in LTE to
the spectrum of WASP-33b. Although our goal is to examine the effects in
emission, we remain consistent with atmospheric codes by generating cross
sections for the absorption case. Atmospheric modelling codes (such as TauREx
3 (Al-Refaie et al., 2021)) ingest absorption data and model the emission case
using these absorption data equivalently by Kirchoff’s law under the
assumption of LTE.

In order to be able to retrieve populations of individual bands of OH in
atmospheric spectroscopic models, individual band cross sections were
generated with normalised populations (state population set to 1) as a sum of
all line intensities for transitions $f\leftarrow i$ from a given vibrational
band $(v^{\prime},v^{\prime\prime})$ (Yurchenko et al., 2018)

$\sigma(\tilde{\nu})_{v^{\prime},v^{\prime\prime}}=\sum_{i,f}f_{if}(\tilde{\nu})\frac{A_{if}^{(v^{\prime},v^{\prime\prime})}N_{\rm
r-v}}{8\pi
c^{2}\tilde{\nu}_{if}^{2}}\left(1-e^{-hc\tilde{\nu}_{if}/kT}\right),$ (1)

where $f_{if}(\tilde{\nu})$ is a line profile,
$A_{if}^{(v^{\prime},v^{\prime\prime})}$ is the corresponding Einstein A
coefficient and $N_{\rm r-v}$ is a population of the lower rotation-vibration
state $(v",i)$.

In order to separate the vibrational and rotational populations in our
analysis, we use the following approximation of a product of these
contributions:

$N_{\rm r-v}=N_{v}N_{\rm rot}$ (2)

where $N_{v}$ is a population of the vibrational (or vibronic) state $v$ and
$N_{\rm rot}$ is the population of the rotational state. For the rotational
population a Boltzmann distribution is assumed with the rotational temperature
$T_{\rm rot}$ the same as structural temperature of the atmosphere as defined
by the atmospheric T-P profile, taken for calculating the scale height. This
approximation is similar to the Treanor distribution where the rotational and
vibrational components are separated as discussed in Treanor et al. (1968).
The Treanor distribution takes the product of two Boltmann distributions; one
in terms of the rotational temperature and one in terms of the vibrational
temperature. We make this assumption for the temperature in Eq. (1) such that
$T=T_{\rm rot}$. The effect of the vibrational population is to scale each of
the band intensities by the population of the state $v$. For example, the
vibrational Boltzmann population is obtained using the formula,

$N_{v}=\frac{e^{-E_{v}/kT_{\rm vib}}}{Q_{v}(T)}$ (3)

where $E_{v}$ is the vibrational energy of the state, $Q_{v}(T)$ is the
temperature-dependent vibrational partition function and $k$ is the Boltzmann
constant. We can relate this back to Fig. 2 which shows this for the
circumstance where $T_{\rm vib}=T_{\rm rot}$, i.e. under the LTE condition. In
this work we aim to retrieve the populations $N_{v}$ directly from the
exoplanetary spectrum as an independent measure of the atmospheric temperature
via $T_{\rm vib}$.

At this step, we use ExoCross to generate cross sections for OH with bespoke
populations $N_{\rm r-v}$ for each of the states defined in the MoLLIST OH
line list. The ExoCross software applies these populations for each state with
the corresponding vibrational quantum number. Individual vibrational
populations can be provided in the ExoCross program in order to isolate each
desired band; in each case the vibrational population $N_{v}$ was set to 1 for
the desired state and 0 for the remaining states and a partition function is
inferred. Using this method we produced seven sets of cross sections for each
of the lower states ($v^{\prime\prime}=0,1,2,3,4,5,6$) of the hydroxyl radical
for the $\Delta v=2$ and $\Delta v=3$ bands; i.e. those where opacity falls
into the NIR window $<1.8$ $\mu$m. The corresponding population of each band
was normalised to 1 for temperatures ranging between 2850 K and 3350 K and
pressures between $1\times 10^{-8}$ bar and $100$ bar which covers the T-P
profile of WASP-33b which follows a Guillot (2010) profile as shown in Fig. 1.
The resulting cross sections are plotted in Fig. 3 for a rotational
temperature of $T_{\rm rot}=3050$ K and $N_{v}=1.0$, where individual spectral
components of the hydroxyl radical are shown; all the bands are overlaid in
the top pane to show their relative intensities and then plotted individually
below for clarity.

Figure 1: The Guillot T-P profile adopted for modelling the atmosphere of
WASP-33b.

Bands with different lower states excitations $v^{\prime\prime}\leq 5$ are
indicated by different colours. The stronger features on the right are from
$\Delta v=2$ bands, while the weaker features on the left are from the $\Delta
v=3$ bands. The labels $v0$, $v1$, $v2$, …indicate the lower vibrational
states 0, 1, 2, …, respectively.

For this case, we see that the strengths of the $\Delta v=2$ bands increase
with $v$ quadratically, using an extension of the double harmonic
approximation to the overtone bands (Clark & Yurchenko, 2021), while the
intensity increase of the $\Delta v=3$ band is approximately cubic (Pastorek
et al., 2022). Due to the relative line strengths, the $\Delta v=3$ sequence
will be far less significant than the $\Delta v=2$ band features until the
$v5$ ($v^{\prime\prime}=5$) and $v6$ ($v^{\prime\prime}=6$) systems. The
weaker $\Delta v=3$ bands are also included as present within our chosen
spectral range. As the $\Delta v=2$ moves away from the window, for
$v^{\prime\prime}>4$, only $\Delta v=3$ bands remain visible.

To come a full circle for the sake of illustration, the original LTE spectrum
at a given $T$ can be then reconstructed from these $N_{v}=1$ vibrational
band-isolated cross sections by scaling the individual band cross sections for
given $T=T_{\rm rot}$ to the actual population number $N_{v}$ in Eq. (3). For
example, the cross sections in Fig. 2 were obtained by scaling the cross
sections from Fig. 3 using $N_{v}$ from Eq. (3). Accordingly, there is an
additional drop in their intensities with increasing $v$.

Figure 2: Individual band contributions for cross sections of OH in LTE with
the NIR ($\Delta v=2$ and $\Delta v=3$) shown for the 3050 K, 0.01 bar point
on the temperature pressure grid.

##  3 Forward Modelling

Forward modelling of the emission spectrum of WASP-33b was performed by taking
the requisite input parameters from Nugroho et al. (2021), listed in Table 1.
These planetary and stellar parameters were used along with the generated
cross sections and the atmospheric modelling code TauREx 3 (Al-Refaie et al.,
2021) to produce atmospheric models. The standard, unfiltered cross sections
were used to generate a forward model of OH in WASP-33b’s atmosphere as shown
in Fig. 4; this was performed with the raw OH cross sections, allowing the
retention of the native high resolution models for high resolution cross
correlation analysis and the later downbinning during the observation
simulations for the low resolution space borne analysis.

Figure 3: NIR cross sections of OH ($\Delta v=2$ and $\Delta v=3$) produced from the completely populated states plotted overlaid in the top pane for scale and individually below, all at a temperature of 3050 K and a pressure of 0.01 bar. (Note that as we progress to the final two cross section sets the $\Delta v=4$ and $\Delta v=5$ bands start to become visible) Parameter | WASP-33b Value
---|---
Planet Radius ($R_{J}$) | 1.679
Planet Mass ($M_{J}$) | 3.266
Stellar Temperature (K) | 7400
Stellar Radius ($R_{\odot}$) | 1.509
Stellar Metallicity (dex) | 0.1
Semi-major Axis (AU) | 0.02558
Period (Days) | 1.21987
OH Volume Mixing Ratio (VMR) Profile Parameters |
Surface VMR |  7.94328x10-7
Middle VMR |  6.30957x10-5
Top VMR |  1x10-10
Surface Pressure (Pa) |  1x107
Middle Pressure (Pa) |  3x103
Top Pressure (Pa) |  1x10-3
Table 1: Planetary and stellar parameters used in forward modelling WASP-33b’s
atmosphere.

##  4 James Webb observation simulations and retrievals

###  4.1 Observation simulations

Simulations of observations of WASP-33b’s atmosphere with JWST using the G140H
configuration of the NIRSpec instrument were conducted using ExoWebb (Edwards
& Stotesbury, 2021) noise simulation along with Pandeia data (Pontoppidan et
al., 2016) incorporated into TauREx’s instrument simulation code. The
simulation is shown in Fig. 4 for a 10 secondary eclipse simulation of
WASP-33b, overlaying the TauREx forward model for comparison. An increasing
number of simulated secondary eclipses yielded simulated observations with
successively decreasing noise in the relative flux values. The 10 secondary
eclipse simulation shows the visibility of OH lines clearly in the wavelength
region above 1.4 $\mu$m.

Figure 4: Observation simulation points and errors of the WASP-33b emission
spectra plotted over the forward model for observations over 10 secondary
eclipses.

###  4.2 Retrievals of vibrational populations

Atmospheric retrievals of the populations of individual bands were performed
with TauREx using a new non-LTE plugin (Wright (in prep.)) on the simulated
observation data using cross sections constructed as detailed in Section 2.
Nested sampling was performed to fit for the log of the mixing ratios for the
cross sections included in each case (log$({\rm OH}n)$ with uniform priors
between -12 and -1) and also the planetary radius ($R_{p}$ with uniform priors
between 1.5 and 2 $R_{J}$).The fit for planetary radius is included to ensure
any degeneracies in fitting are visible. From this a derived mean molecular
weight ($\mu$) is calculated and its distribution displayed alongside the
inferred probability distributions; this is included as a summary
characteristic of the retrieved atmosphere.

The first set of retrievals was conducted using the cross sections
corresponding to the vibrational Boltzmann LTE case as our baseline, where
each band is effectively considered as an individual molecule with the
population pre-defined using the Boltzmann LTE rule in Eq. (3) with $T_{\rm
vib}=T_{\rm rot}$ following the T-P profile discussed in Section 2. This can
be directly correlated with a population of an upper vibrational state
$v^{\prime}$ or equivalently of a lower state $v^{\prime\prime}$, since we
only deal with states differing by $\Delta v=2$ or $\Delta v=3$.

Based on this property, in the following, the vibrational bands in question
$v0$, $v1$, $v2$ etc. as labelled in Fig.3 will be referenced as bands 0, 1, 2
etc., i.e. using the corresponding lower state $v^{\prime\prime}$ vibrational
quantum number.

The mixing ratios for each band associated with $v^{\prime\prime}$ were
considered in two ways: (i) following the three point mixing profile used in
the forward modelling step and (ii) also using a constant mixing ratio. Though
given the altitude of the thermal gradient modelled, the lower and upper
points of the mixing profile cannot be well constrained on account of the
limited pressure region probed by this wavelength, as such we opted to proceed
to perform retrievals with a constant mixing ratio profile.

Figure 5: Retrieval nested sampling posteriors for one simulated secondary
eclipse with bands 0 to 2. Figure 6: Retrieval nested sampling posteriors for
one simulated secondary eclipse with bands 0 to 3. Figure 7: Retrieval nested
sampling posteriors for five simulated secondary eclipses with bands 0 to 3.
Figure 8: Retrieval nested sampling posteriors for ten simulated secondary
eclipses with bands 0 to 3.

Figure 5 shows the retrieved posterior probability distributions with the
three strongest emission bands of OH in NIR with $\Delta v=2$ and $\Delta
v=3$, bands 0, 1, and 2 (i.e. corresponding to the emissions from three states
$v^{\prime}=2,3$ and 4 with the highest populations). For this three band
case, treating each band as a separate species, a single secondary eclipse
observation is sufficient to well constrain the individual band mixing ratios
($10^{-4.53}$, $10^{-4.51}$ and $10^{-4.60}$) to within reasonable agreement
with the ground truth $10^{-4.2}$ (all within $2\sigma$).

When we include the hotter band 3, the lines of which are weaker on account of
the reduced relative populations, additional secondary eclipses are required
to well constrain the band mixing ratios. The posteriors for a single
secondary eclipse retrieval are shown in Fig. 6, while the posteriors for
retrievals conducted on 5 and 10 secondary eclipses are shown in Figs. 7 and 8
respectively. In this final retrieval, the band abundances are neatly
constrained, with little degeneracy.

###  4.3 Spectroscopic temperature fitting from vibrational populations

We then attempted to deduce the temperature spectroscopically from the
individual vibrational populations using the simulated data, again treating
each emission band individually. We used the vibrational band cross sections
with the population normalised to 1 ($N_{v}=1$) introduced in section 2. The
retrieved mixing ratio now acts as a proxy for the band population (in
addition to the information it contains about the overall mixing ratio of OH
in the atmospheric model).

When considering the first three bands in isolation, using data simulated for
five secondary eclipse observations, mixing ratios for each of the individual
bands can be well constrained (see Fig. 9). The crucial difference compared to
a typical retrieval of the mixing ratio for a molecule is that we see a drop
off in retrieved mixing ratio with increasing $v^{\prime\prime}$,
corresponding to the reduced populations in higher bands for OH in LTE at this
temperature. The relative weakness of the lines due to these higher states can
be visualised in Fig. 10 where the individual contributions used in the
retrieval are plotted. The limit of our ability to retrieve such mixing ratios
for less populated bands with a simulated ten secondary eclipses of data is
apparent in Fig. 11 and is well understood by their marginal contribution to
the OH spectra as well as by the limitations of the instrument.

Figure 9: Posteriors for five simulated secondary eclipses with bands 0 to 3,
using artificial completely populated band cross sections. Figure 10:
Individual band contributions of OH in NIR retrieved using the completely
populated band cross sections for 0 to 3 with a simulated five secondary
eclipses, with the plot highlighting the wavelength region with the strongest
lines. Figure 11: Posteriors for ten simulated secondary eclipses with bands 0
to 5, using artificial completely populated band cross sections.

In order to illustrate the limitations of detectability imposed by the low
state populations and the associated limited emissions of very hot bands, in
Fig. 12 we examine the case of 50 secondary eclipse observations. One can see
that even for this idealised case, only one more band can be constrained
properly. For the subsequent two bands, only an upper limit can be given and
the 1$\sigma$ confidence interval for the mixing ratio of these bands spans a
wide, poorly constrained range, as shown in Fig. 13.

Figure 12: Posteriors for the retrieval including bands 0 through 6 for a
simulated 50 secondary eclipses. Figure 13: Mix ratios retrieved for bands 0
through 6 for a simulated 50 secondary eclipses.

We can now use this 50 simulated secondary eclipse retrieval to demonstrate
spectroscopic retrieval of the vibrational temperature $T_{\rm vib}$ via
fitting the retrieved populations with a Boltzmann distribution. The band
mixing ratios retrieved individually using the individual complete band sets
were normalised by their sum. We then employ a simple least squares fitting
routine with the vibrational temperature as our free parameter and the
vibrational populations $N_{v}(T)$ are given by Eq. (3), from which we derive
a best fit vibrational temperature of 2815 K. It should be noted that this
simulated example is not an independent measurement of the thermal conditions,
since the literature T-P profile used in simulation is assumed in the
retrieval. Instead it demonstrates the retrievability of vibrational
populations as a technique which could be used in conjunction with fitting the
T-P profile when conducting such an analysis in the empirical case, with an
unknown temperature profile.

The agreement between the Boltzmann distribution and the retrieved mixing
ratios for the bands, as proxies for the band populations, are plotted in Fig.
14. The bars show two population levels for each vibrational number, with the
values given in table 2; these arise since OH is an open-shell molecule and so
fine structure splitting gives rise to two energies per vibrational number.

###  4.4 non-LTE retrieval

In addition, we consider a hypothetical non-LTE scenario where the vibrational
state populations have been driven away from a Boltzmann distribution, similar
to that shown for the excited bending vibrational mode of H2O by Han et al.
(2022). In our non-LTE scenario, we take the molecules to have been
vibrationally excited to an arbitrary peak of $v=3$ with the remaining state
populations following a normal distribution around this peak, as shown in
table 2. This normally distributed configuration is an example of non-
Boltzmann population; the retrieval result is illustrated in Fig. 15. An
example of a physically occurring non-LTE scenario is the non-LTE vibrational
population of OH in the air glow in the upper layers of the terrestrial
atmosphere (Chen et al., 2021; Noll et al., 2020).

$v$ | Normal population | Boltzmann population spin 1 | Boltzmann population spin 2
---|---|---|---
0 | 0.0044 | 1.0000 | 0.9086
1 | 0.0540 | 0.1614 | 0.1467
2 | 0.2420 | 0.0283 | 0.0258
3 | 0.3990 | 0.0054 | 0.0049
4 | 0.2420 | 0.0011 | 0.0010
5 | 0.0540 | - | -
6 | 0.0044 | - | -
Table 2: Populations corresponding to $v$ states for the normal distribution
(a non-LTE) example and the Boltzmann LTE example for both populations arising
from the spin splitting.  Figure 14: Retrieved mixing ratios for OH compared
with the theoretical Boltzmann distribution used in generating cross sections
for OH in LTE. Retrieved from a simulated 50 secondary eclipses. Figure 15:
Retrieved mixing ratios for OH compared with a simulated spectra generated
using OH with a disequilibrium, normal distribution peaking at $v=3$.

##  5 High resolution cross-correlation approach

A high resolution cross-correlation approach comprises part of this analysis.
While we have demonstrated the theoretical capability to separate individual
bands at the modelling stage and to integrate this into a retrieval framework,
we were motivated by the limited resolution of such data from space-based
observatories to investigate the use of high resolution data from ground-based
facilities. In this case, that is the IRD instrument on the Subaru telescope.
To achieve this we perform high resolution cross-correlation fits for the free
parameters: orbital velocity $K_{\mathrm{P}}$, systemic velocity
$v_{\mathrm{sys}}$, and the scale factor of the model $\alpha$. These fits are
performed for a number of forward models, each of which has a fixed set of
vibrational band contributions.

###  5.1 Observations and Data Reduction

We used three data sets taken after the secondary eclipse using IRD mounted at
the Subaru 8.2-m telescope. The first data set was obtained on October 1, 2020
(N1, PID: S20B-008, PI: Nugroho) and used to detect OH emission, see Nugroho
et al. (2021). While the additional two data sets were taken on August 30,
2021 (N2) and December 18, 2021 (N3) as part of program S21B-111 (PI: Nugroho,
see Table 3 for the observation log). The target was continuously observed in
natural guide star mode without the laser frequency comb with an exposure time
of 300 s per frame and overhead time of 25 s. We discarded the first nine
exposures of N3 which had low S/N due to high humidity. These data sets are
part of the work of Nugroho et al. (in prep.) and we refer the reader there
for full details.

Name of | Date | Observing time | Airmass | Phase | Number | Mean S/N | # SysRem
---|---|---|---|---|---|---|---
the night | (UT) | (UT) | change | coverage | of frames | (at 1 $\mu$m) | iterations
N1 | 2020-10-01 | 10:35-13:42 | 1.12-1.05-1.12 | 0.597-0.700 | 33 | 145 | 3
N2 | 2021-08-30 | 10:25-15:15 | 1.64-1.05-1.08 | 0.569-0.734 | 56 | 156 | 9
N3 | 2021-12-18 | 05:24-10:18 | 1.13-1.05-1.44 | 0.573-0.740 | 55 | 110 | 6
Table 3: Summary of IRD observations

In short, the data were reduced using the IRD pipeline (Kuzuhara et al., 2018,
Kuzuhara et al. in prep.). The extracted spectra consist of 72 spectral orders
ranging from $\approx$9260 to $\approx$17645 Å . We exclude the first eight
and the last spectral orders due to relatively low S/N. In addition, we also
exclude the spectral orders with the median wavelength between 14600-15203 Å
from further analysis due to heavy telluric contamination. We fit the
continuum profile of the spectrum with the highest S/N of each data-set using
the continuum task in IRAF111The Image Reduction and Analysis Facility (IRAF)
is distributed by the US National Optical Astronomy Observatories, operated by
the Association of Universities for Research in Astronomy, Inc., under a
cooperative agreement with the National Science Foundation.. Then we divided
the data sets by their corresponding best-fit continuum profile to normalize
them. The sky emission lines were masked based on Oliva et al. (2015);
Rousselot et al. (2000). The spectra were put into a ”common blaze
function/continuum profile” following the procedure of Merritt et al. (2020).
Any 5 $\sigma$ outliers were masked from the data following Gibson et al.
(2020). Any remaining bad pixels and bad regions were visually identified and
masked out from the data. The final spectra of each data set were then stacked
into an array with three axes: number of spectral orders (58), number of
exposures, and wavelength (2048). Finally, we also masked the pixels that have
a value less than 0.1 from the data to avoid any saturated telluric lines and
estimated the uncertainty of each pixel by calculating the outer product of
the standard deviation along the wavelength and time, then normalized by the
standard deviation of the whole array excluding the masked pixels.

In the wavelength range of IRD, telluric and stellar lines dominate the data.
The signal of the planet is several orders of magnitude weaker than the signal
of the host star. To extract the planet signal using the cross-correlation
technique pioneered by Snellen et al. (2010), the telluric and stellar lines
must be removed from the data. Fortunately, the radial velocity of the planet
ranged from tens to hundreds of km s-1, shifting the planet spectra across
many different pixels, while the telluric and stellar lines were (quasi-)
stationary during the observations. To remove them, first, we divided out the
mean spectra from each data set. We then applied the detrending algorithm
SysRem (Tamuz et al., 2005) independently to each spectral order in each data
set to remove the temporal variations in each wavelength bin. This was first
applied to high-resolution datasets by Birkby et al. (2013). Any wavelength
bin in the residual of each iteration of each spectral order that has a
standard deviation of more than 5 times the standard deviation of the whole
array was masked. We stopped the iteration when the detection strength of OH
is at its maximum (see Table 3 for the number of iterations selected for each
data set).

###  5.2 Cross-correlation and Likelihood Mapping

We cross-correlated the telluric-cleaned data sets with our forward models
Doppler-shifted using linear interpolation over a range of radial velocity. We
used the cross-correlation function of Gibson et al. (2020):

$\mathrm{CCF}(v)=\sum_{i}\frac{f_{i}m_{i}(\nu)}{\sigma_{i}^{2}},$ (4)

where CCF is the cross-correlation function for a given Doppler-shifted
velocity ($\nu$), $f_{i}$ is the mean-subtracted telluric-cleaned data at
$i$th wavelength bin, $\sigma_{i}$ is the uncertainty, and $m_{i}$ is the
mean-subtracted model Doppler-shifted to a radial velocity of $\nu$.

The cross-correlation functions of all of the data sets were then co-added at
the planetary rest frame over a range of orbital velocity ($K_{\mathrm{P}}$)
and systemic velocity ($v_{\mathrm{sys}}$) creating a
$K_{\mathrm{P}}$-$v_{\mathrm{sys}}$ CCF map. We then converted it to an S/N
map by dividing it by the standard deviation of the area $\pm$ 50 km s-1 away
from the expected signal. Based on the previous studies, the planetary signal,
if any, would appear around a $K_{\mathrm{P}}$ of 230 km s-1 and
$v_{\mathrm{sys}}$ of -3 km s-1 (e.g., Nugroho et al., 2017, 2020, 2021; Cont
et al., 2021, 2022; Yan et al., 2021; van Sluijs et al., 2022; Yan et al.,
2022).

To avoid any bias in selecting which area to represent the noise, we followed
Brogi et al. (2022) by fitting a Gaussian function to the CCF distribution of
the noise and used the best-fit standard deviation as the estimated noise
floor.

We also calculated a likelihood map ($\mathcal{L}$) for each model using the
$\beta$-optimised likelihood function following Gibson et al. (2020, see also
).

$\ln\mathcal{L}(v)=-\frac{N}{2}\ln\left[\frac{1}{N}\left(\sum\frac{f_{i}^{2}}{\sigma_{i}^{2}}+\alpha^{2}\sum\frac{m_{i}(\nu)^{2}}{\sigma_{i}^{2}}-2\alpha\mathrm{CCF(\nu)}\right)\right],$
(5)

where $\alpha$ is the scale factor of the model and $N$ is the total number of
pixels.

Previous studies (e.g. Brogi & Line, 2019; Gibson et al., 2022) show that
removing the telluric lines using a detrending method (e.g. airmass
detrending, singular value decomposition, SysRem) alters the planetary
spectrum in the data. As cross-correlation or likelihood mapping is basically
a template matching if the model is processed differently than the observed
data any constraint that relies on the line profile (strength, position,
width, and shape) would be unreliable. To correct this effect, for each
spectral order of each data set, we Doppler-shifted the model template for a
given value of $K_{\mathrm{P}}$ and $v_{\mathrm{sys}}$ as a function of
orbital phase and divided out the mean spectrum to normalize it. We then
applied the preprocessing technique described by Gibson et al. (2022) that
uses the same SysRem basis vectors that were used to remove the telluric lines
and subtracted them out from the normalized Doppler-shifted model array. This
preprocessed model array was then directly used to calculate the likelihood
value. We repeated this process for a range of $K_{\mathrm{P}}$ and
$v_{\mathrm{sys}}$. To account for different scaling of the model to the data,
we calculated the likelihood for a range of $\alpha$ values. Lastly, previous
literature (e.g., Gibson et al., 2020; Nugroho et al., 2020; Gibson et al.,
2022) estimated the significance of the detection by dividing the median
values of the conditional distribution of $\alpha$ at the best-fit parameters
by its standard deviation. Here, we expanded it to produce a significance map
by performing it to all of the combinations of $K_{\mathrm{P}}$ and
$v_{\mathrm{sys}}$.

###  5.3 LTE forward model-all bands

Using our OH LTE forward model, we obtained a strong detection of OH emission
with an S/N of 9.5 at the $K_{\mathrm{P}}$ of 231.8 km s-1 and
$v_{\mathrm{sys}}$ of $-$1.1 km s-1 consistent with the previous detection
(Fig. 16a). With the preprocessed model and likelihood mapping, the same
signal is detected at $>$ 12 $\sigma$ (Fig. 16b) demonstrating the advantage
of the likelihood framework and preprocessing the model over a pure cross-
correlation without preprocessing the model in detecting atmospheric signal.

Figure 16: (a). The $K_{\mathrm{P}}$-$v_{\mathrm{sys}}$ map after combining
the cross-correlation functions of all of the data sets using a non-
preprocessed OH LTE model. The upper panel shows the S/N curve at the
$K_{\mathrm{P}}$ of 231.8 km s-1 while the right panel shows the S/N curve at
the $v_{\mathrm{sys}}$ of -1.1 km s-1. The color-bar shows the S/N of the
$K_{\mathrm{P}}$-$v_{\mathrm{sys}}$ map. (b). The marginalized likelihood
distribution of $K_{\mathrm{P}}$, $v_{\mathrm{sys}}$, and $\alpha$ using a
preprocessed OH LTE model. The red dashed lines show the median value of the
corresponding distribution. The upper right panel is the
$K_{\mathrm{P}}$-$v_{\mathrm{sys}}$ map of the detection from the likelihood.
The color-bar shows the significance of the
$K_{\mathrm{P}}$-$v_{\mathrm{sys}}$ map. This map is produced by dividing the
median value of the conditional likelihood distribution of $\alpha$ by its
uncertainty at every combination of $K_{\mathrm{P}}$ and $v_{\mathrm{sys}}$.

###  5.4 Individual bands

Here, we tried to detect the individual OH band assuming a completely
populated state for each band similar to Fig. 3 Rotational state populations
are taken to be Boltzmann distributed, while the vibrational state populations
are unweighted in order to explore any non-LTE (non-Boltzmann) circumstances.
For the high resolution case we extended our isolated cross section set up to
$v=9$ (which corresponds to a maximum highest state of $v=12$ via the $\Delta
v=3$ band) in order to reflect the potential for higher sensitivty than the
lower resolution modelling. We obtained a strong detection both for bands 0
and 1 at $>$ 6 $\sigma$ consistent with the velocity of the planet detected
using the OH LTE model. Unfortunately, there is no clear evidence for hotter
vibrational bands, although there is a hint of a weak signal of band 2 in the
$K_{\mathrm{P}}$-$v_{\mathrm{sys}}$ significance map (see Fig. 17a).

Figure 17: The $K_{\mathrm{P}}$-$v_{\mathrm{sys}}$ significance map of the
combined data sets ($N_{1}+N_{2}+N_{3}$) for the preprocessed individual band
model (a), unpeeling model from band 0 only to the total 0+1+2+3+4+5+6+7+8+9
bands (b), unpeeling model from band 9 only to the sum of bands 0, 1, 2, 3, 4,
5, 6, 7, 8, and 9 (c). The white dashed lines indicate the location of the OH
LTE model signal as a reference. The color bar shows the significance of the
maps.

###  5.5 Band unpeeling

High-resolution spectroscopy, obtained from ground-based instruments, contains
a substantial noise component due to sources such as the telluric lines from
the Earth’s atmosphere; in comparison the component from the actual
Exoplanet’s atmosphere has a low signal-to-noise ratio. With this limitation,
and since the individual bands in isolation provide even smaller signal
contributions, we obtain inconclusive results when performing the cross-
correlation technique with templates containing isolated vibrational bands.

To mitigate the limited contributions from individual bands we attempt a
spectral ‘unpeeling’ approach; cross-correlation templates were constructed
with bands sequentially added from 0 to 9 (i.e. bands 0, 0+1, 0+1+2, 0+1+2+3,
etc), and also backwards from 9 to 0 (i.e. bands 9, 9+8, 9+8+7, …,
9+8+7+6+5+4+3+2+1+0). We propose the unpeeling approach to include
contributions from a greater number of bands in the model templates. This is
to present an approach which may enhance cross-correlation signal in
circumstances where individual band signals are too weak, while still
performing some band isolation across sequential models as is the goal of this
technique. The combination of a number of the vibrational bands then yielded a
recognizable change in the detection significance when comparing cross-
correlation results between different templates (see Fig. 17b & 17c). This
change is especially significant from 0 to 0+1 consistent with the 0 and 1
states being mainly populated at this temperature (see Section 5.4 and Fig.
17a) and the high detection significance of backwards unpeeling model
9+8+7+6+5+4+3+2+1 (i.e. with band 0 excluded). Meanwhile, the 9+8+7+6+5+4+3+2
backward model (with bands 0 and 1 excluded) does not yield any detection
although visually we can see a possible signal of this model around the
expected velocity of the planet similar to the individual band 2 model in
Section 5.4. This is also supported by the increase of the detection
significance when band 2 is added into the 0+1 band model. Interestingly, the
single band 0 of the unpeeling band model yields a much stronger detection
(9$\sigma$) than the completely populated (unity scaled) $v0$ case
(6.6$\sigma$). This difference in the detection significances arises from a
difference in the model templates, driven by differing line intensities. These
line intensities are governed by the vibrational state populations which in
the case of the completely populated model have been artificially set to
unity, while in the filtered case the populations obey a Boltzmann
distribution but with filtering applied afterwards to select the $v0$ model.

The difficulty of detecting vibrational hot bands is in principle in
accordance with the LTE conditions corresponding to the T-P profile
temperature; the vibrational state populations relate to the detection
significances for each individual band model, with only the first two detected
in line with Boltzmann distributed populations corresponding to a vibrational
temperature consistent with the effective temperature. This difficulty of
detection is also related to the limited NIR region used, which does not
properly cover the contributions of the hot bands $v^{\prime\prime}>5$.

##  6 Discussion

Although the wavelength separation is good for the individual hot bands of
molecules such as OH, the ability to distinguish between the bands in all
situations is not unlimited and sufficient wavelength range coverage is
required to avoid degeneracy. This limitation becomes apparent when attempting
retrievals on spectra for state populations peaking at higher vibrational
numbers. For instance, using the current wavelength coverage of 0.97-1.5
$\mu$m in atmospheric retrievals, it would be difficult to detect the
terrestrial-like airglow population, which is known to peak at the higher
vibrational band $(8,6)$ (Chen et al., 2021). When attempting a retrieval of
this non-LTE scenario, we encountered a substantial overestimation of the
corresponding vibrational band, $(7,5)$ or $(8,6)$, abundances characterised
by a convergence to the top boundary of the prior when fitting via a nested
sampling method. This problem exists for these higher bands when observing
over this wavelength range since the stronger band becomes truncated, to
constrain these bands a wavelength coverage reaching higher wavelengths would
be required.

As the ground truth population for a given state becomes small, the associated
transition lines in the spectrum become weaker. When this weakness becomes
sufficiently small compared to the observation uncertainty (simulated, in this
case), the populations for these states are no longer retrievable. This can be
seen in Fig. 14, where it is characterised by a far wider $1\sigma$ interval
for the retrieved $v=4$ population.

It is theoretically possible to add an additional level of complexity and also
attempt retrieval of fine structure spin splitting state populations, for
which theoretical values are represented by the overlapping bars for each
vibrational number in Fig. 14. OH in this electronic configuration is
characterised by spin splitting of 126.3 cm-1. The quantum number assignments
available in this line list are sufficient to do this theoretically, however
since observational data quality is limited, we leave this for future work.

WASP-33 is a pulsating $\delta$-Scuti variable (Changeat et al., 2022) which
makes it a slightly more complex system; this should not be an issue for the
high resolution analysis as there is no OH in the stellar atmosphere to affect
the cross-correlation results.

Possible further work could explore using this technique for a range of other
molecules (such as TiO or other diatomics due to their good vibrational band
separation) and a range of other data sets; including different wavelength
ranges and target planets. One such candidate could be WASP-18b where OH has
also been detected (Brogi et al., 2022).

##  7 Conclusion

From the analysis conducted here, we can conclude that individual vibrational
state populations can be obtained for OH in NIR as fitted parameters in
exoplanet atmosphere retrieval techniques by including isolated individual
bands as cross sections of individual species. Here we have shown the
viability of this approach for OH, though this generalises to diatomic
molecules more broadly on account of their good vibrational band separation.
Not only is this shown for simulated observational spectra for OH with a
typical LTE Boltzmann distribution across vibrational states, but also for a
non equilibrium case. In addition, we demonstrate that non Boltzmann
distributions are also retrievable: for the non equilibrium case, this
approach can also distinguish where higher states are populated, as in the
example of the $v=3$ peak normal population model. The detection of such a
higher peak, or retrieved populations incongruous with the Boltzmann
distribution, in observed data would be consistent with a non-LTE signature.

To facilitate the individual vibrational band detection, an ‘unpeeling’ model
is suggested, where individual band contributions are added to or subtracted
from the total cross sections. Applying an unpeeling approach to the recent
high-resolution detection of individual bands of OH in WASP-33b, we obtained a
significant difference when band 1 was added to a single band 0 model
illustrating the potential of detection of individual spectroscopic bands.

## Acknowledgements

This research is based [in part] on data collected at the Subaru Telescope,
which is operated by the National Astronomical Observatory of Japan. We are
honored and grateful for the opportunity of observing the Universe from
Maunakea, which has the cultural, historical, and natural significance in
Hawaii. Our data reductions benefited from (PyRAF and) PyFITS that are the
products of the Space Telescope Science Institute, which is operated by AURA
for NASA. We are also grateful to the developers of the Numpy, Scipy,
Matplotlib, Jupyter Notebook, and Astropy packages, which were used
extensively in this work (Virtanen et al., 2020; Hunter, 2007; Kluyver et al.,
2016; Astropy Collaboration et al., 2013; Price-Whelan et al., 2018). Funding:
SW was supported by the STFC UCL Centre for Doctoral Training in Data
Intensive Science (grant number ST/P006736/1) This project received funding
from the European Research Council (ERC) under the European Union’s Horizon
2020 research and innovation programme (grant agreements 758892, ExoAI and
883830, ExoMolHD) and the European Union’s Horizon 2020 COMPET programme
(grant agreement No 776403, ExoplANETS A). Furthermore, we acknowledge funding
by the UK Space Agency and Science and Technology Funding Council (STFC)
grants: ST/R000476/1, ST/K502406/1, ST/P000282/1, ST/P002153/1, ST/S002634/1
and ST/T001836/1. The authors acknowledge the use of the UCL Legion High
Performance Computing Facility (Legion@UCL) and associated support services in
the completion of this work, along with the Cambridge Service for Data Driven
Discovery (CSD3), part of which is operated by the University of Cambridge
Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk).
The DiRAC component of CSD3 was funded by BEIS capital funding via STFC
capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant
ST/R00689X/1. DiRAC is part of the National e-Infrastructure. This work was
supported by JSPS KAKENHI grant Nos. JP22K14092, JP19K14783, JP21H00035,
JP20H00170, 21H04998, and SATELLITE Research from the Astrobiology center
(AB022006). JLB acknowledges funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation program under
grant agreement No 805445.

## Data Availability

The underlying molecular linelist data from the ExoMol project (Tennyson et
al., 2020) are available from the Exomol website (www.exomol.com). Derived
cross section data is available from the corresponding author upon request.

## References

  * Al-Refaie et al. (2021) Al-Refaie, A. F., Changeat, Q., Waldmann, I. P., & Tinetti, G. 2021, ApJ, 917, 37, doi: 10.3847/1538-4357/ac0252
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Birkby et al. (2013) Birkby, J. L., de Kok, R. J., Brogi, M., et al. 2013, MNRAS, 436, L35, doi: 10.1093/mnrasl/slt107
  * Borsa et al. (2021) Borsa, F., Fossati, L., Koskinen, T., Young, M. E., & Shulyak, D. 2021, Nature Astronomy, 6, 226, doi: 10.1038/s41550-021-01544-4
  * Brogi & Line (2019) Brogi, M., & Line, M. R. 2019, ApJ, 157, 114, doi: 10.3847/1538-3881/aaffd3
  * Brogi et al. (2022) Brogi, M., Emeka-Okafor, V., Line, M. R., et al. 2022, arXiv e-prints, arXiv:2209.15548.  https://arxiv.org/abs/2209.15548
  * Brooke et al. (2016) Brooke, J. S. A., Bernath, P. F., Western, C. M., et al. 2016, J. Quant. Spectrosc. Radiat. Transf., 138, 142, doi: 10.1016/j.jqsrt.2015.07.021
  * Changeat et al. (2022) Changeat, Q., Edwards, B., Al-Refaie, A. F., et al. 2022, ApJS, 260, 3, doi: 10.3847/1538-4365/ac5cc2
  * Chen et al. (2021) Chen, Q., Hu, X., Guo, H., & Xie, D. 2021, J. Phys. Chem. Lett., 12, 1822, doi: 10.1021/acs.jpclett.1c00159
  * Clancy et al. (2013) Clancy, R. T., Sandor, B. J., Garcia-Munoz, A., et al. 2013, Icarus, 226, 272, doi: 10.1016/j.icarus.2013.05.035
  * Clark & Yurchenko (2021) Clark, V. H. J., & Yurchenko, S. N. 2021, Phys. Chem. Chem. Phys., 23, 11990, doi: 10.1039/D1CP00839K
  * Cont et al. (2021) Cont, D., Yan, F., Reiners, A., et al. 2021, A&A, 651, A33, doi: 10.1051/0004-6361/202140732
  * Cont et al. (2022) —. 2022, A&A, 657, L2, doi: 10.1051/0004-6361/202142776
  * Edwards & Stotesbury (2021) Edwards, B., & Stotesbury, I. 2021, AJ, 161, 266, doi: 10.3847/1538-3881/abdf4d
  * Fisher & Heng (2019) Fisher, C., & Heng, K. 2019, ApJ, 881, 25, doi: 10.3847/1538-4357/ab29e8
  * Fossati et al. (2021) Fossati, L., Young, M. E., Shulyak, D., et al. 2021, A&A, 653, A52, doi: 10.1051/0004-6361/202140813
  * Gibson et al. (2022) Gibson, N. P., Nugroho, S. K., Lothringer, J., Maguire, C., & Sing, D. K. 2022, MNRAS, 512, 4618, doi: 10.1093/mnras/stac091
  * Gibson et al. (2020) Gibson, N. P., Merritt, S., Nugroho, S. K., et al. 2020, MNRAS, 493, 2215, doi: 10.1093/mnras/staa228
  * Guillot (2010) Guillot, T. 2010, Astronomy and Astrophysics, 520, A27, doi: 10.1051/0004-6361/200913396
  * Han et al. (2022) Han, S., Zhao, B., Conte, R., et al. 2022, J. Phys. Chem. A, 126, 6944, doi: 10.1021/acs.jpca.2c05704
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Schmidt, IOS Press, 87 – 90
  * Kotani et al. (2018) Kotani, T., Tamura, M., Nishikawa, J., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII, ed. C. J. Evans, L. Simard, & H. Takami, 1070211, doi: 10.1117/12.2311836
  * Kuzuhara et al. (2018) Kuzuhara, M., Hirano, T., Kotani, T., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII, ed. C. J. Evans, L. Simard, & H. Takami, 1070260, doi: 10.1117/12.2311832
  * Meinel (1950) Meinel, I. A. B. 1950, ApJ, 111, 555, doi: 10.1086/145296
  * Merritt et al. (2020) Merritt, S. R., Gibson, N. P., Nugroho, S. K., et al. 2020, arXiv e-prints, arXiv:2002.02795.  https://arxiv.org/abs/2002.02795
  * Noll et al. (2020) Noll, S., Winkler, H., Goussev, O., & Proxauf, B. 2020, Atmos. Chem. Phys., 20, 5269, doi: 10.5194/acp-20-5269-2020
  * Nugroho et al. (2020) Nugroho, S. K., Gibson, N. P., de Mooij, E. J. W., et al. 2020, ApJ, 898, L31, doi: 10.3847/2041-8213/aba4b6
  * Nugroho et al. (2017) Nugroho, S. K., Kawahara, H., Masuda, K., et al. 2017, AJ, 154, 221, doi: 10.3847/1538-3881/aa9433
  * Nugroho et al. (2021) Nugroho, S. K., Kawahara, H., Gibson, N. P., et al. 2021, ApJL, 910, L9, doi: 10.3847/2041-8213/abec71
  * Nugroho et al. (2021) Nugroho, S. K., Kawahara, H., Gibson, N. P., et al. 2021, ApJ, 910, L9, doi: 10.3847/2041-8213/abec71
  * Öberg et al. (2011) Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16, doi: 10.1088/2041-8205/743/1/L16
  * Oliva et al. (2015) Oliva, E., Origlia, L., Scuderi, S., et al. 2015, A&A, 581, A47, doi: 10.1051/0004-6361/201526291
  * Pastorek et al. (2022) Pastorek, A., Clark, V. H., Yurchenko, S. N., Ferus, M., & Civiš, S. 2022, Spectra Chimica Acta A, 278, 121322, doi: 10.1016/j.saa.2022.121322
  * Pontoppidan et al. (2016) Pontoppidan, K. M., Pickering, T. E., Laidler, V. G., et al. 2016, in Observatory Operations: Strategies, Processes, and Systems VI, ed. A. B. Peck, R. L. Seaman, & C. R. Benn, Vol. 9910, International Society for Optics and Photonics (SPIE), 381 – 395, doi: 10.1117/12.2231768
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Rousselot et al. (2000) Rousselot, P., Lidman, C., Cuby, J. G., Moreels, G., & Monnet, G. 2000, A&A, 354, 1134
  * Snellen et al. (2010) Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049, doi: 10.1038/nature09111
  * Tamura et al. (2012) Tamura, M., Suto, H., Nishikawa, J., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 84461T, doi: 10.1117/12.925885
  * Tamuz et al. (2005) Tamuz, O., Mazeh, T., & Zucker, S. 2005, MNRAS, 356, 1466, doi: 10.1111/j.1365-2966.2004.08585.x
  * Tennyson et al. (2020) Tennyson, J., Yurchenko, S. N., Al-Refaie, A. F., et al. 2020, J. Quant. Spectrosc. Radiat. Transf., 255, 107228, doi: 10.1016/j.jqsrt.2020.107228
  * Tessenyi et al. (2013) Tessenyi, M., Tinetti, G., Savini, G., & Pascale, E. 2013, Icarus, 226, 1654, doi: https://doi.org/10.1016/j.icarus.2013.08.022
  * Treanor et al. (1968) Treanor, C. E., Rich, J. W., & Rehm, R. G. 1968, J. Chem. Phys., 48, 1798, doi: 10.1063/1.1668914
  * van Sluijs et al. (2022) van Sluijs, L., Birkby, J. L., Lothringer, J., et al. 2022, arXiv e-prints, arXiv:2203.13234.  https://arxiv.org/abs/2203.13234
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: https://doi.org/10.1038/s41592-019-0686-2
  * Wright et al. (2022) Wright, S. O. M., Waldmann, I., & Yurchenko, S. N. 2022, MNRAS, 512, 2911, doi: 10.1093/mnras/stac654
  * Yan et al. (2021) Yan, F., Wyttenbach, A., Casasayas-Barris, N., et al. 2021, A&A, 645, A22, doi: 10.1051/0004-6361/202039302
  * Yan et al. (2022) Yan, F., Pallé, E., Reiners, A., et al. 2022, arXiv e-prints, arXiv:2204.10158.  https://arxiv.org/abs/2204.10158
  * Yurchenko et al. (2018) Yurchenko, S. N., Al-Refaie, A. F., & Tennyson, J. 2018, A&A, 614, A131, doi: 10.1051/0004-6361/201732531