[1]Matthew M. Murphy

[1]Steward Observatory, University of Arizona, Cherry Ave., Tucson, 85705, AZ,
USA

2]Department of Astronomy, University of Wisconsin–Madison, Madison, WI, USA

3]Bay Area Environmental Research Institute, NASA’s Ames Research Center,
Moffet Field, CA, USA

4]Space Science and Astrobiology Division, NASA’s Ames Research Center, Moffet
Field, CA, USA

5]School of Earth and Space Exploration, Arizona State University, Tempe, AZ,
USA

6]Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire
Lagrange, France

7]Department of Astronomy, University of Michigan, Ann Arbor, MI, USA

8]Department of Astronomy and Astrophysics, University of California Santa
Cruz, Santa Cruz, CA, USA

# Evidence for Morning-to-Evening Limb Asymmetry on the Cool Low-density
Exoplanet WASP-107b

<EMAIL_ADDRESS>Thomas G. Beatty  Everett Schlawin  Taylor J. Bell
Michael R. Line  Thomas P. Greene  Vivien Parmentier  Emily Rauscher  Luis
Welbanks  Jonathan J. Fortney  Marcia Rieke  *  [  [  [  [  [  [  [

###### Abstract

Transmission spectroscopy has enabled unprecedented insights into the makeup
of exoplanet atmospheres. A transmission spectrum combines contributions from
a planet’s morning and evening limbs, but these limbs may have different
temperatures, compositions, and aerosol properties due to atmospheric
circulation. High-resolution ground-based observations have detected limb
asymmetry on several ultra-hot ($>$2000 K) exoplanets, but space-based
investigation into limb asymmetry is in its infancy, and limb asymmetry’s
prevalence in the broader exoplanet population remains unexplored. We find
evidence for limb asymmetry on the exoplanet WASP-107b via transmission
spectroscopy from 2.5 to 4.0 micrometers with JWST/NIRCam. This is one of the
first low-resolution space-based measurements of limb asymmetry and is unique
because, at 770 K, WASP-107b is in a relatively cool regime where planetary
terminators are expected to be homogeneous. These observations imply a
difference in temperature and cloud properties between WASP-107b’s limbs,
challenging our models of limb asymmetry in this cooler regime.

The WASP-107 system consists of an active K6 star [1, 2, 3] orbited by at
least two planets. First discovered in 2017 [4], the inner planet WASP-107b is
among the lowest bulk-density exoplanets known to date. With a mass comparable
to Neptune (Mp = 30.5 M⊕ [5]) but a radius comparable to Jupiter (Rp = 11.0 R⊕
[4]), its bulk density of 0.126 g/cm3 is roughly one-fifth that of Saturn, the
lowest density planet in our Solar System. WASP-107b’s density suggests it is
a gaseous planet with a substantial atmosphere of primarily hydrogen and
helium [4, 5]. Further, such low density implies that WASP-107b has a uniquely
large atmospheric scale height, which enhances the strength of molecular
absorption features and makes it a favorable target for transmission
spectroscopy.

A transmission spectrum is measured from starlight that passes through the
planet’s atmosphere at its terminators, which are the dividing lines between
the planet’s day- and nightside. In the transiting geometry, the planetary
limbs encompass its terminators. Tidally synchronized planets have perpetual
day- and nightsides, leading to day-night irradiation contrast and global
atmospheric circulation. Limb heterogeneity develops as this circulation
advects heat in one direction, driving a temperature difference between the
evening terminator, which is downstream of the heat flow, and the morning
terminator, which is upstream [6, 7, 8]. This heterogeneity is generally
expected to be strongest on tidally synchronized planets with temperatures
between $\sim$1200 - 2100 K, with cooler planets having homogeneous limbs [6,
7, 9].

Several models have been used to simulate the signature of limb asymmetry in
transmission [7, 10, 9]. In short, a difference in limb temperatures imparted
by atmospheric circulation primarily causes a difference in scale height (via
temperature) and cloud properties, which together alter the wavelength-
dependent transmission through each limb. For example, with clouds, this
temperature difference can be sufficient that species are gaseous in the
evening limb, but have condensed in the morning limb. If neglected, limb
asymmetry can lead to biased retrievals of atmospheric properties from low-
resolution observations, including JWST [7, 11, 12, 13, c.f. [14]]. Therefore,
observationally characterizing limb asymmetry is vital to ensuring that our
insights to exoplanet atmospheres are bias-free and benchmarking models of
atmospheric circulation and cloud formation. As mentioned, atmospheric models
typically predict limb properties to homogenize below $\sim$1200 K, but this
has not been observationally confirmed. At an equilibrium temperature of 770
K, WASP-107b lies in this cooler regime, and its large atmospheric scale
height makes it a uniquely favorable target for investigating limb asymmetry
at these temperatures.

##  1 Results

###  1.1 Observations

We observed one transit of WASP-107b on January 14th, 2023 using JWST/NIRCam
as part of the MANATEE NIRCam GTO program (Obs.8, JWST-GTO-1185, ref.[15]). We
simultaneously used NIRCam’s F210M filter in the short-wavelength channel and
the F322W2 grism in the long-wavelength channel [16, 17]. For both, this
observation was a continuous time-series of 1,293 integrations at 20.22-second
cadence using the BRIGHT2 readout pattern on the SUBGRISM256 subarray. We
reduced these data using three independent pipelines: Eureka! [18], tshirt
[19], and Pegasus [Beatty, T., et al. in review]. Each pipeline produced
consistent results, discussed in the Methods, validating our results against
systematic biases from the particular reduction method used. We focus
primarily on the results from tshirt due to its more effective 1/f noise
correction [20]. We produced a band-integrated light curve at 2.093 $\pm$
0.102 $\mu$m from F210M, as well as a band-integrated (2.45 – 3.95 $\mu$m)
light curve and thirty (R$\sim$64) spectroscopic light curves from F322W2. Our
F322W2 dynamic spectrum and broad-band light curve are shown in Figure 1.

Figure 1: Dynamic spectrum and broad-band light curve from our JWST/NIRCam
F322W2 transit observation of WASP-107b. The top panel shows our spectroscopic
transit data at full spectral and temporal resolution, with wavelength on the
y-axis and time on the x-axis in terms of the sequence of integrations. Each
column has been normalized by its median value, so the color-coding
illustrates the corresponding relative flux values. In our analysis, we only
use dispersion columns 56 - 1594, or wavelengths from 2.45 - 3.95 $\mu$m, in
order to avoid excess noise visible on either edge of the detector. The bottom
panel shows our broad-band light curve integrated over those wavelengths, with
no additional light curve detrending done. Overall, this light curve displays
no strong non-linear systematic trends or starspot crossing events.

We also collected new and archival transit and radial velocity observations of
WASP-107b using other instruments over a wide wavelength range. Uncertainty in
the planet’s orbital parameters, particularly the time of conjunction, can
bias the signal of limb asymmetry [7, 21, 22]. To combat this, we use these
ancillary data to precisely measure WASP-107b’s orbit before our analysis of
the JWST/NIRCam data. These include radial velocity measurements with the
CORALIE and Keck I/HIRES spectrographs collected by refs.[4, 5]. Also, we
downloaded public data covering three transits with the Transiting Exoplanet
Survey Satellite, and one transit at 4.5 $\mu$m with the Infrared Array Camera
on the Spitzer Space Telescope (Program 13052, PI: M. Werner). We observed one
new transit in the SDSS i-band using the Goodman Spectrograph [23] on the
Southern Astrophysical Research Telescope. Lastly, we obtained the broad-band
(5 – 12 $\mu$m) light curve of one transit using JWST’s Mid-Infrared
Instrument in Low-Resolution Spectroscopy mode (MIRI/LRS), observed in JWST
program 1280 and given to us by PI P.O. Lagage. We describe the reduction and
fitting of these data in the Methods, with particular focus on our derivation
of WASP-107b’s time of conjunction in Section 3.5. The results are tabulated
in Table 1. With these data, we measured WASP-107b’s time of conjunction to a
1-$\sigma$ precision of 0.7 seconds.

###  1.2 Evidence for limb asymmetry

####  1.2.1 Limb transmission spectra

We fit the spectroscopic light curves with asymmetric-limb transit models
using the modeling package catwoman [22, 24]. Described further in Methods
Section 3.4.4, we fit for the planet’s evening and morning limb planet-star
radius ratios at each wavelength, and derived transmission spectra for each
limb. We assume the planet’s rotation axis is aligned to its orbit, discussed
further in Methods Section 3.4.5. We left the time of conjunction ($t_{c}$) as
a free parameter with a Bayesian prior set to its posterior distribution from
fitting our ancillary observations, and $t_{c}$ did not vary from this prior.
The best-fit spectra are shown in Figure 2 and tabulated in Table 2.

Figure 2:  The separate transmission spectra of WASP-107b’s morning and
evening limbs. The top panel shows the transmission spectra of the evening
(orange circles) and morning (blue circles) limbs, derived from our
JWST/NIRCam F322W2 observations. We find the evening and morning spectra
differ by $\sim$250 ppm, with the evening spectrum at larger transit depth
(i.e., that limb is physically larger) at most wavelengths. The strongest
asymmetry appears near 2.7 $\mu$m, where the evening spectrum displays a
strong apparent feature. In contrast, the morning spectrum appears relatively
flat and lacks strong features. We compare these to atmospheric models based
on grid fits to these spectra (Figure 4). The solid gray lines show our best-
fit solution where the evening limb is hotter by $\sim$180 K, and the limbs
have different gray cloud opacity. To demonstrate the individual effects of
temperature and cloud opacity, we show two additional model pairs that are not
best-fit solutions, but exhibit the effect of changing just one parameter in
an extreme way. These are the dashed lines, which share the same cloud opacity
(10-1.58 cm2 g-1) but differ in temperature by $\sim$190 K; and the dotted
lines, which share the same temperature (713 K) but differ in log cloud
opacity by 0.8 cm2 g-1. The lower panels show the difference between the
morning and evening depths in units of their mutual uncertainty, and radii in
units of WASP-107b’s atmospheric scale height. The limb depths are
$>$1-$\sigma$ different at most wavelengths, and over 3-$\sigma$ different
near 2.7 $\mu$m and 3.3 $\mu$m. WASP-107b’s limb radii differ by 0.5 - 2 scale
heights at these wavelengths.

We find that WASP-107b’s morning and evening limbs have significantly ($>$99%
confidence) different transmission spectra, as seen in Figure 2. Across our
bandpass, WASP-107b’s evening limb is generally larger than its morning limb.
The evening limb spectrum is $\sim$250 ppm higher in their pseudo-continuum
regions where the limb depths differ by $\sim$2-$\sigma$, and the
corresponding radii differ by roughly one scale height. The evening spectrum
displays strong features near 2.7 $\mu$m and 3.3 $\mu$m where the asymmetry is
larger, the individual limb depths differ by over 3-$\sigma$, and the
corresponding radii differ by up to two scale heights. In contrast, the
morning limb spectrum appears relatively flat and lacks strong features.

As discussed by ref.[9], limb asymmetry can be more confidently attributed to
an atmospheric signal, rather than systematics or ephemeris error, if larger
asymmetry is observed at wavelengths corresponding to expected features in the
transmission spectrum. The larger asymmetry we see around 2.7 $\mu$m matches
an expected water absorption band, which was detected in WASP-107b’s
atmosphere by previous HST (6.5$\sigma$ significance) [25] and JWST/MIRI
observations ($>$12.5$\sigma$ significance) [26]. The presence of other
molecules in WASP-107b’s atmosphere is being explored in a forthcoming paper
[Welbanks, L., et al. in review]. There is no reason to suggest water would
not be present at both terminators, so the relatively flatter morning limb
spectrum is likely the result of that limb being cooler which reduces the size
of transmission features. As discussed further later, cloudiness likely plays
a role too.

We ran several statistical tests to verify that the limb asymmetry we observe
is real. First, we compared our transmission spectra to a “null” case of a
model atmosphere with uniform-limbs using a $\chi^{2}$ test, which gave
$\chi^{2}/\mathrm{dof}$ = 5.37. This rejects the uniform-limb model at high
($>$99%) confidence. We also perform a paired sample t-test on the spectra,
which returned a p-value of 4 $\times~{}10^{-11}$, and also strongly rejects
our morning and evening spectra having identical means. Both tests assume the
measured depths are uncorrelated across spectral channels, which we believe is
true based on several further tests (Methods Section 3.1). We also verified
that our asymmetric-limb light curve models fit the data better than a
uniform-limb model with the same $t_{c}$ and total transit depth. Our
asymmetric-limb models were always preferred by the reduced $\chi^{2}$, the
Bayesian Information Criterion (BIC) by up to $\Delta$BIC = 10, and the Akaike
Information Criterion (AIC) by up to $\Delta$AIC = 8.1. The Bayes Factor also
preferred the asymmetric-limb model by factors of 1.1 - 137.6, corresponding
to detection significances of 1.15 - 3.62$\sigma$ [27, 28].

Next, we verified that our fits were not artificially injecting limb
asymmetry. We generated thirty synthetic light curves of a uniform-limb planet
with the same cadence, transit depth, and wavelength-dependent light curve
scatter as our real observations. We drew the scatter from a Gaussian
distribution, as our real data residuals are Gaussian based on an Anderson-
Darling test. We fit these simulated observations with catwoman and repeated
for ten realizations, extracting the limb spectra each time (Extended Data
Figure 5). In almost all cases we recovered equal limb depths, and outliers
with differing limb depths were random in wavelength and polarity, unlike our
real observations.

Finally, we checked for the signal of limb asymmetry directly in the data.
Figure 3 shows data at 2.68 $\mu$m, where we see the largest limb asymmetry.
We folded the light curve about the time of conjunction to compare the data in
ingress and egress. If WASP-107b has asymmetric limbs, we would expect to see
differences in the ingress and egress times caused by the larger or small
sides of the planet reaching each contact point earlier or later than for a
symmetric planet. Our observations indicate that the evening limb on WASP-107b
is larger, so we would expect transit ingress and egress to occur late
relative to a uniform-limb transit. In the folded data shown in Figure 3, this
would manifest as a positive difference between the observations taken during
ingress and egress. We see this predicted signal directly in our data (Figure
3 Panel A). Subtracting the data after transit center from that before transit
center, we see a large and consistently positive difference during the
ingress/egress period.

The data also show significant residuals compared to a uniform-limb transit
model during ingress and egress, as predicted by refs.[7, 9]. Panel B of
Figure 3 compares the residuals of our best-fit asymmetric-limb model to a
uniform-limb model with the same transit depth. The residuals from the two
fits match during mid-transit (between the second and third contact points)
and the residuals are consistent with each other out-of-transit. However, the
residuals of the asymmetric- and uniform-limb models do not match during
ingress and egress, where we see the residuals to the uniform-limb model have
a significantly non-zero mean of $198\pm 61$ ppm, while the mean of
asymmetric-limb model’s residuals is $30\pm 61$ ppm. This indicates that the
uniform-limb model does not well fit the observations, and increases our
confidence that that the observed ingress and egress of WASP-107b are indeed
asymmetric.

Figure 3: Signal of limb asymmetry in our JWST/NIRCam data. Here, we focus on
our 2.68 $\mu$m channel which has the largest evening-morning difference
(Figure 2). Panel A shows observed data, while Panel B shows residuals from
light curve models with and without limb asymmetry. We fold each x-axis about
the time of conjunction ($t_{c}$) so that ingress and egress are overlaid and
more easily compared. We zoom in around the ingress/egress. In Panel A, blue
points show the data before $t_{c}$ while gray points show the data after
$t_{c}$. We see a consistently positive difference between the pre- and
post-$t_{c}$ points reaching $\sim$1000 ppm within ingress and egress, but no
difference outside of ingress and egress. This is caused by the larger evening
limb effectively shifting egress further from $t_{c}$ relative to ingress, so
that ingress points at a given distance from $t_{c}$ are at higher relative
flux than the corresponding egress points. If the limbs were uniform, then
ingress and egress would be symmetric about $t_{c}$. In Panel B, gray points
show the residuals of our data with the best-fit catwoman asymmetric-limb
model at this wavelength (from which Figure 2 is derived). The red points are
residuals with a uniform-limb model. The models are identical except within
ingress and egress, where the uniform-limb model residuals are much larger,
reaching a mean values of 198 ppm versus only 30 ppm for the asymmetric-limb
model. The difference between the models (black line) accordingly reaches a
peak of 253 ppm midway through ingress/egress. This large residual shows that
the data is not well fit by a model without limb asymmetry.

####  1.2.2 Ruling out alternate explanations

We can rule out several possible non-atmospheric origins for the limb
asymmetry we infer. First, we investigated the impact of the measurement
uncertainties in WASP-107b’s orbital parameters and the star’s limb-darkening
coefficients. We repeated our light curve fits allowing these to freely vary
in each channel, using Gaussian-shaped priors set by our joint-fit results and
an ATLAS stellar model (Methods Section 3.4.1), as well as fully free limb-
darkening coefficients. This did not change the resulting limb spectra and
only increased the individual depth uncertainties by 1% on average.

Next, we checked the impact of transit timing variations (TTVs). Previous
observations have ruled out TTVs greater than 20 seconds [1], but smaller TTVs
must be considered given our sensitivity to timing precision. Using TTVFaster
[29], we estimate that the peak-to-peak TTV amplitude of WASP-107b due to its
outer companion WASP-107c [5] is $<$0.4 seconds. This is smaller than our
uncertainty on WASP-107b’s $t_{c}$, and our spectra are robust to varying
$t_{c}$ by this amount.

We also considered the timing accuracy of our instruments. Our strongest
constraint on $t_{c}$ comes from the JWST observations. We consulted with JWST
project scientists and verified that the image time-stamps (INT_TIMES) are
accurate within $<$0.5 seconds, and done consistently between instruments. We
also verified that the onboard barycentric correction of these time-stamps
fully accounts for the spacecraft position and velocity. For all JWST
instruments, the INT_TIMES values correspond to the time that the last pixel
along the read direction is read-out. In our case, this occurs 1.347 s after
the first pixel is read-out. In this work, we always used the mid-integration
time, as given in the INT_TIMES table, corrected to the middle of array read-
out and converted to BJD TDB format.

We also checked for excess correlated noise in our JWST/NIRCam data. The
residuals between the data and the best-fit asymmetric-limb models were
Gaussian distributed and passed Anderson-Darling and Shapiro-Wilks tests. From
an Allan variance test, however, the root-mean-square of the residuals did
deviate slightly from pure-white noise behavior, falling off as $N^{-0.4}$
with the number of integrations per bin ($N$) rather than $N^{-0.5}$. We also
see minor fluctuations in the light curves of the detector reference pixels
throughout the observation. Excess red noise may exist in our data, but we
consider it very unlikely that red noise introduced a false limb asymmetry
signal based on the injection-recovery tests we performed (Extended Data
Figure 5).

Lastly, we consider the transit light source effect (TLSE). Stellar surface
heterogeneities can affect the observed transmission spectrum [30, 31], and
spot crossings have been previously reported at optical wavelengths [4, 2, 1,
25] in transit observations of WASP-107b. However, we do not see evidence for
a starspot crossing in our data. Based on TLSE models for a K6 star such as
WASP-107[31], unocculted starspots would change our measured transit depths by
$\sim$50 ppm – much smaller than the observed morning-to-evening differences.

We also considered if a small starspot located right on the stellar limb could
be causing the observed limb-asymmetry signal. If such a spot crossing is the
sole source of the limb asymmetry, it would have to both offset the limb
spectra and change their relative shapes. Due to limb-darkening, starspot
crossings on the stellar limb have an order of magnitude weaker effect in the
light curve than when the spot is occulted mid-transit, but may induce large
($\sim$10-100 s) TTVs [32]. If a limb-spot crossing were affecting our data,
it would also bias the $t_{c}$ we derive from the simultaneous F210M
observation (relative to the ancillary observations), but this is not the case
(Methods Section 3.5).

##  2 Discussion

To date, limb asymmetry has been detected via ground-based high-resolution
transmission spectroscopy on the T${}_{eq}=$ 2150 K WASP-76b [33, 34, 35, 36],
the 2260 K MASCARA-2b [37], and the 2720 K WASP-121b [38, 39] via molecular
abundance and wind speed gradients. WASP-107b joins the 1200 K WASP-39b
[Espinoza, N., et al. in review, Delisle, S., et al. in review] as the only
planets on which limb asymmetry has been detected via space-based low-
resolution transmission spectroscopy, and separate limb transmission spectra
recovered. Our discovery on WASP-107b is uniquely surprising because, at an
equilibrium temperature of only 770 K, WASP-107b is significantly cooler than
these other planets. Such cool atmospheres are not expected to exhibit strong
spatial heterogeneity, as their circulation is expected to efficiently
distribute heat globally, leading to generally low day-night contrast and
homogeneous terminators [e.g. 7, 6, 9]. Our result suggests that either limb
asymmetry is more predominant than current models predict, or WASP-107b is an
outlier.

The comprehensive modeling and the interpretation of our observational results
is beyond the scope of this work, since it will likely require a unified
treatment of 2D effects and clouds in WASP-107b’s atmosphere. Nevertheless, we
did compare our results to separate 1D atmosphere models for the morning and
evening limbs to provide an initial estimate of the atmospheric parameters.

We compared our observations to atmospheric models for the morning and evening
terminator generated using ScCHIMERA. Described further in Methods Section
3.6, we fixed each model to a metallicity of $\sim$10 times the Solar value
and carbon-to-oxygen ratio of 0.35. These are the bestfit values from a
combined retrieval analysis of previous HST/WFC3 and JWST/MIRI observations
[26], and are very consistent with the additional retrievals of Welbanks, et
al. (in review). We include a vertically uniform gray “cloud”. The clouds or
hazes in WASP-107b’s atmosphere are likely more complex than this, but this
simple treatment enables initially probing the relative role of clouds/hazes
without adding more free parameters than is justified. We fit a grid of
temperature and cloud opacity to the evening and morning limb spectra
individually using the $\chi^{2}$-minimization method [40, 41].

We find that WASP-107b’s evening and morning terminators must have
significantly different atmospheric properties, shown in Figure 4. There is
some correlation between the effects of temperature and cloud opacity, but the
evening and morning distributions do not overlap even far beyond their 99%
confidence intervals. The data prefer a cooler morning limb and a hotter
evening limb, with a slight difference in the cloud opacity. The “best-fit”
models, corresponding to the minimum $\chi^{2}$ values, are marked by the
points in Figure 4 and compared to the data in Figure 2. Based on these best-
fit values and 68% confidence intervals for our tshirt reduction, we infer an
evening-morning temperature difference of 180 $\pm$ 36 K. The precise
difference is slightly reduction dependent with the PEGASUS and Eureka!
spectra suggesting differences of 68 $\pm$ 29 K and 137 $\pm$ 20 K,
respectively, but a hotter evening limb is always preferred. More
spectroscopic data will help refine this measurement, explore the difference
in cloudiness in more detail, and investigate whether limb asymmetry affects
other molecules such as photochemically produced SO2 [26, 42].

The evening-morning temperature differences we derive are $\sim$10-23% of
WASP-107b’s equilibrium temperature. Adapting a common metric for quantifying
day-night heat redistribution on hot Jupiters (e.g.,[43]), we compute Ae-m
$\approx$ 0.33-0.64 as the proportional evening-morning energy difference,
where Ae-m $\equiv\left(F_{evening}-F_{morning}\right)/F_{evening}$ and
$F_{i}=\sigma T_{i}^{4}$. Unfortunately, small sample size and differences in
analyses make it difficult to compare WASP-107b in detail to the other planets
on which limb asymmetry has been detected. WASP-76b is the only other planet
with published evening-morning temperature differences; refs.[35, 36] present
large, varied estimates ranging 16 - 35% of WASP-76b’s equilibrium
temperature, or Ae-m $\approx$ 0.38-0.83. Using the day- and nightside
temperatures of WASP-121b from ref.[38] as proxies for its evening and morning
temperatures, we can infer a difference of $>$25% its equilibrium temperature,
and Ae-m $\approx$ 0.65. Therefore, WASP-107b’s limb asymmetry may be
proportionally comparable to those of the much hotter WASP-76b and WASP-121b,
though more data is needed before any population-level trends can be drawn.

To illustrate the individual effects of temperature and cloudiness on the limb
spectra, we show two additional model pairs in Figure 2. The dashed models
share an opacity of 10-1.58 cm2 g-1, midway between the distributions in
Figure 4, and differ in temperature by $\sim$190 K (611 K and 798 K). The
dotted models share a temperature of 713 K, again between the distributions,
and differ in log cloud opacity by 0.8 cm2 g-1. These pairs roughly illustrate
the limits of the distributions in Figure 4 along each axis. Varying
temperature has the strongest effect on the limb spectra, recreating both the
offset and relative amplitude between morning and evening. Varying cloud
opacity can recreate the offset but does not as well explain the difference in
spectral amplitudes. In reality, both temperature and cloud/haze properties
may change spatially and feedback on one another, highlighting the complex
multidimensionality of exoplanet atmospheres. In this work, we demonstrate
that unraveling this multidimensionality, and probing morning-to-evening
differences is possible with JWST. This new capability opens an exciting new
path toward a complete picture of these worlds beyond our own.

Figure 4: Atmospheric model fits to WASP-107b’s evening and morning limb
spectra. Shown are the results of our grid-fits described in Section 2 and 3.6
for the results of all three data reductions tried in this work. We primarily
focused on the results from tshirt (bottom panel) due to its more robust
correction for 1/f noise, but the spread of parameter values between
reductions is instructive for interpreting these values. We fixed the
atmospheric metallicity to $\sim$10x the solar value and carbon-to-oxygen
ratio to 0.35 based on the results of ref.[26] and Welbanks, et al. (in
review), then fit for the redistribution parameter (i.e., terminator
temperature, related as $T=738.31\times\text{redistribution}^{0.25}$) and gray
cloud opacity for each limb individually. Results for the evening and morning
limb are shown by the orange and blue contours, respectively, with contours
drawn at the 68% and 99% confidence intervals. We find that the atmospheric
conditions on each limb are significantly different. There is some correlation
between the temperature and cloud opacity, but this difference in atmospheric
conditions is robust. WASP-107b’s morning limb is generally cooler than its
evening limb, with a best-fit difference in temperature of 180 $\pm$ 36 K from
tshirt, 68 $\pm$ 29 K from PEGASUS, and 137 $\pm$ 20 K from Eureka!. These
fits suggest a difference in cloud properties as well, which should be
explored in more detail with additional observations.

##  3 Methods

###  3.1 JWST/NIRCam Observation

Our transit observation of WASP-107b using JWST/NIRCam was Observation 8 of
program JWST-GTO-1185, which executed on January 14th, 2023. We acquired on
our target WASP-107 using one 0.152 second integration with JWST/NIRCam
imaging in its F335M filter, with RAPID readout on the SUB32TATSGRISM
subarray. After target acquisition, we began our science observation
consisting of simultaneous short-wavelength channel imaging with filter
WLP8/F210M and long-wavelength channel spectroscopy with GRISMR/F322W2 [16,
17]. We used the BRIGHT2 readout pattern on subarray SUBGRISM256, and
collected 1,293 integrations at 20.22 second cadence. Each integration
consisted of seven groups. Our total exposure time was 26,145.53 seconds.

Our reductions of these data using three different pipelines are described in
the following section. From each reduction, we extracted thirty spectroscopic
lightcurves from our F322W2 data in equally-spaced wavelength bins from 2.45
to 3.95 $\mu$m, corresponding to a spectral resolution of approximately 64. We
chose to use this resolution because it offers a good trade-off between
spectral information and transit depth measurement precision, as well as
protects against optical correlations. At this low spectral resolution, we are
binning together $\sim$28 detector resolution elements per channel. As a
result, the point spread function is much smaller than our bin sizing and
there is no optical correlation between channels. We further mitigate possible
bin-to-bin correlations by correcting for 1/f noise, as described in the
following Section.

###  3.2 JWST/NIRCam Data Reduction

####  3.2.1 Reduction with tshirt

We used the tshirt pipeline (https://tshirt.readthedocs.io) to extract
spectroscopic lightcurves of WASP-107b’s transit. As in ref.[44], we modified
the CALDETECTOR1 stage of the jwst pipeline to turn the _uncal files into rate
images with less 1/$f$ noise. We started with JWST pipeline version 1.8.4,
CRDS Version 11.16.16, and CRDS context jwst_1039.pmap. tshirt replaces the
default reference pixel correction with a row-by-row, odd/even by amplifier
(ROEBA) correction [45] using the background pixels, which can reduce the
1/$f$ noise as compared to reference pixels alone [20]. For the odd/even slow-
read direction, we used the reference pixels in the bottom 4 rows. For the
row-by-row fast-read correction, we used all pixels from X=1846 to 2043
(0-based) to estimate the sky background and 1/f noise in that row. For the
jump step of the pipeline, we used a threshold of 6 $\sigma$.

We then continued the rest of the steps in CALDETECTOR1 with the default
parameters. After constructing rate files for each integration, we manually
divided the images by a jwst_nircam_flat_0266.fits imaging flat field and
marked all pixels that have a “DO_NOT_USE” DQ value to NaN. We multiplied all
count rate images by the gain and integration time to estimate the total
number of electrons at the end of the ramp for photon error calculation.
Finally, we cleaned cosmic rays from all images by combining 150 rate images
at a time and marking all pixels that deviate from the median rate image by
more than 20 times the jwst pipeline error for a given rate image for
replacement. The bad pixels are replaced by the linearly interpolated value of
the 150-pixel time series.

We performed a column-by-column background subtraction with a linear fit along
the Y direction to pixels Y=5 to 24 and Y=44 to 65 for all rate images. For
the spectral extraction, we first fit the spectrum row-by-row along the
dispersion direction of integration 647 with a smooth 40-knot univariate
spline with SciPy in order to build a profile function as a function of
wavelength. We used a co-variance weighted extraction [20], assuming a uniform
correlation between pixels of 0.08 and a read noise of 14 $e^{-}$. We
extracted a rectangular aperture from 29 to 39 pixels in the Y direction and 4
to 1747 in the X direction. While this rectangular aperture does not track the
curvature of the trace and could influence the absolute flux and aperture
corrections, we normalized the spectra by the out-of-transit baseline, so this
effect is not expected to affect the derived planet spectra. We used the
profile along the Y direction to estimate the missing flux in pixels that
marked as NaN and which deviate from the average profile by $\geq$30$\sigma$.

We fit the resulting lightcurves produced by tshirt in the same method as for
all reductions, which is described in Section 3.4.4. The best-fit transmission
spectra are compared to those of our other reductions in Figure 6 of the
Extended Data.

####  3.2.2 Reduction with Eureka!

Our reduction of the NIRCam F322W2 data with Eureka! followed the same
procedure as described by ref.[44], with the only differences being package
version numbers. In particular, we used version 0.10.dev15+g80126b56.d20230613
of Eureka![18], CRDS version 11.17.0 and context 1093, and jwst package
version 1.10.2. The exact Eureka! Control Files and Eureka! Parameter Files we
used are available for download (Zenodo link to be created on final
submission). In summary, we cropped the frame to focus on relevant pixels,
corrected for the curvature of the spectral trace, performed background
subtraction per-column and per-integration, and then performed optimal
spectral extraction using the pixels within 9 pixels of the flattened spectral
trace. We then sigma-clipped the observations compared to a boxcar smoothed
version of the observations to remove cosmic rays without removing sharp
astrophysical features like the transit ingress/egress.

We fit the resulting lightcurves produced by Eureka! in the same method as for
all reductions, which is described in Section 3.4.4. The best-fit transmission
spectra are compared to those of our other reductions in Figure 6 of the
Extended Data.

####  3.2.3 Reduction with Pegasus

We also reduced and extracted the transmission spectrum of WASP-107b using the
Pegasus pipeline (https://github.com/TGBeatty/PegasusProject). This pipeline
will be described in more detail in a forthcoming paper [Beatty, T., et al. in
review], and we briefly summarize it here. We began with the rateint files
from the jwst pipeline v1.10.2, using CRDS version 11.17.0 and context 1093.

For each rateint file we first performed a background subtraction step by
fitting a two-dimensional second-order spline to each integration using the
entire 256$\times$2048 rateint images. We masked out image rows 5 to 75 to not
self-subtract light from the WASP-107 system, and then we fit individual
background splines to each of the four amplifier regions in the images. We
then extrapolated the combined background spline over the masked portions near
the star to perform the background subtraction. Visual inspection of the
rateint images also showed that in roughly 5% of the integrations the
reference pixel correction failed for at least one of the amplifier regions,
so after the spline fitting and subtraction we re-ran the reference pixel
correction using hxrg-ref-pixel (https://github.com/JarronL/hxrg_ref_pixels).
This appeared to correct the issue.

We then performed a spectrophotometric extraction on the background-subtracted
images using optimal extraction techniques [46]. We iteratively constructed a
smoothed spatial profile for the F322W2 data using pixel columns from column
22 to 38 (inclusive) and pixel rows from row 5 to 1650 (inclusive), and we
iteratively estimated a variance matrix for this same region starting from the
pixel uncertainties and read-noise values following established techniques
[46]. We then fit a 4th-order polynomial spectral trace to each integration
and extracted fluxes in each wavelength bin using an extraction aperture with
a half-height of seven pixels centered on the estimated trace in each detector
column. In doing so we accounted for partial pixel effects in both the
spectral and spatial directions.

We fit the resulting lightcurves produced by Pegasus in the same method as for
all reductions, which is described in Section 3.4.4. The best-fit transmission
spectra are compared to those of our other reductions in Figure 6 of the
Extended Data.

####  3.2.4 Comparing Reduction Pipelines

We reduced our JWST/NIRCam F322W2 data with three independent reduction
pipelines described above, but ultimately present results from the tshirt
reduction in our main analysis. Here, we compare the results from each
reduction to show that while this choice is justified, it does not make a
difference to our overall result.

As in our primary analysis, we fit the spectroscopic light curves from each
reduction with an asymmetric-limb transit model (see Section 3.4.4 for
details) and derived morning and evening limb transmission spectra for each.
All three reductions lead to the same qualitative result: an evening-limb
spectrum that is generally at higher transit depth and displays stronger
spectral features compared to the relatively flat morning-limb spectrum
(Extended Data Figure 6). Comparing the transit depths in each individual
channel, the morning and evening limb depths are all consistent within
1$\sigma$ between reductions in each individual channel, except for one
outlier at $\sim$3.77 $\mu$m where Eureka! and Pegasus differ by
$\sim$1.08-$\sigma$ (Extended Data Figure 6). There are no channels where
tshirt disagrees with the other two reductions by more than 1-$\sigma$. The
agreement between the evening limb spectra are particularly good, where the
depths differ by less than 0.5 mutual standard deviations in most channels.
The only notable difference between the three reductions’ results is that the
morning limb spectrum derived from tshirt is generally slightly lower than
that of PEGASUS and Eureka! at short wavelengths, though never by more than
1-$\sigma$ as mentioned. A potential reason for this, as well as our
justification of choosing tshirt as our primary focus, is described further
below.

In the main text (Section 1.2), we presented several statistical tests to show
that the morning and evening limb transmission spectra from the tshirt are
significantly different. We repeated these tests for the other two reductions
as well. Between the morning and evening spectra, we compute
$\chi^{2}/\mathrm{dof}$ values of 2.80 and 4.39 for the PEGASUS and Eureka!
reductions, respectively. In both cases, for thirty degrees of freedom, the
value of the survival function for the “null” hypothesis that the morning and
evening spectra are identical is effectively zero. Similarly, we compute
paired t-test p-values of 2 $\times~{}10^{-6}$ and 1.3 $\times~{}10^{-7}$ for
PEGASUS and Eureka!, respectively. In both cases, we again reject a case where
the morning and evening limb spectra have identical means. We also repeated
the grid-fit of model atmosphere spectra over temperature and gray cloud
opacity, described in the Main text and Section 3.6, for the PEGASUS and
Eureka! spectra. The resulting $\chi^{2}$ distributions (e.g., as in Figure 4)
overlapped, and tell a consistent story that the atmospheric properties on
each limb are different. Regardless of which reduction is used, our result
that WASP-107b’s evening and morning limb transmission spectra are
significantly different in the JWST/NIRCam F322W2 bandpass, and indicate a
difference in temperature and cloudiness between terminators, holds true.

We ultimately chose to use tshirt in our primary analysis because, of the
three pipelines, t-shirt performs the most robust correction for 1/f noise.
This 1/f noise is one of, if not the most, dominant systematic noise source in
NIRCam time series data using GRISMR readout [20, 45]. Eureka! and PEGASUS
both adopt very similar techniques for 1/f noise correction, including a
“reference pixel step” where the mean value of each detector amplifier’s
reference (i.e., unexposed) pixels are subtracted from the value of the
exposed pixels. This step is meant to eliminate overall bias offsets which may
be different between amplifiers, as well as odd/even column effects, pre-
amplifier resets, and 1/f noise. However, these reference pixels exhibit noise
and do not sample frequencies greater than $\sim$200 Hz, reducing the
effectiveness of this correction when used alone. Further, this step can
occasionally fail, leaving systematic biases including 1/f noise present in
the data. In these failure cases, PEGASUS then uses a custom routine hxrg-ref-
pixel (linked in Section 3.2.3) to re-correct the data, while Eureka! does no
re-correction. On the other hand, t-shirt replaces this reference pixel step
with row-by-row, odd/even by amplifier (ROEBA) correction, as described above
and in ref.[45]. ROEBA correction better corrects for and reduces 1/f noise in
the data than the default correction [20, 45]. This difference in 1/f noise
correction is most likely the reason for the slight differences between the
spectra derived from tshirt and the other two pipelines, particularly at the
shortest wavelengths of the morning limb spectra (Figure 6). A similar slight
offset in absolute transit depth between tshirt and Eureka! has been seen in
other analyses, and was directly traceable to this 1/f noise correction [44].

As a final comparison, we repeat the grid-fit of terminator temperature (via
the redistribution parameter) and gray cloud opacity, as done for tshirt in
the main text, on the spectra derived from PEGASUS and Eureka!. The results
are shown in Figure 4. We find that, while the specific solution is slightly
different between reductions, the general result remains the same. From all
three reductions, we recover the result that WASP-107b’s evening limb is
hotter than its morning limb, and that there may be a difference in cloudiness
or haziness as well. In all three cases, the evening and morning distributions
do not overlap even well beyond their respective 99% confidence intervals.
Even marginalizing over the cloud opacity, the evening-morning temperature
difference is significant at the respective 68% confidence intervals for all
three reductions. We find that tshirt gives the largest evening-morning
temperature difference of 180 $\pm$ 36 K, PEGASUS gives the smallest
difference of 68 $\pm$ 29 K, and Eureka! is in between with 137 $\pm$ 20 K. Of
these, the difference derived by tshirt has the largest uncertainty.

###  3.3 Ancillary Data Reduction

As mentioned in Section 1, we used additional observations of the WASP-107
system to support the analysis of our JWST/NIRCam observation. The radial
velocity data we used included 31 measurements collected by ref.[4] between
2011 and 2014 using the CORALIE spectrograph on the 1.2-meter Euler-Swiss
telescope, in addition to 60 measurements collected by the California Planet
Search (PI: A. Howard) for ref.[5] between 2017 and 2020 using the High
Resolution Echelle Spectrometer [47] on the Keck I telescope. A full table of
these radial velocity measurements, including uncertainties, is included in
the electronic journal version of ref.[5]. We recorded these tabulated values,
assuming all necessary data reduction and detrending was already performed,
and used them in our analysis.

TESS observed the WASP-107 system in its Sector 10, and captured three full
transits and one partial transit of WASP-107b at a thirty-minute cadence
between March and April 2019. We downloaded the publicly available light
curves of these observations from Mikulski Archive for Space Telescopes, using
the “light curve” data files (file extension .lc) which have known systematics
already removed from the data. The out-of-transit flux values exhibited no
strong outstanding systematic trends. We discarded the partial transit and
only used data around the three full transits. In order to improve
computational efficiency, we only used data within $\pm$0.45 days of each
transit midpoint which is sufficient out-of-transit baseline.

We also observed one transit of WASP-107b using the Goodman High Throughput
Spectrograph [23], in imaging mode with an SDSS-i filter, on SOAR during
program N23A-840705 (PI: M. Murphy). We reduced these observations using the
AstroImageJ pipeline [48], which included comparison star detrending for basic
telluric variation corrections. Despite this initial detrending, the light
curve still displayed a significant, non-linear systematic trend so we first
performed additional airmass detrending. The airmass of WASP-107 varied
quadratically throughout our observation, reaching a minimum near the middle
of our night. As a result, the flux of WASP-107 relative to the comparison
stars was an asymmetric and non-injective function of airmass. That is, the
relative flux at a given airmass took different values depending on whether
the time was before or after the airmass minimum, likely the result of the
telescope structure needing to flip orientation as well as the sky properties
changing throughout the night. We fit two separate quadratic functions of the
form

$f=p_{2}x^{2}+p_{1}x+p_{0}$ (1)

to the data before and after the airmass minimum. Here, $f$ is the relative
flux, $p_{0}$ \- $p_{2}$ are fitting coefficients, $x\equiv z-z_{median}$ for
the corresponding subset of data, and $z$ is the airmass. The best fit
coefficients were $p_{2}$ = 0.2022, $p_{1}$ = -0.0374, and $p_{0}$ = 0.9958
for the pre-minimum set, and $p_{2}$ = -0.0066, $p_{1}$ = 0.01097, and $p_{0}$
= 1.0053 for the post-minimum set. This airmass detrending still did not
remove all systematic trends from the light curve, so we applied a Gaussian
Process (GP) using George [49] to fit the remaining residuals using an
exponential-squared kernel. The best-fit amplitude and scale of this kernel
function were 0.00860 and 0.00057, respectively. The results of these
detrending fits are agnostic to the precise time of conjunction, so we were
confident in performing this detrending prior to fully fitting this light
curve simultaneously with the other data sets. The primary reason for doing so
is computation efficiency, as the GP parameter sampling was too slow to
justify including it in the full routine.

Spitzer/IRAC has observed the WASP-107 system multiple times in both transit
and eclipse. We downloaded the publicly available light curve of one transit
observation of WASP-107b on April 26, 2017 using IRAC’s channel 2 ($\sim$4.5
$\mu$m) observed for Program 13052 (PI: M. Werner) from the Spitzer Heritage
Archive. This program also collected a transit observation using IRAC’s
channel 1 ($\sim$3.6 $\mu$m). We chose to only use the channel 2 dataset
because its bandpass does not overlap with our F322W2 observations, both
reducing duplicate information and providing more robustness against the bias
of limb asymmetry on our time of conjunction measurement, and because
reductions of channel 1 data are more susceptible to strong instrument
systematics that would likely hinder getting a precise tc constraint from the
data. Our reduction of the data generally follows the process of ref.[50].
Starting with the basic calibrated data, we first perform background
subtraction by masking the middle of each image, 10 pixels in from each edge,
then subtracting the median value of the unmasked pixels. To precisely
determine the position of the star in each image, we first performed three
iterations of 5-$\sigma$ clipping on each pixel’s time series to identify bad
pixels, and replaced their values with the median value of their time series.
Then, we used the Howell centroiding technique [51] to determine the position
of WASP-107 in each image and recorded it. Afterward, we returned to the pre-
bad pixel correction images and performed photometric extraction. We extracted
the flux within a 2.3-pixel radius circular aperture centered on WASP-107. We
tested various aperture sizes and chose the one that minimized the resulting
light curve scatter. The light curve showed a strong ramp at the beginning
likely due to detector persistence, so we cut out the first 2,000 points.

Finally, JWST/MIRI observed one transit of WASP-107b in Low Resolution
Spectroscopy (LRS) mode during program JWST-GTO-1280 (PI: P.O. Lagage). We
were given use of the broad-band (5–12 $\mu$m) light curve of this transit
observation for our analysis. The data was reduced using Eureka!, following
the same procedure described in the Eureka! v1 reduction section of ref.[52].
We used a different custom linearity correction than ref.[26], but verified
that this did not ultimately change our result.

###  3.4 Lightcurve Fitting

####  3.4.1 Limb Darkening Treatment

Before fitting any of the data, we computed the expected limb darkening
coefficients for each instrument’s bandpass to use as priors. We used the
online Exoplanet Characterization Toolkit Limb Darkening Calculator [53] with
the Kurucz ATLAS9 stellar model grid, with parameters for WASP-107 of Teff =
4425 K, log$\left(g\right)$ = 4.63 dex, and [Fe/H] = 0.02 dex. We computed the
corresponding model quadratic limb darkening coefficients for each
instrument’s bandpass. The coefficient values pertaining to the ancillary data
fit priors are listed in Table 1, and all values including those used for the
spectroscopic F322W2 light curves are listed in Table 3 of the Extended Data.
When fit freely with no connection between wavelengths, limb darkening
coefficient profiles can take unphysical, non-smooth, non-self-consistent
shapes as a function of wavelength. We use these model coefficients as
Gaussian priors in our fit of the ancillary data to ensure that the
coefficients are physical and self-consistent, while still allowing their
uncertainty to be reflected in the ultimate result. We fix these coefficients
when fitting the spectroscopic F322W2 data to reduce the number of free
parameters, but verified that this decision did not affect our result (Section
1.2.2).

####  3.4.2 Uniform-limb transit model

When fitting the data assuming or forcing the planet to have uniform limbs, we
use the transit modelling package batman [54]. This code models the planetary
disk in transit as a perfect circle, where the radius of this circle varies
only with wavelength. We allowed the orbital parameters to be free for our
ancillary data fit described below, and fixed them to the result of that fit
(Table 1) when analyzing our JWST/NIRCam F322W2 spectroscopic data.

####  3.4.3 Ancillary Data Fits

We fit the seven ancillary transit light curves simultaneously with the radial
velocity data. The transit model we use for the transit observations is
described in Section 3.4.2. Using the Markov Chain Monte Carlo method with the
package emcee, we sampled the time of conjunction $t_{c}$, orbital period
($P$) as $\log_{10}\left(P\right)$, semi-major axis ($a$) as
$\log_{10}\left(a/R_{\star}\right)$, orbital inclination ($i$) as
$\cos\left(i\right)$, eccentricity ($e$) and argument of periastron ($\omega$)
as $\sqrt{e}\sin\left(\omega\right)$ and $\sqrt{e}\cos\left(\omega\right)$,
individual band uniform-limb planet-star radius ratios $R_{p}/R_{\star}$,
quadratic limb darkening coefficients $u_{1}$ and $u_{2}$ for each
instrument/filter, radial velocity semi-amplitude $K$, and radial system
velocity $\gamma$. For each transit light curve except the SOAR data, we also
applied a linear flux vs. time systematics ramp model of the form

$f=a_{1}x+a_{0},$ (2)

where $x\equiv t-t_{median}$, and we also fit for $a_{1}$ and $a_{0}$ of each
individual light curve. We did not apply a linear systematic ramp to the SOAR
data because any such systematic trend would have been removed from our pre-
detrending of that data.

When fitting data through a Bayesian inference method, such as MCMC, it is
common practice to apply priors to fitting parameters based on previous,
independent measurements. For our analysis, we must be careful to select
priors which could not have been affected by unaccounted-for limb asymmetry.
For this reason, we do not enforce any priors on the transit-related
parameters except for the semi-major axis, inclination, and limb darkening
coefficients. In other words, the values of these parameters do not factor
into the prior probability value. Since the semi-major axis and inclination
symmetrically affect the transit light curve on either side of $t_{c}$, these
measurements should not be affected by unaccounted-for limb asymmetry. We
therefore felt safe in applying priors to the semi-major axis and inclination,
based on the independent measurements by ref.[2] using K2 observations.
Similarly, limb darkening is a property of only the star and the instrument
being used, so it would also not be affected. We applied Gaussian priors to
each instrument’s limb darkening coefficients based on an ATLAS stellar model,
as described in Section 3.4.1. For the RV-related parameters, we did not apply
any priors since the only other measurements (refs.[4, 5]) are derived from
the same data we reuse here.

We ran this sampling for 10,500 steps, which was sufficient for all parameters
to converge. The resulting best-fit values and derived parameter values are
listed in Table 1. The data and best-fit models are plotted in Figure 7 of the
Extended Data, and we highlight our constraints on WASP-107b’s orbital
parameters in Figure 8 of the Extended Data.

####  3.4.4 Asymmetric-limb model

When fitting the spectroscopic JWST/NIRCam F322W2 light curves and
investigating for limb asymmetry, we used the transit modeling package
catwoman [22, 24], which models a planet with asymmetric limbs as the
combination of two semi-circles with independent radii. This effectively
splits the uniform disk modeled by batman in two along a central axis. The
orientation of this axis with respect to the planet’s orbital plane is an
adjustable parameter in catwoman, and we fixed this axis to be aligned with
the planet’s orbital axis (see our discussion of this assumption in Section
3.4.5). Again, we fixed the orbital parameters and limb darkening
coefficients, and sampled the time of conjunction, two planet-star radius
ratios per wavelength channel, as well as the slope and intercept of a linear
flux vs. time ramp. As mentioned in the Main Text, we applied a Gaussian prior
to the time of conjunction based on the posterior distribution from our
ancillary data fit. We ran each channel’s MCMC sampling for 10,000 steps,
which was more than sufficient for all parameters to converge.

####  3.4.5 WASP-107b’s Obliquity and Orbital Misalignment

Fits of the current radial velocity measurements of WASP-107 favor WASP-107b
having a small but non-zero eccentricity [5, and this work]. Further, starspot
crossing frequencies [2] and measurements of the Rossiter-McLaughlin effect in
this system [55] suggest that WASP-107b may be on a misaligned orbit. This
eccentricity and misalignment are indicative of dynamical interactions that
occurred in the planet’s past, and may still be operating today. Here, we
discuss the impact of these effects on our experimental design and results.

From our fit of the radial velocity data, we measure WASP-107b’s eccentricity
to be e = 0.05 $\pm$ 0.01 (Table 1), consistent with that of ref.[5] which
analyzed the same data. Neither this eccentricity, or any parameter in our
fit, is affected by the Rossiter-McLaughlin (RM) effect. The RM effect only
affects radial velocity measurements made during the planet’s transit, and all
of the measurements we use were made outside of WASP-107b’s transit. Also, the
measurement of the RM effect in this system by ref.[55] was done by assuming
and fixing WASP-107b’s orbital parameters to those derived by ref.[5], with
whom we are fully consistent. Relatedly, the detailed models of ref.[56] show
that the RM effect does not bias transmission spectra observed by JWST,
especially at the very low spectral resolution that we bin our data to.

In deriving our limb spectra (Figure 2), our catwoman-based transit models
[24, 22] assume that WASP-107b’s rotation axis is aligned to its orbital axis.
Tidal interaction models find that both rotational axis alignment (i.e., zero
rotational obliquity) and synchronization (i.e., equal rotation and orbital
periods) should be inevitable and happens very quickly for short period (P $<$
10 days) gaseous planets [57, 58, 59, 60]. Further, this alignment and
synchronization can occur even if the orbit has not fully circularized [58].
First derived by ref.[61] and rewritten by ref.[62], the timescale $\tau$ for
this alignment can be estimated as

$\tau\approx
0.067\times\left(\frac{Q}{10^{5}}\right)\left(\frac{\omega_{p}M_{p}R_{J}^{3}}{\omega_{J}M_{J}R_{p}^{3}}\right)\left(\frac{P}{10~{}day}\right)^{4}~{}Gyr.$
(3)

Here, $Q$ is the tidal quality factor, and $\omega$, $M$, and $R$ are the
rotation rate, mass, and radius of the planet (subscript $p$) and Jupiter
(subscript $J$). Assuming WASP-107b is synchronously rotating, so that
$\omega_{p}=2.02\times 10^{-6}~{}s^{-1}$, we calculate an alignment timescale
of $\sim$980 years for a Neptune-like $Q$ of $10^{4}$, or $\sim$9800 years for
a Jupiter-like $Q$ of $10^{5}$. This is extremely short in either case,
especially compared to its orbital circularization timescale of $\sim$66 Myr
[5] and the estimated 3.4 - 8.3 Gyr age of the system [1, 5]. Therefore, even
if WASP-107b’s orbit is misaligned to its star’s rotation axis, the planetary
rotation axis should almost certainly be aligned.

###  3.5 Constraining WASP-107b’s time of conjunction

As described in Section 3.4.3, we fit the ancillary transit and radial
velocity observations described in Sections 1.1 and 3.3 together
simultaneously. The goal of this fit was to precisely constrain WASP-107b’s
orbit before fitting our NIRCam/F322W2 spectroscopic data. These orbital
parameters must be known precisely a priori to enable accurately separating
the contribution of each limb to the observed transmission signal. Uncertainty
in the time of conjunction ($t_{c}$), in particular, is degenerate with the
effect of limb asymmetry on the planet’s transit light curve [7, 21, 9, 22].
Using these ancillary data, we were able to strongly constrain each of
WASP-107b’s main orbital parameters (Table 1). We highlight the resulting
posterior distributions of several primary orbital parameters from this fit in
Figure 8 of the Extended Data. Our best-fit values are consistent with
previous measurements from various other datasets [2, 1, 4, 5, 63] and, most
importantly, we achieve very high precision on WASP-107b’s $t_{c}$. At the
epoch of our NIRCam F322W2 observation, we measure WASP-107b’s $t_{c}$ to a
1-$\sigma$ precision of 0.70 seconds. This tight constraint is enabled by the
wide wavelength coverage of these ancillary data, and the unprecedented
combination of high photometric precision and fast integration cadence of the
included JWST data.

The degeneracy between limb asymmetry and $t_{c}$ offsets naturally raises the
question whether our best-fit value may be biased by limb asymmetry in our
ancillary data. Our approach combining multi-wavelength measurements ensures
such bias is eliminated, even if present in individual channels. This multi-
wavelength approach is the current best method based on discussions of this
problem in the literature [9, 22], and has been shown to eliminate any limb-
asymmetry induced bias to the measured $t_{c}$ even when each data-set
individually exhibits limb asymmetry of the same polarity [9]. Ref.[9]
describes the effectiveness of this approach very well, which we summarize
here. This method takes advantage of the fact that limb asymmetry, even if
present at all wavelengths, will be a chromatic effect – the asymmetry will be
different at different wavelengths. When fit simultaneously, each wavelength
will attempt to “pull” $t_{c}$ by different amounts which, since $t_{c}$ is
shared between every wavelength, leads to general disagreement between
channels and a poor fit overall. In effect, the offset in each channel works
against each other. The best solution is the true geometric time of
conjunction, which was proven to be recovered in each test of this method by
ref.[9].

This approach to deriving $t_{c}$ assumes that limb asymmetry varies
chromatically, though this cannot be empirically known a priori. However, the
fact that we derive different morning and evening spectra in one data set
(Figure 2) is proof that the exact same limb asymmetry is not uniformly
present at all wavelengths. For instance, when fitting our broadband NIRCam
F322W2 data by itself assuming a uniform-limb model and allowing $t_{c}$ to
freely vary, we get a best-fit $t_{c}$ that is 6.64 $\pm$ 1.91 seconds later
than the joint-fit value (Figure 9). If the exact same limb asymmetry were
present at all wavelengths and the biases did not cancel out as previously
described, then our joint-fit would have given a time offset by this same
value. If we then applied this offset time in our catwoman fits to the F322W2
spectroscopic data, we would not have observed any significant limb asymmetry.
In this case, the assumed time of conjunction would be biased by an amount
corresponding to the underlying asymmetry, which would have erased the
asymmetry’s signal.

To further evaluate our joint-fit’s effectiveness, we tested fitting each
ancillary data set individually using a uniform-limb model and allowing
$t_{c}$ to freely vary. We then compared the difference between each
individually fit $t_{c}$ and the joint-fit $t_{c}$. As just mentioned, when
fitting the broadband JWST/NIRCam F322W2 data, we get a best-fit $t_{c}$ that
6.64 $\pm$ 1.91 seconds later than the joint-fit result from the broadband
light curve, and we also see chromatic offsets in the best-fit $t_{c}$ shown
in Figure 9 from the spectroscopic data. Both the SOAR and Spitzer/IRAC Ch.2
observations had differences consistent with zero (7.6 $\pm$ 19.7 s and 4.42
$\pm$ 7.39 s, respectively). The TESS data, which overlaps with the SOAR
bandpass, had a best-fit offset of -33.5 $\pm$ 19.0 s. This is negative and
non-zero at over 1$\sigma$ significance, but it is not clear whether such
large offset is real given the sparse 30-minute cadence of the data. The
JWST/NIRCam F210M data, which was observed simultaneously with the NIRCam
F322W2 data with the same cadence, had an offset of -1.11 $\pm$ 1.06 s. This
is also negative but nearly consistent with zero. Finally, the JWST/MIRI LRS
data had a difference of 2.62 $\pm$ 1.28 s which is non-zero by just over
2$\sigma$. Due to their superior data quality, our joint-fit is most effected
by the JWST/NIRCam F210M and JWST/MIRI LRS observations. These individual fits
show that both data sets may also be subject to some limb asymmetry, though
not to the same degree as the JWST/NIRCam F322W2 data, and potentially of
opposite polarity. Our multi-wavelength approach will therefore work as
described to mitigate these individual biases. All together, these fits
underscore the importance of combining data at multiple wavelengths. The mean
offset of these individual fit values, weighted by their uncertainties
($1/\sigma_{t_{c}}^{2}$), is 0.41 s which is well within the joint-fit value’s
uncertainty. Even if we ignore the largest offset, which was the TESS value,
the mean offset is still only 0.47 s. To the best of our current knowledge,
the $t_{c}$ we derive from our joint-fit is as accurate and bias-free as
possible.

The above discussion raises an additional question of whether there is some
systematic error in our JWST/NIRCam F322W2 data that causes the derived
transit time to be different from the joint-fit value, and thus be the cause
of the limb asymmetry we observe. We can directly compare the F322W2 data to
the F210M data as these were observed simultaneously. Both data sets were
reduced using the same method, so there is no difference in how the data time-
stamps are calculated or corrected (e.g., barycentric correction). NIRCam’s
short (F210M) and long (F322W2) channels use separate detectors and the
observatory software does not command synchronization, leading to slight
differences in the exact timestamps between each light curve. On average, the
timestamps of our F322W2 points are 0.41 s later than the corresponding F210M
points. This timestamp offset is within the uncertainty of our joint-fit
derived $t_{c}$, so it cannot explain the 6.64 $\pm$ 1.91 s $t_{c}$ difference
derived from the F322W2 data. To verify this, we tested adding this timestamp
offset to the F210M light curve points, then repeating the individual uniform-
limb fit. The resulting offset was not significantly different (now -0.76
$\pm$ 1.08 s). We got a similar result when instead testing subtracting this
offset from the F322W2 points. Therefore, there is no evidence that any such
systematic errors are present and driving our result.

###  3.6 Forward Modeled Transmission Spectra

We compared our observational results to 1D atmosphere models for the morning
and evening limbs to provide an initial estimate of the atmospheric parameters
on each terminator. The model transmission spectra presented in Figure 2 and
discussed in the main text were generated using the ScCHIMERA model. This
model was first described in ref.[65] with additional updates in refs.[66, 67,
68], but we provide a short summary here. See also ref.[44] for a more
detailed description and most recent updates (including opacity sources). In
short, the ScCHIMERA tool is a self-consistent 1D radiative-convective-
photochemical equilibrium solver (1D-RCPE). Chemistry can be computed in
either thermochemical equilibrium or coupled with the VULCAN chemical kinetics
code [69] to self-consistently compute disequilibrium abundances arising from
photochemistry and vertical mixing. We initially produced a grid of 4830
WASP-107 specific chemical-disequilibrium 1D-RCPE models over a range of
metallicities (23 points from -0.5 - +2.25, where the value is a base-10
logarithm), redistributions (scaling to the incident stellar flux, 14 points
from 0.46 - 1.4 where 1 is full, 2 is dayside), and carbon-to-oxygen ratios
(15 points from 0.1 - 0.75). To fit the spectra we post-process the
atmospheric structure (T-P profile and gas volume mixing ratios) with the
CHIMERA [70] transmission spectrum tool while also adjusting for an opaque
gray, vertically uniform cloud. We parameterize this cloud as a “gray” cross-
section that is added to the gas absorption cross-sections. It can be taken as
the product of a cloud droplet mixing ratio times an absorption cross-section
but does not consider cloud microphysics.

From the grid of models described above, we used an atmospheric metallicity of
10 times the solar value and carbon-to-oxygen ratio of 0.35 based on the best-
fit retrieval values on HST/WFC3 and JWST/MIRI observations of this planet
from ref.[26]. These values are also consistent with the additional retrievals
in Welbanks, et al. (in review). We expect the metallicity and carbon-to-
oxygen ratio to be globally uniform quantities, and do not vary them between
the models for each limb. In this work, we do not seek to make specific,
precise constraints on global properties, but rather develop an understanding
for what may be causing the observed difference between the evening and
morning limbs. To do this, we performed a grid fit over the heat
redistribution parameter, over its full range, and the gray cloud opacity,
ranging from 10-2.58 to 10-0.78 cm2 g-1. This heat redistribution parameter is
effectively the temperature at the terminator, related by
$T=738.31\times\text{redistribution}^{0.25}$. We computed the $\chi^{2}$ value
between the grid model and the observed morning and evening spectra
individually, following the $\chi^{2}$-minimization technique of refs.[40,
41].

The results of our grid fit are shown in Figure 4 and discussed in the main
text. The best-fit model spectra are shown along with the observed limb
spectra in Figure 2. In Figure 4, we show the $\chi^{2}$ distributions of the
morning-limb fit in blue and the evening-limb fit in orange via contours drawn
at the 68% (i.e., 1-$\sigma$) and 99% (i.e., 3-$\sigma$) confidence intervals.
The distributions for each limb do not overlap even well beyond their 99%
confidence intervals, showing that the atmospheric properties on WASP-107b’s
terminators are significantly different. The morning-limb distribution is
against the lower edge of our temperature grid with a minimum $\chi^{2}$ at
611 K (redistribution parameter of 0.47), and is below 630 K (667 K) at 68%
(99%) confidence. The morning-limb uniform gray cloud opacity is well
constrained between 10-1.98 and 10-1.58 cm2 g-1 at 99% confidence, though as
we have noted this is a simple initial exploration and the true nature of
clouds or hazes on WASP-107b are likely more complex, and beyond the scope of
this work. The evening-limb distribution is bimodel, with a cooler/clearer
solution (T $\approx$ 732 K, $\kappa_{cloud}\approx 10^{-0.78}$ cm2 g-1) and a
hotter/cloudier solution (T $\approx$ 760 K, $\kappa_{cloud}\approx
10^{-1.58}$ cm2 g-1). In either case, these fits suggest that the evening-limb
is significantly hotter than the morning-limb, and that there is a slight
difference in the cloudiness as well.

Highlighted by the apparent bimodality of our evening-limb distribution, we
see a correlation between the temperature and cloud opacity. This makes sense
as both generally affect the amplitude of spectral features. To get a better
sense for the individual effect of each parameter, we show two additional
pairs of models in Figure 2 – one with fixed temperature but varied cloud
opacity, and vice-versa. As discussed in the Main text, we find that the
temperature has the strongest effect on both the relative shapes of the
evening and morning spectra, as well as their relative offset. In our models,
the cloud opacity (at least within this range of interest) also affects the
offset, but does not as strongly change the shape of each spectrum. This
precise behavior is obviously subject to our parametrization of these clouds
as a vertically uniform gray opacity source. For example, more detailed cloud
or haze models, especially focused on specific species, may affect the shapes
differently or even induce relative slopes across the spectra. Disentangling
these details will require additional observations over a wider wavelength
coverage as well as much more complex multi-dimensional modeling.
Nevertheless, our results show that recovering the separate properties of a
transiting exoplanet’s morning and evening terminators is possible, and that
WASP-107b exhibits significant morning-to-evening asymmetry.

## Acknowledgements

The authors would like to acknowledge Néstor Espinoza and Karl Misselt for
assistance in verifying the reliability of JWST’s instrument time-stamps,
Pierre-Olivier Lagage for providing the JWST/MIRI LRS data observed in JWST
program 1280, Charles Beichman for suggesting we observe this planet as part
of our JWST guaranteed time program, Kevin Hardegree-Ullman for help with
using the package TTVFaster, and Néstor Espinoza for helpful discussion about
this manuscript. T.P.G. and T.J.B. acknowledge funding from NASA in WBS
411672.07.05.05.03.02. M.M., E. S., and M.R. acknowledge funding from NASA
Goddard Spaceflight Center via NASA contract NAS5-02105. This work benefited
from the 2023 Exoplanet Summer Program in the Other Worlds Laboratory (OWL) at
the University of California, Santa Cruz, a program funded by the Heising-
Simons Foundation and NASA.

## Declarations

  * •

Funding: T.P.G. and T.J.B. acknowledge funding from NASA in WBS
411672.07.05.05.03.02. M.M., E.S., and M.R. acknowledge funding from NASA
Goddard Spaceflight Center via NASA contract NAS5-02105.

  * •

Conflict of interest/Competing interests: None

  * •

Ethics approval: Not applicable

  * •

Consent to participate: Not applicable.

  * •

Consent for publication: Not applicable.

  * •

Availability of data and materials: The JWST/NIRCam and JWST/MIRI data used in
this analysis have completed their one-year proprietary period and are
currently available in the Mikulski Archive for Space Telescopes (MAST). The
other data used in this analysis are also publicly available in MAST, or can
be provided upon request.

  * •

Code availability: The three reduction pipelines used in this work (Eureka!,
t-shirt, and Pegasus) are all either currently open-source or are planned to
be made open-source. The codes of each are hosted on Github, and links for
each may be found in the Methods section.

  * •

Authors’ contributions: M.M. led this analysis including formulating the idea,
performing the data and supporting analyses, and was the primary author of the
paper. T.G.B. contributed the Pegasus reduction of the NIRCam data, E.S.
contributed the t-shirt reduction, and T.J.B. contributed the Eureka!
reduction of the NIRCam F322W2 spectra and the broadband reduction of the
MIRI/LRS data. T.J.B., E.S., and T.G.B. also contributed to the interpretation
of our results and the writing of the paper. M.L. contributed the modeling
analysis using ScCHIMERA. T.P.G. selected the object, devised the observing
strategy, and contributed to the analysis. V.P., E.R., J.F., M.L., and L.W.
each contributed to the methodology and interpretation of the results. M.R.
designed the JWST/NIRCam instrument used for the observations, and provided
the observing time necessary for completing these observations.

##  Appendix A Extended Data

Table 1: Priors for fitting the ancillary radial velocity and transit observations of WASP-107b, as well as the best-fit parameter values from these fits. The fitted time of conjunction corresponds to the actual epoch of our JWST/NIRCam F210M / F322W2 transit observation. Note that we use the convention for $\omega$ used by the RadVel code [64].  Parameter | Units | Value |
---|---|---|---
Priors — |  |  |
$t_{c}$ | Time of Conjunction [BJD${}_{\text{TDB}}$] | None |
P | Orbital Period [day] | None |
a / R⋆ | Semi-major axis | 18.164 $\pm$ 0.037 [2] |
i | Inclination [deg] | 89.8 $\pm$ 0.2 [2] |
e | Eccentricity | None |
$\omega$ | Argument of Periastron [deg] | None |
$\left(u_{1},u_{2}\right)_{\text{TESS}}$ | Quadratic LD coefficients, TESS | $\left(0.541\pm 0.005,0.106\pm 0.007\right)$ |
$\left(u_{1},u_{2}\right)_{\text{i-band}}$ | Quadratic LD coefficients, i-band | $\left(0.535\pm 0.001,0.120\pm 0.001\right)$ |
$\left(u_{1},u_{2}\right)_{\text{F210M}}$ | Quadratic LD coefficients, NIRCam F210M | $\left(0.13\pm 0.01,0.26\pm 0.01\right)$ |
$\left(u_{1},u_{2}\right)_{\text{IRAC}}$ | Quadratic LD coefficients, IRAC Ch.2 | $\left(0.10\pm 0.01,0.12\pm 0.01\right)$ |
$\left(u_{1},u_{2}\right)_{\text{MIRI}}$ | Quadratic LD coefficients, MIRI broad-band | $\left(0.076\pm 0.001,0.089\pm 0.001\right)$ |
Astrophysical Parameter Best-fit Values — |  |  |
$t_{c}$ | Time of Conjunction [BJD${}_{\text{TDB}}$] | 2459958.747244 $\pm$ 8.2 $\times$ 10-6 |
$\log_{10}\left(P\right)$ | log Orbital period [days] | 0.75750893 $\pm$ 2.9 $\times$ 10-8 |
$\log_{10}\left(a/R_{\star}\right)$ | log semi-major axis [R⋆] | 1.25639 $\pm$ 0.0024 |
$\cos~{}i$ | Cosine of orbital inclination | 0.0074 $\pm$ 0.0005 |
$R_{p}/R_{\star}$ | Radius ratio, TESS band | 0.1421 $\pm$ 0.0005 |
$R_{p}/R_{\star}$ | Radius ratio, i-band (SOAR/Goodman) | 0.1435 $\pm$ 0.0003 |
$R_{p}/R_{\star}$ | Radius ratio, JWST/NIRCam F210M | 0.143982 $\pm$ 0.00006 |
$R_{p}/R_{\star}$ | Radius ratio, Spitzer/IRAC Ch.2 | 0.1407 $\pm$ 0.0004 |
$R_{p}/R_{\star}$ | Radius ratio, JWST/MIRI broadband | 0.14368 $\pm$ 0.00004 |
$K$ | RV semi-amplitude [m/s] | 13.74 $\pm$ 0.31 |
$\gamma$ | System velocity [km/s] | -0.0022 $\pm$ 0.0002 |
$\sqrt{e}\cos\left(\omega\right)$ |  | 0.234 $\pm$ 0.024 |
$\sqrt{e}\sin\left(\omega\right)$ |  | -0.009 $\pm$ 0.024 |
Derived Parameter Values — |  |  |
$P$ | Orbital period [days] | 5.72148722 $\pm$ 3 $\times$ $10^{-7}$ |
$a/R_{\star}$ | Semi-major axis [R⋆] | 18.05 $\pm$ 0.1 |
$i$ | Orbital inclination [deg] | 89.57 $\pm$ 0.03 |
$e$ | Eccentricity | 0.05 $\pm$ 0.01 |
$\omega$ | Argument of Periastron [deg] | -2.3 $\pm$ 6.1 |
|  |  |

Table 2: Best-fit limb transit depths from our fit to WASP-107b’s JWST/NIRCam F322W2 transit light curves. The combined limb depth is the sum of the evening and morning limb depths, accounting for their covariance. Note that we used constant wavelength spacing for our bins. Wavelength [$\mu$m] | Evening limb depth [ppm] | Morning limb depth [ppm] | Total transit depth [ppm] | Limb depths covariance
---|---|---|---|---
2.475 $\pm$ 0.025 | 10397 $\pm$ 126 | 10077 $\pm$ 125 | 20474 $\pm$ 49 | -1.5 $\times$ 10-8
2.525 $\pm$ 0.025 | 10280 $\pm$ 122 | 10400 $\pm$ 121 | 20681 $\pm$ 39 | -1.4 $\times$ 10-8
2.575 $\pm$ 0.025 | 10546 $\pm$ 111 | 10197 $\pm$ 111 | 20744 $\pm$ 43 | -1.2 $\times$ 10-8
2.625 $\pm$ 0.025 | 10504 $\pm$ 105 | 10233 $\pm$ 105 | 20737 $\pm$ 45 | -1.0 $\times$ 10-8
2.675 $\pm$ 0.025 | 10837 $\pm$ 97 | 10247 $\pm$ 96 | 21084 $\pm$ 40 | -8.6 $\times$ 10-9
2.725 $\pm$ 0.025 | 10798 $\pm$ 104 | 10313 $\pm$ 103 | 21112 $\pm$ 39 | -1.0 $\times$ 10-8
2.775 $\pm$ 0.025 | 10722 $\pm$ 101 | 10245 $\pm$ 99 | 20968 $\pm$ 42 | -9.2 $\times$ 10-9
2.825 $\pm$ 0.025 | 10665 $\pm$ 106 | 10293 $\pm$ 104 | 20959 $\pm$ 42 | -1.0 $\times$ 10-8
2.875 $\pm$ 0.025 | 10541 $\pm$ 114 | 10276 $\pm$ 115 | 20817 $\pm$ 47 | -1.2 $\times$ 10-8
2.925 $\pm$ 0.025 | 10522 $\pm$ 103 | 10275 $\pm$ 103 | 20798 $\pm$ 37 | -1.0 $\times$ 10-8
2.975 $\pm$ 0.025 | 10498 $\pm$ 98 | 10231 $\pm$ 96 | 20729 $\pm$ 34 | -8.9 $\times$ 10-9
3.025 $\pm$ 0.025 | 10490 $\pm$ 100 | 10197 $\pm$ 99 | 20687 $\pm$ 36 | -9.4 $\times$ 10-9
3.075 $\pm$ 0.025 | 10405 $\pm$ 97 | 10251 $\pm$ 94 | 20656 $\pm$ 41 | -8.4 $\times$ 10-9
3.125 $\pm$ 0.025 | 10486 $\pm$ 102 | 10070 $\pm$ 101 | 20557 $\pm$ 32 | -9.9 $\times$ 10-9
3.175 $\pm$ 0.025 | 10348 $\pm$ 88 | 10150 $\pm$ 88 | 20499 $\pm$ 36 | -7.2 $\times$ 10-9
3.225 $\pm$ 0.025 | 10375 $\pm$ 90 | 10133 $\pm$ 90 | 20508 $\pm$ 35 | -7.6 $\times$ 10-9
3.275 $\pm$ 0.025 | 10283 $\pm$ 91 | 10222 $\pm$ 89 | 20505 $\pm$ 34 | -7.6 $\times$ 10-9
3.325 $\pm$ 0.025 | 10546 $\pm$ 88 | 10017 $\pm$ 87 | 20563 $\pm$ 34 | -7.2 $\times$ 10-9
3.375 $\pm$ 0.025 | 10417 $\pm$ 91 | 10089 $\pm$ 90 | 20506 $\pm$ 35 | -7.6 $\times$ 10-9
3.425 $\pm$ 0.025 | 10232 $\pm$ 97 | 10128 $\pm$ 96 | 20360 $\pm$ 40 | -8.6 $\times$ 10-9
3.475 $\pm$ 0.025 | 10326 $\pm$ 91 | 10079 $\pm$ 91 | 20406 $\pm$ 38 | -7.7 $\times$ 10-9
3.525 $\pm$ 0.025 | 10188 $\pm$ 94 | 10130 $\pm$ 94 | 20319 $\pm$ 38 | -8.2 $\times$ 10-9
3.575 $\pm$ 0.025 | 10339 $\pm$ 95 | 9924 $\pm$ 95 | 20264 $\pm$ 34 | -8.6 $\times$ 10-9
3.625 $\pm$ 0.025 | 10284 $\pm$ 98 | 10035 $\pm$ 98 | 20319 $\pm$ 43 | -8.7 $\times$ 10-9
3.675 $\pm$ 0.025 | 10263 $\pm$ 111 | 9908 $\pm$ 112 | 20172 $\pm$ 48 | -1.1 $\times$ 10-8
3.725 $\pm$ 0.025 | 10205 $\pm$ 99 | 10079 $\pm$ 98 | 20284 $\pm$ 40 | -9.0 $\times$ 10-9
3.775 $\pm$ 0.025 | 10334 $\pm$ 101 | 9867 $\pm$ 100 | 20201 $\pm$ 38 | -9.5 $\times$ 10-9
3.825 $\pm$ 0.025 | 10222 $\pm$ 102 | 10040 $\pm$ 100 | 20262 $\pm$ 44 | -9.4 $\times$ 10-9
3.875 $\pm$ 0.025 | 10317 $\pm$ 110 | 9912 $\pm$ 109 | 20230 $\pm$ 34 | -1.1 $\times$ 10-8
3.925 $\pm$ 0.025 | 10280 $\pm$ 85 | 9970 $\pm$ 84 | 20251 $\pm$ 35 | -6.6 $\times$ 10-9
|  |  |  |

Figure 5:  Uniform-limb injection test results. To better interpret our observed limb spectra (Figure 2), we tested what we would derive if there were no limb asymmetry at all. As described in Section 1.2, we generated several realizations of synthetic data with the same Gaussian light curve scatter and cadence of our real observations, based around a transit model with identical morning and evening limbs. Then, we fit these synthetic data with catwoman as we did with the real data. As in our presentation of the real data in Figure 2, we show the best-fit evening-limb spectrum in orange and the best-fit morning-limb spectrum in blue. The top four panels show the results of four individual tests. The bottom left panel shows all ten tests overlaid on top of each other, and the bottom right panel shows the average limb depth in each channel. In almost all cases, we recover the absence of limb asymmetry in each channel and the best-fit evening and morning transit depths are consistent with one another within 1-$\sigma$. Due to noise in the data, there are occasional outliers where the depths differ by more than 1-$\sigma$, but these occur randomly and have no consistent or regular pattern with wavelength. Also, these random variations switch polarity at random, further suggesting they cannot be real.  Figure 6:  Comparison of the transmission spectra determined for WASP-107b’s morning and evening limbs from each of our reduction methods. Top: the morning and evening limb transmission spectra derived from each reduction of our JWST/NIRCam F322W2 observation. The results from t-shirt are shown as circles, Pegasus as squares, and Eureka! as triangles. Each reduction produces morphologically consistent spectra for both limbs. Bottom: the difference, in terms of the mutual standard deviation, between the limb transit depths derived by each reduction. Those between the evening-limb depths are shown in the upper panel, and those between the morning-limb depths in the lower panel. The exact limb transit depths in each channel vary slightly between reductions, but are all consistent within 1-$\sigma$ except for one outlier. At shorter wavelengths, the morning-limb spectrum from t-shirt is slightly lower than that of the other two reductions due to t-shirt’s improved 1/f noise correction, though this difference is always less than 1-$\sigma$. In our main analysis, we ultimately chose to use the t-shirt reduction due to this improved 1/f noise correction.  Table 3: Quadratic limb darkening coefficients used in this analysis. We computed all coefficients using the Kurucz ATLAS9 stellar model grid available through the online Exoplanet Characterization Toolkit.  Instrument | Wavelength(s) [$\mu$m] | u1 | u2 | Prior type
---|---|---|---|---
TESS | 0.78 $\pm$ 0.21 | 0.541 $\pm$ 0.005 | 0.106 $\pm$ 0.007 | Gaussian
SOAR/Goodman | 0.77 $\pm$ 0.10 | 0.535 $\pm$ 0.001 | 0.120 $\pm$ 0.001 | Gaussian
JWST/NIRCam F210M | 2.09 $\pm$ 0.20 | 0.13 $\pm$ 0.01 | 0.26 $\pm$ 0.01 | Gaussian
JWST/NIRCam F322W2 | 2.475 $\pm$ 0.025 | 0.140 $\pm$ 0.012 | 0.221 $\pm$ 0.017 | Fixed
| 2.525 $\pm$ 0.025 | 0.141 $\pm$ 0.012 | 0.215 $\pm$ 0.016 | Fixed
| 2.575 $\pm$ 0.025 | 0.136 $\pm$ 0.012 | 0.214 $\pm$ 0.016 | Fixed
| 2.625 $\pm$ 0.025 | 0.133 $\pm$ 0.012 | 0.211 $\pm$ 0.016 | Fixed
| 2.675 $\pm$ 0.025 | 0.130 $\pm$ 0.012 | 0.209 $\pm$ 0.016 | Fixed
| 2.725 $\pm$ 0.025 | 0.127 $\pm$ 0.012 | 0.208 $\pm$ 0.016 | Fixed
| 2.775 $\pm$ 0.025 | 0.125 $\pm$ 0.012 | 0.205 $\pm$ 0.016 | Fixed
| 2.825 $\pm$ 0.025 | 0.123 $\pm$ 0.012 | 0.202 $\pm$ 0.016 | Fixed
| 2.875 $\pm$ 0.025 | 0.122 $\pm$ 0.012 | 0.198 $\pm$ 0.016 | Fixed
| 2.925 $\pm$ 0.025 | 0.122 $\pm$ 0.011 | 0.195 $\pm$ 0.015 | Fixed
| 2.975 $\pm$ 0.025 | 0.121 $\pm$ 0.011 | 0.192 $\pm$ 0.015 | Fixed
| 3.025 $\pm$ 0.025 | 0.120 $\pm$ 0.011 | 0.189 $\pm$ 0.015 | Fixed
| 3.075 $\pm$ 0.025 | 0.119 $\pm$ 0.011 | 0.186 $\pm$ 0.015 | Fixed
| 3.125 $\pm$ 0.025 | 0.118 $\pm$ 0.011 | 0.183 $\pm$ 0.015 | Fixed
| 3.175 $\pm$ 0.025 | 0.118 $\pm$ 0.011 | 0.180 $\pm$ 0.014 | Fixed
| 3.225 $\pm$ 0.025 | 0.118 $\pm$ 0.011 | 0.177 $\pm$ 0.014 | Fixed
| 3.275 $\pm$ 0.025 | 0.115 $\pm$ 0.010 | 0.176 $\pm$ 0.014 | Fixed
| 3.325 $\pm$ 0.025 | 0.115 $\pm$ 0.010 | 0.173 $\pm$ 0.014 | Fixed
| 3.375 $\pm$ 0.025 | 0.115 $\pm$ 0.010 | 0.170 $\pm$ 0.014 | Fixed
| 3.425 $\pm$ 0.025 | 0.114 $\pm$ 0.010 | 0.167 $\pm$ 0.013 | Fixed
| 3.475 $\pm$ 0.025 | 0.114 $\pm$ 0.010 | 0.165 $\pm$ 0.013 | Fixed
| 3.525 $\pm$ 0.025 | 0.114 $\pm$ 0.010 | 0.163 $\pm$ 0.013 | Fixed
| 3.575 $\pm$ 0.025 | 0.113 $\pm$ 0.010 | 0.161 $\pm$ 0.013 | Fixed
| 3.625 $\pm$ 0.025 | 0.112 $\pm$ 0.010 | 0.159 $\pm$ 0.013 | Fixed
| 3.675 $\pm$ 0.025 | 0.110 $\pm$ 0.009 | 0.157 $\pm$ 0.013 | Fixed
| 3.725 $\pm$ 0.025 | 0.109 $\pm$ 0.009 | 0.155 $\pm$ 0.012 | Fixed
| 3.775 $\pm$ 0.025 | 0.109 $\pm$ 0.009 | 0.153 $\pm$ 0.012 | Fixed
| 3.825 $\pm$ 0.025 | 0.108 $\pm$ 0.009 | 0.151 $\pm$ 0.012 | Fixed
| 3.875 $\pm$ 0.025 | 0.108 $\pm$ 0.009 | 0.150 $\pm$ 0.012 | Fixed
| 3.925 $\pm$ 0.025 | 0.107 $\pm$ 0.009 | 0.148 $\pm$ 0.012 | Fixed
Spitzer/IRAC | 4.49 $\pm$ 0.50 | 0.10 $\pm$ 0.01 | 0.12 $\pm$ 0.01 | Gaussian
JWST/MIRI | 5– 12 | 0.076 $\pm$ 0.001 | 0.089 $\pm$ 0.001 | Gaussian
|  |  |  |
Figure 7: Gallery of ancillary transit and radial velocity data used to
support our analysis. We collected archival transit observations from TESS and
Spitzer/IRAC, observed new data with SOAR/Goodman and JWST/NIRCam F210M, and
received broadband JWST/MIRI LRS data from JWST program 1280 from P.I. Pierre-
Olivier Lagage. The data for TESS are shown as the white points. The data for
the other transit observations are shown as the gray points, with temporally
binned data as the white points. The best-fit models, including the transit
model and systematics models combined, are shown in green. We also collected
archival radial velocity measurements observed from Keck/HIRES and CORALIE [4,
5]. These data are shown phase-folded in the bottom-right panel, and the best-
fit model is shown in green. We used these transit and radial velocity data to
constrain WASP-107b’s orbital parameters for use in fitting our JWST/NIRCam
F322W2 spectroscopic light curves. In particular, with these data we
constrained WASP-107b’s time of conjunction, at the epoch of our NIRCam
observation, to a 1-$\sigma$ precision of 0.7 seconds. We see no significant
evidence for starspot crossings during any of our transit observations. Figure
8:  Posterior distributions of WASP-107b’s primary orbital parameters from our
ancillary data joint-fit. As described in the Methods, we simultaneously fit
transit observations from TESS, SOAR/Goodman, JWST/NIRCam F210M, Spitzer/IRAC,
and JWST/MIRI LRS, as well as radial velocity observations from Keck/HIRES and
CORALIE [4, 5]. These data and best-fit models are shown in Figure 7, and here
we show the posterior distributions for the main planetary orbital parameters.
Each parameter was well constrained. The primary goal of this fit was to
measure the time of conjunction precisely, shown in the left-most panel, which
we did to a 1-$\sigma$ precision of only 0.70 seconds. Figure 9: Apparent
chromatic offsets in WASP-107b’s time of conjunction with wavelength due to
limb asymmetry. Shown are the best-fit times of conjunction as a function of
wavelength when fitting our JWST/NIRCam F322W2 data by itself. We fit the data
using a uniform-limb transit model, described in Section 3.4.2, fitting only
for $t_{c}$, the planet-star radius ratio, and a linear flux vs. time ramp
model. The black points are the result from fitting each of our thirty
spectroscopic light curves. The purple line is the result from fitting the
broadband light curve, with the shaded region representing the 1-$\sigma$
uncertainty. We show the difference between these results and the joint-fit
time of conjunction derived by fitting other multi-wavelength data sets
simultaneously, whose uncertainty is shown as the green shaded region. The
observed difference is due to the significant evening-morning limb asymmetry
on WASP-107b at these wavelengths.

## References

  *   * [1] Močnik, T., Hellier, C., Anderson, D. R., Clark, B. J. M. & Southworth, J.  Starspots on WASP-107 and pulsations of WASP-118.  _MNRAS_ 469, 1622–1629 (2017).
  * [2] Dai, F. & Winn, J. N.  The Oblique Orbit of WASP-107b from K2 Photometry.  _AJ_ 153, 205 (2017).
  * [3] Hejazi, N. _et al._ Elemental Abundances of the Super-Neptune WASP-107b’s Host Star Using High-resolution, Near-infrared Spectroscopy.  _ApJ_ 949, 79 (2023).
  * [4] Anderson, D. R. _et al._ The discoveries of WASP-91b, WASP-105b and WASP-107b: Two warm Jupiters and a planet in the transition region between ice giants and gas giants.  _A &A_ 604, A110 (2017).
  * [5] Piaulet, C. _et al._ WASP-107b’s Density Is Even Lower: A Case Study for the Physics of Planetary Gas Envelope Accretion and Orbital Migration.  _AJ_ 161, 70 (2021).
  * [6] Kataria, T. _et al._ The Atmospheric Circulation of a Nine-hot-Jupiter Sample: Probing Circulation and Chemistry over a Wide Phase Space.  _ApJ_ 821, 9 (2016).
  * [7] Line, M. R. & Parmentier, V.  The Influence of Nonuniform Cloud Cover on Transit Transmission Spectra.  _ApJ_ 820, 78 (2016).
  * [8] Schwartz, J. C., Kashner, Z., Jovmir, D. & Cowan, N. B.  Phase Offsets and the Energy Budgets of Hot Jupiters.  _ApJ_ 850, 154 (2017).
  * [9] Powell, D. _et al._ Transit Signatures of Inhomogeneous Clouds on Hot Jupiters: Insights from Microphysical Cloud Modeling.  _ApJ_ 887, 170 (2019).
  * [10] Powell, D., Zhang, X., Gao, P. & Parmentier, V.  Formation of Silicate and Titanium Clouds on Hot Jupiters.  _ApJ_ 860, 18 (2018).
  * [11] Feng, Y. K. _et al._ The Impact of Non-uniform Thermal Structure on the Interpretation of Exoplanet Emission Spectra.  _ApJ_ 829, 52 (2016).
  * [12] Caldas, A. _et al._ Effects of a fully 3D atmospheric structure on exoplanet transmission spectra: retrieval biases due to day-night temperature gradients.  _A &A_ 623, A161 (2019).
  * [13] Welbanks, L. & Madhusudhan, N.  On Atmospheric Retrievals of Exoplanets with Inhomogeneous Terminators.  _ApJ_ 933, 79 (2022).
  * [14] MacDonald, R. J., Goyal, J. M. & Lewis, N. K.  Why Is it So Cold in Here? Explaining the Cold Temperatures Retrieved from Transmission Spectra of Exoplanet Atmospheres.  _ApJ_ 893, L43 (2020).
  * [15] Schlawin, E., Greene, T. P., Line, M., Fortney, J. J. & Rieke, M.  Clear and Cloudy Exoplanet Forecasts for JWST: Maps, Retrieved Composition, and Constraints on Formation with MIRI and NIRCam.  _AJ_ 156, 40 (2018).
  * [16] Greene, T. _et al._ Oschmann, J., Jacobus M., Clampin, M. C. & MacEwen, H. A. (eds) _NIRCam: development and testing of the JWST near-infrared camera_.  (eds Oschmann, J., Jacobus M., Clampin, M. C. & MacEwen, H. A.) _Space Telescopes and Instrumentation 2010: Optical, Infrared, and Millimeter Wave_ , Vol. 7731 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , 77310C (2010).
  * [17] Rieke, M. J. _et al._ Performance of NIRCam on JWST in Flight.  _PASP_ 135, 028001 (2023).
  * [18] Bell, T. _et al._ Eureka!: An End-to-End Pipeline for JWST Time-Series Observations.  _The Journal of Open Source Software_ 7, 4503 (2022).
  * [19] Schlawin, E. & Glidic, K.  tshirt (2022).  URL https://github.com/eas342/tshirt.
  * [20] Schlawin, E. _et al._ JWST Noise Floor. I. Random Error Sources in JWST NIRCam Time Series.  _AJ_ 160, 231 (2020).
  * [21] von Paris, P., Gratier, P., Bordé, P., Leconte, J. & Selsis, F.  Inferring asymmetric limb cloudiness on exoplanets from transit light curves.  _A &A_ 589, A52 (2016).
  * [22] Espinoza, N. & Jones, K.  Constraining Mornings and Evenings on Distant Worlds: A new Semianalytical Approach and Prospects with Transmission Spectroscopy.  _AJ_ 162, 165 (2021).
  * [23] Clemens, J. C., Crain, J. A. & Anderson, R.  Moorwood, A. F. M. & Iye, M. (eds) _The Goodman spectrograph_.  (eds Moorwood, A. F. M. & Iye, M.) _Ground-based Instrumentation for Astronomy_ , Vol. 5492 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , 331–340 (2004).
  * [24] Jones, K. & Espinoza, N.  catwoman: A transit modelling python package for asymmetric light curves.  _Journal of Open Source Software_ 7, 2382 (2022).  URL https://doi.org/10.21105/joss.02382.
  * [25] Kreidberg, L., Line, M. R., Thorngren, D., Morley, C. V. & Stevenson, K. B.  Water, High-altitude Condensates, and Possible Methane Depletion in the Atmosphere of the Warm Super-Neptune WASP-107b.  _ApJ_ 858, L6 (2018).
  * [26] Dyrek, A. _et al._ SO2, silicate clouds, but no CH4 detected in a warm Neptune.  _Nature_ 625, 51–54 (2024).
  * [27] Benneke, B. & Seager, S.  How to Distinguish between Cloudy Mini-Neptunes and Water/Volatile-dominated Super-Earths.  _ApJ_ 778, 153 (2013).
  * [28] Welbanks, L. & Madhusudhan, N.  Aurora: A Generalized Retrieval Framework for Exoplanetary Transmission Spectra.  _ApJ_ 913, 114 (2021).
  * [29] Agol, E. & Deck, K.  TTVFaster: First order eccentricity transit timing variations (TTVs).  Astrophysics Source Code Library, record ascl:1604.012 (2016).  1604.012.
  * [30] Rackham, B. V., Apai, D. & Giampapa, M. S.  The Transit Light Source Effect: False Spectral Features and Incorrect Densities for M-dwarf Transiting Planets.  _ApJ_ 853, 122 (2018).
  * [31] Rackham, B. V., Apai, D. & Giampapa, M. S.  The Transit Light Source Effect. II. The Impact of Stellar Heterogeneity on Transmission Spectra of Planets Orbiting Broadly Sun-like Stars.  _AJ_ 157, 96 (2019).
  * [32] Oshagh, M. _et al._ Effect of stellar spots on high-precision transit light-curve.  _A &A_ 556, A19 (2013).
  * [33] Ehrenreich, D. _et al._ Nightside condensation of iron in an ultrahot giant exoplanet.  _Nature_ 580, 597–601 (2020).
  * [34] Kesseli, A. Y. & Snellen, I. A. G.  Confirmation of Asymmetric Iron Absorption in WASP-76b with HARPS.  _ApJ_ 908, L17 (2021).
  * [35] Gandhi, S., Kesseli, A., Snellen, I., Brogi, M., Wardenier, J.P., Parmentier, V., Welbanks, L., & Savel, A.B.,  Spatially resolving the terminator: variation of Fe, temperature, and winds in WASP-76 b across planetary limbs and orbital phase.  _MNRAS_ (2022).
  * [36] Maguire, C., Gibson, N.P., Nugroho, S.K., Fortune, M., Ramkumar, S., Gandhi, S., & de Mooij, E.,  High-resolution atmospheric retrievals of WASP-76b transmission spectroscopy with ESPRESSO: Monitoring limb asymmetries across multiple transits.  _A &A_ (2024).
  * [37] Hoeijmakers, H. J. _et al._ High-resolution transmission spectroscopy of MASCARA-2 b with EXPRES.  _A &A_ 641, A120 (2020).
  * [38] Bourrier, V. _et al._ Hot Exoplanet Atmospheres Resolved with Transit Spectroscopy (HEARTS). III. Atmospheric structure of the misaligned ultra-hot Jupiter WASP-121b.  _A &A_ 635, A205 (2020).
  * [39] Borsa, F. _et al._ Atmospheric Rossiter-McLaughlin effect and transmission spectroscopy of WASP-121b with ESPRESSO.  _A &A_ 645, A24 (2021).
  * [40] Avni, Y.  Energy spectra of X-ray clusters of galaxies.  _ApJ_ 210, 642–646 (1976).
  * [41] Wall, J. V. & Jenkins, C. R.  _Practical Statistics for Astronomers_ (Cambridge University Press, 2012).
  * [42] Tsai, S.-M., Moses, J. I., Powell, D. & Lee, E. K. H.  Day-Night Transport-induced Chemistry and Clouds on WASP-39b: Gas-phase Composition.  _ApJ_ 959, L30 (2023).
  * [43] Perez-Becker, D. & Showman, A.P.,  ATMOSPHERIC HEAT REDISTRIBUTION ON HOT JUPITERS.  _ApJ_ 776, 134 (2013).
  * [44] Bell, T. J. _et al._ Methane throughout the atmosphere of the warm exoplanet WASP-80b.  _Nature_ 623, 709–712 (2023).
  * [45] Schlawin, E. _et al._ JWST NIRCam Defocused Imaging: Photometric Stability Performance and How It Can Sense Mirror Tilts.  _PASP_ 135, 018001 (2023).
  * [46] Horne, K.  An optimal extraction algorithm for CCD spectroscopy.  _PASP_ 98, 609–617 (1986).
  * [47] Vogt, S. S. _et al._ Crawford, D. L. & Craine, E. R. (eds) _HIRES: the high-resolution echelle spectrometer on the Keck 10-m Telescope_.  (eds Crawford, D. L. & Craine, E. R.) _Instrumentation in Astronomy VIII_ , Vol. 2198 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , 362 (1994).
  * [48] Collins, K. A., Kielkopf, J. F., Stassun, K. G. & Hessman, F. V.  ASTROIMAGEJ: IMAGE PROCESSING AND PHOTOMETRIC EXTRACTION FOR ULTRA-PRECISE ASTRONOMICAL LIGHT CURVES.  _The Astronomical Journal_ 153, 77 (2017).  URL https://doi.org/10.3847%2F1538-3881%2F153%2F2%2F77.
  * [49] Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W. & O’Neil, M.  Fast Direct Methods for Gaussian Processes.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 38, 252 (2015).
  * [50] Beatty, T. G. _et al._ A significant overluminosity in the transiting brown dwarf CWW 89ab.  _The Astronomical Journal_ 156, 168 (2018).  URL https://doi.org/10.3847/1538-3881/aad697.
  * [51] Howell, S. B.  _Handbook of CCD Astronomy_ 2 edn.  Cambridge Observing Handbooks for Research Astronomers (Cambridge University Press, 2006).
  * [52] Bell, T. J. _et al._ Nightside clouds and disequilibrium chemistry on the hot Jupiter WASP-43b.  _arXiv e-prints_ arXiv:2401.13027 (2024).
  * [53] Bourque, M. _et al._ The exoplanet characterization toolkit (exoctk) (2021).  URL https://doi.org/10.5281/zenodo.4556063.
  * [54] Kreidberg, L.  batman: BAsic Transit Model cAlculatioN in Python.  _PASP_ 127, 1161 (2015).
  * [55] Rubenzahl, R. A. _et al._ The TESS-Keck Survey. IV. A Retrograde, Polar Orbit for the Ultra-low-density, Hot Super-Neptune WASP-107b.  _AJ_ 161, 119 (2021).
  * [56] Carteret, Yann, Bourrier, Vincent & Dethier, William.  The rossiter-mclaughlin effect and exoplanet transits: A delicate association at medium and low spectral resolution.  _A &A_ 683, A63 (2024).  URL https://doi.org/10.1051/0004-6361/202348328.
  * [57] Hut, P.  Tidal evolution in close binary systems.  _A &A_ 99, 126–140 (1981).
  * [58] Peale, S. J.  Origin and Evolution of the Natural Satellites.  _ARA &A_ 37, 533–602 (1999).
  * [59] Fabrycky, D. C., Johnson, E. T. & Goodman, J.  Cassini States with Dissipation: Why Obliquity Tides Cannot Inflate Hot Jupiters.  _ApJ_ 665, 754–766 (2007).
  * [60] Correia, A. C. M. & Laskar, J.  Tidal evolution of exoplanets (2010).  1009.1352.
  * [61] Guillot, T., Burrows, A., Hubbard, W. B., Lunine, J. I. & Saumon, D.  Giant Planets at Small Orbital Distances.  _ApJ_ 459, L35 (1996).
  * [62] Rauscher, E., Cowan, N. B. & Luger, R.  Warm Jupiters Beyond the Tidal Synchronization Limit may Exhibit a Wide Range of Secondary Eclipse Depths.  _AJ_ 165, 261 (2023).
  * [63] Kokori, A. _et al._ ExoClock Project. II. A Large-scale Integrated Study with 180 Updated Exoplanet Ephemerides.  _ApJS_ 258, 40 (2022).
  * [64] Fulton, B. J., Petigura, E. A., Blunt, S. & Sinukoff, E.  RadVel: The Radial Velocity Modeling Toolkit.  _PASP_ 130, 044504 (2018).
  * [65] Piskorz, D. _et al._ Ground- and Space-based Detection of the Thermal Emission Spectrum of the Transiting Hot Jupiter KELT-2Ab.  _AJ_ 156, 133 (2018).
  * [66] Gharib-Nezhad, E. & Line, M. R.  The Influence of H2O Pressure Broadening in High-metallicity Exoplanet Atmospheres.  _ApJ_ 872, 27 (2019).
  * [67] Mansfield, M. _et al._ A unique hot Jupiter spectral sequence with evidence for compositional diversity.  _Nature Astronomy_ 5, 1224–1232 (2021).
  * [68] Iyer, A. R., Line, M. R., Muirhead, P. S., Fortney, J. J. & Gharib-Nezhad, E.  The SPHINX M-dwarf Spectral Grid. I. Benchmarking New Model Atmospheres to Derive Fundamental M-dwarf Properties.  _ApJ_ 944, 41 (2023).
  * [69] Tsai, S.-M. _et al._ VULCAN: An Open-source, Validated Chemical Kinetics Python Code for Exoplanetary Atmospheres.  _ApJS_ 228, 20 (2017).
  * [70] Line, M. R. _et al._ A Systematic Retrieval Analysis of Secondary Eclipse Spectra. I. A Comparison of Atmospheric Retrieval Techniques.  _ApJ_ 775, 137 (2013).