# Carbon monoxide emission lines reveal an inverted atmosphere in the ultra
hot Jupiter WASP-33 b and indicate an eastward hot spot

Lennart van Sluijs,1,2 Jayne L. Birkby,1,2,3 Joshua Lothringer,4 Elspeth K. H.
Lee,5 Ian J. M. Crossfield,6Vivien Parmentier,7 Matteo Brogi,8,9,10 Craig
Kulesa,11 Don McCarthy,11 Keith Powell,12David Charbonneau3
1Department of Astrophysics, University of Oxford, Oxford, OX1 3RH, United
Kingdom
2Anton Pannekoek Institute for Astronomy, University of Amsterdam, Amsterdam,
1098 XH, The Netherlands
3Center for Astrophysics | Harvard $\&$ Smithsonian, 60 Garden Street, Cambridge, MA, 02138, USA
4Department of Physics, Utah Valley University, Orem, UT 84058, USA
5Center for Space and Habitability, University of Bern, Gesellschaftsstrasse
6, CH-3012 Bern, Switzerland
6Department of Physics and Astronomy, University of Kansas, Lawrence, KS
66045, USA
7Atmospheric, Oceanic $\&$ Planetary Physics, Department of Physics,
University of Oxford, Oxford OX1 3PU, UK
8Department of Physics, University of Warwick, Coventry CV4 7AL, UK
9INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, I-10025, Pino
Torinese, Italy
10Centre for Exoplanets and Habitability, University of Warwick, Gibbet Hill
Road, Coventry, CV4 7AL, UK
11Steward Observatory, University of Arizona, Tucson, AZ, USA
12MMT Observatory, University of Arizona, Tucson, AZ, USA E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We report the first detection of CO emission lines at high spectral resolution
in the day-side infrared thermal spectrum of an exoplanet. These emission
lines, found in the atmosphere of the ultra hot Jupiter WASP-33 b, provide
unambiguous evidence of its thermal inversion layer. Using spectra from the
MMT Exoplanet Atmosphere Survey (MEASURE, $R\sim 15,000$), covering pre- and
post-eclipse orbital phases ($0.33<\phi<0.73$), we performed a cross-
correlation analysis with 1D PHOENIX model atmospheres to detect CO at S/N=7.9
at $v_{sys}=0.15^{+0.64}_{-0.65}$ km/s and $K_{p}=229.5^{+1.1}_{-1.0}$ km/s.
However, using the framework of Cross-Correlation-to-log-Likelihood mapping,
we further find that the spectral line depths, as probed by the scaling
parameter, change with phase: the line contrast is larger after the eclipse
than before. We then use the general circulation model SPARC/MITgcm post-
processed by the 3D gCMCRT radiative transfer code and interpret this
variation as due to an eastward-shifted hot spot. Before the eclipse, when the
hot spot is facing Earth, the thermal profiles are shallower, leading to a
smaller line depth despite greater overall flux. After the eclipse, the
western part of the day-side is facing Earth, where the thermal profiles are
much steeper, leading to larger line depth despite less overall flux. We thus
demonstrate that even relatively moderate resolution spectra can be used to
understand the 3D nature of close-in exoplanets, if assessed within the log-
likelihood framework, and that resolution can be traded for photon collecting
power when the induced Doppler-shift is sufficiently large. We highlight that
CO in ultra hot Jupiters is a good probe of their thermal structure and
corresponding dynamics, and does not suffer from stellar activity unlike some
atomic species, such as iron, that also appear in the hot host star spectrum.

###### keywords:

planets and satellites: atmospheres – planets and satellites: fundamental
parameters – techniques: spectroscopic

††pubyear: 2022††pagerange: Carbon monoxide emission lines reveal an inverted
atmosphere in the ultra hot Jupiter WASP-33 b and indicate an eastward hot
spot–A

##  1 Introduction

Ultra hot Jupiters (UHJs) provide the unique opportunity to study gaseous
planets in a strongly irradiated environment to learn about their composition
and global circulation patterns (e.g. for a review see Showman et al., 2020;
Fortney et al., 2021). Their high day-side temperatures above 2200 K provide a
unique window to directly detect volatile species in their vapor phase (e.g.
Visscher et al., 2010; Parmentier et al., 2018; Hoeijmakers et al., 2018;
Lothringer et al., 2018; Kitzmann et al., 2018; Hoeijmakers et al., 2019;
Ehrenreich et al., 2020; Merritt et al., 2021; Cont et al., 2021a). Due to
their tidally locked rotation, they have an extreme day-to-night-side
temperature contrast (e.g. Knutson et al., 2007). Consequently, global
circulation models (GCMs) predict strong global jets and hot spots (Showman et
al., 2009; Menou & Rauscher, 2009; Parmentier et al., 2018; Tan & Komacek,
2019). This combined with their relatively high star-to-planet contrast ratio,
makes them ideal targets for atmospheric characterisation with current
observing facilities (e.g. Kreidberg, 2018; Birkby, 2018).

Hydrodynamic simulations and theoretical calculations predict that hot spots
occur in ultra hot Jupiters with an eastward offset from the sub-stellar point
due to global wind circulation patterns (e.g. Showman & Guillot, 2002; Dobbs-
Dixon & Lin, 2008; Menou & Rauscher, 2009; Dobbs-Dixon et al., 2010; Rauscher
& Menou, 2010; Showman & Polvani, 2011; Perez-Becker & Showman, 2013; Debras
et al., 2020). Most observations support this eastward offset hot spot
prediction (Harrington et al., 2006; Cowan et al., 2007; Knutson et al., 2007,
2009; Charbonneau et al., 2008; Swain et al., 2009; Crossfield et al., 2010;
Wong et al., 2016), but for some UHJs westward offsets have been observed as
well (Armstrong et al., 2016; Dang et al., 2018; Jackson et al., 2019; Bell et
al., 2019; von Essen et al., 2020; Herman et al., 2022). Several mechanisms
including cloud asymmetries, asynchronous rotation and magnetohydrodynamical
effects have been proposed to explain these westward offset hot spots (e.g.
Hindle et al., 2021).

Thermal inversions of the pressure-temperature profiles (P-T profiles) of UHJs
were at first predicted due to the the strong optical/UV absorbing molecules
TiO and VO (e.g. Hubeny et al., 2003). Later work found absorption by other
atomic and molecular species, in addition to TiO and VO, can also result in
thermal inversions (Lothringer et al., 2018; Arcangeli et al., 2018; Gandhi &
Madhusudhan, 2019). This is due to a combination of short-wavelength stellar
irradiation around early type stars and absorption of these wavelengths by
continuous opacity sources, metal atoms, metal hydrides and SiO. Analysis of
Spitzer and Hubble Space Telescope (HST) secondary eclipse observations
support this scenario as they detect emission signatures suggestive of an
inverted stratosphere for a handful of UHJs (e.g. Deming et al., 2012; Haynes
et al., 2015; Evans et al., 2016, 2017; Sheppard et al., 2017; Arcangeli et
al., 2018; Kreidberg et al., 2018; Mansfield et al., 2018; Nugroho et al.,
2020a; Yan et al., 2020; Garhart et al., 2020; Baxter et al., 2021).

One well-studied UHJ for which a thermal inversion was predicted by
atmospheric modeling is WASP-33 b. This planet orbits every 1.22 days around
an A5 star (Collier Cameron et al., 2010). The first observational evidence of
a thermal inversion of WASP-33 b’s atmospheric profile was excess infrared
emission observed by Spitzer from secondary eclipse observations and (Deming
et al., 2012; Zhang et al., 2018) and NIR low-resolution HST/Wide Field Camera
3 spectra, possibly due to TiO (Haynes et al., 2015). Furthermore, (Zhang et
al., 2018) find a phase angle offset by $-12.8\pm 5.8\degree$ in the Spitzer
3.6$\micron$-band, indicative of an eastward hospot. On the contrary, (von
Essen et al., 2020) find a $28.7\pm 7.1\degree$ phase angle offset in the
optical with TESS, suggesting a westward offset hot spot instead, although
they mention host star variability may introduce a spurious westward offset.
Recently, Herman et al. (2022) report day-to-night brightness contrast
variations and a $22\pm 12\degree$ westward phase offset from their detection
of Fe-I emission lines.

For these low resolution observations, stellar pulsations of the $\delta$
Scuti pulsating host star induce variability of the stellar continuum (Herrero
et al., 2011), limiting observational constraints as they induce a quasi-
sinusoidal trend in transit light curve observations (e.g. Deming et al.,
2012; de Mooij et al., 2013; Garhart et al., 2020; von Essen et al., 2019,
2020). Furthermore, stellar pulsations also affect high resolution spectral
observations targeting any molecule present in both the stellar and
exoplanetary atmosphere due to changes in the spectral line shape (Herrero et
al., 2011). For WASP-33 b this has been observed for the detection of neutral
iron lines by Nugroho et al. (2020b); Herman et al. (2022).

Pursuits to follow-up the tentative low-resolution detection of TiO in the
atmosphere of WASP-33 b using high resolution spectra have lead to
inconclusive results so far. Nugroho et al. (2017) detected the the optical
TiO molecular signature using the High Dispersion Spectrograph on the Subaru
Telescope. This detection was later challenged by a re-analysis of the same
data by Serindag et al. (2021) who unexpectedly found a slightly weaker signal
using the updated TiO ExoMol TOTO line list. Furthermore, Herman et al. (2020)
analysed emission and transmission spectra from the Echelle
SpectroPolarimetric Device for the Observation of Stars (ESPaDOnS) on the
Canada-France-Hawaii Telescope and the High Resolution Echelle Spectrometer
(HIRES) on the Keck telescope, but are unable corroborate their observations
with a thermal inversion due to TiO. There is thus a lack of a current
consensus regarding the observational evidence for TiO in the atmosphere of
WASP-33 b.

Despite the elusive results for TiO, several new detections suggest WASP-33 b
must have a thermally inverted atmosphere. Recently, emission signatures of
FeI (Nugroho et al., 2020b; Cont et al., 2021a), neutral Si (Cont et al.,
2021b) and OH (Nugroho et al., 2021) have been detected. Nugroho et al. (2021)
also marginally detect $\textrm{H}_{2}\textrm{O}$, consistent with the
theoretical predictions that $\rm{H}_{2}O$ molecules dissociate into OH and
$\textrm{H}^{-}$ at the temperature regime of UHJs. Other atmospheric
detections for WASP-33 b include Ca II (Yan et al., 2019), the hydrogen
H$\alpha$, H$\beta$ and H$\gamma$ Balmer lines (Cauley et al., 2021; Yan et
al., 2021; Borsa et al., 2021), evidence for AlO (von Essen et al., 2019), and
lastly, Kesseli et al. (2020) place upper limits on the Volume Mixing Ratio
(VMR) for FeH based on their null detection. These detections are thus in line
with predictions by Lothringer et al. (2018); Gandhi & Madhusudhan (2019) of
thermally inverted UHJ atmospheres, regardless of the presence of TiO or VO.

In this paper we report CO emission lines in the atmosphere of WASP-33 b using
observations with the MMT in Arizona, USA, equipped with the ARizona Infrared
imager and Echelle Spectrograph (ARIES). It is the first detection of CO
emission lines detected with high resolution cross correlation spectroscopy
(HRCCS) utilizing the large Doppler-shift induced by the planet’s orbital
motion, and despite the lower resolution of the spectra presented here in
comparison to other HRCCS observations, it also reveals information about the
planet’s atmospheric dynamics. The observations and the subsequent data
reduction are discussed in Section 2. Our models and the methods used to
characterise the exoplanet’s atmosphere are described in Section 3. The
results are presented in Section 4 and discussed in Section 5. Our main
conclusions and recommendations for future work are summarised in Section 6.

##  2 Observations and Data Reduction

###  2.1 Observations

We observed the ultra hot Jupiter WASP-33 b with the ARizona Infrared imager
and Echelle Spectrograph (ARIES; McCarthy et al. 1998; Sarlot et al. 1999) in
combination with the f/15 adaptive secondary mirror and adaptive optics system
at the 6.5-m MMT Observatory on Mt Hopkins in Arizona, USA. The adaptive
secondary mirror provides a low thermal background, while augmenting the total
throughput of the instrument. ARIES can provide spectral observations using
both long slit and echelle spectrograph modes. It can observe in the 1-5
$\micron$ range and at spectral resolutions of 2,000-30,000. The observations
of WASP-33 b presented in this work are the first of a larger survey called
MEASURE (MMT Exoplanet Atmosphere Survey) which contains observations of a
diverse set of eleven exoplanets with a wide range of temperatures, masses and
radii.

For these observations of WASP-33 b, we obtained echelle spectra in the
1.37-2.56 $\micron$ wavelength range during three half nights in October 2016
using the ARIES/MMT combination. In total 211 spectra were obtained with an
exposure time of $300$ s per frame at a theoretical instrumental resolution of
$R=30,000$ using the $1^{\prime\prime}\times 0.2^{\prime\prime}$ slit and the
f/5.6 camera mode. The simultaneously operated ARIES imager was used for in-
slit guiding, with occasional manual offsets made to correct any sustained
drift during an exposure not corrected by the AO guiding system.

