# A Combined Ground-based and JWST Atmospheric Retrieval Analysis: Both IGRINS
and NIRSpec Agree The Atmosphere of WASP-77A b is Metal-Poor

Peter C. B. Smith  School of Earth and Space Exploration, Arizona State
University, Tempe, AZ 85287  Michael R. Line  School of Earth and Space
Exploration, Arizona State University, Tempe, AZ 85287  Jacob L. Bean
Department of Astronomy and Astrophysics, The University of Chicago, Chicago,
IL 60637  Matteo Brogi  Department of Physics, University of Turin, Via Pietro
Giuria 1, I-10125, Turin, Italy  INAF – Osservatorio Astrofisico di Torino,
Via Osservatorio 20, I-10025, Pino Torinese, Italy  Prune August  Department
of Space Research and Technology, Technical University of Denmark, Kongens
Lyngby, Denmark  Luis Welbanks  School of Earth and Space Exploration, Arizona
State University, Tempe, AZ 85287  NHFP Sagan Fellow  Jean-Michel Desert
Anton Pannekoek Institute of Astronomy, University of Amsterdam, Amsterdam,
Netherlands  Jonathan Lunine  Department of Astronomy, Cornell University,
Ithaca, NY 14850  Jorge Sanchez  School of Earth and Space Exploration,
Arizona State University, Tempe, AZ 85287  Megan Mansfield  Steward
Observatory, University of Arizona, Tucson, AZ 85719  NHFP Sagan Fellow
Lorenzo Pino  INAF – Osservatorio Astrofisico di Arcetri, Largo Enrico Fermi
5, 50125 Firenze, Italy  Emily Rauscher  Department of Astronomy, University
of Michigan, Ann Arbor, MI 48109  Eliza Kempton  Department of Astronomy,
University of Maryland, College Park, MD 20742  Joseph Zalesky  Department of
Astronomy, University of Texas at Austin, Austin, TX 78712  Martin Fowler

###### Abstract

Ground-based, high-resolution and space-based, low-resolution spectroscopy are
the two main avenues through which transiting exoplanet atmospheres are
studied. Both methods provide unique strengths and shortcomings, and combining
the two can be a powerful probe into an exoplanet’s atmosphere. Within a joint
atmospheric retrieval framework, we combined JWST NIRSpec/G395H secondary
eclipse spectra and Gemini South/IGRINS pre- and post-eclipse thermal eclipse
observations of the hot Jupiter WASP-77A b. Our inferences from the IGRINS and
NIRSpec data sets are consistent with each other, and combining the two allows
us to measure the gas abundances of H2O and CO as well as the vertical thermal
structure with higher precision than either data set provided individually. We
confirm WASP-77A b’s subsolar metallicty ([(C+O)/H]=-0.61${}^{+0.10}_{-0.09}$)
and solar C/O ratio (C/O = 0.57${}^{+0.06}_{-0.06}$). The two types of data
are complementary, and our abundance inferences are mostly driven by the
IGRINS data while inference of the thermal structure is driven by the NIRSpec
data. Our ability to draw inferences from the post-eclipse IGRINS data is
highly sensitive to the number of singular values removed in the detrending
process, potentially due to high and variable humidity. We also search for
signatures for atmospheric dynamics in the IGRINS data and find that
propagated ephemeris error can manifest as both an orbital eccentricity or a
strong equatorial jet. Neither are detected when using more up-to-date
ephemerides. However, we find moderate evidence of thermal inhomogeneity and
measure a cooler nightside that presents itself in the later phases after
secondary eclipse.

planets and satellites: atmospheres

††facilities: Gemini South (IGRINS), JWST (NIRSpec)††software: maptlotlib
(Hunter, 2007), numpy (van der Walt et al., 2011), pymultinest (Buchner,
2016), scipy (Virtanen et al., 2019)

##  1 Introduction

A wealth of information is contained in the atmospheres of exoplanets.
Measuring their compositions and thermal structures at a population level
provides insight into a planet’s formation pathways (Öberg et al., 2011;
Madhusudhan et al., 2017) and atmospheric chemical and physical processes
(Fortney et al., 2008; Parmentier et al., 2018). One of the major goals of
exoplanet science is to synthesize this information and uncover population-
level trends in compositional diagnostics such as metallicity and the carbon-
to-oxygen ratio (e.g., Sing et al., 2016; Welbanks et al., 2019; Baxter et
al., 2020; Mansfield et al., 2021). Spectroscopic campaigns with the Hubble
and Spitzer Space Telescopes (hereafter HST and Spitzer) aimed to measure the
molecules H2O and CO, which make up a significant fraction of metals in hot
Jupiters, as well as CO2, whose strong absorption can provide another
indicator of atmospheric metallicity. Applying Bayesian atmospheric parameter
estimation tools (“retrievals”) to these datasets enabled constraints on the
water abundances in several planets (e.g., Kreidberg et al., 2014; Haynes et
al., 2015; Welbanks et al., 2019). However, due to limited wavelength
coverage, similar measurements of CO and CO2 remained elusive, and estimates
of metallicity and C/O were based largely on these H2O constraints.

An alternative to space-based spectroscopy is High Resolution (R =
$\lambda/\Delta\lambda>$20000) Cross-Correlation Spectroscopy (HRCCS) using
ground-based telescopes. HRCCS leverages the planetary signal’s Doppler shift
due to its orbital motion to disentangle the planetary atmospheric emission or
transmission lines from telluric and stellar lines. HRCCS has been used to
detect several C- and O-bearing species (e.g., de Kok et al., 2013; Birkby et
al., 2013; Hawker et al., 2018) in addition to numerous refractory metals
(e.g., Hoeijmakers et al., 2019; Gandhi et al., 2023; Pelletier et al., 2023)
in hot- and ultra-hot Jupiters. HRCCS also provides probes of atmospheric
dynamics and has been used to measure wind and rotation speeds (Snellen et
al., 2010; Brogi et al., 2016; Gandhi et al., 2022). Brogi & Line (2019) and
Gibson et al. (2020) demonstrated that absolute abundance and temperature
profile constraints could be retrieved from these datasets, placing HRCCS
observations on a similar footing as space-based transit spectoscopy for
quantitative estimation. Using simulated HST/WFC3 and CRIRES data, Brogi &
Line (2019) also showed that that high- and low-resolution data can be
combined in a joint retrieval process. The two types of data are
complementary: low-resolution data contains continuum information and probes
high pressures, and high-resolution data is sensitive to molecular line shapes
and probes lower pressures where the strongest line cores are located. Brogi &
Line (2019) predicted that such combinations would enable more precise
estimates of atmospheric properties than either data set could provide alone.
Using real data this has indeed been the case in the literature, but there are
only a handful of studies that have done such analysis (Gandhi et al., 2019;
Kasper et al., 2023; Boucher et al., 2023).

Line et al. (2021) recently demonstrated the HRCCS capabilities of IGRINS
(Immersion GRating INfrared Spectrometer, R$\sim$45,000, 1.45-2.6 micron) on
Gemini-South by applying the HRCCS retrieval framework to pre-secondary-
eclipse (0.325 $<\varphi<$ 0.47, $\varphi=$0.5 being secondary eclipse) data
of the hot Jupiter WASP-77A b
($T_{\mathrm{eq}}=1700\mathrm{K},R_{P}=1.21R_{J},P=1.36$ day; Maxted et al.,
2013). Owing to the instrument’s stability and large wavelength coverage, Line
et al. (2021) placed precise constraints on the planet’s thermal structure and
the abundances of H2O and CO ($\pm$0.1-0.2 dex). From the latter, they
inferred a sub-solar metallicity ([(C+O)/H] = -0.48${}^{+0.15}_{-0.13}$) and
near-solar C/O (0.59$\pm$0.08), indicative of more diverse formation pathways
commonly predicted in the literature (such as Madhusudhan et al. (2014);
Mordasini et al. (2016); Khorshid et al. (2021)).

Mansfield et al. (2022) investigated the atmosphere of WASP-77A b at low
resolution using HST/WFC3 (1.1-1.7 micron, R $\sim$ 70) and the Spitzer/IRAC
channels centered at 3.6 and 4.5 microns. Using a similar free retrieval
framework to Line et al. (2021), they were unable to obtain any informative
composition constraints, only placing a subsolar lower limit on the
metallicity broadly consistent with the IGRINS results. Again, WASP-77A b was
recently observed in eclipse using NIRSpec/G395H on the James Webb Space
Telescope (JWST, August et al., 2023). Bounded constraints could be placed on
the abundances of H2O and CO and an upper limit was placed on the abundance of
CO2 when applying a free retrieval, showing a clear improvement over WFC3 and
allowing the confirmation of the low metallicity measured with IGRINS by Line
et al. (2021).

While previous combined high- and low-resolution retrievals have only used HST
and Spitzer data (Gandhi et al., 2019; Kasper et al., 2023; Boucher et al.,
2023), no such analysis has been attempted with data taken with JWST. Both
IGRINS and NIRSpec have been shown to provide stringent estimates on
atmospheric composition and thermal structure, and combining the unique
strengths of both instruments could provide powerful probes into transiting
giant atmospheres. However, such a combination may also be challenged in ways
the previous high- and low- resolution combinations were not due to the fact
that the two instruments probe different altitudes and hence potentially
different gas abundances. Thus, a test of its feasibility and effectiveness is
needed. In this paper, we will perform a combined ground-based and JWST
retrieval analysis using the IGRINS and NIRSpec data of WASP-77A b.

We first present two additional nights of IGRINS data covering the post-
eclipse phases of WASP-77A b’s orbit in Section 2. In the subsequent sections,
we search for additional information that may be gained by incorporating these
new data through molecular detection via cross-correlation in Section 3,
searching for signatures of atmospheric dynamics in Section 4, and atmospheric
retrieval in Section 5. The results of combining the IGRINS and NIRSpec data
are also presented in Section 5. We place WASP-77A b’s atmosphere into context
and discuss the synergies between NIRSpec and IGRINS in Section 6 before
concluding in Section 7.

##  2 Observations & Data Reduction

WASP-77A System Properties
---
Stellar Properties |
Spec. Type | G8Vb
$R_{\star}$ | 0.955 $R_{\odot}^{a}$
$M_{\star}$ |  $1.002\pm 0.045$ $M_{\odot}^{b}$
$T_{\mathrm{eff}}$ | 5605 $K^{b}$
log$g$ |  $4.33\pm$ 0.08b
$M_{K}$ | $8.405\pm\ 0.031^{b}$
$\gamma$ |  $1.6845\pm$ 0.0004 km s${}^{{-1}^{b}}$
$\mathrm{[Fe/H]}$ | $0.00\pm\ 0.11^{b}$
$\mathrm{[C/H]}$ | $-0.02\pm 0.05^{d}$
| $0.10\pm 0.09^{e}$
| $-0.04\pm 0.06^{f}$
$\mathrm{[O/H]}$ | $0.06\pm 0.07^{d}$
| $0.23\pm 0.02^{e}$
| $-0.14\pm 0.06^{f}$
$\mathrm{C/O}$ | $0.46\pm 0.09^{d}$
| ${0.44^{+0.07}_{-0.08}}^{e}$
| $0.59\pm 0.08^{f}$
Planet Properties |
$R_{P}$ |  $1.21\pm\ 0.02$ $R_{J}^{b}$
$M_{P}$ |  $1.76\pm\ 0.06$ $M_{J}^{b}$
$T_{eq}$ | 1740 Kb
$K_{P}$ | 192 $\pm$ 4.5 km s-1
$e$ | ${0.0074^{+0.007}_{-0.005}}^{c}$
$\omega_{\star}$ | ${-166^{+66}_{-75}}^{\circ,{c}}$
$T_{C}$ |  $2455870.44977\pm 0.00014$ $\mathrm{BJD}^{b}$
|  $2457420.88439^{+0.00080}_{-0.00085}$ $\ \mathrm{BJD}^{c}$
$P$ |  $1.3600309\pm\ 0.0000020$ $\mathrm{day}^{b}$
|  $1.36002854\pm 0.00000062$ $\mathrm{day}^{c}$
$a$ |  $0.02405\pm\ 0.00036$ $\mathrm{AU}^{b}$
Table 1: Relevant stellar and planet parameters for WASP-77A b and its star.
$K_{P}$ is calculated by assuming a circular orbit with the reported semi-
major axis and period: $K_{P}=2\pi a/P$. References: a: Bonomo et al. (2017);
b Maxted et al. (2013); c: Cortés-Zuleta et al. (2020); d: Polanski et al.
(2022); e: Reggiani et al. (2022); f: Kolecki & Wang (2022).

The high-resolution data were taken over the course of three half nights in
December 2020 as part of observing program GS-2020B-Q-249 (P.I. J. Zalesky)
using the IGRINS instrument on Gemini South (Park et al., 2014; Mace et al.,
2018). The data taken on December 14, 2020, covering the pre-eclipse phases,
were previously presented and analyzed in Line et al. (2021). Data taken on
December 6 and December 21, each covering post-eclipse phases, are presented
here for the first time111The reduced spectral matrices and observing
conditions for all three nights as well as the spectral templates used for the
subsequent cross-correlation analysis are publicly available in a Zenodo
repository linked here.. For each night, we took a continuous sequence of 70s
exposures using an AB-BA nodding pattern (140s per AB pair, hereafter referred
to as “frames”), consistently achieving high ($\sim$200) SNR per pixel per
frame in most orders (Table 2). The primary star has a companion, WASP-77 B,
3.3” away ($\sim$ 10 slit widths) with a position angle 150∘ east of north.
The slit was rotated from the default 90∘ to 60∘ to avoid contamination when
nodding. Observing conditions on each night are summarized in Figure 1.

Date | Orbital Phase | Start (BJD) | End (BJD) | # Frames | Med. SNR H | Med. SNR K | Air Mass
---|---|---|---|---|---|---|---
12/06/2020 | 0.535-0.605 | 2459189.61326 | 2459189.75917 | 39 | 195 | 185 | 1.09 $\rightarrow$1.86
12/14/2020 | 0.325-0.47 | 2459197.53017 | 2459197.74095 | 79 | 205 | 190 | 1.12 $\rightarrow$ 1.09
|  |  |  |  |  |  |  $\rightarrow$ 1.82
12/21/2020 | 0.535-0.628 | 2459204.61064 | 2459204.74449 | 48 | 175 | 165 | 1.15 $\rightarrow$ 2.66
Table 2: Details for each of the three observing sequences taken with IGRINS.
Figure 1: Observing conditions on each of the three nights data were taken.
Figure 2: Observed phases. Transit is at $\phi$ = 0 and secondary eclipse is
at $\phi$=0.5. The colors of the sectors indicate whether the planet signal is
blue- or red-shifted during those phases. Not to scale.

