Is the Atmosphere of the Ultra-hot Jupiter WASP-121 b Variable?

We present a comprehensive analysis of the Hubble Space Telescope observations of the atmosphere of WASP-121 b, an ultra-hot Jupiter. After reducing the transit, eclipse, and phase-curve observations with a uniform methodology and addressing the biases from instrument systematics, sophisticated atmospheric retrievals are used to extract robust constraints on the thermal structure, chemistry, and cloud properties of the atmosphere. Our analysis shows that the observations are consistent with a strong thermal inversion beginning at ∼104 Pa on the dayside, solar to subsolar metallicity Z (i.e., −0.77<log(Z)<0.05 ), and supersolar C/O ratio (i.e., 0.59 < C/O < 0.87). More importantly, utilizing the high signal-to-noise ratio and repeated observations of the planet, we identify the following unambiguous time-varying signals in the data: (i) a shift of the putative hotspot offset between the two phase curves and (ii) varying spectral signatures in the transits and eclipses. By simulating the global dynamics of WASP-121 b's atmosphere at high resolution, we show that the identified signals are consistent with quasiperiodic weather patterns, hence atmospheric variability, with signatures at the level probed by the observations (∼5% to ∼10%) that change on a timescale of ∼5 planet days; in the simulations, the weather patterns arise from the formation and movement of storms and fronts, causing hot (as well as cold) patches of atmosphere to deform, separate, and mix in time.

1. INTRODUCTION Spectroscopic measurements of transiting exoplanets have provided a wealth of "snapshot" information about the thermal structure, chemistry, and cloud properties of exoplanet atmospheres (e.g., Charbonneau et al. 2002;Tinetti et al. 2007;Kreidberg et al. 2014a,b;Swain et al. 2008b,a;Stevenson et al. 2014;Madhusudhan et al. 2014;Line et al. 2016;Sing et al. 2016;Tsiaras et al. 2018;Tsiaras et al. 2019;Benneke et al. 2019;Mansfield et al. 2021;Edwards et al. 2020;Skaf et al. 2020;Pluriel et al. 2020;Mugnai et al. 2021;Line et al. 2021;Changeat et al. 2022;Edwards et al. 2023; JWST Transiting Exoplanet Community Early Release Science Team et al. 2023).However, temporally varying information has yet to be unambiguously obtained by observations.This is partly because, prior to the recently launched James Webb Space Telescope (JWST), exoplanet atmospheres have generally been studied with a single observation whose spectral feature signal-to-noise ratio (S/N) is too low.In an attempt to reduce the noise, the current standard practice is to average the signal from different observations; however, the averaging removes any temporal variability that may be captured.On the other hand, when a single observation can achieve a high enough S/N, the observation of the planet is generally not repeated-due to observing time constraints.
With the Spitzer telescope, a number of studies have analyzed repeated measurements of individual transiting exoplanets via photometric multi-epoch measurements of eclipses.Many of these studies did not detect atmospheric variability below a certain level, due to the quality of the data (e.g., Agol et al. 2010;Crossfield et al. 2012;Ingalls et al. 2016;Morello et al. 2016;Kilpatrick et al. 2020;Murphy et al. 2022).Others, however, have suggested time-dependent shifts in phase-curve offsets for at least three exoplanets-HAT-P-7 b, WASP-12 b, and Kepler-76 b-using either the Kepler or Spitzer telescopes (Armstrong et al. 2016;Bell & Cowan 2018;Jackson et al. 2019;Wilson et al. 2021;Ouyang et al. 2023).The latter studies have speculated that such changes might be due to varying cloud structures, but conclusive interpretations of the datasets have remained elusive (Bell et al. 2019;Lally & Vanderburg 2022;Wong et al. 2022).Hence, presently there exists no unambiguous detection of atmospheric variability in the atmospheres of transiting exoplanets.
In contrast, atmospheric variability is commonly reported for non-transiting exoplanets which are characterized by high-contrast imaging (e.g., Artigau et al. 2009;Biller et al. 2015;Metchev et al. 2015;Biller 2017;Manjavacas et al. 2019;Vos et al. 2022).Among them, the ∼11-19 Jupiter-mass planet VHS 1256-1257 b, which has recently been observed by the JWST-NirSpec and JWST-MIRI instruments as part of the Early Release Science program (Miles et al. 2023), exhibits one of the largest amplitudes of observed atmospheric variability; for example, a periodic brightness change of up to 38% with a period of ∼15 hours is reported by Zhou et al. (2022).Many other planetary mass companions exhibit similar levels of variability.
Recently, intriguing possibility of time-variability for the ultra-hot Jupiter WASP-121 b has been reported in two studies, Wilson et al. (2021) and Ouyang et al. (2023).The former study compares spectra from a ground-based observation using Gemini-GMOS and a Hubble Space Telescope (HST) observation using HST-STIS and finds differences in the two spectra, which could be associated with a temporal variability.The latter study uses ground-based data from SOAR-GHTS and finds a spectra that also does not match that of the previous HST observations.These studies associate the observed differences with the presence of enhanced scattering slope in the case of GMOS, which could be explained by clouds or hazes, and varying abundances of molecular TiO/VO in the case of GHTS.Additionally, phase-curve observations with the HST (Mikal-Evans et al. 2022), Spitzer (Morello et al. 2023), and JWST-NirSpec G395H (Mikal-Evans et al. 2023) also report different phase-curve characteristics (i.e., "hotspot" offset and shape), which could be indicative of moving hot regions in the atmosphere.However, the variability inferred in these works again relies on combining the constraints from different instruments and/or observing conditions.
On the atmospheric dynamics modeling side, many hot Jupiter simulations in the past have suggested the presence of a single, stationary hot region eastward of the substellar point-particularly between the ∼10 4 Pa to ∼10 3 Pa pressure levels (e.g., Cooper & Showman 2006;Dobbs-Dixon et al. 2010;Parmentier et al. 2018;Komacek & Showman 2020).However, at high-resolution, highly-dynamic, variably-shaped (and often multiple) hot regions emerge instead (Skinner & Cho 2021;Cho et al. 2021;Skinner & Cho 2022;Skinner et al. 2022).In these simulations, a long-lived giant storm-pair forms near the substellar point, drifts initially toward one of the terminators, then rapidly translate westward thereafter-traversing the nightside and ultimately breaking up or dissipating near the eastern terminator; this cycle is quasi-periodic.Throughout each cycle, hot (as well as cold) patches of air are chaotically mixed over large areas and distances by the storms and sharp fronts around them.Similar mixing due to storms and fronts has initially been predicted in highresolution simulations (with a different initial condition than in the above studies) by Cho et al. (2003), who suggested that weather patterns would lead to potentially observable variability on hot exoplanets.
In this backdrop, we present here results from an indepth study of WASP-121 b atmosphere-focusing on its variability.WASP-121 b is one of the best targets for atmospheric characterization because it is characterized by a high S/N and has been observed multiple times.It has been observed four times with the HST Wide Field Camera 3 Grism 141 (WFC3-G141): one transit in June 2016, one eclipse in November 2016, and two phase-curves in March 2018 and February 2019.Significantly, we utilize the entirety of these observations in this work.Previous analyses from a combination of facilities-HST, TESS, Spitzer, JWST and groundbased facilities (Tsiaras et al. 2018;Evans et al. 2018;Mikal-Evans et al. 2020;Sing et al. 2019;Cabot et al. 2020;Ben-Yami et al. 2020;Hoeijmakers et al. 2020;Borsa et al. 2021;Daylan et al. 2021;Mikal-Evans et al. 2022;Changeat et al. 2022;Gibson et al. 2022;Azevedo Silva et al. 2022;Mikal-Evans et al. 2023)-have detected the presence of water vapour, absorbers of visible radiation (VO and TiO), hydrogen ions (H − ), and atomic species (Ba, Ca, Cr, Fe, H, K, Li, Mg, Na, V, and Sr); but, atmospheric variability has not yet been detected.
The outline of this paper is as follows.Our basic methodology and codes are presented in Section 2. The reduction procedure for the construction of a consistent set of spectra from the raw observations is described in Section 3.Then, the results from our state-of-the-art atmospheric retrievals, which use the newly-developed "1.5D phase-curve retrieval" models (Changeat & Al-Refaie 2020;Changeat et al. 2021) and one-dimensional (1D) models, are presented in Section 4 and Section 5, respectively; the 1.5D models enable mapping of the atmospheric properties (i.e.chemistry, cloud, and thermal structure) as a function of longitude using the entire phase-curve data, while the 1D models are used to extract information from each of the recovered spectrum individually.In Section 6, we present the findings from the highest resolution, three-dimensional (3D) dynamics simulations of WASP-121 b atmosphere to date.Finally, discussion and conclusion are presented in Section 7. Additional material is included in an Appendix as well.