An overview of our observed orbital phase coverage111Orbital phases have been
calculated using the ephemerides in Table 4, which contains all relevant
WASP-33 system parameters used throughout this work. is shown in Fig. 1 and
our observational parameters and data quality are summarised in Table 1. The
layout of the 26 echelle orders of ARIES is shown in the top-left panel of
Fig. 2. The three observing nights cover both pre- and post-eclipse orbital
phases, but the post-eclipse data covers significantly more of the orbit,
almost to quadrature, where the planet shows half its day-side and half its
night-side to the Earth.

At the beginning and end of each half night, we observed a set of calibration
images. This included a set of dark frames obtained using a blank in the
filter wheel and the ARIES entrance covered to prevent any light entering the
instrument. We observed flat-fields with the grating in, using an incandescent
light bulb arranged such that its light reflected off an aluminumized board
into the dichroic of ARIES. We also observed a thorium-argon lamp in a similar
manner at the start of the nights to focus the spectrograph, but used the
observed telluric absorption lines for simultaneous wavelength calibration.

Date | Total | No. | Phase | S/N∗
---|---|---|---|---
(UTC) | exposure | spectra | Coverage |
| time (min) |  |  |
2016 Oct 15 | 220 | 44 | 0.331-0.469 | 91
2016 Oct 19 | 465 | 93 | 0.461-0.731 | 26
2016 Oct 20 | 370 | 74 | 0.337-0.554 | 71

∗This quantity refers to the median S/N per wavelength bin
as measured from the extracted and normalised spectral time series
of the 23rd spectral order around the $2.3\micron$ CO emission region.

Table 1: Observational parameters of the WASP-33 data set for each night.
Figure 1: Orbital phase ($\phi$) coverage for each observing night of the
WASP-33 data set. As can be seen, all nights contain some spectra taken during
full eclipse of WASP-33 b. The blue dotted line indicates up to which phases
we have symmetric phase coverage for both the pre- and post eclipse data.

###  2.2 Data reduction: pre-processing

Since this is the first time the MMT/ARIES combination has been used for
exoplanet atmosphere characterisation, we developed a purpose-built end-to-end
data reduction pipeline that we describe here first. This new pre-processing
pipeline will be made publicly available on the author’s
GitHub222https://github.com/lennartvansluijs. As ARIES is currently being
upgraded into MAPS (Morzinski et al., 2020), the instrument is no longer
available in this exact form, but our data reduction pipeline uses a
generalized approach which may be useful to reduction of the instrument when
it returns. The procedure is described below, and is applied in the same way
to all three nights of data:

  1. 1.

Cross-talk correction: ARIES experiences cross-talk between the four quadrants
of its detector which needs to be corrected. A bright pixel in one quadrant
produces distinct, cross-shaped, negative shadows at the same position in all
other three quadrants on the detector. This effect is apparent in all dark,
flat, science frames.

We use the C-based corquad routine provided by the ARIES team333See
http://66.194.178.32/~rfinn/pisces.html to correct for the cross-talk. A
Python-based version is available at: https://github.com/jordan-stone/ARIES..
Three input parameters define the convolution kernel size for each detector
quadrant. The parameter space is explored using a grid-based search to find
the best cross-talk convolution kernels, which we define as those that
minimize the standard deviation in a $10\times 10$ box around a visually
identified prominent shadow feature. We apply the corquad routine using the
best convolution kernels parameters to all dark, flat and science frames,
before any further data reduction.

  2. 2.

Dark correction: we need to remove hot pixels from individual darks and
combine them into a master dark. We identify hot pixels as $\geq
3\sigma$-outliers. Like the shadow cross-talk features, they are cross-shaped
and we replace them by the median of their neighbours outside of a cross-
shaped footing in a $11\times 11$ box around each hot pixel. We median combine
all frames of identical exposure time into a master dark. We subtract the
appropriate master dark from each science and flat frame.

  3. 3.

Flat correction: all flats contain fringes arising from Fabry-Pérot
interference between the optical elements of the instrument. Fringes in the
flat field introduce modulation in the continuum that would be erroneously
divided into the science flat during flat field correction, reducing the
sensitivity of our later procedures. To remove the fringes from the flats, we
developed a modified version of the flat fringe correction by Stone et al.
(2014). The procedure is shown in Fig. 2. We use the getthem routine from the
publicly available Python package CERES444https://github.com/rabrahm/ceres
(Brahm et al., 2017) to fit the echelle traces as 4th order polynomials using
a 10 pixel wide aperture. We use the solution of all 26 spectral orders to
calculate a 2D polynomial transform to dewarp the flat frame. We fit each row
by a 7th order polynomial to create a dewarped illumination model, essentially
representing the dewarped blaze function. This model is subtracted from the
dewarped flat frame to reveal the dewarped flat fringes. Most fringes have
frequencies $<0.025\ \rm{pixel}^{-1}$. We apply a 6th order Butterworth high-
pass filter with a cutoff-frequency $f_{\rm{c}}=0.025\ \rm{pixel}^{-1}$ to
create a dewarped fringe model. We then use the inverse 2D polynomial
transform to warp the fringes model. We remove the fringes from the observed
flat frame by subtracting the warped fringes model. It is key this subtraction
is done in the observed coordinate frame such that the observed flat was not
interpolated by a transform operation, thus preservering any intra-pixel
systematics.

After defringing all flats, we perform a hot pixel correction similar to the
one done for the dark frames and create a master flat. The ARIES detector has
a central column of dead pixels. We correct this column by replacing it by the
median of the two neighbouring columns in the master flat and science frames.
To create a master flat hot pixel map, we first normalize the master flat by
dividing the illumination model and identifying pixels with values outside of
the [0.5, 1.5] range, as hot pixels. We combine the master flat hot pixel map
with the master dark hot pixel map and replace all master flat hot pixels by
the fit from the illumination model. For each science frame we correct each
flagged hot pixel by their closest neighbours in the dispersion direction of
the spectrograph. Finally, the master flat is divided out of all of the
science frames. The science frames still contain fringes at this stage, but
these will be removed by applying a high-pass-filter to the extracted spectral
time series in the post-processing stage of the pipeline.

  4. 4.

Spectral order extraction: after the flat correction we can extract the
spectral orders from the detector images. We cannot directly fit the detector
echelle traces for the science frames as these traces are often disrupted due
to telluric absorption lines. Instead we fit the echelle traces using the
CERES retrace routine. This treats the initial set of master flat traces as an
aperture mask and cross-correlates it with each science frame to determine the
drift in the cross-dispersion direction. We use the CERES Marsh-algorithm
(Marsh, 1989) routine with a 5 pixel aperture radius and clipping $>$5$\sigma$
outliers. The final output of the pre-processing part of the data reduction
pipeline is a spectral time series for each of the 26 spectral orders for each
of the 3 observing nights of the system (see example for one order in the top
panel of Fig. 3).

Figure 2: Flat fringe correction illustrated in six steps: (1) the observed
flat frame with the fitted echelle traces overlaid in red (numbered from blue-
to-red). (2) the flat frame after applying a 2D polynomial transform to
straighten the echelle traces. (3) the 4th order polynomial row-by-row
illumination model. (4) fringes in the dewarped flat frame. (5) row-by-row
low-pass filter fringe model. (6) fringe model after applying inverse 2D
polynomial transform.

###  2.3 Data reduction: post-processing

From the spectral time series we would like to extract the planet spectrum,
but it is faint and buried in the noise. Thus we first need to remove the
telluric and stellar lines. An overview of all steps of the post-processing
procedure for a single spectral order covering the most line dense part of the
CO spectrum is shown in Fig. 3. We perform the following steps to post-process
the extracted spectral time series for each spectral order in a homogeneous
manner:

  1. 1.

Bad pixel/column correction (see Fig. 3, panel 2): in each spectral time
series, some bad pixel columns or detector artifacts may still be present in
the spectral time series. To reveal these anomalies, we normalise each spectra
by the median of each row, followed by subtracting the mean of each columns.

We apply a Laplacian of the Gaussian algorithm, also known as a Blob Detection
Algorithm, to detect outliers (e.g. see Kong et al., 2013), where we identify
features of the size of a single pixel. These outliers are corrected by
interpolating in the dispersion direction over the closest nearby good pixels.
We ignore any convolution boundary effects here as they only affect the edge
pixels and these will be clipped during the post-processing stage of the
pipeline.

Bad columns in the spectral time series are identified by taking the median of
each spectral channel and applying a median-filter using a three pixel width
and subtract this from the original. Any 5$\sigma$ outliers in the residual
time series are considered bad columns and are corrected by interpolation from
their nearest neighbouring good columns. We iteratively perform the bad pixel
and column correction procedure five times to ensure all outliers are
corrected.

  2. 2.

Alignment of spectra (see Fig. 3, panel 2): instrumental flexure and changes
in air mass cause sub-pixel drift of the spectra across the detector. We need
to align the spectra to correct for this effect. We use the first spectrum of
every night as a reference spectrum. This reference spectrum is cross-
correlated with each frame. We fit a Gaussian to this the cross-correlation
function (CCF) to determine the drift on a sub-pixel level for each frame. We
then shift each spectra using linear interpolation to align them with the
reference spectra.

  3. 3.

Telluric wavelength calibration: after alignment, we perform wavelength
calibration using the observed telluric lines. We generate a telluric model
using ATRAN555https://atran.arc.nasa.gov/cgi-bin/atran/atran.cgi (Lord, 1992).
For each aligned spectral order we cross-match by-eye telluric lines in the
model to the corresponding telluric lines in the median of our observed
spectra. We robustly fit a polynomial to all these points for each order to
obtain a wavelength solution. Robust here refers to the fact that we
iteratively clip $\geq 5\sigma$ outliers and refitting a polynomial until no
more outliers are found. By default we fit 3rd order polynomials, but for some
orders where there are few absorption lines towards the edges of the detector,
we fitted 2nd order polynomials instead, to avoid massive deviation of the
wavelength solution towards the edges. All these fits were manually inspected
to ensure a proper wavelength solution was obtained. For some orders the
wavelength solutions were poorly constrained due to strong residual science
fringes and/or weak telluric absorption features. For these reasons we
excluded orders 7-9 and 11 from the rest of the analysis.

  4. 4.

Throughput correction (see Fig. 3, panel 3): throughput variations can be
caused by variations in airmass and misalignment of the target on the slit. We
correct for these variations by dividing by the mean of the brightest 50
pixels per spectrum following the procedure described in Brogi & Line (2019).

  5. 5.

Telluric removal (see Fig. 3, panel 4-5): quasi-stationary trends mostly due
to the telluric lines are still dominating the spectra. To remove these,
Principal Component Analysis (PCA) is used to decompose the spectral time
series. In our implementation, we perform a Singular Value Decomposition (SVD)
(similar to de Kok et al., 2013; Line et al., 2021). Visual inspection of
other ARIES observations from other systems indicates that using seven
Principal Components works well to balance telluric removal and retrieve an
injected exoplanet signal at a high S/N (as described in Section 3.2.1), so we
adopt this also for WASP-33 b. However, we emphasize that we keep the number
PCA iterations fixed for all spectral orders and all observing nights to avoid
optimisation of order/night specific systematic effects and noise, which can
produce artificially large S/N (e.g. Cabot et al. 2020; Spring et al. 2022).

Remaining low-frequency residuals persist due to several effects: (1)
throughput variations due to small offsets of the target position on the
spectrographs slit entrance; (2) residual airmass variations; (3) echelle
traces drift on the detector due to the changing gravity vector during the
observations, causing instrumental flexure; and (4) uncorrected science frame
fringes. We apply a 6th order Butterworth high-pass filter with a cutoff
frequency of $0.02\ \rm{pixel}^{-1}$ to remove most of these trends.

  6. 6.

Masking (see Fig. 3, panel 6): remaining bad columns in the residual time
series are masked. This is done by identifying $>3\sigma$ outliers in the
residual matrix and flagging columns with more than five outliers. We combine
flagged columns if there are more than two columns within a five pixel-wide
sliding window. Additionally, the columns in a 50 pixel window from the edge
of the detector are masked by default.

The final result of the post-processing pipeline is shown in panel 6 of Fig.
3.

Figure 3: Example of post-processing procedure of a spectral time series for
order 23 for the first observing night. First panel: The extracted spectral
time series from the processed detector image is shown in the top panel.
Second panel: after bad pixel/column correction and alignment with respect to
the first spectrum. Third panel: normalised by the median of the 50 brightest
pixels per spectrum. Fourth panel: after removing the first 7 components using
Principal Component Analysis (PCA). Fifth panel: After applying a high-pass-
filter with cutoff-frequency of 0.02 $\mathrm{pixel}^{-1}$. Sixth panel: After
masking the columns with the worst residual structure. Seventh panel: same as
the previous panel, but with a $\times$10 PHOENIX self-consistent model
injected into the spectral time series between panels 2 and 3.