The raw 1D spectra were extracted by the IGRINS facility team using the IGRINS
Pipeline Package (PLP; Lee & Gullikson, 2016; Mace et al., 2018). While the
PLP provides an initial wavelength solution, this solution can drift by up to
0.2 pixels (0.46 km s-1) over the course of an observing sequence (e.g., see
Brogi et al., 2023). To correct these small misalignments, for each sequence,
we re-fit each frame to the last frame in the sequence to place the entire
sequence on a common wavelength grid. We choose to align to the last frame in
the sequence as this is the closest in time to when the wavelength calibrator
is used by the PLP. We also normalize the counts in each frame, discard 8
orders due to low throughput and/or strong telluric contamination, and discard
200 pixels on the low throughput edges of each order. After this cleaning
process, we place the data into data cubes of shape $N_{\mathrm{orders}}\times
N_{\mathrm{frames}}\times N_{\mathrm{pixels}}$.

To detrend each raw spectral sequence, we apply a singular value decomposition
(SVD) to each order from a data cube using numpy.lingalg.svd (de Kok et al.,
2013; Brogi & Line, 2019; Line et al., 2021; Brogi et al., 2023). By visual
inspection, there are no remaining telluric artifacts after the first 4
principal components are removed. Neither the cross-correlation nor the
retrieval analyses change whether we choose 4, 6 or 8 principal components,
and for our initial analysis we remove 4 for all three sequences in order to
be consistent with the Line et al. (2021) analysis of the pre-eclipse data.
However, we find that removing only 3 is best for the post-eclipse nights (see
Section 6.5). For a given sequence, we save two matrices – the data cube
recomposed with the first 4 principal components removed and a scaling matrix
of the data cube recomposed using only the first 4 – as well as an array of
the times of each frame in BJD. In the subsequent sections, we will refer to
the sequence taken on 12/14/2020 as the pre-eclipse data. The sequences taken
on 12/06/2020 and 12/21/2020 will be referred to together as the post-eclipse
data, and in all analysis except the CCF trail, we sum their cross-correlation
coefficients and log-likelihoods together and treat this sum as if from a
single sequence.

##  3 Molecular Detection via Cross-Correlation

Figure 3: Cross-correlation maps for the pre- and post-eclipse sequences as
well as from all 3 nights of IGRINS data combined. Each column is the
resultant map from cross-correlating using a different atmospheric template.
Below each map is a cross section of both the CCF SNR (red) and
$\Delta$log$\mathcal{L}$ (divided by its standard deviation, blue) at the
highest likelihood $K_{P}$.

As an initial check of the quality of the data and to determine the strength
of the planetary signal, we cross-correlate a model spectrum with the post-SVD
data. We used a solar composition ([X/H] = 0; C/O = 0.55) 1D radiative-
convective-thermoequilibrium model (1D-RCTE) using the ScCHIMERA modeling
framework as described in Arcangeli et al. (2018); Piskorz et al. (2018);
Mansfield et al. (2022). The ScCHIMERA model provides the dayside average 1D
pressure-temperature (P-T) profile as well as gas abundance profiles. For a
hot Jupiter, C and O, in the form of H2O and CO, will contain the majority of
metals in the atmosphere. Therefore, H2O and CO, both of which have many lines
in the H and K bands, are of particular interest to detect.

A high-resolution (R=250,000) thermal emission model spectrum is computed by
passing the 1D-RCTE atmospheric structure through a GPU-accelerated version of
the atmospheric forward modeling code CHIMERA (Line et al., 2013; Brogi &
Line, 2019). The model spectrum is then convolved with the appropriate
equatorial rotation kernel (vsin$i$=4.2 km s-1) and a Gaussian intrumental
profile at the IGRINS nominal resolution. Template spectra that include only
H2O and CO individually (in addition to continuum opacities) are also computed
using the same P-T and abundance profiles as output by the solar composition
1D-RCTE model but with all other gas abundances set to zero. We include the
most recent ExoMol (Tennyson et al., 2020), HITEMP (Rothman et al., 2010), and
HITRAN opacities222H2-H2/He CIA cross-sections from Karman et al. (2019); H2O
line list from Polyansky et al. (2018) and absorption cross-sections computed
via the process described in Gharib-Nezhad et al. (2021); the CO isotopologue
cross-sections are from Li et al. (2015); CH4 from Hargreaves et al. (2020),
H2S from Azzam et al. (2016); NH3 from Coles et al. (2019); and HCN from
Barber et al. (2014). from H2-H2 and H2-He CIA, H2O, 12CO, 13CO, CH4, H2S,
NH3, and HCN.

We convert from $F_{P}$ to $F_{P}/F_{\star}$ by dividing by a smoothed
(Gaussian over 200 pixels) PHOENIX stellar spectrum (Husser et al., 2013) at
the appropriate $T_{\mathrm{eff}}$ and log$g$. Because the SVD process can
modify and stretch the planet signal, before each cross-correlation or
likelihood evaluation we follow the model injection procedure outlined in
Brogi & Line (2019), Line et al. (2021), and Brogi et al. (2023) to similarly
modify the model spectrum.

The planetary lines Doppler shift at each orbital phase as the planet orbits
the star, and the total line-of-sight velocity at a given time $t$ is:

$V_{\mathrm{LOS}}(t)=\gamma+V_{\mathrm{bary}}(t)+V_{P}(t)+V_{\mathrm{sys}}$
(1)

where $\gamma$ is the star-planet system’s radial velocity,
$V_{\mathrm{bary}}(t)$ is the Solar System barycentric radial velocity in the
observatory’s frame (via radial_velocity_correction from astropy), $V_{P}(t)$
is the planet’s velocity in the star-planet system’s barycentric frame, and
$V_{\mathrm{sys}}$ is an additive term to account for any systematic offset.
Assuming a circular orbit (however, see Section 4) $V_{P}(t)$ becomes
$K_{P}\sin[2\pi\times\varphi(t)]$, where $K_{P}$ is the planet’s radial
velocity semi-amplitude and $\varphi$ is its orbital phase. As is typical in
the literature (e.g., Birkby, 2018; Brogi & Line, 2019; Line et al., 2021), we
Doppler shift and cross-correlate the model spectrum on a 2D grid of possible
$K_{P}$ and $V_{\mathrm{sys}}$ values, creating a 2D map of correlation
coefficients. We then median subtract this map, find the standard deviation of
a 3-sigma-clipped copy, and divide the map by this standard deviation to
determine the detection signal-to-noise ratio of planetary absorption or
emission lines (Kasper et al., 2021, 2023). In addition to calculating
correlation coefficients, we use the likelihood function from Brogi & Line
(2019) to calculate a log-likelihood based detection map. Figure 3 shows the
cross-correlation (CCF) detection maps for the atmospheric signal from the
above mentioned model.

Similar to Line et al. (2021), we achieve a strong detection of absorption
lines from the planet’s thermal emission with a signal-to-noise-ratio (SNR) of
9.4 when cross-correlating the full atmospheric model with the pre-eclipse
data alone (Figure 3, top left). We note that this is different from Line et
al. (2021), who report a SNR of 12.3 due to their different method of
normalizing the CCF map. This and the many other different normalizing methods
in the literature demonstrates that the CCF SNR alone is not necessarily a
robust metric for planet signal strength. Using more statistically robust
methods, such as a Welch t-test (Brogi et al., 2012) or the log-likelihood
approach, may be more appropriate. Nonetheless, to keep in line with the
literature, we list here the CCF SNR values, but for each model-data
combination we also show the detection significances via the log-likelihood
maps in Figure 3, all of which are comparable to the CCF SNR.

Searching for the individual gases of interest, we detect H2O with a SNR of
9.2 and CO with a SNR of 3.6. We also use a spectral template without H2O and
CO to assess the presence of additional species still in the model like HCN,
CH4, and NH3. We are unable to detect these species in the CCF alone (Fig. 3,
top right), which is unsurprising because they have significantly lower
abundances in the model compared to H2O and CO.

Absorption features are also detected in both post-eclipse data sets. Summing
over both post-eclipse nights, the full-atmospheric SNR is 8 while H2O and CO
are detected with SNRs of 7.7 and 3.4, respectively (Figure 3, middle row).
Combining all three nights increases the full atmosphere SNR to 10.9. H2O, and
CO have SNRs of 10.7 and 5.0, respectively (Figure 3, bottom row). We note
that at this relatively high SNR, many of the structures in the CCF maps are
aliases of the planet signal or wings of the central CCF peak. Therefore, it
is difficult to truly estimate the CCF noise, and the actual detection
significance is likely underestimated. When cross-correlating the spectral
template with H2O and CO removed for both the post-eclipse data and all nights
combined, the remaining gases (CH4, NH3, HCN, 13CO, and H2S) are not detected
to any significance.

We also note that the correlated streaks in the CCF maps in Figure 3 lean in
opposite directions between pre- and post-eclipse phases, indicative that we
are indeed seeing the planet signal as it moves away and then back toward the
line-of-sight. We also search for the planetary signal’s line-of-sight Doppler
trail with orbital phase (Figure 4). If the planetary signal is truly present,
the CCF should trace out the predicted radial velocity with orbital phase
given the literature reported $K_{P}$ and $V_{\mathrm{sys}}$. Figure 4 shows
this is indeed the case. The fact that all correlation coefficients along the
planet’s path are positive also indicates that we are detecting absorption
features, not emission features, as would be expected for a hot Jupiter with
no thermal inversion in the infrared photosphere.

Figure 4: Cross correlation trails tracking the planet signal throughout the
phases captured in this program. The same 1D equilibrium model used to create
the 2D CCF maps in Figure 3 was cross correlated with each individual frame at
a range of line-of-sight velocities (horizontal axis). For visual clarity, the
two post-eclipse sequences have been combined and pairs of frames have been
binned together to boost the visible CCF. The dashed lines are the expected
planet radial velocity per the best-fit values of KP and Vsys from Section 4,
offset by 35 km s-1 on either side of the actual CCF trail for clarity. The
horizontal dashed lines indicate secondary eclipse, when the planet is
occulted by the star. The white space indicates phases at which no data were
taken.

##  4 Searching for Signatures of Atmospheric Dynamics

Figure 5: Posterior distributions for the radial velocity semi-amplitude,
$K_{P}$, and systemic offset velocity, $V_{\mathrm{sys}}$, with different
assumptions about WASP-77A b’s orbit as described in Section 4. Each label
states the period reference and eccentricity. T16 stands for Turner et al.
(2016) and CZ20 stands for Cortés-Zuleta et al. (2020).

The sub-km s-1 velocity precision from current high resolution spectroscopy
can enable sensitivity to the main atmospheric dynamical features in hot- or
ultra-hot Jupiters such as day-night winds or equatorial jets (Flowers et al.,
2019; Beltz et al., 2021). Phase-resolved CCFs can provide a powerful tool to
probe these dynamics by tracking the planet signal’s velocity in time
(Ehrenreich et al., 2020; Gandhi et al., 2022; Pino et al., 2022). In order to
test for velocity consistency between the pre- and post-eclipse sequences and
search for possible signatures of dynamics, we use the log-likelihood
formalism with the PyMultinest sampler (Feroz et al., 2009; Buchner, 2016) to
obtain quantitative constraints on $K_{P}$ and Vsys for a given atmospheric
template spectrum. We also include a multiplicative scale factor, $a$, that
can stretch the model spectrum to account for any line amplitude mis-matches
either due to a model inaccuracy or the SVD process. If multiple nights are
considered at once, we include a separate $a$ for each individual night. We
fit for $K_{P}$, $V_{\mathrm{sys}}$, and log${}_{10}(a)$ using the full solar
composition, H2O only, and CO only models from Section 3 to search for any
velocity offsets between gases as were found in Brogi et al. (2023) or between
pre- and post-eclipse sequences as were found in Pino et al. (2022).

The planet’s time-resolved velocity depends on its orbital phase, typically
calculated by dividing the time elapsed since a measured transit midpoint by
the period: $\varphi(t)=(t-T_{C})/P$. Using up-to-date ephemerides is crucial
– initially we used the midtransit time and period from WASP-77A b’s entry in
the Transiting Exoplanet Catalogue (TEPCat, Southworth, 2011, midtransit time
from Bonomo et al., 2017 and period from Turner et al., 2016) and found the
pre- and post-eclipse $K_{P}$’s to be inconsistent by 3$\sigma$ (Figure 5, top
left). To assess goodness-of-fit provided by each subset of the IGRINS data,
we fit Gaussian profiles to each frame of the CCF orbit trail to measure the
velocity at each frame time. To build signal, we took the average of bins of 5
frames weighted by the individual CCF amplitudes and treated the standard
error of that average as the 1$\sigma$ uncertainty on the velocity in that
time bin. As can be seen in Figure 6, the best fit $K_{P}$ and
$V_{\mathrm{sys}}$ from using both sequences provides a somewhat adequate fit
($\chi^{2}/N$=1.27, N=32), but the individual sequences fail to predict each
other (pre-eclipse and post-eclipse giving $\chi^{2}/N$’s of 4.25 and 2.97,
respectively).

Such velocity asymmetries can arise from an orbital eccentricity or
atmospheric dynamics, so we expand our velocity model to include either an
eccentricity or an equatorial jet (the details of both are described in
Appendix A). While the best fit eccentric model achieves consistency between
pre- and post-eclipse and provides a better fit ($\chi^{2}/N$=1.02 using both
sequences, Figure 6, second panel), we can only place an upper limit on the
eccentricity and there is insufficient evidence to favor a freely eccentric
orbit over a circular one. The jet model provides an even better fit
($\chi^{2}/N$=1.01, Figure 6, third panel) but “detects” a supersonic westward
jet at 3$\sigma$ (by comparing Bayes factors), which is physically
implausible.