METHODOLOGY
The WASP-121 b data we have analyzed in this work is obtained using HST WFC3-G141.This data is publicly available at the Mikulski Archive for Space Telescopes (MAST) 1 .Importantly, we have chosen to not include the observations from other telescopes or instruments, since the combination of multiple instruments is known to produce incompatible results (Yip et al. 2020(Yip et al. , 2021;;Edwards et al. 2023).The analyzed data is from one eclipse, one transit, and two phase-curves-comprising a total observing time of about 90 hours; each phasecurve observation contains two eclipses and one transit.For consistency, we have used the same data reduction pipeline, Iraclis, and have adopted identical assumptions for all of the observations.The individual eclipse and transit events which are not from the phasecurve observations have been previously extracted using Iraclis by Changeat et al. (2022) [hereafter C22], and the phase-curve data has been recently analyzed by Mikal-Evans et al. (2022) [hereafter ME22].However, the phase-curves have not been extracted using Iraclis.Therefore, we have re-analyzed the data with Iraclis in order to ensure consistency of treatment with that by C22.
Equipped with a consistently treated set of transit, eclipse, and phase-curve spectra, we use a suite of atmospheric retrieval codes and a high-resolution atmospheric dynamics model code which is extensively tested and validated specifically for hot exoplanet simulations.Our aim here is to perform a robust extraction of the thermal structures and chemical abundance profiles in WASP-121 b's atmosphere, which allows us to estimate planetary formation markers and investigate potential time variability.Broadly, our methodology can be grouped into four main activities, or parts: Part 1: extracting a consistent set of WFC3-G141 light curves for the transit, eclipse, and phase-curve data using Iraclis; fitting the light curves for the phase-curve data, using the PoP (Pipeline of Pipes) code (see description in Materials and Methods); and, testing various assumptions to model the instrument systematics in the literature.The two observations, while close, present differences that could come from either atmospheric variability or instrument systematics.Note that the second-order coefficients from the sine phase-curve models are fixed to the best-fit values from the combined fit (see Figure B3).
Part 3: analyzing the transit and eclipse data using 1D retrievals; and, incorporating the constraints (e.g., chemical parameters) obtained in Part 2 to reduce the degeneracies between temperature and chemistry so that observations can be analyzed individually rather than just in sum.
Part 4: performing high-resolution, global atmospheric dynamics simulations with the psudospectral code BoB (Built On Beowolf) (e.g., Skinner & Cho 2021;Polichtchouk et al. 2014), suitably optimized and forced with T -p profiles obtained in Part 2; and, interpreting the observed variability informed by these simulations.
A more detailed description of each part is provided in the Materials and Methods section of the Appendix.

Combined white light curve correction
Performing the extraction from the raw full-frame images of the two observed phase-curves with Iraclis, we obtain very similar results to ME22 (a comparison is provided in Figure A2).We then correct and fit the reduced light curves with PoP (see description in Materials and Methods).Here, the two observations are fitted together, sharing orbital (mid-transit time t mid , inclination i, and semi-major axis a) and model (planet-to-star radius ratios R p /R s , and phase-coefficients: C 0 , C 1 , C 2 , C 3 , C 4 ) parameters, as well as the parameters for the short-term HST systematics.The parameters for the long-term HST systematics are not shared to accommodate for the six different observation segments.Fitting the white light curves, we explore the effects of different short and long-term HST systematics on our results.
For the short-term ramps, we have attempted two models: a simple exponential (e.g., as in Tsiaras et al. (2016b) [hereafter 2-param IS short ]) and a double exponential (e.g., as in de Wit et al. (2018) and Mikal-Evans et al. (2022) [hereafter 4-param IS short ]).For the long-term systematics, three options are possible: linear, quadratic, and hybrid (e.g., quadratic for the first segments of each visit and linear for the others, as in ME22).The comparison of these runs can be found in Figures B1 and B2.Overall, we conclude that assuming simple or double exponential short-term ramp does not change our results for this dataset.For consistency with ME22, we therefore adopt the double exponential ramp model for the remainder of the study; for the long-term ramp, however, the corrected phase-curve observations depend on the assumption of linear, hybrid, or quadratic option.
Comparing the Bayesian evidence (E), a linear, hybrid, and quadratic correction give log-values, ln(E), of 5733, 5797, and 5803, respectively.While the quadratic case gives a slightly higher log-evidence, such an assumption is too flexible for the second segment (i.e., the data around transit), generating artificial nightside flux and causing large degeneracies (see posterior distribution in Figure B2).We therefore adopt the more conservative approach and adopt a hybrid long-term ramp, as was done in ME22.The final recovered white light curve is corrected from the instrument systematics and compared with the results of ME22 in Figure B3.The residuals, in particular, demonstrate good agreement with the previous literature results.

Individual white light curve corrections
To investigate the variability of the atmosphere and the instrument systematics in the available phase-curves observations, we have reproduced our white light curve analysis for each visit individually.Due to the lower amount of information, we choose to fix the orbital and second-order phase-curve parameters (C 3 and C 4 ) to the ones of the combined fit.For those runs, we have also employed the hybrid trend for the long-term ramps and the double exponential model for the detector shortterm systematics.Figure 1 shows the corrected white light curves for the individual fits; Figure B4 shows the corresponding posterior distributions.
The recovered phase-curves in Figure 1 exhibit clear differences.For instance, the first-order phase of the sinusoidal model is larger in the 2018 visit (0.19 rad) than in the 2019 visit (0.03 rad), or a variation in white transit depth is found.These differences could come from a combination of atmospheric variability (e.g., movement of hot/cold regions or changing cloud coverage)-a possibility we explore further in Section 6-and instrument systematic effects.In any case, these results suggest that combining HST phase-curve observations may not be straightforward.

Spectral light curve fitting
The spectral extraction is done using two different binning strategies, at "low" (i.e., as in ME22) and "medium" (i.e., as in C22) resolutions; see Materials and Methods for more details.Figures B5 and B6 show the corrected spectral light curves for the two cases, respectively.Inspecting the light curves, both strategies lead to similar residuals.We therefore moved on to the extraction of the spectra from the corrected light curves.Using 16 temporal bins (about 1.5h of observation), the final extracted spectra compared to that of ME22 are shown in Figure B7.Overall, the spectra agree very well-except for phase 0.07, where the reductions of our paper are consistent with each other but show larger flux than for phase 0.05 of ME22; this is despite the similar flux for phase 0.93, when compared with phase 0.95 of ME22.This slight difference could arise from the fact that ME22 includes the in-transit planetary flux (after removal of the transit signal) for bins 0.07 and 0.93 (those bins labeled 0.05 and 0.95 in ME22 span 3h), while we chose to ignore the in-transit planetary flux for consistency and simplicity.

PHASE-CURVE ATMOSPHERIC RETRIEVALS
One of the goals of this study is to robustly characterize the bulk properties of WASP-121 b atmosphere using the combined phase-curve data.As described in Materials and Methods, we use the 1.5D atmospheric retrievals to interpret the observed data.2This retrieval technique analyzes the two phase-curve observations using a single unified atmospheric model (e.g., a single likelihood); hence, it efficiently exploits all the information content-as opposed to the more traditional 1D retrieval performed individually for each phase (see, e.g., Stevenson et al. 2017).Since the "hotspot" offset D HS and (angular) size A HS are difficult to constrain from HST data alone (Changeat et al. 2021), we have tested different combinations and found that (D HS , A HS ) = (30 • , 50 • ) leads to the highest Bayesian evidence.Therefore, we here focus on this case.
Figure 2 (and the associated animation) shows the best-fit spectra and recovered thermal structure from 1.5D retrievals of the low-resolution data.Here, we obtain consistent T -p profiles by reproducing the retrievals on the "Medium" resolution spectra as well as those from ME22 (see Figures C1 and C2), demonstrating that this information is independent of the data reduction.The chemical parameters (i.e., metallicity and C/O ratio) are slightly dependent on the spectral resolution (especially in terms of precision, see Figure 3), which could be due to reduced correlation between the spectral channels at lower resolution, or information dilution occurring during the fit due to the lower S/N of the light curves at higher resolution.Since our "Low" resolution reduction is consistent with ME22, we focus the rest of our discussion on this case; nevertheless, full posterior distributions for all three retrievals are provided in Figure C3.
Importantly, due to the resolving power of phase-curve data, our retrievals allow precise thermal and chemical estimates at different locations in the planet's atmospheres to be obtained.Importantly, for ultra-hot Jupiters, the presence of absorption features for water, refractory species (TiO, VO, and FeH), and hydrogen ions in WFC3 allows to break the degeneracies between metallicity and C/O ratio (see also Changeat 2022).Given our retrieval assumptions, we find a strong thermal inversion on the dayside with the hottest region (here labeled hotspot in Figure 2) being ∼300 K hotter than the rest of the dayside between 10 5 Pa and 10 3 Pa (see red and orange profiles).
The thermal inversion could be caused by two different mechanisms.One mechanism is the production of energy at high altitudes by the presence of refractory molecules and H − (the latter from H 2 thermal dissocia-1000 2000 3000 4000 5000 Temperature (K) Recovered temperature-pressure (T -p) profiles (left) and best-fit spectra (right) for the phases from 0.05 (blue) to 0.5 (red), obtained from the phase-curve atmospheric retrieval.In the T -p plot, the shaded regions correspond to one and three sigma confidence regions (dark to light, respectively).The radiative contribution function is also shown in dashed line, colored for each region: hotspot (red), dayside (orange), and nightside (blue).These retrievals show good agreement with the observed data and demonstrate a strong dayside thermal inversion, with the presence of a hotter region (e.g.hotspot).The best-fit T − p profiles (solid lines, left) are used to thermally force the atmospheric dynamics simulations.This figure is accompanied by a 15 s video, available online at the journal, showing the evolution of WASP-121 b emission (from blue to red) and the corresponding thermal structure as a function of phase.As the planet moves from transit to eclipse, absorption features in the data are replaced by emission features.These spectral variations enable the characterization of the thermal structure and chemistry across WASP-121 b's atmosphere.tion): the temperature in the inversion region of the dayside is indeed hot enough to dissociate most moleculesincluding water and even more stable volatiles (CO and CO 2 ) and refractory molecules (FeH, TiO, and VO), along with H 2 (see the retrieved chemical profiles of Figure C4).At lower pressures, the atmosphere could partially be ionized, with an increased abundance of free electrons creating a continuum H − opacity, as suggested for other similar ultra-hot Jupiters (Edwards et al. 2020;Pluriel et al. 2020;Changeat & Edwards 2021;Changeat 2022).Another possible mechanism is heat deposition of breaking or saturating planetary and gravity waves launched from the atmospheric region below (e.g., Watkins & Cho 2010;Cho et al. 2015).Both mechanisms likely contribute to the observed thermal inversion layer.The retrievals we performed include a "gray" cloud model (i.e., constant opacity cloud deck); however, large cloud patches are not favored by the data (see Figure C3) despite the temperatures at the nightside being potentially suitable for silicate cloud formation (Powell et al. 2018;Gao et al. 2021).
Comparing the results from all the reductions (see also recovered T -p profiles and posteriors for other hotspot characteristics in Figures C5 and C6), we conclude that the data is consistent with a solar to slightly sub-solar metallicity Z and a super-solar C/O ratio.For the lowresolution case, we estimate the mean chemical characteristics of the planet to be log(Z) = −0.19+0.11  −0.13 and C/O = 0.80 +0.03 −0.05 .A more conservative estimate, encompassing the uncertainties from all the reductions and retrievals tested in this work is −0.77< log(Z) < 0.05 and 0.59 < C/O < 0.87.As a byproduct, this enables us to also speculate on the formation history of this planet; specifically, the obtained Z and the super-solar C/O ratio are suggestive of an early formation via significant gas accretion (i.e., without significant planetesimal pollution) and beyond the snowline of the proto-planetary disk (e.g., Öberg et al. 2011;Mordasini et al. 2016;Madhusudhan et al. 2017;Brewer et al. 2017;Cridland et al. 2019;Shibata et al. 2020;Turrini et al. 2021;Pacetti et al. 2022).
To further support the above conclusions, we have performed additional sensitivity tests, which are shown in Figure C7.We apply ±3σ departures, where σ is the retrieved uncertainty on the modified parameter, to the best-fit Z and C/O from the low-resolution fit and compare the simulated spectra with the observations.Introducing these departures leads to spectra that do not properly explain the observations, confirming the magnitude of the recovered uncertainties for the atmospheric chemistry of WASP-121 b.
Due to the high constraints on the mean atmospheric properties of this planet obtained via the phase-curve data, the parameters extracted at this stage (e.g., thermal profiles and chemistry) can serve as priors for the subsequent parts of our analysis.In particular, assuming that Z and C/O ratio remains spatially homogeneous and constant in time allows us to re-use the retrieved values to analyze the transits and eclipse data individually; the thermal profile, chemistry, and cloud properties extracted from individual transit and eclipse observations are known to be much more degenerate on their own.Additionally, the recovered thermal profiles provide important information for dynamics calculations.For example, they can be re-introduced as an observation-driven forcing to enable physically realistic and case-specific simulations.

TRANSIT AND ECLIPSE ATMOSPHERIC RETRIEVALS
As mentioned, C22 has previously reduced the individual transit and eclipse datasets with the Iraclis pipeline; therefore, we make use of the spectra from that work directly.Already interesting differences appear in the transit spectra, although the eclipse spectra look more alike (Figure 4).For both transits and eclipses, we perform 1D retrievals using the standard TauREx3 models.We have first attempted to retrieve all the free parameters of the models without particular priors; but, as expected for HST data, the degeneracies between thermal structure, chemistry, and cloud properties were difficult to break from individual HST transit/eclipse spectra: we could not extract a consistent picture.However, since Z and C/O ratio are expected to remain time-independent3 and have better constraints from the phase-curve data, we have decided to re-inject this information from the 1.5D retrievals.Therefore, the chemistry is fixed to the median value from the retrieval on the low-resolution spectra; this has allowed to obtain a consistent fit of the spectra from all the visits (Figure D1).For completeness, the posterior distributions are presented in Figures D2 and D3 for the transits and eclipses, respectively.
For the transits, the three individual observations show the presence of covering hazes or clouds.Specifically, the transit spectra captured during the two phasecurves (blue and green) are fully cloudy (i.e., consistent with featureless), while the observation obtained in 2016 Wavelength ( m)  6 and 7).For example, during the transits, the observed temporal variations could be interpreted as a formation of intermittent clouds and/or hazes; during the eclipse, the observed variations may be due to subtle changes in the thermal structure of the atmosphere-induced by motions of hot/cold regions from the planet's atmospheric dynamics.
(orange) shows clear spectral modulations from water vapor.Note that the red end of this spectrum, however, cannot be fit properly using the chemical equilibrium assumption; however, free chemistry retrievals, which use H 2 O, VO, and H − opacities beyond equilibrium, could achieve a better match.Nevertheless, Transit 1 (orange) is not consistent with a featureless spectrum, as shown by ∆ ln(E) = 18.3, and possesses strong water vapor absorption features.Transit 2 (blue) displays an interesting multi-modal solution.The spectrum is best explained by either very high temperatures (T ≈ 3000 K) at the terminator, forcing dissociation of the main molecules and leaving a flat contribution from the H − continuum, or a more cloudy/hazy atmosphere with a slight slope towards the blue end of the spectrum.Given the un-physical nature of the high-temperature solution, we suggest that this second observation is consistent with the presence of hazes, especially as a lower dimension featureless fit achieves a similar Bayesian evidence, ∆ ln(E) = 0.5.For Transit 3 (green), we find the observation consistent with a flat spectrum, which would be well explained by clouds and/or hazes, given ∆ ln(E) = −0.6.In transit, stellar activity (i.e., unocculted stellar spots and faculae) can can cause important spectral variations between repeated observations (Rackham et al. 2018;Thompson et al. 2023).However, long-term monitoring campaign from the ground (Delrez et al. 2016;Evans et al. 2018) suggests that WASP-121 is a very quiet and stable star.If confirmed, these results would indicate a transient formation of cloud/haze structures at the terminator of WASP-121 b.
We note that gray clouds, which were only considered for the night-side, were not recovered by the retrievals on the phase-curve data.Many reasons could explain this result: i) the emitted signal in phase-curve does not probe the same altitudes as the transit data, ii) the observed clouds are poorly described by the gray cloud model, or iii) the clouds are located around the terminator region and not covering large patches of WASP-121 b.Additionally, unambiguously determining the presence of gray clouds from emission data is more challenging due to the degeneracies with the vertical temperature profile and the shorter geometrical path length through the atmosphere.
For the eclipse data, the spectral differences are much smaller and difficult to infer by visual inspection of the spectra.We have conducted atmospheric retrievals and extracted the thermal structure for the five different eclipses individually (see left panel of Figure 5).The thermal profiles are overall consistent across the five observations (i.e., similar thermally inverted structure), but we find variations in the mean temperature of 310 K when averaging the profiles over the 10 5 Pa to 10 3 Pa region.Importantly, this range is much larger than the average one-sigma uncertainty of the profiles in the same region (the averaged standard deviation is 108 K).For instance, in Figure 5, the thermal profiles extracted from Eclipse 2 (blue) and Eclipse 4 (red) are not consistent within the retrieved uncertainties.Similar to the phasecurve data, the observed differences in eclipse could be attributed to temporal variations of hot/cold regions in the planet's dayside and/or a changing thermal structure of the substellar point.While instrument systematics could remain (see Section above), the observed differences in hot region offset from the phase-curves, cloud coverage from the transits, as well as dayside thermal structure from the eclipses are all plausibly explained by the presence of atmospheric temporal variations.The observed differences are in fact expected from a theoretical atmospheric dynamics viewpoint due to the intense stellar heating contrast from the planet's parent star WASP-121.To investigate more precisely the possible origins of atmospheric temporal variations on the planet WASP-121 b and verify if they can affect our data to observable levels, we model its atmosphere with high-resolution dynamics simulations, to which we now turn.

