11institutetext:  Department of Astronomy, Stockholm University, AlbaNova
University Center, 10691 Stockholm, Sweden 22institutetext:  Weltraumforschung
und Planetologie, Physikalisches Institut, Universität Bern,
Gesellschaftsstrasse 6, 3012 Bern, Switzerland 33institutetext:  Center for
Space and Habitability, Universität Bern, Gesellschaftsstrasse 6, 3012 Bern,
Switzerland 44institutetext:  Observatoire astronomique de l’Université de
Genève, Chemin Pegasi 51, 1290 Versoix, Switzerland 55institutetext:  Jet
Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91011,
USA 66institutetext:  Faculty of Physics, Ludwig Maximilian University,
Scheinerstrasse 1, D-81679, Munich, Bavaria, Germany 77institutetext:  ARTORG
Center for Biomedical Engineering Research, University of Bern, Murtenstrasse
50, CH-3008, Bern, Switzerland 88institutetext:  University College London,
Department of Physics & Astronomy, Gower St, London, WC1E 6BT, United Kingdom
99institutetext:  University of Warwick, Department of Physics, Astronomy &
Astrophysics Group, Coventry CV4 7AL, United Kingdom 1010institutetext:
Department of Space, Earth and Environment, Chalmers University of Technology,
Onsala Space Observatory, 439 92 Onsala, Sweden 1111institutetext:  Department
of Astronomy and Astrophysics, The University of Chicago, Chicago, IL 60637,
USA 1212institutetext:  Centre Vie dans l’Univers, Faculté des sciences de
l’Université de Genève, Quai Ernest-Ansermet 30, 1205 Geneva, Switzerland
1313institutetext:  Leiden Observatory, University of Leiden, PO Box 9513,
2300 RA Leiden, The Netherlands 1414institutetext:  Division of Geological and
Planetary Sciences, California Institute of Technology, Pasadena, CA 91125,
USA 1515institutetext:  Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3
0HE, UK

# JWST reveals a rapid and strong day side variability of 55 Cancri e