A more recent ephemeris analysis of WASP-77A b lists the eccentricity as
$e=0.0074^{+0.007}_{-0.005}$. Fixing the eccentricity to this value once again
achieves consistency between pre- and post-eclipse phases and is favored over
the circular orbit by 4.2$\sigma$ (Figure 5, lower left). However, when using
the more up-to-date midtransit time from Cortés-Zuleta et al. (2020), the
$K_{P}$ values become inconsistent again and this eccentricity is disfavored
over a circular orbit, indicating 0.0074 is too large of an eccentricity.
Again, there is insufficient evidence for an eccentricity measured on our own,
and a circular orbit perfectly fits the data ($\chi^{2}$/N=1.00, Figure 6) and
yields consistent $K_{P}$ values between sequences (Figure 5, upper right).
Fitting for a jet speed with the updated ephemerides is also disfavored over a
simply circular orbit by 2.5$\sigma$.

This points to the velocity asymmetries arising from propagated ephemeris
error. Indeed, if we add a time correction parameter to our initial circular
orbit analysis with the TEPCat ephemerides, we are able to measure that
WASP-77A b appears 8.53${}^{+2.41}_{-2.39}$ minutes behind in phase, which
translates to an initial measured period error of 0.21$\pm$0.06 seconds.
Including this offset achieves consistent velocities between pre- and post-
eclipse and is favored over the initial circular orbit fit by 3.6$\sigma$
(Figure 5, bottom right). Fitting for such a correction with the Cortés-Zuleta
et al. (2020) midtransit time and period yields a value of 0 minutes to one
part in a thousand. Therefore, we can conclude that WASP-77A b’s orbit is
effectively circular, but propagated ephemeris error from the TEPCat values
induced an effective eccentricity. It should be noted that had the true
eccentricity been as “large” as 0.0074, we would have been sensitive to its
effects.

In the future, the impact of propagated ephemeris error can be mitigated by
using the most recent measured midtransit time before a given observation when
calculating orbital phase. This may seem obvious, but many planets only have a
few published ephemerides. In the likely scenario that published ephemerides
measured recently before a given HRCCS observation don’t exist, we recommend
taking advantage of the several campaigns monitoring known exoplanet host
stars including with the Transiting Exoplanet Survey Satellite (e.g., Wong et
al., 2020) and citizen science projects like ExoClock (Kokori et al., 2022)
and Exoplanet Watch (Zellem et al., 2020). Searching for entries in the
American Association of Variable Star Observers (AAVSO)
database333https://www.aavso.org/, we discovered that a transit of WASP-77A b
was observed with the 6-inch MicroObservatory and subsequently analyzed
through the Exoplanet Watch citizen science program in late October of 2020.
Estimating once again WASP-77A b’s orbital velocity using the Turner et al.
(2016) period with this midtransit time, the pre- and post-eclipse $K_{P}$’s
are consistent.

Besides line positions, line shapes are also affected by atmospheric dynamics.
Assuming a circular orbit, we added an average line full width at half maximum
(FWHM) parameter. This is the FWHM in pixels of a Gaussian kernel we convolved
the 1D-RCTE model spectrum with instead of the nominal instrument profile and
rotation kernels. For the pre- and post-eclipse sequences we were able to
measure this FWHM to be 9.34${}^{+0.59}_{-0.49}$ and 8.39${}^{+0.76}_{-0.68}$
pixels, respectively, which is consistent with what we would expect from a
combination of the instrument profile (4.5 pixels) and tidally locked
planetary rotation (4.52 km s-1, $\sim$7 pixels). There is nothing to suggest
that dynamics or anything else beyond these two sources are affecting the line
shape.

Finally, we find no significant velocity offsets between different gases.
Therefore, we detect no signatures of atmospheric dynamics that may bias a 1D
atmospheric retrieval analysis. Such analysis is presented in the following
section.

Figure 6: Measured velocity from the CCF signal in the planet literature rest
frame (diamonds) and best fits from various assumptions about WASP-77A b’s
orbit. The purple lines are the estimated velocity trail using both IGRINS
sequences, and the red and blue lines are estimates from each sequence alone.
Pre-eclipse is red and post-eclipse is blue. The reduced $\chi^{2}$ from each
best fit is also listed in its corresponding color. Only a combined pre+post
fit appears for the jet model because after that initial analysis it was clear
the jet was not a viable path of inquiry. Note that some of the data points
appear to be at the same phase because the CCF frames were binned in BJD time
and have been phase folded for visualization.

##  5 Estimating Gas Abundances and Vertical Thermal Structure

To estimate WASP-77A b’s atmospheric composition and vertical thermal
structure, we apply the same retrieval framework described in Line et al.
(2021) and Brogi et al. (2023). Using PyMultinest paired with a GPU-
accelerated version of CHIMERA (Line et al., 2013, 2021), we estimate the
constant-with-altitude abundances of the same gases mentioned in Section 3,
the 6-parameter analytic pressure-temperature (P-T) profile described in
Madhusudhan & Seager (2009), and the velocities and scale factor estimated in
Section 4. As in Line et al. (2021), we indirectly estimate the 13CO/12CO
isotopologue ratio by estimating directly [13CO/12CO] = log(13CO/12CO)Planet
\- log(13CO/12CO)Earth, where (13CO/12CO)Earth = 1/89 (Meibom et al., 2007).
For example, a retrieved [13CO/12CO] of 0.5 corresponds to a 13CO/12CO ratio
of 1/(89-0.5) = 1/28. Table 3 lists each parameter and their uniform prior
bounds as well as derived quantities such as the C/O ratio.

For retrievals on the NIRSpec/G395H data, the model spectra are computed at
R=100,000 over those wavelengths ($\sim$2-5 micron). We then assume a top hat
profile to bin the spectra down to the same wavelength bin widths as the data,
which has an average R$\sim$250\. We perform retrievals with the following
data combinations: IGRINS pre-eclipse only, IGRINS post-eclipse nights
combined, all IGRINS data combined, NIRSpec eclipse only, and IGRINS pre-
eclipse and NIRSpec eclipse combined. Table 3 lists each parameter, its prior,
and its posterior median value for each retrieval. When considering both
IGRINS and NIRSpec data simultaneously, we add the log-likelihood of the low-
resolution data ($-\frac{1}{2}\chi^{2}$) to the IGRINS log-likelihood
following Brogi & Line (2019). For the retrievals using only IGRINS data we
use 500 live points, and we use 2500 for any retrieval including the NIRSpec
data. This is because in our initial, exploratory retrieval analyses, the
NIRSpec retrieval results varied with model resolution and number of live
points. This variance asymptotically decreased as we increased both to the
values adopted here. This behavior was not observed for the IGRINS-only
retrievals. In the following subsections, we summarize the dependence of gas
abundance and P-T profile estimates on the inclusion of different nights of
IGRINS data (Figure 7) and the NIRSpec data (Figure 8).

Parameter | Prior | Pre-Eclipse | Post-Eclipse | Pre+Post | NIRSpec | Pre+NIRSpec
---|---|---|---|---|---|---
$\log_{10}(n_{\mathrm{H_{2}O}})$ |  $\mathcal{U}$(-12, 0) | -3.97${}^{+0.11}_{-0.10}$ | >-3.89 | -3.80${}^{+0.16}_{-0.12}$ | -3.80${}^{+0.34}_{-0.28}$ | -4.02${}^{+0.07}_{-0.06}$
$\log_{10}(n_{\mathrm{CO}})$ |  $\mathcal{U}$(-12, 0) | -3.81${}^{+0.19}_{-0.17}$ | >-4.28 | -3.60${}^{+0.22}_{-0.20}$ | -3.73${}^{+0.55}_{-0.52}$ | -3.91${}^{+0.13}_{-0.13}$
$\log_{10}(n_{\mathrm{CH_{4}}})$ |  $\mathcal{U}$(-12, 0) | <-5.97 | <-1.77 | <-5.57 | - | <-6.93
$\log_{10}(n_{\mathrm{H_{2}S}})$ |  $\mathcal{U}$(-12, 0) | <-4.45 | <-1.50 | <-4.36 | - | <-5.07
$\log_{10}(n_{\mathrm{NH_{3}}})$ |  $\mathcal{U}$(-12, 0) | <-5.88 | <-3.72 | <-6.01 | - | <-6.03
$\log_{10}(n_{\mathrm{HCN}})$ |  $\mathcal{U}$(-12, 0) | <-5.69 | <-2.43 | <-5.46 | - | <-.612
$[^{13}\mathrm{CO}/^{12}\mathrm{CO}]_{\oplus}$ |  $\mathcal{U}$(-5, 5) | 0.63${}^{+0.30}_{-0.34}$ | <0.70 | 0.26${}^{+0.33}_{-0.43}$ | <1.78 | 0.33${}^{+0.35}_{-0.71}$
$\log_{10}(n_{\mathrm{CO_{2}}})$ |  $\mathcal{U}$(-12, 0) | - | - | - | <-6.34 | <-7.25
$T_{0}$ [K] |  $\mathcal{U}$(500,2500) | 1470${}^{+210}_{-370}$ | 1480${}^{+110}_{-110}$ | 1460${}^{+180}_{-280}$ | 1300${}^{+80}_{-50}$ | 1400${}^{+30}_{-40}$
$\log_{10}P_{1}$ [$\log_{10}$ bar] |  $\mathcal{U}$(-5.5,2.5) | >-4.25 | -1.78${}^{+0.81}_{-0.61}$ | >-2.14 | >-1.37 | 0.69${}^{+0.22}_{-0.24}$
$\log_{10}P_{2}$ [$\log_{10}$ bar] |  $\mathcal{U}$(-5.5,2.5) | Unc. | Unc. | Unc. | Unc. | Unc.
$\log_{10}P_{3}$ [$\log_{10}$ bar] |  $\mathcal{U}$(-2,2) | >-0.83 | Unc. | -0.59${}^{+0.24}_{-0.18}$ | <-1.16 | >-0.51
$\alpha_{1}$ |  $\mathcal{U}$(0.02,2) | >0.31 | >0.29 | >0.29 | 0.48${}^{+0.04}_{-0.03}$ | 0.60${}^{+0.04}_{-0.05}$
$\alpha_{2}$ |  $\mathcal{U}$(0.02,2) | Unc. | Unc. | Unc. | Unc. | >0.03
$dK_{P}$ [km/s] |  $\mathcal{U}$(-20,20) | -1.26${}^{+0.85}_{-0.85}$ | -9.62${}^{+2.02}_{-2.07}$ | -5.75${}^{+0.27}_{-0.27}$ | - | -1.29${}^{+0.75}_{-0.78}$
$dV_{\mathrm{sys}}$ [km/s] |  $\mathcal{U}$(-20,20) | -5.27${}^{+0.51}_{-0.53}$ | -3.84${}^{+0.86}_{-0.86}$ | -2.46${}^{+0.15}_{-0.16}$ | - | -5.27${}^{+0.48}_{-0.47}$
$\log(a_{1})$ |  $\mathcal{U}$(-2,2) | - | 0.80${}^{+0.37}_{-0.44}$ | 0.04${}^{+0.53}_{-0.33}$ | - | -
$\log(a_{2})$ |  $\mathcal{U}$(-2,2) | 0.15${}^{+0.52}_{-0.53}$ | 0.73${}^{+0.38}_{-0.43}$ | 0.14${}^{+0.54}_{-0.32}$ | - | 0.03${}^{+0.04}_{-0.04}$
$\log(a_{3})$ |  $\mathcal{U}$(-2,2) | - | - | -0.05${}^{+0.55}_{-0.33}$ | - | -
$\log(a_{4})$ |  $\mathcal{U}$(-2,2) | - | - | - | 0.01${}^{+0.02}_{-0.03}$ | 0.03${}^{+0.01}_{-0.02}$
C/O |  | 0.59${}^{+0.07}_{-0..07}$ | 0.58${}^{+0.16}_{-0.19}$ | 0.61${}^{+0.08}_{-0.07}$ | 0.54${}^{+0.15}_{-0.16}$ | 0.57${}^{+0.06}_{-0.06}$
[(C+O)/H]⊙ |  | -0.53${}^{+0.16}_{-0.13}$ | >-0.68 | -0.32${}^{+0.20}_{-0.17}$ | -0.42${}^{+0.49}_{-0.42}$ | -0.61${}^{+0.11}_{-0.10}$
[O/H] |  | -0.54${}^{+0.15}_{-0.12}$ | >-0.62 | -0.34${}^{+0.19}_{-0.16}$ | -0.41${}^{+0.45}_{-0.38}$ | -0.61${}^{+0.10}_{-0.09}$
[C/H] |  | -0.51${}^{+0.15}_{-0.12}$ | >-0.97 | -0.29${}^{+0.22}_{-0.20}$ | -0.43${}^{+0.55}_{-0.52}$ | -0.60${}^{+0.13}_{-0.13}$
${}^{13}\mathrm{CO}/^{12}\mathrm{CO}$ |  | 1/21${}^{+25}_{-11}$ | <1/18 | 1/50${}^{+127}_{-27}$ | <1/1.48 | 1/31${}^{+92}_{-18}$
Table 3: Free parameters and their retrieved values for each retrieval. Values
listed with uncertainties are bounded constraints, while values with a $>$ or
$<$ are 3-$\sigma$ lower or upper limits, respectively. Entries with “Unc.”
are unconstrained, and those with “-” were not included as free parameters for
that specific retrieval. The parameters below the horizontal line are derived
and not retrieved for directly. The scale factors, $a_{i}$, correspond to each
specific night/data set. 1, 2, 3, and 4 are IGRINS nights 12/6, 12/14, 12/21,
and the low resolution data, in that order. Figure 7: Marginalized posterior
distributions of relevant gas abundances as well as the P-T profile from the
three different subsets of IGRINS data: pre-eclipse (red), post-eclipse
(blue), and all three nights (purple). The vertical P-T panel shows the median
posterior profile (solid lines) and the shaded regions are the 1$\sigma$
confidence intervals. Also plotted is the $\tau$=1 spectrum assuming the best
fit parameters from the pre+post retrieval.

###  5.1 Pre-Eclipse Data Only

