# Multiple Clues for Dayside Aerosols and Temperature Gradients in WASP-69 b
from a Panchromatic JWST Emission Spectrum

Everett Schlawin Steward Observatory, 933 North Cherry Avenue, Tucson, AZ
85721, USA  Sagnick Mukherjee Department of Astronomy and Astrophysics,
University of California, Santa Cruz, CA, USA  Kazumasa Ohno Division of
Science, National Astronomical Observatory of Japan, Tokyo, Japan  Department
of Astronomy and Astrophysics, University of California, Santa Cruz, CA, USA
Taylor Bell Bay Area Environmental Research Institute, NASA’s Ames Research
Center, Moffett Field, CA 94035, USA  Space Science and Astrobiology Division,
NASA’s Ames Research Center, Moffett Field, CA 94035, USA  Thomas G. Beatty
Department of Astronomy, University of Wisconsin–Madison, Madison, WI 53703,
USA  Thomas P. Greene Space Science and Astrobiology Division, NASA’s Ames
Research Center, Moffett Field, CA, USA  Michael Line School of Earth and
Space Exploration, Arizona State University, Tempe, AZ, USA  Ryan C. Challener
Department of Astronomy, Cornell University, 122 Sciences Drive, Ithaca, NY
14853, USA  Vivien Parmentier Université de la Côte d’Azur, Observatoire de la
Côte d’Azur, CNRS, Laboratoire Lagrange, France  Jonathan J. Fortney
Department of Astronomy and Astrophysics, University of California, Santa
Cruz, CA, USA  Emily Rauscher Department of Astronomy, University of Michigan,
1085 S. University Ave., Ann Arbor, MI 48109, USA  Lindsey Wiser School of
Earth and Space Exploration, Arizona State University, Tempe, AZ, USA  Luis
Welbanks School of Earth and Space Exploration, Arizona State University,
Tempe, AZ, USA  Matthew Murphy Steward Observatory, 933 North Cherry Avenue,
Tucson, AZ 85721, USA  Isaac Edelman Bay Area Environmental Research
Institute, NASA’s Ames Research Center, Moffett Field, CA 94035, USA  Natasha
Batalha Space Science and Astrobiology Division, NASA’s Ames Research Center,
Moffett Field, CA, USA  Sarah E. Moran Lunar and Planetary Laboratory,
University of Arizona, Tucson, AZ 85721, USA  Nishil Mehta Université de la
Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange,
France  Marcia Rieke Steward Observatory, 933 North Cherry Avenue, Tucson, AZ
85721, USA

(Accepted May 31, 2024)

###### Abstract

WASP-69 b is a hot, inflated, Saturn-mass planet (0.26 MJup) with a zero-
albedo equilibrium temperature of 963 K. Here, we report the JWST 2 to 12 µm
emission spectrum of the planet consisting of two eclipses observed with
NIRCam grism time series and one eclipse observed with MIRI LRS. The emission
spectrum shows absorption features of water vapor, carbon dioxide and carbon
monoxide, but no strong evidence for methane. WASP-69 b’s emission spectrum is
poorly fit by cloud-free homogeneous models. We find three possible model
scenarios for the planet: 1) a Scattering Model that raises the brightness at
short wavelengths with a free Geometric Albedo parameter 2) a Cloud Layer
model that includes high altitude silicate aerosols to moderate long
wavelength emission and 3) a Two-Region model that includes significant
dayside inhomogeneity and cloud opacity with two different temperature-
pressure profiles. In all cases, aerosols are needed to fit the spectrum of
the planet. The Scattering model requires an unexpectedly high Geometric
Albedo of 0.64. Our atmospheric retrievals indicate inefficient redistribution
of heat and an inhomogeneous dayside distribution, which is tentatively
supported by MIRI LRS broadband eclipse maps that show a central concentration
of brightness. Our more plausible models (2 and 3) retrieve chemical
abundances enriched in heavy elements relative to solar composition by
6$\times$ to 14$\times$ solar and a C/O ratio of 0.65 to 0.94, whereas the
less plausible highly reflective scenario (1) retrieves a slightly lower
metallicity and lower C/O ratio.

stars: atmospheres — stars: individual (WASP-69 b (catalog ))

††facilities: JWST(NIRCam), JWST(MIRI)††software: astropy (Astropy
Collaboration et al., 2013), photutils v0.3 (Bradley et al., 2016), matplotlib
(Hunter, 2007), numpy (van der Walt et al., 2011), scipy (Virtanen et al.,
2020), starry (Luger et al., 2019), batman (Kreidberg, 2015), ThERESA
(Challener & Rauscher, 2022), MC3 (Cubillos et al., 2017), pymc3 (Salvatier et
al., 2016), celerite2 (Foreman-Mackey, 2018), emcee (Foreman-Mackey et al.,
2013), webbpsf (Perrin et al., 2014),

##  1 Introduction

###  1.1 The Value of Spectra of Moderately Irradiated Hot Jupiters

JWST is opening up a new window into giant irradiated planet atmospheres to
reveal their chemical, physical and dynamical processes (e.g. Ahrer et al.,
2023; Alderson et al., 2022; Rustamkulov et al., 2022; Feinstein et al., 2023;
Bean et al., 2023; Grant et al., 2023; Xue et al., 2023). Emission spectra and
phase curves are particularly valuable for understanding the temperature
structure, composition and 3D effects in the atmospheres of hot planets (Bell
et al., 2023a; Coulombe et al., 2023a). Most JWST studies thus far have probed
either cooler, smaller planets (GJ 1214 b Teq=600 K, Kempton et al. (2023);
Gao et al. (2023); LHS 475 b Teq= 590 K, Lustig-Yaeger & Fu et al. (2023); GJ
486 b Teq=700 K, Moran & Stevenson et al. (2023)) or more highly irradiated
Jupiters (WASP-39 b Teq=1170 K (JWST Transiting Exoplanet Community Early
Release Science Team et al., 2023), WASP-43 b Teq =1400 K (Bell et al.,
2023a), HD 149026 b Teq=1700 K (Bean et al., 2023), HD 209458 b Teq=1450 K
(Xue et al., 2023)), WASP-77 A b Teq=1700 K (August et al., 2023). However,
the study of giant planets in the 800-1000 K temperature range is potentially
a valuable window in the chemical equilibrium of the methane (CH4) and ammonia
(NH3) molecules in exoplanet atmospheres, due to both photochemistry and
mixing from deeper layers (Fortney et al., 2020).

Methane and ammonia are expected to be the dominant Carbon-bearing and
Nitrogen-bearing molecules below 900 K to 1000 K for hot Jupiter atmospheres
observed at 0.1 bar (Moses et al., 2013) so their abundances are critical for
evaluating the level of Carbon and Nitrogen in exoplanet atmospheres. Recently
methane has been definitively detected for the first time via low resolution
space-based spectroscopy with JWST, which opens up this molecule as a tracer
of dynamical and chemical processes in exoplanets. Methane was found in the
moderately-irradiated Jupiter, WASP-80 b ($T_{eq}=825K$), (Bell et al., 2023c)
as well as the lower temperature sub-Neptunes K2-18 b T${}_{eq}=270K$
(Madhusudhan et al., 2023) and TOI-270 d T${}_{eq}=354$K (Benneke et al.,
2024; Holmberg & Madhusudhan, 2024). The moderately irradiated Jupiter
HAT-P-18 b (Teq=850 K) may also have methane, as found by retrieval modeling,
but there are no strong, high signal-to-noise absorption features in the SOSS
bandpass from 0.6 µm to 2.8 µm (Fu et al., 2022) and a subsequent re-analysis
did not find evidence for methane in HAT-P-18 b with the same data (Fournier-
Tondreau et al., 2023). Additional spectra of irradiated giant planets in the
800-1000 K temperature regime will be valuable probes of the atmospheric
physics. This temperature regime efficiently covers the transition between
CH4-dominated versus CO-dominated reservoirs of carbon in planet atmospheres.
The presence of CH4 in planets below $\sim$950 K is also predicted to supply
the building blocks for photochemical hazes in planets that can explain some
of the muted atmospheric features below 950 K as observed by HST transmission
spectra (Gao et al., 2020).

###  1.2 WASP-69 b

WASP-69 b is an inflated Saturn mass transiting planet (0.26 $M_{Jup}$, 1.06
$R_{Jup}$) orbiting a K-type star (4715 K, 0.83 M⊙, 0.81 R⊙, Age $\sim$ 2
Gyr), with a zero-albedo full redistribution equilibrium temperature of 963 K
(Anderson et al., 2014). The host star shows evidence for metal enrichment
with metallicities reported from [Fe/H]= 0.15 $\pm$ 0.08 (Anderson et al.,
2014) to 0.29 $\pm$ 0.04 (Sousa et al., 2021). The equilibrium temperature
regime of WASP-69 b is where CH4 can potentially become a dominant carbon-
bearing molecule in chemical equilibrium (Moses et al., 2013), but the cooling
history of a planet and vertical mixing can potentially alter the composition
of the planet from equilibrium expectations (Fortney et al., 2020).

Atoms and molecules have been previously detected in WASP-69 b’s atmosphere
and give initial insights into its composition and aerosols. Searches for Na
in the atmosphere of WASP-69 b with high resolution spectroscopy have yielded
a range of results, but most of the studies find that Na is significantly
detected in the Na D2 line but less significantly detected in the Na D1 line,
potentially from obscuration by atmospheric aerosols (Casasayas-Barris et al.,
2017; Deibert et al., 2019; Langeveld et al., 2022; Khalafinejad et al.,
2021).111 Specifically, Casasayas-Barris et al. (2017) report that the Na D2
line was detected at a 5$\sigma$, while the Na D1 was not detected. Casasayas-
Barris et al. (2017) find that the Na D2 has an excess transit depth of 0.53
$\pm$0.14% for a a 1.5 Å bandpass and a line contrast of 5.8 $\pm$ 0.3% from a
Gaussian fit the transit profile. The slightly different oscillator strengths
of the two lines makes it possible for Na D1 to be more obscured (Huitson et
al., 2012). Deibert et al. (2019), by contrast, do not find significant
absorption by WASP-69 b with the GRACES instrument for either the Na D2 or Na
D1 lines. Langeveld et al. (2022) find that WASP-69 b has a significant Na
line contrast of 3.28 $\pm$ 0.70% for the Na D2 line and 1.26 $\pm$ 0.61 % for
the Na D1 line. Muted but significant detections of both the Na D2 and Na D1
lines were detected by Khalafinejad et al. (2021). The ratio of the two lines
Na D2/ Na D1 of 2.5 $\pm$ 0.7 suggests that there is aerosol opacity that
muted the Na D1 line (Khalafinejad et al., 2021). Thus, the majority of Na
line absorption studies at high spectral resolution mostly suggest an
atmosphere with significant sodium but the Na D1 line is at least partially
obscured by aerosols as compared to Na D2. WASP-69 b’s Na detections have also
been theorized to be potentially influenced by a volcanically active moon (Oza
et al., 2019); so far no moons have been detected (Narang et al., 2023).

Other atmospheric features of WASP-69 b found with transmission spectroscopy
include its scattering slope, tentative evidence for TiO, Carbon-bearing and
Oxygen-bearing molecules and an escaping outflow, all described below. Murgas
et al. (2020) find a spectral slope in the GTC OSIRIS optical transmission
spectrum that is consistent with Rayleigh scattering from hydrogen. However,
the spectral slope is also consistent with stellar activity and the transit
light source effect (e.g. Rackham et al., 2018). Estrela et al. (2021) find an
even stronger spectral slope than Murgas et al. (2020) with HST’s Space
Telescope Imaging Spectrograph (STIS). The large transit depths down to
$\sim$0.4 µm are best explained by aerosols at high altitudes ($\mu$bar
pressure levels) (Estrela et al., 2021). There is further evidence of Raleigh
scattering in ground-based spectroscopy from 0.5 to 0.9 µm with tentative
evidence of TiO in the spectrum (Ouyang et al., 2023). Extended Helium has
also been detected in WASP-69 b with a signal-to-noise ratio of 18 that is
blue-shifted by several km/s (Nortmann et al., 2018). This was confirmed with
narrowband photometry, constraining the mass loss rate to be 10-4 to 10-3
MJup/Gyr (Vissapragada et al., 2020). The envelope escaping WASP-69 b extends
several radii and is confined to a tail like structure that extends up to 7.5
planetary radii behind the planet’s direction of motion (Tyler et al., 2023).
However, the inferred mass loss rate from the observed escaping tail ($\sim
10^{-3}$ MJup/Gyr) is not a significant fraction of the overall mass (0.26
MJup) to measurably alter the deeper atmosphere.

In transit transmission, 5 molecules were detected with high resolution cross
correlation (Guilluy et al., 2022): CO, CH4, NH3, H2O and C2H2. In the same
analysis, the detection of H2O and CH4 and the absence of CO2 disfavored
models larger than 10$\times$ solar enrichment of heavy elements. If the
atmosphere is assumed to be solar composition, the presence of the 5 detected
molecules favors models that have high C/O ratios, perhaps larger than 1.0.
The presence of C2H2 in the terminator also implies some disequilibrium
processes.

Baxter et al. (2020) and Wallack et al. (2019) used secondary eclipse
measurements at the Spitzer 3.6 µm and 4.5 µm bandpasses to constrain the
composition and temperature profile of WASP-69 b. Baxter et al. (2020) found
that WASP-69 b and other planets below $\lesssim 1660$K lack temperature
inversions, ie. their gas temperatures decrease as a function of altitude. The
spectral slope of the Spitzer 3.6 µm and 4.5 µm eclipse depths was also fit
with chemically equilibrated atmospheric models. The low spectral slope
observed in the planet ($\lesssim 20ppm/\micron$) favors a high atmospheric
metallicity of the planet above $\sim 30\times$ solar (Wallack et al., 2019).
However, the two photometric eclipse depths cannot uniquely constrain the
composition of WASP-69 b due to the correlation with the planet’s temperature-
pressure profile, aerosols, abundance ratios and possible disequilibrium
effects (Wallack et al., 2019). JWST observations covering a wide range of
wavelengths that encompass many molecular features and probe different
pressure levels can better constrain the abundances, dynamics, chemical
processes and aerosols in WASP-69 b and similar planets’ atmospheres.

####  1.2.1 Paper Outline

In Section 2, we describe the observations of 3 eclipses spanning 2 µm to 12
µm. In Section 3, we describe three independent analyses to extract the
spectrum of WASP-69 b from the raw data and also describe the absolute flux of
the star as compared to a stellar model. In Section 4, we present the planet’s
dayside spectrum and its overall properties. We introduce several types of
models, the results of atmospheric retrievals and why they must incorporate
aerosols and 3D effects in Section 5. We find some evidence for an
inhomogeneous dayside temperature distribution, which we investigate
independently with broadband eclipse mapping in Section 6. We discuss the
energy budget of the planet’s dayside and the properties of atmospheric
aerosols and synthesize WASP-69 b’s inferred composition across all of our
models in Section 7. Finally, we conclude in Section 8 that WASP-69 b is
chemically enriched relative to solar composition, has dayside aerosols and
temperature gradients away from the substellar point.

##  2 Observations

Figure 1: Dynamic Spectra for WASP-69 b’s show the eclipse as a horizontal
dark band. The eclipse depth (as seen by the darkness of the band) increases
from short wavelengths toward long wavelengths and also shows a dip where
there is strong CO2 atmospheric absorption near 4.3 µm. The dynamic spectra
are the extracted spectra for each detector column (NIRCam) and row (MIRI) at
full spectral and full time resolution without binning. The top panels of each
dynamic spectrum show the average spectrum in electrons for one integration.