J. A. Patel${}^{\href https://orcid.org/0000-0001-5644-6624}$ 11
<EMAIL_ADDRESS>A. Brandeker 11 D. Kitzmann${}^{\href
https://orcid.org/0000-0003-4269-3311}$ 2233 D. J. M. Petit dit de la Roche 44
A. Bello-Arufe 55 K. Heng 66778899 E. Meier Valdés 33 C. M. Persson 1010 M.
Zhang 1111 B.-O. Demory 3322 V. Bourrier 44 A. Deline 44 D. Ehrenreich 441212
M. Fridlund 13131010 R. Hu${}^{\href https://orcid.org/0000-0003-2215-8485}$
551414 M. Lendl 44 A. V. Oza 141455 Y. Alibert 2233 M. J. Hooton 1515

(Received XX; accepted XX)

###### Abstract

Context. The nature of the close-in rocky planet 55 Cnc e is puzzling despite
having been observed extensively. Its optical and infrared occultation depths
show temporal variability, in addition to a phase curve variability observed
in the optical.

Aims. We wish to explore the possibility that the variability originates from
the planet being in a 3:2 spin-orbit resonance, thus showing different sides
during occultations. We proposed and were awarded Cycle 1 time at the James
Webb Space Telescope (JWST) to test this hypothesis.

Methods. JWST/NIRCam observed five occultations (secondary eclipses), of which
four were observed within a week, of the planet simultaneously at 2.1 and 4.5
$\mu$m. While the former gives band-integrated photometry, the latter provides
a spectrum between 3.9–5.0 $\mu$m.

Results. We find that the occultation depths in both bandpasses are highly
variable and change between a non-detection ($-5\pm 6$ ppm and $7\pm 9$ ppm)
to $96\pm 8$ ppm and $119^{+34}_{-19}$ ppm at 2.1 $\mu$m and 4.5 $\mu$m,
respectively. Interestingly, the variations in both bandpasses are not
correlated and do not support the 3:2 spin-orbit resonance explanation. The
measured brightness temperature at 4.5 $\mu$m varies between 873–2256 K and is
lower than the expected dayside temperature of bare rock with no heat re-
distribution (2500 K) which is indicative of an atmosphere. Our atmospheric
retrieval analysis of occultation depth spectra at 4.5 $\mu$m finds that
different visits statistically favour various atmospheric scenarios including
a thin outgassed CO/CO2 atmosphere and a silicate rock vapour atmosphere. Some
visits even support a flat line model.

Conclusions. The observed variability could be explained by stochastic
outgassing of CO/CO2, which is also hinted by retrievals. Alternatively, the
variability, observed at both 2.1 and 4.5 $\mu$m, could be the result of a
circumstellar patchy dust torus generated by volcanism on the planet.

###### Key Words.:

techniques: spectroscopic – techniques: photometric – planets and satellites:
atmospheres – planets and satellites: terrestrial planets – planets and
satellites: individual: 55 Cnc e

##  1 Introduction

Ultra-short-period planets (USPs) provide a unique opportunity to study
planets in extreme environments that have no counterparts in our Solar system
(see Winn et al., 2018, for a review). Many USPs are consistent with a bare
rock composition, while some of them might have a secondary metal-rich
atmosphere or a disintegrating surface (e.g., Brogi et al., 2012; Kreidberg et
al., 2019; Zieba et al., 2022). Being in an orbit around the nearby ($d=12.6$
pc), bright naked eye star 55 Cancri ($V=5.95$ mag), 55 Cancri e (hereafter 55
Cnc e) is one of the best targets for investigating the nature of a USP. Out
of the five known planets in the system, planet e is the only one transiting
the star.

55 Cnc e was discovered by McArthur et al. (2004) with an orbital period of
$\sim 2.8$ d, which was later found to be an alias of the true $0.74$ d period
(Dawson & Fabrycky, 2010). This was confirmed by the detection of planetary
transits in the optical and infrared (IR) independently (Winn et al., 2011;
Demory et al., 2011), enabling its radius measurement. Together with mass
estimates derived from radial velocity measurements, the earlier works
attempted to constrain the internal structure of the planet and found that the
planetary density was consistent with either a purely rocky planet, a rocky
planet with a thick super-critical water envelope, or carbon-rich interior
with no envelope (Demory et al., 2011; Gillon et al., 2012; Madhusudhan et
al., 2012). More recently, Bourrier et al. (2018a) refined planetary mass
($8.3\,M_{\oplus}$) and radius ($1.88\,R_{\oplus}$) using radial velocity data
and HST/STIS transit observations. Their internal structure modelling, based
on these updated mass-radius measurements, suggests a rocky planet surrounded
by a heavyweight (high mean molecular weight) atmosphere. A low mean molecular
weight, or lightweight, atmosphere on the planet is not possible because of
intense radiation from its host star. Atmospheric escape simulations also
imply that lightweight atmospheres (made of H, He) would not survive on 55 Cnc
e for a long time period (e.g., Gillon et al., 2012; Salz et al., 2016;
Bourrier et al., 2018a; Zhang et al., 2021). Other attempts to model the
internal structure of the planet (e.g., Dorn et al., 2017; Lopez, 2017; Crida
et al., 2018) indicate a rocky interior with a gas or water envelope.

Soon after its discovery, Demory et al. (2012) used Spitzer to detect thermal
emission from 55 Cnc e and determined its day-side temperature to be around
2300 K. Demory et al. (2016a) constructed a temperature map of the planet
using Spitzer/IRAC phase curve measurements at 4.5 $\mu$m. They calculated the
average day-side temperature to be around 2350 K with a maximum $\sim 2700$ K.
Curiously, the hottest location of the planet was found to be shifted by $\sim
41^{\circ}$ to the east compared to the sub-stellar point, indicating a strong
heat re-distribution. On the other hand, the day-night temperature difference
was found to be as large as 1300 K, a sign of inefficient heat transport to
the night side. These conflicting results led Demory et al. (2016a) to
speculate that perhaps efficient heat transport is only happening on the day
side of the planet by a thick atmosphere, or alternatively, a molten lava flow
is responsible for the heat transport. The inefficiency of energy transport to
the night side could be due to gases becoming cold enough to condense.
Similarly, a lava stream could be hindered by the surface solidifying at the
night side. Angelo & Hu (2017) re-analysed the phase-curve data and confirmed
the findings of Demory et al. (2016a). Their physical model of the phase curve
allowed them to show that the radiative and advective time scales must be of
the same order to reproduce the observed phase curve. This disfavours the lava
ocean scenario since a lava flow would have a too large advective time scale
(e.g., Kite et al., 2016) to be an efficient heat transporter (however
interior dynamics models of the planet, in some cases, exhibits a mantle
super-plume away from the sub-stellar point, which can potentially interact
with the lava ocean and increase its temperature at the location of the plume,
mimicking hot-spot offset; Meier et al., 2023). Angelo & Hu (2017) further
propose that a CO or N2 dominated atmosphere on the day side could explain the
phase curve. This claim was corroborated by a 3D global circulation model
climate model by Hammond & Pierrehumbert (2017) that could potentially
describe the observations, assuming a H2 \+ N2 dominated atmosphere with a
trace source of opacity at 4.5 $\mu$m (such as CO2 or H2O), coupled with the
with a presence of night-side clouds. A recent re-reduction and re-analysis of
the Spitzer phase curve by Mercier et al. (2022) yielded an even larger day-
night temperature difference with a smaller phase offset, more consistent with
a poor heat transport typically found on USPs.

The heavyweight atmosphere on the planet, which was implied by the Spitzer
phase curve, climate modelling and mass-radius constraints, is challenging to
detect. Numerous observations have tried but failed to detect any atmosphere
on the planet The singular claim of detection of gas on 55 Cnc e comes from
Tsiaras et al. (2016) who identified HCN in the atmosphere using HST/WFC3
transit observations. However, subsequent observations using high-resolution
spectroscopy from the ground could not reproduce the detection of HCN (Deibert
et al., 2021). Furthermore, the transit observation of 55 Cnc e in the Ly
$\alpha$ band by Ehrenreich et al. (2012) resulted in a non-detection,
suggesting the absence of an extended H upper atmosphere. This was supported
by the non-detection of He in the upper atmosphere by Zhang et al. (2021). A
lack of H and He in the atmosphere could mean that both gases escaped if they
were initially accreted from the disk. In addition to this, several studies
attempted but could not detect other atmospheric species such as H2O, TiO,
NH3, C2H2, Fe, Ca, Mg, K, Na, H (Ridden-Harper et al., 2016; Esteves et al.,
2017; Jindal et al., 2020; Tabernero et al., 2020; Deibert et al., 2021; Keles
et al., 2022; Rasmussen et al., 2023). These non-detections mean that either
those species are absent from the atmosphere or only present at very low
volume mixing ratios if the mean molecular weight of the atmosphere is not too
high to be detected by the transit observations. Another possibility is that
the atmosphere of the planet is cloudy (Mahapatra et al., 2017).

The IR observations of 55 Cnc e in emission posed another challenge for
understanding the behaviour of the planet. Demory et al. (2016b) monitored the
occultation depths of 55 Cnc e with Spitzer at 4.5 $\mu$m during 2012–2013 and
found a variable occultation depth ranging from 47 ppm to 176 ppm. This
translates into a corresponding change of brightness temperature of 1370 K to
2530 K. Variability was also observed in the optical bandpass by MOST
(Microvariability and Oscillations of STars) that discovered significant
changes in phase curves over several seasons (Winn et al., 2011; Dragomir et
al., 2014; Sulis et al., 2019). While the optical observations with MOST found
a significant phase curve amplitude, the secondary occultation remained
undetected. More recently, CHEOPS (CHaracterising ExOPlanet Satellite)
extensively observed 55 Cnc (Morris et al., 2021; Demory et al., 2023; Meier
Valdés et al., 2023) in the optical (G band) and confirmed significant
variability not only in phase amplitude but also in phase offset and
occultation depth, where the occultation depths at some epochs were consistent
with zero. TESS (the Transiting Exoplanet Survey Satellite) also observed 55
Cnc and found a hint of weak variability in occultation depths over three
observing sectors (Meier Valdés et al., 2022). In contrast to the variability
of the occultation depths, no optical or IR variability has been observed in
the transit depths (e.g. Meier Valdés et al., 2023; Bourrier et al., 2018a).

Multiple studies in the literature propose various hypotheses to explain the
observed variability of the occultation depth of 55 Cnc e in the optical and
IR. Demory et al. (2016b) suggested that plumes from volcanic outgassing on
the dayside could explain the observed variability in emission. Assuming an
Earth-like composition for the interior, it can release gases such as CO or
CO2 that are a significant source of opacity around 4.5 $\mu$m. Gas plumes
evolving at different atmospheric pressure levels could be inferred as varying
temperatures during occultation observations in the IR. Given that the
variability was observed throughout the optical and IR, it was suggested by
Morris et al. (2021) that a circumstellar inhomogeneous dusty torus could
provide a variable source of opacity. Meier Valdés et al. (2023) studied the
dusty torus scenario in detail and concluded that such a torus made up of
certain species of a narrow range of particle sizes could indeed reproduce the
level of observed variability in the optical. However, a dusty torus should
extent out to its Hill sphere and, if opaque, is inconsistent with the
observed transit depths (Heng, 2023). Heng (2023) argued that a thin,
transient outgassed atmosphere is consistent not only with the observed
optical and infrared occultation depths, but also provides a plausible
explanation for their variability. Tian & Heng (2024) demonstrated that CO-CO2
atmospheres are outgassed under a broad range of conditions (surface
pressures, oxygen fugacity, temperatures).

Since 55 Cnc e is in a very close-in orbit around its host star, Folsom et al.
(2020) showed that the planet’s orbit is inside the stellar Alfvén surface.
This means that star-planet interactions (SPI) are plausible for the system,
potentially causing variability-inducing star spots. Bourrier et al. (2018b)
proposed coronal rain, a kind of SPI, as a reason for the variability in
chromospheric lines that they observed with HST (see also Sulis et al., 2019).
Morris et al. (2021) ruled out star spot creation by the planet as a plausible
mechanism to explain the optical variability observed by CHEOPS but this does
not prohibit other possible forms of SPIs, such as coronal rain.

Although multiple hypotheses have been provided to describe the thermal phase
curve and variability from 55 Cnc e, each has difficulties in fully explaining
all observed features. The observations with the James Webb Space Telescope
(JWST), presented here, were in part motivated by exploring an alternative
hypothesis that the planet rotates at an asynchronous rate to its orbit,
potentially explaining both the hot-spot shift into the afternoon and the
rapid orbit-to-orbit variability. The idea and the observations motivated by
it are presented in Section 2 followed by results in Section 3. We show the
results from atmospheric retrieval analysis in Section 3.2. Finally, we
interpret the results from our observations and present our conclusions in
Section 4 and 5, respectively. Details of the data analysis methods used are
put into Appendix A.

Table 1: Observation log and wide band occultation depths Visit | Prog. | Start date | End date | Parity | Occultation | Occultation | Brightness | Brightness
---|---|---|---|---|---|---|---|---
| ID |  |  |  | depth at | depth at | temp. at | temp. at
|  |  |  |  | 2.1 $\mu$m (ppm) | 4.5 $\mu$m (ppm) | 2.1 $\mu$m (K) | 4.5 $\mu$m (K)
1 | 2084 | 2022-11-18 14:40:17 | 2022-11-18 19:15:53 | even | $47.4^{+21.0}_{-15.5}$ | $7.0^{+8.8}_{-8.8}$ | $2417^{+335}_{-287}$ | $873^{+167}_{-187}$
2 | 2084 | 2022-11-20 19:43:08 | 2022-11-21 00:18:44 | odd | $-5.1^{+5.5}_{-6.0}$ | $65.2^{+22.3}_{-42.2}$ | $1247^{+190}_{-245}$ | $1716^{+230}_{-315}$
3 | 2084 | 2022-11-23 00:43:57 | 2022-11-23 05:19:33 | even | $37.3^{+4.7}_{-4.6}$ | $101.4^{+17.1}_{-32.4}$ | $2234^{+86}_{-88}$ | $2078^{+172}_{-342}$
4 | 1952 | 2022-11-24 11:38:15 | 2022-11-24 17:28:41 | even | $36.8^{+27.7}_{-32.9}$ | $119.2^{+34.0}_{-19.0}$ | $2302^{+413}_{-807}$ | $2256^{+330}_{-188}$
5 | 2084 | 2023-04-24 11:57:03 | 2023-04-24 16:32:36 | odd | $95.9^{+8.1}_{-7.9}$ | $95.4^{+13.5}_{-16.8}$ | $3138^{+107}_{-107}$ | $2016^{+137}_{-179}$

##  2 Asynchronous rotation scenario for 55 Cnc e, Observations and Methods

###  2.1 55 Cnc e in 3:2 spin-orbit resonance

The planet 55 Cnc e orbits its host star in about 17.7 h with a semi-major
axis of 0.015 AU (Bourrier et al., 2018a). When a planet is orbiting this
close to its host star, it is usually assumed to be in a tidally locked
synchronous spin-orbit configuration because of strong tidal forces. However,
if the planet is part of a multi-planetary system, gravitational interactions
with the other planets can perturb the planet from its synchronous 1:1 spin-
orbit configuration. Rodríguez et al. (2012) simulated the tidal evolution of
the orbit of 55 Cnc e and showed that there is a reasonable likelihood for the
planet to be trapped in an asynchronous spin-orbit resonance, with the 3:2
spin-orbit resonance being the most likely after 1:1 synchronous rotation.
Asynchronous rotation can thus not be ruled out for 55 Cnc e. The consequence
is that the planet would show different faces to the star during the orbit.
This in turn means that the hottest point on the planet would not necessarily
be the sub-stellar point. Just as on Earth the hottest time of the day is in
the afternoon and not at noon, so could thermal inertia on 55 Cnc e shift its
hottest spot to the afternoon (east). The thermal inertia could, like on
Earth, be provided by the atmosphere. In the case of a bare rock, thermal
inertia could be provided by the heating, melting and evaporation of the rock
in the morning with subsequent condensation and crystallisation in the
afternoon. Quantitative models of these scenarios are sensitive to the
detailed assumptions of the mass and composition of the atmosphere that, in
turn, depends on the material equation of state. Using simplified models,
Brandeker (2019) showed that the observations up until then could indeed be
explained by using reasonable assumptions on the physical properties of the
planet, meaning that the asynchronous rotation scenario could not be excluded.

Assuming that the planet is rotating asynchronously in the most probable 3:2
spin-orbit resonance, the planet will show the same face only at every second
occultation instead of showing the same face every time. That means the two
opposite sides will be seen during consecutive occultations. If there are
semi-stable surface features, e.g., due to volcanic activity, on different
sides of the planet they will show up differently during alternate
occultations. In this case, the observed occultation depths would be expected
to highly correlate with the occultation number over a short period, while
this correlation could be broken over a longer time scale due to surface
changes. The variability in occultation depths observed by Demory et al.
(2016b) can then be attributed simply to the planet showing different faces
during occultations. Notably, Tamburo et al. (2018), who confirmed the Spitzer
variability of occultation depths, found the variability to be well fitted by
a sinusoidal with a period as short as 2 days, but discarded this solution as
being unphysical. However, if the planet is indeed in a 3:2 spin-orbit
resonance, it is expected that the period of variability should be equivalent
to the synodic period ($\sim 35.5$ hr), close to the period of 2 days.

To further test this intriguing hypothesis of asynchronous rotation and
simultaneously sensitively measure potential atmospheric signatures, we
designed an observation programme for JWST as detailed in the next section.

###  2.2 Observations

If the planet is indeed in a 3:2 spin-orbit resonance, it will show two
opposite sides in consecutive occultations. Assuming the planetary surface to
evolve slowly, we would then expect every second consecutive occultation to be
strongly correlated. Enumerating the occultations by orbit number, we thus
requested two “odd” and two “even” occultations within a short time-
constrained span of two weeks, to rule out significant surface evolution
within that time. Since 55 Cnc is a very bright IR target ($K=4$ mag),
avoiding saturation while observing it with JWST is challenging. From pre-
launch estimates, our options were essentially limited to a grism time-series
mode of the Near Infrared Camera (NIRCam). The proposal was awarded time in
JWST Cycle 1 as GO 2084 (Brandeker et al., 2021). The observation log is
provided in Table 1. Due to technical difficulties, only three occultations of
the programme were observed within the time constraint of two weeks; the
fourth was postponed until five months later. Fortunately, a different
programme (GO 1952, Hu et al., 2021) that also targeted 55 Cnc had an
occultation observed in the same instrument mode and within the same first
week (Hu et al., 2024). In the following, we thus present an analysis of all
five visits.

NIRCam offers simultaneous observations in short-wave (SW) and long-wave (LW)
channels at 0.6–2.3 $\mu$m and 2.4–5.0 $\mu$m, respectively. The SW channel
allows the use of a weak lens with a filter providing photometric monitoring
of the target, while the LW channel provides a spectroscopic mode using a
grism and a filter. Our observations in the LW channel used the F444W filter
with GRISMR element and RAPID readout mode. On the other hand, the WLP4/F212N2
weak lens/filter with RAPID readout mode was used in the SW channel. Both
channels employed the SUBGRISM64 subarray that has 2048 columns and 64 rows.
This gave us spectroscopic data between 3.9–5 $\mu$m (centred at around 4.5
$\mu$m) in the LW channel (or, 4.5 $\mu$m channel) and one single photometric
data point in a narrow-band (2.3%) bandpass at 2.12 $\mu$m from the SW channel
(also referred to as the 2.1 $\mu$m channel). Given the brightness of the host
star, we chose two groups per integration with a total integration time of
about 1.03 s.

We used five independent pipelines to reduce and analyse the spectroscopic
data at 4.5 $\mu$m and two different pipelines to analyse the short-wave
photometric data. The details of these methods are described in Appendix A.

###  2.3 Retrieval model and atmospheric scenarios

We chose two representative independent reductions of occultation depth
spectra, from stark and HANSOLO pipelines, to perform atmospheric retrieval.
Both reductions differ in their treatment of correlated noise and thus produce
slightly different results, which was the reason for choosing two different
reductions for retrieval (see Appendix A for more details).

To interpret the observational data, we used the open-source HELIOS-r2
atmospheric retrieval code (Kitzmann et al., 2020), which uses the nested
sampling algorithm (Skilling, 2004) implemented in the MultiNest library
(Feroz & Hobson, 2008). For the atmospheric characterisation, we tested four
different models with a varying level of complexity. The simplest model tries
to fit the observational data with a flat line, while the second one assumes
the planet to emit like a pure blackbody of temperature $T_{\mathrm{bb}}$.
Since observations by, for example, Ehrenreich et al. (2012) and Zhang et al.
(2021) rule out the presence of a thick primordial hydrogen-helium atmosphere,
a potential atmosphere has to be secondary in nature. There are two essential
pathways to create a secondary atmosphere for a hot planet such as 55 Cnc e.
The atmosphere can either be dominated by outgassing from the planetary
interior (e.g., Tian & Heng, 2024) or be created through evaporation of mantle
material, or a combination thereof. Thus, for the two atmospheric scenarios,
we assumed a secondary atmosphere with outgassed carbon monoxide (CO)/carbon
dioxide (CO2) (e.g. Heng, 2023) or an atmosphere produced by an evaporating
mantle with a bulk silicate earth composition that is composed of silicon
oxide (SiO), silicon dioxide (SiO2), and magnesium oxide (MgO) (Zilinskas et
al., 2022).

Nested sampling allows Occam’s Razor (of Ockham, 1495) to be enforced via the
calculation of the Bayesian evidence (or marginalised likelihood function),
see Trotta (2008, 2017). In practice, this allows us to favour simpler
explanations for some of the data (e.g. flat line or blackbody function). To
provide good constraints on the Bayesian evidence values, within MultiNest we
used 5000 live points (Feroz & Hobson, 2008) for each retrieval calculation.
Increasing this value further did not alter the resulting evidence values to a
significant degree.

The atmosphere was considered to be isothermal with the surface pressure
$p_{\mathrm{surf}}$ as a free parameter in the retrieval model. The atmosphere
and surface were allowed to have their own distinct temperatures,
$T_{\mathrm{atm}}$ and $T_{\mathrm{surf}}$, respectively.

The cross sections of CO, CO2, SiO, SiO2, and MgO were taken from Li et al.
(2015), Yurchenko et al. (2020), Yurchenko et al. (2022), Owens et al. (2020),
and Li et al. (2019), respectively. All temperature and pressure-dependent
cross sections were calculated with the open-source opacity calculator
HELIOS-K (Grimm & Heng, 2015; Grimm et al., 2021).

Table 2: Summary of retrieval parameters and prior distributions used for the retrieval models. Parameter | Prior
---|---
| Type | Value
Flat line |  |
Occultation depth | uniform |  $0$ ppm – $200$ ppm
Blackbody |  |
$d_{wl}$ | Gaussian | see Table 1
$R_{p}/R_{*}$ | Gaussian | $0.0182\pm 0.0002$
$T_{\mathrm{bb}}$ | uniform |  $300$ K – $3000$ K
Atmosphere |  |
$d_{\mathrm{wl}}$ | Gaussian | see Table 1
$R_{p}/R_{*}$ | Gaussian | $0.0182\pm 0.0002$
$p_{\mathrm{surf}}$ | log-uniform |  $10^{-10}$ bar – $500$ bar
$T_{\mathrm{surf}}$ | uniform | 300 K – 3000 K
$T_{\mathrm{atm}}$ | uniform | 300 K – 3000 K
$\xi_{j}$ | uniform | $10^{-10}\leq x_{j}\leq 1$
Figure 1: Detrended occultation light curves from the short-wave photometric
channel (2.1 $\mu$m, left panel) and white-light light curves from the long-
wave channel (4.5 $\mu$m, right panel). Only binned data points are shown
here. The darker and lighter shades of the points depict even and odd orbital
number parity, respectively. The dates and occultation depth (median and
68-percentile confidence intervals) of the visits are indicated above each
plot. The best-fitted models and models computed from randomly selected
posteriors to show the model uncertainties are plotted with thick and thin
lines. Figure 2: Observed wide band occultation depths in long-wave (in
orange) and in short-wave (in blue) channels. The depths are plotted as a
function of epoch number starting from the first visit.

The atmospheric composition in the retrieval model was described through a
centred-log-ratio prior which allows a more optimised sampling of the
parameter space when the dominant background gas is not known (Benneke &
Seager, 2012). For a given mixture of $n$ gases, the centred-log-ratio
conversion (clr) for the mixing ration $x_{j}$ of a given molecule $j$ in the
mixture is given by

$\xi_{j}=\mathrm{clr}(x_{j})=\ln\frac{x_{j}}{g(\mathbf{x})}\ ,$ (1)

where $g(\mathbf{x})$ is the geometric mean of all mixing ratios $\mathbf{x}$:

$g(\mathbf{x})=\left(\prod_{j=1}^{n}x_{j}\right)^{1/n}\ .$ (2)

Due to the constraint that

$\sum_{j=1}^{n}x_{j}=1\quad\text{or}\quad\sum_{j=1}^{n}\xi_{j}=0\ ,$ (3)

only $n-1$ free parameters are needed in the retrieval. We used uniform priors
to produce $\xi_{j}$ values subject to the constraints that
$\min\left(\mathbf{x}\right)=10^{-10}$ and $\max\left(\mathbf{x}\right)=1$,
see Benneke & Seager (2012) for details. We note that the prior boundaries for
$\xi_{j}$ depend on the number of molecules in the retrieval and the chosen
value of the smallest allowed mixing ratio.

For the retrieval of the data from the stark reduction, we performed the
calculations on the relative occultation depths. Thus, for these calculations,
we needed to add an additional free parameter to the retrieval: the white-
light occultation depth $d_{\mathrm{wl}}$. For these, we used Gaussian priors
with the values provided in Table 1. Since HANSOLO reduction provides absolute
occultation depths this additional parameter was not needed. Additionally, we
binned the data provided by stark which uses the instrument’s native
resolution to about 30 spectral bins.

All retrieval parameters for the different models are summarised in Table 2.
The empirically calibrated stellar spectrum of 55 Cnc from Crossfield (2012)
was used to transform the emission spectra calculated by the retrieval model
to wavelength-dependent occultation depths.

##  3 Results

###  3.1 Wide band occultation depths

Table 3: Overview of retrieval results for the stark and HANSOLO reductions. Boldface indicates the statistically preferred models. Model | Visit 1 | Visit 2 | Visit 3 | Visit 4 | Visit 5
---|---|---|---|---|---
| $\ln\mathcal{Z}$ | $B$ | $\ln\mathcal{Z}$ | $B$ | $\ln\mathcal{Z}$ | $B$ | $\ln\mathcal{Z}$ | $B$ | $\ln\mathcal{Z}$ | $B$
stark |  |  |  |  |  |  |  |  |  |
Flat line | $-169.98$ | $e^{32.3}$ | $-148.08$ | $e^{14.1}$ | $-154.10$ | $e^{16.2}$ | $-147.70$ | $e^{11.8}$ | $\mathbf{-135.03}$ | -
Blackbody | $-159.53$ | $e^{21.8}$ | $-134.26$ | 1.3 | $-154.33$ | $e^{16.4}$ | $\mathbf{-135.90}$ | - | $-140.10$ | 159.9
CO, CO2 | $\mathbf{-137.72}$ | - | $\mathbf{-133.96}$ | - | $-147.66$ | $e^{9.8}$ | $-135.96$ | 1.1 | $-137.48$ | 11.6
SiO, SiO2, MgO | $-139.56$ | 6.3 | $-135.17$ | 3.4 | $\mathbf{-137.90}$ | - | $-136.71$ | 2.3 | $-141.01$ | $e^{6.0}$
HANSOLO |  |  |  |  |  |  |  |  |  |
Flat line | $-115.19$ | $9.5$ | $-109.72$ | 12.2 | $-143.00$ | 27.2 | $-129.95$ | 1.7 | $-134.64$ | 1.3
Blackbody | $\mathbf{-112.94}$ | - | $\mathbf{-107.22}$ | - | $\mathbf{-139.68}$ | - | $\mathbf{-129.41}$ | - | $\mathbf{-134.41}$ | -
CO, CO2 | $-113.66$ | 2.1 | $-108.23$ | 2.7 | $-139.72$ | 1.0 | $-130.10$ | 2.0 | $-134.97$ | 1.8
SiO, SiO2, MgO | $-114.06$ | 3.0 | $-108.35$ | 3.1 | $-140.39$ | 2.2 | $-130.36$ | 2.6 | $-135.43$ | 2.5

We used six pipelines to reduce and fit our JWST/NIRCam dataset. The methods
are described in detail in Appendix A. Here we present results from our
primary analysis from the stark pipeline (Appendix A.1). A summary of our
results, along with the observation log, is tabulated in Table 1.

Our main finding is the strong variability in occultation depths. The white-
light occultation depths (computed by fitting an occultation model to the
band-averaged occultation time series) at 4.5 $\mu$m are highly variable even
during the short time scale of a week (Table 1). During the time span of 6
days (8 planetary orbits), the measured occultation depths at 4.5 $\mu$m
continuously increased from basically non-detection in Visit 1 ($7\pm 9$ ppm)
to $119^{+34}_{-19}$ ppm in Visit 4. The occultation depth from our final
visit (Visit 5), observed 5 months after the other visits, is $\sim 95\pm 16$
ppm and consistent with the depths from Visits 3 and 4 but differs
significantly from the depths from Visit 1 and 2. Fig. 2 shows occultation
depths as a function of time, illustrating this point. The best-fitted
occultation models along with the de-trended data are shown in Fig. 1 for all
visits.

We used an empirically calibrated stellar spectrum of 55 Cnc from Crossfield
(2012), stellar and planetary parameters from Bourrier et al. (2018a) and the
NIRCam response function111http://svo2.cab.inta-csic.es/theory/fps/ to compute
brightness temperatures using the measured white-light occultation depths at
4.5 $\mu$m. As shown in Table 1, the brightness temperature changes
significantly from 873 K to 2256 K within a week. Notably, the brightness
temperature almost doubled from Visit 1 to 2, i.e., only after three planetary
orbits.

Similarly, the 2.1 $\mu$m channel occultation depths are also variable. Within
a week, the 2.1 $\mu$m occultation depths remained almost constant at around
40 ppm for Visit 1, 3 and 4, while we found a non-detection of occultation for
Visit 2 which was observed between Visit 1 and 3 (see, Fig. 2). However, the
final observation that was taken 5 months later (Visit 5) shows a
significantly higher occultation depth of $96\pm 8$ ppm, which is almost equal
to the depth observed at 4.5 $\mu$m in the same epoch. The corresponding
brightness temperatures varies significantly between 1247 K and 3138 K (see,
Table 1). Interestingly, there is no correlation between the occultation depth
variability observed at 2.1 and 4.5 $\mu$m (Fig. 2). Fig. 1 present the de-
trended SW data with best-fitted models.

The variability, plotted in Fig. 2, is clearly not correlated with the parity
of the orbit number. Occultation depths are also variable between occultations
from orbits of the same parity, e.g., in even (Visit 1, 3, 4) or odd (Visit 2,
5) visits. The rapid variability thus cannot be explained by simply
alternating between two sides of the planet. This does not rule out that the
planet rotates asynchronously but means that an explanation for the rapid
variability has to be found elsewhere.

All visits showed various degrees of significant correlated noise of unknown
origins, in both the 2.1 and 4.5 $\mu$m channels. The leftover correlated
noise can be seen in Fig. 1 and are also quantified in the Allan deviation
plots in Fig. 13. We perform an injection-retrieval test to estimate proper
uncertainties on occultation depths in the presence of correlated noise (see,
Sec. A.1.1). We report uncertainties from this analysis in Table 1. We,
however, found that various methods to account for correlated noise could
somewhat change the results of occultation depths and emission spectra (see,
Appendix A for more details).

###  3.2 Occultation depth spectra and atmospheric retrieval

We computed the relative occultation depth spectra as outlined in Appendix A.1
using the stark reduction and absolute occultation depth spectra from HANSOLO
pipeline as described in Appendix A.4. Since different methods to handle the
correlated noise could lead to different results, we chose to perform
atmospheric retrieval analysis on results from two pipelines, stark and
HANSOLO, which use two representative techniques to deal with the correlated
noise (see, Appendix A for details). The occultation depth spectra, shown in
Fig. 11, are also variable from visit to visit and do not show any consistent
spectral features.

Figure 3: Posterior spectra for all model scenarios and visits. The left
column shows predicted occultation depths in the shortwave channel. The black
data points indicate the observed value, while diamonds represent the
retrieval results for the HANSOLO reduction and squares refer to the outcome
for the stark. The vertical error bars represent the 1-$\sigma$ confidence
intervals. The middle column shows the posterior spectra for stark, while the
column on the right-hand side displays the corresponding results for HANSOLO.
Solid lines refer to the median spectra from the posterior sample, while the
shaded areas correspond to the 1-$\sigma$ intervals. We note that the
retrievals for the stark reductions were made for relative occultation depths,
i.e. the mean occultation depths in the middle column are close to zero.
Figure 4: Posterior distributions of the free parameters for the first visit,
representing the CO/CO2-atmosphere scenario. Results are shown for the stark
reduction. We note that $\xi_{\mathrm{CO_{2}}}$ is not a free parameter in the
retrieval but was calculated during a postprocess procedure following the
requirement that in each posterior sample the sum of all $\xi$ values must be
zero. Figure 5: Posterior distributions of the free parameters for the third
visit, representing the SiO/SiO2/MgO-atmosphere scenario. Results are shown
for the stark reduction. We note that $\xi_{\mathrm{MgO}}$ is not a free
parameter in the retrieval but was calculated during a postprocess procedure
following the requirement that in each posterior sample, the sum of all $\xi$
values must be zero.

####  3.2.1 Summary of the retrieval results

The retrieval results for the two different reductions across all five visits
and for the four different model scenarios described in Sect. 2.3 are
summarised in Table 3. The table shows the resulting Bayesian evidence values
$\ln\mathcal{Z}$ and the Bayes factors $B$ with respect to the models with the
highest likelihood value. The former are marked in bold for every visit. Fig.
3 additionally shows the posterior spectra for all models, visits, and
reductions. The detailed posterior distributions for all atmospheric
retrievals can be found in Fig. 4 & 5, as well as in Appendix C.

The results presented in Table 3 suggest that for the HANSOLO reduction, the
planetary blackbody model is always the preferred model. This is likely caused
by the relatively large errors of the reduction that results in the retrieval
favouring a simpler model as can clearly be noticed in the spectra shown in
the right column of Fig. 3.

However, for most visits, the preference for the simple blackbody model is not
statistically significant. The more complex atmospheric scenarios usually have
a Bayes factor of less than three, which suggests that they are essentially
equally likely. For the first three visits, a flat-line fit to the measured
spectrum is effectively ruled out by the Bayesian evidence. The last two
visits, on the other hand, can be fit with any of the four models. There seems
to be little statistical preference for any of the different modelling
scenarios.

The results for the stark reduction show a much broader range of different
models that are statistically preferred. As suggested by Table 3, the first
visit strongly prefers a CO/CO2 atmosphere, the second visit can be explained
by either a CO/CO2 atmosphere, a planetary blackbody, or a siliciate vapour
atmosphere, while the third model overwhelmingly prefers the SiO/SiO2
scenario. The fourth visit is consistent with a planetary blackbody spectrum,
as well as an atmosphere with CO & CO2, or SiO, SiO2 & MgO. Finally, the last
visit strongly prefers a flat-line model.

####  3.2.2 Detailed posterior distributions

Detailed posterior distributions for the preferred model from the stark
reduction of Visit 1 (CO & CO2) and Visit 3 (SiO, SiO2 & MgO), where
atmospheric models are favoured, are shown in Fig. 4 & 5. The posterior
distributions for the first visit reveal a bimodal distribution for the
surface pressure $p_{\mathrm{surf}}$ and the abundances of CO and CO2. As the
two-dimensional correlation plots suggest, the surface pressure has a solution
with a very low value of about $10^{-6.5}$ bar that is dominated by CO in
composition, as well as a higher-pressure mode at about $10^{-3}$ bar that
contains mostly CO2. For comparison, if the outgassing flux were to be
balanced by flux-limited atmospheric escape then the implied surface pressure
is $\sim 10^{-7}$ bar (Heng, 2023). At about 2000 K, the atmosphere
temperature is much warmer than the retrieved temperature for the surface. It
is also important to note, that the posterior distribution for the white-light
occultation depths $d_{\mathrm{wl}}$ is shifted from its prior value of $7\pm
9$ ppm, though they are both still within their 1-$\sigma$ intervals.

The posterior distribution for the SiO/SiO2/MgO model shown in Fig. 5 for the
third visit, on the other hand, exhibits a unimodal pressure distribution with
a median value of about 0.1 bar. Here, the atmosphere is clearly dominated by
SiO2 with only an upper limit for SiO and essentially no constraints on MgO.
The posterior spectra shown in Fig. 3 clearly show the drop-off in the
occulation depth near a wavelength of 4.8 $\mu$m caused by SiO2. Just like in
the previous CO/CO2 scenario for Visit 1, the retrieved atmosphere temperature
is again much higher than the one of the surface.

####  3.2.3 Blackbody temperatures

Figure 6: Retrieved temperatures for all five visits using the blackbody
model. Top panel: results for the stark reduction, bottom panel: HANSOLO.

The resulting posterior distributions of the blackbody temperature models are
shown in Fig. 6 for all visits and the two different reductions. In the case
of the HANSOLO reduction, the blackbody is always the preferred model
according to the Bayesian evidence, though, as previously mentioned, this
preference is statistically not very significant.

As the distributions depicted in the figure suggest, the temperatures
retrieved from the HANSOLO observational data are found in two different
clusters. A low-temperature mode near 750 K is found for Visits 1 and 2 and a
second one at about 1200 K to 1300 K for the other three visits. The
temperatures are quite well constrained with 1-$\sigma$ intervals usually in
the range of about $\pm 100$ K, despite the rather large errors on the
observational data points (see Fig. 3).

For stark, the temperatures are clustered much closer together around a mean
temperature of 1500 K. In comparison to the HANSOLO reduction, however, these
temperatures are less well constrained with 1-$\sigma$ intervals typically
covering a range of several 100 K. This is likely caused by the white-light
occultation depths that are directly correlated with these temperatures.
Following Table 1, they have in general quite large associated errors that
translate into less well-constrained temperatures.

####  3.2.4 Surface pressures

Figure 7: Surface pressure posterior distribution from the CO/CO2 model for
all five visits. Top panel: results for the stark reduction, bottom panel:
HANSOLO.

Figure 8: Surface pressure posterior distribution from the SiO/SiO2/MgO model
for all five visits. Top panel: results for the stark reduction, bottom panel:
HANSOLO.

For the two model scenarios that involve atmospheres, we also retrieve the
surface pressure. For the CO/CO2 model, the corresponding posterior
distributions are shown in Fig. 7, while those for the SiO/SiO2/MgO scenario
are shown in Fig. 8.

In general, the HANSOLO reduction only weakly constrains the surface pressure
with posteriors that usually cover the entire prior range of the pressure from
$10^{-10}$ bar to 500 bar. The posterior distributions seem to be essentially
bimodal for almost every visit, with a very low-pressure mode and a high-
pressure one. These more or less unconstrained pressures are the result of the
rather large errors of the observational data from the HANSOLO reduction.
Those make it difficult to provide good constraints for actual atmospheric
models.

For the stark reduction, the results are more diverse. Some visits seem to
result in very well-constrained surface pressures. This includes Visits 1 and
5 for the CO/CO2 model (see upper panel of Fig. 7) and Visits 1 and 3 for the
SiO/SiO2/MgO case (see upper panel of Fig. 8).

Other visits show the same behaviour as for the HANSOLO reduction: rather
unconstrained surface pressures with usually a bimodal posterior distribution.
Even though not very visible in Fig. 7, the posterior distribution for Visit 1
is also bimodal in shape, with a smaller, high-pressure mode of an atmosphere
dominated by CO2, as discussed above.

We note that our retrieved surface pressures differ from the one reported by
Hu et al. (2024), which corresponds to our Visit 4 and is based on the JWST
program by Hu et al. (2021). However, given that even the two reductions of
the same data in our study produce different results regarding the atmospheric
properties, this is not too surprising. Furthermore, Hu et al. (2024) employed
a different retrieval approach. This includes not using the white-light
eclipse depths of the NIRCam data, imposing a lower limit on the surface
temperature and allowing for a non-radiatively interacting background gas.
Especially the latter assumption will affect the posterior distributions of
the surface pressure.

####  3.2.5 Surface and atmosphere temperatures

Figure 9: Posterior distributions for the atmosphere (top) and surface
temperatures (bottom). The distributions are shown for the CO/CO2 model and
the stark reduction.

For the CO/CO2 model, we present the posteriors for the surface and atmosphere
temperatures in Fig. 9. As discussed in Sect. 2.3, we allow these two
temperatures to have distinct values. We only present the posteriors for the
stark reduction since, as shown above, the HANSOLO one does not provide good
constraints on the atmospheric properties.

Just like the surface pressure, the temperatures are rather well-constrained
for some visits, such as the surface temperatures for Visits 4 and 5.
Observational data from other visits yield much broader distributions, such as
Visit 2, some of which also seem to possess a bimodal shape or only provide
upper limits.

Visit 1 is the only case where the atmosphere seems to have a distinctly
higher temperature than the surface. For other visits, this trend is less
clear. For example, Visit 5 yields a very high surface temperature but the
atmospheric one is less well-constrained and only seems to provide an upper
limit that is roughly equal to the surface temperature. In the case of Visit
3, this situation is reversed. Here, the atmosphere temperature is constrained
with a median value of roughly 1400 K, while the surface temperature only has
an upper limit of about the same value.

##  4 Interpretation of observations

As mentioned in Section 2.1, if the variability in the emission from the
planet is caused by the planet showing different faces during consecutive
occultations, we would expect the occultation depth to be correlated with the
orbit number. However, Fig. 2, which plots the occultation depths as a
function of orbit number, shows that this is not the case. This means that the
observations give no support for a 3:2 spin-orbit resonance being the root
cause for the variability. It is still possible that the planet is trapped in
some higher-order spin-orbit resonance, but to show this by establishing a
pattern would require many more occultation observations than we currently
have.

There are several hypotheses that could potentially explain the full or part
of the observations. We outline two such models in the subsections below: a
transient outgassing atmosphere model and a circumstellar material supported
by the volcanism model. Moreover, the NIRCam data also constrain the presence
of spectral features from a mineral atmosphere resulting from a purported lava
ocean, as described in Section. 4.1 below.

###  4.1 Constraints on silicate atmosphere on 55 Cnc e

Being in proximity to its host, the substellar temperature on 55 Cnc e can
reach ¿ 2000 K. The surface of the planet at such a high temperature is
expected to be molten if there is no atmosphere on the planet. A molten
surface on the planet could then produce a thin rock vapour atmosphere on the
planet. Zilinskas et al. (2022) recently calculated self-consistent models of
outgassed atmospheres for all USPs at the time. They solved the radiative
transfer equations along with equilibrium chemistry models for the outgassed
atmosphere to compute temperature-pressure profile and emission spectra. They
showed that gases such as SiO, SiO2, Na, MgO etc. are some of the main
constituents of these outgassed atmospheres. Their models for 55 Cnc e222All
models are publicly available at https://github.com/zmantas/LavaPlanets are
shown in Fig. 10 overplotted with our observations. The models assume bulk
silicate (oxidised) Earth (BSE) composition for the planet with unevolved and
evolved surface with 80% outgassed efficiency (evolved BSE composition).

Figure 10: Theoretical models of evaporating lava atmospheres for 55 Cnc e
from Zilinskas et al. (2022). Two models are for bulk-silicate composition (in
yellow) and for evolved bulk-silicate composition (in purple). Also
overplotted are photometric occultation depths from 2.1 $\mu$m channel (in
blue) and white-light occultation depths for 4.5 $\mu$m channel (in maroon).
The blue and maroon points are slightly spread in wavelength near their
bandpasses to avoid overlap. We show the two bandpasses corresponding to both
of these channels. The black points show the predicted occultation depths for
both NIRCam bandpasses with different shapes representing bulk-silicate
(stars) and evolved bulk-silicate (squares) compositions.

All of their models with different outgassing efficiencies predict occultation
depths of 70–80 and 145–150 ppm for NIRCam 2.1 and 4.5 $\mu$m channels. As
depicted in Fig. 10 these values are larger compared to our observations. Some
occultation depths are, however, consistent with models at 1–3$\sigma$. One
occultation depth at 2.1 $\mu$m in Visit 5 produces a larger depth compared to
the models. This hints towards a lack of short-wave absorbers such as SiO
and/or SiO2 from the atmosphere that are responsible for thermal inversion
and, in turn, larger occultation depths in NIRCam bandpasses. Indeed, only one
visit (Visit 3) favoured the SiO/SiO2/MgO model in the retrieval analysis. The
band-averaged occultation depth for this visit at 4.5 $\mu$m agrees with the
model prediction (145 ppm for BSE case) at 2.4$\sigma$. However, the short-
wave occultation depth in this visit is inconsistent with the model prediction
at 7$\sigma$. We here note that Hu et al. (2024) found that the occultation
depths in the MIRI bandpass are significantly lower than what is predicted by
Zilinskas et al. (2022) models and thus do not support the presence of the
silicate-rich atmosphere.

At the same time, lower occultation depths in the NIRCam bandpasses could
imply the presence of a gaseous species that have opacity sources in our
NIRCam bandpasses. Alternatively, the lower occultation depths, translated
into lower brightness temperatures, suggest a thick atmosphere with a strong
heat re-distribution (e.g., Hammond & Pierrehumbert, 2017). The estimated
dayside brightness temperatures (see, Table 1) at 4.5 $\mu$m (Table 1) in all
visits are smaller than the expected dayside temperature333Computed using
$T_{\textnormal{day}}=T_{\star}\sqrt{\frac{R_{\star}}{a}}(1-A_{B})^{1/4}f^{1/4}$,
while using zero bond albedo and the heat re-distribution factor $f=2/3$, for
a bare rock with no heat re-distribution (Burrows, 2014; Koll et al., 2019).
of 2537 K indicating the presence of heat transfer. In either case, our
observations seem to indicate the existence of volatiles in the atmosphere of
55 Cnc e. However, it is still challenging to explain the very large
occultation depth (and, thus, hot brightness temperature — 3138 K; see, Table
1) observed at 2.1 $\mu$m in Visit 5.

###  4.2 Constraints on an outgassed secondary atmosphere

Heng (2023) previously suggested that a transient, outgassed secondary
atmosphere is capable of simultaneously explaining the observed variability of
55 Cnc e in both the optical/visible and infrared range of wavelengths.
Specifically, atmospheres of several tens of bars of pure carbon monoxide (CO)
are capable of producing occultation depths of about 21 ppm in the CHEOPS and
TESS bandpasses, which are consistent with most of the occultation depths
measured by CHEOPS (Meier Valdés et al., 2023) and TESS (Meier Valdés et al.,
2022). However, a change in atmospheric surface pressure of several tens of
bars through loss processes or outgassing over the observed variability time
scale in the CHEOPS data is difficult to explain. Such outgassed atmospheres
are incapable of producing occultation depths as high as $\approx 40$–50 ppm,
which were measured thrice in Fig. 3 of Meier Valdés et al. (2023). Similarly,
they cannot produce phase variations as high as 110 ppm as measured by MOST
(Sulis et al., 2019). It cannot be ruled out that these anomalously high
occultation depths are associated with stellar activity.

For the first data reduction (stark), the outgassed atmosphere with CO and CO2
is associated with the highest Bayesian evidence in Visits 1 & 2\. Bayesian
model comparison does not disfavour this interpretation of Visit 4 as well.
Fig. 4 shows the interpretation of the spectrum from Visit 1 using a CO+CO2
atmosphere. For Visit 3, a silicate-vapour atmosphere is strongly preferred
over an outgassed atmosphere (with the logarithm of the Bayes factor being
9.8; Fig. 5). For the more conservative second data reduction (HANSOLO), the
retrieval associated with the highest Bayesian evidence is a blackbody curve
over all 5 visits.

The simplest interpretation of the spectra is using a blackbody curve, which
is consistent with the data in Visits 2 and 4 of the stark reduction and all 5
visits of the HANSOLO reduction. Fig. 6 shows the posterior distributions of
the blackbody temperature. For Visits 2 and 4 of the stark reduction, the
blackbody temperature is broadly between 1500 K and 2000 K. Note that a
blackbody curve does not automatically imply that one is probing a bare rocky
surface, since an optically thick, isothermal atmosphere may also produce a
blackbody curve (Heng, 2023). For the HANSOLO reductions, the blackbody
temperature is about 750 K for Visits 1 and 2 and increases to about 1250 K
for Visits 3, 4 and 5 over a period of about 2.2 days (between Visits 2 and
3). Such a duration is not inconsistent with the radiative timescale, which is
under an Earth day for $\sim 1$ bar atmospheres (Heng, 2023). If 55 Cnc e has
a bare rocky surface and negligible albedo, then its temperature would be the
equilibrium temperature of about 2000 K. If we take these blackbody
temperatures (750 K and 1250 K) seriously, then it implies that the
observations are not probing a bare rocky surface that has reached a steady
state with the stellar instellation, unless one assumes implausibly high
surface albedos.

If we focus on the interpretation of the spectra using CO-CO2 atmospheres,
then Figs. 7 and 9 show the posterior distributions of surface pressures,
atmospheric temperatures and surface pressures. For the HANSOLO data
reductions, the surface pressure is unconstrained. For Visits 1, 2 and 4 of
the stark reduction, the inferred surface pressure is $\sim 1$ $\mu$bar. The
surface temperature is $\sim 1000$ K, which is only possible if the surface
has not come to radiative equilibrium with the stellar instellation because of
the presence of an atmosphere. The atmospheric temperature jumps from $\sim
2000$ K to $\sim 2500$ K to $\sim 1500$ K from Visits 1 to 2 to 3. While this
is not implausible because of the short radiative timescales, we do not have a
mechanism to explain how and why this happens.

###  4.3 Can a circumstellar inhomogeneous dusty torus explain variability?

Two of our observations, Visit 1 at 4.5 $\mu$m and Visit 2 at 2.1 $\mu$m, show
occultation depths that are consistent with zero at 1-$\sigma$. These non-
detections are challenging to explain with any kind of atmospheric phenomena.
Moreover, the occultation depths observed at 2.1 $\mu$m and 4.5 $\mu$m are not
correlated with each other (Fig. 2), which potentially hints towards different
origins of variability in different wavelength channels.

A grey absorber could explain the optical and 2.1 $\mu$m channel variability.
A natural candidate for this grey absorber is a circumstellar dust torus
(Sulis et al., 2019; Meier Valdés et al., 2023). The progenitor of the dusty
torus could be the volcanism on 55 Cnc e developed by the extreme tidal
heating akin to Io (e.g., Oza et al., 2019b; Gebek & Oza, 2020). The most
common gases from volcanism seen on the Earth, Io, and Venus - e.g., SO2, CO2,
generate a tenuous atmosphere on the planet. Volcanism, supported by
significant tidal heating, is expected to expel a prodigious quantity of dust
grains into the upper atmosphere, which ultimately escape the planet’s
gravitational sphere of influence due to impinging stellar ions. Upon escape,
such a mechanism may eventually generate a patchy, circumstellar dust torus,
which has been shown to be sufficiently opaque in visible light to produce
optical variability (Meier Valdés et al., 2023). Volcanic gases are additional
non-trivial sources of opacity in our NIRCam 4.5 $\mu$m channel. Analytical
models showed that an optically thin (e.g., Gebek & Oza, 2020) SO2 atmosphere
with a range of pressures can produce the IR variability observed with
Spitzer. Since the Spitzer/IRAC bandpass at 4.5 $\mu$m and our NIRCam/F444W
bandpass have a large overlap in wavelength, it remains a possibility that a
similar thin SO2 (or, any other volcanic gases such as CO2 which also absorbs
at 4.5 $\mu$m) atmosphere with several tens of $\mu$bar could explain the
observed variability in our NIRCam dataset. To evaluate this idea in detail is
however beyond the scope of the present work and instead planned for an
upcoming publication (Oza et al., in prep.).

The variability at 2.1 $\mu$m is difficult to explain with a thin atmosphere
consisting volcanic gases such as SO2 or CO2 since they do not have
significant opacity in the 2.1 $\mu$m bandpass. Instead, the dust grains
present in the torus could be a cause of this variability, which was also
hypothesised by Meier Valdés et al. (2023). If the grain size is larger than
0.3 $\mu$m from the size range of 0.1–0.7 $\mu$m discussed in Morris et al.
(2021) and Meier Valdés et al. (2023), the particles will be opaque in the 2.1
$\mu$m channel, but transparent in the 4.5 $\mu$m channel. Although many
Earth-like dust species do not survive long enough in the circumstellar
environment, dust made of quartz, silicon carbide and graphite can survive a
significant fraction of an orbit to generate a patchy torus (Meier Valdés et
al., 2023). Following the same formalism from Meier Valdés et al. (2023), the
mass loss needed to account for the maximum change in occultation depth (95.9
ppm, in visit 5) 2.5–5.7 $\times$ 106 kg s-1 is within a factor of two of the
maximum escape rate derived by CHEOPS, reported to be as large as $\sim$ 2.9
$\times$ 106 kg s-1 (Meier Valdés et al., 2023). If the particle size is
larger than 0.7 $\mu$m, they can, in principle, even explain the variability
at 4.5 $\mu$m channel. However, the non-correlation of occultation depths at
2.1 $\mu$m and 4.5 $\mu$m channels suggests that although the two sources may
be linked, they are indeed distinct absorbers, e.g, grains and gas at 2.1 and
4.5 $\mu$m, respectively, as mentioned above. However, the effect of the dust
torus on the transit observations is yet to be found observationally. In
particular, if the dust escape happens during a transit event, dust could
float in the Hill sphere of the planet or form a comet-like tail (e.g., Brogi
et al., 2012). Both processes should affect the transit light curve in the
form of a significantly large transit depth and an asymmetric transit shape,
respectively, unless dust very quickly leaves the vicinity of the planet.

It is unknown what escape mechanism is currently operating at 55 Cnc e, and
therefore more phase curve observations, especially at shorter wavelengths
where Si in the dust have emission lines, are needed to monitor the
variability. Multiple phase curves would scan the whole circumstellar region
over time to determine the location of the dusty torus and how it evolves,
helping in a better understanding of the escape mechanism and thus
variability. However, based on its close proximity several mechanisms
including canonical photoevaporation and boil-off (Affolter et al., 2023) are
able to reproduce the estimated escape rate. For close-in rocky bodies like 55
Cnc e, more energetic plasma escape mechanisms including ion-neutral
interactions such as atmospheric sputtering (Oza et al., 2019b; Meyer zu
Westram et al., 2024), which, similar to Io, drive a feedback process sourced
by the melting and degassing of the rocky body itself via induction-heating
(Lanza, 2021) and two body tidal-heating (Oza et al., 2019a; Quick et al.,
2020; Charnoz et al., 2021).

The aforementioned escape mechanisms are source-limited by geological activity
and expected to vary on orbital timescales in phase-curve observations (Meyer
zu Westram et al., 2024). Source-limited implies that the escape rate is
ultimately limited by the outgassing rate below the escape layer, such that if
the supply rate were zero, escape would not occur. Effectively, the discussed
energetic escape mechanisms naturally generate extended neutral and grain
clouds that provide a toroidal opacity source in the circumstellar
environment.

###  4.4 Can stellar activity cause the occultation depth variability?

Stellar activity can, in principle, cause the occultation depth variability of
55 Cnc e. Demory et al. (2023) checked whether stellar granulation could
explain the optical occultation depth variability found with CHEOPS. They,
however, rejected stellar activity as a source of variability due to very low
occultation depths in some visits and their detection of a sinusoidal temporal
trend of the variability. Furthermore, the photometric monitoring of the star
for about 11 years in the optical from the ground revealed a photometric
variability of 0.006 mag which is too small to explain the $\sim$ 50 ppm
occultation depth variability observed with CHEOPS (Fischer et al., 2008;
Demory et al., 2023). The stellar activity signal is expected to decrease at
longer wavelengths. This means that it is challenging to explain IR
variability with the photometric variation of mmag level observed by Fischer
et al. (2008) in the optical. Moreover, the activity has to happen every
instance during the short time window around the occultation, which is
improbable. In any case, the inflation of uncertainties with the injection-
retrieval method accounts for any noise, including the correlated noise. The
fact that the maximum difference in the occultation depths is significant even
with inflated uncertainties suggests that the origin of the occultation depth
variability is not related to the star.

##  5 Conclusions

We obtained time on JWST/NIRCam to study the dayside emission variability of
55 Cnc e (GO 2084: PI Brandeker and GO 1952: PI Hu). In particular, we test
the hypothesis that 55 Cnc e is in a 3:2 spin-orbit resonance, thus showing
different faces at every occultation and thereby explaining the observed
dayside variability and also the hot-spot displacement from the sub-stellar
location. The prediction was that this would result in occultation depths
highly correlated with their orbital number parity, at least over short time
scales.

We observed five occultations of 55 Cnc e in two wavelength bands, or
channels, a spectroscopic band at 4.5 $\mu$m and a single photometric band at
2.1 $\mu$m. Four of them are observed within a week, i.e., in the duration of
eight planetary orbits, while the last was observed after five months. We
analysed the data using six different pipelines. Our main finding is that the
occultation depths change strongly, from a non-detection to 100 ppm, and
rapidly (within a week). The variability is however not observed to correlate
with the occultation number parity, implying that a planet 3:2 spin-orbit
resonance is not the reason for its variability. The variability is observed
in both 2.1 and 4.5 $\mu$m channels, but is curiously not correlated between
channels. The estimated brightness temperature at 4.5 $\mu$m varies between
873 K – 2256 K. These values are less than the predicted dayside temperature
in case of zero heat re-distribution and zero albedo, 2537 K, which hints at
the presence of a planetary atmosphere enabling the heat re-distribution.

The spectroscopic data at 4.5 $\mu$m is affected by correlated noise of
unknown origin. Although the results from different reductions overall agree
well with each other, there are several differences in white-light occultation
depths and emission spectra that can be attributed to different treatments of
correlated noise. We select two representative reductions, stark and HANSOLO,
to perform atmospheric retrieval. Our atmospheric retrieval was performed
using two simple atmospheric models containing an isothermal atmosphere made
up of either CO/CO2 or SiO/SiO2/MgO. Additionally, we also tested a blackbody
model and a flat line model with no atmospheric features. Retrievals performed
with HANSOLO results mainly favour a blackbody model owing to larger errorbars
on the occultation depths. However, other models with CO/CO2 or SiO/SiO2/MgO
were not discarded either, statistically. The retrievals with stark prefer
CO/CO2 atmospheres in at least two visits, SiO/SiO2/MgO atmosphere in one
visit and blackbody and flat line models in the remaining two visits. The
CO/CO2 atmosphere could be generated from outgassing of the surface (e.g.,
Heng, 2023). The outgassing could be stochastic and thus can potentially
explain the variability. As already advocated by Heng (2023), simultaneous
observations in the optical and infrared are needed to corroborate (or refute)
the presence of a transient outgassed CO/CO2 atmosphere.

The occultation depth variability in the 2.1 $\mu$m channel, especially its
uncorrelated behaviour with its 4.5 $\mu$m channel counterpart, is challenging
to explain with a simple atmospheric model. It is possible that the
variability seen at 2.1 $\mu$m and that at 4.5 $\mu$m have different origins.
A circumstellar inhomogeneous cloud of dust could potentially describe the
variability at 2.1 $\mu$m. Volcanism induced by extreme tidal heating of 55
Cnc e could be a natural source of dust in the atmosphere of the planet which
would eventually escape the planet and generate a patchy dusty torus in the
circumstellar environment. The presence of dust in the circumstellar
environment could also be helpful in the interpretation of several non-
detection of occultation depths found in our observations as it could hide our
view of the planet. More observations at shorter wavelengths, e.g., in
ultraviolet, would help to more strongly constrain the presence of a
circumstellar patchy dust torus. Simultaneous observations in near and mid-IR
around 4 and 8 $\mu$m where volcanic gases CO2/SO2 have opacity would be
helpful in constraining their presence. Such multiple observations in the
optical and IR would not only constrain the presence of a circumstellar dust
torus and atmosphere on the planet but also probe how these components evolve
with time, essentially distinguishing both scenarios discussed in this work.

While we do find a hint of an atmosphere on the planet in at least some
visits, corroborating Hu et al. (2024), the simple picture of a static
atmosphere cannot explain all observational features. A more complex model,
including an outgassed atmosphere, circumstellar material, and perhaps
dynamical processes in the atmosphere, would probably be needed to explain the
entire range of observations. Moreover, given the strong variability of the
system, simultaneous multi-wavelength observations would go a long way to
distinguish between possible explanations and help probe the true nature of 55
Cnc e.

###### Acknowledgements.

We would like to thank an anonymous referee for their detailed referee report
and suggestions which significantly improved the manuscript. JAP acknowledges
Néstor Espinoza for discussing the peculiarities of JWST data analysis. JAP
would like to thank Ludmila Carone for an insightful dialogue about
theoretical models of the planet. JAP and ABr were supported by the Swedish
National Space Agency (SNSA). The contributions of DP and ML have been carried
out within the framework of the NCCR PlanetS supported by the Swiss National
Science Foundation under grants 51NF40$\\_$182901 and 51NF40$\\_$205606\. DP
and ML also acknowledge support of the Swiss National Science Foundation under
grant number PCEFP2_194576. EMV acknowledges support from the Centre for Space
and Habitability (CSH). This work has been carried out within the framework of
the National Centre of Competence in Research PlanetS supported by the Swiss
National Science Foundation under grant 51NF40_205606. EMV acknowledges the
financial support of the SNSF. This project has received funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (project Spice Dune, grant agreement No
947634, and Four Aces; grant agreement No 724427). ADe and DEh have received
funding from the Swiss National Science Foundation for project 200021_200726.
This work has also been carried out within the framework of the National
Centre of Competence in Research PlanetS supported by the Swiss National
Science Foundation under grant 51NF40_205606. This research has made use of
the Spanish Virtual Observatory (https://svo.cab.inta-csic.es) project funded
by MCIN/AEI/10.13039/501100011033/ through grant PID2020-112949GB-I00. CMP and
MF gratefully acknowledge the support of the SNSA (DNR 65/19, 177/19). BOD
acknowledges support from the Swiss State Secretariat for Education, Research
and Innovation (SERI) under contract number MB22.00046. Part of this research
was carried out at the Jet Propulsion Laboratory, California Institute of
Technology, under a contract with the National Aeronautics and Space
Administration (80NM0018D0004). Part of the High Performance Computing
resources used in this investigation were provided by funding from the JPL
Information and Technology Solutions Directorate. Finally, we thank ERASMUS
student Charlotte Zimmermann for her contributions to the initial studies of
this work.

## References

  * Affolter et al. (2023) Affolter, L., Mordasini, C., Oza, A. V., Kubyshkina, D., & Fossati, L. 2023, A&A, 676, A119
  * Ahrer et al. (2023) Ahrer, E.-M., Stevenson, K. B., Mansfield, M., et al. 2023, Nature, 614, 653
  * Angelo & Hu (2017) Angelo, I. & Hu, R. 2017, AJ, 154, 232
  * Bell et al. (2022) Bell, T., Ahrer, E.-M., Brande, J., et al. 2022, The Journal of Open Source Software, 7, 4503
  * Benneke & Seager (2012) Benneke, B. & Seager, S. 2012, ApJ, 753, 100
  * Bourrier et al. (2018a) Bourrier, V., Dumusque, X., Dorn, C., et al. 2018a, A&A, 619, A1
  * Bourrier et al. (2018b) Bourrier, V., Ehrenreich, D., Lecavelier des Etangs, A., et al. 2018b, A&A, 615, A117
  * Bradley et al. (2023) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2023, astropy/photutils: 1.8.0
  * Brandeker (2019) Brandeker, A. 2019, in AAS/Division for Extreme Solar Systems Abstracts, Vol. 51, AAS/Division for Extreme Solar Systems Abstracts, 311.07
  * Brandeker et al. (2021) Brandeker, A., Alibert, Y., Bourrier, V., et al. 2021, Is it raining lava in the evening on 55 Cancri e?, JWST Proposal. Cycle 1, ID. #2084
  * Brogi et al. (2012) Brogi, M., Keller, C. U., de Juan Ovelar, M., et al. 2012, A&A, 545, L5
  * Bruntt et al. (2008) Bruntt, H., De Cat, P., & Aerts, C. 2008, A&A, 478, 487
  * Burrows (2014) Burrows, A. S. 2014, Proceedings of the National Academy of Science, 111, 12601
  * Bushouse et al. (2023) Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2023, JWST Calibration Pipeline
  * Castelli & Kurucz (2004) Castelli, F. & Kurucz, R. L. 2004, astro-ph/0405087 [astro-ph/0405087]
  * Charnoz et al. (2021) Charnoz, S., Sossi, P. A., Lee, Y.-N., et al. 2021, Icarus, 364, 114451
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102
  * Crida et al. (2018) Crida, A., Ligi, R., Dorn, C., & Lebreton, Y. 2018, ApJ, 860, 122
  * Crossfield (2012) Crossfield, I. J. M. 2012, A&A, 545, A97
  * Dawson & Fabrycky (2010) Dawson, R. I. & Fabrycky, D. C. 2010, ApJ, 722, 937
  * Deibert et al. (2021) Deibert, E. K., de Mooij, E. J. W., Jayawardhana, R., et al. 2021, AJ, 161, 209
  * Demory et al. (2016a) Demory, B.-O., Gillon, M., de Wit, J., et al. 2016a, Nature, 532, 207
  * Demory et al. (2011) Demory, B. O., Gillon, M., Deming, D., et al. 2011, A&A, 533, A114
  * Demory et al. (2016b) Demory, B.-O., Gillon, M., Madhusudhan, N., & Queloz, D. 2016b, MNRAS, 455, 2018
  * Demory et al. (2012) Demory, B.-O., Gillon, M., Seager, S., et al. 2012, ApJ, 751, L28
  * Demory et al. (2023) Demory, B. O., Sulis, S., Meier Valdés, E., et al. 2023, A&A, 669, A64
  * Dorn et al. (2017) Dorn, C., Hinkel, N. R., & Venturini, J. 2017, A&A, 597, A38
  * Doyle et al. (2014) Doyle, A. P., Davies, G. R., Smalley, B., Chaplin, W. J., & Elsworth, Y. 2014, MNRAS, 444, 3592
  * Dragomir et al. (2014) Dragomir, D., Matthews, J. M., Winn, J. N., & Rowe, J. F. 2014, in Formation, Detection, and Characterization of Extrasolar Habitable Planets, ed. N. Haghighipour, Vol. 293, 52–57
  * Ehrenreich et al. (2012) Ehrenreich, D., Bourrier, V., Bonfils, X., et al. 2012, A&A, 547, A18
  * Espinoza (2022) Espinoza, N. 2022, TransitSpectroscopy
  * Espinoza et al. (2019) Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262
  * Espinoza et al. (2023) Espinoza, N., Úbeda, L., Birkmann, S. M., et al. 2023, PASP, 135, 018002
  * Esteves et al. (2017) Esteves, L. J., de Mooij, E. J. W., Jayawardhana, R., Watson, C., & de Kok, R. 2017, AJ, 153, 268
  * Feroz & Hobson (2008) Feroz, F. & Hobson, M. P. 2008, MNRAS, 384, 449
  * Fischer et al. (2008) Fischer, D. A., Marcy, G. W., Butler, R. P., et al. 2008, ApJ, 675, 790
  * Folsom et al. (2020) Folsom, C. P., Ó Fionnagáin, D., Fossati, L., et al. 2020, A&A, 633, A48
  * Gebek & Oza (2020) Gebek, A. & Oza, A. V. 2020, MNRAS, 497, 5271
  * Gillon et al. (2012) Gillon, M., Demory, B. O., Benneke, B., et al. 2012, A&A, 539, A28
  * Grimm & Heng (2015) Grimm, S. L. & Heng, K. 2015, ApJ, 808, 182
  * Grimm et al. (2021) Grimm, S. L., Malik, M., Kitzmann, D., et al. 2021, ApJS, 253, 30
  * Hammond & Pierrehumbert (2017) Hammond, M. & Pierrehumbert, R. T. 2017, ApJ, 849, 152
  * Heng (2023) Heng, K. 2023, ApJ, 956, L20
  * Hu et al. (2024) Hu, R., Bello-Arufe, A., Zhang, M., et al. 2024, arXiv e-prints, arXiv:2405.04744
  * Hu et al. (2021) Hu, R., Brandeker, A., Damiano, M., et al. 2021, Determining the Atmospheric Composition of the Super-Earth 55 Cancri e, JWST Proposal. Cycle 1, ID. #1952
  * Husser et al. (2013) Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6
  * Jindal et al. (2020) Jindal, A., de Mooij, E. J. W., Jayawardhana, R., et al. 2020, AJ, 160, 101
  * Keles et al. (2022) Keles, E., Mallonn, M., Kitzmann, D., et al. 2022, MNRAS, 513, 1544
  * Kirk et al. (2024) Kirk, J., Stevenson, K. B., Fu, G., et al. 2024, AJ, 167, 90
  * Kite et al. (2016) Kite, E. S., Fegley, Bruce, J., Schaefer, L., & Gaidos, E. 2016, ApJ, 828, 80
  * Kitzmann et al. (2020) Kitzmann, D., Heng, K., Oreshenko, M., et al. 2020, ApJ, 890, 174
  * Koll et al. (2019) Koll, D. D. B., Malik, M., Mansfield, M., et al. 2019, ApJ, 886, 140
  * Kreidberg (2015) Kreidberg, L. 2015, PASP, 127, 1161
  * Kreidberg et al. (2019) Kreidberg, L., Koll, D. D. B., Morley, C., et al. 2019, Nature, 573, 87
  * Kurucz (1993) Kurucz, R. L. 1993, VizieR Online Data Catalog, VI/39
  * Kurucz (2013) Kurucz, R. L. 2013, ATLAS12: Opacity sampling model atmosphere program, Astrophysics Source Code Library
  * Lanza (2021) Lanza, A. F. 2021, A&A, 653, A112
  * Lendl et al. (2017) Lendl, M., Cubillos, P. E., Hagelberg, J., et al. 2017, A&A, 606, A18
  * Lendl et al. (2016) Lendl, M., Delrez, L., Gillon, M., et al. 2016, A&A, 587, A67
  * Li et al. (2015) Li, G., Gordon, I. E., Rothman, L. S., et al. 2015, ApJS, 216, 15
  * Li et al. (2019) Li, H. Y., Tennyson, J., & Yurchenko, S. N. 2019, MNRAS, 486, 2351
  * Lopez (2017) Lopez, E. D. 2017, MNRAS, 472, 245
  * Madhusudhan et al. (2012) Madhusudhan, N., Lee, K. K. M., & Mousis, O. 2012, ApJ, 759, L40
  * Mahapatra et al. (2017) Mahapatra, G., Helling, C., & Miguel, Y. 2017, MNRAS, 472, 447
  * Marsh (1989) Marsh, T. R. 1989, PASP, 101, 1032
  * Mayor et al. (2003) Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20
  * McArthur et al. (2004) McArthur, B. E., Endl, M., Cochran, W. D., et al. 2004, ApJ, 614, L81
  * Meier et al. (2023) Meier, T. G., Bower, D. J., Lichtenberg, T., Hammond, M., & Tackley, P. J. 2023, A&A, 678, A29
  * Meier Valdés et al. (2023) Meier Valdés, E. A., Morris, B. M., Demory, B. O., et al. 2023, A&A, 677, A112
  * Meier Valdés et al. (2022) Meier Valdés, E. A., Morris, B. M., Wells, R. D., Schanche, N., & Demory, B. O. 2022, A&A, 663, A95
  * Mercier et al. (2022) Mercier, S. J., Dang, L., Gass, A., Cowan, N. B., & Bell, T. J. 2022, AJ, 164, 204
  * Meyer zu Westram et al. (2024) Meyer zu Westram, M., Oza, A. V., & Galli, A. 2024, Journal of Geophysical Research (Planets), 129, e2023JE007935
  * Morris et al. (2021) Morris, B. M., Delrez, L., Brandeker, A., et al. 2021, A&A, 653, A173
  * of Ockham (1495) of Ockham, W. 1495, Quaestiones et decisiones in quattuor libros Sententiarum Petri Lombardi: Centilogium theologicum (Johannes Trechsel)
  * Owens et al. (2020) Owens, A., Conway, E. K., Tennyson, J., & Yurchenko, S. N. 2020, MNRAS, 495, 1927
  * Oza et al. (2019a) Oza, A., Bower, D. J., Demory, B.-O., et al. 2019a, in EPSC-DPS Joint Meeting 2019, Vol. 2019, EPSC–DPS2019–1714
  * Oza et al. (2019b) Oza, A. V., Johnson, R. E., Lellouch, E., et al. 2019b, ApJ, 885, 168
  * Persson et al. (2018) Persson, C. M., Fridlund, M., Barragán, O., et al. 2018, A&A, 618, A33
  * Piskunov & Valenti (2017) Piskunov, N. & Valenti, J. A. 2017, A&A, 597, A16
  * Quick et al. (2020) Quick, L. C., Roberge, A., Mlinar, A. B., & Hedman, M. M. 2020, PASP, 132, 084402
  * Rasmussen et al. (2023) Rasmussen, K. C., Currie, M. H., Hagee, C., et al. 2023, AJ, 166, 155
  * Ridden-Harper et al. (2016) Ridden-Harper, A. R., Snellen, I. A. G., Keller, C. U., et al. 2016, A&A, 593, A129
  * Rodríguez et al. (2012) Rodríguez, A., Callegari, N., Michtchenko, T. A., & Hussmann, H. 2012, MNRAS, 427, 2239
  * Salz et al. (2016) Salz, M., Czesla, S., Schneider, P. C., & Schmitt, J. H. M. M. 2016, A&A, 586, A75
  * Schlawin et al. (2023) Schlawin, E., Beatty, T., Brooks, B., et al. 2023, PASP, 135, 018001
  * Skilling (2004) Skilling, J. 2004, in 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, American Institute of Physics Conference Series, Vol. 735, American Institute of Physics Conference Series, ed. R. Fischer, R. Preuss, & U. V. Toussaint, 395–405
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132
  * Sulis et al. (2019) Sulis, S., Dragomir, D., Lendl, M., et al. 2019, A&A, 631, A129
  * Tabernero et al. (2020) Tabernero, H. M., Allende Prieto, C., Zapatero Osorio, M. R., et al. 2020, MNRAS, 498, 4222
  * Tamburo et al. (2018) Tamburo, P., Mandell, A., Deming, D., & Garhart, E. 2018, AJ, 155, 221
  * Tian & Heng (2024) Tian, M. & Heng, K. 2024, ApJ, 963, 157
  * Trotta (2008) Trotta, R. 2008, Contemporary Physics, 49, 71
  * Trotta (2017) Trotta, R. 2017, arXiv e-prints, arXiv:1701.01467
  * Tsiaras et al. (2016) Tsiaras, A., Rocchetto, M., Waldmann, I. P., et al. 2016, ApJ, 820, 99
  * Valenti & Piskunov (1996) Valenti, J. A. & Piskunov, N. 1996, A&AS, 118, 595
  * van Dokkum (2001) van Dokkum, P. G. 2001, PASP, 113, 1420
  * van Leeuwen (2007) van Leeuwen, F. 2007, A&A, 474, 653
  * Vines & Jenkins (2022) Vines, J. I. & Jenkins, J. S. 2022, arXiv e-prints, arXiv:2204.03769
  * von Braun et al. (2011) von Braun, K., Boyajian, T. S., ten Brummelaar, T. A., et al. 2011, ApJ, 740, 49
  * Winn et al. (2011) Winn, J. N., Matthews, J. M., Dawson, R. I., et al. 2011, ApJ, 737, L18
  * Winn et al. (2018) Winn, J. N., Sanchis-Ojeda, R., & Rappaport, S. 2018, New A Rev., 83, 37
  * Yurchenko et al. (2020) Yurchenko, S. N., Mellor, T. M., Freedman, R. S., & Tennyson, J. 2020, MNRAS, 496, 5282
  * Yurchenko et al. (2022) Yurchenko, S. N., Tennyson, J., Syme, A.-M., et al. 2022, MNRAS, 510, 903
  * Zhang et al. (2021) Zhang, M., Knutson, H. A., Wang, L., et al. 2021, AJ, 161, 181
  * Zieba et al. (2022) Zieba, S., Zilinskas, M., Kreidberg, L., et al. 2022, A&A, 664, A79
  * Zilinskas et al. (2022) Zilinskas, M., van Buchem, C. P. A., Miguel, Y., et al. 2022, A&A, 661, A126

##  Appendix A Data analysis methods

This section details six independent methods of analysing the JWST/NIRCam
data. In Table 4, we summarise the white-light occultation depths between
about 4 and 5 $\mu$m (see, below for exact wavelength range for different
methods) and photometric occultation depths at 2.1 $\mu$m. Figure 11 compares
the relative occultation depth spectra for all visits from different methods.
It can be seen from Figure 11 and Table 4 that the results obtained with
various independent analysis methods overall agree with each other, however,
there are some differences which could be attributed to the different handling
of correlated noise in the data. For example, HANSOLO reduction uses Gaussian
processes (GP) to model the correlated noise and thus produces results, white-
light and spectroscopic occultation depths, that are the most distinct from
the rest of the methods. On the other hand, reduction methods from, e.g.,
stark, inflate errorbars on occultation depths to account for correlated
noise. We use results from HANSOLO and stark as two representative methods in
our atmospheric retrieval analysis and interpretation. We describe each
analysis method below.

Table 4: Comparison of white-light and photometric occultation depths from
different methods

Visit | stark | Eureka! | Eureka! | HANSOLO | transitspectroscopy | SPARTA
---|---|---|---|---|---|---
| (ppm) | R1 (ppm) | R2 (ppm) | (ppm) | (ppm) | (ppm)
White-light occultation depths |  |  |  |  |  |
Visit 1 (Nov 18, 2022) | $7.0^{+8.8}_{-8.8}$ | $-$ | $49.2^{+12.4}_{-12.3}$ | $2.6^{+14.1}_{-2.6}$ | $15.9^{+11.6}_{-11.4}$ | $52.1^{+11.1}_{-10.3}$
Visit 2 (Nov 20, 2022) | $65.2^{+22.3}_{-42.2}$ | $-$ | $85.1^{+9.6}_{-9.8}$ | $6.4^{+31.0}_{-6.1}$ | $52.2^{+11.3}_{-11.5}$ | $79.0^{+10.0}_{-9.5}$
Visit 3 (Nov 23, 2022) | $101.4^{+17.1}_{-32.4}$ | $-$ | $130.9^{+10.3}_{-11.3}$ | $112.1^{+28.4}_{-31.9}$ | $141.9^{+11.5}_{-12.0}$ | $119.1^{+10.8}_{-10.3}$
Visit 4 (Nov 24, 2022) | $119.2^{+34.0}_{-19.0}$ | $-$ | $134.1^{+9.6}_{-9.5}$ | $37.8^{+28.8}_{-24.1}$ | $115.5^{+8.9}_{-8.9}$ | $82.9^{+18.0}_{-18.3}$
Visit 5 (Apr 24, 2023) | $95.4^{+13.5}_{-16.8}$ | $-$ | $106.7^{+9.2}_{-11.7}$ | $73.5^{+21.3}_{-21.4}$ | $98.6^{+11.0}_{-10.8}$ | $95.9^{+11.3}_{-10.1}$
Photometric occultation depths |  |  |  |  |  |
Visit 1 (Nov 18, 2022) | $47.4^{+21.0}_{-15.5}$ | $42.8^{+4.9}_{-4.7}$ | $-$ | $-$ | $-$ | $-$
Visit 2 (Nov 20, 2022) | $-5.1^{+5.5}_{-6.0}$ | $-9.8^{+5.6}_{-6.0}$ | $-$ | $-$ | $-$ | $-$
Visit 3 (Nov 23, 2022) | $37.3^{+4.7}_{-4.6}$ | $28.2^{+5.5}_{-5.6}$ | $-$ | $-$ | $-$ | $-$
Visit 4 (Nov 24, 2022) | $36.8^{+27.7}_{-32.9}$ | $39.5^{+6.0}_{-5.6}$ | $-$ | $-$ | $-$ | $-$
Visit 5 (Apr 24, 2023) | $95.9^{+8.1}_{-7.9}$ | $92.4^{+5.9}_{-5.5}$ | $-$ | $-$ | $-$ | $-$

444The uncertainties are 68 percentile of the corresponding posterior
distribution. Visit 4 is the archival observation from Hu et al. (2024).

Figure 11: Comparison of occultation depth spectra for all observations from
different methods: (Left) Relative occultation depth spectra from stark
(baseline spectra, in orange), Eureka! (in blue) and transitspectroscopy (in
green), and absolute occultation depth spectra minus white-light depth for
SPARTA (in purple). (Right) stark relative occultation depth spectra (in
orange) and HANSOLO absolute occultation depth spectra minus white-light depth
(in grey).

###  A.1 stark

As described in Section 2.2, the observations were carried out using NIRCam
grism timeseries observing mode, which has two channels, a long-wave (LW, 4.5
$\mu$m channel) spectroscopic channel and a short-wave (SW, 2.1 $\mu$m
channel) photometric channel. We analysed both datasets with our pipeline.

Figure 12: Raw photometric light curves from the short-wave channel at 2.1
$\mu$m (in blue) and raw white-light light curves from the long-wave channel
at 4.5 $\mu$m (in orange) for Visit 1 to 3 and 5 (GO 2084, in the top panel)
and for Visit 4 (GO 1952, bottom panel). A darker and lighter shade of colours
depicts the even and odd parity of the observations. The darker points on the
top of the main data show the binned data points.

Figure 13: Allan deviation plots of residuals from photometric light curve
analysis from 2.1 $\mu$m (SW) channel (left panel, in blue) and 4.5 $\mu$m
(LW) channel white-light light curve analysis (right panel, in orange). Figure
14: Posteriors of occultation depths from injection-retrieval exercise (see,
text) for 2.1 $\mu$m (SW) channel (the top row, in blue) and 4.5 $\mu$m (LW)
channel (the bottom row, in orange). The dashed and dotted vertical lines are
injected and retrieved – a median of the posteriors – values of occultation
depths, respectively. The median and 68-percentile confidence intervals of the
posterior are written on the top of the plots.

####  A.1.1 Long-wave data analysis

We downloaded uncalibrated data files (uncal files) from the MAST archive and
used the official jwst pipeline to produce calibrated from them. We ran Stage
1 of the jwst pipeline with some modifications to the uncal files. The main
change in Stage 1 is that we skipped the dark current step and jump step. This
is justified because the dark current level in NIRCam detectors is low.
Furthermore, since our observations were carried out using only two groups per
integration, the jump step would become obsolete. Once we have rateints data
from Stage 1 processing, we replace all NaN values in data and error arrays
with average values of their neighbouring pixels. We add these pixels to the
default bad-pixel map generated by the jwst pipeline. We performed a column-
by-column and row-by-row background subtraction to reduce $1/f$ noise from the
data. In this process, we subtracted a median of background pixels from each
row while we fitted a line to the column background pixels and subtracted the
estimated background from each column pixel. We then searched for cosmic ray
events in the data file by comparing each frame with a median frame. We
replaced all detected events with the mean of neighbouring pixels. However, we
added these events to the bad-pixel map in the end. We did not run Stage 2 of
the jwst pipeline because it does not change the science images.

Once we have corrected timeseries data, we used an open-source package
stark555https://stark-package.readthedocs.io/en/latest/ to extract spectra.
stark fits one and two-dimensional splines to the spectral data to find a
robust estimate of PSF which can later be used to extract the spectrum. Before
spectral extraction, we computed the location of the spectral trace using the
centre-of-flux method. We found that the location of the trace on the detector
remains extremely stable and varies only within 0.03 pixels. To estimate the
stellar spectrum, we first need to compute the stellar PSF, which we did by
fitting splines to the data. As a first approximation, we assume that the PSF
does not change with wavelength and with time, so we fitted a 1D spline to the
data as a function of distance from trace (known as pixel coordinates). This
is a poor assumption because while the PSF stays constant in time, it varies
significantly with wavelength. We improved our PSF estimate by fitting a 2D
spline to the data as a function of pixel coordinates and wavelength. This
robust PSF is then used to find stellar timeseries spectra. We used aperture
half-widths of 9 and 2 pixels to fit PSF and extract spectra, respectively. We
ran this procedure iteratively. At the end of each iteration, we subtracted
the median static residual noise from the raw data. The median static noise is
defined as a median difference between data and synthetic images constructed
using stellar PSF and spectra. Only two iterations were sufficient to find
robust stellar spectra. We compute the white-light light curve by taking a
weighted average of light curves in all spectroscopic channels between 3.8612
and 4.9771 $\mu$m. The raw white-light light curves for all visits are shown
in Figure 12.

Now that we have generated light curves we can fit an occultation model to the
data. The light curves show a strong ramp in the beginning of each visit (see,
Figure 12), so we discarded the first 35 min of the data before the analysis.
In the light curve analysis, we fixed all planetary parameters except
occultation depth to their values from the literature (Bourrier et al. 2018a;
Meier Valdés et al. 2022). We used a wide uniform prior between -500 to 500
ppm to the occultation depth parameter. We analysed white-light light curves
from all five visits together. We used juliet (Espinoza et al. 2019) to fit an
occultation model to the data, which uses an occultation model from batman
(Kreidberg 2015) and samples posteriors using dynesty (Speagle 2020). In
addition to the planetary model, we added linear and quadratic polynomials in
time to correct for long-term trends seen in the light curve. The best-fitted
values of white-light occultation depths are tabulated in Table 4. We could
not, however, model hour-long correlated noise (see, e.g., Fig. 1), with this
simple polynomial model. This is also evident from the Allan deviation plots,
shown in Fig. 13, of residuals that show additional noise at larger bin sizes.
The presence of uncorrected correlated noise means that the uncertainties
found on the occultation depths are underestimated. We could not determine the
origin of this noise: we searched engineering data but could not find any
parameter that correlates with the noise, pointing towards a possible
astrophysical origin. However, recent transit observations of a bright star
(GJ 341, K = 5.6 mag, Kirk et al. 2024) with the same observing mode also show
a similar noise as our dataset (see, their Fig. 2). So, the correlated noise
could be a previously unknown systematics of the instrument. We looked at the
2D spectral data at the group level to further test this possibility.
Generally, the data from the first and last groups are discarded as they could
be unreliable. We cannot do this since our dataset has only two groups. We
took the 2D spectral data for both groups independently and extracted spectral
timeseries from them in exactly the same manner described earlier. We finally
computed and analysed white-light light curves from both groups. We found that
the correlated noise similar to the integration level light curve is also
present at “group level” white-light light curves. This suggests that the
correlated noise does not originate from unreliable first and last groups (see
also our companion paper for more details, Patel & Brandeker, in prep).

We perform injection-retrieval tests on the white-light light curves to
estimate proper uncertainties on the occultation depths in the presence of
correlated noise. We first subtract the normalised planetary signal from the
raw white-light light curve keeping the long-term trend and the correlated
noise as it is in the data. We next produced 1000 realisations of light curves
by injecting an occultation signal at random times in the data. The depth of
the signal is equal to the median value from the full light curve analysis
presented earlier. In this process, we made sure that the full signal remained
inside the data. We fit a full model, consisting of an occultation model and
polynomial – linear and quadratic – trend, using juliet to each of the
realisations. We build a posterior of occultation depth using randomly
selected samples from the posteriors of occultation depth in each realisation.
These posteriors, shown in Fig. 14 for all visits, are clearly not Gaussian
for most of the cases illustrating the effect of correlated noise. A
68-percentile confidence interval of this posterior should be more
representative of uncertainties on white-light occultation depths. In the
cases where the uncertainties obtained this way were smaller than the “white”
uncertainties from the light curve analysis, we choose to report the larger
value.

The correlated noise is also present in the spectroscopic light curves of each
column. We first boosted the estimated errors of the spectroscopic light
curves and the white-light light curve according to the scatter in the light
curves. Then we divided spectroscopic light curves from each column with the
white-light light curve to remove the correlated noise from the spectroscopic
data. This mostly removed correlated noise from the spectroscopic light
curves. Finally, we computed relative occultation depths as
$1-(F_{\mathrm{in}}/F_{\mathrm{out}})$, where $F_{\mathrm{in}}$ and
$F_{\mathrm{out}}$ are the flux inside and outside of the occultation
duration, respectively. Before computing this, we made sure that the baseline
before and after the occultation signal was the same. Note that we compute
relative occultation depths at the native resolution of the instrument before
binning them to a lower resolution. This method minimizes the impact of any
leftover 1/f noise in the data (see, e.g., Espinoza et al. 2023).

####  A.1.2 Short-wave data analysis

The Stage 1 processing of 2.1 $\mu$m channel uncal files was mostly done in
the same way as for the 4.5 $\mu$m channel uncal files described above. The
main difference is that here we only perform a row-by-row background
subtraction. The short-wave PSF spreads to almost all pixel ranges along the
column so that there are very few background pixels along the column making it
impossible to perform background subtraction along columns.

Once we got rateints data, we performed simple aperture photometry to 2.1
$\mu$m channel data to obtain a photometric light curve. Before doing this, we
computed the centroids of the PSF using the centre-of-flux method. We then
computed a growth function – flux inside an aperture as a function of
increasing aperture radius – to optimally select an aperture radius. We find
that the growth function flattens out at around 45 pixel radius that we
eventually used in our analysis. We adapted the
photutils666https://photutils.readthedocs.io/en/stable/index.html (Bradley et
al. 2023) package to compute aperture photometry. photutils simply calculates
the total flux inside the aperture. Since we already did a row-by-row
background subtraction we did not perform another sky annulus subtraction.
Uncorrected short-wave photometric light curves are plotted in Figures 12.

We fitted an occultation model to thus-obtained SW light curves in almost the
same manner as for the occultation model fitting of LW white-light light
curves. The instrumental model used here was different from what was used in
the LW case. Here we used a linear polynomial in time and PSF centroids as
decorrelation vectors. Additionally, light curves from two of our visits
(Visits 1 and 4) show abrupt flux jumps analogous to what was found in
Schlawin et al. (2023) (see, Figure 12). These flux jumps may or may not be
caused by mirror tilting events as described in Schlawin et al. (2023) — a
thorough investigation of the origin of these jumps is ongoing (see also our
companion work Patel & Brandeker, in prep.). Here we model these flux jumps
using multiple step functions; since the jumps are abrupt and affect only a
few integrations, it is fairly easy to set the boundaries of step functions.
For certainty, we masked all integrations near jumps, which is safe because
the masked integrations consist of only a few per cent of the total number of
data points and none of these are near the ingress or egress. Another source
of noise in the SW light curves is the high-frequency periodic noise possibly
caused by the thermal cycling of heaters in the Integrated Science Instrument
Module on JWST (see, Espinoza et al. 2023). This is clearly visible in the
power spectrum of the light curve as a peak period near 3.8 min in all visits.
We performed a principal component analysis (PCA) of the PSF time series to
see if we could capture this noise as a principal component (PC) or not.
Indeed, one of the first PCs in all visits show a periodic pattern with a
period of about 3.8 min. While we are uncertain about the origin of this
noise, we simply use this PC as a decorrelation vector in our light curve
analysis.

In summary, our total model fitted to the SW light curve includes an
occultation model, linear models in time, PSF centroids and a PC. Step
functions were also included as decorrelation vectors in Visits 1 and 4. We
used juliet to fit the light curve data. The best-fitted occultation depths
can be found in Table 4. These data are also affected by a correlated noise
that we could not model using our simple model. This is also evident from the
Allan deviation of the residuals shown in Figure 13. We performed injection-
retrieval tests similar to the LW data analysis described in Appendix A.1.1 to
properly estimate the uncertainties on the occultation depths.

###  A.2 Eureka! — Reduction 1

Here we provide an independent reduction of the short-wave (SW) observations
of NIRCam. To reduce the nrca1 uncal files we used Eureka! (version
0.11.dev276+g4e12d23d, Bell et al. 2022) pipeline. Stage 1 consists of running
default jwst detector processing steps, but we skip the saturation step. On
stage 2 we only correct for the flat field. On Stage 3, we crop the full array
to a window between pixels 1400 and 2000 in the $x$-axis and between pixels 1
and 64 in the $y$-axis. We also mask pixels flagged as bad quality and reject
outliers above 7$\sigma$ along time axis. We interpolate bad pixels with a
linear function and perform row-by-row background subtraction and $1/f$ noise
correction. Aperture photometry is performed using a circular 40 pixel radius
aperture. We subtract the background region with an annulus with an inner edge
of 45 pixels and an outer edge of 60 pixels. Finally, Stage 4 uses the
calibrated files to produce the light-curve. Visit 1 and 4 exhibit strong
discontinuities, dividing the light-curve into five and six clearly defined
segments, respectively. To correct the discontinuities, first, we mask the
occultation. To flatten the light-curve, we fit a linear function to each
segment and then fit an occultation model with exoplanet in a Hamiltonian
Monte Carlo algorithm with PyMC3. The rest of the visits did not exhibit such
discontinuities and thus we fit only one linear function in time. The
resulting occultation depths are shown in Table 4. Compared to the stark
reduction and analysis, all occultation depths are consistent within
1$\sigma$.

###  A.3 Eureka! — Reduction 2

We produced an independent reduction of the NIRCam spectra using the jwst
(version 1.12.5, Bushouse et al. 2023) and Eureka! (version 0.9, Bell et al.
2022) pipelines, including purpose-built steps that we describe here. Starting
from the uncalibrated raw data, we ran the default jwst detector processing
steps up to (and including) the dark current step. Prior to the ramp fitting
step, we subtracted from each row the median of the left-most 650 pixels in
the corresponding row and group. By using these unilluminated pixels as a
reference of the level of noise added during readout, this helps reduce 1/f
noise. We then applied the remaining jwst calibration steps.

We ran the resulting calibrated files through Eureka!. We extracted columns
850 through 1945 and discarded the reference pixels. To straighten the trace,
we vertically slid each detector column by an integer number of pixels. We
performed background subtraction using the average value of each column,
rejecting $7\sigma$ outliers and excluding a window with a half-width of 15
pixels centred on the trace. Constructing the spatial profile from the median
frame, we performed optimal extraction on a region centred on the source and
with a half-width of 5 pixels. We generated 30 spectroscopic light curves
between 3.9365 and 4.9265 $\mu$m, each spanning 0.033 $\mu$m. In each light
curve, we discarded values farther than 4$\sigma$ from the mean of a sliding
window.

The flux in the light curves follows a downward trend with time, and they show
significant time-correlated noise. After trimming the initial 20 min of data,
where the ramp is the steepest, we modelled the white light curve in each
visit as the product of an exponential ramp, a linear polynomial and a batman
occultation model, where the occultation depth acted as a free parameter. The
fits included an estimated error multiplier to match the scatter in the
residuals. We assumed a circular orbit, and fixed the orbital period and mid-
transit time to the values in Zhang et al. (2021), and planet radius, orbital
inclination and scaled semi-major axis to those reported by Bourrier et al.
(2018a). For each visit, we also calculated the relative occultation depths
following the methodology outlined in Appendix A.1.1.

###  A.4 HANSOLO

The HANSOLO (atmospHeric trANsmission SpectrOscopy anaLysis cOde) pipeline was
originally developed to analyse ground-based transmission spectra observed
with 8m-class telescopes, but has been adapted to also enable its use on
NIRCam data (Lendl et al. 2016, 2017; Ahrer et al. 2023). HANSOLO takes
calibrated rateints outputs of the jwst pipeline Stage 1 as input.

We used the LACOSMIC algorithm (van Dokkum 2001) to remove cosmic ray effects
from the two-dimensional images and identified the spectral trace by using a
Moffat function fit to each column. The sky background was calculated on a
column-by-column basis by calculating a linear trend in the column background,
which was defined as at least 20 pixels away from the centre of the spectral
trace. This linear trend was then subtracted from the whole column. We
extracted the spectrum by summing over an aperture with a half-width of 4
pixels.

Consistent with the other reductions, we generated a white light curve and 30
spectroscopic light curves from which we clipped the first 35 min to remove
the worst of the ramp that is present in all the data. For each light curve we
applied a $5\sigma$ outlier rejection filter. We used the light curve and RV
fitting code CONAN to fit the white light curves with an occultation model and
a GP with a 3/2 Matern kernel to account for both the remaining ramp and the
correlated red noise. We leave the occultation depth and the GP parameters
(amplitude, lengthscale and a white noise factor) as free parameters and fix
all orbital parameters to the literature values found by Bourrier et al.
(2018a). The white light occultation depths are presented in Table 4. We then
calculate the common mode for each visit by removing the fitted occultation
from the white light curve and divide the common mode out of the spectroscopic
light curves. Since the spectroscopic light curves still show some correlated
noise even with the common mode removed, we then fit each spectroscopic light
curve individually in the same way as the white light curves, with the orbital
parameters held fixed and the occultation depth and GP parameters as free
parameters. The resulting emission spectra are shown in Figure 11.

###  A.5 transitspectroscopy

We take the corrected timeseries data from stark long-wave analysis and use an
open-source tool transitspectroscopy (Espinoza
2022)777https://github.com/nespinoza/transitspectroscopy. We first use a
centre of flux method to find the location of trace on the detector. We used
the optimal extraction algorithm from Marsh (1989) to extract 1D stellar
spectra from the timeseries data. In this procedure, we used an aperture half-
width of 3 pixels. The optimal extraction naturally clips all outliers not
identified by the pipeline. We masked all such 10$\sigma$ outliers. White-
light light curves for each visit were computed by taking a weighted average
of spectroscopic light curves between 3.8612 and 4.9771 $\mu$m.

We used juliet to fit the occultation model to the white-light light curve
data. In addition to the occultation model (from batman, Kreidberg 2015), our
full model includes linear, quadratic and cubic polynomials to model a long-
term decreasing trend. We also added white noise to the errors on the flux. We
fixed all planetary parameters except occultation depth from the literature
(Bourrier et al. 2018a; Meier Valdés et al. 2022). The median and
68-percentile confidence intervals for the best-fitted occultation depths are
tabulated in Table 4. We also determined relative occultation depth spectra
using the procedure described in Appendix A.1.1 and plotted in Fig. 11.

###  A.6 SPARTA

Our SPARTA reduction is very similar to that used in Hu et al. (2024), which
analyzed the one occultation observed by GO 1952 (PI Hu). The steps that we
used to go from the uncalibrated files to the spectroscopic light curves are
identical. In stage 1, we perform superbias subtraction, reference pixel
subtraction, non-linearity correction, dark subtraction, and up-the-ramp
fitting (which amounted to subtracting the two reads since we only have two).
In stage 2, we remove the background, which also removes some of the 1/f noise
because we perform row-by-row subtraction in addition to column-by-column
subtraction. In stage 3, we perform sum extraction with a window half-width of
2 pixels, obtaining spectroscopic light curves.

Using emcee, we fit the white light curve with a model that has the
occultation time and occultation depth as astrophysical free parameters, while
the light curve normalization factor, exponential ramp amplitude and
timescale, x and y linear correlation parameters, linear slope with time, and
error inflation multiple are free systematics parameters. We save the
systematics model corresponding to the best fit to the white light curve. To
fit the spectroscopic light curves, we first divide each light curve by the
aforementioned systematics model, and then fit the result with a model that
includes every parameter in the white light curve fit except the occultation
time (which we fix to the white light value).

Table 5: Spectroscopic parameters for 55 Cnc.

Method | $T_{\mathrm{eff}}$ | $\log g_{\star}$ | [Fe/H] | [Ca/H] | [Mg/H] | [Na/H] | $V\sin i$
---|---|---|---|---|---|---|---
| (K) | (dex) | (cgs) | (cgs) | (cgs) | (cgs) | (km s-1)
SME | $5234\pm 55$ | $4.33\pm 0.05$ | $+0.31\pm 0.05$ | $+0.33\pm 0.05$ | $+0.44\pm 0.12$ | $+0.60\pm 0.11$ | $2.0\pm 0.7$
astroARIADNE$a$$a$footnotetext:  | $5269\pm 46$ | $4.34\pm 0.07$ | $+0.34\pm 0.07$ | … | … | … | …

888 $a$$a$footnotetext: Posteriors from the SED modelling.

##  Appendix B Properties of the star

###  B.1 Observed stellar spectrum

We produced rateints files from uncalibrated data using the jwst pipeline
using the same procedure as described in Appendix A.1.1. We then ran Stage 2
of the jwst pipeline with some modifications, namely skipping the flat
fielding and extract1d steps, to produce calibrated spectrum files. This was
followed by correcting data and error files for NaN and cosmic rays as
described in Appendix A.1.1. Despite being classified as a point source by the
jwst pipeline, the physical unit of calibrated 2D spectrum data is given as
MJy/sr. We converted the units to Jy using the pixel area quoted in a header
file of calints data products from Stage 2 of the jwst pipeline. We finally
extracted the spectrum using stark as described in Appendix A.1.1. We
extracted a timeseries of spectra from part of the data from our most recent
visit, Visit 5. A median spectrum of these timeseries spectra is plotted in
Figure 15 and compared with the Crossfield (2012) empirical spectrum and black
body spectrum. We found that similar to Hu et al. (2024), the NIRCam observed
spectrum is discrepant with the Crossfield (2012) empirical spectrum. We think
that this may be because of improper photometric correction for bright stars
provided by the jwst pipeline. Furthermore, Hu et al. (2024) found that their
MIRI observed spectrum agrees very well with Crossfield (2012) spectrum. Here,
we use the Crossfield (2012) spectrum in our atmospheric retrieval analysis.

Figure 15: The observed stellar spectrum with NIRCam/JWST (in yellow) is shown
with Crossfield (2012) empirical spectrum and a blackbody at 5269 K.

###  B.2 Stellar parameters from modelling

We modelled 85 publically available spectra from the High Accuracy Radial
velocity Planet Searcher (HARPS; Mayor et al. 2003) spectrograph with a
resolution of 115 000. The spectra were co-added and modelled with
Spectroscopy Made Easy999http://www.stsci.edu/~valenti/sme.html (SME; Valenti
& Piskunov 1996; Piskunov & Valenti 2017) version 5.2.2 and the stellar
atmosphere grid Atlas12 (Kurucz 2013). SME computes synthetic spectra and
adjusts the chosen free parameters based on comparison with the observed
spectrum. We modelled one parameter at a time, utilising spectral features
sensitive to different photospheric parameters and iterating until all
parameters converged. Throughout the modelling, we held the macro- and micro-
turbulent velocities, $V_{\rm mac}$ and $V_{\rm mic}$, fixed at 2.7 km s-1
(Doyle et al. 2014) and 0.95 km s-1 (Bruntt et al. 2008). A description of the
modeling procedure is detailed in Persson et al. (2018). The results are
listed in Table 5.

The stellar radius was modelled with the SED fitting software
astroARIADNE101010https://github.com/jvines/astroARIADNE (Vines & Jenkins
2022) using priors from SME and photometry from the Johnson $B$ and $V$
magnitudes (APASS), $GG_{\rm BP}G_{\rm RP}$ (DR3), $JHK_{S}$ magnitudes
(2MASS), WISE W1-W2, and the Gaia DR3 parallax. We utilized three different
atmospheric model grids from Phoenix v2 (Husser et al. 2013), Castelli &
Kurucz (2004), and Kurucz (1993). The final radius was computed with Bayesian
Model Averaging and was found to be $0.953\pm 0.011$ $R_{\odot}$. The
luminosity is $0.63\pm 0.02$ $L_{\odot}$, and the visual extinction is
consistent with zero ($0.03\pm 0.03$). We derived a stellar mass of
$0.639^{+0.021}_{-0.020}$ $M_{\odot}$ interpolating the MIST (Choi et al.
2016) isochrones with astroARIADNE. Our results are very close to previous
results; von Braun et al. (2011) derive a stellar radius of $0.943\pm 0.010$
$R_{\odot}$ based on interferometric measurements and the parallax from van
Leeuwen (2007). Updating this calculation with the Gaia DR3 parallax, this
radius becomes $0.962\pm 0.010$ $R_{\odot}$ in good agreement with our
results.

##  Appendix C Detailed retrieval posterior distributions

In this appendix we present all posterior distributions from our retrieval
calculations for the CO/CO2 and SiO/SiO2/MgO cases. The posterior
distributions are ordered in chronological order and shown for the stark and
HANSOLO reductions. Due to the fact that for the HANSOLO reduction, the
retrievals are performed on absolute occultation depths, the posterior
distributions do not include the white-light occultation depths parameter
$d_{\mathrm{wl}}$.

It is also important to note that the depicted centre-log-ratio posterior
$\xi_{j}$ for the last molecule is not a free parameter in the retrieval as
mentioned in Sect. 2.3. Instead, we calculated the corresponding posterior
distribution following the requirement that for each posterior sample, the sum
of all $\xi$ values must be zero.

For Visits 1 and 3, the posterior distributions are already shown in Figs. 4 &
5 in the main text and are not repeated here. The corresponding posterior
spectra for the posteriors are shown in Fig. 3.

Figure 16: Posterior distributions of the free parameters for the first visit,
representing the CO/CO2-atmosphere scenario. Results are shown for the HANSOLO
reduction. The corresponding distributions for the stark reduction are
depicted in Fig. 7. We note that $\xi_{\mathrm{CO_{2}}}$ is not a free
parameter in the retrieval but was calculated during a postprocess procedure
following the requirement that in each posterior sample the sum of all $\xi$
values must be zero.

Figure 17: Posterior distributions of the free parameters for Visit 2,
representing the CO/CO2-atmosphere scenario. Results are shown for the stark
(left) and HANSOLO (right) reductions.

Figure 18: Posterior distributions of the free parameters for Visit 3,
representing the CO/CO2-atmosphere scenario. Results are shown for the stark
(left) and HANSOLO (right) reductions.

Figure 19: Posterior distributions of the free parameters for Visit 4,
representing the CO/CO2-atmosphere scenario. Results are shown for the stark
(left) and HANSOLO (right) reductions.

Figure 20: Posterior distributions of the free parameters for Visit 5,
representing the CO/CO2-atmosphere scenario. Results are shown for the stark
(left) and HANSOLO (right) reductions.

Figure 21: Posterior distributions of the free parameters for Visit 1,
representing the SiO/SiO2/MgO-atmosphere scenario. Results are shown for the
stark (left) and HANSOLO (right) reductions.

Figure 22: Posterior distributions of the free parameters for Visit 2,
representing the SiO/SiO2/MgO-atmosphere scenario. Results are shown for the
stark (left) and HANSOLO (right) reductions. Figure 23: Posterior
distributions of the free parameters for Visit 3, representing the
SiO/SiO2/MgO-atmosphere scenario. Results are shown for the HANSOLO reduction.
The corresponding distributions for the stark reduction are depicted in Fig.
5.

Figure 24: Posterior distributions of the free parameters for Visit 4,
representing the SiO/SiO2/MgO-atmosphere scenario. Results are shown for the
stark (left) and HANSOLO (right) reductions.

Figure 25: Posterior distributions of the free parameters for Visit 5,
representing the SiO/SiO2/MgO-atmosphere scenario. Results are shown for the
stark (left) and HANSOLO (right) reductions.