DYNAMICS MODELING
We simulate the dynamics of WASP-121 b atmosphere with the BoB code at "T682L50" resolution, where T = 682 is the triangular truncation wavenumber (i.e., number of total and zonal modes each in the spherical harmonics) and L = 50 is the number of vertical layers (uniformly space in p); see Sec A.5, as well as Skinner & Cho (2021) and Polichtchouk et al. (2014), for detailed descriptions of the numerical model and simulation parameters.The use of BoB at this resolution-to directly guide the retrieval interpretation with numerical robust-ness and verisimilitude-is a significant feature of this study.The simulations are performed to obtain a broad idea and insights into the variability plausible on planets like WASP-121 b, when the flow is adequately resolved: it has recently been shown explicitly that hot-exoplanet simulations are not converged if the resolution employed is much below that in this work (Skinner & Cho 2021).As in most past studies, the atmosphere initially at rest is set in motion via a thermal relaxation to T -p profiles.
Here the forcing profiles are obtained from the retrievals described above (e.g., Fig. 2) and prescribed.Profiles at many different times, which have deviated from the prescribed one due to the nonlinear atmospheric motion, are compiled in Fig. 5 (right panel): they should be compared with observations (left panel).
Figure 6 shows temperature maps at p = 10 5 Pa from three widely separated times in the simulation.The maps demonstrate the highly variable nature of the planet's temperature field at a pressure level from which observable flux would originate.Snapshots of the temperature at two different pressure levels, p = 5×10 3 Pa and p = 10 5 Pa, show the vertical distribution, which is strongly barotropic-i.e,vertically aligned (Figure E1); a full movie of this simulation is provided in the Supplementary Materials.We also show chemical species maps, which are simply post-processed using the instantaneous density and temperature distributions from the simulation (Figures E2-E4) at p = 1 bar = 10 5 Pa and p = 50 mbar = 5×10 3 Pa (left two columns and right two columns, respectively), at t = 49 days and t = 62 days (left and right columns, respectively, for each p-level); distributions of the main relevant molecules (H 2 , H, H − , e − , H 2 O, CO, CO 2 , CH 4 , TiO, VO, and FeH) are shown.Not surprisingly, variations (vertical, horizontal, and temporal) induced by WASP-121 b's atmospheric dynamics impacts the chemistry as well as temperature, the latter of which is discussed more in detail below.Note that here a simple chemical equilibrium assumption is made, which may not be valid everywhere in the modeled domain-particularly if the reaction time is comparable to the advection time.
As can be seen in the figures, the complex motion of the atmosphere-and the organized structures thereincause hot and cold regions to chaotically mix in time.Generally, the hottest region periodically forms slightly eastward of the substellar point, but it always moves away from the point of emergence.Depending on when the atmosphere is observed, the hottest region can even be located on the west side of the substellar pointeither sequestered in a long-lived storm or generated near hyperbolic flow points between the storms (middle and left panels in Figure E1, respectively).
Interestingly, the dynamical behavior here is reminiscent of that of WASP-96 b in the p ≳ 104 Pa vertical region.Skinner et al. (2022) have recently reported a quasi-periodic generation of giant cyclonic storms moving away westward from the point of emergence.This is due to "deep heating" (i.e.strong heating at ∼10 5 Pa), which may be experienced by some hot Jupiters.Here the similarity in behavior is likely due to the morphologically similar T -p profiles in the aforementioned region on WASP-121 b and WASP-96 b.The strength of the dayside-nightside difference is greater on WASP-121 b, but the profiles drive-and steer-the dynamics in a qualitatively similar way to WASP-96 b.
The similar dynamical behavior also leads to qualitatively similar disk-averaged flux signatures at p = 10 5 Pa, for example; cf. Figure 7 with Figure 4 in Skinner et al. (2022) 4 .In both figures, the fluxes exhibit quasiperiodic variations, with excess flux at the substellar ("Day side" in Figure 7) and east terminator regions.Hence, the flux is persistently "shifted to the east of the substellar point"-when averaged over the disk.It is important to note that, when not averaged, the hottest regions are rarely located near the equator and often situated westward and at higher latitude of the substellar point (middle and left panels in Figure 6, respectively); the regions are not always vertically aligned either (cf.temperature maps from two pressure levels at Day 49 in Figure E2).Movies of the simulation show both clearly.The movies also show that the timescale of the variability is ∼5 planet days with sharp flux changes occurring in much shorter (≲ 0.5 day) windows-as was reported for WASP-96 b-like atmospheres by Skinner et al. (2022).The above behavior overall may also explain why most infrared observations to date-except by Dang et al. (2018), Bell et al. (2019), and Morello et al. (2023)have reported only "eastward-shifted hotspots".
As expected, the magnitude of the variability in the model domain varies with the p-level: it is greatest at p ∼ 10 5 Pa.Above p ≈ 2×10 4 Pa, the dayside maintains a strong thermal inversion associated with an atmosphere that nevertheless remains highly variable.Given our model assumptions, this altitude-dependent behavior originates from the greater intensity of stellar irradiation on WASP-121 b compared to that of a typical hot Jupiter, leading to a much shorter thermal relaxation time, as well as the presence of the visible light absorbers H − , arising from the dissociation of H 2 (Figure E2).On WASP-121 b, while the stellar irradiation can still penetrate as deep as 10 5 Pa (as in a typical hot Jupiter), the upper layers of the atmosphere likely maintain much stronger dayside-nightside temperature gradient overall because of the shorter relaxation time (which is ∝ T −4 : Andrews et al. 1987;Salby 1996;Cho et al. 2008).The thermal profile at the substellar point shows a short period (∼5 days) temporal variation of the order of ±200 K, which could be captured by the HST observations.
To round out our investigation, we compare the temperature profiles and fluxes obtained in our dynamics simulations with the information obtained from the data.In Figure 5, we have already shown the agreement between the extracted temperature profiles from the observed eclipses (left panel) and the typical profiles resulting from simulations at the substellar point (right panel).Significantly, the predicted and observed spread of the temperature profiles at the substellar point in time is very similar, suggesting a qualitative agreement between observation and theory.From the simulations, we find that the three-σ altitude-averaged temperature range is 680 K between p ∈ [10 5 , 10 3 ] Pa (e.g., assuming a normally distributed temporal spread of the T -p profiles, this value means that altitude averaged T -p profiles 340 K hotter or cooler than the mean should be  considered outliers).In comparison, we find from the retrievals that the five eclipses have an average temperature range of 311 K, when averaged over the same pressure domain range.Despite the possibility of our data being affected by small HST systematics still remaining and the difficulty of directly comparing inherently different models (i.e, 3D vs 1D), the scale of the temperature variations from the eclipse observations and from the theory is compatible.2023) respectively.In addition, we present a comparison of the modeled spectra in Figure 8, which shows that variability of up to 10% in the observed flux is compatible with the simulations.Such a level of flux variability agrees with theoretical works on other objects (Skinner et al. 2022), and is within the capabilities of HST-provided that its instrument systematics are kept under control.Additionally, this level of variability should easily be captured by repeated observations of similar planets from JWST and Ariel.
Note that we obtain a lower level of variability when performing full radiative transfer calculation for the presented frames in Figure 8.This could be because of a selection bias (with the chosen frames) and/or because of the baroclinicity (vertical non-alignment) of the atmosphere, which smooth the variability when the flux is integrated over a large pressure range.Moreover, we find that the variability is wavelength dependent, which we believe is caused by the changes in the chemistry.In particular, in the case of WASP-121 b, H 2 thermally dissociates faster than H 2 O; hence, we expect larger variability signals for wavelengths between 1.0 µm and 1.3 µm.