Constraints from the pre-eclipse data alone, shown in red in Figure 7, are
identical within the statistical noise of the nested sampler to those
presented in Line et al. (2021). We place bounded constraints444We refer to
marginalized posterior distributions that are not against the prior bounds as
“bounded” and report the median value and 1$\sigma$ quantiles on either side
of this median. For posterior distributions against the upper (lower) prior,
we report 3$\sigma$ lower (upper) limits. Posteriors against both prior bounds
are “unconstrained”. on the the gas abundance of both H2O and CO to within 0.2
dex as well as a bounded constraint on the 13CO/12CO isotopologue ratio at
1/22${}^{+14}_{-10}$. This ratio is enriched compared to the terrestrial
standard (1/89; Meibom et al., 2007) and local ISM (1/70; Wilson & Rood,
1994). We can only place upper limits on the other gases in the model, all of
which are expected to be much less abundant than H2O and CO following
expectations from equilibrium chemistry. The median retrieved P-T profile is
non-inverted and temperature is monotonically increasing with pressure. These
gas abundances are physically plausible under assumptions of both equilibrium
and disequilibrium chemistry, and the P-T profile is consistent with an
atmosphere with efficient day-to-night heat transport and/or day-night cold
trapping (Line et al., 2021). All bounded values are within 1$\sigma$ of their
reported values in Line et al. (2021). For comparison to NIRSpec, we also
performed a retrieval on the pre-eclipse data with the abundance of CO2 as an
added parameter while leaving out the rest of the trace gases besides H2O and
CO. Unsurprisingly, only an upper limit is placed on the CO2 abundance.

###  5.2 Including the Post-Eclipse Data

The retrieval constraints from the post-eclipse data alone are summarized in
blue in Figure 7. Initially, when removing the first 4 principal components,
we struggled to make any informative inferences and could only place lower
limits on the abundances of H2O and CO. Adjusting to remove only the first 3
principal components makes a drastic improvement in our inference capabilities
and enable us to place bounded constraints on both H2O and CO consistent with
the pre-eclipse values albeit with slightly lower precision. We can only place
an upper limit on the 13CO/12CO ratio, but it is consistent with the retrieved
pre-eclipse value. The top-of-atmosphere temperature is well constrained, but
the rest of the retrieved P-T profile is poorly constrained and near-
isothermal. The slightly less stringent constraints compared to the pre-
eclipse retrieval may be due to the detrending process but also due to a
changing P-T profile with phase that the 1D model struggles to capture. Both
are discussed further in Section 6.

When we combine all three nights of IGRINS data, shown in purple in Figure 7,
the constraints on H2O and CO are consistent with both the pre- and post-
eclipse sequences individually. The median retrieved values are between those
from the pre- and post-eclipse sequences, weighted more in favor of the pre-
eclipse sequence, which appears to be driving these inferences. These
constraints are not more precise than those provided by the pre-eclipse data,
suggesting that the post-eclipse data contributed little new compositional
information.

###  5.3 NIRSpec

Figure 8: Similar to Figure 7, marginalized posterior distributions of
relevant gas abundances as well as the P-T profile from the pre-eclipse IGRINS
data alone (red), NIRSpec alone (blue), and both combined (purple). From the
confidence intervals of these parameters, it is evident that in the combined
retrieval, the composition inferences are driven by the high resolution data
while the P-T inferences are driven by the low resolution data.

We apply the same retrieval framework to the NIRSpec/G395H secondary eclipse
data presented in August et al. (2023). For initial exploratory retrievals
with lower model spectral resolution and fewer live points, we include all of
the gases used in the high-resolution retrievals in addition to CO2 and SO2,
prominent absorbers over the 4-5 micron region. Just as with the IGRINS
retrievals, we could only place bounded constraints on H2O and CO and upper
limits on all the other gases. The upper limits on the other gases did not get
more stringent than the IGRINS limits, and for the final, fiducial retrieval
with 2500 live points and R=100,000, we only include H2O, 12CO, 13CO, and CO2
in addition to the continuum opacities.

The results of this retrieval are summarized in blue in Figure 8. The
constraints both on the composition and vertical thermal structure are
consistent with IGRINS. The constraints on the H2O and CO abundances are
slightly less precise ($\sim 3\times$) than what we achieved with IGRINS. We
can only place an upper limit on the 13CO/12CO ratio, but it is consistent
with the measured IGRINS value. The P-T profile is more precisely constrained
in the upper atmosphere due to the lower pressures probed at these longer
wavelengths. Only an upper limit can be placed on the abundance of CO2, and
its inclusion as a model parameter is only favored by 1.8$\sigma$ (Bayes
factor 2). Within the context of theoretical predictions that the abundance of
CO2 is highly sensitive to atmospheric metallicity (e.g., Lodders & Fegley,
2002; Zahnle et al., 2009; Moses et al., 2013), the absence of a clear CO2
absorption feature around 4.5 microns further qualitatively confirms the sub-
solar metallicity of WASP-77A b’s atmosphere as previously found with the
IGRINS data (in contrast to published JWST results showing CO2, e.g., JWST
Transiting Exoplanet Community Early Release Science Team et al., 2023; Bean
et al., 2023).

###  5.4 IGRINS+NIRSpec

Figure 9: Top: Best fit model spectrum from the combined IGRINS+NIRSpec
retrieval at R=500,000 (gray), converted to brightness temperature and
smoothed to R=1000 (black) for visual clarity. The shaded regions illustrate
contributions and light absorbed from H2O, CO, and CO2 individually. The
NIRSpec data is shown as gold diamonds and has been further binned down to
R$\sim$85 for visual clarity. While the planet spectrum is not visible in the
actual IGRINS data, to give a sense of its wavelength coverage and spectral
resolution, the best fit has been further convolved to the IGRINS resolution
(R=45,000) and interpolated onto the instrument’s wavelength grid, alternating
colors with echelle order. Goodness of fit with both data sets is discussed in
detail in Section 6.3. Also shown for comparison are the previously published
HST/WFC3 eclipse measurements presented in Mansfield et al. (2022) but not
considered in this study (circles). Bottom Left: Cross-correlation residuals
with the IGRINS data after the best-fit spectrum has been divided out. A 2D
detection map was calculated as in Section 3, and shown here is a horizontal
slice at the planet’s $K_{P}$. The lack of a signal at rest velocity indicates
the best fit spectrum is an adequate fit to the true planet signal. This
process is described more in Section 6.3. Bottom Right: Residuals from
subtracting the binned NIRSpec data from the best fit binned onto the same
wavelength grid (black squares in top panel).

Here, we combine the IGRINS and NIRSpec data into a single retrieval. Because
the post-eclipse data did not provide new information, we only use the pre-
eclipse data in the interest of computational efficiency (high resolution
retrievals take days to weeks to complete). The results of this combined
retrieval are shown in purple in Figure 8. The inferences of the planet’s
atmosphere are consistent with both the IGRINS and NIRSpec individual
analyses, and the constraints on both the composition and vertical thermal
structure are more precise than either data set provided alone (composition by
about 30% compared to pre-eclipse only and P-T profile by about 50% compared
to NIRSpec alone). It is apparent that the composition constraints are driven
by the IGRINS data, whereas the NIRSpec observations provide more precise
contraints on the P-T profile. This is not surprising as the low resolution
data retains continuum information while the high resolution data is more
sensitive to gas abundances via individual molecular line shapes and ratios.
Similar to the NIRSpec-only retrieval, we can only place an upper limit on the
abundance of CO2, but a bounded constraint is placed on the 13CO/12CO ratio
consistent with the IGRINS-only retrieval. The best fit model spectrum is
shown in Figure 9, showing remarkable agreement with both data sets across
this wide wavelength range (1.5-5 micron).

##  6 Discussion

###  6.1 Bulk Composition and Placing WASP-77A b into Context

A major goal of exoplanet science is to tie back the composition of planetary
atmospheres to formation pathways via diagnostics like bulk metallicity and
the C/O ratio. These diagnostics can be derived from our inferred gas
abundances like so:

⊙ $\displaystyle=$ (2)
$\displaystyle\mathrm{log}_{10}\Bigg{[}\frac{n_{\mathrm{H2O}}+2n_{\mathrm{CO}}+n_{\mathrm{CH4}}+n_{\mathrm{HCN}}+3n_{\mathrm{CO_{2}}}}{2n_{\mathrm{H2}}[(n_{\mathrm{O}}+n_{\mathrm{C}})/n_{\mathrm{H}}]_{\odot}}\Bigg{]}$

and

$\mathrm{C/O}=\frac{n_{\mathrm{C}}}{n_{\mathrm{O}}}=\frac{n_{\mathrm{CO}}+n_{\mathrm{CH4}}+n_{\mathrm{HCN}}+n_{\mathrm{CO2}}}{n_{\mathrm{H2O}}+n_{\mathrm{CO}}+2n_{\mathrm{CO_{2}}}}$
(3)

where $n_{i}$ is the abundance of gas $i$. We assume the solar values from
Asplund et al. (2009). The median values of these derived quantities from each
retrieval set up are listed in Table 3. For the combined IGRINS+NIRSpec
analysis, these are [(C+O)/H] = -0.61${}^{+0.11}_{-0.10}$, C/O =
0.57${}^{+0.06}_{-0.06}$, [O/H] = -0.61${}^{+0.10}_{-0.09}$, and [C/H] =
-0.60${}^{+0.13}_{-0.13}$. These values are consistent with August et al.
(2023)’s values for [M/H] (-0.91${}^{+0.24}_{-0.16}$) and C/O
(0.36${}^{+0.10}_{-0.09}$) within $\sim$2$\sigma$. The slight differences can
be attributed to the many differences between the two analyses including
modeling assumptions, model resolution, number of live points, and the
modeling code itself. Slight differences in interpretations between modeling
codes is a known occurrence for other JWST data sets (Taylor et al., 2023,
e.g.,), and model synthesis in the context of JWST exoplanet retrievals is an
ongoing effort beyond the scope of this paper.

Line et al. (2021) interpreted WASP-77A b’s metallicity to be sub-stellar as
previous studies had measured WASP-77A’s metallicity (via [Fe/H]) to be
consistent with solar (0.00 $\pm$ 0.11, -0.10${}^{+0.10}_{-0.11}$ Maxted et
al., 2013; Cortés-Zuleta et al., 2020, respectively). More recent studies
measuring the stellar [O/H] and [C/H] have since placed WASP-77A’s [(C+O)/H]
at slightly super-solar (0.32 $\pm$0.04, 0.33$\pm$0.09; Polanski et al., 2022;
Reggiani et al., 2022), and the qualitative interpretation of the planet’s
sub-stellar metallicity does not change. Polanski et al. (2022) and Reggiani
et al. (2022) also both measured WASP-77A’s C/O ratio to be slightly sub-solar
(0.46 $\pm$0.09 and 0.44${}^{+0.07}_{-0.08}$, respectively). While our median
retrieved C/O ratio for the planet could be interpreted then as super-solar,
it is consistent with these new stellar values within about 2$\sigma$. On the
other hand, Kolecki & Wang (2022) measure WASP-77A’s C/O ratio to be $0.59\pm
0.08$, in which case the the planet’s C/O ratio would be almost exactly
stellar. Furthermore, if we account for partial sequestration of the total O
inventory due to rainout of refractory condensates (Burrows & Sharp,
1999)555$n_{O,true}=n_{O,observed}+3.28\times n_{Si}$. For $n_{Si}$ we assumed
the Si/O ratio is the same as the star using [O/H] and [Si/H] from Polanski et
al. (2022)., which is plausible as discussed in Section 6.6, the C/O ratio
drops to 0.52$\pm$0.06, which is consistent with the stellar values from
Polanski et al. (2022) and Reggiani et al. (2022).

While we present here extremely precise composition estimates, linking a
planet’s composition to its formation history is not trivial and depends
heavily on planet and disk modeling assumptions. As stated in Line et al.
(2021), the combination of a substellar metallicity and stellar C/O ratio is
not a common prediction from many planet formation theories (Öberg & Bergin,
2016; Madhusudhan et al., 2017; Khorshid et al., 2021, e.g.,), but it is not
implausible. One possible pathway under this scenario is formation via pebble
accretion interior to the H2O ice line from C-depleted gas with little-to-no
planetesimal pollution or core erosion (Madhusudhan et al., 2017). Formation
beyond the CO2 ice line is also possible if the planet migrated after disk-
dissipation regardless of if the planet’s final C/O is stellar (Madhusudhan et
al., 2014; Schneider & Bitsch, 2021) or super-stellar (Reggiani et al., 2022).
Schneider & Bitsch (2021) further predict that low metal content is more
common for planets that formed at large distances in the presence of pebble
evaporation. However, low metallicities combined with stellar C/O ratios are
rare in these models and depend on a low disk viscosity. Mousis et al. (2019)
construct a model in which devolatilization of amorphous ice from pebbles
enriches the gas in Jupiter’s feeding zone, which enhances the C/O ratio while
keeping the overall metallicity modest, but this model is tuned to the
specific problem of Jupiter’s composition.

These formation pathways are supported by more recent, targeted studies –
Bitsch et al. (2022) and Khorshid et al. (2023) run formation simulations
specifically for WASP-77A b, and both infer that the planet most likely formed
beyond the CO2 ice line and migrated late. The enrichment of 13CO in WASP-77A
b’s atmosphere may also be an indicator of formation at large distances. Zhang
et al. (2021) suggest a similar enrichment in the atmosphere of the planetary
mass companion YSES-1 b may be due to a combination of less efficient 12CO
self shielding and more efficient isotope exchange reactions at large orbital
separations in the protoplanetary disk. However, this is speculative and no
standard model linking isotope ratios with exoplanet formation has been
developed in the literature. Ultimately, these measurements of WASP-77A b’s
composition will be more meaningful in the future in the context of a larger
sample size of similar measurements with which we can test the wide range of
planet formation theories.

Figure 10: WASP-77A b in context of other hot- and ultra-hot Jupiters with
both measured [O/H] and [C/H] values. These include HD 209458 b (Gandhi et
al., 2019, VLT/CRIRES and HST/WFC3), $\tau$ Bootis b (Pelletier et al., 2021,
CFHT/SPIRou), KELT-20 b (Kasper et al., 2023, GN/MAROON-X and HST/WFC3), HD
149026 b (Bean et al., 2023, JWST/NIRCam), and WASP-96 b (Taylor et al., 2023,
JWST/NIRISS). Note the first two studies also combined high- and low-
resolution data to obtain these measurements. The values and uncertainties
here for WASP-77A b are from the combined IGRINS+NIRSpec retrieval (turquoise
diamond), which were mostly driven by the IGRINS data, but the NIRSpec-only
constraints are also shown for comparison (denoted with a hexagon). Also
included is Jupiter (Atreya et al., 2016) to illustrate the high precision
achieved on WASP-77A b’s composition with the combined IGRINS and NIRSpec
analysis, as well as regions in this parameter space associated with several
broad predicted formation pathways as summarized in Reggiani et al. (2022).