WASP-69 b was observed as part of the MIRI And NIRCam Assay for the
Transmission and Emission of Exoplanets (MANATEE), which combines observations
from GTO program 1177 (Observation 3) and GTO program 1185 (Observation 6 and
7) to build up near-infrared to mid-infrared spectra on a wide variety of
planets. Table 1 contains a summary of observations.

For NIRCam, we used the readout pattern with the most number of reads per
integration to mitigate the effects of 1/f noise (Schlawin et al., 2020) while
still keeping data volume manageable222With BRIGHT2, our NIRCam observations
had data volume excess above the sustainable rate (0.87 MB/s) with an excess
of 7 (9) GB for the F322W (F444W) filters, respectively, but below the medium
volume excess threshold of 15 GB. With the usually-favored RAPID readout
pattern (1 frame per group), the observations would exceed the medium
threshold with an excess of 24 GB., which led us to the BRIGHT2 readout
pattern, which has no skipped frames and 2 co-added frames per group. We
selected the F210M filter and WLP8 pupil with the short wavelength arm of
NIRCam to measure the the eclipse at 2.1 $\mu$m for both NIRCam observations.
For MIRI, we used the FASTR1 read pattern with the slit-less Low Resolution
Spectrometer (LRS) mode. The 3 separate observations spanned from 372 to 373
minutes in duration, or 2.8 times the eclipse duration for the planet. This
long baseline ensured significant detector settling time for persistence
effects like charge trapping (Zhou et al., 2017) and MIRI detector upward and
downward ramps (Bell et al., 2023b; Bouwman et al., 2023) to settle. The
observation details were specified before Cycle 1 measurements of these
settling constants, but the on-orbit measurements of exoplanet lightcurves
show that the time constants for the charge trapping ramp-up are short for the
near-infrared instruments that employ HgCdTe H2RG detectors ($\lesssim$ 15
min) (e.g. Ahrer et al., 2023; Feinstein et al., 2023; Espinoza et al., 2022;
Schlawin et al., 2023b), so future near-infrared observations may succeed with
shorter pre-transit or pre-eclipse baselines.

The pixel-level lightcurves for the three observations are shown in Figure 1
along with the detected counts in electrons per integration. These dynamic
spectra are normalized by the median spectrum per integration to show the
change as a function of time and wavelength. The eclipse as the planet goes
behind the star is clearly detected in all three observations at this pixel-
level resolution, as seen by the dark horizontal band across the dynamic
spectra. As expected, the eclipse depth rises from the shortest to longest
wavelengths because the planet has a cooler temperature than the star. Also
visible in the dynamic spectra is the shallower eclipse depth near 4.3 µm due
to carbon dioxide in the planet’s atmosphere.

Some systematic correlated noise is also visible in the dynamic spectra in
both the NIRCam and MIRI data. For NIRCam, horizontal stripes and bands are
visible due to the 1/f noise present in the detector, which arises from the
electronics that read the detector in the horizontal (wavelength) direction
for NIRCam (Schlawin et al., 2020; Ahrer et al., 2023). 1/f noise appears as
horizontal striping in the individual detector groups that manifests as
horizontal striping in the dynamic spectra for the two NIRCam observations.
The 1/f noise was already mitigated by row-by-row subtraction using background
pixels (e.g. Schlawin et al., 2023b), as described in Section 3. However, the
row-by-row subtraction is only possible with background pixels of the right-
most amplifier for the F322W2 NIRCam observation and the left two amplifiers
of the F444W NIRCam observation, when they are oriented in the Data Management
System orientation that is used in MAST. 1/f noise that is present in the
middle two amplifiers can be mitigated by subtracting the median of rows of
the other amplifiers, but this can only subtract 1/f noise that is common to
the shared SIDECAR ASIC device and not noise that is unique to each
amplifier’s p-type field-effect transistor (PFET) (Schlawin et al., 2020).
Therefore, residuals 1/f noise is visible in Figure 1, mostly at wavelengths
shorter than 3.89 µm for the the F322W2 observation 6 and longer than 4.09 µm
for observation 7. Also noticeable in the dynamic spectra is a downward flux
trend, most prominently in the NIRCam F322W2 spectrum near 3.1 µm. As will be
described in Section 3, we further mitigate and account for 1/f noise with co-
variance weighted extraction, fitting excess noise and model the time-
dependent trend with either a polynomial or the detector housing temperature.

For MIRI, there are also two systematic effects visible in Figure 1. There is
a ramp visible in the first several hundred integrations which is typical of
MIRI/LRS observations, especially at the shortest wavelengths (e.g., Bell et
al., 2023b, 2024; Bouwman et al., 2023). There is also an odd-even effect in
alternating spectral-axis pixels that appears as vertical striping in the
dynamic spectrum; this is caused by the simultaneous resetting of pairs of
rows (in MIRI’s 90-degree rotated reference frame), and such odd-even striping
is also typical of MIRI/LRS observations (e.g., Bell et al., 2024). This odd-
even effect, the 390 Hz noise seen in the data from some MIRI subarrays
(Bouwman et al. 2023; Bell et al. 2024), and other still-unknown noise sources
result in excessively noisy MIRI/LRS spectra at high-resolutions which is
significantly improved by spectral binning (Bell et al., 2024; Welbanks & et
al., 2023). As a result, we used wider bin sizes (0.25 $\mu$m full width) to
average over these effects and give more reliable uncertainty estimates.

Table 1: Summary of WASP-69 b observations UT Exp Date | Prog | Obs # | Instrument | SW Filter | Spectroscopic | Spec Wave Range | Nint | Ngroup | Readout
---|---|---|---|---|---|---|---|---|---
YYYY-mm-dd |  |  |  |  | Description | µm |  |  | Pattern
2022-10-25 | 1185 | 6 | NIRCam | F210M | F322W2 | 2.4 - 4.01 | 2367 | 3 | BRIGHT2
2023-06-06 | 1185 | 7 | NIRCam | F210M | F444W | 3.88 - 4.98 | 1507 | 5 | BRIGHT2
2023-05-02 | 1177 | 3 | MIRI |  $\cdots$ | LRS | 5.0 - 14 | 8276 | 16 | FASTR1
Table 2: Summary of WASP-69 b Properties and Priors

Parameter | Prior/Existing Value | Prior/Existing Source | TESS Posterior (this work)
---|---|---|---
t0 | 2455748.83344 $\pm$ 1.8e-4 | Casasayas-Barris et al. (2017) | 2455748.83345 $\pm$ 1.8e-4
P | 3.8681382 $\pm$ 1.7e-6 d | Casasayas-Barris et al. (2017) | 3.86813942 $\pm$ 1.8e-7 d
$i$ | 86.71∘ $\pm$ 0.20 | Casasayas-Barris et al. (2017) | 86.631∘ $\pm$ 0.029
a/R∗ | 12.0 $\pm$ 0.46 | Casasayas-Barris et al. (2017) | 11.919 $\pm$ 0.05
Rp/R∗ | 0.1336 $\pm$ 0.005 | Casasayas-Barris et al. (2017) | 0.1285 $\pm$ 0.00025
e | 0.0 | fixed | fixed
R∗ | 0.813 $\pm$ 0.028 R⊙ | Anderson et al. (2014) |  $\cdots$
log g∗ (cgs) | 4.535 $\pm$ 0.023 | Anderson et al. (2014) |  $\cdots$
Distance | 50.29 $\pm$ 0.04 pc | Gaia Collaboration et al. (2023) |  $\cdots$

Note. — For the spectroscopic analysis, we use the prior values from
Casasayas-Barris et al. (2017), which come from Anderson et al. (2014). We
also fit the publicly available TESS data and report the posteriors (right) to
improve the orbital parameter precisions for eclipse mapping.

##  3 Data Analysis and Extraction

###  3.1 Reduction and Lightcurve Extraction

####  3.1.1 tshirt

We extract spectra in a similar manner as in previous JWST data analyzed with
tshirt and published (Ahrer et al., 2023; Bell et al., 2023c; Welbanks & et
al., 2023), but with the exception that we allow for the small curvature and
tilt of the trace (about 4 pixels across the spectrum with NIRCam, but
negligibly for MIRI LRS). For NIRCam data, we begin with the _uncal.fits stage
0 data products from MAST and process them with a modified version of the jwst
pipeline (Bushouse et al., 2023) version 1.10.2, except for the F210M
photometry on Observation 6, which was extracted soon after observation with
an older jwst version 1.6.0. We ran the jwst pipeline with CRDS context
jwst_1077.pmap for Program 1185 observation 6’s spectroscopy and CRDS context
jwst_1093.pmap for observation 7’s spectroscopy, but checked that the used
reference files were the same for both observations’ contexts. For the
photometry, we used jwst_1009.pmap for both observations. We also used tshirt
to extract the spectrum of the MIRI LRS data in Program 1177 observation 3.
For MIRI LRS data, we used SDP version 2023_4a, jwst version 1.13.4, CRDS
version 11.17.15 and CRDS context jwst_1225.pmap. For reference on different
version numbers, we have found negligible differences in extracted tshirt
photometry and spectroscopy $\lesssim$ 10 ppm depths between JWST version
numbers and reference file contexts from JWST versions 1.8 to 1.13 and pmap
contexts 1039 to 1188.