###  2.4 Observed instrument performance

As this is the first time ARIES/MMT is used to characterize an exoplanet
atmosphere, we compare the observed performance to the theoretical
instrumental performance. We measured the instrumental throughput, resolving
power and Precipitable Water Vapour (PWV) directly from the normalised
spectral time series. The instrumental throughput is measured directly from
the total S/N for each spectra. To measure the observed resolving power and
PWV we largely follow procedure described in Chiavassa & Brogi (2019), where
we cross-correlate a telluric model with the normalised spectral time series
and evaluate the log-likelihood defined by Zucker (2003). This method is
sensitive to continuum variations and first we correct for them the following
way: (1) a median filter removes any possible outliers (2) each spectrum is
binned and the maximum of each bin is calculated (3) we fit a low 3rd-order
polynomial to these local maxima and (4) we divide out this polynomial fit
from each spectra. The Telfit666https://github.com/kgullikson88/Telluric-
Fitter code is used to compute telluric spectra (Gullikson et al., 2014). We
fix the airmass to the logged airmass for each frame which we obtain from the
raw FITS-header information, instead of keeping it a free parameter as done by
Chiavassa & Brogi (2019). The Markov Chain Monte Carlo (MCMC)
PyMultiNest777https://johannesbuchner.github.io/PyMultiNest/index.html#citing-
pymultinest (Buchner et al., 2014a) package is used to constrain $PWV$ and
$R$.

The measured relative throughput, instrumental resolving power and PWV are
shown in Fig. 4. A key results is that the measured spectral resolving power
is $R=14,000-16,000$, notably lower than the nominal specification for ARIES
of $R=30,000$. We hypothesize this large offset may be due to sub-optimal
focus of the spectrograph resulting in consistently lower observed spectral
resolution. The other fluctuations of spectral resolving power in the course
of each night could be due to a combination of target misalignment on the slit
entrance, instrumental instability and occasional opening of the AO loop. This
analysis also highlights the importance of measuring the spectral resolving
power from the observations when using HRCCS, as our cross-correlation
templates will need to be convolved to the appropriate spectral resolution
later on.

Figure 4: Relative throughput, spectral resolving power and Precipitable Water
Vapour (PWV) for each observing night, as measured from order 23 which is
centered around the CO 2.3 $\micron$ spectral feature. The dashed lines
indicate the median resolving power used to convolve the spectral templates
for each night. The overall consistently reduced resolution of ARIES is
hypothesised as sub-optimal focusing of the spectrograph.

##  3 Methods

As we use the HRCCS technique to characterise the atmosphere of WASP-33 b, we
require forward models for the cross-correlation templates i.e. synthetic
spectra of WASP-33 b, which we discuss in Section 3.1, followed by a
description of the HRCCS methods used in Section 3.2.

###  3.1 Models for WASP-33 b

We use several atmospheric modeling codes in this work to address the
structure, chemistry, and dynamics of the atmosphere of WASP-33 b. To identify
absorption and emission features initially, we used PHOENIX atmospheric
models. PHOENIX is a general purpose atmosphere code, well-tested on objects
from cool planets to hot stars (e.g. Hauschildt et al., 1997, 1999; Barman et
al., 2001; Barman et al., 2011; Allard et al., 2011; Lothringer et al., 2018;
Lothringer & Barman, 2019). We used both self-consistent atmosphere models in
Section 3.1.1 in thermodynamic equilibrium, and a grid of models where the
temperature structure and/or the molecular abundances have been manually
varied in Section 3.1.2. We then further use the Monte Carlo Radiative
Transfer code gCMCRT (Lee et al., 2019, 2021) to model 3D effects and orbital
phase-dependent variations of the observed spectra in Section 3.1.3.

####  3.1.1 Self-consistent PHOENIX models

The self-consistent models used here are similar to the extremely irradiated
hot Jupiter models presented in Lothringer et al. (2018). The self-consistent
atmosphere models solve the atmosphere structure and composition iteratively,
calculating the spectrum at each iteration and comparing it to radiative
equilibrium. For our models the radiative-convective boundary is usually below
the highest pressure we model and no convective adjustments were made in order
to increase speed of convergence. Temperatures are then modified to approach
radiative equilibrium after which the chemistry is calculated based on
chemical equilibrium and the spectrum is re-calculated. This process is
repeated until the temperature corrections are small, indicating that the
atmosphere is in radiative and chemical equilibrium. Each spectrum is
calculated from 10 to $10^{7}$ nm using planetary parameters by Haynes et al.
(2015) (as summarised in Table 1 of Lothringer et al., 2018) on a grid with
varying atmospheric metallicity and a heat re-radiation efficiency $f$ (where
$f$ is the parameter as defined in Madhusudhan & Seager, 2009). For the
metallicity, we explore a grid of models at 0.1$\times$, 1$\times$ and
10$\times$ Solar metallicity. For the heat re-radiation factors $f$ we explore
the values: $f=1/4$ for day-side and night-side heat redistribution, $f=1/2$
for day-side-only heat redistribution and $f=2/3$ for instantaneous heat re-
radiation.

These self-consistent models include opacities for TiO and VO, but as
explained in Lothringer et al. (2018), TiO and VO are not the sole cause of
temperature inversions in extremely irradiated hot Jupiters. Temperature
increase at these pressures is likely due to the absorption of short-
wavelength radiation by absorbers like atomic Fe and the lack of coolants such
as $\rm{H}_{2}\rm{O}$. Such hot temperatures result in the thermal
dissociation of most molecules in the atmosphere including H2 and H2O (see Fig
15). The abundance of CO, which exhibits the strongest molecular bond, is also
reduced, but to a lesser extent compared to other potential molecular opacity
sources. The resulting spectrum, shown in Fig 5, predicts CO emission lines
with signatures of H2O and other molecules muted due to thermal dissociation
combined with the isothermal lower atmosphere.

Figure 5: WASP-33 b 1D PHOENIX self-consistent synthetic emission spectra.
Left panels: PHOENIX self-consistent spectra as a function of the metallicity
and re-radiation factor $f$: $f=1/4$ for full heat-redistribution, $f=1/2$ for
day-side-only heat-redistribution and $f=2/3$ for instantaneous re-radiation.
Different metallicity models have been offset by 0.75 along the y-axis with
respect to the 0.1$\times$Solar model. The ARIES spectral order wavelength
ranges included are indicated by the gray bars. Right panel: corresponding
PHOENIX self-consistent P-T profiles.

####  3.1.2 PHOENIX models with Modified Structure and/or Abundances

In addition to the self-consistent modeling, we ran a grid of models with
varying temperature structures and atmospheric metallicity to further explore
the parameter space. Again, we explore metallicities of 0.1$\times$,
1$\times$, 10$\times$ and 100$\times$ Solar. For the P-T profiles, we use the
parameterisation from Madhusudhan & Seager (2009), but modified to have an
isothermal upper atmosphere. These simple structures used five parameters to
describe the temperature throughout the atmosphere: $P_{1}$, $P_{2}$, $T_{1}$,
$T_{2}$ and $\alpha_{2}$. $P_{1}$ describes the pressure at the tropopause,
above which the temperature was set to $T_{1}$. The temperature then varies
until it reaches temperature $T_{2}$ at pressure $P_{2}$, as set by the
gradient of the inversion $\alpha_{2}$. We explored a wide range of P-T
structures with upper atmosphere temperatures from 2000 to about 10,100 K,
shown in Fig. 16. In total we ran 576 forward models.

We also generated a small subset of these models to explore how removing a
single opacity source affects the significance of the planet detection i.e we
re-run the best matching model without H2O, without CO, and without OH. We do
this in order to identify potential absorbers or emitters in the atmosphere,
while keeping the P-T structure fixed. Line lists for the investigated opacity
sources were taken from: CO (Goorvitch, 1994), OH (Barber et al., 2006) and
H2O (et al., 2009).

####  3.1.3 GCM and post-processing modelling

For the 3D GCM models, we use output from a SPARC/MITgcm model (Showman et
al., 2009) of WASP-33 b with similar set-up to the simulations performed in
the Parmentier et al. (2018) study for UHJs. The GCM output is then post-
processed using the gCMCRT 3D RT code (Lee et al., 2021) to provide high-
resolution emission spectra that account for variation in the spectra
throughout the orbital phases of the observation. We use line-list data from
various sources, namely: OH (Hargreaves et al., 2019), H2O (Polyansky et al.,
2018), CH4 (Hargreaves et al., 2020), CO (Li et al., 2015), CO2 (Yurchenko et
al., 2020), NH3 (Coles et al., 2019), HCN (Barber et al., 2014). Cross-
sections were calculated using the HELIOS-K opacity code (Grimm et al., 2021)
and interpolated between 1.27-2.66 $\mu$m at a resolution of R = 100,000.

The resulting synthetic high-resolution emission spectra include the Doppler
shifting of spectral lines due to winds and planetary rotation towards the
line of sight for each phase. We include all molecular species listed above as
well as CIA H${2}$-H-He pair continuum opacities from Karman et al. (2019) and
H- from John (1988).

The resulting GCM phase-dependent model and P-T profiles and the gCMCRT output
synthetic spectra are shown in the bottom-left and right panels of Fig. 6. The
GCM predicts an eastward hot spot at $+\degree$ east from the substellar
point, which results in an asymmetry between P-T profiles with an equal
longitudinal distance west or east as seen from the sub-stellar point. As the
planet rotates in the same direction as it orbits its host star, we see more
of the hot spot region pre-eclipse compared to post-eclipse. This results in
phase-dependent variations in both the continuum and relative line strengths
(see bottom-left panel of Fig. 6). As it is computationally expensive to run a
full GCM at all observed orbital phases, we run the full GCM at 9 phases
uniformly sampled over the phase coverage during our three observing nights.
For spectra observed at intermediate phases, we linearly interpolate these
spectra in phase.

To enable comparison and cross-check between the 3D GCM and the 1D PHOENIX
model atmospheres, we also created a set of GCM day-side-only models i.e.
fixed at orbital phase $\phi=0.5$ (see top-left panel of Fig. 6). These
templates were created by simulating the emission assuming continuum opacity
without the effects of Doppler shifting. We also again explore the impact of
any single absorber or emitter on the planet signal by creating day-side-only
templates with only a single opacity source included i.e. only CO, only H2O
and only OH. CO is most prolific in spectral lines in the reddest parts of the
ARIES spectra, while OH dominates its shorter wavelengths, and water persists
throughout.

Figure 6: WASP-33 b GCM (Global Circulation Model) synthetic emission spectra.
Top-left panel: GCM day-side-only (i.e. fixed at phase $\phi=0.5$) spectra
with all opacity sources included or only specific molecules included. The
different models have been offset by 0.5 along the y-axis with respect to the
spectrum at the bottom. The ARIES spectral order window wavelength ranges
included are indicated by the gray bars. Bottom-left panel: GCM spectra at
different orbital phases. Planck black body curves are plotted at a range of
temperatures. The GCM predicts hotter models pre-eclipse compared to post-
eclipse due to an hotspot approximately $+8\degree$ eastward seen from the
substellar point. Right panel: GCM P-T profiles as a function of longitude,
where $0^{\circ}$ is defined as the substellar point and increases towards the
eastward direction as seen from the substellar point.

###  3.2 High-resolution Cross Correlation Spectroscopy (HRCCS)

The HRCCS technique has been developed in recent years to enable not simply
detection of molecular species, i.e. the S/N-method, but into a framework that
enables statistically rigorous assessment of atmospheric models that may match
the data e.g. Cross-Correlation-to-log-Likelihood (CC-to-log(L)) mapping (e.g.
Brogi & Line, 2019; Gibson et al., 2020). Both methods have different uses,
and we use both in this work to characterise the composition, structure, and
dynamics of the WASP-33 b atmosphere, as described below.

####  3.2.1 S/N method

The S/N-method cross-correlates each observed spectrum with a spectral
template and the S/N is calculated from the maximum and standard deviation of
the corresponding Cross-Correlation Function (CCF). First the template model
is scaled to the expected star-planet contrast ratio using

${F_{\rm{model,scaled}}=\frac{F_{\rm{model}}}{B_{\rm{\star}}}\frac{{R_{\rm{p}}}^{2}}{{R_{\star}}^{2}}},$
(1)