Such inferences about a planet’s formation history were more difficult in the
era of the Hubble Space Telescope. WFC3 was mainly sensitive to H2O and not CO
or CO2, so attempting to measure the complete O or C inventory in a hot
Jupiter was impossible. With the advent of high resolution retrievals (Brogi &
Line, 2019; Gibson et al., 2020), the abundances of both H2O and CO, and
indirectly [O/H] and [C/H], have been measured in a handful of hot- and ultra-
hot Jupiters in the past few years (Gandhi et al., 2019; Line et al., 2021;
Pelletier et al., 2021; Kasper et al., 2023). The launch of JWST has increased
our capabilities to do so and already the sample of such measurements has
grown (Bean et al., 2023; August et al., 2023; Taylor et al., 2023).
Nonetheless, the current dearth of planets in which these two have been
reliably estimated is stark (Figure 10), highlighting the long path ahead for
characterizing the transiting giant planet population as a whole.

In the landscape of other transiting giants for which [O/H] and [C/H] have
both been measured in some capacity (i.e., constrained), IGRINS has provided
some of the most precise measurements to date, on par with measurements of
[O/H] and [C/H] for planets in the solar system. However, as we demonstrated
here and as apparent in the case of WASP-18 b (Brogi et al., 2023; Coulombe et
al., 2023), JWST has a better grasp on vertical thermal structure.
Measurements of both composition and climate go hand-in-hand and are necessary
to interpret each other – thus there is high value in combining high- and low-
resolution data as the capabilities of both are expanded.

###  6.2 Comparison to Mansfield et al. 2022

WASP-77A b was previously observed in eclipse by Mansfield et al. (2022) using
HST/WFC3 and Spitzer’s IRAC channels centered at 3.6 and 4.5 microns. Using a
similar “free” chemistry prescription and the same analytic P-T profile
parameterization as used in this paper, only a rough lower limit could be
placed on the planet’s atmospheric metallicity. However, through a 1D-RCTE
model grid search, a moderately super-solar metallicity was inferred,
apparently in tension with the sub-solar metalicity preferred by the IGRINS
data. A frequentist $\chi^{2}$ statistic could not reliably place a clear
preference for either case, with the low metallicity IGRINS best-fit from Line
et al. (2021) giving a $\chi^{2}_{\nu}$ of 1.32 and the super-solar
metallicity best grid fit yielding 1.24. Qualitatively, the WFC3 data fall
well within the Line et al. (2021)’s posterior distribution of model spectra,
indicating that while the two data sets are consistent, the WFC3 data may be
insufficient, either in quality or wavelength coverage, to reliably place
precise constraints on gas abundances.

When incorporating the WFC3 data in their NIRSpec analysis, August et al.
(2023) inferred a higher, solar metallicity than from the NIRSpec data alone.
However, again frequentist metrics could not place a clear preference between
the solar and sub-solar metallicity best fit models, and those authors suggest
the WFC3 data may be unreliable. In initial exploratory retrieval analyses, we
tested combining the WFC3 data with the IGRINS data and found negligible
differences in the inferences or precision on the atmospheric composition.
Because the data did not appear to contribute new information, we chose not to
include the WFC3 and Spitzer data in the combined IGRINS and NIRSpec analyses.
In regards to pressures probed and estimating the vertical thermal structure,
NIRSpec has access to the same pressures as WFC3 (the role this plays in
inferences of the P-T profile is discussed further in Section 6.3), so little
vertical thermal information was lost by excluding the WFC3 data.

Our combined IGRINS+NIRSpec best fit spectrum is broadly consistent with the
WFC3 data (Figure 9). However, the spectral slope starts to diverge in the
reddest data points and the reduced chi-square statistic with the WFC3 data is
$\chi^{2}/N$ = 1.92. Investigating whether any specific points are driving
inferences toward a higher or lower metallicity, such as with a leave-one-out
cross validation analysis (Welbanks et al., 2023), may elucidate why both
IGRINS and NIRSpec are yielding different results than WFC3. However, such an
analysis is beyond the scope of this paper and the application of leave-one-
out cross validation has not yet been validated for use on high-resolution
data. That NIRSpec has the ability to give meaningful constraints on
individual gas abundances when WFC3 struggled to do so for the same planet
speaks to the exceptional quality of JWST data and the dramatic improvements
in the capabilities of space-based exoplanet spectroscopy in little over a
year’s time.

###  6.3 Comparing the Predictive Power of IGRINS and NIRSpec

The joint IGRINS and NIRSpec retrieval presented in Section 5.4 shows that
these two datasets are a powerful combination for atmospheric model parameter
inference. Estimates from this retrieval are more precise than either data set
provided alone, and it appears that the composition inferences were primarily
driven by IGRINS while the P-T inferences were driven by NIRSpec. In this
subsection, we will identity in what specific ways the two datasets are
complementary and how they support the shortcomings of each other.

Figure 11: Residual $\Delta$log-likelihood maps from cross-correlating the
solar composition 1D-RCTE template with the IGRINS data after the best fit
model spectra from each of the three main retrievals were divided out. The
maps have been median subtracted, and in contrast to the CCF maps in Figure 3,
there are no peaks near the expected planet velocity, indicating the best fit
spectra matched the true underlying planet signal within the precision of the
IGRINS data.

To compare how well each retrieval predicts the IGRINS data, we divide out the
best fit model spectrum from the IGRINS, NIRSpec, and IGRINS+NIRSpec
retrievals and then calculate cross-correlation maps using the same solar
composition model as in Section 3. If a given model is a good fit, there
should be no significant peak at the planet $K_{P}$ and Vsys, and the maps
will only show the cross-correlation signal of residual noise. This is indeed
the case for all three models (Figure 11). There is no significant difference
between the residual maps themselves, indicating that each of the three best
fit spectra predicted the IGRINS data equally well to within the quality of
the IGRINS data.

Figure 12: Top: The median spectrum from 2000 random draws from the IGRINS
pre-eclipse retrieval (red line) and 1$\sigma$ distribution about the median
(red shaded region). The median spectrum from 2000 draws from the NIRSpec
retrieval is in blue, and the 1$\sigma$ distribution is so tightly around this
median it is not visible in this figure. The NIRSpec data is plotted in
diamonds, with the color indicating the difference in the median log-
likelihood While the predictions from the IGRINS retrieval are largely
consistent with the NIRSpec data, it struggles to predict the data around 4
microns the most. Middle: Similar to the top panel, but now in red is the
median spectrum resulting from combining random draws from the IGRINS gas
posteriors but the NIRSpec P-T posteriors. Bottom: Draws from the combined
IGRINS+NIRSpec retrieval. Similar to the NIRSpec retrieval, the 1$\sigma$
contours are not visible.

To compare how well the IGRINS and NIRSpec retrievals predict the NIRSpec
data, we post-process 2000 random draws from each posterior distribution to
the NIRSpec wavelengths and calculate the by-point log-likelihood for each
draw. Unsurprisingly, the NIRSpec posteriors better predict the data with a
median $\chi^{2}/N=0.77\pm 0.02(N=160)$ and a best fit $\chi^{2}/N=0.72$. The
best fit spectrum from the IGRINS-only retrieval adequately fits the data
($\chi^{2}/N=1.72$), and the distribution of the draws is largely consistent
with the data as well, albeit with a much wider range of $\chi^{2}/N$ ( median
$\chi^{2}/N=7.56\pm 9.93$). Taking the difference in the medians of the by-
point log-likelihood distributions, we can see that the points IGRINS
struggles the most to predict are around 4 microns (dark blue points in the
top panel of Figure 12), at which it tends to underpredict the flux.

This region is outside the CO2 feature at 4.5 microns as well as wavelength
regions in which NIRSpec probes higher altitudes than IGRINS. Plotting the
joint IGRINS+NIRSpec best fit photosphere against the data, we can see that
these wavelengths probe $\sim$100-10-2 bar. Comparing the median retrieved P-T
profiles from IGRINS and NIRSpec in Figure 8, we can see that NIRSpec
estimates the atmosphere to be hotter, and therefore brighter, than IGRINS
does at these pressures. When we again post-process draws from the IGRINS
posterior but replace the P-T draws with those from the NIRSpec posterior, the
fits to the data improve. The median $\chi^{2}$/N reduces to 5.2 and the by-
point differences in log$\mathcal{L}$ reduce as well (middle panel of Figure
12). Therefore, we can conclude that NIRSpec’s ability to probe deep into the
atmosphere and more accurately measure the temperature at those pressures is
what gave it the edge over IGRINS in this particular wavelength region.

The joint retrieval combining the IGRINS and NIRSpec data predicts the NIRSpec
data just as well as the NIRSpec data alone without a significant improvement.
The median $\chi^{2}/N$ is 0.78$\pm$0.03, and the best fit $\chi^{2}/N$ is
0.75. Additionally, the by-point differences in log$\mathcal{L}$ tend to be
negligible (bottom panel of Figure 12). Between the NIRSpec-only and combined
retrievals, the precision on the gas abundances increased by a factor of
$\sim$4-5, while the P-T constraints remained largely the same. We can then
conclude that, at least in this low metallicity case with weak molecular
features, fits to the NIRSpec data are more sensitive to the P-T profile than
the composition, and additional information about the gas abundances did not
improve the fit. This further highlights the greater sensitivity to
composition of IGRINS for the molecules it is sensitive to, which does not
include CO2, which NIRSpec was able to place a more stringent upper limit on.

The slight improvements to the precision on the gas abundances between the
IGRINS-only and combined retrievals can be attributed to the decreased
uncertainty on the P-T profile at pressures less than $\sim$ 1 mbar. The gas
abundances are slightly anti-correlated with the P-T parameters $T_{0}$,
log$P_{1}$, and $\alpha_{1}$, and between the IGRINS-only and combined
retrieval the confidence intervals on all three decreased significantly. As
mentioned above, the constraints on the P-T profile are mostly driven by the
NIRSpec data, and the certainty on the P-T profile it provided propagated back
to decreased uncertainty on the gas abundances themselves.

Ultimately, the combined analysis did not qualitatively change the previous
interpretations of WASP-77A b’s atmosphere by Line et al. (2021) and August et
al. (2023). However, the excellent agreement between IGRINS and NIRspec
greatly increases the confidence in the accuracy of these interpretations as
well as those from the combined analysis itself. The most significant
improvement is the more precise constraints on the P-T profile. This can
enable comparisons to predictions from global circulation models to
interrogate assumptions about the distribution of heat in WASP-77A b’s
atmosphere, but such comparisons are beyond the scope of this paper.

###  6.4 The Effect of Ephemeris Error on Atmospheric Inferences

As shown in Section 4, propagated ephemeris error can have a non-negligible
impact on measurements of an exoplanet’s velocity. As the field of high
resolution exoplanet spectroscopy moves toward the measurement of atmospheric
dynamics through the measurement of wind speeds (Gandhi et al., 2022; Pino et
al., 2022) and velocity offsets between gases (Brogi et al., 2023), it is
crucial to avoid false positives from small eccentricities or ephemeris error.
While our “detection” of a superrotating westward jet is physically
implausible and certainly not reality, it is conceivable that a similar error
in ephemeris could lead to the false detection of a more believable dynamical
result for another hot- or ultra-hot Jupiter atmosphere.

The degeneracy between eccentricity and dynamics has been discussed before by
Pino et al. (2022), who also found asymmetries in their retrieved $K_{P}$
values between pre- and post-eclipse phases for the ultra-hot Jupiter KELT-9
b. When fitting for an eccentricity, they obtain a value much higher than the
upper limit derived from previous photometric studies and conclude that the
orbit is only effectively non-circular. Instead, Pino et al. (2022) suggest
the anomalous Doppler shifts in the planet signal are indeed from rotation and
winds. As they note, the degeneracy between eccentricity and atmospheric
dynamics is difficult to break. In the two cases of WASP-77A b and KELT-9 b,
this degeneracy was able to be broken using physical intuition and prior
information, but more ambiguous cases are bound to occur in the future. For a
simple 1D analysis of an atmosphere like this paper, attempting to break the
degeneracy between ephemeris error, winds, and orbital eccentricity may be
beyond the scope for such a study. If one measures a velocity asymmetry
between pre- and post-eclispe phases, it may be adequate to simply give each
sequence separate $K_{P}$ and $V_{\mathrm{sys}}$ values.

We test whether the effective deviation from a circular orbit due to
propagated ephemeris error is enough to affect the planet signal-to-noise
ratio, and we repeat the calculation of CCF and $\Delta$log$\mathcal{L}$ maps
as described in Section 3. Assuming either a circular or eccentric orbit and
using the midtransit time and period values reported by Cortés-Zuleta et al.
(2020) to determine the phase, we find no significant difference in the CCF
SNR and $\Delta$log$\mathcal{L}$ values compared to when we assumed a circular
orbit and the Bonomo et al. (2017) midtransit time. We also repeat the three-
night IGRINS retrieval performed in Section 5.2 with the Cortés-Zuleta et al.
(2020) period and midtransit time in order to test whether the decrease in
precision on the gas abundances compared to the pre-eclipse retrieval was due
to the planet signal being imperfectly summed along its velocity. The
inferences and precision were unaffected.

###  6.5 Sensitivity of the Post-Eclipse Data to the Number of Principal
Components Removed

Figure 13: Marginalized posterior distributions for the abundances of H2O and
CO as well as posterior draws of the P-T profile for when removing 3 principal
components (red) and 4 principal components (blue) from the post-eclipse data.
For comparison, the same marginalized posteriors for the pre-eclipse data when
12 principal components are removed.

As noted in Section 5.2, the post-eclipse data are highly sensitive to the
number of principal components (PCs) removed in the detrending process, and we
chose to remove 3 for our main reported retrieval results on this sequence.
The CCF SNR is stronger after only removing 3 PCs (9.5, compared to 8.0 when
removing 4 PCs and 9.4 for pre-eclipse, Figure 14), but we discounted removing
so few principal components in our initial analysis because in some orders
telluric artifacts were still visible by eye in the post-SVD matrices.
However, no obvious spurious telluric peaks appear in the CCF map. The choice
to initially remove 4 PCs was also motivated as an attempt to be consistent
with the previous analysis by Line et al. (2021) on the pre-eclipse data.
However, as can be seen in Figure 13, which compares the retrieval results
from the 3 PC and 4 PC cases in red and blue, respectively, removing no more
than 3 PCs is required to get any informative constraints from the post-
eclipse data. While this conforms to intuition that removing more principal
components removes more of the planet signal, this sensitivity is concerning,
especially considering there was no indication of such a dramatic shift in
quality from the small improvement in the CCF SNR, and repeating retrievals
several times to test the number of principal components to remove is a time
intensive process that many studies are likely to forgo.