Our main modification to the jwst pipeline for NIRCam is to replace the
reference pixel correction with a row-by-row, odd-even by amplifier (ROEBA)
correction to mitigate 1/f noise in the photometry and spectroscopy (e.g.
Schlawin et al., 2023b). For MIRI LRS data, we do not use the ROEBA
correction. We also skip the dark current step for NIRCam because it adds
noise for subarray modes. We also set the jump sigma rejection threshold for
cosmic rays to 6 $\sigma$ for NIRCam and 7 $\sigma$ for MIRI and stop
processing the jwst pipeline after the _rateints.fits products. We manually
divide by the NIRCam spectral images by the flat field reference file for
imaging in CRDS (jwst_nircam_flat_0313.fits for the F444W filter and
jwst_nircam_flat_0266.fits for the F322W2 filter). The jwst pipeline was run
on NIRCam data with the JWST Time Series Observation Wrapper (jtow) version
0.1.4 (https://github.com/eas342/jtow.git) and the spectra and photometry were
extracted with the Time Series Helper and Integration Reduction Tool tshirt
version 0.3 (https://github.com/eas342/tshirt). For MIRI LRS data, we used an
updated tshirt version 0.4 that includes wavelength calibration of MIRI LRS.
We do not apply a flat field correction to the MIRI LRS data.

For the photometric F210M data, we use a circular aperture of 92 pixels and a
background annulus from 93 to 129 pixels, as shown in Appendix C. For
spectroscopy, we fit the traces of the spectrum for one image and use this to
define the aperture and background for all images. For the trace centroids, we
fit the profile with a Gaussian in each column and then fit the centroids with
a 3rd order polynomial with iterative 3 $\sigma$ clipping of outlier points.
For the long wavelength extraction, we next subtract the background using a
linear fit along the Y (spatial) direction for all the pixels that are more
than 7 pixels away from the source (rounded to the nearest whole pixel) and
also located in a rectangular region from pixels Y=5 to 65. We use a
covariance-weighted extraction (Schlawin et al., 2020) for the non-NaN pixels
that are within 5 pixels from the source trace, assuming a read noise of 14
$e^{-}$ and a correlation of 0.08 between pixels. We use the interpolated PSF
to estimate the missing flux in the spatial direction for pixels that are
marked as NaN or outliers more than 30$\sigma$ from a spline fit and multiply
the flux to account for the missing pixels’ fractional flux. The extracted
average spectrum in ($e^{-}$) and spectroscopic lightcurves at the pixel
resolution are shown in Figure 1. We next bin the spectra in wavelength for
lightcurve analysis to approximately 0.010 µm wide wavelength bins, but
rounded to the nearest whole pixel.

####  3.1.2 Eureka!

Our Eureka! NIRCam and MIRI reductions used version 0.10 of the Eureka!
pipeline (Bell et al., 2022); jwst package version 1.10.2 (Bushouse et al.,
2023); CRDS version 11.17.0; and CRDS context 1094 for NIRCam/F322W2 and
NIRCam/F444W, and 1097 for MIRI/LRS. Our NIRCam reduction methods closely
follow those used in previous Eureka! NIRCam spectroscopy analyses (Ahrer et
al., 2023; Bell et al., 2023c; Welbanks & et al., 2023), and our MIRI/LRS
reduction method generally followed the Eureka! v1 method described by Bell et
al. (2023a) and the Eureka! MIRI/LRS analysis of Welbanks & et al. (2023).
Eureka! Control Files and Eureka! Parameter Files that can be used to
reproduce this work are available for download
(https://doi.org/10.5281/zenodo.11168833; Schlawin et al. 2024), but the
important parameters are summarized below.

Eureka!’s Stages 1 and 2 make use of the jwst pipeline’s (Bushouse et al.,
2023) Stages 1 and 2 but allow for some modifications. For both NIRCam
spectroscopy datasets, the only change to Stage 1 was increasing the jump
rejection threshold to 6 to avoid excessive false positives when ramp fitting.
For MIRI/LRS, we similarly increased the jump rejection threshold to 7 and
also turned on the firstframe and lasstframe steps (which skip the first and
last frames of each integration) as these frames are subject to excessive
noise and generally increase the final noise level in the data. Finally,
during the Stage 1 processing of the MIRI/LRS data we also applied Eureka!’s
newly developed 390 Hz noise removal and group-level background subtraction
technique described in detail by Welbanks & et al. (2023); in general, this
step removes the wavelength correlated noise in raw MIRI/LRS SLITLESSPRISM
subarray data first pointed out by Bouwman et al. (2023). Our Stage 2
processing followed the jwst standard processing with the exception of turning
off the photom and extract1d steps which are not required for time-series
observations. Later, we re-analyzed the data for checking the stellar absolute
flux against a model.

Our Stage 3 processing of the NIRCam data closely followed the Eureka! NIRCam
reductions of Welbanks & et al. (2023). In summary, we cropped the frames to
only include the relevant pixels, masked outlier pixels, corrected for the
curvature of the spectral trace, performed a linear column-by-column
background subtraction using pixels $>$13 px away from the source center,
computed the position and width of the source for later use when fitting the
lightcurves, and performed optimal spectral extraction (Horne, 1986)
considering only the pixels within 5 pixels of the source center and using a
cleaned median integration to compute our variance-weighted spatial profile.
Our Stage 3 processing of the MIRI data also closely followed the Eureka! MIRI
reduction of Welbanks & et al. (2023) which differs only slightly from the
NIRCam reduction steps. In particular, the differences from the NIRCam
reduction steps include rotating the MIRI/LRS data to have wavelength
increasing toward the right, manually applying an estimated gain of 3.1 e/DN
(Bell et al., 2023b, a), using pixel indices 11–61 (excluding pixels within 9
px of the spectral trace) to perform a column-by-column background
subtraction, and only using pixels within 3 px of the source center when
performing optimal spectral extraction.

In Stage 4, we spectrally binned the data and sigma-clipped temporal outliers.
For the NIRCam data we used 0.01 $\mu$m spectral bins across 2.45–3.95 $\mu$m
for F322W2 and 3.89–4.97 $\mu$m for F444W, while for the MIRI data we used
coarser 0.25 $\mu$m spectral bins from 5–12 $\mu$m following the
recommendation of Bell et al. (2023a). We did not see evidence of the
“shadowed region effect” described by Bell et al. (2023a), and therefore
include the 10.6–11.8 $\mu$m data in our reduction. For all three
observations, we performed 4-sigma clipping comparing each point to a smoothed
version of the data computed with a 20-integration wide boxcar filter which
helped to removed otherwise missed temporal outliers while ensuring not to
mask the eclipse ingress or egress.

####  3.1.3 Pegasus

We also reduced the NIRCam observations using the Pegasus pipeline
(https://github.com/TGBeatty/PegasusProject). Our Pegasus reduction began with
background subtraction on the rateint files provided by version 1.10.2 of the
jwst pipeline using CRDS version 11.17.0. We began by applying a basic
background subtraction to each rateint file by fitting a two-dimensional,
second-order spline to each integration, covering the full 256$\times$2048
rateint images. For the spline fitting, we masked out image rows 5 through 75
to prevent the self-subtraction of light from the WASP-69 system. We then
performed a single round of $3\,\sigma$ clipping on the unmasked portions of
the image. Next, we fit individual splines to each amplifier region on the
image using the unmasked, unclipped pixel values with a median box size of 20
pixels. This per-amplifier spline fitting was necessary to eliminate residual
bias differences across the amplifier areas. We extrapolated the combined
background spline for the whole image over the masked portions near WASP-69
and subtracted it from the original image values. Visual inspection of the
rateint images showed that in roughly 2% 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). Finally, we
attempted to remove some of the red-noise caused by the NIRCam readout
electronics which is present along detector rows, by calculating the robust
mean of each row using pixels from column 1800 onwards for the F322W2 images
and up to column 600 for F444W (both chosen to avoid light from WASP-69), and
then subtracting this mean from each row. Visually, this removed most of the
horizontal banding typical for NIRCam grismr images.

We then extracted spectroscopic lightcurves from our background-subtracted
images. To do so, we fit the spectral trace using a fourth-order polynomial
and then used optimal extraction to measure the 1D spectrum in each image. We
performed three rounds of iterative profile estimation for the optimal
extraction routine, after which we judged the profile fit to have converged.
Using the resultant 1D spectra, we extracted a broadband lightcurve from 2.45
$\mu$m to 3.95 $\mu$m at F322W2 and from 3.89 $\mu$m to 4.97 $\mu$m at F444W.
For the spectroscopic lightcurves we subdivided each of these wavelength
regions into individual 0.01µm-wide spectral channels. During this extraction
process, we linearly interpolated over each spectral column to account for
partial-pixel effects in the 0.01µm wavelength bins.

###  3.2 Lightcurve Fitting

####  3.2.1 tshirt

We fit the lightcurves with a starry (Luger et al., 2019) lightcurve model,
assuming a uniform intensity map characterized by an amplitude term only. We
use the probabilistic programming suite pymc3 (Salvatier et al., 2016) and the
pymc3-ext tools from Foreman-Mackey et al. (2021) to calculate the posterior
distributions of the model with No U-Turns sampling. For the orbital parameter
priors, we use the values from Casasayas-Barris et al. (2017) listed in Table
2. For NIRCam, we bin the time series to a resolution of 300 equally spaced
time bins at a cadence of 1.24 minutes to speed up model evaluation and also
provide an empirical estimate of noise. We estimate the errors by taking the
standard error in the man of all integrations within each 1.24 minute bin and
then taking the median of all bins’ standard deviations and adopting this
across all time bins. edit1For MIRI data, we trim the first 1000 integrations
(ie. 45 minutes) to discard much of the initial ramp-up behavior (Bell et al.,
2023b) and bin the data to 100 equally spaced time bins.

We have found that the measured flux of many different targets anti-correlates
with NIRCam’s detector housing temperature. Figure 18 in Section C shows the
Focal Plane Housing Temperature and the out-of-transit and out-of-eclipse
fluxes of many different exoplanet observations and how they anti-correlate.
While the mechanism is not exactly understood, we can model the flux as a
linear function of the focal plane housing temperature deviation. For both the
short wavelength photometry and the long wavelength grism spectroscopy, we
collect the detector housing temperature from the MAST engineering database.
The short wavelength and long wavelength temperatures are accessible by
querying the IGDP_NRC_A_T_SWFPAH1 and IGDP_NRC_A_T_LWFPAH1 telemetry
mnemonics, respectively. We find the temperature deviation by subtracting all
temperatures by the median. The temperature sensors have significant noise, so
we smooth both the short wavelength and long wavelength housing temperatures
first by fitting them with a 5th order polynomial and using this as an
interpolation function for the temperature at each integration. Finally, we
model the NIRCam flux as a linear function of the smoothed temperature
deviation.

In addition to the drift with focal plane housing temperature for NIRCam, we
include a linear slope in time to account for other drifts such as stellar
rotation. We assume a prior slope of 0$\pm$1.6%/hr for the linear drift. For
the MIRI data, we also include an exponential baseline function to fit for the
ramp behaviors of MIRI LRS lightcurves (Bell et al., 2023b). For the
exponential ramp, we assume the behavior decays exponentially from the first
integration with a timescale prior with a lognormal prior that has a timescale
of 1 minute with a wide geometric mean of 2 and amplitude of 0.1%.
Additionally, a strong level of correlated noise is visible in both NIRCam
observations using the same F210M photomeric filter, as seen in Figure 2.
Therefore, for the short wavelength F210M lightcurves only, we also model the
time series with a Gaussian Process (e.g. Gibson et al., 2012). We use
celerite2 code (Foreman-Mackey, 2018) with a stochastically-driven damped
harmonic oscillator kernel that has a fixed quality factor of Q=0.25 and a
log-normal prior on the period of the oscillator $\rho$ that has a geometric
mean equal to the duration of the observation.

We iteratively sigma clipped any fluxes that deviated by the model by more
than 5 $\sigma$ and iterated 2 times to maximize the a priori probability
using the starry, pymc3 and pymc3-ext software suites. The resulting maximum a
priori lightcurve solutions are shown in Figure 2. This maximum a prior
solution was used as an initial set of parameters for the Hamiltonian No
UTurns Sampler (NUTS). The sampling was tuned with 3000 steps and samples for
another 3000 steps with two different chains.

Figure 2: Top Half-Panels: Broadband eclipse Lightcurves with a Best-Fit
Model. Bottom Half-Panels: Residuals of the model.

We first found the posterior distribution of the broadband eclipse lightcurve.
We then fix the semi-major axis, inclination, period, transit time (linked to
the eclipse time) at the broadband fit value for all spectroscopic lightcurve
fits. We allow the linear baseline and FPAH temperature deviation terms to be
free parameters with the same priors as applied to the broadband.

We note that the short wavelength F210M lightcurves exhibited significant
correlated noise, as seen in Figure 2, upper left, especially for Observation
7. Ignoring this noise and fitting with polynomial de-trending vectors results
in significantly discrepant eclipse depths by 3.6$\sigma$. However, when we
use a GP fit to both observations, we find agreement to within 0.7$\sigma$.
The correlated noise seen in these observations is larger than found for
HAT-P-14 b with the same defocused photometry mode (Schlawin et al., 2023b).
The increased noise for WASP-69 b may be related to activity in the host star,
which has been found to have spots and a rotation period of $\sim$24 days
(Khalafinejad et al., 2021; Anderson et al., 2014). Another possibility is
JWST primary mirror flexures because the defocused photometry is highly
sensitive to JWST’s wavefront (e.g. Schlawin et al., 2023b; McElwain et al.,
2023). The longer wavelength spectroscopic observations (grism time series and
MIRI LRS) show less of this correlated noise from either of these effects
because stellar spot contrast decreases as function of wavelength and these
modes are not as sensitive to JWST’s wavefront.

The resulting NIRCam emission spectrum of WASP-69 b is shown in Figure 3.
There is a rise from 2.5 to 5.0 µm expected for a $\sim$963 K zero-albedo
equilibrium temperature planet orbiting a 4715 K star. Additionally, there is
a deep CO2 absorption feature visible at 4.3 µm, but no obvious CH4 feature at
3.3 µm. A comparison with the Eureka! and Pegasus reductions (described below)
is also shown in Figure 3, which show consistency within 1 $\sigma$ errors,
but some subtle 30-50 ppm offsets between reduction methods.

Figure 3: Left: Eclipse spectrum with three independent extraction and
lightcurve fitting methods on a common $\sim$0.010 µm-spaced grid. All spectra
agree within the $\sim$50-80 ppm errorbars. Right: Eclipse Spectrum with MIRI
with two independent extraction and lightcurve fitting methods .

####  3.2.2 Eureka!

As was done with tshirt, our Eureka! lightcurve fitting also used a starry
(Luger et al., 2019) astrophysical model assuming a uniform map. For the
orbital parameter priors, we fixed all parameters except $t_{0}$ to the values
from Casasayas-Barris et al. (2017). For $t_{0}$, we used a Normal prior based
on the posterior of Casasayas-Barris et al. (2017) when fitting the MIRI/LRS
broadband (5–12 $\mu$m) lightcurve which resulted in an inferred
$t_{0}=2455748.83473\pm 0.00012$; we then adopted this new value and fixed
$t_{0}$ for the spectroscopic fits to all lightcurves. We also assumed a
stellar radius of 0.813 $R_{\odot}$ (Casasayas-Barris et al., 2017) combined
with $a/R_{*}$ to account for the light travel delay when computing the
expected time of eclipse. For our systematic noise model, we used a polynomial
in time (quadratic for NIRCam and linear for MIRI), a linear decorrelation
against the spatial position and PSF-width of the spectral trace, and for MIRI
an exponential ramp following the recommendations of Bell et al. (2023a). We
also removed the first 800 integrations from the MIRI observations to reduce
the impact of the strong initial exponential ramp. We also fitted an error
inflation parameter which multiples our estimated uncertainties by a constant
factor; this ensures a reduced chi-squared of 1.0 and avoids issues caused by
incorrectly estimated gain values in Stage 3. For all our systematic noise
models, we used minimally informative priors.

We used pymc3’s No U-Turns Sampler (Salvatier et al., 2016) to sample our
posteriors using two independent chains with a target acceptance rate of 0.85.
For the MIRI spectra we used 3000 tuning steps and 1500 posterior samples,
while for the NIRCam spectra we used 4000 tuning steps and 3000 posterior
samples. For all fits, the Gelman-Rubin statistic (Gelman & Rubin, 1992) was
at or below 1.01, ensuring the chains had converged. We then used the 16th,
50th, and 84th percentiles of the posterior samples to estimate the best-fit
parameter values and their corresponding uncertainties. The fits to our NIRCam
data showed minimal residual red noise, and the white noise in our residuals
was $\sim$20–30% above the estimated photon limit for F322W2 and $\sim$10–20%
above the estimated limit for F444W. Meanwhile, the MIRI data did exhibit
significant residual red noise in some channels (especially between 6–8
$\mu$m), and for MIRI the white noise level in our residuals ranged from
$\sim$20–30% above the estimated photon limit from 6–10 $\mu$m with a gradual
increase to 80% above the limit longward of 10 $\mu$m and a steep increase to
140% above the limit shortward of 6 $\mu$m (the cause of which could not be
identified). To account for the impact of unmodelled red noise on our final
spectra, we used the $\beta$ error inflation method developed by Winn et al.
(2008); in particular, we computed $\beta$ at timescales spanning 22–32
minutes (within 5 minutes of WASP-69b’s $\sim$27 minute ingress/egress
duration) and multiplied our emission spectra uncertainties by the computed
$\beta$ value. For most spectroscopic channels this inflated the uncertainties
by $<$20% but significantly increased the uncertainties of three channels
(6.125, 6.875, and 7.75 $\mu$m).

####  3.2.3 Pegasus

Our Pegasus analysis used a BATMAN eclipse model (Kreidberg, 2015) to fit the
spectroscopic lightcurves extracted via the Pegasus pipeline. We fixed the
orbital parameters of WASP-69b to those measured in Casasayas-Barris et al.
(2017), which left the free parameters in our spectroscopic lightcurve fitting
to be the secondary eclipse depth and the slope and normalization of a
background linear trend. We did not impose a prior on any of these parameters.
We fit each spectral channel individually.

To fit the spectroscopic lightcurves we performed an initial likelihood
maximization using a Nelder-Mead sampler followed by MCMC likelihood sampling.
We used the maximum likelihood point identified by the Nelder-Mead
maximization as the starting locus about which we initialized the MCMC chains.
To perform the MCMC runs, we used the emcee Python package (Foreman-Mackey et
al., 2013) using 12 walkers with a 2,000-step burn-in and then a 4,000-step
production run for each spectral channel. We checked that the MCMC had
converged by verifying that the Gelman-Rubin statistic was below 1.1 for each
parameter in each spectral channel.

We additionally checked the goodness-of-fit and statistical properties of our
eclipse modeling in each spectral channel. We did so by first verifying that
the average of the per-point flux uncertainties in each channel’s lightcurve
matched the standard deviation of the residuals to the best fit eclipse model.
We also computed the Anderson-Darling statistic for each channel’s lightcurve
residuals to check that the residuals themselves appeared Gaussian. We did not
find statistically significant non-Gaussianity in the residuals to our
spectroscopic lightcurve fits.

###  3.3 Comparison of Extracted Spectra

A comparison of our three NIRCam reductions on a common $\sim$0.010µm
wavelength grid is visible in Figure 3. The tshirt, Eureka! and Pegasus
reductions all agree well within the error bars of 50-80 ppm. Even though each
spectral extraction treated 1/f noise and there were differences in treatment
of the systematic trends, the three reductions give very similar absorption
features, slopes and overall depths. A moving average of the observations
showed agreement to better than the 50 ppm level. Given the close agreement
between reductions, we proceed with the tshirt reduction, which had the most
1/f noise mitigations and smallest broadband out-of-eclipse lightcurve
scatter. We also did two independent reductions of the MIRI LRS spectrum shown
in Figure 3. Again, the agreement for both reductions is broadly consistent
with only 2 points differing by 2.6 and 2.8 times the 1$\sigma$ uncertainties
of the Eureka! reduction. These two wavelengths do not correspond to any
strong opacity sources in our models. We use the Eureka! data for our modeling
because it has a longer heritage of analysis with MIRI LRS data (e.g. Bell et
al., 2024; Kempton et al., 2023).

###  3.4 Absolutely Calibrated Stellar Spectrum

Figure 4: Top: Calibrated absolute flux for 3 observations of the host star
WASP-69 A and a 4750 K BOSZ model (Bohlin et al., 2017). Additionally,
photometry from 2MASS KS and WISE bands 1-3 is shown as data points with error
bars. The black circles are synthetic photometry of the BOSZ model for those
same photometric bands. Bottom: The ratio of the data over the model shows
closest agreement in the F322W2 filter (where a direct ratio with a calibrator
was possible) and some 1-3% discrepancies where stellar CO absorption is
present.. No circumstellar debris disks are present above $\gtrsim 3\%$ flux
levels that could dilute the planet’s flux significantly.

We also create an absolutely calibrated stellar spectrum across the NIRCam and
MIRI wavelengths to find the intrinsic dayside planet flux as well as verify
if a stellar model provides an accurate estimate for the star. We first
calculated calibrated stellar spectra and compared them to stellar models that
most closely matched the stellar parameters from Anderson et al. (2014):
T∗,eff = 4700 $\pm$ 50 K, log(g∗)=4.5$\pm$0.15, [Fe/H] = 0.15 $\pm$ 0.08,
R=0.813 R⊙ and a Gaia DR3 distance of 50.29 $\pm$ 0.04 pc (Gaia Collaboration
et al., 2023). For the BOSZ model (Bohlin et al., 2017), the parameters are
T=4750 K, log(g)=4.5 and [Fe/H]=0.0.

For the NIRCam observations, we extracted stellar spectra of the solar-analog
calibrator GSPC P-330E with the same extraction parameters as WASP-69 b to
minimize errors due to aperture corrections and/or field-dependence
differences with the NIRCam Wide Field Slitless calibration observations. We
used calibration program 1076 observations 1 and 2 for the F322W2 and F444W
filters, respectively. We used the CALSPEC model for GSPC P-330E
(p330e_mod_006.fits) to convert from DN/s to physical flux units by dividing
the stellar model by the observed count rates.

Finally, we multiplied this calibration factor by the observed count rate to
calculate the calibrated NIRCam stellar spectrum of WASP-69 b. Program 1076
observation 2 was take at a different field point that was separated by 94
pixels in the detector X direction from the science observations, so it
necessitated interpolating the F444W calibrator to the same wavelengths, which
was a shift of about 0.092 µm. An improved model for P-330E is being prepared
(Rieke & et al., 2023) that modifies the CALSPEC model by up to 1.5% at 2 µm
but this is below the level of systematic errors, such as non-linearities in
our measured absolute fluxes. We also compared the F322W2 calibration to the
value when extracting the spectrum with the jwst pipeline with the photometry
step applied, an imaging flat field jwst_nircam_flat_0610.fits, aperture
correction jwst_nircam_apcorr_0003.fits, and pixel area
jwst_nircam_area_0054.fits from CRDS and find agreement to within 3% with the
jwst pipeline correction being 3% higher than the method with a ratio to GSPC
P-330E.

For the MIRI observations, we re-ran Eureka!’s Stage 2 with the photometric
calibration step turned on and then used Eureka!’s Stage 3 optimal extraction
on the _cal.fits calibrated data product using the exact same procedure we
used to produce the MIRI lightcurve. The Eureka! optimal extraction allowed
for better interpolation over bad pixels than a simple box extraction.
Additionally, we analyzed the MIRI/LRS SLITLESS calibration observations of HD
167060, HD 106252, and HD 37962 using the exact same methods as for our
WASP-69b science data to improve the flux-calibration of our observations. We
describe the use of these standard star observations below.

We compared the measured absolute fluxes of the star from before eclipse to
stellar models, as shown in Figure 4. While the measured fluxes from before
eclipse include the planet emission, the planet contributes less than 0.3% to
the total flux so we ignore it at the level of calibration and systematic
errors (several %). We use the stellar modeled flux at the photosphere of the
star (i.e. $\pi$ times the intensity) and multiply it by the square of the
ratio of the stellar radius (0.813 R⊙, Anderson et al., 2014) to the Gaia DR3
distance (50.2871 pc, Gaia Collaboration et al., 2023). We were able to
confirm that a BOSZ model accurately predicted the spectrum to a 5% level, but
required an additional multiplicative factor to best match the data, which is
1.0266. In other words, the combined flux multiplicative factor was
1.36399$\times 10^{-19}$ times the modeled flux at the stellar photospheric
“surface”, which is $\pi$ times the modeled intensity. The multiplicative
offset of 1.027 can either be due to calibration flux uncertainties or
physical system parameter differences such as the distance to the system and
radius. Given 2% uncertainties in the flux calibration of the NIRCam slitless
grism mode and MIRI LRS mode (STScI, 2016, accessed on 2024-02-01), this
multiplication factor is within the calibration errors of both modes. The
strength of the CO absorption feature near 4.5 µm is also slightly deeper by
1-3% in the model as compared to the measured flux, either due to model
systematics or differences in stellar abundances. We also used our extracted
and calibrated spectra for HD 167060, HD 106252 and HD 37962 to check for
systematic calibration errors. We used the average calibration factor (from Jy
to DN/s) for these three stars using CALSPEC models and used this as an
alternative calibration factor for WASP-69 b’s calibration. In this re-
normalized spectrum, the drop from 5-7 µm is still present, the residuals
between 10-12 µm are smaller and the overall spectrum is shifted upward by 5%.
Given that the stellar model does not show any large ($>3\%$) discrepancies
with the calibrated spectrum, we can rule out the possibility that debris
disks can significantly dilute the spectrum of the planet at the measured
precisions in $\sigma_{F_{\rm p}/F_{*}}$ (3-16%, where $F_{\rm p}$ is the flux
of the planet and $F_{*}$ is the flux of the star and $\sigma$ represents the
uncertainty on the ratio of the two fluxes at the binned resolution).

We also checked our absolutely calibrated spectrum against 2MASS (Skrutskie et
al., 2006) and WISE (Wright et al., 2010) photometry. The 2MASS and WISE
photometry show slightly elevated fluxes as compared with the NIRCam grism
spectra. It is possible that there are residual non-linearity effects or the
brighter-fatter effect (Plazas et al., 2018) that degrades the linearity
correction of individual pixels. A general analysis of the absolute
calibration well below JWST requirements of 10% accuracy for spectroscopy
(Gordon et al., 2022) is beyond the scope of this work. However, given that
there are no large deviations from the stellar model, no significant evidence
for dust excesses, a small discontinuity between NIRCam F444W and MIRI LRS and
future improvements to the JWST flux calibration ahead (Gordon et al., 2022),
we proceed in using the stellar model and multiply it by $F_{\rm p}/F_{*}$ to
derive the planet’s dayside flux.

##  4 General Properties of WASP-69 b’s Emission Spectrum

Figure 5 shows the combined Eureka! MIRI and tshirt NIRCam emission spectrum,
with NIRCam data binned to lower resolution for visualization purposes. The
bins for visualization are about 0.05 µm for NIRCam, as compared to the
original lightcurve bins of about 0.01 µm with the MIRI LRS binning of 0.25 µm
already sufficiently low enough for visualization purposes. Figure 5 also
shows some representative blackbody planet spectra. For these blackbody
spectra, we multiply the Planck function by the planet-to-star radius ratio
squared from Table 2 and divide this by the stellar intensity from our BOSZ
stellar model (Bohlin et al., 2017) with a stellar effective temperature of
4750 K, log(g)=4.5 log(cm s-2), described in Section 3.4. The measured
spectrum clearly is not well fit by a blackbody and shows absorption features
and crosses from high to low blackbody temperature from short to long
wavelengths.

We also show the Spitzer secondary eclipse values for the 3.6 µm and 4.5 µm
bandpasses in Figure 5 (421 $\pm$ 29, 463 $\pm$ 39; Wallack et al., 2019) and
that they are shallower than the NIRCam spectroscopy. We calculate the
synthetic photometric eclipse depths using our JWST emission spectra and a
photometric response curve interpolated to a common grid to compare the
results between JWST and Spitzer. We find synthesized eclipse depths of 542
$\pm$ 7 ppm and 644 $\pm$ 7 ppm for the 3.6 µm and 4.5 µm bandpasses
respectively, assuming that each pixel is statistically independent. In both
bands, the JWST NIRCam eclipses are significantly deeper than the Spitzer
values (Wallack et al., 2019) of 421 $\pm$ 29 ppm and 463 $\pm$ 39 ppm for the
3.6 µm and 4.5 µm bands, respectively. For context, a previous comparison with
the Spitzer IRAC 3.6 µm eclipse depth and the synthetic NIRCam IRAC 3.6 µm
eclipse depth for WASP-39 b indicated consistency to within the 1 $\sigma$
uncertainty (Ahrer et al., 2023); however the Spitzer uncertainty on this
transit depth is 176 ppm in Sing et al. (2016), so the 120 and 180 ppm offsets
we find in WASP-69 b would not be significant for WASP-39 b. Bean et al.
(2023) find a 4.4$\sigma$ discrepancy with the 3.6 µm eclipse depth from Zhang
et al. (2018), ie 84 ppm $\pm$ 19 ppm for the hot Jupiter HD 149026 b. Some of
the difference we find between NIRCam and Spitzer IRAC for WASP-69 b could be
due to red noise in the Spitzer lightcurves, which exhibited deviations from
$1/\sqrt{N}$ statistics (where the noise for $N$ averaged data points is
$1/\sqrt{N}$ times the non-averaged data). This was observed for the second
visits in both the 3.6 µm and 4.5 µm filters (Wallack et al., 2019). Another
possibility is stellar or planet variability and/or a different set of
systematics with JWST NIRCam versus Spitzer IRAC.

Given the large dynamic range of the eclipse depth (Fp/F∗) across our full
spectrum from 2 µm to 12 µm, it hard to see how individual gases affect the
spectrum and see where atmospheric models perform well. Instead, it is
preferable for visualization purposes to compare the brightness temperature of
the planet with the brightness temperature of atmospheric models because its
dynamic range is much smaller. We calculate the brightness temperature
spectrum for WASP-69 b in the following manner: First, we average the
intensity of the star from the BOSZ stellar model described in Section 4
across each wavelength bin where the intensity is evaluated in wavelength
space. Second, we multiply this stellar intensity by the planet-to-star flux
ratio and by the square of the star-to-planet radius ratio to evaluate the
planet intensity. Third, we invert the Planck function in wavelength space to
calculate the associated blackbody temperature that corresponds to the planet
intensity. Fourth, we propagate intensity error intervals to brightness
temperature error intervals by linearizing the inverse Planck function (with a
first order Taylor expansion) and evaluating the width of the temperature
uncertainty at its median intensity. The brightness temperature (as well as
the binning) are used for illustrative purposes to better view the measured
spectrum and model over a wide dynamic range, but the original lightcurve
wavelength bins (0.01 µm bins, rounded to the nearest pixel for NIRCam
spectroscopy, a filter bandpass from F210M and 0.25 µm for MIRI LRS) were used
to fit the models and perform retrievals using the eclipse depth (Fp/F∗).

The brightness temperature spectrum of WASP-69 b is shown in Figure 5 (second
panel from the top and bottom panel) for the binned version of the spectrum
for visual clarity. Additionally, we show the relative opacities for abundant
and observable gases at the temperature of WASP-69 b in the second from the
bottom panel. The brightness temperature spectrum more clearly shows the
spectral features of H2O near 2.8 µm and 6.7 to 8.0 µm, and CO2 near 4.3 µm.
The red edge of the CO2 feature also shows excess absorption due to CO near
4.7 µm. CH4 is expected to be highly abundant in chemical equilibrium at the
zero albedo full redistribution equilibrium temperature of WASP-69 b (963 K)
for near-solar composition at 0.1 bar (Moses et al., 2013). However, our data
do not show any strong CH4 features near 3.3 µm nor 7.7 µm.

Figure 5: The emission spectrum of WASP-69 b shows significant absorption
features from molecules in its atmosphere and a trend from high to low
brightness temperature from 2.1 µm to 12µm. Top: Binned emission spectrum of
WASP-69 b in terms of planet-to-star flux for the NIRCam and MIRI data,
colored by their respective filters. The full higher resolution NIRCam data
(which are used for atmospheric retrieval) are shown as points and the binned
spectra are shown for illustrative purposes as lines with error bars. The MIRI
full resolution is the same as the binned one due to its intrinsically lower
resolution. Synthetic photometry is calculated for the IRAC 3.6 and 4.5
channels (gray circles) indicates larger eclipse depths than Spitzer published
measurements (purple points with error bars Wallack et al., 2019). Blackbody
planetary spectra for example temperatures are shown as black lines. Middle:
The spectrum is poorly fit by a Clear 1D Homogeneous model (dashed blue line)
where it over-predicts the flux at MIRI wavelengths (5-12 µm), over-predicts
the 3.9 µm flux and under-predicts the shortest wavelengths (2.5-3.2 µm). It
is necessary to include the effects of clouds (thin orange line) as well as
dayside inhomogeneities (thick green line) to fit the spectrum of the planet.
Bottom: The brightness temperature spectrum shows absorption signatures of
H2O, CO2 and CO. The relative pressure level where each molecule has an
optical depth of 1 is shown for each model in the style of Mukherjee et al.
(2024). We plot 2 models (described in Section 5, which can alternatively fit
the spectrum of WASP-69 b with either significant reflection (Scattering
Model) or a high altitude cloud (Cloud Layer Model).

##  5 Atmospheric Retrievals

###  5.1 Initial Fits with 1D Atmospheric Models for a Homogeneous Dayside

We interpret the measured emission spectrum with theoretical atmospheric
models beginning with the simplest assumptions. We first attempted to fit the
spectrum under the assumption that the disk-averaged emission from a planetary
dayside can be represented by a single 1D atmospheric column, ie. a
homogeneous dayside:

  1. 1.

We initially do not include any aerosols in this model. We term this model the
“Clear One-Region Model.”

As will be described below, the “Clear One-Region Model” provides a poor fit
to the data, so that is why we also include a cloud for the dayside following
results of the transmission spectra that show evidence for aerosols (e.g.
Khalafinejad et al., 2021; Estrela et al., 2021).

  1. 2.

In a second fit we term the “Cloudy One-Region Model”, we assume a vertically
uniform gray cloud opacity.

We computed a grid of self-consistent radiative-convective equilibrium
temperature-pressure (TP) profiles using the Extrasolar Giant Planet (EGP)
code (Marley & McKay, 1999; Fortney et al., 2005; Marley & Robinson, 2015;
Thorngren et al., 2019) for possible combinations of atmospheric metallicity
of [M/H]=0.0, 0.5, 1.0, 1.5 and 2.0, C/O ratio of 0.25$\times$, 0.5$\times$,
1.0$\times$ and 2.0$\times$ solar C / O value, and the heat redistribution
factor of $f_{\rm redist}=0.5,0.6,0.7,0.8,0.9,1.0$. The redistribution factor
is $f_{\rm redist}=0.5$ in the case of full redistribution of absorbed stellar
radiation over a solid angle of $4\pi$, whereas $f_{\rm redist}=1.0$ is
dayside-only redistribution of absorbed stellar radiation to a solid angle of
$2\pi$ (ie. a hemisphere).333The Dayside Effective Temperature is
$T_{day}=T_{*}\sqrt{R_{*}/a}(1-A_{\rm B})^{1/4}(0.5f_{\rm redis})^{1/4}$,
where T∗ is the stellar effective temperature, $R_{*}$ is the stellar radius,
$a$ is the semi-major axis and $A_{\rm B}$ is the Bond Albedo. We then use the
open-source radiative transfer code CHIMERA (Line et al., 2013a) to compute
the emission spectrum from the atmospheric RCE grid, where the molecular
abundances at each pressure level are computed by a precomputed equilibrium
chemistry table implemented by CHIMERA in default. We used PyMultinest
(Buchner et al., 2014), a python implementation of the Multinest tool (Feroz
et al., 2009), to obtain posterior distributions for [M/H], C/O, and $f_{\rm
redist}$. The temperature at each pressure level for sub-grid point values of
${\rm[M/H]}$, C/O, and $f_{\rm redist}$ are interpolated from the RCE grid
with the RegularGridInterpolator in Python scipy library, as conducted in Bell
et al. (2023c). We assumed a uniform prior for each parameter and set the
Nested Sampling live points to $500$.

As shown in the second from the top panel of Figure 5, the Clear One-Region
model tested here is unable to explain the emission spectrum of WASP-69b. The
best-fit spectrum has a $\tilde{\chi}^{2}$ of 14.5. The difficulty originates
from the difference in brightness temperature between the NIRCam and MIRI
bandpass. The data reveal the average brightness temperature of $\sim
1050\leavevmode\nobreak\ {\rm K}$ at $2.1$–$4\leavevmode\nobreak\
{\rm{\mu}m}$, which is much hotter than the temperature of $\sim
950\leavevmode\nobreak\ {\rm K}$ at $>5\leavevmode\nobreak\ {\rm{\mu}m}$. The
Clear One-Region model underestimates the brightness temperature at
$2.1$–$3.3\leavevmode\nobreak\ {\rm{\mu}m}$ and overestimates it at
$>5\leavevmode\nobreak\ {\rm{\mu}m}$. This result causes a dilemma for the
one-region model in fitting the data: the atmosphere needs to be hotter to fit
the planet’s emission at $2.1$–$4$ µm but it further worsens the discrepancy
at $>5\leavevmode\nobreak\ $µm. The addition of a gray cloud opacity in the
“Cloudy One-Region Model” shown in Figure 5 helps lower the flux at long
wavelengths and at 3.9 µm but still cannot achieve a high enough brightness
temperature from 2.1 to 3.3 µm and simultaneously low enough brightness
temperature from 6 to 9 µm. In the following sections, we overcome the
difficulty found here by investigating more detailed atmospheric properties.

###  5.2 Models With Additional Parameters and Complexity

Given that our initial fit does not match the observed dayside spectrum of
WASP-69 b well, we added additional complexity to the models. We considered a
variety assumptions about aerosols, reflected light, the temperature-pressure
profiles and dayside inhomogeneities. We used two different radiative transfer
codes to simulate the emission spectra of WASP-69 b that each take as inputs
temperature-pressure profiles and chemical abundances: PICASO (Batalha et al.,
2019; Mukherjee et al., 2023) and CHIMERA (Line et al., 2013b). The
independent sets of radiative transfer models provide robustness to the
results to a particular radiative transfer implementation. We include the
following types of reflection assumptions, aerosols and 3D effects, which are
described in more detail later in this section:

  1. 3.

A “Cloudy Two-region model” that includes two different temperature-pressure
profiles and sets of cloud properties between two geographically distinct
regions. Scattering is not significant. (7 free parameters)

  2. 4.

A “Scattering” model that includes a wavelength-independent geometric albedo
parameter as a free “knob” to the model, which can approximate the effects of
aerosols without a specific assumption about their composition or sizes. This
adds relatively more flux at short wavelengths than long wavelengths because
the thermal contribution decreases at short wavelengths but the scattering
term remains constant. (10 free parameters)

  3. 5.

A “Cloud Layer” model that includes a distribution of silicate condensates
that can emit at a different brightness temperature than the underlying
molecular emission. There is no free parameter for the albedo like Model 4,
but there is a modest level of scattering due to molecules and silicate
condensates. (14 free parameters)

For model 3 we use the Extrasolar Giant Planet (EGP) code (Marley & McKay,
1999; Fortney et al., 2005; Marley & Robinson, 2015; Thorngren et al., 2019)
with radiative transfer calculated by the CHIMERA code (e.g. Line et al.,
2013b). For models 4 and 5 we use the PICASO code (e.g. Batalha et al., 2019;
Mukherjee et al., 2023). We introduce more detailed model descriptions in what
follows.

####  5.2.1 Cloudy Two-Region Model

It has been known that modeling exoplanet daysides with a single TP profile,
as adopted in the One-Region models, has an inability to fit the emission
spectrum if there is a strong temperature contrast in the dayside (Feng et
al., 2016; Taylor et al., 2020). In fact, the mock retrieval of Feng et al.
(2016) demonstrated that a single TP profile cannot produce bright emission at
$\sim 2$–$3\leavevmode\nobreak\ {\rm{\mu}m}$ produced by a dayside with strong
temperature contrast (see their Figure 4), which is reminiscent of the
difficulty found in Section 5.1. Fitting a homogeneous dayside model to planet
with an inhomogeneous dayside can also result in biased abundances and
spurious molecular detections (Feng et al., 2016; Taylor et al., 2020).

To account for the dayside inhomogeneity, we introduce the Two-Region model
that computes the emission spectrum using 2 TP profiles. The setup of the
model is largely the same as the one-region model presented in Section 5.1. We
utilize the grid of radiative-convective equilibrium TP profiles used in the
one-region models (see Section 5.1), although we have extended the heat
redistribution factor to $f_{\rm redist}=2$ because the dayside now has a
local region that is hotter than the dayside average. We split the dayside
into hot and cool regions and compute the the emission spectrum from each
region using the CHIMERA. The observable spectrum is then computed as

$F_{\rm obs}=x_{\rm hot}F_{\rm p}(f_{\rm redist,hot},\kappa_{\rm
cld,hot})+(1-x_{\rm hot})F_{\rm p}(f_{\rm redist,cold},\kappa_{\rm
cld,cold}),$ (1)

where $x_{\rm hot}$ is the areal fraction of the hotter region of the dayside,
and we have assigned different values of the heat redistribution factors
$f_{\rm redist,hot}$ and $f_{\rm redist,cold}$ for each hot and cool regions.
$\kappa_{\rm{cld,hot}}$ and $\kappa_{\rm{cld,cold}}$ are the gray cloud
opacities of the hot region and cold region, respectively. Since the cloud
properties can also be very different at each region due to the distinct TP
profiles, we also independently assign the gray cloud opacity for hot and cold
regions. We note that the gray cloud is treated as an purely isotropic
scattering opacity source in CHIMERA. The atmospheric metallicity and C/O
ratio are expected to be horizontally uniform and thus are assumed to be the
same in the hot and cool regions. As in the One-Region models, we use
PyMultinest with 400 live points to obtain the posterior distribution of 5
climate parameters: [M/H], C/O, $f_{\rm redist,hot}$, $f_{\rm redist,cold}$,
and $x_{\rm hot}$ and 2 cloud parameters: $\kappa_{\rm redist,hot}$,
$\kappa_{\rm redist,cold}$ for a total of 7 parameters. We note that the model
includes a gray cloud for computing the spectrum, but ignores the radiative
effects of clouds on TP profiles. For the priors on the temperature of the
cold region, we restrict $f_{\rm redist,cold}$ to be between 0.5 and 1.0 (ie
between full and and no heat re-distribution). For the hot region, we set a
wider prior on $f_{\rm redist,hot}$, from 0.5 and 2.0. The upper bound
corresponds to the substellar point temperature in pure radiative equilibrium.
This is higher than the no redistribution case discussed above, because the
latter corresponds to the averaged dayside temperature, whereas the $f_{\rm
redist,hot}=2$ bound corresponds to the maximum possible local temperature.

####  5.2.2 Scattering Model

We use the open-source PICASO model (Batalha et al., 2019; Mukherjee et al.,
2021) to model the planet’s atmosphere within a 1D Bayesian retrieval
framework. The temperature-pressure profile is assumed to be in an analytic
form of Line et al. (2013b), with five free parameters: the equilibrium
temperature Teq (simply a parameter and not necessarily consistent with the
planet’s equilibrium temperature), log(g), Planck mean infrared opacity
$\kappa_{\rm IR}$, internal temperature $T_{\rm int}$ and the relative
fraction of the second visible stream in the two-stream approximation
$\alpha$. We assume thermochemical equilibrium throughout the atmosphere in
our retrieval framework. For a given $T(P)$ profile, atmospheric metallicity,
and C/O ratio, the FASTCHEM chemical equilibrium model (Stock et al., 2022) is
used to generate the abundance profiles of gases like CO2, CO, CH4, NH3, H2O,
etc. The solar composition elemental abundances from Lodders et al. (2009) are
scaled for different metallicities and C/O ratios within the retrieval. The
C/O ratio for a given metallicity is varied by scaling the C and O elemental
abundances such that the total C+O remains unaltered. In addition to
thermochemical equilibrium gases, we also retrieve on a constant (with
altitude) SO2 abundance to capture possibility of photochemically produced SO2
(Tsai et al., 2022). This SO2 abundance is not self-consistent with the
chemical equilibrium chemical profiles, but it gives a clue whether
photochemical calculations such as with VULCAN (Tsai et al., 2021) may be
necessary. These chemical abundances along with the $T(P)$ profile are then
used to generate the 1D emission spectrum of the planet with PICASO. To
account for the horizontal temperature contrast, the model also adds a
dilution parameter ($s_{\mathrm{dilute}}$), which is the fractional area of
the planet’s dayside that is emitted by this 1D model. This smaller emitting
area can compensate for the 2D structure of the planet and mitigate the biases
in retrievals (Taylor et al., 2020).

In the scattering model retrieval setup, we include a constant reflected light
term within our retrieval model in addition to the thermal component to fit
the eclipse spectrum. The planet-to-star flux ratio spectrum in this case is
defined as,

$\dfrac{F_{\rm planet}(\lambda)}{F_{\rm star}(\lambda)}=A_{\rm g}\dfrac{R_{\rm
p}^{2}}{a^{2}}+\dfrac{F_{\rm planet,thermal}(\lambda)}{F_{\rm star}(\lambda)}$
(2)

where $A_{\rm g}$ is the geometric albedo of the planet, Rp is the planet
radius, and $a$ is the semi-major axis. We choose the $A_{\rm g}$ term to be
wavelength-independent to simply asses the amount of starlight reflection that
would be needed to explain the near-infrared data, even though such a
wavelength-independent reflected light term is not very realistic and we do
not include the effect that scattering can have on radiative transfer. For the
prior on the geometric albedo, we put in a strict upper limit of 0.67 on
$A_{\rm g}$. We use this upper limit because if the $A_{\rm g}$ is greater
than 0.67, then the Bond Albedo of the planet would be higher than 1 if the
planet is assumed to be Lambertian. We also note that this retrieval setup can
violate energy conservation if the thermal component plus the reflected light
component are larger than the total incident energy on the planet. We sample
the posterior distribution with the Dynasty code (Speagle, 2020).

####  5.2.3 Cloud Layer Model

We use a second model setup with the same PICASO framework, equilibrium
chemistry, parameters for the metallicity, C/O ratio, and analytic Line et al.
(2013b) Temperature-Pressure profile but with a cloud deck. We also use the
Dynasty sampler (Speagle, 2020). Our cloud deck setup attempts to fit the data
without an arbitrary free Geometric Albedo that allows high values of
reflected light. Instead, we simulate the (small) contribution to scattering
light from molecular Raleigh scattering and condensate clouds and the (larger)
effects on radiative transfer of planet thermal emission that the clouds can
have. We include a parametric form of condensate cloud deck in this retrieval
setup, which can be used with any cloud species’ optical properties. The cloud
is parameterized assuming a log-normal cloud particle size distribution, where
$r_{\rm mean}$ is the geometric mean (and also the median) of the particle
size distribution and $\sigma$ is the geometric standard deviation the
particle size distribution as in Ackerman & Marley (2001) Equation 9. We
assume that $r_{\rm mean}$ and $\sigma$ remain the same for all cloudy
atmospheric layers and calculate the layer-by-layer cloud optical depth
$\tau_{\rm{cld}}$, asymmetry parameter $g_{0}$, and single scattering albedo
$\rm{w}_{0}$ using Mie scattering calculations. These layer-by-layer Mie
properties are then scaled using the normalization factor–$ndz$. This lets us
calculate the cloud Mie properties at the base of the cloud deck. We fit for a
base pressure – $P_{\rm base}$ deeper than which the cloud optical depth is
zero. At pressures smaller than $P_{\rm base}$, the layer-by-layer optical
depth of the cloud deck is calculated using $r_{\rm mean}$ and $ndz$ is scaled
using the scaling factor,

$f=e^{-{f_{\rm sed}}z/H}$ (3)

where $z$ is altitude relative to the cloud deck base. We set the reference
altitude $z$ arbitrarily and the scale height $H$ to be a constant throughout
the atmosphere and use $f_{\rm sed}$ as a free fitting parameter. We fit for
the five free cloud parameters $log(r_{\rm mean})$, $\sigma$, $P_{\rm base}$,
${f_{\rm sed}}$, and $ndz$ in this setup. We also include a dilution parameter
($s_{\mathrm{dilute}}$) to approximate the 3D effects of day-to-terminator
temperature contrasts, as for the Scattering model described in Section 5.2.2.

We apply our cloud layer retrieval setup using optical properties of Enstatite
(MgSiO3) (Scott & Duley, 1996). We initially tried other sulfide compositions
such as Na2S but found unrealistically high temperatures that deviated
significantly from radiative-convective equilibrium grids.

###  5.3 Retrieval Results

####  5.3.1 Two-Region Retrieval

Figure 6: Retrieved Posterior Distribution for the Two-Region Model. The
physical constraints of the model give tight bounds on the metallicity and C/O
ratio of the planet. The hotter region fills about 68% of the dayside in this
retrieval.  Figure 7: Left: The Two-Region Model can fit the spectrum with a
hot temperature-pressure profile that covers 68% of the dayside area and a
cooler cloudier region that contributes negligibly to the planet flux. The
binned data are shown as gray diamonds with error bars. The individual spectra
from the hot component (gray dotted line) and cold component (gray dashed
line) are combined via Equation 1 to give a best fit spectrum (thick cyan
line). A model with no reflected light and only thermal flux (thin red line)
and the hot component’s emission only (dotted orange line) show that reflected
light and the cold component’s emission are both negligible compared to the
thermal emission from the planet’s hot region. This also means that the
dilution parameter used in the Cloudy Layer and Scattering models should give
a very similar results as two separate regions. Right: Our measured spectrum
shows signatures of H2O, CO2 and CO molecules as well as clouds to both
suppress the emission from the cold region and the 4 µm and 9–11 µm emission
from the hot region. We remove one molecule at at time and plot the change in
spectrum as a shaded region spanning from the best fit model (cyan line) to
the molecule-removed and cloud-removed spectra (color-shaded regions). The
binned data data (gray diamonds) are the same as on the left plot. We note
that removing H2O can decrease the flux as compared to the best-fit model
because it increases the single scattering albedo and thus decreases the
emissivity of the photosphere.

The introduction of a dayside inhomogeneity greatly improves the model ability
to explain the data. Figure 5 shows the median emission spectrum of the two-
region model. While One-Region models over-predicts the long wavelengths and
under-predict the short wavelengths, the Two-Region model better fits the
overall shape of the spectrum (see Figure 5 and 7). This result demonstrates
that it is critical to accounting for the temperature contrast on the day
side, as suggested by Feng et al. (2016) and Taylor et al. (2020). Figure 6
shows the posterior distributions of the two-region model. We retrieved a
super-solar metallicity of [M/H] = 1.01${}^{+0.08}_{-0.06}$ and moderately
super-solar C/O ratio of $0.79^{+0.03}_{-0.03}$. It should be noted that the
uncertainties of the parameters retrieved here are likely underestimated, as
our two-region model relies on the TP profiles interpolated from the 1D RCE
grid, which has less flexibility to change the TP profile compared to
parameterized TP profiles. It is also interesting to mention that the cloud
opacity retrieved for the cold region is much higher than that for the hot
region, which indicates a strong cloudiness contrast on the dayside.

To further elaborate on the results of the two-region model, the left panel of
Figure 7 shows the median emission spectrum along with individual emission
spectra from the hot and cold regions. The data can be well fitted by the hot
region with an area fraction of $68\%$. A comparison between the median
spectrum with the diluted hot region spectrum further reveals that the
emission from the cold region has negligible impacts on the total thermal
emission from the dayside, which is attributed to thick clouds in cold
regions. Thus, our result indicates that the overall shape of the observed
spectrum is controlled by the emission from the hot regions where it is hotter
and less cloudy as compared to the cold regions. We also note that the median
spectrum is invariant if we omit the reflected light component, indicating
that the reflected light is negligible in the Two-Region model. Given that the
cold region’s flux is negligible due to thermal scattering in the atmosphere,
the temperature on the cold region largely unconstrained by the data. The
$f_{\rm redist,cold}$ parameter pushes up against our prior of 0.5 (full
redistribution of heat) as shown in Figure 6. Thus, we only find an upper
limit on the redistribution factor and the data do not probe the deeper layers
as will be discussed in the contribution functions described in Section 7.3.
That being said, the temperatures implied by $f_{\rm redist}$ (1220$\pm$12 K
in the hot region and $\sim$1030 K in the cold region) are roughly consistent
but slightly below the temperature contrast in a General Circulation Model for
WASP-69 b in a cloudless simulation (Mehta & Parmentier, 2024).

The observed spectrum shows the absorption features of several molecules. As
shown in the left panel of 7, the spectrum clearly shows the absorption
feature of CO2 at ${\sim}4.3\leavevmode\nobreak\ {\rm{\mu}m}$ that is a strong
indicator of a high metallicity atmosphere (e.g., Moses et al., 2011; JWST
Transiting Exoplanet Community Early Release Science Team et al., 2023). The
spectrum also indicates the presence of CO that is needed to explain the
relatively low flux at ${\sim}4.6\leavevmode\nobreak\ {\rm{\mu}m}$. The H2O
absorption and cloud opacity have comparable impacts and control the overall
shape of the spectrum. On the other hand, the spectrum does not show any
noticeable CH4 features. As will be discussed in Section 7.2 , the atmospheric
TP profile inferred from the observed spectrum leads to an equilibrium CH4
abundance of $\lesssim{10}^{-6}$ at $P<{10}^{-1}\leavevmode\nobreak\ {\rm
bar}$, which is orders of magnitude lower than the abundances of CO and H2O.
This makes the CH4 feature unnoticeable even without the disequilibrium
quenching from the deep hot atmosphere. However, this would not be obvious
under the assumption of a homogeneous model that is closer to the zero-albedo
full-redistribution equilibrium temperature of 963 K for WASP-69 b, which
would contain more significant CH4 absorption.

The observed spectrum indicates the presence of clouds in both the hot and the
cold regions. The hot region needs to be veiled by a moderate amount of
clouds; otherwise, the spectrum around $4\leavevmode\nobreak\ {\rm{\mu}m}$
becomes too bright to explain the observation. The cold regions should be
veiled by clouds much thicker than those in the hot regions to suppress the
emission from the cold region. The model fails to explain the faint emission
at $5\leavevmode\nobreak\ {\rm{\mu}m}$ if the cold regions have clear
atmospheres. We will discuss the potential cloud compositions in the hot and
cold regions in Section 7.3.

Our two-region model still struggles to explain the sudden increase of the
observed emission at ${\gtrsim}9\leavevmode\nobreak\ {\rm{\mu}m}$. Since
removing clouds from the hot region leads to thermal emission comparable to
the observed value at ${\gtrsim}9\leavevmode\nobreak\ {\rm{\mu}m}$, the data
potentially indicate that clouds in the hot region begin to be transparent at
${\sim}9\leavevmode\nobreak\ {\rm{\mu}m}$ through Mie scattering, which cannot
be modeled by our gray scattering clouds in the two-region model. If this is
true, the clouds in the hot region would mainly consist of ${\sim}(9/2\pi)$=
1.4 µm cloud particles. Another intriguing possibility is the sudden decrease
in the cloud’s single scattering albedo at ${>}9\leavevmode\nobreak\
{\rm{\mu}m}$, which we will further discuss in Section 7.3.

Figure 8: Posterior densities of retrieved atmospheric parameters from PICASO
for the scattering model (red) and the cloud layer model (blue). The modeled
parameters (from left to right and down are the atmospheric log metallicity
relative to solar composition ([M/H]), absolute carbon-to-oxygen ratio (C/O),
geometric albedo ($A_{g}$), scattering model only), dayside area dilution
($s_{\mathrm{dilute}}$), Five Temperature-Pressure Profile parameters from
from Line et al. (2013b): the $T_{eq}$ parameter (not consistent with the
planet’s equilibrium temperature, described in the text), $log(g)$,
$\kappa_{IR}$, the internal heat flux temperature $T_{int}$,$\alpha$, and
finally a free-floating SOs abundance (SO2). Stair-step 2D representations of
these posteriors can be found in Figures 16 and 17.

####  5.3.2 Scattering Model Retrieval

We show the posterior distribution of the Scattering Model spectrum in the
bottom panel of Figure 5. The posterior distributions on different parameters
like atmospheric metallicity, C/O ratio, dilution factor, and geometric albedo
obtained from this retrieval setup is shown in the panels of Figure 8 in blue.
This retrieval setup constrains the atmospheric metallicity of the planet to
be supersolar at with [M/H] = +0.75$\pm$0.15 dex and the absolute C/O ratio at
0.42$\pm$0.16. The dilution parameter is significantly below 1.0, indicating
that the emission is dominated by a central region of higher brightness
approximated by this 1D model. There is also a peak in the posterior abundance
of SO2 near a mixing ratio of $10^{-6}$, which likely improves the fit near
7.5 µm, but it has a long tail extending to lower mixing ratios. Given that
there is no significant absorption feature in the brightness temperature
spectrum at 7.7 µm we do not consider this significant evidence for SO2.

In our scattering model setup, we find that the retrieved geometric albedo
pushes right up against our prior of 0.67 in order to fit the spectra at short
wavelengths, as shown in Figure 8. The retrieved very high reflected component
is much higher than the albedo inferred in other hot Jupiter atmospheres which
are, for example, 0.096 $\pm$ 0.02 for HD 209458 b (Brandeker et al., 2022)
and 0.076 $\pm$ 0.02 for HD 189733 b (Krenn et al., 2023), as measured by
CHEOPS secondary eclipse, with a range of 0.02 to 0.27 for 1500 K to 1700 K a
sample of hot Jupiters (Adams et al., 2022). Recently, an ultra-hot (1980 K)
Neptune LTT 9779 b was found to have a large geometric albedo of
0.80${}^{+0.10}_{-0.17}$ from its CHEOPS optical secondary eclipse depth
(Hoyer et al., 2023), and even at 1980 K, the thermal contribution contributed
negligibly ($<$10) ppm for a 115 ppm eclipse. However, a Neptune-mass planet
like LTT 9779 b may have a different composition than a warm Jupiter like
WASP-69 b so it is unknown if Jupiter-mass planets can attain such high
geometric albedos. The $T(P)$ profile constrained from this retrieval is shown
along with the abundances of key molecular absorbers of CH4, CO, CO2 and, H2O
using blue in Figure 9.

####  5.3.3 Cloud Layer Retrieval

The fit to the brightness temperature spectra obtained with this setup is
shown in red in Figure 5. The posterior distributions of the parameters are
shown in Figure 8. This setup estimates the planet’s atmospheric metallicity
to be +0.96$\pm$0.2 above Solar and the absolute C/O ratio to be 0.75$\pm$0.1,
as seen in Figure 8. Thus, the compositional constraints depend on the
assumptions about aerosols and the level of reflection from WASP-69 b. As with
the Scattering retrieval described in Section 5.3.2, there is a peak in the
SO2 posterior, but it is at very low abundances where the 7.7 µm feature does
not show up significantly in the spectrum. Our constraints on the layer-by-
layer cloud optical depth, asymmetry parameter, and single scattering albedo
at a wavelength of 10 µm along with the constrained $T(P)$ profile is shown in
Figure 9. The retrieval prefers a silicate cloud deck near $\sim$ $10^{-4.5}$
bar pressure level extending up to $10^{-6}$ bar with a peak cloud optical
depth close to 0.1. The retrieved particle size is very small with a median
radius of 10-6 cm.

Figure 9 shows the retrieved temperature-pressure profiles for the Cloud Layer
model. The temperature-pressure profile crosses the saturation pressure curves
for Na2S and MgSiO3 (enstatite). This indicates that Na2S can form a cloud at
a similar altitude as the inferred cloud deck base $10^{-4.5}$ bar. However,
Na2S does not have the right optical properties to work in our Cloud Layer
model, which requires a high absorption coefficient at MIRI wavelengths ($5$
to $12$µm) as compared to short wavelengths. MgSiO3, on the other hand, has
the right optical properties (e.g. Taylor et al., 2021) to explain the low
($\sim$930 K) brightness temperatures at the MIRI wavelengths without
significantly absorbing the NIRCam wavelengths $2$ to $5$ µm. However, our
retrieved TP profiles for even the warmest Cloud Layer model does not cross
the MgSiO3 condensation curve shown in Figure 9, except at very low altitudes
(10 bars or deeper), so MgSiO3 particles would likely rain down from
observable pressures. The enstatite particles would need to be lofted to the
retrieved 10-4.5 to $10^{-6}$ bar pressure level from below 10 bars. This
would require an extreme level of vertical mixing (such as inferred for
WASP-107 b to loft high altitude silicate cloud particles (Dyrek et al.,
2024)).

The normalized contribution per wavelength is plotted in Figure 10, which
shows that the cloud dominates the emission at longer wavelengths (thus
explaining the lower brightness temperature) while gas near 10-1.5 bar
dominates the emission at short wavelengths (thus explaining the higher
brightness temperature). MgSiO3 clouds produce this wavelength dependence
because the absorption coefficient of these silicates is much larger at
wavelengths longer than 4 µm than at short wavelengths (e.g. Taylor et al.,
2021). The cloud emission at long wavelengths has a lower brightness
temperature because it is emitted by a cooler upper layer. The gas emission at
short wavelengths has a higher brightness temperatures because it comes from a
warmer layer below the cloud deck. We also include reflected light from
Rayleigh scattering by molecular gas and silicate clouds, but their Albedo is
negligible ($\sim$10 ppm or A${}_{g}\approx 0.08$) compared to the thermal
component or the high geometric Albedo inferred from the Scattering Model
(A${}_{g}\approx$0.64) described in Section 5.2.2.

Figure 9: Temperature-Pressure Profile (upper left) and cloud properties for
the silicate cloud model . The TP profile for the Cloud Layer model crosses
the condensation curve for Na2S (dashed gray line) at very high altitudes and
MgSiO3 (dashed pink line) at very deep altitudes, while the Scattering model
crosses neither of these cloud candidate’s condensation curves. While the
Cloud Layer crosses these two condensation curves, Na2S has inefficient
nucleation rates from theoretical models (Gao & Benneke, 2018) and MgSiO3
would require extreme verticle mixing to be lofted to the high altitudes
needed to match the observed brightness spectrum. The upper right 3 panels
show the optical properties of the retrieved silicate clouds: optical depth
$\tau_{cld}$, asymmetry $g_{0}$ and single scattering albedo $w_{0}$ at a
reference wavelength of 10 µm. Figure 10: Contribution functions for the Two
Region model (top row) and the Cloud Layer model (bottom plot). The individual
regions’ contributions are shown for the Two Region model (top left two plots)
as well as the combined full dayside (right plot). A high altitude ($10^{-4}$
to $10^{-5}$ bar) cloud dominate the longer wavelength emission (with smaller
brightness temperatures) while the warmer lower layers dominate at short
wavelengths (with larger brightness temperatures) in both models. The
posterior retrieved Temperature-Pressure Profile for the CHIMERA model (bottom
middle) is shown as a red curve with the temperature at the top axis.

####  5.3.4 Summary of Model Retrieval Results

Table 3: Summary of retrieved atmospheric metallicity and C/O ratio from each model setup. The most plausible scenarios are the Two-Region (Cloudy) and Cloud Layer retrievals Model | Retrieval Code | Dayside Geometry | TP Profile | Reflected Light | Metallicity [M/H] | C/O ratio | $\tilde{\chi}^{2}$ | log$Z$
---|---|---|---|---|---|---|---|---
One-Region (Clear) | CHIMERA | Homogeneous | RCE Grid | Molecular | $1.30^{+0.02}_{-0.02}$ | $0.11^{+0.02}_{-0.01}$ | 14.5 | $-588$
One-Region (Cloudy) | CHIMERA | Homogeneous | RCE Grid | Molecular | $0.96^{+0.03}_{-0.04}$ | $0.65^{+0.03}_{-0.03}$ | 3.07 | $-439$
Two-Region (Cloudy) | CHIMERA | 2 TP | RCE Grid | Molecular | $1.01^{+0.08}_{-0.06}$ | $0.79^{+0.04}_{-0.04}$ | 1.79 | $-255$
Scattering | PICASO | Diluted | Line et al. (2013b) | Constant AG | 0.75${}^{+0.15}_{-0.15}$ | 0.42${}^{+0.16}_{-0.16}$ | 2.00 | $-253$
Cloud Layer | PICASO | Diluted | Line et al. (2013b) | Silicate | 0.96${}^{+0.20}_{-0.17}$ | 0.75${}^{+0.19}_{-0.10}$ | 1.65 | $-256$

The 5 models considered in this work are summarized in Table 3, with the major
assumptions, derived abundances, reduced chi-squared ($\tilde{\chi}^{2}$) and
Bayesian Evidence (log(Z)). The 1D homogeneous models that were originally
explored in Section 5.1 under-predict the short wavelength planet flux and
over-predicts the long wavelength planet flux thus giving large
$\tilde{\chi}^{2}$ values. The difference in log Bayesian evidence between the
1D Cloudy Homogeneous model and the Two-Region model is 184, corresponding to
a 19 $\sigma$ difference. The three models that best fit the data (with a
$\tilde{\chi}^{2}$ of 2.0 or less) are the the Scattering Model, the Two-
Region Model and the Cloud Layer Model with $\tilde{\chi}^{2}$ of 2.0, 1.79
and 1.65 respectively. The Bayesian evidence favors the Scattering model, but
the extremely high Geometric Albedo of 0.64 is likely unrealistic. The Two-
Region model and Cloud Layer models, by contrast have negligible stellar
reflected light, as shown in Figure 7. While we regard the Scattering Model as
less plausible, it is still consistent with existing JWST data and can be used
to assess the robustness of other conclusions like the composition of WASP-69
b’s atmosphere inferred from the Cloud Layer and Two-Region Models. We further
discuss the synthesized inferences about WASP-69 b’s cloud properties and
atmospheric composition across all models in Section 7.

##  6 Eclipse Mapping

Figure 11: MIRI LRS broadband lightcurves and Residuals for Lightcurve Fits to
the Uniform Map as well as the model difference between non-uniform map and
the uniform map. The non-uniform map allows for a gradient from the sub-
stellar point to the terminator through the Y2,0 and Y1,0 spherical harmonics.
The Bayesian Information Criterion favors the non-uniform map (BIC=184.1
versus BIC=191.1), but higher signal to noise is needed to robustly map
WASP-69 b to high confidence . Figure 12: Posterior maps for the non-uniform
map fits to WASP-69 b’s dayside. The maps only include the Y2,0 and Y1,0
spherical harmonics to sense temperature gradients, but do not reveal location
information. The mean map (Left) and Error (Right) favor a dayside brightness
distribution that is centrally concentrated and inefficient at redistributing
heat. The plus sign shows the peak brightness, which is forced to be at the
substellar point by the limited number of spherical harmonics. .

Motivated by the inhomogeneous dayside inferred by both the Two-Region
retrieval and the Scattering Model’s dilution parameter, we investigated
whether the lightcurves show further and independent evidence of an
inhomogeneous dayside. Eclipse mapping (e.g. Williams et al., 2006; Rauscher
et al., 2007; de Wit et al., 2012; Majeau et al., 2012; Coulombe et al.,
2023b), can use the star as a spatial scanner to map some non-uniform features
of the planet. We used the broadband MIRI lightcurve, which had the highest
signal-to-noise eclipse signal to map the dayside of the planet.

First, it is necessary to use the highest available precision orbital
parameters to decrease the correlations of these parameters with maps. This is
achieved by fitting the publicly available TESS lightcurve from sector 55
taken between 2022-08-05 and 2022-09-01. The posterior parameters we find are
listed in Table 2.

We model the MIRI LRS lightcuve a with starry (Luger et al., 2019) and include
a systematic trend for the lightcurves. The lightcurve is modeled as

$F(t)=F_{a}(t)\left(1+Ce^{-(t-t_{\rm{min}})/\tau}\right)(A+Bx),$ (4)

where $F_{a}(t)$ is the astrophysical variation modeled by starry, $t$ is the
time, $t_{\rm{min}}$ is the minimum (ie. start) time, $\tau$ is an exponential
time constant, $C$ is the exponential amplitude constant and $A$ and $B$ are
polynomial baseline terms. $x$ is the scaled time

$x=2(t-t_{mid})/(t_{max}-t_{min}),$ (5)

where tmid is the mid-time.

We fit the broadband MIRI LRS lightcurve with a uniform map, a spherical
harmonic degree 1 map and a spherical harmonic degree 2 map and a fourth fit
with the spherical harmonics corresponding to dayside brightness gradients
only. We sample the posterior with No-U-Turns sampling with pymc3 (Salvatier
et al., 2016) and the pymc3-ext tools from Foreman-Mackey et al. (2021). All
of our model fits results show that there is excess noise as compared to the
photon and read noise estimate. The best-fit (maximum a-priori) spherical
harmonic degree 2 lightcurve model has a $\tilde{\chi}^{2}$ = 1.05, which is
smaller than the Uniform model, which has $\tilde{\chi}^{2}$ = 1.11, where
$\tilde{\chi}^{2}$ is the reduced-chi-squared metric. We included an error
inflation factor that was fit as a hyper-parameter for the uniform fit and
then that same uncertainty is used in cross-model comparisons. We find that
the theoretical photon and read noise error should be 33 ppm, whereas the
inflated error is 56 ppm. The ratio of the inflated to theoretical error is on
the higher end as compared with other MIRI LRS lightcurves, but less than the
commissioning target L169-9 b (Bouwman et al., 2023).= We also test whether
the extra degrees of freedom in the spherical harmonic degree 2 lightcurve are
justified with the Bayesian Information Criterion (BIC) and find that it is
214.5 for the spherical harmonic degree 2 versus 191.1 for a uniform map,
which favors uniform map due to its 8 fewer free parameters. We repeated an
eclipse mapping fit but with only two free parameters in the spherical
harmonic map besides the overall amplitude: Y2,0, Y1,0, which create a
gradient from the sub-stellar point to the terminator - see Schlawin et al.
(2023a) Figure 5 or Luger et al. (2019) Figure 1, for example. This time the
non-uniform map has $\tilde{\chi}^{2}$ = 0.99 and a BIC=184.0, so both
statistics favor a non-uniform map with dayside temperature gradients. The
spherical harmonic degree 1 map, on the other hand, is not statistically
preferred over the Uniform map.

To better visualize the residuals, we also show the difference between the two
best-fit models and the residuals for a uniform eclipse map fit in Figure 11,
bottom. The lightcurve for the non-Uniform explains some of the 50-100 ppm
deviations in the residuals at ingress and egress. While the evidence for a
non-Uniform map is still at a similar level as the excess noise in the data,
the MIRI LRS broadband lightcurve supports the possibility that a non-Uniform
dayside model may explain the spectrum.

Figure 12 shows the posterior map distributions for the spherical harmonic
degree 2 fits to the lightcurves. The maps suggest significant day-to-
terminator contrasts on the planet. We also analyzed the lightcurve
independently with the eigenmapping method (Rauscher et al., 2018) using the
same systematic model as described in Equation 4 using the ThERESA mapping
code in its single-wavelength implementation (Challener & Rauscher, 2022). The
BIC preferred 2 eigenmap components up to 2nd order spherical harmonics. The
ThERESA maps also significantly favored a non-Uniform dayside map over a
uniform one by $\Delta$BIC of 53.6. The excess noise in the lightcurves beyond
photon and read noise means that caution is warranted when interpreting the
maps. Transit observations (planned for JWST Cycle 2 GO program 3712 and Cycle
3 GO program 5924) and a full phase curve will better measure if this high
contrast is indeed needed for the planet because the transmission spectrum
will constrain the properties around the terminator (ie longitudes of $\pm$ 90
degrees). The black points in Figure 12 are the locations of peak brightness
for random sample posterior draws from the pymc3 No U-Turns sampling of the
lightcurve. The location of peak brightness is consistent with the substellar
point (5$\pm 14^{\circ}$,-5$\pm 24^{\circ}$). The large variance in the peak
brightness longitude is due in part to the baseline we fit to the data in
Equation 4. Non-flat baselines can decrease the precision of the m=1 spherical
harmonic terms, which constrain the peak brightness longitude (Schlawin et
al., 2023a). The ThERESA analysis of the same lightcurve indicates an eastward
hotspot shift of $\sim 20^{\circ}$ with less uncertainty in the ThERESA maps
due to the fewer number of mapping terms and free parameters, but an improved
orbital ephemeris is needed to better constrain the hotspot shift. This
ephemeris will be improved with the upcoming JWST transit observations of
WASP-69 b. As with confirming the central concentration of brightness from the
eclipse maps, a full phase curve would better constrain the longitude of
WASP-69 b’s hotspot.

##  7 Discussion

###  7.1 Energy Balance

Figure 13: The dayside temperature from our JWST spectrum indicates that heat
is not fully re-circulated around the planet. The allowed ranges of dayside
temperature from full redistribution of heat around the planet (solid line) to
zero redistribution of heat (dotted line) for a given Bond Albedo are bounded
within the gray region. The integrated thermal flux from the planet
extrapolated to 0 µm and $\infty\leavevmode\nobreak\ \micron$ is calculated
for our better-fitting models (thick red line, an orange line and and thin
blue line). The high geometric albedo of the Scattering model (blue) which
pushes against the prior of 0.67 strongly indicates a high Bond albedo so we
show values below 0.4 with a dashed line. The Scattering Model barely
conserves energy at the highest Geometric Albedo and would result in a very
extreme dayside-to-nightside temperature contrast..

The total planet emission as measured over the NIRCam and MIRI wavelength
ranges can be used to estimate the dayside temperature, which constrains the
Bond albedo and the heat recirculation efficiency of the planet. For the
energy budget, we assume that internal heat flux is negligible for this
circularized hot Jupiter planet, which is also supported by our internal
temperature $T_{\rm int}$ of less than 200 K shown in Figure 8. We calculate
the dayside effective temperature of the planet for the Scattering model,
Cloud Layer Model and Two-Region models, which were the best fits in terms of
$\tilde{\chi}^{2}$. We integrate the thermal component of a best-fit model
spectrum from 0 to $\infty$, which requires extrapolation beyond our measured
wavelengths from 2.0 µm to 12 µm. We then find the temperature of a spherical
homogeneous blackbody radiator the size of the planet that has the same flux.

This dayside effective temperature can be compared to the incoming radiation
from the host star to derive a Bond albedo. However, the day/night heat
recirculation efficiency cannot be measured with only an eclipse lightcurve of
WASP-69 b. A full phase curve would be needed to measure the global effective
temperature, as has been done for GJ 1214 b (Kempton et al., 2023). We
therefore show the combined constraints on both the Bond albedo ($A_{B}$) and
heat redistribution efficiency in Figure 13. We mark the the allowed dayside
temperatures in the gray region in Figure 13, which are bounded by full heat
redistribution ($f=0.5$ or $\varepsilon=1.0$ in the parameterization of Cowan
& Agol (2011)) and no heat redistribution ($f=\frac{4}{3}$ or
$\varepsilon=0$).

We examine the constraints on the redistribution efficiency, model-by-model
from their dayside temperature. The Scattering Model spans a wide range of
possible redistribituion efficiencies that depend on the Bond Albedo. However,
the high constant-wavelength Geometric Albedo inferred from the Scattering
model 0.64 or higher indicates that the Bond Albedo is also likely large. For
a Lambertian sphere that scatters equally at all wavelengths, the Geometric
Albedo is $\frac{2}{3}$ times the Bond Albedo (Heng et al., 2021). In all four
of the Solar System’s giant planets, the infrared geometric albedo is less
than the Bond Albedo (de Pater & Lissauer, 2001). If the Bond albedo of
WASP-69 b is indeed larger than 0.64 (ie the right hand side of the blue line
in Figure 13, the redistribution efficiency would have to be near $f=1.33$
(ie. no redistribution of heat). However, it is possible (with a preference
for back-scattering starlight), that the geometric Albedo at some wavelengths
exceeds the Bond albedo. In the Solar System giant planets, the visual
Geometric Albedo is on average $\sim\frac{2}{3}$ times the Bond Albedo.
Therefore, we mark all Bond Albedos below 0.4 with a dashed line for the
Scattering model. Given the likely high Bond Albedo for the planet, we
conclude that the Scattering Model’s dayside temperature implies an
inefficient heat redistribution for the planet (with $f>0.87$ or
$\varepsilon<0.56$ in the parameterization of Cowan & Agol (2011)). The Bond
Albedo can be much smaller (near 0.0) for the Two Region and Cloud Layer
models shown in Figure 13 (ie the left hand side of the Figure). However, the
dayside temperatures inferred from these models are significantly warmer, so
even for a Bond Albedo of 0.05, the heat redistribution must be inefficient
($f>0.7$ or equivalently $\varepsilon<0.76$ ) for the Two Region Model and
($f>0.8$ or equivalently $\varepsilon<0.64$) for the Cloud Layer Model. Thus,
in all of the models we fit to the dayside spectrum of WASP-69 b, full
redistribution of heat is ruled out and some significant day-to-night
temperatures are expected for the planet.

Dynamical models can help inform whether inefficient heat redistribution and
high Bond Albedo (thinner blue line in Figure 13) are plausible for WASP-69 b.
In models of zero albedo equilibrium temperatures of 1061 K with no clouds or
hazes, circulation is efficient with temperature contrasts of less than 100 K
between the day and night sides at 1 mbar pressures, but clouds can
significantly enhance the day-to-night temperature contrast at the same
pressure to 500 K as well as alter the photosphere altitude (Roman et al.,
2021). Therefore, depending on the clouds in an atmosphere, there can be a
large variety of day-to-night contrasts and recirculation efficiencies
possible for a planet. The phase curve amplitudes of similar temperature gas
giant planets in short orbits measured by the Spitzer space telescope have
also shown a wide variety of longitudinal brightness distributions are
possible at 1100-1500 K equilibrium temperatures, potentially influenced by
clouds, with some cooler planets like WASP-34 b showing significant day/night
contrasts (May et al., 2022). WASP-43 b also shows evidence for nightside
clouds and a day-to-night brightness temperature that varies between 1520 and
860 K for the two hemispheres (Bell et al., 2024). This is similar for our
emerging picture of WASP-69 b, that there are significant aerosols in the
atmosphere that can enhance the day-to-night contrast of the planet.
Furthermore, the inferred dayside inhomogeneity from the Two Region model
(Section 5.3.1), dilution factor in the Scattering model (Section 5.3.2) and
tentative results of eclipse mapping (Section 6), all point to a central
concentration of high temperature gas as compared to the limbs of the planet,
which would be expected for inefficient heat redistribution.

###  7.2 WASP-69 b’s Atmospheric Composition

All of the models explored in this work give posterior distributions that are
enriched in heavy elements as compared to solar composition. This largely
comes from the constraints on the dominant molecular absorbers in the observed
2 µm to 12 µm region: H2O, CO and CO2. Even the lowest metallicity retrieval
from the Scattering model has [M/H]= 0.75$\pm$ 0.15. The two models which give
the lowest $\tilde{\chi}^{2}$ and (the Cloud Layer and Two-Region models) give
consistent compositional constraints within 1$\sigma$, in addition to showing
similar cloud effects described in Section 7.3. Despite different assumptions
about the 3D geometry of the planet and the temperature-pressure profile(s),
the Cloud Layer and Two-Region retrievals overlap well in inferred chemical
enrichment over solar of [M/H]=0.96${}^{+0.20}_{-0.17}$ and
[M/H]=1.01${}^{+0.08}_{-0.06}$, respectively. Similarly, the posterior C/O
ratios overlap with values of 0.75${}^{+0.19}_{-0.10}$ and
0.79${}^{+0.04}_{-0.04}$ respectively. The Cloud Layer model has broader
compositional constraints because of the more flexible temperature-pressure
profile whereas the Two-Region model uses fixed profiles based on radiative-
convective equilibrium calculations. We did not include radiative feedback
from the clouds in the Two-Region model so the very small standard deviations
in the posterior metallicity and C/O ratio likely under-estimate model
uncertainty. Therefore, we recommend adopting a metallicity of 6$\times$ to
14$\times$ and a C/O ratio of 0.65 to 0.94 (from the Cloud Layer model) for
1$\sigma$ (68%) posterior intervals for modeling and comparative planetology.

We also note that the Scattering Model is not completely ruled out (and has
the largest log(Z) or Bayesian evidence in Table 3) so WASP-69 b’s composition
may be less metal-enriched than the Two-Region and Cloud Layer retrievals
indicate. While we consider this a less likely scenario, it is worth
considering the possibility that WASP-69 b’s atmosphere has a metallicity of
4$\times$ to 8$\times$ solar and a C/O ratio of 0.26 to 0.58 for 1$\sigma$
(68%) posterior intervals. The Scattering Model can be ruled out or confirmed
by measuring the optical reflection from WASP-69 b, which diverges strongly
from the Cloud Layer and Two-Region models below wavelengths of 2 µm, as shown
in Figure 5. CHEOPS (archived), WFC3 UVIS (planned) or JWST NIRISS SOSS
lightcurves of WASP-69 b can help answer this question at short wavelengths.
Until those data are analyzed and published, we carry forward these two
possible compositions for WASP-69 b.

WASP-69 b’s host star also shows enrichment relative to solar composition with
[Fe/H]=0.144$\pm$0.077 (Anderson et al., 2014) and another determination of
[Fe/H]= 0.29 $\pm$ 0.04 (Sousa et al., 2021). Machine-learning
characterization of the spectrum indicates [Fe/H]=0.41 $\pm$ 0.03, [C/H] =
0.27 $\pm$ 0.05 and [O/H]=0.18 $\pm$ 0.07 (Polanski et al., 2022). If indeed
WASP-69 A has a super-solar C/O ratio, this could explain the slightly super-
solar C/O ratio for WASP-69 b using the Cloud Layer and Two Region models.
More characterization of the star should be performed to better constrain its
metallicity and heavy element enrichment. If the star is super-solar in heavy
element enrichment of [Fe/H]=0.29 $\pm$ 0.04 from Sousa et al. (2021), WASP-69
b may only be $\sim$0.7 dex (ie. 5$\times$) above stellar composition. This
would place WASP-69 b more similar in chemical enrichment to Jupiter
($\sim$5$\times$) than Saturn ($\sim$10$\times$) relative to the Sun (Atreya
et al., 2022).

The brightness temperature spectrum plotted in Figure 5 does not show any
obvious signs of a 7.7 µm SO2 feature, that can be very prominent due to
photochemical production of SO2 (Dyrek et al., 2024). The Two Region model
fits the data well in terms of both Bayesian Evidence and $\tilde{\chi}^{2}$
and this model assumes chemical equilibrium and thus essentially no SO2. The
Cloud Layer retrieval has a long tail in the posterior mixing ratio down 10-10
and the Cloud Layer model retrieves a mixing ratio of 10-8, where the feature
size becomes small. Thus, we do not see the significant abundances of SO2
inferred in the hot Jupiter WASP-39 b of 10-5.2 to 10-5.7 (Rustamkulov et al.,
2022; Alderson et al., 2022). One reason could be the moderately high C/O
ratio of 0.65 to 0.94 inferred for WASP-69 b in the Cloud Layer and Two Region
models, which decreases the expected SO2 abundances from photochemical models
(Tsai et al., 2022).

Our observations are potentially sensitive to absorption features from CH4
near 3.3 µm and 7.8 µm, but there are no such features in the emission
spectrum of WASP-69 b, as shown in Figure 5. This is well modeled in our
Scattering, Cloud Layer and Two-Region models, which assume chemical
equilibrium. At the high inferred temperatures of these models ($\gtrsim$1000
K at 10-3 bar), the CH4 abundance is low enough that H2O absorption dominates
and obscures any CH4 features. Note in Figure 5 (second panel from the bottom)
that CH4 doesn’t become optically thick until pressures deeper than 10-1 bar,
which is deeper than where H2O becomes optically thick. In cooler planets,
such as WASP-80 b (850 K), methane presents a strong feature that was observed
in absorption (Bell et al., 2023c). If the terminator and nightside of WASP-69
b is significantly cooler than the dayside, as inferred by our retrievals, CH4
will become a more significant absorber in WASP-69 b’s transmission spectrum.
Our inferred moderately large C/O ratio of 0.65 to 0.94 from the Cloud Layer
and Two Region models would significantly enhance the CH4 abundance in
chemical equilibrium. Continued analysis of existing spectra (Guilluy et al.,
2022; Khalafinejad et al., 2021) and future JWST measurements (e.g. GO
programs 3712 and 5924) will provide better insight into the chemical
processes shaping WASP-69 b’s atmosphere and CH4, now that the dayside
abundances have been constrained in this work.

One other consideration in our abundance inferences from atmospheric
retrievals is that there may be systematic offsets in the JWST emission
spectra. We do not expect significant non-linearity effects in our data
because the difference between the in-eclipse and the out-of-eclipse (0.1% to
0.25% in this case for MIRI, 0.02 to 0.08% for NIRCam) do not appreciably
change the detector well filling factor. However, we did experiment with a
retrieval that allows for the MIRI LRS data spectrum to have an offset
relative to the NIRCam data. We used the Cloud Layer model with an MIRI offset
with an allowed range of -500 ppm to +500 ppm. This retrieval results in a
bimodal posterior distribution with offsets of -80 ppm and +40 ppm that have
significantly different abundances. The -80 ppm offset resulted in a similar
metallicity as the Cloud Layer, Scattering and Two Region models ([M/H]
$\approx$ 1.0, with a moderately lower C/O ratio (C/O$\approx 0.6$). The
higher MIRI offset resulted in a very extended cloud and an order magnitude
higher metallicity. Therefore, if our MIRI data have a systematic offset at
all wavelengths, the inferred composition of the planet can change.
Fortunately, we do not find significant evidence for offsets in our
independent data reductions shown in Figure 3, with a mean MIRI LRS offset
between Eureka! and tshirt of 2 ppm, despite larger differences at individual
data points.

###  7.3 Aerosol Properties

Both the Cloud Layer and Two-Region Models show evidence for clouds that most
significantly affect the long wavelength radiation, as seen in Figure 10. The
Cloud Layer model has a combination of this infrared-absorbing cloud deck
located at 10-4.5 to 10-6 bar and gaseous deeper emission from 10-1 to 10-3
bar with a homogeneous (but slightly diluted) 1D model. The Two-Region model
instead fits the data with separate temperature-pressure profiles and gray
scattering clouds so that the hot region is dominated by the molecular
features near 10-1 to 10-3 bars, and cloud scattering suppresses the flux at
4.0 µm and 7 to 9 µm. The cold region of the Two-Region model has negligible
emission because of the combination of cold temperature and thermal scattering
from clouds (e.g. Taylor et al., 2021) that strongly suppresses the planet
emission. In either case, the averaged dayside properties are best explained
with the inclusion of clouds, and the clouds help reduce the flux at long
wavelengths. The next best model in terms of $\tilde{\chi}^{2}$ the Scattering
model, also would require aerosols to produce the significant reflection
inferred by the geometric albedo parameter A${}_{G}\approx$0.64. The Raleigh
scattering of molecules only provides a geometric albedo of A${}_{G}\lesssim
0.1$, so this would require reflective aerosols in the atmosphere of WASP-69 b
to achieve the inferred geometric albedo. Thus, in all three of the best
models with a $\tilde{\chi}^{2}$ $\lesssim 2.0$, we infer the presence of
aerosols in the dayside of WASP-69 b.

####  7.3.1 Potential Mineral Cloud Compositions

We now discuss the possible compositions of aerosols that could exist in
WASP-69 b’s dayside, first considering cloud formation from commonly assumed
condensed minerals (e.g. Zhang, 2020; Gao et al., 2021). As shown in Figure 9,
the retrieved temperature-pressure profiles from the Cloud Layer model crosses
the the saturation vapor pressure curves for Na2S and MgSiO3 and the
Scattering model crosses the curve for Na2S within 1 $\sigma$. While the Cloud
Layer model uses MgSiO3 (enstatite) optical properties (Scott & Duley, 1996)
to fit the observed dayside spectrum of WASP-69 b, the inferred cloud that
extends from 10-4.5 to 10-6 bar is far higher than the $\sim$10 bar pressure
where MgSiO3 cloud is expected to form. The situation is the same for other
silicates like quartz (SiO2), which condenses at a similar temperature (Grant
et al., 2023) and forsterite (Mg2SiO4), which condenses at even higher
temperatures (e.g., Morley et al., 2012; Wakeford et al., 2017). This
situation is reminiscent of WASP-107b, for which a recent JWST-MIRI
transmission spectrum indicates the presence of silicate clouds at high
altitudes based on a 9-11 µm feature despite its cool equilibrium temperature
of 740 K. At these low temperatures, silicate clouds would be expected to rain
down to the deep atmosphere (Dyrek et al., 2024). An extreme lofting
mechanism, such as vigorous vertical mixing (Gao & Benneke, 2018; Ormel & Min,
2019) and/or nonspherical cloud particles with high porosity (Ohno et al.,
2020; Samra et al., 2020), would be needed to elevate these particles from the
$\sim 10\leavevmode\nobreak\ {\rm bar}$ to the
10-4.5–10${}^{-6}\leavevmode\nobreak\ {\rm bar}$ pressure levels inferred by
our model retrievals for aerosols in WASP-69 b.

Besides silicate clouds, salt clouds such as KCl have been predicted to form
efficiently thanks to their low surface energy which drastically enhance the
nucleation rate (Gao & Benneke, 2018; Lee et al., 2018; Gao et al., 2020).
However, as shown in Figure 9, the atmospheric temperature profile inferred
from the dayside emission spectrum is too hot to form KCl clouds. Instead, KCl
would have to be formed in the cooler terminator and nightside of the planet
and horizontally transported to the dayside.

As discussed in Section 5.3.3, the Cloud Layer model crosses the saturation
vapor pressure curve of Na2S at a similar pressure level as the inferred
cloud. Nucleation theory suggests that Na2S, ZnS, and MnS nucleate very
inefficiently due to their high surface energy, leading to the predction that
those sulfide clouds are absent in exoplanet atmospheres (Gao & Benneke, 2018;
Gao et al., 2020). However, one should keep in mind that the nucleation theory
involves many uncertainties, such as the contact angles of the nucleating
embryo on the condensation nuclei, and Na2S clouds may form if there are
condensation nuclei suitable for Na2S nucleation (for discussion, see Arfaux &
Lavvas, 2023).

To investigate the possibility of Na2S clouds, we attempted to fit the
spectrum of WASP-69 b with an alternative model with the same configuration as
the Cloud Layer model described in Section 5.2.3 but with Na2S’s optical
properties (Montaner et al., 1979). This Na2S Cloud Layer model required a
temperature-pressure profile with about twice the temperature from a
radiative-convective equilibrium calculation and thus was not a physically
plausible solution. One other challenge for Na2S is that WASP-69 b’s
transmission spectrum shows Na atomic features (Casasayas-Barris et al., 2017;
Langeveld et al., 2022; Khalafinejad et al., 2021), which should be suppressed
for significant levels of Na2S cloud formation. Thus overall, there are a few
candidate cloud particle compositions that could form directly on WASP-69 b’s
dayside at high altitudes but do not nucleate efficiently (Na2S and MnS) and
there are more readily-formed cloud particle candidates (KCl and silicates)
that could form in colder regions or deep in the dayside, respectively, and
transported to the observed dayside pressures. Photochemical haze particles
can also form in conditions like atmosphere of WASP-69 b but are expected to
create a temperature inversion not observed in WASP-69 b’s dayside. We discuss
photochemical hazes in Section 7.3.2.

We next consider atmospheric features of candidate cloud particles that may
explain some of the features in WASP-69 b’s dayside spectrum that are not
fully explained by the Two Region, Scattering and Cloud Layer models.
Inspection of Figure 5 reveals that the MIRI LRS spectrum shows a jump in
brightness temperature at a wavelength of 9.4 µm. This is not the same
wavelength as the known instrument artifact, the shadow-region effect, from
10.6 to 11.8 µm (Bell et al., 2023a), nor is there any evidence for a shadow-
region effect in the lightcurves. Two possibilities for this model-data
disagreement are either an absorption feature from WASP-69 b near 8.6 µm
against a rising continuum or a change in the optical properties of aerosols
at 9.4 µm. Evidence for SiO2 (Quartz) has been found in the hot exoplanet
WASP-17 b from its 8.6 µm absorption feature (Grant et al., 2023). WASP-69 b
also has an apparent 8.6 µm absorption feature compared to the Scattering and
Cloud Layer models as visible in Figure 5. Unfortunately, SiO2 is condensed at
temperatures above our inferred temperature-pressure profile for WASP-69 b,
except at deep pressures more than 10 bar, so it would require extreme
vertical mixing to reach the high altitudes near 10-4.5, as with the enstatite
and forsterite discussed above. Another possibility is that the scattering
properties of a cloud (like in the Two-Region model) change abruptly with
wavelength near 9.4 µm. The single scattering albedo of MgSiO3 also drops off
steeply from near 1.0 at wavelengths of 8 µm. If another type of particle
dropped in its single scattering albedo near 9.4 µm, we expect that it would
fit the spectrum of WASP-69 b well. Thus, we encourage future studies of
particles that can exist in solid form on the dayside temperatures and
pressures of WASP-69 b and also have an 8.6 µm absorption feature or changes
in single scattering albedo near 9.4 µm.

####  7.3.2 A Lack of Photochemical Haze?

It has been suggested that hydrocarbon photochemical hazes tend to produce
temperature inversions in the atmospheres of exoplanets (Morley et al., 2015;
Lavvas & Arfaux, 2021; Arfaux & Lavvas, 2022; Steinrueck et al., 2023). Our
observed emission spectrum disfavors the presence of such a temperature
inversion, at least in the hot region that controls the spectral shape,
because the CO2 feature unambiguously appears in absorption rather than in
emission. We also independently computed an emission spectrum by taking into
account the haze’s radiative feedback using the TP profiles computed by the
EGP code coupled with a two-moment microphysical haze model (Ohno & Kawashima,
2020) and found that the addition of hazes tends to worsen the model fit for
haze optical constants tested. Namely, we used Titan-like haze compositions
(Khare et al., 1984) and soot (Lavvas & Koskinen, 2017) . The lack of
hydrocarbon hazes could be consistent with the apparent lack of CH4 in the
dayside spectrum and several modeling studies that predicted the decline in
the haze abundance at $T_{\rm eq}\gtrsim 950\leavevmode\nobreak\ {\rm K}$ due
to the conversion of CH4 into CO (Morley et al., 2015; Kawashima & Ikoma,
2019; Gao et al., 2020). On the other hand, WASP-69 b is known to show a steep
spectral slope in optical wavelengths (Murgas et al., 2020; Estrela et al.,
2021), and our results may challenge previous studies that attribute the
optical slope to photochemical hazes (e.g., Lavvas & Koskinen, 2017; Ohno &
Kawashima, 2020). However, we note that the impacts of photochemical hazes on
the TP profile should be sensitive to the actual optical properties of
exoplanetary hazes, which have been highly uncertain to date. Further
understanding of the haze optical properties, such as by experimental studies
(Corrales et al., 2023; He et al., 2023), would be needed to quantitatively
investigate to what extent photochemical hazes might exist in the atmosphere
of WASP-69b.

##  8 Conclusions

We analyzed a 2 µm to 12 µm JWST emission spectrum of WASP-69 b and find the
following:

  * •

We observe molecular features in absorption, so the temperature profile
decreases with altitude

  * •

We observe features of of CO2, CO and H2O in the atmosphere, but no features
of CH4.

  * •

The shortest wavelengths (less than 3 µm) observed in NIRCam have unexpectedly
high brightness temperatures above 1000 K as compared to the longest
wavelengths (more than 5 µm)

  * •

There is some marginal evidence of an inhomogeneous dayside for WASP-69 b from
the MIRI broadband eclipse maps.

We model the emission spectrum with a variety of assumptions about the
atmosphere and find the following:

  * •

Cloudless 1D homogeneous chemical equilibrium models cannot fit our observed
spectrum well and imply that additional complexity is needed to explain
WASP-69 b’s dayside spectrum.

  * •

Aerosols are present in WASP-69 b’s dayside atmosphere because they are needed
in our models to fit the shape of the planet’s 2 µm to 12 µm spectrum,
especially the high brightness temperature at short wavelengths as compared to
long wavelengths and the lack of a strong peak in brightness near 4.0 µm.

  * •

We find three models that can explain WASP-69 b’s spectrum 1) a Scattering
Model with a free parameter for the geometric albedo, possibly due to aerosols
2) A Cloud Layer Model that has a high altitude silicate cloud layer and 3) A
Two-Region model that has an inhomogeneous dayside that is a combination warm
and cool components that both contain gray aerosols. The Scattering model
requires unexpectedly high albedos of 0.64, which is less plausible.

  * •

The abundances of CO2, CO and H2O point to super-solar concentrations of heavy
elements in the atmosphere. Considering the close agrement between the Cloud
Layer and Two-Region Models, the metallicity of the atmosphere is enriched by
6$\times$ to 14$\times$ compared to solar composition and the C/O ratio is
0.65 to 0.94, but we cannot completely rule out the Scattering model, which
has an atmospheric metallicity of 4$\times$ to 8$\times$ and a C/O ratio of