with planet radius $R_{\rm{p}}$, stellar radius $R_{\rm{star}}$ and assuming a
stellar black body flux $B_{\rm{\star}}$ at the effective stellar temperature
(using values in Table 4). A stellar black body is a reasonable assumption for
the host star WASP-33 given its hot effective temperature which results in few
spectral stellar absorption lines in the wavelength region of interest.
Previous works indicate that it is important to match the line shape correctly
for HRCCS (e.g. Spring et al., 2022), hence before cross-correlating, we
convolve each synthetic spectrum of WASP-33 b to the median observed resolving
power for each night (see Fig. 4), using a Gaussian instrumental profile.
Since we applied a 6th-order Butterworth high-pass filter at a cutoff-
frequency of $0.02\ \rm{pixel}^{-1}$ to our data, we also apply it to the
scaled model template, where we adjust the cutoff-frequency to account for the
difference in sampling rate between the model and data wavelength grid. We
calculate the correlation coefficient as defined by Brogi & Line (2019) for
each observed residual for each spectral order for all observing nights with
the scaled template. We explore a radial velocity range of [-500, 500] km/s in
steps of $\Delta v=5\ \rm{km/s}$, close to the theoretical instrumental
resolution of ARIES. This results in a CCF time series for each spectral
order, we will refer to this as the CCF-matrix (see Figure 7). To robustly
detect if the emission source is indeed the planet, it is common practice to
shift each CCF-matrix using a grid of trial velocities into the proposed
planet rest frame velocity (see right panel of Figure 7):

${v_{\rm{trial}}=v_{\rm{bary}}+v_{\rm{sys}}+K_{\rm{p}}\sin{\phi_{\rm{p}}}},$
(2)

with $v_{\rm{bary}}$ the barycentric velocity, $v_{\rm{sys}}$ the trial system
velocity, trial $K_{\rm{p}}$ the semi-major amplitude of the planet’s orbital
velocity and $\phi_{\rm{p}}$ the planet’s orbital phase, where we have assumed
a circular Keplerian orbit. The Python package
Barycorrpy888https://github.com/shbhuk/barycorrpy (Wright & Eastman, 2014;
Kanodia & Wright, 2018) is used to calculate $v_{\rm{bary}}$ at the observed
times. We explore a $31\times 31$ grid in $(v_{\rm{sys}},K_{\rm{p}})$-space
around the expected values of ($v_{\rm{sys}}$, $K_{\rm{p}}$) = (-0.3, 230.9)
km/s as previously found in the detection of OH by Nugroho et al. (2020b)
within the range of $\pm 100$ km/s for $v_{\rm{sys}}$ and $\pm 150$ km/s for
$K_{\rm{p}}$ on a 31x31 grid. For each trial velocity on the grid, we shift
the CCF-matrix to $v_{\rm{trial}}$ and calculate the combined CCF by taking
the mean along the time axis. The S/N is determined by calculating the maximum
of the combined CCF and dividing it by the standard deviation of the CCF.
Those exposures taken in between the first and last contact of the eclipse of
WASP-33 b (see Fig. 1) are excluded from the S/N calculation to avoid dilution
of the combined CCF signal.

Figure 7: Cross-Correlation Functions (CCFs) for the best-modified P-T profile
detected at a S/N of 7.9$\sigma$. Top-left: The CCF for each frame including
all three observing nights. Frames are sorted by orbital phase but are not
equally spaced, and the blue dotted line marks the post-eclipse frame with
symmetric phase to the earliest pre-eclipse frame (see Figure 1). The
start/end of ingress/egress are indicated by the white dashed lines labeled
TI-TIV for the four contact points. The two red lines are centered around the
planet radial velocity trail. Top-right: Same as the top-left panel, but
shifted to the planet rest frame. Bottom-left: the mean CCF along the y-axis.
Frames where the planet is not fully visible (in between TI-TIV) were excluded
when computing the CCF to prevent dilution of the planet signal. Bottom-right:
same as the bottom-right panel, but shifted to the planet’s rest frame. Figure
8: Weights per spectral order after cross-correlation with the best self-
consistent PHOENIX model used in our S/N-method adopted. Weights are
calculated from the S/N of the difference between the injected and observed
CCF, indicated by the miniature CCF above each order. Orders 7-9 and 11 were
excluded from the analysis in the wavelength calibration step. Orders 23-25
around the CO 2.3$\micron$ spectral feature are naturally assigned the highest
weights.

Not all spectral orders are equally sensitive when cross-correlating with the
template model. Lower sensitivity is expected for the bluer ARIES orders due
to poorer data quality due to lower throughput and stronger spectral fringes
and weaker spectral features in the models. To account for both, we use an
order-specific weighting scheme in the S/N-method to combine the CCF-matrices,
based on the data quality and model. To calculate the weights we follow Spring
et al. (2022), where we first inject the model (before PCA cleaning) into the
observations using the expected model strength and orbital parameters (based
on Nugroho et al. (2020b)). We create a CCF-matrix for each injected order,
and then subtract its corresponding observed (no injection) CCF-matrix. For
each order we calculate the S/N of this CCF difference. For each order $n$,
weights $w_{n}$ are assigned proportional to this S/N and normalised such that
$\sum^{N}_{n=1}w_{n}=1$, where $N$ is the total number of spectral orders
included. An example of the weighting scheme is shown in Fig. 8. As expected,
using the self-consistent PHOENIX model template at $\times$1 Solar
metallicity with a re-radiation factor of $f=1/2$ as an example, the reddest
orders around the $\sim$2.3 micron densely-packed CO spectral feature region
contribute the highest weights due to the numerous strong spectral lines and
higher S/N of the observed spectra.

####  3.2.2 CC-to-log(L) mapping and model scaling parameter

While the S/N-method is widely used in literature and thus a convenient way to
compare the significance of our detection with previous works, it is
inherently normalised by definition and therefore insensitive to scaling of
the model or observed spectra. This is why there was a need to implement a
weighting scheme to combine multiple spectral orders in the S/N-method. The
CC-to-log L Bayesian Framework described by Brogi & Line (2019) however does
not suffer from this. CC-to-log(L) mapping maximizes the log-likelihood as:

${\log{L}=-\frac{N}{2}\log{\big{(}s_{g}^{2}-2R+s_{f}^{2}}\big{)}}$ (3)

with $s_{\rm{g}}$ and $s_{\rm{f}}$ the variance of model and data,
respectively, $R$ the cross-covariance and $N$ the number of spectral
channels. When combining multiple nights and spectral orders, we can compute
the total sum of their log-likelihood values. There is thus no need to weight
each spectral order as this is accounted for in the CC-to-log(L) mapping
(Brogi & Line, 2019). Moreover, CC-to-log(L) mapping gives us a framework to
extract uncertainties on the measured and derived quantities from the
modelling.

Our implementation uses
Pymultinest999https://johannesbuchner.github.io/PyMultiNest/ (Buchner et al.,
2014b), a generic Python package connected to MultiNest, a Bayesian inference
tool101010https://github.com/JohannesBuchner/MultiNest (Feroz & Hobson, 2008;
Feroz et al., 2009; Feroz et al., 2019), for all our parameter estimations.
Importantly, following Brogi & Line (2019), we also introduce one additional
free parameter, the scaling factor $a$, which allows us to compensate for any
potential unknown scaling of the model template, but ideally should be
retrieved at $a=1$. We multiply the template by $a$ after applying the high-
pass filter. The multidimensional parameter space is explored using the priors
in Table 2.

Parameter | Symbol | Prior
---|---|---
System velocity offset | $\Delta v_{\rm{sys}}$ | Uniform(-50,50) km/s
Orbital velocity offset | $\Delta K_{\rm{p}}$ | Uniform(-50,50) km/s
Scaling parameter | $a$ | LogUniform(-2,2)
Table 2: Priors used in CC-to-log(L) framework. Velocity offsets are given
relative to the previous WASP-33 b velocity values of ($v_{\rm{sys}}$,
$K_{\rm{p}}$) = (-0.3, 230.9) km/s as reported by Nugroho et al. (2021).

To statistically compare the models amongst each other we apply Wilks’ Theorem
(Wilks, 1938). In practical terms, this theorem states that the difference in
$\log(L)$-values between two models with $n$ free parameters follows a
$\chi^{2}$-distribution with $n$ degrees of freedom (for a concise statistical
description see Pino et al., 2020). The confidence interval p-value can now be
calculated from the corresponding $\chi^{2}$-distribution and converted to a
familiar $\sigma$-value. This way, by its construction, we can only compare
models relative to the best trial model.

##  4 Results

We now present the results using both the S/N-method and the CC-to-log(L)
mapping for our four different modelling suites described in Section 3.1. We
detect the atmosphere of WASP-33 b and measure its properties to varying
significance and confidence for the best matching model in each modelling
suite. The key results are summarised in Table 3. The best-matching model from
each suite results in the detection the Doppler-shifting spectrum of the
planet, as demonstrated by the radial velocity trail in the CCF-matrices shown
in Figures 7 and 17. Notably, the trail disappears at the calculated start
time of secondary eclipse ingress and reappears at the end of its egress,
confirming that the detected signal is associated with the planet’s spectrum
as it is obscured by the star during secondary eclipse.

Model | Highest S/N | $\Delta\log(a)$
---|---|---
|  | (symmetric phase coverage)
PHOENIX (1D) modified P-T structure models | 7.9 | $0.21^{+0.03}_{-0.03}$
PHOENIX (1D) self-consistent models | 7.6 | $0.24^{+0.03}_{-0.03}$
GCM (3D) day-side-only models (fixed orbital phase $\phi=0.5$) | 7.1 | $0.27^{+0.03}_{-0.04}$
GCM (3D) phase-dependent models | 7.1 | $0.13^{+0.03}_{-0.03}$
Table 3: Overview of the highest S/N-ratio of the sets of models explored in
this work and their corresponding scaling parameter $a$ from the CC-to-log(L)
mapping. The final column shows $\Delta\log(a)$ i.e. the difference between
the mean of the pre-eclipse (Nights 1 & 3) scaling parameter and that for the
post-eclipse data subset for which we have symmetric phase coverage. The
$\log(a)$ scaling parameter has a value of zero when the model is a good
description of the data. The phase-dependent GCM gives the smallest
discrepancy in $\log(a)$, indicating that it is accounting for the phase-
dependence of the scaling parameter.

###  4.1 Thermal structure and composition

The model that results in the highest S/N detection is a modified PHOENIX
model that includes all opacity sources, with the highest likelihood
corresponding to at a modified P-T profile with $T_{3}=2000\textrm{K}$,
$P_{1}=10^{-3}\ \textrm{bar}$, $P_{3}=10^{-6}\ \textrm{bar}$,
$\alpha_{2}=0.17$ and a metallicity of $10\times$ Solar. As shown in the top
row of Fig. 9 this model gives $S/N=7.9\sigma$ when all nights are combined,
and is higher during pre-eclipse orbital phases than post-eclipse. The
correspondingly similar results for the PHOENIX self-consistent models are
shown in Fig. 18. Although the S/N of the detection varies with phase, this
alone is not sufficient evidence to interpret it as due to a planetary feature
e.g. an offset hot spot, day-to-night variations in brightness, abundances, or
thermal structure. However, we discuss in Section 5.3 the additional evidence
supplied by the CC-to-log(L) mapping scaling parameter that could support such
interpretations.

Figure 9: S/N as a function of the system velocity $v_{\rm{sys}}$ and orbital
velocity $K_{\rm{p}}$. The gray plus symbol marks the maximal S/N location.
The dashed white grid indicates the expected location at ($v_{\rm{sys}}$,
$K_{\rm{p}}$) = (-0.3, 230.9) km/s based on the OH detection by Nugroho et al.
(2021). Our detection is in agreement with their ($v_{\rm{sys}}$,
$K_{\rm{p}}$)-values for all nights combined. Top row: the highest S/N-ratio
modified P-T structure PHOENIX model detected at 7.9$\sigma$. Bottom row: the
phase-dependent GCM model detected at 7.1$\sigma$.

The top-left panel of Fig. 10 shows the results from the log(L)-to-CC mapping
for the best-matching modified PHOENIX model, giving the posterior
distributions for the systemic velocity, orbital velocity, and the scaling
parameter. The planet is detected at $v_{\rm{sys}}=0.15^{+0.64}_{-0.65}\
\rm{km/s}$ and $K_{\rm{p}}=229.53^{+1.11}_{-1.02}\ \rm{km/s}$ with this
highest S/N model, in agreement with previous literature e.g. ($v_{\rm{sys}}$,
$K_{\rm{p}}$) = [-0.3, 230.9] km/s (Nugroho et al., 2021). The scale factor
posterior distribution, which would have $\log(a)=0$ for a perfectly matched
model, is systematically offset to smaller values indicating that the model
does not fully encapsulate all of the physics and chemistry to describe the
atmosphere. Similar is seen in the top right panel of Fig. 10 for the best-
matching self-consistent PHOENIX model. However, all nights are in agreement
with their ($v_{\rm{sys}}$, $K_{\rm{p}}$)-values within 2$\sigma$. We discuss
the variation of $\log(a)$ in later sections.