DISCUSSION AND CONCLUSION
The extreme atmospheric conditions of WASP-121 b make it an ideal laboratory to test our understanding of physical and chemical processes in the atmospheres of exoplanets.Until now, detecting and studying weather patterns on exoplanets have remained elusive because of the lack of either adequate S/N or repeated, crossverifiable observations that can be directly compared.However, multiple, comparable observations with the HST for WASP-121 b exist-and now permit progress to be made.Moreover, truly high-resolution, numerically converged simulations also permit those observations to be interpreted with some fidelity, since the governing equations are now much more accurately solved than have been in the past.The uniform treatment in the analysis of observation, together with informed guidance from numerically accurate and validated simulations, are the salient features of this study: without at least both, variability cannot be confidently assessed at present.Here, bolstered by state-of-the-art simulations, we identify spectroscopic variability in the HST observations of WASP-121 b.
While some caution must still be exercised in interpreting HST data, given the well-known high level of systematics and the large number of assumptions required in reducing spectroscopic observations, we demonstrate here a strong potential evidence for variability associated with weather on WASP-121 b.Here the weather Figure 8. Spectroscopic eclipse flux obtained at three different times, t = {49, 62, 137} days, when the dynamics calculations are post-processed (top), and spectroscopic differences between two times, t = 62 days and t = 49 days.In the top panel, the eclipse spectra obtained for Observation 3 (green) and Observation 4 (purple) are also shown for reference.The simulation shows that the atmospheric variability should be wavelength dependent since, e.g., H − is highly sensitive to thermal dissociation (impacting the short wavelengths) while H2O remains more chemically stable.
is inferred from the following: i) the movement of the peak emission in two phase-curves, ii) the changing depth of the water feature in three transits, and iii) the variable retrieved thermal profiles in five eclipses.On WASP-121 b, the large dayside-nightside temperature gradient-which is not necessarily fixed in space and time5 -is expected to power (as well as steer) dynamical, thermal, and chemical processes.These include vortex instability, gravity wave and front generation, thermal dissociation, chemistry changes, and potential cloud/haze formation (e.g., silicate clouds), whose consequences are observable.Other studies using groundbased data have also suggested variable atmospheric conditions on WASP-121 b (Wilson et al. 2021;Ouyang et al. 2023).
A new finding in this study is that high-resolution dynamics simulations forced by T (p) information retrieved from observations show that ultra-hot Jupiters, such as WASP-121 b, likely have hot regions that are generally situated slightly eastward of the substellar point when disk-averaged-but whose actual shape and location are markedly variable in time.The variability is particularly visible in the modeled region, p ∈ [10 3 , 10 5 ] Pa, for which we have given the lion's share of focus in this study.This is in stark contrast with nearly all past hot Jupiter simulations, which show a fixed location and shape for a singular hot region; well-resolved simulations consistently indicate otherwise (e.g.Cho et al. 2003;Skinner & Cho 2021).More specifically, our simulations indicate that variability for WASP-121 b should be ∼5% to ∼10% in the disk-averaged flux with a frequency of ∼5 days.Hence, the signatures generated by these quasi-periodic weather patterns should be detectable with current as well as future instruments-e.g., by the JWST (Greene et al. 2016) and Ariel (Tinetti et al. 2021) telescopes-if repeated, high-quality observations are obtained.For the reduction and extraction of the spatially scanned spectroscopic images, we use the dedicated and publicly available pipeline Iraclis (Tsiaras et al. 2016b(Tsiaras et al. ,c, 2018)).The individual transit and eclipse events have already been extracted in the population study of C22, so we obtain those outputs from this study.For the phase-curves, however, the data has not previously been extracted with Iraclis, so we perform all the extraction steps described in Tsiaras et al. (2019) as implemented in C22.Those steps consisted in zero-read subtraction, reference-pixels correction, non-linearity correction, dark current subtraction, gain conversion, sky background subtraction, calibration, flat-field correction, and bad-pixels/cosmic-rays correction.We then use Iraclis to extract the white (1.088-1.68µm) and spectral light curves from the reduced images, taking into account the geometric distortions caused by the design of the WFC3 tilted detector.The two observed phase-curves had re-acquisition events, separating the observations in three distinct segments (see Figure A1).Re-acquisitions cause displacements of the images onto the detector and induce larger systematics at the beginning of each event.Pixels with higher flux rates are affected more, causing wavelength-dependent long-term ramps that require an additional treatment.