Figure 14: Top: Cross-correlation map for the post eclipse data calculated as
described in Section 3 when removing the first 4 principal components from the
data. Bottom: Similar cross-correlation map but when removing only the first 3
principal components.

Investigating the fraction of total variance accounted for by each principal
component, we can see that on average, about 1% more variance is projected
onto a given principal component for the post-eclipse nights compared to the
same component from the pre-eclipse night (Figure 15). This is not surprising
as the number of principal components a matrix has is equal to its rank, which
in this case is the number of frames in each sequence, of which both of the
post-eclipse nights had fewer. Therefore, a higher fraction of variance per
component is necessary. Indeed, if we crop the pre-eclipse sequence to the
same amount of frames, the variance contained per component matches that of
the post-eclipse nights (green in Figure 15).

Figure 15: Cumulative variance projected onto each principal component for
each night of IGRINS data. Because the number of principal components is
determined by the rank of a matrix, sequences with fewer frames have fewer
principal components, and therefore a higher fraction of variance is projected
onto a given component compared to a shorted sequence.

About 0.67% more total variance in contained in the first 4 principal
components of the post-eclipse nights than the first 4 of the pre-eclipse
night. It takes the first 12 pre-eclipse principal components to account for
the same amount of variance. Based on injections of the IGRINS best-fit model,
the contributed variance of the planet signal itself is less than this
difference, and it is possible more of the planet signal was removed in the
post-eclipse SVD process due to this difference. To determine whether the
increased fraction of variance removed significantly affected the planet
signal, we performed a retrieval on both the cropped pre-eclipse sequence with
the first 4 principal components removed and the full sequence with the first
12 principal components removed. For the cropped sequence, both the
constraints on the composition and P-T profile were largely unchanged. The
confidence intervals on these quantities were slightly larger compared to the
original retrieval but still smaller than those from the post-eclipse
retrievals. The composition constraints from the 12-PC-removed case were also
similar, and we can conclude that the higher fraction of variance removed per
principal component is not an issue.

For the 12-PC-removed case, the median retrieved P-T profile was more
isothermal, similar to the median post-eclipse profile (Figure 13, top row).
This suggests the possibility that in the 4 PC case, the poorly constrained
gas abundances and isothermal P-T profile are separate issues affecting the
post-eclipse data. Some of the P-T parameters, especially $\alpha_{1}$, are
correlated with the multiplicative scale factor $a$, which is constrained
worse in the pre-eclipse 12 PC retrieval. Compared to the 4 PC retrieval, the
1$\sigma$ confidence intervals on log${}_{10}a$ are about 2 times larger, and
the marginalized posterior distribution is not approximately Gaussian as it
was before. It is likely that after the removal of the first 12 principal
components, we have less information about the overall line contrasts, which
is why $a$ and by extension the P-T profile are more poorly constrained, but
information about relative line contrasts is preserved, allowing the absolute
gas abundances to be estimated to similar precision as the 4 PC retrieval.

While the 3 PC retrieval on the post-eclipse data yields significant
improvement over the 4 PC case and places bounded constraints on the H2O and
CO gas abundances consistent with the pre-eclipse constraints and only a
factor of $\sim$2 worse precision (Figure 13), the P-T profile is still near-
isothermal. The similarity to the 12 PC pre-eclipse retrieval (bounded gas
constraints but poorly constrained P-T profile) while more variance still
being removed compared to the pre-eclipse SVD process suggests that
information about the absolute line contrasts and therefore the P-T profile
could be the “first to go” in the SVD process. However, this may also be a
result of our 1D modeling framework failing to capture 3D atmospheric effects,
such as a changing avgerage P-T profile with visible longitudes (discussed
further in the next subsection). Telluric artifacts are prevalent in all
orders of the the post-SVD matrices when only the first 2 PCs are removed, so
we did not attempt a retrieval for a 2 PC removed case.

Next, we investigate the observing conditions of each night and how well the
components of the SVD capture them. Similar to de Kok et al. (2013), we can
identify correlations with recorded observing conditions, such as the air mass
and humidity, with the left singular values (LSVs) in each order. For each
sequence and each order, we identify linear and quadratic correlations between
the median count value (essentially the continuum level), air mass, and
humidity with the first 3 or 4 LSVs (Figure 16), indicating that the removal
of these first few singular values is indeed removing these from the data.
However, the specific singular vectors that correlate to these quantities
changes order-by-order. This implies that each order may have a unique number
of singular values necessary to remove in the detrending process.

For most orders, the first one or two eigenvectors (right singular vectors)
appear consistent across all three nights. However, in many orders, there
appears to be either an “extra” eigenvector on 12/06 that does not exist for
12/14 and 12/21, or 12/06 has similar eigenvectors but in a different order
(Figure 17). The extra eigenvectors displace what would have been otherwise
common eigenvectors with the other nights to lower singular values. For
example, for the order shown in Figure 17, the third eigenvector on 12/06
closely matches the second on both 12/14 and 12/21 but has been “pushed back”
by that night’s own unique second eigenvector.

Figure 16: The first four left singular vectors for order 26 (1.69-1.71
micron) from the sequence taken on 12/21 compared to the continuum level, air
mass, and humidity. Each column corresponds to the $n$th singular value,
indicated at the top. For each left singular vector, we search for both linear
(blue) and quadratic (red) correlations with these. The correlation
coefficients for linear and quadratic fits are shown, with a value $\geq$0.5
being considered a significant correlation. In this particular order, it
appears that effects from humidity are captured in the first two singular
vectors, whereas the continuum takes the first three and there still remains a
correlation with air mass in the fourth. Figure 17: The first four
eigenvectors (columns) from the SVD on order 26 for each night (rows). In this
particular echelle order, all three nights have a common first eigenvector,
while 12/06 appears to have an “extra” eigenvector associated with the second
singular value. The third eigenvector on 12/06 resembles the second
eigenvectors on 12/14 and 12/21.

Likewise, in orders heavily impacted by telluric features, the eigenvectors on
12/06 and 12/14 closely match and it is 12/21 that can have an “extra”
eigenvector associated with the second or third singular value. Again, the
next eigenvector of 12/21 will match closely with eigenvectors of 12/06 and
12/14 but ranked differently by the singular values. We believe these “extra”
eigenvectors are the cause of the post-eclipse data’s sensitivity to number of
principal components removed. It appears that whenever an order has one of
these extra eigenvectors, the additional information removed contains a large
part of the planet signal.

The trade off between which night and which order this extra vector occurs in
appears to be related to the humidity. The humidity on 12/14 was relatively
stable for the majority of the sequence and hence does not have this issue,
while on 12/06 the humidity was both higher and more variable by 40%. On 12/21
the humidity was the lowest of the three nights but it was also more variable
than 12/14 by 25%. It appears that with high humidity, 12/06 “needs” the extra
eigenvector and the first three singular values in orders with few telluric
lines to remove the effects of humidity, while in orders with many telluric
lines it only needs the first two components and eigenvectors to do so. The
opposite is true on 12/21, where in the case of low humidity, the third and
extra eigenvector is needed to remove the effects of humidity in orders with
heavy telluric contamination while for the other orders only the first two are
needed.

Because the extra eigenvector switches between orders and between nights,
selecting the number of principal components to remove both by-night and by-
order may be warranted in the future. Additionally, to obtain a more stable
set of eigenvectors to remove, one could perform SVD on all three sequences
concatenated together or some other large matrix of IGRINS data and then
remove them via multilinear regression such as in Lafarga et al. (2023). This
approach might be more successful and robust in future analyses.

For such a more finely tuned approach to selecting the number of principal
components to remove, it would be appropriate to develop a quantitative metric
to indicate how many is “enough.” This is not trivial – CCF signals can be
spuriously over-optimized to particular models or velocities (see e.g., Cabot
et al., 2019; Cheverall et al., 2023), and searching for correlations between
the LSVs and physical components of the data requires a comprehensive record
of such components. Air mass and humidity appear to be prevalent quantities
captured by the SVD in this data, but information about e.g., the seeing or
dew point were not available yet these may still impact the data.

####  6.5.1 Salvaging the Post-eclipse Data By Combining with NIRSpec

Figure 18: Marginalized posterior distributions of relevant parameters from
retrievals on the post-eclipse data with 4 principal components removed (red),
NIRSpec (blue), and the two combined (purple). Only lower limits would be
placed on the abundances of H2O and CO from the post-eclipse data alone. By
combining it with the NIRSpec data, we are able to use what gas abundance
information there was in that data and gain more precise estimates for these
values than afforded by NIRSpec alone.

The retrieval on the post-eclipse data with 4 principal components removed
struggled to produce informative results in large part due to the inability to
break the degeneracy between the gas abundances and certain parameters for the
analytic P-T profile. In principal, if information about the temperature in
the deep atmosphere was regained and broke this degeneracy, informative
inferences about the composition could be made from the post-eclipse data in
the hypothetical scenario in which removing less principal components did not
improve the signal. The NIRSpec data contains such information, so we combined
it with the 4 PC-removed post-eclipse data in a single retrieval similar to
the combination with the pre-eclipse data in Section 5.4. This retrieval is
able to place bounded constraints on H2O and CO to precision slightly better
than achieved with the NIRSpec data alone but not quite at the precision the
pre-eclipse data provided (this exercise is summarized in Figure 18). This
test confirms that once the abundance and P-T degeneracy was broken, we could
extract gas abundance information from the post-eclipse data that improved
upon the constraints possible with NIRSpec alone. Therefore, even if a high
resolution data set seems fruitless, combining it with low-resolution data can
unlock previously inaccessible information.

###  6.6 Evidence for Thermal Inhomogeneity

Figure 19: Left: Marginalized posterior distributions of the H2O and CO gas
abundances from retrievals on the full post-eclipse sequence (red), the first
half only (blue), and the second half only (purple). Right: Median retrieved
day- and night-side P-T profiles from the 2 P-T retrieval in Section 6.6.

As can be seen in the CCF orbit trail in Section 3, the post-eclipse signal
disappears around phase 0.57. On both post-eclipse nights, the air mass
increased over the course of the night, and the subsequent decrease in total
flux of the planet may be reflected in the orbit trail. However, air mass is
similarly increasing with time during the pre-eclipse night, yet the planet
signal appears to be robust for the entire pre-eclipse sequence, including
phases further away from eclipse (i.e., $\sim$100% dayside visibility) than
those covered in the post-eclispe data. Therefore, the strength of the CCF
trail with time may instead be reflecting a physical change in the planet
signal itself.

The morning terminator of the planet is rotating into view during the post-
eclipse sequences. If WASP-77A b has efficient heat redistribution, which the
retrieved P-T profile is consistent with (Line et al., 2021), the hotspot
should be offset to the east (evening) and the morning (western) terminator
should be cooler than the evening terminator, which is visible in pre-eclipse
phases. There is precedence for hotspot offsets affecting HRCCS data, such as
Herman et al. (2022) or van Sluijs et al. (2023), both of which inferred
thermal inhomogeneities via a phase dependence of line contrasts.

The decrease in CCF signal may simply be from the hotspot rotating out of
view, or it may also be due to patchy cloud coverage. While the day sides of
hot Jupiters are too hot for clouds to form, it is expected that clouds will
form on their night sides before blowing over and vaporizing again on the day
side (Parmentier et al., 2016). WASP-77A b’s night side would be around the
same temperature as the L/T-dwarf transition, which is marked by the
appearance of cloud species such as forsterite and enstatite. WASP-77A b’s UV
transit depth is also consistent with silicate clouds (Turner et al., 2016).
Therefore, it is possible that clouds form on the night side and
heterogeneously cover the morning limb and other parts of the western
hemisphere and are affecting the IGRINS data.

As a first order test of whether we are sensitive to longitudinal thermal
variations, we split the post-eclipse sequence into two halves demarcated at
phase 0.57 and performed retrievals on each half individually. The constraints
from the first half are virtually identical to those from the entire sequence,
while for the second half only a lower limit can be placed on the abundance of
H2O and CO is entirely unconstrained (Figure 19). We can conclude that the
full sequence constraints were driven almost entirely by the first half alone.
This suggests we are indeed probing thermal gradients in the post-eclipse
sequence and a longitudinally varying treatment of the atmosphere is
warranted.

To test both (1) whether the increasing fraction of the cooler morning
terminator and night side as we see different longitudes is “diluting” the
observed thermal emission spectrum and (2) whether clouds on these cooler
parts of the planet are further obscuring the deep atmosphere, we repeat the
retrieval analysis with a simplified 2 P-T atmosphere (such as in e.g., Feng
et al., 2016). The total outgoing spectrum is a phase dependent linear
combination of a day and night side spectrum:

$F(t)=k(t)F_{\mathrm{day}}+(1-k(t))F_{\mathrm{night}}$ (4)

where $k$ is the fraction of the day side that is visible at a given phase
angle $\alpha$:

$k(t)=\frac{1}{2}(1+\cos\alpha(t));\ \alpha(t)=2\pi\varphi(t)+\pi$ (5)

Both the day and night side have an individual free P-T profile and a gray
cloud opacity while we enforce that the global bulk composition of H2O and CO
remains constant and therefore are each single free parameters for the entire
planet (see e.g., Cooper & Showman, 2006 for a discussion on the global
chemical homogenization on hot Jupiters). With this framework, we are unable
to measure a nightside cloud opacity, but we are able to measure two separate
P-T profiles, with the night side profile cooler than the dayside by several
standard deviations. The inclusion of the 2 P-T profiles is moderately favored
over the baseline 1 P-T profile retrieval with a Bayes factor of 16.54
(2.88$\sigma$). Further analysis of the three dimensional nature of WASP-77A
b’s atmosphere is beyond the scope of this paper, but our 1D retrieval results
are likely not significantly biased by the night side because the
multiplicative scale factor $a$ can adequately account for dilution by thermal
inhomogeneities (Taylor et al., 2020). Nonetheless, our 1D retrieval models
failing to capture underlying 3D effects might be one of several factors
(along with lower SNR and the issues with detrending) contributing to the
slightly less precise inferences offered by the post-eclipse data than
compared to the pre-eclipse sequence.

###  6.7 Conflicting Evidence for the Presence of 13CO?