Figure 10: Marginalized distributions from the CC-to-log(L) framework based on
Brogi & Line (2019). Results are shown for the best-matching PHOENIX modified
P-T structure model (top-left panel), the best-matching self-consistent
PHOENIX model (top-right panel), the GCM day-side only model (phase fixed at
$\phi=0.5$; bottom-left panel) and the GCM phase-dependent model (bottom-right
panel). The contours indicate the 1$\sigma$, 2$\sigma$ and 3$\sigma$
confidence intervals. All frames with phases during the eclipse have been
excluded from this analysis to prevent dilution of the planetary signal.
Results are shown for individual observing nights, all nights combined, and
for the subset of post-eclipse phases that are symmetrically matched with the
pre-eclipse phases (see Fig. 1). The dashed black grid indicates the expected
location at ($v_{\rm{sys}}$, $K_{\rm{p}}$) = (-0.3, 230.9) km/s based on the
OH detection by Nugroho et al. (2021).

We compare the best-matching modified P-T structure PHOENIX model with the
rest of the models in the suite, as well as that of the best-matching self-
consistent PHOENIX model, to give confidence intervals with respect to the
alternative P-T profiles, fixed at 10 $\times$ Solar metallicity. This was
calculated using Wilks’ Theorem (with 7 free parameters), and is shown in Fig.
11. The best-matching modified P-T structure is slightly cooler than the best-
matching self-consistent PHOENIX model, and inverted P-T profiles with an
upper and lower atmosphere temperature difference close to the best-matching
model are favoured. Both can be understood by the fact that we are sensitive
to the the relative line strength with respect to the continuum, which is set
by the upper and lower atmosphere temperature contrast (or equivalently the
strength of the inversion layer). On the contrary, we are less sensitive to
the absolute temperature due to loss of the continuum information in the HRCCS
processing, resulting in a range absolute lower and upper atmosphere
temperatures inside the calculated 1$\sigma$-confidence interval.

Figure 11: Confidence intervals for the PHOENIX modified structure models at
$10\times$ Solar metallicity. These are for all observing nights combined and
were calculated using Wilks’ Theorem. The P-T profile of the best modified
structure model, the one with the highest log-likelihood, is plotted in dark
green and has the following parameterisation: $T_{3}=2000\ \textrm{K}$,
$P_{1}=10^{-3}\ \textrm{bar}$, $P_{3}=1\ \textrm{bar}$, $\alpha_{2}=0.17$. All
confidence intervals here are with respect to this best model. The best self-
consistent PHOENIX model’s P-T profile, at $1\times$ Solar metalictiy and with
a heat reradiation factor $f=1/2$, is plotted in red for comparison.

To determine if any one particular species was contributing the majority of
the detected signal, we fix the P-T profile to the best-matching modified
PHOENIX model (all opacities included), and then run the S/N-analysis for each
individual night without CO, $\rm{H}_{2}O$ and OH, respectively. The results
are shown in Fig. 12. Exclusion of CO significantly impacts the S/N on every
night, and drops from $7.9\sigma$ to $4.3\sigma$ for all nights combined. On
the other hand, exclusion of OH and $\rm{H}_{2}O$ only marginally reduces the
measured S/N, indicating that we are not sensitive to these species at their
abundance with this data. However, on the first night of pre-eclipse data, we
still detect the planet signal to relatively high S/N, despite the exclusion
of CO (the corresponding S/N matrices are shown in Fig. 19). This suggests
that other species in the model are summing to give a detection, but we do not
identify these individual species. The data on the third night, which covers
the same orbital phase, is of lower quality and may explain why a similar
feature is not seen during this night.

Figure 12: S/N results for the models with different opacity sources included
or excluded. Left panel: results from the best modified P-T structure PHOENIX
model for the cases of all opacity sources included and with a single opacity
source excluded, respectively models with OH, $\rm{H}_{2}O$ and CO excluded.
All these models fix the P-T profile to the best modified structure PHOENIX
model where all opacity sources were included. Only the models without CO
result in a significant drop of the S/N when compared to the all included
case. Right panel: Results from the GCM day-side-only opacity templates for
the cases of all opacity sources included and respectively with only CO, only
$\rm{H}_{2}O$ and only OH included. From all opacity sources, only CO remains
robustly detected.

While the modified PHOENIX models allow us to explore the structure and
chemical species of the atmosphere of WASP-33 b, the self-consistent PHOENIX
models allow us to assess its heat-redistribution efficient and metallicity.
The highest log-likelihood in the PHOENIX self-consistent model suite
corresponds to a day-side-only heat-redistribution efficiency $f=0.5$ and a
metallicity of 1$\times$ Solar. Confidence intervals comparing the self-
consistent PHOENIX models computed using Wilks’ Theorem (with 4 free
parameters) are shown in Fig. 13, where we use a visualisation similar to that
recently shown by Giacobbe et al. (2021). Models with high re-radiation
factors $f\geq 1/2$ and high metallicity $\geq 1\times$ Solar are favoured and
within 1$\sigma$ of our best model. We note that the preference for the
modified PHOENIX models is not surprising given the additional free parameters
in the model and the larger number of models explored in the suite.

Figure 13: Confidence intervals for the PHOENIX self-consistent model grid for
all observing nights combined using Wilks’ Theorem. The best model, the one
with the highest log-likelihood, is indicated by the star-marker. By
definition of Wilks’ Theorem, all confidence intervals here are with respect
to this best model. Models with $f\geq 1/2$ and metallicity $\geq 1\times$
Solar are all within 1$\sigma$.

###  4.2 Atmospheric dynamics and longitudinal variation

The post-eclipse phase coverage on Night 2 is longer than the pre-eclipse
coverage on Nights 1 and 3. To enable a direct comparison, we split the post-
eclipse Night 2 spectra into a subset of data that contains only the symmetric
phases corresponding to the pre-eclipse coverage. A comparison of results from
Night 2 and this subset of data with the different model suites are shown in
Fig. 10 by the grey and green contours respectively, and in the right column
of the S/N-maps in Fig. 9. The S/N for the subset data is reduced, but only by
$\Delta$S/N=0.7 at most, indicating that the data containing more night side
hemisphere contributes less to the overall S/N of the detection. The green
contours for the subset of phases are also tighter and better constrained in
all cases in the $v_{sys}-K_{p}$ plots, again indicating a better match with
the model when the night hemisphere contributions in the data set are reduced.
As shown in Table 3, the scaling parameter from the CC-to-Log(L) mapping
differs for the symmetric pre- and post-eclipse phases, indicating a phase-
dependence in the observed spectra of the planet due to asymmetry in the disk
of the planet caused by e.g. an offset hot spot. To qualify this further, we
use the GCM models described in Section 3.1.3.

To first check for general consistency between the GCM models and the results
from the PHOENIX 1D models, we perform our two analysis methods using the day-
side-only (i.e. no phase dependence) 3D GCM model with all opacities included.
This is detected at $S/N=7.1\sigma$, which is slightly lower the PHOENIX
models, and found at the expected ($v_{\rm{sys}}$, $K_{\rm{p}}$)-values for
all nights combined (see bottom row of Fig. 18 and the corresponding cross-
correlation trail is shown in the middle row of Fig. 17). The lower S/N is not
necessarily surprising as we explore only a single GCM day-side-only model vs.
many PHOENIX models, and we know that the GCM lacks some of the important
optical absorbers in its energy balance e.g. atomic iron, such that it cannot
produce as strong temperature inversions. We can also assess in the GCM day-
side-only model if any single species contributes the majority of the
detection. The S/N results for single species models are shown in the right
panel of Fig. 12 and the corresponding S/N-plots are shown in Fig. 20. Only CO
has significant S/N at $6.6\sigma$. The GCM H2O model has $S/N=3.3\sigma$,
which is insufficient for a robust detection in HRCCS as previous work
demonstrated spurious signals at a SNR $\leq$4$\sigma$ often persist Cabot et
al. (2020); Spring et al. (2022), even though it appears at the expected
planet velocity. OH is not detected. The CC-to-log(L) analysis shown in Fig.
21 further highlights the lack of robust detection of any species by CO with
this model, in general agreement with the PHOENIX models.

Our main goal is to understand possible phase dependence in $\log(a)$. The
bottom-left panel of Fig. 10 shows the results from the CC-to-log(L) framework
for the GCM day-side-only model. The results are qualitatively similar to the
the PHOENIX models, again with offsets in $\log(a)$ for the different phase
ranges, but with larger error margins. We hypothesize these larger error
contours may be due to the broadening of the lines by planet rotation in the
GCM model and its overall lower S/N detection. The broad agreement with the
trend for offsets in $\log(a)$ is confirmed by the day-side-only GCM and thus
we proceed to allow the GCM to have phase dependence and determine if this
resolves the $\log(a)$ scaling discrepancy.

The GCM phase-dependent model provides a different cross-correlation template
for each observed spectrum at each phase, where the spectrum corresponds to a
P-T profile that a combination of all the profiles from different longitudes
visible across the planet disk at that time. The GCM phase-dependent model is
similarly detected at $S/N=7.1\sigma$ at the expected ($v_{\rm{sys}}$,
$K_{\rm{p}}$)-location for all night combined (the corresponding CCF-matrix
and $K_{p}-v_{sys}$ maps are shown in the bottom rows of Fig. 17 and Fig. 9,
respectively). However, the results from the CC-to-log(L) mapping shown in the
bottom-right panel Fig. 10 are notably different for the scaling parameter.
The phase-dependent GCM appears to resolve the majority of the $\log(a)$
discrepancy between the pre- and post-eclipse symmetric phases, and brings all
the data sets into broad agreement within $2\sigma$. We discuss this further
in Section 5.3. Although the GCM gives $\log(a)=0$ closer to 1 than the
PHOENIX 1D models, indicating that the inclusion of 3D dynamical effects such
as hot spots and orbital phase brightness variation are a better match to the
data, it does not fully describe the atmosphere of WASP-33 b. Additional
optical absorbers not in the GCM, such as Fe, Fe+ and SiO would enable greater
thermal inversion (Lothringer et al., 2020) likely leading to even deeper
lines. While missing physics in the GCM may account for some of the
$log(a)<0$, it may be that the data reduction process also impacts the overall
line shape as discussed in the next section.

##  5 Discussion

###  5.1 Thermal structure and composition

The detection of CO emission lines with the modified PHOENIX P-T structures
provide unambiguous evidence for a thermal inversion in the atmosphere of
WASP-33 b. A non-inverted atmosphere would have exhibited absorption lines
instead. This is in agreement with previous observations of other atomic and
molecular lines in emission in the atmosphere of WASP-33 b (e.g. Deming et
al., 2012; Haynes et al., 2015; Zhang et al., 2018; Nugroho et al., 2020b;
Cont et al., 2021a, b; Herman et al., 2022) and expected from atmospheric
modelling of UHJ atmospheres (e.g. Lothringer et al., 2018; Gandhi &
Madhusudhan, 2019). Although, we do not explore non-inverted or isothermal P-T
structures, our highest log-likelihood P-T structure (see Fig. 11) shows
temperature contrasts with well constrained upper- and lower limits when
including a thermal inversion layer.

From all molecular templates investigated in this work, only CO has resulted
in a robust detection. This is further supported by the CC-to-log(L) results
shown in Fig. 21 where only CO has similar constraints on the velocities
compared to when all opacity sources are included:
$v_{\rm{sys}}=0.75^{+0.81}_{-0.79}\ \rm{km/s}$,
$K_{\rm{p}}=228.74^{+1.26}_{-1.29}\ \rm{km/s}$. Furthermore, only the
exclusion of CO had a significant impact on the detection strength of the
planet signal. It does remain somewhat puzzling however that we still detect
the best modified PHOENIX model without CO at a $\textrm{S/N}=5.8$ in the
first night of pre-eclipse data. This suggests we still detect the combined
set of all other opacity sources, besides CO, during the first night, even
though we cannot attribute these to $\rm{H}_{2}O$ or OH alone, and we do not
expect significant contributions from iron lines in the ARIES wavelength range
either. We do not detect the model without CO during the third night, despite
covering similar orbital phases. This may be explained by the better observing
conditions during the first observing night relative to the third observing
night (see S/N measurements in Table 1). Although there have been many
detections of CO in absorption (Snellen et al., 2010; Brogi et al., 2012;
Rodler et al., 2012; de Kok et al., 2013; Rodler & et al., 2013; Lockwood et
al., 2014; Brogi et al., 2014, 2016b), this is the first robust detection of
CO emission lines using the HRCCS technique, utilizing only the large Doppler-
shift induced by the planet’s orbital motion.

Previous works have indicated both water (Haynes et al., 2015) and OH (Nugroho
et al., 2021) in the atmosphere of WASP-33 b, but we do not make a robust
detection of these molecules in the ARIES data. To demonstrate that the very
weak GCM $\rm{H}_{2}O$ signal in our data is caused by non-planetary residuals
e.g. tellurics, we ran the CC-to-log(L) analysis using only frames obtained
when the planet is not visible i.e. during the full eclipse of WASP-33 b.
Spurious signals due to telluric or stellar residuals should persist whereas
signals originating from WASP-33 b should disappear or move away from the
planet’s expected ($v_{\rm{sys}}$, $K_{\rm{p}}$)-location. For $\rm{H}_{2}O$
we found that a signal is still retrieved at the same offset ($v_{\rm{sys}}$,
$K_{\rm{p}}$)-position albeit with larger errors, while the CO signal was not
obtained at the expected planet ($v_{\rm{sys}}$, $K_{\rm{p}}$)-position
anymore. This supports the interpretation that any $\rm{H}_{2}O$ signal in the
ARIES data is caused by non-planetary residuals.