A.1. Data reduction with Iraclis
The extracted white light curves for both visits are compared to the ones obtained by ME22 in Figure A2.While completely independent, both reductions lead to very similar results.To extract the spectral light curves and study the impact of spectral binning on our conclusions, we explore two strategies: 1) A medium spectral binning.To ensure consistency, we test the same binning as other Iraclis reductions (Tsiaras et al. 2018;Changeat et al. 2022;Edwards et al. 2023), adopting an HST template composed of 18 wavelength bins (see: Tsiaras et al. 2019).
2) The low spectral binning.For better comparisons with ME22, we test their proposed 12 bins strategy.2022a).While TauREx was originally designed for atmospheric retrievals (Waldmann et al. 2015b,a), its last version provides generic classes that can easily be used to solve optimization problems outside of its original scope.
In terms of code structure, PoP combines an observation pipeline with a scientific model pipeline.Here, pipelines refer to a series of units representing individual transformation steps.In the context of this study, the "observation" pipeline has a unit to load the Iraclis reduced light curves and a sequence of other units to correct for the HST instrument systematics.The "model" pipeline contains an idealized light curve model.Note that the "model" pipeline can also possess additional transformation steps to convert the model's outputs to the "observation" pipeline format.During optimization, the free parameters are marginalized over likelihood for both the observation and the model pipelines.This structure makes rapid modifications of the pipeline easy and creates a clear separation between astrophysical models and instrument models.
Model Pipeline: In the case of WASP-121 b, we model the transit and eclipse events using the Pylightcurve package (Tsiaras et al. 2016a) as done in past works employing Iraclis, such as C22.The transit light curve drop (LC T ) and the normalized eclipse light curve drop (LC E ) from Pylightcurve are combined with a description of phase-curve variations (LC S ).Following the standard practice in the literature for HST (see e.g., Mikal-Evans et al. 2022)  ations using a first or second-order sinusoidal function: (1) where Φ is the orbital phase, C 0 is the minimum of the nightside flux, C 1 is the maximum of the dayside flux, and C 2 is the phase-curve offset.C 3 and C 4 correspond to the second-order terms of the sine phase variations.As in Dang et al. (2018), the full phase-curve model M P C is constructed by: In transit, we computed the limb-darkening coefficients with the Claret (2000) law using the stellar ATLAS models (Kurucz 1970;Howarth 2011;Espinoza & Jordán 2015) included as part of the ExoTETHyS package (Morello et al. 2020a;Morello et al. 2020b).Those coefficients were fixed during the light curve fits.
Observation Pipeline: With HST, the instrument systematics typically consists of a long-term trend affecting each visit and short-term ramps affecting each orbit.For the long-term trend (IS long ), the standard practice is to assume a linear or a quadratic behavior (Tsiaras et al. 2016b;Kreidberg et al. 2018;Arcangeli et al. 2019).In the case of WASP-121 b, since re-acquisition events had occurred, we were required to fit the long-term trends separately for each segment of the observations (here segments are noted with the index i with i ∈ {1, 2, 3}).The corresponding long-term trend is: where A i 0 is the linear coefficient of the segment i, A i 1 the quadratic coefficient of the segment i, and N i is a normalization factor for the segment i.The time t s refers to the time at the beginning of each observation segment.For the detector short-term orbital ramps (IS short ), previous studies have modeled its behavior using a combination of exponentials.Most studies discard the first frame of each orbit due to the larger systematics.Additionally, most methods also discard the first orbit, which is usually affected by a much larger and distinct exponential ramp.When discarding the first orbit and the first frame of each orbit, a simple exponential common to all the visits can be used (Kreidberg et al. 2014a;Stevenson et al. 2014;Tsiaras et al. 2016b): where B 0 and B 1 are the exponential factors, common to all visits in each observation and t i o is the time at the beginning of each orbit j.More recently, an evolution of this approach for the short-term ramps of HST has emerged (Zhou et al. 2017;de Wit et al. 2018).Their physically motivated model combines two exponential functions and is able to recover the first orbit.It is given by: where R is the correcting function mainly acting on the first orbit and given by: B 2 and B 3 are the exponential factors of R and t v is the time at the beginning of the visit.
In this paper, we explore the impact of various instrument systematics assumptions on the recovered data.
The approach, however, is the same in all cases and followed C22 for consistency.We start by fitting the white light curve with the goal to infer the free parameters of the system that are wavelength independenti.e., the mid-transit time (t mid ), the inclination (i), and the semi-major-axis (a).The other orbital parameters are fixed to literature values (Bourrier et al. 2019), as they have better constraints from complementary observations.The fit is conducted with the nested sampling algorithm MultiNest (Feroz et al. 2009;Buchner et al. 2014), using 1024 live points and an evidence tolerance of 0.5.An initial estimate of the observational noise is made by Iraclis (see Tsiaras et al. 2016b).However, we account for additional systematic noise or other unaccounted effects by re-scaling the uncertainties according to the RMS of the residuals, and perform a second fit using the updated uncertainties.In practice, this step has little effect for those WASP-121 b observations as the RMS of the residuals was very close to 1 at the first stage.Since nested sampling computes accurate Bayesian log-evidence, ln(E), we note that the Bayes Factor can be used to perform model selection.
To fit the spectral light curves, we employ a dividewhite strategy similar to Kreidberg et al.  2022) and we divide each spectral light curve by the corresponding white light curve.Since the detector ramps are correlated between wavelengths, this step essentially removes the short-term ramps and reduces the baseline drift in the data.As such, the observation pipeline for the spectral light curve only contains the long-term trend correction IS long (i.e., IS short is not needed).To match the divide-white observations, the model pipeline is also modified with an additional step.This step normalizes the modeled spectral light curves, dividing by the median white light curve obtained during the prior fit.Each spectral light curve is then fitted separately with this model, keeping t mid , i and a fixed to the best-fit value from the white.As with the white light curves, the errors are re-scaled to match the RMS of the residuals.
Finally, we extract the phase-curve spectra from the binned corrected light curves using 16 different phase bins (Φ) of equal dimensions 0.05, giving us: Φ ∈ {0.07, 0.12, 0.17, 0.23, 0.28, 0.32, 0.38, 0.43, 0.57, 0.62, 0.68, 0.72, 0.78, 0.82, 0.88, 0.93}.Except for Φ ∈ {0.07, 0.93}, this matches the bins adopted in ME22.For those bins the planetary flux around Φ = 0.0 is not included as it is blended during the transit, and a consistent bin size of 0.05 is used instead of 0.1 in ME22.Additionally, for this binning step, the finite integration time is not corrected for as the expected photometric distortions are below 5 ppm in the 0.05 intervals (Morello et al. 2022).For the transit spectra, self-blend (i.e, contamination by planetary emission) is not corrected for as the effect only affects the transit data by a few ppm in the case of WASP-121 b (Morello et al. 2019;Mikal-Evans et al. 2022;Morello et al. 2021).