The measurement of the 13CO/12CO isotopologue ratio first presented in Line et
al. (2021)’s analysis of the pre-eclipse IGRINS data was an intriguing result.
This was the first such measurement in a transiting exoplanet atmosphere, and
isotope ratios can potentially provide another avenue for shedding light on a
planet’s formation history (Pontoppidan et al., 2014). We are able to recreate
this measurement in the pre-eclipse data, and the inclusion of 13CO is favored
by 4.2$\sigma$, but we are unable to measure the isotope ratio from the post-
eclipse data alone and inclusion of 13CO is neither statistically favored or
disfavored.

To test the sensitivity of the post-eclipse data to the isotopologue, we
divide out of the data the best fit model from the 3-night IGRINS retrieval,
then inject the same model but with the pre-eclipse isotope ratio back into
the data at an offset $V_{\mathrm{sys}}$, and perform a retrieval again. We
are still unable to place a bounded measurement on the isotopologue ratio,
indicating that the quality of the post-eclipse data is preventing us from
confirming the pre-eclipse measurement.

Similarly, from the NIRSpec data we are unable to measure the isotopologue
ratio and instead place an upper limit consistent with the pre-eclipse IGRINS
data at 1/1.48. To test NIRSpec’s sensitivity to the presence of 13CO, we
perform a retrieval on the best-fit model from the three-night IGRINS
retrieval. We post-process the model spectrum to the NIRSpec wavelength range,
bin onto the data wavelength grid, then add a noise instance based on the
error bars. The input value for 13CO/12CO was 1/22, and we again can only
place an upper limit at 1/5.88, similar to the retrieval on the real NIRSpec
data. Therefore, we conclude that the NIRSpec data is also not sensitive
enough to definitively rule out the presence of 13CO.

##  7 Conclusions

We have presented the first combined exoplanet retrieval analysis using both
ground-based and JWST data as well as two new nights of IGRINS data covering
the post-eclipse phases of WASP-77A b. Our findings are as follows:

  * •

In addition to the IGRINS pre-eclipse day side thermal emission data first
presented in Line et al. (2021), we present here two additional nights
capturing the planet’s day side emission in the post-eclipse phases. Using
traditional cross-correlation methods, we detect H2O and CO absorption
features in these additional data.

  * •

We found no signatures of atmospheric dynamics and WASP-77A b’s orbit is
effectively circular. However, propagated ephemeris error manifested as
measurements of both an effective eccentricity and an implausibly strong
westward equatorial jet. Using updated midtransit times decreases the chances
of false positive detections of dynamics.

  * •

The additional post-eclipse IGRINS data is highly sensitive to the number of
principal components removed in the SVD detrending process. Removing 4, as was
done on the pre-eclipse data, removed too much of the planet signal to make
informative inferences about WASP-77A b’s atmosphere. Only by removing 3 could
informative inferences on its composition be made. Investigation of the
individual singular vectors produced by the SVD indicates that high and
variable humidity played a large role in what information was projected onto a
given principal component in a given order. Future observations would benefit
from low or less variable humidity or a by-order choice of number of principal
components to remove.

  * •

We find moderate evidence for thermal inhomogeneity in the post-eclipse data
and are able to retrieve both a hot dayside and cooler nightside P-T profile.

  * •

The IGRINS and NIRSpec data are in excellent agreement with each other,
lending much more confidence in the analysis of each. We reproduce the sub-
solar metallicity as measured by Line et al. (2021) and August et al. (2023)
and the solar C/O ratio as measured by Line et al. (2021). Recent studies
suggest this combination may indicate formation beyond the CO2 ice line, but
formation within the H2O ice line is also possible within different formation
modeling frameworks.

  * •

Combining the pre-eclipse IGRINS data with the NIRSpec/G395H data allowed us
to more stringently constrain both the composition and vertical thermal
structure of WASP-77A b than possible with either data set alone. IGRINS is
more sensitive to absolute gas abundances, while NIRSpec is more sensitive to
the P-T profile.

  * •

Our inferences from the IGRINS data alone largely predict the NIRSpec data
well. Wavelengths in which IGRINS struggled to predict the NIRSpec data
coincided with pressures both instruments probe, indicating that NIRSpec is
able to estimate the P-T profile more accurately due to the preservation of
continuum information in low resolution data.

  * •

Neither the post-eclipse or NIRSpec data are sensitive enough to the 13CO/12CO
isotopologue ratio to either refute or confirm the measurement made with the
pre-eclipse data.

As shown by its increased sensitivity to absolute gas abundances and ability
to search for atmospheric dynamics, high resolution instruments like IGRINS
are still relevant and necessary for exoplanet atmospheric science even with
the advent of JWST. Combinations with low resolution transit spectroscopy can
provide more powerful probes of exoplanetary atmospheres than either method
can provide alone thanks to the complementary nature of the physical
quantities each method is sensitive to. High resolution observations in
emission are essentially partial spectroscopic phase curves, and as shown in
this paper, there is potential to uncover three dimensional information about
the thermochemical structure of hot Jupiters. However, thermochemical
gradients of hot and ultra-hot Jupiters are largely unexplored in high
resolution thermal emission studies. Dynamics, which are measurable with
HRCCS, and thermal structure, measurable with low-resolution instruments, are
intrinsically and physically linked, but the utility of combining high- and
low-resolution data for probing thermochemical gradients remains to be seen.

## Acknowledgements

M.R.L. and J.L.B. acknowledge support for this work from NSF grant
AST-2307177. MM acknowledges support provided by NASA through the NASA Hubble
Fellowship grant HST-HF2-51485.001-A awarded by the Space Telescope Science
Institute, which is operated by AURA, Inc., for NASA, under contract
NAS5-26555. Based on observations obtained at the international Gemini
Observatory, a program of NSF’s NOIRLab, which is managed by the Association
of Universities for Research in Astronomy (AURA) under a cooperative agreement
with the National Science Foundation on behalf of the Gemini Observatory
partnership: the National Science Foundation (United States), National
Research Council (Canada), Agencia Nacional de Investigación y Desarrollo
(Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina),
Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and
Korea Astronomy and Space Science Institute (Republic of Korea). This work
used the Immersion Grating Infrared Spectrometer (IGRINS) that was developed
under a collaboration between the University of Texas at Austin and the Korea
Astronomy and Space Science Institute (KASI) with the financial support of the
US National Science Foundation 27 under grants AST-1229522 and AST-1702267, of
the University of Texas at Austin, and of the Korean GMT Project of KASI. This
publication makes use of data products from Exoplanet Watch, a citizen science
project managed by NASA’s Jet Propulsion Laboratory on behalf of NASA’s
Universe of Learning. This work is supported by NASA under award number
NNX16AC65A to the Space Telescope Science Institute, in partnership with
Caltech/IPAC, Center for Astrophysics—Harvard & Smithsonian, and NASA Jet
Propulsion Laboratory.