The non-detection of OH in with ARIES/MMT is not necessarily surprising either
when compared to the OH detection by Nugroho et al. (2021) using Subaru/IRD
data. The bluer wavelength region (0.97-1.75 $\micron$) covered by Subaru/IRD
with respect to ARIES/MMT contains stronger OH emission lines and Subaru/IRD’s
has higher spectral resolution ($R=70,000$). Thus, the lack of an OH detection
with ARIES can likely be explained by the poorer data quality of the spectral
orders covering the OH line dense regions compared to CO line dense regions.
The synthetic GCM spectrum in the left-top panel of Fig. 6 shows OH line dense
regions in the wavelength range of the bluer spectral orders and order 26 at
the far red-end. As explained in Section 3.2.1, lower weights are assigned to
the bluer ARIES orders which have a lower throughput and contain residual
fringing. The reddest order 26 falls at the edge of the detector and is
consequently of poorer data quality. Lastly, telluric absorption is strong in
some OH line dense spectral orders.

Due to the high dissociation temperature for CO, the abundances of CO are
approximately constant in- and outside of the hotter and cooler regions of the
planet, at the pressures probed by HRCCS (Lodders & Fegley, 2002). This makes
CO a particularly interesting molecule to target for the characterisation of
the thermal structures of UHJs, especially compared to atomic species or other
molecules such as $\rm{OH}$, $\rm{H}_{2}O$ or TiO which dissociate on the day-
side, resulting in an additional phase-dependent chemical gradient,
complicating the interpretation of the line contrast. We therefore advocate
for CO as a good probe of the thermal structure and inversions of UHJs
compared to $\rm{H}_{2}O$ in HRCCS.

###  5.2 Robustness of CO as a temperature tracer in the presence of stellar
activity

The CCFs inside the planet radial velocity trail in Fig. 7 are distinct and
easily seen even by eye. We do not see contaminating signatures of the
$\delta$ Scuti pulsations. On the contrary, the contamination by stellar
pulsations can be seen clearly in the cross-correlation analyses neutral iron
emission lines (Nugroho et al., 2020a; Herman et al., 2022). This is not
exclusive to WASP-33 b; any pulsating star with opacity sources present in
both the star and planet will contaminate HRCCS. This contamination limits the
available effective orbital phase coverage, as frames obtained at phases close
to secondary-eclipse have to be either heavily processed to remove the stellar
pulsations (Johnson et al., 2015; Temple et al., 2017; van Sluijs et al.,
2019) or disregarded (Nugroho et al., 2020a; Herman et al., 2022; Spring et
al., 2022). In contrast to atomic species often present in both the planet and
star, CO is dissociated in the hot stellar atmosphere of early-type A/F-stars.
These type of stars are also more likely to be pulsating stars due to their
location within the Hertzprung-Russel diagram’s instability strip (e.g.
Gautschy & Saio, 1996). This highlights a further advantage of using CO in the
NIR to probe the atmospheric temperature structure and inversions of UHJ
around hot pulsating stars.

###  5.3 Atmospheric dynamics and longitudinal variation

Figure 14: Line contrast (flux - mean(flux)) of a single CO line (2.314-2.316
$\micron$) in the GCM phase-dependent model as function of orbital phase
($\phi$) coverage. The phase-dependent variations show shallower lines pre-
eclipse compared to post-eclipse for symmetric phases around the secondary
eclipse. This is due to a shallower average temperature gradient in the pre-
eclipse phases compared to post-eclipse.

Hot Jupiter thermal structure is expected to be asymmetric, with the hottest
point being shifted away from the sub-stellar point. This shift arises from
the competition between heat transport by winds of planetary-scale waves and
the radiative cooling of the parcel of gas (e.g. Showman & Guillot, 2002;
Perez-Becker & Showman, 2013). The vast majority of observations and models
feature a shift of the hot spot in the direction of the planet rotation,
usually defined as an eastward shift. In such a case, the pre-eclipse phases
are more dominated by the hot spot than the post-eclipse phases, due to the
planet rotation during the eclipse.

Previous evidence for a hot spot on WASP-33 b comes from NIR and optical
photometric phase curve observations, where eastward hot spot results in
negative phase curve offset. Zhang et al. (2018) report the first evidence of
an eastward hot spot on WASP-33 b, measuring a phase curve offset of $-12.8\pm
5.8\degree$ in the Spitzer $3.6\micron$-band. Conversely, von Essen et al.
(2020) found a $+28.7\pm 7.1\degree$ westward phase curve offset in the
optical TESS light curve, and a range of theories including
magnetohydrodynamic effects, non-synchronous rotation, and clouds have been
invoked to explain such westward offset hot spots (e.g. Dang et al., 2018;
Hindle et al., 2021). Recent work by Herman et al. (2022) further report a
$+22\pm 8\degree$ westward offset based on the phase offset of neutral iron
thermal emission lines at high spectral resolution. Our CO emission line
detection also supports a phase-dependence in the planet spectrum. However,
the GCM that we use to model it includes an eastward (not westward) hot spot,
which is oriented towards us during the pre-eclipse phases of WASP-33 b.

It is important to note that the phase-dependence detected in this work and in
Herman et al. (2022) is found via a model scaling parameter ($a$ in this work,
and a similar $A_{p}(\phi)$ in equation 4 of the other study). Both studies
detect that a larger scaling is required to best model the observations just
after eclipse. Herman et al. (2022) interpret the larger $A_{p}(\phi)$ post-
eclipse as an overall brightness variation, indicating more flux west from the
sub-stellar point and hence a westward hot spot. This is the typical
interpretation for photometric phase curve measurements, which are sensitive
to absolute flux as a function of phase. Instead, we highlight that because
the HRCCS we use in this work divides out the continuum, that the scaling
parameter is a measure only of line contrast i.e. that the lines in the planet
spectrum are larger just after eclipse, as shown in Figure 14.

The phase-dependent line contrast can then be understood in the context of our
GCM model. First, the eastward hot spot results in a hotter continuum and thus
greater overall brightness for pre-eclipse phases as shown in the bottom-left
panel of Fig. 6. But, since our post-processing of the high resolution spectra
removes the continuum information we argue we are insensitive to these
variations. Fig. 6 then shows that the GCM predicts thermal profiles that are
shallower in the eastward-shifted hot spot and steeper on the western part of
the day-side. Overall, this leads to a disk-integrated temperature gradient
that is shallower during pre-eclipse when it contains the eastward hot spot
(where our line-of-sight is centered on longitudes $\sim 17-61\degree$ east of
the sub-stellar point), and steeper during post-eclipse when it contains more
of the western day-side hemisphere (where our line-of-sight is centered on
longitudes $\sim 277-341\degree$ east of the sub-stellar point). Consequently,
the average line contrast associated with these temperature gradients will be
shallower during pre-eclipse and stronger during post-eclipse. We therefore
conclude that the data can be well described by an eastward offset hot spot.
We do not find significant evidence for a relative blue-shift pre-eclipse
compared to post-eclipse that would further indicate eastward hot spot. Such
net Doppler shifts are predicted of $\sim$1-3 km/s as the hot spot rotates in
and out of the observers view over while orbiting its host star (Zhang et al.,
2017). However, our velocity sampling is $\sim$4-6 km/s for the best self-
consistent PHOENIX model for the individual observing nights. Thus we cannot
resolve such Doppler shifts at our spectral resolution. We note that when
including all orbital phases out to quadrature in the post-eclipse night, the
scaling parameter prefers shallower lines (see grey contours in Fig. 10),
despite the strong line contrasts shown in Fig. 14. While these later phases
give only a small addition to the overall S/N, the change in line contrast may
be explained by missing physics in the GCM e.g. clouds close to quadrature
phases that flatten out the continuum, reducing the overall line contrast.

We further find a systematic bias towards small scaling parameter
$\log_{10}{a}<0$ in all models. Rather than missing physics, we consider that
the data cleaning steps of the HRCCS can affect the line shape and strength of
planetary signal. PCA likely degrades more of the planetary signal at lower
resolution, where the lines spread across multiple pixels. Applying PCA to our
template in every iteration of our CC-to-log(L) PyMultinest implementation is
computationally expensive. As an alternative, we investigated these effects by
injecting the model at $a=1$ into the frames observed during the eclipse of
WASP-33 b. We still found a bias towards $a$-values $<1$ for the recovered
signals, albeit much closer to one. This result supports our hypothesis that
our data cleaning steps cause at least part of the bias towards lower scale
factors. Given that we apply the data processing in a homogeneous way to all
data sets, and while we do not give large weight to the absolute values of
$\log(a)$, we still give weight to the interpretation of the relative
differences. Novel new model-filtering techniques such as the on introduced by
Gibson et al. (2022) may alleviate these issues in future work.

###  5.4 Low resolution limit for HRCCS technique

For the HRCCS technique in the photon-limited regime, to first order the S/N
is given by:

${S/N_{\rm{planet}}=(\frac{S_{\rm{p}}}{S_{\rm{\star}}})\ S/N_{\rm{star}}\
\sqrt{N_{\rm{lines}}}},$ (4)

where $(\frac{S_{\rm{p}}}{S_{\rm{\star}}})$ is the planet-to-star contrast
ratio, $\textrm{S/N}_{\rm{star}}$ the stellar S/N-ratio, and $N_{\rm{lines}}$
the number of spectral lines resolved (Snellen et al., 2015; Birkby, 2018).
For the close-in UHJs, we cannot be aided by High Contrast Imaging to
spatially resolve the exoplanet and increase our planet-to-star contrast
ratio. We thus rely solely on disentangling the exoplanet’s lines from it’s
the Doppler shift induced by its orbital motion. We can thus increase our
exoplanet’s S/N in two ways: Firstly, we can increasing our photon collecting
power such that $\textrm{S/N}_{\rm{star}}$ increases or secondly increase the
total number of spectral lines observed by either expanding our instantaneous
wavelength coverage or our spectral resolving power.

Since S/N is proportional to $\sqrt{N_{\rm{lines}}}$ it is crucial to have a
sufficiently high spectral resolution to resolve more lines. Conservatively,
Birkby (2018) placed a lower limit of $R=25,000$ on the use of the HRCCS
method based on the detection of water vapour in the atmosphere of $\tau$ Böo
with pre-upgraded NIRSPEC1.0/Keck combination reported by (Lockwood et al.,
2014). Other evidence for a detection of the thermal spectrum of HD 88133 b,
also with NIRSPEC1.0/Keck was reported by Piskorz et al. (2016) using a multi-
epoch approach, but later on disputed (Buzard et al., 2021). Our detection
with ARIES/MMT pushes this to a lower limit at $R=14,000-16,000$.

Our detection of WASP-33 b at this lower spectral resolution provides proof-
of-concept HRCCS is still possible at lower resolving power given a
sufficiently large orbital Doppler shift. Although less and diluted spectral
lines are observed at a lower spectral resolution, more photons are collected
per resolution element. Therefore one can trade off spectral resolving power
for photon collecting power. This is particularly interesting for short period
exoplanets where we do not need a high spectral resolving power to resolve
their large orbital Doppler shifts. Examples include other UHJs (for example
KELT-9b) or terrestrial lava worlds (for example CoRoT-7b) as both have
orbital velocities of similar order of magnitude as WASP-33 b.

##  6 Conclusions

In this work we presented the first high resolution spectroscopy observations
from the NIR spectrograph ARIES mounted behind the MMT targeting WASP-33b. We
develop public pre-processing pipeline for ARIES to extract a spectral time
series from the raw detector images. We use the HRCCS technique with PCA to
characterise WASP-33b’s atmosphere using both the S/N-method and CC-to-log(L)
mapping. The primary conclusions of this work are the following:

  1. 1.

We present the first robust detection (7.9$\sigma$) of CO emission lines using
HRCCS in an exoplanet atmosphere. The detection of CO emission lines in the
atmosphere of WASP-33 b is unambiguous evidence of an inverted P-T profile.

  2. 2.

The best PHOENIX self-consistent atmospheric models prefer day-side-only heat
redistribution.

  3. 3.

The modified P-T structure PHOENIX models show we are sensitive to the
temperature gradient of the inversion layer, but less sensitive to the
absolute temperatures of the lower atmosphere.

  4. 4.