A.3. Phase-curve retrievals
We analyze the information contained in the phasecurve data using atmospheric retrievals.We employ the TauREx 3.1 framework (Al-Refaie et al. 2021, 2022a), following the previously established methodology detailed in Changeat & Al-Refaie (2020); Changeat et al. (2021); Changeat (2022).We use the 1.5D phase-curve model, which is specifically designed to handle this type of observation, to simultaneously fit all the spectra in a Bayesian retrieval framework.The 1.5D model is composed of three different regions, referred to as hotspot, dayside, and nightside.Each region possesses independent properties allowing us to resolve large-scale atmospheric features from the data.The contribution of each region to the emitted flux at each phase is computed using a quadrature integration scheme (Changeat & Al-Refaie 2020).For this study, the structure of the planet is defined by 90 layers equally spaced in log space from p ∈ [10 6 , 10 −1 ] Pa.To first-order, such a model accounts for the main one-dimensional biases discussed in Feng et al. ( 2016); Taylor et al. (2020); Changeat et al. (2021).
As described in Changeat et al. (2021), the hotspot region is parameterized by two free parameters (hotspot size, A HS , and hotspot location, D HS ).We have first attempted to retrieve those parameters but due to the large degeneracies between those parameters, the thermal structure, and the chemistry, this led to un-physical solutions (a similar behaviour was found in Changeat et al. 2021).Therefore, we have decided to fix those parameters and instead explored fixed values.For the hotspot size, we test A HS ∈ {30, 50} degree cases.For the hotspot location, we test D HS ∈ {10, 20, 30, 40} degree eastward-shifted cases (those choices are also motivated by the results in ME22).
To parameterize the thermal structure, the temperature-pressure (T -p) profiles are created by linearly interpolating between T -p nodes (the pros and cons of such description are discussed in detail in the Appendix D of Changeat et al. 2021).For the hotspot region and dayside, we use seven nodes at fixed pressures (p ∈ {10 6 , 10 5 , 10 4 , 10 3 , 100, 10, 0.1} Pa), while for the nightside, since the information content is reduced due to the lower planetary emission, we choose to only use five nodes (p ∈ {10 6 , 10 5 , 10 3 , 10, 0.1} Pa).For the chemistry, we employ the GGChem (Woitke et al. 2018) chemical equilibrium code in its TauREx 3.1 plugin and we couple the two only free parameters (metallicity and C/O ratio) between the three different regions.Previous works have shown the importance of dis-equilibrium processes for hot Jupiters (Moses et al. 2011;Drummond et al. 2020;Venot et al. 2020Venot et al. , 2012;;Al-Refaie et al. 2022b), however, given the temperatures of WASP-121 b, chemical reactions should be fast, favoring chemical equilibrium for the investigated species (Parmentier et al. 2018;Kitzmann et al. 2018) (Abel et al. 2011(Abel et al. , 2012;;Fletcher et al. 2018) and Rayleigh scattering (Cox 2015).A fully opaque cloud deck (referred here as gray clouds) was also included on the night side of the planet but we were not able to find evidence for clouds from this phase-curve data.The parameter space of the model is explored using the nested sampling optimizer MultiNest (Feroz et al. 2009;Buchner et al. 2014), with 500 live points and an evidence tolerance of 0.5.To explore the parameter space, priors were chosen to be non-informative (i.e., uniform priors).Specifically, for the temperature points, we explored the space from T ∈ [300, 6000] K.For the chemistry, the metallicity was explored in log space from Z ∈ [0.1, 100] times solar, while the C/O ratio was explored from C/O ∈ [0.1, 2].From this retrieval analysis, we were able to extract averaged chemical properties of WASP-121 b as well as the thermal structure of the three considered regions.

A.4. Transit and eclipse retrievals
To evaluate the variability of WASP-121 b's atmosphere, we also analyze each eclipse and transit spectra individually, using 1D atmospheric retrievals with Tau-REx 3.1.Previous studies have shown the difficulty of extracting reliable constraints from single HST visits (Changeat et al. 2020) due to degeneracies between chemical abundances and thermal properties.To break those degeneracies, we use our most accurate estimate of the time-independent parameters from our phasecurve retrieval as priors for our individual retrievals.Using the low-resolution results, the metallicity Z of the atmosphere is therefore fixed at log(Z) = -0.19,while the carbon-to-oxygen ratio C/O is fixed at C/O = 0.80.Note that this simplification is not expected to always be correct, for instance cloud condensation can locally (and temporally) change the C/O ratio of the gas phase, which we do not model here (i.e., GGChem is used without condensation).However, without additional knowledge or constraints on condensates in WASP-121 b, this remains a reasonable and necessary assumption.
Transits: Since HST transits probe a narrow pressure range, and because it is mainly affected by the atmospheric scale height (Rocchetto et al. 2016), we consider a simple isothermal profile with a unique free parameter T for those observations (spectra shown in Figure D1).As with the phase-curve data, the temperature is explored with uniform priors (T ∈ [300, 6000] K).However, transit is much more sensitive to clouds than eclipse or phase-curve, so we have used a more complex representation of clouds from Lee et al. (2013), for which particle size and mixing ratio were fitted.This cloud model was favored by the Bayesian evidence compared to the gray-cloud case.
Eclipses: As planetary radius is known to be degenerate with temperatures in HST eclipses (Edwards et al. 2020;Pluriel et al. 2020), we fix this parameter to the literature value.For each spectrum (see Figure D1), we retrieve a thermal profile with a similar parameterization to the phase-curve case.Namely, the T -p profile is parameterized by linearly interpolating between seven freely moving T -p nodes.The pressure of each node is fixed to log-spaced values of pressures (i.e., p ∈ {10 6 , 10 5 , 10 4 , 10 3 , 100, 10, 0.1} Pa) and we retrieve the temperature of each point individually using the same uniform non-informative priors.We refer to Changeat et al. (2021); Rowland et al. (2023) for a more complete discussion on thermal structure parameterizations and their trade-offs.A simpler 3-point thermal profile (where the pressure levels of each node are left free) was also tested, which did not change our overall conclusions.For the eclipse spectra, we also decided to run a simple blackbody planet fit, which served as our comparison baseline.Following this procedure, and because of the additional chemistry priors from our phasecurve analysis, we have obtained well-defined thermal structures for each of the five eclipses.

A.5. Dynamics modeling
We model the atmospheric dynamics of WASP-121 b using the pseudospectral dynamical core, BoB (e.g., Rivier et al. 2002;Scott et al. 2004;Polichtchouk et al. 2014;Skinner & Cho 2021;Skinner & Cho 2022).The core has been outfitted and set up in Skinner & Cho (2022), Skinner &Cho (2021), andCho et al. (2021) especially for high-resolution hot Jupiter simulations; and, we refer the reader to those works for a more complete description of the code and governing equations solved.However, for the readers' convenience, we provide a brief summary of the key features of the model here -especially as they pertain to modeling of WASP-121 b atmosphere.Other numerical models have been used to study hot Jupiter atmospheric dynamics in the past (e.g.Showman & Guillot 2002;Cho et al. 2003;Cho et al. 2008;Cho 2008;Dobbs-Dixon et al. 2010;Thrastarson & Cho 2010;Polichtchouk et al. 2014;Heng et al. 2014;Mayne et al. 2014;Mendonça et al. 2018;Parmentier et al. 2018, and references therein), and we direct the reader to consult those works for instructive context.
BoB calculates the large-scale dynamics of the atmosphere by numerically solving the traditional and hydrostatic primitive equations in (longitude, latitude, pressure) = (λ, ϕ, p) coordinates an in vorticity-divergencepotential-temperature formulation.In the vertical (p) direction, BoB employs a second-order finite differencing scheme with free-slip boundary conditions at the top and bottom pressure surfaces.The present study builds upon extensive testing and validation of BoB, including simulations at high-resolution and under numerically challenging conditions resembling those found on WASP-121b (Polichtchouk & Cho 2012;Polichtchouk et al. 2014;Cho et al. 2015;Skinner & Cho 2021).
For the physical setup of WASP-121 b we use the parameters in Delrez et al. (2016).To simulate irradiation from the planet's host star we implement an "idealized" thermal forcing using the Newtonian relaxation scheme which accelerates the initially resting atmosphere (u = 0) towards specified hotspot and nightside equilibrium T -p profiles that are obtained from the phase-curve retrievals of the planet (see Fig. 2).The initial temperature is the average of the hot spot and nightside equilibrium temperatures.Due to the planet's close proximity to its star, the effect of its spherical geometry on stellar irradiation deposition is accounted for using a cosine profile to graduate the hot spot temperature from the substellar point to the terminators.The radiative cooling time τ r (p) is computed from the initial temperature profile following Cho et al. (2008); τ r (p) is approximately linear in log(p), ranging from O(10 6 )s at p = 10 5 Pa to O(10 2 )s at p = 10 3 Pa.
For the numerical setup, we use high horizontal and high vertical resolutions: T682 and L50 respectively.In the former, 'T' denotes the highest wavenumber retained in the spherical harmonic expansion (the truncation wavenumber) is n t = 682 and the latter 'L' denotes the number of vertical levels which are distributed lin-early in p-coordinate.By "high-resolution" we mean our simulations are above the minimum resolution required for numerical convergence of flow solutions on hot, tidally synchronized exoplanets (Skinner & Cho 2021).Note that high vertical resolution is also necessary to ensure the baroclinic region of the dayside and hotspot temperature profiles (p ≲ 10 4 Pa) is well represented in p.A small timestep size of ∆t = 6 seconds is used concomitantly with the fine grid spacing to ensure flows at the maximum sound speed (i.e., c s ∼ 4880 m/s in the hottest region of the atmosphere) are well captured.Hence, the simulations maintain a Courant-Friedrichs-Lewy (CFL: Courant et al. 1928) condition of well below unity.Simulations are time-integrated for 200 planet days, to ensure they reach a quasi-stable state long after their initial ramp-up period of ∼40 WASP-121 b days.
For the domain size, we model the expected p range (from p t = 10 3 Pa to p b = 10 5 Pa) from which radiation originates, as indicated by HST data (see Figure 2 and Figure 5).We have verified that the model equations are valid for this region by confirming that the retrieved thermal forcing profiles are stably stratified (i.e. that the Brunt-Väisälä frequency N = (g/ρ) (dρ/dz) is real).This is shown in Fig. C2.Below this (p > 10 5 ), the nightside temperature profile exhibits a jump in stratification.Although this could be due to increased uncertainty on the retrieved thermal structure in this region.
Finally, for numerical dissipation we use a high-order hyper-viscosity ∇ 16 with a small corresponding artificial viscosity coefficient of ν 16 = 10 48 .This prevents excessive kinetic energy removal from small-scale flows and hence ensures the dynamics of large-scale flows are well represented (see e.g.Cho & Polvani 1996;Skinner & Cho 2021).In addition, a very weak Robert-Asselin filter with coefficient ϵ = 0.02 is used to filter the additional computational mode arising from the models leapfrog time integration scheme (Asselin 1972;Thrastarson & Cho 2010).Besides this, no other drags are applied as these can coerce the flows to dynamically unphysical states (Polichtchouk et al. 2014).The simulations are allowed to evolve freely under thermal forcing from the retrieved temperature profiles.
For the post-processing of the three-dimensional atmospheric dynamics simulations, we employ two different approaches for evaluating the evolution of WASP-121,b from an observational perspective.The first is a one-dimensional time-series analysis of flux emitted by a single layer at the mid-region of our computational domain (0.1 bar) in order to evaluate the qualitative behavior of the planet's weather over hundreds of planet days.Here the black-body emission is calculated from BoB temperature maps, centered on key regions of interest (i.e., the substellar point, eastern terminator, western terminator, and antistellar point).While this approach does not account for the entire domains contribution towards atmospheric variability, it enables the key dynamical processes in the atmosphere to be isolated and studied in detail.
For second phase of our post processing, we link the outputs of BoB with the TauREx3 library to simulate observables and produce detailed 3D chemical maps of the planet's atmosphere.Due to the large size of the BoB simulation grid cube we isolate frames from the BoB calculations with significant spatial temperature differences that are likely observable for this analysis (e.g., t = 50 and t = 62 days).First, we derive the chemical maps over the entire (λ, ϕ, p) = (2048, 1024, 50) (i.e., grid cube) by performing calculations with the GGChem (Woitke et al. 2018) chemical equilibrium code using the T -p values from each grid square of the BoB simulations.We maintain fixed values for the metallicity and C/O ratio based on the median values obtained during our atmospheric retrieval exploration (log(Z) = -0.19 and C/O = 0.80).The flux emitted by the planet is then computed using the TauREx3 plane-parallel radiative transfer model, modified for our three-dimensional grid and to account for the changing viewing angles.That is, for each column of the computational grid, the flux is propagated upwards from 0.5µm to 50µm at resolution R = 15,000 and then summed by taking into account the viewing angle of each grid element.In this calculation, the same opacities as during the retrievals were used: molecular absorption via Exo-Mol cross-sections, H − opacity, Collision Induced Absorption (CIA) and Rayleigh scattering.We computed the planetary flux in those frames as if the planet were observed in eclipse (i.e., phase 0.5).

B. COMPLEMENTARY FIGURES TO THE SECTION 3
This appendix contains the complementary figures to the main article Section 3, Figure B1

Figure 1 .
Figure1.Corrected white light curves when the two WASP-121 b phase-curve visits are fitted individually.The two observations, while close, present differences that could come from either atmospheric variability or instrument systematics.Note that the second-order coefficients from the sine phase-curve models are fixed to the best-fit values from the combined fit (see FigureB3).
Figure2.Recovered temperature-pressure (T -p) profiles (left) and best-fit spectra (right) for the phases from 0.05 (blue) to 0.5 (red), obtained from the phase-curve atmospheric retrieval.In the T -p plot, the shaded regions correspond to one and three sigma confidence regions (dark to light, respectively).The radiative contribution function is also shown in dashed line, colored for each region: hotspot (red), dayside (orange), and nightside (blue).These retrievals show good agreement with the observed data and demonstrate a strong dayside thermal inversion, with the presence of a hotter region (e.g.hotspot).The best-fit T − p profiles (solid lines, left) are used to thermally force the atmospheric dynamics simulations.This figure is accompanied by a 15 s video, available online at the journal, showing the evolution of WASP-121 b emission (from blue to red) and the corresponding thermal structure as a function of phase.As the planet moves from transit to eclipse, absorption features in the data are replaced by emission features.These spectral variations enable the characterization of the thermal structure and chemistry across WASP-121 b's atmosphere.

Figure 3 .
Figure 3. Posterior distribution of the chemistry retrieved from the phase-curve data using the "Low" resolution spectra (Left), the "Medium" resolution spectra (Middle) and the spectra from ME22 (Right).The constraints obtained on metallicity Z and C/O ratio are consistent with a solar to slightly sub-solar metallicity and a super-solar C/O ratio, indicating the formation of the planet likely occurred beyond the snow-line.The blue line indicates the solar values for Z and C/O ratio.The full posterior distributions for the low-resolution retrieval are available in Figure C3

Figure 5 .
Figure5.One-dimensional (1D) thermal structure recovered by our retrieval analysis of the five eclipse observations with one sigma confidence region (left), and T -p profiles from multiple times (t ∈ [40, 185] days) at the substellar point from a threedimensional (3D) atmospheric dynamics simulation (right).The magnitude of variability in p ∈ [10 5 , 10 3 ] Pa is ∼300 K, which is consistent with the variation predicted by the 3D simulation.Dashed gray lines show the vertical extent of the atmosphere modeled by the simulations in this study.Note, while these profiles are not like-for-like comparable because the retrieved thermal structure is global and substellar temperature predictions are local.

Figure 6 .
Figure 6.Instantaneous spatial temperature maps T (λ, ϕ, p, t), where λ is the longitude, ϕ is the latitude, p is the pressure, and t is the time); p = 10 5 Pa level at t = {49, 62, 137} planetary days after the start of the simulation are shown.The maps demonstrate that very different temperature distributions arise at different times.Maps are in spherical orthographic projection, where the white circle marks the substellar point.The fields result from simulations which are thermally forced by the T -p profiles in Fig 2.The planet's "hotspot" is highly variable not only in location and time but also in shape.These storms are planetary-scale, coherent, and exhibit repetitive behaviors, leading to high amplitude and identifiable periodic signature in the planet's flux.