## References

  * Arcangeli et al. (2018) Arcangeli, J., Désert, J.-M., Line, M. R., et al. 2018, ApJ, 855, L30, doi: 10.3847/2041-8213/aab272
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222
  * Atreya et al. (2016) Atreya, S. K., Crida, A., Guillot, T., et al. 2016, arXiv e-prints, arXiv:1606.04510.  https://arxiv.org/abs/1606.04510
  * August et al. (2023) August, P. C., Bean, J. L., Zhang, M., et al. 2023, ApJ, 953, L24, doi: 10.3847/2041-8213/ace828
  * Azzam et al. (2016) Azzam, A. A. A., Tennyson, J., Yurchenko, S. N., & Naumenko, O. V. 2016, MNRAS, 460, 4063, doi: 10.1093/mnras/stw1133
  * Barber et al. (2014) Barber, R. J., Strange, J. K., Hill, C., et al. 2014, MNRAS, 437, 1828, doi: 10.1093/mnras/stt2011
  * Baxter et al. (2020) Baxter, C., Désert, J.-M., Parmentier, V., et al. 2020, A&A, 639, A36, doi: 10.1051/0004-6361/201937394
  * Bean et al. (2023) Bean, J. L., Xue, Q., August, P. C., et al. 2023, Nature, 618, 43, doi: 10.1038/s41586-023-05984-y
  * Beltz et al. (2021) Beltz, H., Rauscher, E., Brogi, M., & Kempton, E. M. R. 2021, AJ, 161, 1, doi: 10.3847/1538-3881/abb67b
  * Birkby (2018) Birkby, J. L. 2018, arXiv e-prints, arXiv:1806.04617.  https://arxiv.org/abs/1806.04617
  * Birkby et al. (2013) Birkby, J. L., de Kok, R. J., Brogi, M., et al. 2013, MNRAS, 436, L35, doi: 10.1093/mnrasl/slt107
  * Bitsch et al. (2022) Bitsch, B., Schneider, A. D., & Kreidberg, L. 2022, A&A, 665, A138, doi: 10.1051/0004-6361/202243345
  * Bonomo et al. (2017) Bonomo, A. S., Desidera, S., Benatti, S., et al. 2017, A&A, 602, A107, doi: 10.1051/0004-6361/201629882
  * Boucher et al. (2023) Boucher, A., Lafreniére, D., Pelletier, S., et al. 2023, MNRAS, 522, 5062, doi: 10.1093/mnras/stad1247
  * Brogi et al. (2016) Brogi, M., de Kok, R. J., Albrecht, S., et al. 2016, ApJ, 817, 106, doi: 10.3847/0004-637X/817/2/106
  * Brogi & Line (2019) Brogi, M., & Line, M. R. 2019, AJ, 157, 114, doi: 10.3847/1538-3881/aaffd3
  * Brogi et al. (2012) Brogi, M., Snellen, I. A. G., de Kok, R. J., et al. 2012, Nature, 486, 502, doi: 10.1038/nature11161
  * Brogi et al. (2023) Brogi, M., Emeka-Okafor, V., Line, M. R., et al. 2023, AJ, 165, 91, doi: 10.3847/1538-3881/acaf5c
  * Buchner (2016) Buchner, J. 2016, PyMultiNest: Python interface for MultiNest.  http://ascl.net/1606.005
  * Burrows & Sharp (1999) Burrows, A., & Sharp, C. M. 1999, ApJ, 512, 843, doi: 10.1086/306811
  * Cabot et al. (2019) Cabot, S. H. C., Madhusudhan, N., Hawker, G. A., & Gandhi, S. 2019, MNRAS, 482, 4422, doi: 10.1093/mnras/sty2994
  * Cheverall et al. (2023) Cheverall, C. J., Madhusudhan, N., & Holmberg, M. 2023, MNRAS, 522, 661, doi: 10.1093/mnras/stad648
  * Coles et al. (2019) Coles, P. A., Yurchenko, S. N., & Tennyson, J. 2019, MNRAS, 490, 4638, doi: 10.1093/mnras/stz2778
  * Cooper & Showman (2006) Cooper, C. S., & Showman, A. P. 2006, ApJ, 649, 1048, doi: 10.1086/506312
  * Cortés-Zuleta et al. (2020) Cortés-Zuleta, P., Rojo, P., Wang, S., et al. 2020, A&A, 636, A98, doi: 10.1051/0004-6361/201936279
  * Coulombe et al. (2023) Coulombe, L.-P., Benneke, B., Challener, R., et al. 2023, Nature, 620, 292, doi: 10.1038/s41586-023-06230-1
  * de Kok et al. (2013) de Kok, R. J., Brogi, M., Snellen, I. A. G., et al. 2013, A&A, 554, A82, doi: 10.1051/0004-6361/201321381
  * Eastman et al. (2013) Eastman, J., Gaudi, B. S., & Agol, E. 2013, PASP, 125, 83, doi: 10.1086/669497
  * Ehrenreich et al. (2020) Ehrenreich, D., Lovis, C., Allart, R., et al. 2020, Nature, 580, 597, doi: 10.1038/s41586-020-2107-1
  * Feng et al. (2016) Feng, Y. K., Line, M. R., Fortney, J. J., et al. 2016, ApJ, 829, 52, doi: 10.3847/0004-637X/829/1/52
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601, doi: 10.1111/j.1365-2966.2009.14548.x
  * Flowers et al. (2019) Flowers, E., Brogi, M., Rauscher, E., Kempton, E. M. R., & Chiavassa, A. 2019, AJ, 157, 209, doi: 10.3847/1538-3881/ab164c
  * Fortney et al. (2008) Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419, doi: 10.1086/528370
  * Fulton et al. (2018) Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504, doi: 10.1088/1538-3873/aaaaa8
  * Gandhi et al. (2022) Gandhi, S., Kesseli, A., Snellen, I., et al. 2022, MNRAS, 515, 749, doi: 10.1093/mnras/stac1744
  * Gandhi et al. (2019) Gandhi, S., Madhusudhan, N., Hawker, G., & Piette, A. 2019, AJ, 158, 228, doi: 10.3847/1538-3881/ab4efc
  * Gandhi et al. (2023) Gandhi, S., Kesseli, A., Zhang, Y., et al. 2023, AJ, 165, 242, doi: 10.3847/1538-3881/accd65
  * Gharib-Nezhad et al. (2021) Gharib-Nezhad, E., Iyer, A. R., Line, M. R., et al. 2021, ApJS, 254, 34, doi: 10.3847/1538-4365/abf504
  * Gibson et al. (2020) Gibson, N. P., Merritt, S., Nugroho, S. K., et al. 2020, MNRAS, 493, 2215, doi: 10.1093/mnras/staa228
  * Hargreaves et al. (2020) Hargreaves, R. J., Gordon, I. E., Rey, M., et al. 2020, ApJS, 247, 55, doi: 10.3847/1538-4365/ab7a1a
  * Hawker et al. (2018) Hawker, G. A., Madhusudhan, N., Cabot, S. H. C., & Gandhi, S. 2018, ApJ, 863, L11, doi: 10.3847/2041-8213/aac49d
  * Haynes et al. (2015) Haynes, K., Mandell, A. M., Madhusudhan, N., Deming, D., & Knutson, H. 2015, ApJ, 806, 146, doi: 10.1088/0004-637X/806/2/146
  * Herman et al. (2022) Herman, M. K., de Mooij, E. J. W., Nugroho, S. K., Gibson, N. P., & Jayawardhana, R. 2022, AJ, 163, 248, doi: 10.3847/1538-3881/ac5f4d
  * Hoeijmakers et al. (2019) Hoeijmakers, H. J., Ehrenreich, D., Kitzmann, D., et al. 2019, A&A, 627, A165, doi: 10.1051/0004-6361/201935089
  * Householder & Weiss (2022) Householder, A., & Weiss, L. 2022, arXiv e-prints, arXiv:2212.06966, doi: 10.48550/arXiv.2212.06966
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Husser et al. (2013) Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6, doi: 10.1051/0004-6361/201219058
  * JWST Transiting Exoplanet Community Early Release Science Team et al. (2023) JWST Transiting Exoplanet Community Early Release Science Team, Ahrer, E.-M., Alderson, L., et al. 2023, Nature, 614, 649, doi: 10.1038/s41586-022-05269-w
  * Karman et al. (2019) Karman, T., Gordon, I. E., van der Avoird, A., et al. 2019, Icarus, 328, 160, doi: 10.1016/j.icarus.2019.02.034
  * Kasper et al. (2021) Kasper, D., Bean, J. L., Line, M. R., et al. 2021, ApJ, 921, L18, doi: 10.3847/2041-8213/ac30e1
  * Kasper et al. (2023) —. 2023, AJ, 165, 7, doi: 10.3847/1538-3881/ac9f40
  * Khorshid et al. (2023) Khorshid, N., Min, M., & Désert, J. M. 2023, A&A, 675, A95, doi: 10.1051/0004-6361/202245469
  * Khorshid et al. (2021) Khorshid, N., Min, M., Désert, J. M., Woitke, P., & Dominik, C. 2021, arXiv e-prints, arXiv:2111.00279.  https://arxiv.org/abs/2111.00279
  * Kokori et al. (2022) Kokori, A., Tsiaras, A., Edwards, B., et al. 2022, Experimental Astronomy, 53, 547, doi: 10.1007/s10686-020-09696-3
  * Kolecki & Wang (2022) Kolecki, J. R., & Wang, J. 2022, AJ, 164, 87, doi: 10.3847/1538-3881/ac7de3
  * Kreidberg et al. (2014) Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, ApJ, 793, L27, doi: 10.1088/2041-8205/793/2/L27
  * Lafarga et al. (2023) Lafarga, M., Brogi, M., Gandhi, S., et al. 2023, MNRAS, 521, 1233, doi: 10.1093/mnras/stad480
  * Lee & Gullikson (2016) Lee, J.-J., & Gullikson, K. 2016, Plp: V2.1 Alpha 3, v2.1-alpha.3, Zenodo, doi: 10.5281/zenodo.56067
  * Li et al. (2015) Li, G., Gordon, I. E., Rothman, L. S., et al. 2015, ApJS, 216, 15, doi: 10.1088/0067-0049/216/1/15
  * Line et al. (2013) Line, M. R., Wolf, A. S., Zhang, X., et al. 2013, ApJ, 775, 137, doi: 10.1088/0004-637X/775/2/137
  * Line et al. (2021) Line, M. R., Brogi, M., Bean, J. L., et al. 2021, Nature, 598, 580, doi: 10.1038/s41586-021-03912-6
  * Lodders & Fegley (2002) Lodders, K., & Fegley, B. 2002, Icarus, 155, 393, doi: 10.1006/icar.2001.6740
  * Mace et al. (2018) Mace, G., Sokal, K., Lee, J.-J., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII, ed. C. J. Evans, L. Simard, & H. Takami, 107020Q, doi: 10.1117/12.2312345
  * Madhusudhan et al. (2014) Madhusudhan, N., Amin, M. A., & Kennedy, G. M. 2014, ApJ, 794, L12, doi: 10.1088/2041-8205/794/1/L12
  * Madhusudhan et al. (2017) Madhusudhan, N., Bitsch, B., Johansen, A., & Eriksson, L. 2017, MNRAS, 469, 4102, doi: 10.1093/mnras/stx1139
  * Madhusudhan & Seager (2009) Madhusudhan, N., & Seager, S. 2009, ApJ, 707, 24, doi: 10.1088/0004-637X/707/1/24
  * Mansfield et al. (2021) Mansfield, M., Line, M. R., Bean, J. L., et al. 2021, Nature Astronomy, 5, 1224, doi: 10.1038/s41550-021-01455-4
  * Mansfield et al. (2022) Mansfield, M., Wiser, L., Stevenson, K. B., et al. 2022, AJ, 163, 261, doi: 10.3847/1538-3881/ac658f
  * Maxted et al. (2013) Maxted, P. F. L., Anderson, D. R., Collier Cameron, A., et al. 2013, PASP, 125, 48, doi: 10.1086/669231
  * Meibom et al. (2007) Meibom, A., Krot, A. N., Robert, F., et al. 2007, ApJ, 656, L33, doi: 10.1086/512052
  * Mordasini et al. (2016) Mordasini, C., van Boekel, R., Mollière, P., Henning, T., & Benneke, B. 2016, ApJ, 832, 41, doi: 10.3847/0004-637X/832/1/41
  * Moses et al. (2013) Moses, J. I., Line, M. R., Visscher, C., et al. 2013, ApJ, 777, 34, doi: 10.1088/0004-637X/777/1/34
  * Mousis et al. (2019) Mousis, O., Ronnet, T., & Lunine, J. I. 2019, ApJ, 875, 9, doi: 10.3847/1538-4357/ab0a72
  * Öberg & Bergin (2016) Öberg, K. I., & Bergin, E. A. 2016, ApJ, 831, L19, doi: 10.3847/2041-8205/831/2/L19
  * Öberg et al. (2011) Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16, doi: 10.1088/2041-8205/743/1/L16
  * Park et al. (2014) Park, C., Jaffe, D. T., Yuk, I.-S., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V, ed. S. K. Ramsay, I. S. McLean, & H. Takami, 91471D, doi: 10.1117/12.2056431
  * Parmentier et al. (2016) Parmentier, V., Fortney, J. J., Showman, A. P., Morley, C., & Marley, M. S. 2016, ApJ, 828, 22, doi: 10.3847/0004-637X/828/1/22
  * Parmentier et al. (2018) Parmentier, V., Line, M. R., Bean, J. L., et al. 2018, A&A, 617, A110, doi: 10.1051/0004-6361/201833059
  * Pelletier et al. (2021) Pelletier, S., Benneke, B., Darveau-Bernier, A., et al. 2021, AJ, 162, 73, doi: 10.3847/1538-3881/ac0428
  * Pelletier et al. (2023) Pelletier, S., Benneke, B., Ali-Dib, M., et al. 2023, arXiv e-prints, arXiv:2306.08739, doi: 10.48550/arXiv.2306.08739
  * Pino et al. (2022) Pino, L., Brogi, M., Désert, J. M., et al. 2022, A&A, 668, A176, doi: 10.1051/0004-6361/202244593
  * Piskorz et al. (2018) Piskorz, D., Buzard, C., Line, M. R., et al. 2018, AJ, 156, 133, doi: 10.3847/1538-3881/aad781
  * Polanski et al. (2022) Polanski, A. S., Crossfield, I. J. M., Howard, A. W., Isaacson, H., & Rice, M. 2022, Research Notes of the American Astronomical Society, 6, 155, doi: 10.3847/2515-5172/ac8676
  * Polyansky et al. (2018) Polyansky, O. L., Kyuberis, A. A., Zobov, N. F., et al. 2018, MNRAS, 480, 2597, doi: 10.1093/mnras/sty1877
  * Pontoppidan et al. (2014) Pontoppidan, K. M., Salyk, C., Bergin, E. A., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 363–385, doi: 10.2458/azu_uapress_9780816531240-ch016
  * Reggiani et al. (2022) Reggiani, H., Schlaufman, K. C., Healy, B. F., Lothringer, J. D., & Sing, D. K. 2022, AJ, 163, 159, doi: 10.3847/1538-3881/ac4d9f
  * Rothman et al. (2010) Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010, J. Quant. Spec. Radiat. Transf., 111, 2139, doi: 10.1016/j.jqsrt.2010.05.001
  * Schneider & Bitsch (2021) Schneider, A. D., & Bitsch, B. 2021, A&A, 654, A71, doi: 10.1051/0004-6361/202039640
  * Sing et al. (2016) Sing, D. K., Fortney, J. J., Nikolov, N., et al. 2016, Nature, 529, 59, doi: 10.1038/nature16068
  * Snellen et al. (2010) Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049, doi: 10.1038/nature09111
  * Southworth (2011) Southworth, J. 2011, MNRAS, 417, 2166, doi: 10.1111/j.1365-2966.2011.19399.x
  * Taylor et al. (2020) Taylor, J., Parmentier, V., Irwin, P. G. J., et al. 2020, MNRAS, 493, 4342, doi: 10.1093/mnras/staa552
  * Taylor et al. (2023) Taylor, J., Radica, M., Welbanks, L., et al. 2023, MNRAS, 524, 817, doi: 10.1093/mnras/stad1547
  * Tennyson et al. (2020) Tennyson, J., Yurchenko, S. N., Al-Refaie, A. F., et al. 2020, J. Quant. Spec. Radiat. Transf., 255, 107228, doi: 10.1016/j.jqsrt.2020.107228
  * Turner et al. (2016) Turner, J. D., Pearson, K. A., Biddle, L. I., et al. 2016, MNRAS, 459, 789, doi: 10.1093/mnras/stw574
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * van Sluijs et al. (2023) van Sluijs, L., Birkby, J. L., Lothringer, J., et al. 2023, MNRAS, 522, 2145, doi: 10.1093/mnras/stad1103
  * Virtanen et al. (2019) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2019, arXiv e-prints, arXiv:1907.10121.  https://arxiv.org/abs/1907.10121
  * Welbanks et al. (2019) Welbanks, L., Madhusudhan, N., Allard, N. F., et al. 2019, ApJ, 887, L20, doi: 10.3847/2041-8213/ab5a89
  * Welbanks et al. (2023) Welbanks, L., McGill, P., Line, M., & Madhusudhan, N. 2023, AJ, 165, 112, doi: 10.3847/1538-3881/acab67
  * Wilson & Rood (1994) Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191, doi: 10.1146/annurev.aa.32.090194.001203
  * Wong et al. (2020) Wong, I., Shporer, A., Daylan, T., et al. 2020, AJ, 160, 155, doi: 10.3847/1538-3881/ababad
  * Zahnle et al. (2009) Zahnle, K., Marley, M. S., Freedman, R. S., Lodders, K., & Fortney, J. J. 2009, ApJ, 701, L20, doi: 10.1088/0004-637X/701/1/L20
  * Zellem et al. (2020) Zellem, R. T., Pearson, K. A., Blaser, E., et al. 2020, PASP, 132, 054401, doi: 10.1088/1538-3873/ab7ee7
  * Zhang et al. (2021) Zhang, Y., Snellen, I. A. G., Bohn, A. J., et al. 2021, Nature, 595, 370, doi: 10.1038/s41586-021-03616-x

##  Appendix A Eccentricity and Wind Models

In Section 4 we introduced two alternate versions of the planet velocity,
$V_{P}(t)$ to use in Eqn. 1. The first is an eccentric orbit, for which
$V_{P}(t)$ becomes:

$V_{P}(t)=K_{P}[\cos(\nu(t)+\omega_{P})+e\cos(\omega_{P})]$ (A1)

where $\nu$(t) is the true anomaly, $\omega_{P}$ is the argument of periastron
of the planet, and $e$ is the eccentricity. $\nu(t)$ is defined as

$\nu(t)=E(t)+2\arctan\frac{\beta\sin E(t)}{1-\beta\cos E(t)};\
\beta=\frac{e}{1+\sqrt{1-e^{2}}}$ (A2)

where $E(t)$ is the eccentric anomaly, which we use the Newton-Raphson method
to solve for from the mean anomaly $M(t)$:

$M(t)=\frac{2\pi}{P}(t-T_{P})=E(t)-e\sin E(t)$ (A3)

where $T_{P}$ is the time of periastron.

$\omega_{P}$ should not be confused with $\omega_{\star}$, the argument of
periastron of the star666The many RV codes in the literature are inconsistent
between which $\omega$ is used and the sign of line-of-sight velocity in Eqn.
A1. For example, Cortés-Zuleta et al. (2020) use ExoFast (Eastman et al.,
2013) and report $\omega_{\star}$, so we add 180∘ to this value to get
$\omega_{P}$. This subtlety caused the authors many headaches. See Householder
& Weiss (2022) for a more detailed overview of this problem. – the two are
different by 180∘. To be consistent with Cortés-Zuleta et al. (2020), we use
$\omega_{\star}$ as the input value for our RV code and convert to
$\omega_{P}$ for Eqn. A1. It is common in the literature to not fit directly
for $e$ and $\omega_{\star}$ but instead for the two quantities
$A=\sqrt{e}\cos\omega_{\star}$ and $B=\sqrt{e}\sin\omega_{\star}$ (see e.g.,
Eastman et al., 2013 and Fulton et al., 2018 for discussion on choice of RV
parameterization). $e$ and $\omega_{\star}$ are obtained from the two via:

$e=A^{2}+B^{2};\ \omega_{\star}=\arctan\frac{B}{A}$ (A4)

For both of these we set a uniform prior from -1 to 1.

For the jet model, we find the disk averaged line-of-sight velocity from a
function of both longitude $\theta$ and latitude $\phi$ on the visible disk:

$v_{\mathrm{LOS}}(\theta,\phi)=u\sin\theta+v_{\mathrm{rot}}\sin\theta\cos\phi$
(A5)

where $u$ is the jet speed and $v_{\mathrm{rot}}$ is the equatorial rotation
velocity, which we set to 4.52 km s-1 assuming tidal locking. At a given
orbital phase, we calculate $v_{\mathrm{LOS}}$ on a grid of of disk longitudes
and latitudes overlapping with the visible dayside then take the weighted
average where the weights are

$\mu(\theta,\phi)=\cos\theta\cos\phi$ (A6)

and the average line-of-sight velocity is then:

$\langle
v_{\mathrm{LOS}}\rangle=\frac{\sum_{i,j}\mu(\theta_{i},\phi_{j})v_{\mathrm{LOS}}(\theta_{i},\phi_{j})}{\sum_{i,j}\mu(\theta_{i},\phi_{j})}$
(A7)

This is added to the velocity from a circular orbit:

$V_{P}(t)=K_{P}\sin[2\pi\varphi(t)]+\langle
v_{\mathrm{LOS}}(\varphi(t))\rangle$ (A8)

##  Appendix B Corner Plots

The associated corner plots for all of the velocity inferences and retrieval
analyses are available in a public Zenodo repository linked here. Note: by
default our plotting routine lists the marginalized posterior medians and
1$\sigma$ confidence intervals, even if the posterior distribution is against
a prior bound. The utility of corner plots is to qualitatively inspect the
correlation between model parameters. For quantitative estimates, refer to
those listed in Table 3 or in the text.