Context. The nature of the closein rocky planet 55 Cnce is puzzling despite it 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 32 spinorbit resonance and 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. JWSTNIRCam Near Infrared Camera observed five occultations secondary eclipses of the planet  of which four were observed within a week  simultaneously at 2.1 and 4.5 m. While the former gives bandintegrated photometry the latter provides a spectrum between 3.95.0 m. Results. We find that the occultation depths in both bandpasses are highly variable and change between a nondetection 5  6 ppm and 7  9 ppm to 96  8 ppm and 119SUB19SUBSUP34SUP ppm at 2.1 m and 4.5 m respectively. Interestingly the variations in both bandpasses are not correlated and do not support the 32 spinorbit resonance explanation. The measured brightness temperature at 4.5 m varies between 8732256 K and is lower than the expected dayside temperature of bare rock with no heat redistribution 2500 K which is indicative of an atmosphere. Our atmospheric retrieval analysis of occultation depth spectra at 4.5 m finds that different visits statistically favour various atmospheric scenarios including a thin outgassed COCOSUB2SUB 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 COCOSUB2SUB which is also hinted at by retrievals. Alternatively the variability observed at both 2.1 and 4.5 m could be the result of a circumstellar patchy dust torus generated by volcanism on the planet.
['techniques: photometric', 'techniques: spectroscopic', 'planets and satellites: atmospheres', 'planets and satellites: terrestrial planets', 'planets and satellites: individual: 55 Cnc e', 'Astrophysics - Earth and Planetary Astrophysics']
JWST reveals the rapid and strong dayside variability of 55 Cancri e
{'JWST reveals a rapid and strong day side variability of 55 Cancri e': "J. A. Patel 1 ⋆ , A. Brandeker 1 , D. Kitzmann 2 , 3 , D. J. M. Petit dit de la Roche 4 , A. Bello-Arufe 5 , K. Heng 6 , 7 , 8 , 9 , E. Meier Valdés 3 , C. M. Persson 10 , M. Zhang 11 , B.-O. Demory 3 , 2 , V. Bourrier 4 , A. Deline 4 , D. Ehrenreich 4 , 12 ,  \nM. Fridlund 13 , 10 , R. Hu 5 , 14 , M. Lendl 4 , A. V. Oza 14 , 5 , Y. Alibert 2 , 3 , M. J. Hooton 15  \n- 1 Department of Astronomy, Stockholm University, AlbaNova University Center, 10691 Stockholm, Sweden\n- 2 Weltraumforschung und Planetologie, Physikalisches Institut, Universität Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland\n- 3 Center for Space and Habitability, Universität Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland\n- 4 Observatoire astronomique de l'Université de Genève, Chemin Pegasi 51, 1290 Versoix, Switzerland\n- 5 Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91011, USA\n- 6 Faculty of Physics, Ludwig Maximilian University, Scheinerstrasse 1, D-81679, Munich, Bavaria, Germany\n- 7 ARTORG Center for Biomedical Engineering Research, University of Bern, Murtenstrasse 50, CH-3008, Bern, Switzerland\n- 8 University College London, Department of Physics & Astronomy, Gower St, London, WC1E 6BT, United Kingdom\n- 9 University of Warwick, Department of Physics, Astronomy & Astrophysics Group, Coventry CV4 7AL, United Kingdom\n- 10 Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, 439 92 Onsala, Sweden\n- 11 Department of Astronomy and Astrophysics, The University of Chicago, Chicago, IL 60637, USA\n- 12 Centre Vie dans l'Univers, Faculté des sciences de l'Université de Genève, Quai Ernest-Ansermet 30, 1205 Geneva, Switzerland\n- 13 Leiden Observatory, University of Leiden, PO Box 9513, 2300 RA Leiden, The Netherlands\n- 14 Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125, USA\n- 15 Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK  \nReceived 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.  \nAims. We wish to explore the possibility that the variability originates from the planet being in a 3:2 spin-orbit resonance, thus showing di ff erent sides during occultations. We proposed and were awarded Cycle 1 time at the James Webb Space Telescope (JWST) to test this hypothesis.  \nMethods. 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 µ m. While the former gives band-integrated photometry, the latter provides a spectrum between 3.95.0 µ m.  \nResults. We find that the occultation depths in both bandpasses are highly variable and change between a non-detection ( -5 ± 6 ppm and 7 ± 9 ppm) to 96 ± 8 ppm and 119 + 34 -19 ppm at 2.1 µ mand 4.5 µ 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 µ mvaries 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 µ m finds that di ff erent 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.  \nConclusions. 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 µ m, could be the result of a circumstellar patchy dust torus generated by volcanism on the planet.  \nKey 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  \n( 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.  \n55 Cnc e was discovered by McArthur et al. (2004) with an orbital period of ∼ 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 es-  \nmates 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 supercritical 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 ⊕ ) and radius (1 . 88 R ⊕ ) 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.  \nSoon after its discovery, Demory et al. (2012) used Spitzer to detect thermal emission from 55 Cnc e and determined its dayside 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 µ m. They calculated the average day-side temperature to be around 2350 K with a maximum ∼ 2700 K. Curiously, the hottest location of the planet was found to be shifted by ∼ 41 · to the east compared to the sub-stellar point, indicating a strong heat re-distribution. On the other hand, the day-night temperature di ff erence was found to be as large as 1300 K, a sign of ine ffi cient heat transport to the night side. These conflicting results led Demory et al. (2016a) to speculate that perhaps e ffi cient 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 ine ffi ciency 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 e ffi cient 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 o ff set; 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 µ 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 di ff erence with a smaller phase o ff set, more consistent with a poor heat transport typically found on USPs.  \nThe heavyweight atmosphere on the planet, which was implied by the Spitzer phase curve, climate modelling and massradius constraints, is challenging to detect. Numerous observa-  \ntions 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 α 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).  \nThe 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 µ 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 o ff set 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).  \nMultiple 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 µ m. Gas plumes evolving at di ff erent 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  \nits 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-CO 2 atmospheres are outgassed under a broad range of conditions (surface pressures, oxygen fugacity, temperatures).  \nSince 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.  \nAlthough multiple hypotheses have been provided to describe the thermal phase curve and variability from 55 Cnc e, each has di ffi culties 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.", '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 spinorbit 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 di ff erent 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  \nof 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.  \nAssuming 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 di ff erent sides of the planet they will show up di ff erently 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 di ff erent 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 ( ∼ 35 . 5 hr), close to the period of 2 days.  \nTo 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 timeconstrained 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 di ffi culties, 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 di ff erent 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.  \nNIRCam o ff ers simultaneous observations in short-wave (SW) and long-wave (LW) channels at 0.6-2.3 µ m and 2.45.0 µ 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  \nTable 1. Observation log and wide band occultation depths  \n3.9-5 µ m (centred at around 4.5 µ m) in the LW channel (or, 4.5 µ m channel) and one single photometric data point in a narrow-band (2.3%) bandpass at 2.12 µ m from the SW channel (also referred to as the 2.1 µ 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.  \nWe used five independent pipelines to reduce and analyse the spectroscopic data at 4.5 µ m and two di ff erent 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 di ff er in their treatment of correlated noise and thus produce slightly di ff erent results, which was the reason for choosing two di ff erent reductions for retrieval (see Appendix A for more details).  \nTo 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 di ff erent 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 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 (CO 2 ) (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 (SiO 2 ), and magnesium oxide (MgO) (Zilinskas et al. 2022).  \nNested 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 re-  \nTable 2. Summary of retrieval parameters and prior distributions used for the retrieval models.  \ntrieval calculation. Increasing this value further did not alter the resulting evidence values to a significant degree.  \nThe atmosphere was considered to be isothermal with the surface pressure p surf as a free parameter in the retrieval model. The atmosphere and surface were allowed to have their own distinct temperatures, T atm and T surf, respectively.  \nThe cross sections of CO, CO 2 , SiO, SiO 2 , 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).  \nThe 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 xj of a given molecule j in the mixture is given by  \nξ j = clr( xj ) = ln xj g ( x ) , (1)  \nwhere g ( x ) is the geometric mean of all mixing ratios x :  \ng ( x ) =         n Y j = 1 xj         1 / n . (2)  \nDue to the constraint that  \nn X j = 1 xj = 1 or n X j = 1 ξ j = 0 , (3)  \nFig. 1. Detrended occultation light curves from the short-wave photometric channel (2.1 µ m, left panel) and white-light light curves from the long-wave channel (4.5 µ 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.  \n<!-- image -->  \nFig. 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.  \n<!-- image -->  \nonly n -1 free parameters are needed in the retrieval. We used uniform priors to produce ξ j values subject to the constraints that min( x ) = 10 -10 and max ( x ) = 1, see Benneke & Seager (2012)  \nfor details. We note that the prior boundaries for ξ j depend on the number of molecules in the retrieval and the chosen value of the smallest allowed mixing ratio.  \nFor 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 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.  \nAll retrieval parameters for the di ff erent 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 wavelengthdependent occultation depths.", '3.1. Wide band occultation depths': 'We used six pipelines to reduce and fit our JWST / NIRCam dataset. The methods are described in detail in Appendix A. Here  \nTable 3. Overview of retrieval results for the stark and HANSOLO reductions. Boldface indicates the statistically preferred models.  \nwe 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.  \nOur 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 µ 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 µ m continuously increased from basically non-detection in Visit 1 (7 ± 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 ∼ 95 ± 16 ppm and consistent with the depths from Visits 3 and 4 but di ff ers 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.  \nWe 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 function 1 to compute brightness temperatures using the measured white-light occultation depths at 4.5 µ 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.  \nSimilarly, the 2.1 µ m channel occultation depths are also variable. Within a week, the 2.1 µ moccultation 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 ± 8 ppm, which is almost equal to the depth observed at 4.5 µ min 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 µ m (Fig. 2). Fig. 1 present the de-trended SW data with best-fitted models.  \nThe 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.  \nAll visits showed various degrees of significant correlated noise of unknown origins, in both the 2.1 and 4.5 µ m channels. The leftover correlated noise can be seen in Fig. 1 and are also quantified in the Allan deviation plots in Fig. A.3. 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 di ff erent methods to handle the correlated noise could lead to di ff erent 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. A.1, are also variable from visit to visit and do not show any consistent spectral features.', '3.2.1. Summary of the retrieval results': 'The retrieval results for the two di ff erent reductions across all five visits and for the four di ff erent model scenarios described in Sect. 2.3 are summarised in Table 3. The table shows the resulting Bayesian evidence values ln 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.  \nThe 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.  \nHowever, 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,  \nFig. 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σ 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σ 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.  \n<!-- image -->  \non the other hand, can be fit with any of the four models. There seems to be little statistical preference for any of the di ff erent modelling scenarios.  \nThe results for the stark reduction show a much broader range of di ff erent models that are statistically preferred. As suggested by Table 3, the first visit strongly prefers a CO / CO 2 atmosphere, the second visit can be explained by either a CO / CO 2 atmosphere, a planetary blackbody, or a siliciate vapour atmosphere, while the third model overwhelmingly prefers the SiO / SiO 2 scenario. The fourth visit is consistent with a planetary blackbody spectrum, as well as an atmosphere with CO & CO 2 , or SiO, SiO 2 &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 & CO 2 ) and Visit 3 (SiO, SiO 2 & 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 surf and the abundances of CO and CO 2 . 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 CO 2 . For comparison, if the outgassing flux were to be balanced by flux-limited atmospheric escape then the implied surface pressure is ∼ 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 wl  \nFig. 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 ξ 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 ξ values must be zero.  \n<!-- image -->  \nis shifted from its prior value of 7 ± 9 ppm, though they are both still within their 1σ intervals.  \nThe posterior distribution for the SiO / SiO 2 / 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 SiO 2 with only an upper limit for SiO and essentially no constraints on MgO. The posterior spectra shown in Fig. 3 clearly show the drop-o ff in the occulation depth near a wavelength of 4.8 µ m caused by SiO 2 . Just like in the previous CO / CO 2 scenario for Visit 1, the retrieved atmosphere temperature is again much higher than the one of the surface.', '3.2.3. Blackbody temperatures': 'The resulting posterior distributions of the blackbody temperature models are shown in Fig. 6 for all visits and the two di ff erent 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.  \nAs the distributions depicted in the figure suggest, the temperatures retrieved from the HANSOLO observational data are found in two di ff erent clusters. A low-temperature mode near  \n750 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σ intervals usually in the range of about ± 100 K, despite the rather large errors on the observational data points (see Fig. 3).  \nFor 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σ 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': 'For the two model scenarios that involve atmospheres, we also retrieve the surface pressure. For the CO / CO 2 model, the corresponding posterior distributions are shown in Fig. 7, while those for the SiO / SiO 2 / MgO scenario are shown in Fig. 8.  \nIn 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  \nFig. 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 ξ 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 ξ values must be zero.  \n<!-- image -->  \nvisit, 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 di ffi cult to provide good constraints for actual atmospheric models.  \nFor 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 / CO 2 model (see upper panel of Fig. 7) and Visits 1 and 3 for the SiO / SiO 2 / MgO case (see upper panel of Fig. 8).  \nOther 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 CO 2 , as discussed above.  \nWe note that our retrieved surface pressures di ff er 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 di ff erent results regarding the atmospheric properties, this is not too surprising. Furthermore, Hu et al. (2024) employed a di ff erent 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 a ff ect the posterior distributions of the surface pressure.  \nFig. 6. Retrieved temperatures for all five visits using the blackbody model. Top panel: results for the stark reduction, bottom panel: HANSOLO .  \n<!-- image -->  \nFig. 7. Surface pressure posterior distribution from the CO / CO2 model for all five visits. Top panel: results for the stark reduction, bottom panel: HANSOLO .  \n<!-- image -->  \nFig. 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 .  \n<!-- image -->  \nFig. 9. Posterior distributions for the atmosphere (top) and surface temperatures (bottom). The distributions are shown for the CO / CO2 model and the stark reduction.  \n<!-- image -->', '3.2.5. Surface and atmosphere temperatures': 'For the CO / CO 2 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.  \nJust 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.  \nVisit 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 di ff erent 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.  \nThere 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 e 2 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 e ffi ciency (evolved BSE composition).  \nFig. 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 µ m channel (in blue) and white-light occultation depths for 4.5 µ 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 di ff erent shapes representing bulk-silicate (stars) and evolved bulksilicate (squares) compositions.  \n<!-- image -->  \nAll of their models with di ff erent outgassing e ffi ciencies predict occultation depths of 70-80 and 145-150 ppm for NIRCam 2.1 and 4.5 µ mchannels. As depicted in Fig. 10 these values are larger compared to our observations. Some occultation depths are, however, consistent with models at 1-3 σ . One occultation depth at 2.1 µ 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 / SiO 2 / MgO model in the retrieval analysis. The bandaveraged occultation depth for this visit at 4.5 µ m agrees with the model prediction (145 ppm for BSE case) at 2.4 σ . However, the short-wave occultation depth in this visit is inconsistent with the model prediction at 7 σ . 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.  \nAt 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 redistribution (e.g., Hammond & Pierrehumbert 2017). The estimated dayside brightness temperatures (see, Table 1) at 4.5 µ m (Table 1) in all visits are smaller than the expected dayside temperature 3 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 µ min 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 di ffi cult to explain. Such outgassed atmospheres are incapable of producing occultation depths as high as ≈ 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.  \nFor the first data reduction ( stark ), the outgassed atmosphere with CO and CO 2 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.  \nThe 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 ∼ 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.  \nIf 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 ∼ 1 µ bar. The surface temperature is ∼ 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 ∼ 2000 K to ∼ 2500 K to ∼ 1500 K from Visits 1 to 2 to 3. While this is not implausible because of the short radiative  \ntimescales, 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 µ mand Visit 2 at 2.1 µ m, show occultation depths that are consistent with zero at 1σ . These non-detections are challenging to explain with any kind of atmospheric phenomena. Moreover, the occultation depths observed at 2.1 µ m and 4.5 µ m are not correlated with each other (Fig. 2), which potentially hints towards di ff erent origins of variability in di ff erent wavelength channels.  \nA grey absorber could explain the optical and 2.1 µ 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 su ffi ciently 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 µ mchannel. 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 / IRACbandpass at 4.5 µ mandour 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 CO 2 which also absorbs at 4.5 µ m) atmosphere with several tens of µ 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. ).  \nThe variability at 2.1 µ mis di ffi cult to explain with a thin atmosphere consisting volcanic gases such as SO 2 or CO 2 since they do not have significant opacity in the 2.1 µ 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 µ mfrom the size range of 0.1-0.7 µ mdiscussed in Morris et al. (2021) and Meier Valdés et al. (2023), the particles will be opaque in the 2.1 µ m channel, but transparent in the 4.5 µ mchannel. 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 × 10 6 kg s -1 is within a factor of two of the maximum escape rate derived by CHEOPS, reported to be as large as ∼ 2.9 × 10 6 kg s -1 (Meier Valdés et al. 2023). If the particle size is larger than 0.7 µ m, they can, in principle, even explain the variability at 4.5 µ m channel. However, the non-correlation of occultation depths at 2.1 µ m and 4.5 µ 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 µ m, respectively, as mentioned above. However, the e ff ect of the dust torus on the  \ntransit 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 a ff ect 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.  \nIt 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-o ff (A ff olter 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 tidalheating (Oza et al. 2019a; Quick et al. 2020; Charnoz et al. 2021).  \nThe 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. E ff ectively, 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 ∼ 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 di ff erence 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 GO1952: PI Hu). In particular, we test the hypothesis that 55 Cnc e is in a 3:2 spin-orbit resonance, thus showing di ff erent 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.  \nWe observed five occultations of 55 Cnc e in two wavelength bands, or channels, a spectroscopic band at 4.5 µ m and a single photometric band at 2.1 µ 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 di ff erent 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 µ m channels, but is curiously not correlated between channels. The estimated brightness temperature at 4.5 µ mvaries 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.  \nThe spectroscopic data at 4.5 µ m is a ff ected by correlated noise of unknown origin. Although the results from di ff erent reductions overall agree well with each other, there are several di ff erences in white-light occultation depths and emission spectra that can be attributed to di ff erent 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 / CO 2 or SiO / SiO 2 / 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 / CO 2 or SiO / SiO 2 / MgO were not discarded either, statistically. The retrievals with stark prefer CO / CO 2 atmospheres in at least two visits, SiO / SiO 2 / MgO atmosphere in one visit and blackbody and flat line models in the remaining two visits. The CO / CO 2 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 / CO 2 atmosphere.  \nThe occultation depth variability in the 2.1 µ mchannel, especially its uncorrelated behaviour with its 4.5 µ mchannel counterpart, is challenging to explain with a simple atmospheric model. It is possible that the variability seen at 2.1 µ mand that at 4.5 µ m have di ff erent origins. A circumstellar inhomogeneous cloud of dust could potentially describe the variability at 2.1 µ 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  \nobservations 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 µ m where volcanic gases CO 2 / SO 2 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.  \nWhile 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.  \nAcknowledgements. 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 S pice D une , grant agreement No 947634, and F our A ces ; 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': '- A ff olter, L., Mordasini, C., Oza, A. V., Kubyshkina, D., & Fossati, L. 2023, A&A, 676, A119\n- Ahrer, E.-M., Stevenson, K. B., Mansfield, M., et al. 2023, Nature, 614, 653 Angelo, I. & Hu, R. 2017, AJ, 154, 232\n- Bell, T., Ahrer, E.-M., Brande, J., et al. 2022, The Journal of Open Source Software, 7, 4503\n- Benneke, B. & Seager, S. 2012, ApJ, 753, 100\n- Bourrier, V., Dumusque, X., Dorn, C., et al. 2018a, A&A, 619, A1\n- Bourrier, V., Ehrenreich, D., Lecavelier des Etangs, A., et al. 2018b, A&A, 615, A117\n- Bradley, L., Sip"ocz, B., Robitaille, T., et al. 2023, astropy / photutils: 1.8.0\n- Brandeker, A. 2019, in AAS / Division for Extreme Solar Systems Abstracts, Vol. 51, AAS / Division for Extreme Solar Systems Abstracts, 311.07\n- 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  \nArticle number, page 14 of 32  \nBrogi, M., Keller, C. U., de Juan Ovelar, M., et al. 2012, A&A, 545, L5 Bruntt, H., De Cat, P., & Aerts, C. 2008, A&A, 478, 487 Burrows, A. S. 2014, Proceedings of the National Academy of Science, 111,  \n- 12601\n- Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2023, JWST Calibration Pipeline\n- Castelli, F. & Kurucz, R. L. 2004, astro-ph / 0405087 [ astro-ph/0405087 ]\n- Charnoz, S., Sossi, P. A., Lee, Y.-N., et al. 2021, Icarus, 364, 114451\n- Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102\n- Crida, A., Ligi, R., Dorn, C., & Lebreton, Y. 2018, ApJ, 860, 122\n- Crossfield, I. J. M. 2012, A&A, 545, A97\n- Dawson, R. I. & Fabrycky, D. C. 2010, ApJ, 722, 937\n- Deibert, E. K., de Mooij, E. J. W., Jayawardhana, R., et al. 2021, AJ, 161, 209\n- Demory, B.-O., Gillon, M., de Wit, J., et al. 2016a, Nature, 532, 207\n- Demory, B. O., Gillon, M., Deming, D., et al. 2011, A&A, 533, A114\n- Demory, B.-O., Gillon, M., Madhusudhan, N., & Queloz, D. 2016b, MNRAS, 455, 2018\n- Demory, B.-O., Gillon, M., Seager, S., et al. 2012, ApJ, 751, L28  \nDemory, B. O., Sulis, S., Meier Valdés, E., et al. 2023, A&A, 669, A64  \n- Dorn, C., Hinkel, N. R., & Venturini, J. 2017, A&A, 597, A38\n- Doyle, A. P., Davies, G. R., Smalley, B., Chaplin, W. J., & Elsworth, Y. 2014, MNRAS, 444, 3592\n- 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\n- Ehrenreich, D., Bourrier, V., Bonfils, X., et al. 2012, A&A, 547, A18\n- Espinoza, N. 2022, TransitSpectroscopy\n- Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262\n- Espinoza, N., Úbeda, L., Birkmann, S. M., et al. 2023, PASP, 135, 018002\n- Esteves, L. J., de Mooij, E. J. W., Jayawardhana, R., Watson, C., & de Kok, R. 2017, AJ, 153, 268\n- Feroz, F. & Hobson, M. P. 2008, MNRAS, 384, 449\n- Fischer, D. A., Marcy, G. W., Butler, R. P., et al. 2008, ApJ, 675, 790\n- Folsom, C. P., Ó Fionnagáin, D., Fossati, L., et al. 2020, A&A, 633, A48\n- Gebek, A. & Oza, A. V. 2020, MNRAS, 497, 5271\n- Gillon, M., Demory, B. O., Benneke, B., et al. 2012, A&A, 539, A28\n- Grimm, S. L. & Heng, K. 2015, ApJ, 808, 182\n- Grimm, S. L., Malik, M., Kitzmann, D., et al. 2021, ApJS, 253, 30\n- Hammond, M. & Pierrehumbert, R. T. 2017, ApJ, 849, 152\n- Heng, K. 2023, ApJ, 956, L20\n- Hu, R., Bello-Arufe, A., Zhang, M., et al. 2024, arXiv e-prints, arXiv:2405.04744\n- 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  \nHusser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6  \nJindal, A., de Mooij, E. J. W., Jayawardhana, R., et al. 2020, AJ, 160, 101  \nKeles, E., Mallonn, M., Kitzmann, D., et al. 2022, MNRAS, 513, 1544  \n- Kirk, J., Stevenson, K. B., Fu, G., et al. 2024, AJ, 167, 90\n- Kite, E. S., Fegley, Bruce, J., Schaefer, L., & Gaidos, E. 2016, ApJ, 828, 80  \nKitzmann, D., Heng, K., Oreshenko, M., et al. 2020, ApJ, 890, 174  \nKoll, D. D. B., Malik, M., Mansfield, M., et al. 2019, ApJ, 886, 140  \nKreidberg, L. 2015, PASP, 127, 1161  \n- Kreidberg, L., Koll, D. D. B., Morley, C., et al. 2019, Nature, 573, 87\n- Kurucz, R. L. 1993, VizieR Online Data Catalog, VI / 39\n- Kurucz, R. L. 2013, ATLAS12: Opacity sampling model atmosphere program, Astrophysics Source Code Library\n- Lanza, A. F. 2021, A&A, 653, A112  \nLendl, M., Cubillos, P. E., Hagelberg, J., et al. 2017, A&A, 606, A18  \n- Lendl, M., Delrez, L., Gillon, M., et al. 2016, A&A, 587, A67\n- Li, G., Gordon, I. E., Rothman, L. S., et al. 2015, ApJS, 216, 15\n- Li, H. Y., Tennyson, J., & Yurchenko, S. N. 2019, MNRAS, 486, 2351 Lopez, E. D. 2017, MNRAS, 472, 245\n- Madhusudhan, N., Lee, K. K. M., & Mousis, O. 2012, ApJ, 759, L40\n- Mahapatra, G., Helling, C., & Miguel, Y. 2017, MNRAS, 472, 447 Marsh, T. R. 1989, PASP, 101, 1032\n- Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20\n- McArthur, B. E., Endl, M., Cochran, W. D., et al. 2004, ApJ, 614, L81\n- Meier, T. G., Bower, D. J., Lichtenberg, T., Hammond, M., & Tackley, P. J. 2023, A&A, 678, A29\n- Meier Valdés, E. A., Morris, B. M., Demory, B. O., et al. 2023, A&A, 677, A112 Meier Valdés, E. A., Morris, B. M., Wells, R. D., Schanche, N., & Demory, B. O. 2022, A&A, 663, A95\n- Mercier, S. J., Dang, L., Gass, A., Cowan, N. B., & Bell, T. J. 2022, AJ, 164, 204 Meyer zu Westram, M., Oza, A. V., & Galli, A. 2024, Journal of Geophysical Research (Planets), 129, e2023JE007935  \nMorris, B. M., Delrez, L., Brandeker, A., et al. 2021, A&A, 653, A173  \n- of Ockham, W. 1495, Quaestiones et decisiones in quattuor libros Sententiarum Petri Lombardi: Centilogium theologicum (Johannes Trechsel)\n- Owens, A., Conway, E. K., Tennyson, J., & Yurchenko, S. N. 2020, MNRAS, 495, 1927  \nOza, A., Bower, D. J., Demory, B.-O., et al. 2019a, in EPSC-DPS Joint Meeting  \n2019, Vol. 2019, EPSC-DPS2019-1714 Oza, A. V., Johnson, R. E., Lellouch, E., et al. 2019b, ApJ, 885, 168 Persson, C. M., Fridlund, M., Barragán, O., et al. 2018, A&A, 618, A33 Piskunov, N. & Valenti, J. A. 2017, A&A, 597, A16 Quick, L. C., Roberge, A., Mlinar, A. B., & Hedman, M. M. 2020, PASP, 132, 084402 Rasmussen, K. C., Currie, M. H., Hagee, C., et al. 2023, AJ, 166, 155 Ridden-Harper, A. R., Snellen, I. A. G., Keller, C. U., et al. 2016, A&A, 593, A129 Rodríguez, A., Callegari, N., Michtchenko, T. A., & Hussmann, H. 2012, MNRAS, 427, 2239 Salz, M., Czesla, S., Schneider, P. C., & Schmitt, J. H. M. M. 2016, A&A, 586, A75 Schlawin, E., Beatty, T., Brooks, B., et al. 2023, PASP, 135, 018001 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, J. S. 2020, MNRAS, 493, 3132 Sulis, S., Dragomir, D., Lendl, M., et al. 2019, A&A, 631, A129 Tabernero, H. M., Allende Prieto, C., Zapatero Osorio, M. R., et al. 2020, MNRAS, 498, 4222 Tamburo, P., Mandell, A., Deming, D., & Garhart, E. 2018, AJ, 155, 221 Tian, M. & Heng, K. 2024, ApJ, 963, 157 Trotta, R. 2008, Contemporary Physics, 49, 71 Trotta, R. 2017, arXiv e-prints, arXiv:1701.01467 Tsiaras, A., Rocchetto, M., Waldmann, I. P., et al. 2016, ApJ, 820, 99 Valenti, J. A. & Piskunov, N. 1996, A&AS, 118, 595 van Dokkum, P. G. 2001, PASP, 113, 1420 van Leeuwen, F. 2007, A&A, 474, 653 Vines, J. I. & Jenkins, J. S. 2022, arXiv e-prints, arXiv:2204.03769 von Braun, K., Boyajian, T. S., ten Brummelaar, T. A., et al. 2011, ApJ, 740, 49 Winn, J. N., Matthews, J. M., Dawson, R. I., et al. 2011, ApJ, 737, L18 Winn, J. N., Sanchis-Ojeda, R., & Rappaport, S. 2018, New A Rev., 83, 37 Yurchenko, S. N., Mellor, T. M., Freedman, R. S., & Tennyson, J. 2020, MNRAS, 496, 5282  \nYurchenko, S. N., Tennyson, J., Syme, A.-M., et al. 2022, MNRAS, 510, 903 Zhang, M., Knutson, H. A., Wang, L., et al. 2021, AJ, 161, 181 Zieba, S., Zilinskas, M., Kreidberg, L., et al. 2022, A&A, 664, A79 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 A.1, we summarise the whitelight occultation depths between about 4 and 5 µ m (see, below for exact wavelength range for di ff erent methods) and photometric occultation depths at 2.1 µ m. Figure A.1 compares the relative occultation depth spectra for all visits from di ff erent methods. It can be seen from Figure A.1 and Table A.1 that the results obtained with various independent analysis methods overall agree with each other, however, there are some di ff erences which could be attributed to the di ff erent 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.', 'Appendix 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 µ mchannel) spectroscopic channel and a short-wave (SW, 2.1 µ mchannel) photometric channel. We analysed both datasets with our pipeline.', 'Appendix A.1.1: Long-wave data analysis': "We downloaded uncalibrated data files ( uncal files) from the MAST archive and used the o ffi cial 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-byrow 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.  \nOnce we have corrected timeseries data, we used an opensource package stark 4 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  \nonly 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 di ff erence between data and synthetic images constructed using stellar PSF and spectra. Only two iterations were su ffi cient to find robust stellar spectra. We compute the whitelight light curve by taking a weighted average of light curves in all spectroscopic channels between 3.8612 and 4.9771 µ m. The raw white-light light curves for all visits are shown in Figure A.2.  \nNow 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 A.2), 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 A.1. 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. A.3, 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' whitelight 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 ).  \nWe perform injection-retrieval tests on the white-light light curves to estimate proper uncertainties on the occultation depths  \nTable A.1. Comparison of white-light and photometric occultation depths from di ff erent methods  \nNotes. The uncertainties are 68 percentile of the corresponding posterior distribution. Visit 4 is the archival observation from Hu et al. (2024).  \nFig. A.1. Comparison of occultation depth spectra for all observations from di ff erent 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).  \n<!-- image -->", '2.1 um channel': "<!-- image -->  \n<!-- image -->  \nFig. A.2. Raw photometric light curves from the short-wave channel at 2.1 µ m (in blue) and raw white-light light curves from the long-wave channel at 4.5 µ 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.  \n<!-- image -->  \nFig. A.4. Posteriors of occultation depths from injection-retrieval exercise (see, text) for 2.1 µ m (SW) channel (the top row, in blue) and 4.5 µ 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.  \n<!-- image -->  \nFig. A.3. Allan deviation plots of residuals from photometric light curve analysis from 2.1 µ m(SW) channel (left panel, in blue) and 4.5 µ m(LW) channel white-light light curve analysis (right panel, in orange).  \n<!-- image -->  \nOccultation depth  \n<!-- image -->  \nin 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 realisa-  \ntion. These posteriors, shown in Fig. A.4 for all visits, are clearly not Gaussian for most of the cases illustrating the e ff ect 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.  \nThe 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 whitelight 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 in / F out), where F in and F 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).", 'Appendix A.1.2: Short-wave data analysis': 'The Stage 1 processing of 2.1 µ m channel uncal files was mostly done in the same way as for the 4.5 µ m channel uncal files described above. The main di ff erence 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.  \nOnce we got rateints data, we performed simple aperture photometry to 2.1 µ 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 photutils 5 (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 A.2.  \nWe 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 di ff erent 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 A.2). 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 a ff ect only a few integrations, it is fairly easy to set the boundaries of step functions. For certainty, we masked all  \nintegrations 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.  \nIn 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 A.1. These data are also a ff ected 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 A.3. 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.', 'Appendix 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 σ 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 A.1. Compared to the stark reduction and analysis, all occultation depths are consistent within 1 σ .', 'Appendix 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  \nthe 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.  \nWe 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 σ 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 µ m, each spanning 0.033 µ m. In each light curve, we discarded values farther than 4 σ from the mean of a sliding window.  \nThe 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 semimajor 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.', 'Appendix A.4: HANSOLO': 'The HANSOLO (atmospHeric trANsmission SpectrOscopy anaLysis cOde) pipeline was originally developed to analyse groundbased 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.  \nWe used the LACOSMIC algorithm (van Dokkum 2001) to remove cosmic ray e ff ects from the two-dimensional images and identified the spectral trace by using a Mo ff at function fit to each column. The sky background was calculated on a column-bycolumn 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.  \nConsistent 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 σ 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 A.1. 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  \nwith 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 A.1.', 'Appendix A.5: transitspectroscopy': 'We take the corrected timeseries data from stark long-wave analysis and use an open-source tool transitspectroscopy (Espinoza 2022) 6 . 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 σ 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 µ m.  \nWe used juliet to fit the occultation model to the whitelight 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 bestfitted occultation depths are tabulated in Table A.1. We also determined relative occultation depth spectra using the procedure described in Appendix A.1.1 and plotted in Fig. A.1.', 'Appendix 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-theramp 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-byrow 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.  \nUsing 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).  \nTable A.2. Spectroscopic parameters for 55 Cnc.Notes. ( a ) Posteriors from the SED modelling.  \nFig. B.1. The observed stellar spectrum with NIRCam / JWST (in yellow) is shown with Crossfield (2012) empirical spectrum and a blackbody at 5269 K.  \n<!-- image -->', 'Appendix 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 B.1 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.', 'Appendix 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 Easy 7  \n(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 di ff erent photospheric parameters and iterating until all parameters converged. Throughout the modelling, we held the macro- and micro-turbulent velocities, V mac and V 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 A.6.  \nThe stellar radius was modelled with the SED fitting software astroARIADNE 8 (Vines & Jenkins 2022) using priors from SME and photometry from the Johnson B and V magnitudes (APASS), GG BP G RP (DR3), JHKS magnitudes ( 2MASS ), WISE W1-W2, and the Gaia DR3 parallax. We utilized three di ff erent 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 ± 0 . 011 R ⊙ . The luminosity is 0 . 63 ± 0 . 02 L ⊙ , and the visual extinction is consistent with zero (0 . 03 ± 0 . 03). We derived a stellar mass of 0 . 639 + 0 . 021 -0 . 020 M ⊙ 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 ± 0 . 010 R ⊙ based on interferometric measurements and the parallax from van Leeuwen (2007). Updating this calculation with the Gaia DR3 parallax, this radius becomes 0 . 962 ± 0 . 010 R ⊙ 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 / CO 2 and SiO / SiO 2 / 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 wl.  \nIt is also important to note that the depicted centre-log-ratio posterior ξ 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 ξ values must be zero.  \nFor 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.  \nFig. C.1. 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 ξ 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 ξ values must be zero.  \n<!-- image -->  \nTs (K)  \nTa (K)  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nArticle number, page 28 of 32  \n<!-- image -->  \nFig  \n<!-- image -->  \nFig. C.8. 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.  \n<!-- image -->  \nFig  \n<!-- image -->  \nFig  \n<!-- image -->  \nArticle number, page 32 of 32'}