Figure 7 .
Figure7.Disk-averaged black-body flux variations at 1.3µm using Planck's law for a single layer in the middle of our modeled domain (10 4 Pa), and centered on four regions: dayside (substellar point), eastern terminator, western terminator, and nightside (antistellar point).The flux is normalized by the mean value of the dayside.The fluxes are quasi-periodic on periods of ∼5 days, with variations of ∼5 to ∼10%.The normalized flux for the west terminator is generally lower than that of the east terminator, indicating that a phase-curve observation of this planet is more likely to have an eastward-shifted phase offsets.The level of variability is compatible with the uncertainties of our HST observations; however, note that the HST instrument systematics makes absolute flux measurements at such p-level highly uncertain (see, e.g.,Yip et al. 2021; Edwards et al. 2023).
As explained in the Materials and Method section of the Appendix, for three time frames, t = {49, 62, 137} days, we have post-processed the simulation outputs to obtain the emitted spectral flux.Figures E2-E4 demonstrate the wide range of physical and chemical conditions that could appear on WASP-121 b, which could strengthen the changes in cloud coverage and refractory chemistry claimed by Wilson et al. (2021) and Ouyang et al. (

Figure A1 .
Figure A1.Displacement of the raw images onto the detector for the 2019 observation.The three segments, separated by the re-acquisition events, are clearly visible on those diagnostics.

A. 2 .
Light-curve analysis: the PoP pipeline To fit the light curves, we have developed a new pipeline PoP 10 (Pipeline of Pipes), leveraging the flexibility of the TauREx 3.1 framework (Al-Refaie et al.

Figure A2 .
Figure A2.White light curves for the two WASP-121 b phase-curve visits obtained with our Iraclis extraction and compared to the ones in ME22.The two visits were offset vertically (5×10 6 e − ) and in time (by -257 orbital periods) for clarity.The observation re-acquisitions (Re-acq) are also annotated.
(2014b); Tsiaras et al. (2016b); Mikal-Evans et al. ( Figure B1.Comparison of recovered white light curves, fitted with the PoP when using different assumptions for the detector and long-term ramp models.We show in solid line the median of each tested model, and for our preferred set of assumption, we show the resulting observations.Red: Two parameters exponential ramp model fromTsiaras et al. (2016b) and hybrid long-term ramp; Blue: Four parameters exponential ramp model from deWit et al. (2018)  and hybrid long-term ramp; Green: Two parameters exponential ramp model and linear long-term ramp for each segment; Orange: Two parameters exponential ramp model and quadratic long-term ramp for each segment.The hybrid long-term ramp consists in a quadratic trend for the first segment of each visit and a linear ramp for the remaining four segments, as done in ME22.This, as well as FigureB2, demonstrate that the recovered light curve depends on the assumption for the long-term ramp.

Figure B2 .
Figure B2.Posterior distributions of the white light curve fits with PoP using four different detector ramp models (see Figure B1).Red: Two parameter exponential ramp model from Tsiaras et al. (2016b); Blue: Four parameters exponential ramp model from de Wit et al. (2018).Green: Two parameters exponential ramp model and linear long-term ramp for each segment; Orange: Two parameters exponential ramp model and quadratic long-term ramp for each segment.The recovered orbital parameters are independent of the detector short-term ramp, but they are impacted by the choice of the long-term trend.

Figure B3 .
Figure B3.Corrected white light curves when the two WASP-121 b phase-curve visits are fitted together.We show our results with the Iraclis extraction and compared to the ones in Mikal-Evans et al. (2022).Top: Corrected light curves zoomed in the eclipses; Middle: Full corrected light curves; Bottom: Residuals between our best-fit light curve model and corrected data with red and blue for the Iraclis reductions and green for the reduced data in Mikal-Evans et al. (2022).

Figure B4 .
Figure B4.Posterior distributions of the white light curve fits with PoP when fitting the two phase-curve visits independently.The recovered parameters show differences in both the phase-curve model and the instrument systematics.

Figure B5 .
Figure B5.Corrected spectral light curves (left) and corresponding residuals (right) for the low spectral binning.This is the same binning employed in Mikal-Evans et al. (2022).Higher binning resolution fits (e.g.medium resolution) are also performed and presented in Figure B6.

Figure B7 .
Figure B7.Emission spectra of WASP-121 b, obtained at different phases from various reduction methods.Green: spectra at low resolution from ME22.Red: spectra at low resolution obtained using the 4-short instrument systematic model.Blue: spectra at medium resolution obtained using the 4-short instrument systematic model.Orange: spectra at medium resolution obtained using the 2-short instrument systematic model.All the reductions are consistent with each other.Note that phases 0.07 and 0.93 do not exist in ME22; hence, we instead plot their phase 0.05 and 0.095.

Figure C1 .
Figure C1.Temperature -pressure profiles (T − p) obtained by the 1.5D retrievals on the spectra from different reductions.Left panel: nightside.Middle panel: dayside.Right panel: hotspot.We show the extent of the radiative contribution function for the low-resolution retrieval with dashed lines.The retrievals on those different reductions are consistent and provide a similar picture.

Figure C2 .
Figure C2.Profiles of Brunt-Väisälä frequency squared N 2 = (g/ρ) (dρ/dz) for the retrieved thermal profiles in Fig.C1.Dotted lines at p = 10 1 and p = 10 5 Pa show the upper and lower boundary of the GCM model.In this region, flows are stably stratified (N 2 ≥ 0); hence, they satisfy the hydrostatic balance approximation of the model equations.For p > 10 6 , the nightside, hotspot and initial profiles exhibit a jump in stratification.

Figure C3 .
Figure C3.Posterior distributions for the atmosphere of WASP-121 b obtained by the 1.5D retrievals on different reduction methods.Green: spectra at low resolution from ME22.Red: spectra at low resolution obtained using the 4-short instrument systematic model.Blue: spectra at medium resolution obtained using the 4-short instrument systematic model.

Figure C5 .
Figure C5.Temperature -pressure profiles (T − p) obtained by the 1.5D retrievals on the low-resolution spectra by varying the hotspot size (AHS) and offset (DHS).Left panel: nightside.Middle panel: dayside.Right panel: hotspot.While the thermal structure are different, the conclusions on the thermal structure of WASP-121 b from those runs would be the same, independently from the hotspot assumptions.We also note that one model (AHS = 50 and DHS = 30) has a significantly higher Bayesian evidence.The models obtained: ln(E) = 1550.3for AHS = 30 and DHS = 30, ln(E) = 1510.6for AHS = 50 and DHS = 10, ln(E) = 1538.8for AHS = 50 and DHS = 20, ln(E) = 1554.5 for AHS = 50 and DHS = 30, and ln(E) = 1537.5 for AHS = 50 and DHS = 40.

Figure C6 .Figure C7 .Figure D1 .
Figure C6.Posterior distributions for the atmosphere of WASP-121 b obtained on the low-resolution spectra by varying the hotspot size (AHS) and offset (DHS).The color code is the same as in Figure C5.Independently from the hotspot assumptions, all those retrievals of WASP-121 b phase-curve data have a super-solar C/O with 0.62 < C/O < 1.11, while the retrieved metallicity is between C/O with -1.27 < log(Z) < 0.77.For all cases, we do not find evidence of a fully opaque nightside cloud deck.One model (AHS = 50 and DHS = 30) has a significantly higher Bayesian evidence.The models obtained: ln(E) = 1550.3for AHS = 30 and DHS = 30, ln(E) = 1510.6for AHS = 50 and DHS = 10, ln(E) = 1538.8for AHS = 50 and DHS = 20, ln(E) = 1554.5 for AHS = 50 and DHS = 30, and ln(E) = 1537.5 for AHS = 50 and DHS = 40.

Figure D2 .
Figure D2.Posterior distributions obtained for the three transit fits.Yellow: Transit 1. Blue: Transit 2. Green: Transit 3. Note that the chemistry parameters of the equilibrium model GGChem in those fits are fixed to the median values obtained by the 1.5D phase-curve retrieval (log Z = -0.19,C/O = 0.80).Transit 1 shows an atmosphere with moderate hazes but absorption at high altitude from water.Transit 2 shows multimodal solutions involving either a fully ionized atmosphere with un-physically high temperatures, or a cloudy/hazy atmosphere.Transit 3 displays an atmosphere with opaque clouds.

Figure D3 .
Figure D3.Posterior distributions obtained for the five eclipse fits.Yellow: Eclipse 1. Blue: Eclipse 2. Green: Eclipse 3. Red: Eclipse 4. Purple: Eclipse 5. Note that the chemistry parameters of the equilibrium model GGChem in those fits are fixed to the median values obtained by the 1.5D phase-curve retrieval (log Z = -0.19,C/O = 0.80).The five eclipses are consistent with similarly inverted thermal structures, however the posterior distributions are not the same in the middle of the atmosphere (T1 to T3).Large-scale atmospheric variability could creates those observable departures in the thermal profiles.

Figure E2 .
Figure E2.WASP-121 b chemical maps centered around the sub-stellar point for H2, H, H − and e − at two different times (t = 49 and t = 62 days) and two pressure levels (p = 10 5 Pa and p = 5×10 3 Pa).Those maps are obtained by post-processing the temperature fields (top row) from the BoB dynamics calculations with the TauREx library.

Figure E3 .
Figure E3.WASP-121 b chemical maps centered around the sub-stellar point for H2O, CO, CH4 and CO2 at two different times (t = 49 and t = 62 days) and two pressure levels (p = 10 5 and p = 5 × 10 3 Pa).Those maps are obtained by post-processing the temperature fields (top row) obtained from the BoB dynamics calculations with the TauREx library.

Figure E4 .
Figure E4.WASP-121 b chemical maps centered around the sub-stellar point for TiO, VO and FeH at two different times (t = 49 and t = 62 days) and two pressure levels (p = 10 5 and p = 5 × 10 3 Pa).Those maps are obtained by post-processing the temperature fields (top row) obtained from the BoB dynamics calculations with the TauREx library.
Figure 4. Examples of observed transit (left) and eclipse (right) spectra of WASP-121 b analyzed in this work.The spectra shown are corrected for vertical offsets to highlight the differences in spectral shapes.Variability in the planet's weather patterns could create such variations in the spectroscopic data; this is strongly suggested in high-resolution simulations carried out for this planet in this work (see, e.g., Figures , we model the phase-curve vari- . The radiative transfer includes absorption from the main expected opacity sources via ExoMol line lists (Tennyson & Yurchenko 2012; Chubb et al. 2021), namely: H 2 O (Polyansky et al. 2018), CO (Li et al. 2015), CO 2 (Yurchenko et al. 2020), CH 4 (Yurchenko et al. 2017), TiO (McKemmish et al. 2019), VO (McKemmish et al. 2016), FeH (Bernath 2020), and H − (John 1988; Edwards et al. 2020).We also consider Collision Induced Absorption (CIA) by H 2 -H 2 and H 2 -He pairs