The CC-to-log(L) mapping reveals a phase-dependent scaling factor, which shows
that the planet spectral line contrast is greater just after secondary eclipse
than before. We demonstrate that the P-T profiles from a GCM that includes an
eastward offset hot spot can explain the majority of the phase-dependence in
$\log(a)$. We emphasise that inclusion of a scaling parameter in the CC-to-
log(L) mapping can reveal the longitudinal phase-dependent thermal structure
of the atmosphere, even in spectra that lack the resolving power required to
detect the additional induced Doppler shift caused by the offset hot spot.

  5. 5.

We find no evidence of stellar pulsations affecting our detection of CO, in
contrast to previous detections of atomic emission lines in the optical. We
advocate that CO is therefore an advantageous and good probe of the thermal
atmospheric structure of UHJs at pressures seen with HRCCS.

  6. 6.

At a measured resolving power of $R\sim 15,000$, this detection of CO in
emission pushes the HRCCS technique to its lowest resolution limits when using
only the large Doppler-shift induced by the exoplanet’s orbital motion to
disentangle its spectrum from its host star and tellurics. Thus spectral
resolving power can be traded for photon collecting power in cases where the
induced orbital Doppler shift is sufficient move the planet spectrum over
multiple pixels such as expected for close-in terrestrial lava planets and
ultra hot Jupiters.

ARIES is currently being upgraded into the MMT: AO exoPlanet characterization
System (MAPS)111111More information can be found here:
https://www.as.arizona.edu/~ktmorz/maps.html. It will have an updated AO
system, broader instantaneous wavelength coverage and modes at a higher
spectral resolution. A modified version of the ARIES data reduction pipeline
presented in this work may prove helpful in analysing future output from MAPS
as well.

Finally, future low-resolution observations of WASP-33 b with the James Webb
Space Telescope can be combined to with ground-based high-resolution
observations to mitigate the effects of stellar pulsations and give tighter
constraints on the P-T profile and chemical abundances of $\rm{OH}$,
$\rm{H}_{2}O$, $\rm{CO}$, and $\rm{CO}_{2}$ (Beichman & Greene, 2018). UHJ
have the additional benefit that refractory elements can be measured as well
(Lothringer et al., 2021). Combined this will enable calculation of abundance
ratios which can in turn be used to constrain different planet formation
pathways (e.g. Öberg et al., 2011; Madhusudhan, 2012; Cridland et al., 2016;
Cridland et al., 2017, 2019; Khorshid et al., 2021). Finally, combining phase-
dependent low resolution data, for example the recently observed diurnal
variations in the atmosphere of WASP-121b (Mikal-Evans et al., 2022), and high
resolution data using the CC-to-log(L) (Brogi et al., 2016a; Brogi & Line,
2019; Gibson et al., 2020, e.g.) mapping will greatly improve our
understanding of the longitudinal chemical and thermal phase-dependence.

## Acknowledgements

JLB extends special thanks to all the MMT Observatory staff, particularly the
telescope and AO operators, for their assistance and valuable advice in
carrying out the MMT Exoplanet Atmosphere Survey. This research is part of the
exoZoo project that has received funding from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation program
under grant agreement No 805445. Observations reported here were obtained at
the MMT Observatory, a joint facility of the Smithsonian Institution and the
University of Arizona. This work was performed in part under contract with the
Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship
Program executed by the NASA Exoplanet Science Institute.

## Data Availability

All the original data will be made publicly available in the author’s Github
repository.

## References

  * et al. (2009) L. S., et al., 2009, J. Quant. Spectrosc. Radiative Transfer, 110, 533
  * Allard et al. (2011) Allard F., Homeier D., Freytag B., 2011, in Johns-Krull C., Browning M. K., West A. A., eds, Astronomical Society of the Pacific Conference Series Vol. 448, 16th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun. p. 91 (arXiv:1011.5405)
  * Arcangeli et al. (2018) Arcangeli J., et al., 2018, ApJ, 855, L30
  * Armstrong et al. (2016) Armstrong D. J., de Mooij E., Barstow J., Osborn H. P., Blake J., Saniee N. F., 2016, Nature Astronomy, 1, 0004
  * Barber et al. (2006) Barber R., Tennyson J., Harris G., Tolchenov R., 2006, Monthly Notices of the Royal Astronomical Society, 368, 1087
  * Barber et al. (2014) Barber R. J., Strange J. K., Hill C., Polyansky O. L., Mellau G. C., Yurchenko S. N., Tennyson J., 2014, MNRAS, 437, 1828
  * Barman et al. (2001) Barman T. S., Hauschildt P. H., Allard F., 2001, ApJ, 556, 885
  * Barman et al. (2011) Barman T. S., Macintosh B., Konopacky Q. M., Marois C., 2011, ApJ, 733, 65
  * Baxter et al. (2021) Baxter C., et al., 2021, A&A, 648, A127
  * Beichman & Greene (2018) Beichman C. A., Greene T. P., 2018, Observing Exoplanets with the James Webb Space Telescope. Springer International Publishing, p. 85, doi:10.1007/978-3-319-55333-7_85
  * Bell et al. (2019) Bell T. J., et al., 2019, MNRAS, 489, 1995
  * Birkby (2018) Birkby J. L., 2018, arXiv e-prints, p. arXiv:1806.04617
  * Borsa et al. (2021) Borsa F., et al., 2021, A&A, 653, A104
  * Brahm et al. (2017) Brahm R., Jordán A., Espinoza N., 2017, PASP, 129, 034002
  * Brogi & Line (2019) Brogi M., Line M. R., 2019, AJ, 157, 114
  * Brogi et al. (2012) Brogi M., Snellen I. A. G., de Kok R. J., Albrecht S., Birkby J., de Mooij E. J. W., 2012, Nature, 486, 502
  * Brogi et al. (2014) Brogi M., de Kok R. J., Birkby J. L., Schwarz H., Snellen I. A. G., 2014, A&A, 565, A124
  * Brogi et al. (2016a) Brogi M., Line M., Bean J., Désert J.-M., Schwarz H., 2016a, preprint,  (arXiv:1612.07008)
  * Brogi et al. (2016b) Brogi M., de Kok R. J., Albrecht S., Snellen I. A. G., Birkby J. L., Schwarz H., 2016b, ApJ, 817, 106
  * Buchner et al. (2014a) Buchner J., et al., 2014a, A&A, 564, A125
  * Buchner et al. (2014b) Buchner J., et al., 2014b, A&A, 564, A125
  * Buzard et al. (2021) Buzard C., Piskorz D., Lockwood A. C., Blake G., Barman T. S., Benneke B., Bender C. F., Carr J. S., 2021, arXiv e-prints, p. arXiv:2109.13275
  * Cabot et al. (2020) Cabot S. H. C., Madhusudhan N., Welbanks L., Piette A., Gandhi S., 2020, MNRAS, 494, 363
  * Cauley et al. (2021) Cauley P. W., Wang J., Shkolnik E. L., Ilyin I., Strassmeier K. G., Redfield S., Jensen A. G., 2021, AJ, 161, 152
  * Chakrabarty & Sengupta (2019) Chakrabarty A., Sengupta S., 2019, AJ, 158, 39
  * Charbonneau et al. (2008) Charbonneau D., Knutson H. A., Barman T., Allen L. E., Mayor M., Megeath S. T., Queloz D., Udry S., 2008, ApJ, 686, 1341
  * Chiavassa & Brogi (2019) Chiavassa A., Brogi M., 2019, A&A, 631, A100
  * Coles et al. (2019) Coles P. A., Yurchenko S. N., Tennyson J., 2019, MNRAS, 490, 4638
  * Collier Cameron et al. (2010) Collier Cameron A., et al., 2010, MNRAS, 407, 507
  * Cont et al. (2021a) Cont D., et al., 2021a, arXiv e-prints, p. arXiv:2105.10230
  * Cont et al. (2021b) Cont D., et al., 2021b, arXiv e-prints, p. arXiv:2112.10461
  * Cowan et al. (2007) Cowan N. B., Agol E., Charbonneau D., 2007, MNRAS, 379, 641
  * Cridland et al. (2016) Cridland A. J., Pudritz R. E., Alessi M., 2016, MNRAS, 461, 3274
  * Cridland et al. (2017) Cridland A. J., Pudritz R. E., Birnstiel T., Cleeves L. I., Bergin E. A., 2017, MNRAS, 469, 3910
  * Cridland et al. (2019) Cridland A. J., van Dishoeck E. F., Alessi M., Pudritz R. E., 2019, A&A, 632, A63
  * Crossfield et al. (2010) Crossfield I. J. M., Hansen B. M. S., Harrington J., Cho J. Y. K., Deming D., Menou K., Seager S., 2010, ApJ, 723, 1436
  * Dang et al. (2018) Dang L., et al., 2018, Nature Astronomy, 2, 220
  * Debras et al. (2020) Debras F., Mayne N., Baraffe I., Jaupart E., Mourier P., Laibe G., Goffrey T., Thuburn J., 2020, A&A, 633, A2
  * Deming et al. (2012) Deming D., et al., 2012, ApJ, 754, 106
  * Dobbs-Dixon & Lin (2008) Dobbs-Dixon I., Lin D. N. C., 2008, ApJ, 673, 513
  * Dobbs-Dixon et al. (2010) Dobbs-Dixon I., Cumming A., Lin D. N. C., 2010, ApJ, 710, 1395
  * Ehrenreich et al. (2020) Ehrenreich D., et al., 2020, Nature, 580, 597
  * Evans et al. (2016) Evans T. M., et al., 2016, ApJ, 822, L4
  * Evans et al. (2017) Evans T. M., et al., 2017, Nature, 548, 58
  * Feroz & Hobson (2008) Feroz F., Hobson M. P., 2008, MNRAS, 384, 449
  * Feroz et al. (2009) Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601
  * Feroz et al. (2019) Feroz F., Hobson M. P., Cameron E., Pettitt A. N., 2019, The Open Journal of Astrophysics, 2, 10
  * Fortney et al. (2021) Fortney J. J., Dawson R. I., Komacek T. D., 2021, Journal of Geophysical Research (Planets), 126, e06629
  * Gandhi & Madhusudhan (2019) Gandhi S., Madhusudhan N., 2019, MNRAS, 485, 5817
  * Garhart et al. (2020) Garhart E., et al., 2020, AJ, 159, 137
  * Gautschy & Saio (1996) Gautschy A., Saio H., 1996, ARA&A, 34, 551
  * Giacobbe et al. (2021) Giacobbe P., et al., 2021, Nature, 592, 205
  * Gibson et al. (2020) Gibson N. P., et al., 2020, MNRAS, 493, 2215
  * Gibson et al. (2022) Gibson N. P., Nugroho S. K., Lothringer J., Maguire C., Sing D. K., 2022, arXiv e-prints, p. arXiv:2201.04025
  * Goorvitch (1994) Goorvitch D., 1994, The Astrophysical Journal Supplement Series, 95, 535
  * Grimm et al. (2021) Grimm S. L., et al., 2021, ApJS, 253, 30
  * Gullikson et al. (2014) Gullikson K., Dodson-Robinson S., Kraus A., 2014, AJ, 148, 53
  * Hargreaves et al. (2019) Hargreaves R., Gordon I., Kochanov R., Rothman L., 2019, in EPSC-DPS Joint Meeting 2019. pp EPSC–DPS2019–919
  * Hargreaves et al. (2020) Hargreaves R. J., Gordon I. E., Rey M., Nikitin A. V., Tyuterev V. G., Kochanov R. V., Rothman L. S., 2020, ApJS, 247, 55
  * Harrington et al. (2006) Harrington J., Hansen B. M., Luszcz S. H., Seager S., Deming D., Menou K., Cho J. Y. K., Richardson L. J., 2006, Science, 314, 623
  * Hauschildt et al. (1997) Hauschildt P. H., Baron E., Allard F., 1997, ApJ, 483, 390
  * Hauschildt et al. (1999) Hauschildt P. H., Allard F., Baron E., 1999, ApJ, 512, 377
  * Haynes et al. (2015) Haynes K., Mandell A. M., Madhusudhan N., Deming D., Knutson H., 2015, ApJ, 806, 146
  * Herman et al. (2020) Herman M. K., de Mooij E. J. W., Jayawardhana R., Brogi M., 2020, AJ, 160, 93
  * Herman et al. (2022) Herman M. K., de Mooij E. J. W., Nugroho S. K., Gibson N. P., Jayawardhana R., 2022, arXiv e-prints, p. arXiv:2203.11188
  * Herrero et al. (2011) Herrero E., Morales J. C., Ribas I., Naves R., 2011, A&A, 526, L10
  * Hindle et al. (2021) Hindle A. W., Bushby P. J., Rogers T. M., 2021, ApJ, 922, 176
  * Hoeijmakers et al. (2018) Hoeijmakers H. J., et al., 2018, Nature, 560, 453
  * Hoeijmakers et al. (2019) Hoeijmakers H. J., et al., 2019, A&A, 627, A165
  * Hubeny et al. (2003) Hubeny I., Burrows A., Sudarsky D., 2003, ApJ, 594, 1011
  * Jackson et al. (2019) Jackson B., Adams E., Sandidge W., Kreyche S., Briggs J., 2019, AJ, 157, 239
  * John (1988) John T. L., 1988, A&A, 193, 189
  * Johnson et al. (2015) Johnson M. C., Cochran W. D., Collier Cameron A., Bayliss D., 2015, ApJ, 810, L23
  * Kanodia & Wright (2018) Kanodia S., Wright J., 2018, Research Notes of the AAS, 2, 4
  * Karman et al. (2019) Karman T., et al., 2019, Icarus, 328, 160
  * Kesseli et al. (2020) Kesseli A. Y., Snellen I. A. G., Alonso-Floriano F. J., Mollière P., Serindag D. B., 2020, AJ, 160, 228
  * Khorshid et al. (2021) Khorshid N., Min M., Désert J. M., Woitke P., Dominik C., 2021, arXiv e-prints, p. arXiv:2111.00279
  * Kitzmann et al. (2018) Kitzmann D., et al., 2018, ApJ, 863, 183
  * Knutson et al. (2007) Knutson H. A., et al., 2007, Nature, 447, 183
  * Knutson et al. (2009) Knutson H. A., et al., 2009, ApJ, 690, 822
  * Kong et al. (2013) Kong H., Akakin H. C., Sarma S. E., 2013, IEEE Transactions on Cybernetics, 43, 1719
  * Kreidberg (2018) Kreidberg L., 2018, Handbook of Exoplanets, p. 2083–2105
  * Kreidberg et al. (2018) Kreidberg L., et al., 2018, AJ, 156, 17
  * Lee et al. (2019) Lee E., Taylor J., Grimm S. L., Baudino J.-L., Garland R., Irwin P. G. J., Wood K., 2019, MNRAS, 487, 2082
  * Lee et al. (2021) Lee E. K. H., et al., 2021, arXiv e-prints, p. arXiv:2110.15640
  * Lehmann et al. (2015) Lehmann H., Guenther E., Sebastian D., Döllinger M., Hartmann M., Mkrtichian D. E., 2015, A&A, 578, L4
  * Li et al. (2015) Li G., Gordon I. E., Rothman L. S., Tan Y., Hu S.-M., Kassi S., Campargue A., Medvedev E. S., 2015, The Astrophysical Journal Supplement Series, 216, 15
  * Line et al. (2021) Line M. R., et al., 2021, Nature, 598, 580
  * Lockwood et al. (2014) Lockwood A. C., Johnson J. A., Bender C. F., Carr J. S., Barman T., Richert A. J. W., Blake G. A., 2014, ApJ, 783, L29
  * Lodders & Fegley (2002) Lodders K., Fegley B., 2002, Icarus, 155, 393
  * Lord (1992) Lord S. D., 1992, in NASA Technical Memorandum 103957.
  * Lothringer & Barman (2019) Lothringer J. D., Barman T., 2019, ApJ, 876, 69
  * Lothringer et al. (2018) Lothringer J. D., Barman T., Koskinen T., 2018, ApJ, 866, 27
  * Lothringer et al. (2020) Lothringer J. D., Fu G., Sing D. K., Barman T. S., 2020, ApJ, 898, L14
  * Lothringer et al. (2021) Lothringer J. D., Rustamkulov Z., Sing D. K., Gibson N. P., Wilson J., Schlaufman K. C., 2021, ApJ, 914, 12
  * Madhusudhan (2012) Madhusudhan N., 2012, ApJ, 758, 36
  * Madhusudhan & Seager (2009) Madhusudhan N., Seager S., 2009, ApJ, 707, 24
  * Mansfield et al. (2018) Mansfield M., et al., 2018, AJ, 156, 10
  * Marsh (1989) Marsh T. R., 1989, PASP, 101, 1032
  * McCarthy et al. (1998) McCarthy D. W., Burge J. H., Angel J. R. P., Ge J., Sarlot R. J., Fitz-Patrick B. C., Hinz J. L., 1998, in Fowler A. M., ed., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 3354, Infrared Astronomical Instrumentation. pp 750–754, doi:10.1117/12.317340
  * Menou & Rauscher (2009) Menou K., Rauscher E., 2009, ApJ, 700, 887
  * Merritt et al. (2021) Merritt S. R., et al., 2021, MNRAS, 506, 3853
  * Mikal-Evans et al. (2022) Mikal-Evans T., et al., 2022, Nature Astronomy,
  * Morzinski et al. (2020) Morzinski K. M., et al., 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series. p. 114481L, doi:10.1117/12.2563178
  * Nugroho et al. (2017) Nugroho S. K., Kawahara H., Masuda K., Hirano T., Kotani T., Tajitsu A., 2017, AJ, 154, 221
  * Nugroho et al. (2020a) Nugroho S. K., Gibson N. P., de Mooij E. J. W., Watson C. A., Kawahara H., Merritt S., 2020a, MNRAS, 496, 504
  * Nugroho et al. (2020b) Nugroho S. K., Gibson N. P., de Mooij E. J. W., Herman M. K., Watson C. A., Kawahara H., Merritt S. R., 2020b, ApJ, 898, L31
  * Nugroho et al. (2021) Nugroho S. K., et al., 2021, ApJ, 910, L9
  * Öberg et al. (2011) Öberg K. I., Murray-Clay R., Bergin E. A., 2011, ApJ, 743, L16
  * Parmentier et al. (2018) Parmentier V., et al., 2018, A&A, 617, A110
  * Perez-Becker & Showman (2013) Perez-Becker D., Showman A. P., 2013, ApJ, 776, 134
  * Pino et al. (2020) Pino L., et al., 2020, ApJ, 894, L27
  * Piskorz et al. (2016) Piskorz D., et al., 2016, ApJ, 832, 131
  * Polyansky et al. (2018) Polyansky O. L., Kyuberis A. A., Zobov N. F., Tennyson J., Yurchenko S. N., Lodi L., 2018, MNRAS, 480, 2597
  * Rauscher & Menou (2010) Rauscher E., Menou K., 2010, ApJ, 714, 1334
  * Rodler & et al. (2013) Rodler F., et al. 2013, MNRAS, 432, 1980
  * Rodler et al. (2012) Rodler F., Lopez-Morales M., Ribas I., 2012, ApJ, 753, L25
  * Sarlot et al. (1999) Sarlot R. J., McCarthy D. W., Burge J. H., Ge J., 1999, in Fischer R. E., Smith W. J., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 3779, Current Developments in Optical Design and Optical Engineering VIII. pp 274–283, doi:10.1117/12.368218
  * Serindag et al. (2021) Serindag D. B., Nugroho S. K., Mollière P., de Mooij E. J. W., Gibson N. P., Snellen I. A. G., 2021, A&A, 645, A90
  * Sheppard et al. (2017) Sheppard K. B., Mandell A. M., Tamburo P., Gandhi S., Pinhas A., Madhusudhan N., Deming D., 2017, ApJ, 850, L32
  * Showman & Guillot (2002) Showman A. P., Guillot T., 2002, A&A, 385, 166
  * Showman & Polvani (2011) Showman A. P., Polvani L. M., 2011, ApJ, 738, 71
  * Showman et al. (2009) Showman A. P., Fortney J. J., Lian Y., Marley M. S., Freedman R. S., Knutson H. A., Charbonneau D., 2009, ApJ, 699, 564
  * Showman et al. (2020) Showman A. P., Tan X., Parmentier V., 2020, Space Sci. Rev., 216, 139
  * Smith et al. (2011) Smith A. M. S., Anderson D. R., Skillen I., Collier Cameron A., Smalley B., 2011, MNRAS, 416, 2096
  * Snellen et al. (2010) Snellen I. A. G., de Kok R. J., de Mooij E. J. W., Albrecht S., 2010, Nature, 465, 1049
  * Snellen et al. (2015) Snellen I., et al., 2015, A&A, 576, A59
  * Spring et al. (2022) Spring E. F., et al., 2022, arXiv e-prints, p. arXiv:2201.03600
  * Stone et al. (2014) Stone J. M., Eisner J. A., Salyk C., Kulesa C., McCarthy D., 2014, ApJ, 792, 56
  * Swain et al. (2009) Swain M. R., Vasisht G., Tinetti G., Bouwman J., Chen P., Yung Y., Deming D., Deroo P., 2009, ApJ, 690, L114
  * Tan & Komacek (2019) Tan X., Komacek T. D., 2019, ApJ, 886, 26
  * Temple et al. (2017) Temple L. Y., et al., 2017, MNRAS, 471, 2743
  * Visscher et al. (2010) Visscher C., Lodders K., Fegley Bruce J., 2010, ApJ, 716, 1060
  * Wilks (1938) Wilks S. S., 1938, The Annals of Mathematical Statistics, 9, 60
  * Wong et al. (2016) Wong I., et al., 2016, ApJ, 823, 122
  * Wright & Eastman (2014) Wright J. T., Eastman J. D., 2014, PASP, 126, 838
  * Yan et al. (2019) Yan F., et al., 2019, A&A, 632, A69
  * Yan et al. (2020) Yan F., et al., 2020, A&A, 640, L5
  * Yan et al. (2021) Yan F., et al., 2021, A&A, 645, A22
  * Yurchenko et al. (2020) Yurchenko S. N., Mellor T. M., Freedman R. S., Tennyson J., 2020, MNRAS, 496, 5282
  * Zhang et al. (2017) Zhang J., Kempton E. M. R., Rauscher E., 2017, ApJ, 851, 84
  * Zhang et al. (2018) Zhang M., et al., 2018, AJ, 155, 83
  * Zucker (2003) Zucker S., 2003, MNRAS, 342, 1291
  * de Kok et al. (2013) de Kok R. J., Brogi M., Snellen I. A. G., Birkby J., Albrecht S., de Mooij E. J. W., 2013, A&A, 554, A82
  * de Mooij et al. (2013) de Mooij E. J. W., Brogi M., de Kok R. J., Snellen I. A. G., Kenworthy M. A., Karjalainen R., 2013, A&A, 550, A54
  * van Sluijs et al. (2019) van Sluijs L., et al., 2019, A&A, 626, A97
  * von Essen et al. (2019) von Essen C., Mallonn M., Welbanks L., Madhusudhan N., Pinhas A., Bouy H., Weis Hansen P., 2019, A&A, 622, A71
  * von Essen et al. (2020) von Essen C., Mallonn M., Borre C. C., Antoci V., Stassun K. G., Khalafinejad S., Tautvaišienė G., 2020, A&A, 639, A34

##  Appendix A Additional tables and figures

For transparency, completeness and reproducibility of this work, we include
some additional figures and tables here.

Quantity | Value | Unit | Reference
---|---|---|---
Stellar mass | 1.495 | $\rm{M}_{\odot}$ | Collier Cameron et al. (2010)
Stellar effective temperature | 7400 | K | Collier Cameron et al. (2010)
Stellar radius | 1.444 | $\rm{R}_{\sun}$ | Collier Cameron et al. (2010)
Primary transit time | 2454590.17936 | BJD | Smith et al. (2011)
Eccentricity | 0.0 | - | Smith et al. (2011)
Semi-major axis | 0.02558 | au | Smith et al. (2011)
Orbital period | 1.21986967 | days | Smith et al. (2011)
Orbital inclination | 87.7 | degree | Lehmann et al. (2015)
Planetary radius | 1.603 | $\rm{R}_{\rm{Jup}}$ | Lehmann et al. (2015)
Impact parameter | 0.21 | - | Chakrabarty & Sengupta (2019)
Expected system velocity | -0.3 | km/s | Nugroho et al. (2020b)
Expected orbital velocity | 230.9 | km/s | Nugroho et al. (2020b)
Table 4: An overview of relevant WASP-33 system parameters used in this work.
Figure 15: Mixing ratios in chemical equilibrium as functions of pressure for
different species in the self-consistent PHOENIX modek atmospheres of WASP-33
b. This is for the day-side-only heat redistribution model $f=0.5$ at $\times
1$ Solar metallicity. Figure 16: Overview of all modified PHOENIX P-T
structures explored in this work.

Figure 17: Similar to Fig. 7, but for the highest SNR PHOENIX self-consistent
model, GCM day-side-only model (middle panel) and GCM phase-dependent model
(bottom panel). Figure 18: Same as Fig. 9 but shown for the best PHOENIX self-
consistent model (top row), best GCM day-side-only model (bottom row). Figure
19: Similar to Fig. 9, but for the PHOENIX modified P-T structure models with
all opacity sources included except CO. The detection of $5.7\sigma$ for the
first night suggests we are still sensitive to some of the other opacity
sources in the atmosphere of WASP-33 b. Figure 20: Similar to Fig. 9, but for
the GCM day-side-only templates with respectively: all opacity sources, only
CO, only $\rm{H}_{2}O$ and only OH included. Amongst the opacity sources, CO
is the only molecule that is robustly detected at 6.6$\sigma$. Figure 21: Same
as the panels in Fig. 10, but results shown for the GCM day-side-only models
for different opacity sources included. All results are shown for all nights
combined.