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

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; Swain et al. 2008, 2009; Kreidberg et al. 2014a, 2014b; Madhusudhan et al. 2014; Stevenson et al. 2014; Line et al. 2016; Sing et al. 2016; Tsiaras et al. 2018; Benneke et al. 2019; Tsiaras et al. 2019; Edwards et al. 2020; Pluriel et al. 2020; Skaf et al. 2020; Line et al. 2021; Mansfield et al. 2021; Mugnai et al. 2021; Changeat et al. 2022, hereafter C22; 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 the observing time constraints.

With the Spitzer telescope, a number of studies have analyzed the repeated measurements of individual transiting exoplanets via photometric multiepoch 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. 2023). 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 data sets 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 nontransiting 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-1257b, 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 hr is reported by Zhou et al. (2022). Many other planetary mass companions exhibit similar levels of variability.

Recently, an 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, hereafter ME22), 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⁴ to ∼10³ 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 (Cho et al. 2021; Skinner & Cho 2021, 2022; Skinner et al. 2023). In these simulations, a long-lived giant storm pair forms near the substellar point, drifts initially toward one of the terminators, and then rapidly translate westward thereafter—traversing the nightside and ultimately breaking up or dissipating near the eastern terminator; this cycle is quasiperiodic. 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 high-resolution simulations (with a different initial condition than in the above studies) by Cho et al. (2003), who suggested that the weather patterns would lead to a potentially observable variability on hot exoplanets.

In this backdrop, we present here results from an in-depth 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 2016 June, one eclipse in 2016 November, and two phase curves in 2018 March and 2019 February. Significantly, we utilize the entirety of these observations in this work. Previous analyses from a combination of facilities—HST, TESS, Spitzer, JWST, and ground-based facilities (Evans et al. 2018; Tsiaras et al. 2018; Sing et al. 2019; Ben-Yami et al. 2020; Cabot et al. 2020; Hoeijmakers et al. 2020; Mikal-Evans et al. 2020; Borsa et al. 2021; Daylan et al. 2021; Azevedo Silva et al. 2022; C22; Gibson et al. 2022; ME22; Mikal-Evans et al. 2023)—have detected the presence of water vapor, 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, the discussion and conclusion are presented in Section 7. Additional material is included in Appendix A as well.

The WASP-121 b data we have analyzed in this work are obtained using HST WFC3-G141. These data are publicly available at the Mikulski Archive for Space Telescopes (MAST). ¹⁷ 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, 2021; Edwards et al. 2023). The analyzed data are from one eclipse, one transit, and two phase curves—comprising a total observing time of about 90 hr; each phase-curve 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 that are not from the phase-curve observations have been previously extracted using Iraclis by C22, and the phase-curve data have been recently analyzed by ME22. However, the phase curves have not been extracted using Iraclis. Therefore, we have reanalyzed 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:

• 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 Appendix A); and, testing various assumptions to model the instrument systematics in the literature;
• 2.  
  analyzing the recovered phase-curve data set with the 1.5D retrievals developed by Changeat & Al-Refaie (2020), Changeat et al. (2021) in the TauREx3.1 framework (Al-Refaie et al. 2021, 2022a), permitting time-independent, global properties (e.g., mean metallicity Z and C/O ratio as well as mean thermal profiles at different longitudes) to be extracted;
• 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;
• 4.  
  performing high-resolution, global atmospheric dynamics simulations with the pseudospectral code Built On Beowolf (BoB; e.g., Polichtchouk et al. 2014; Skinner & Cho 2021), 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 Appendix A.

3.1. Combined White Light-curve Correction

3.2. Individual White Light-curve Corrections

To investigate the variability of the atmosphere and the instrument systematics in the available phase-curve 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₃ and C₄) 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 short-term systematics. Figure 1 shows the corrected white light curves for the individual fits; Figure B4 shows the corresponding posterior distributions.

Zoom In Zoom Out Reset image size

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.

3.3. Spectral Light-curve Fitting

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 Appendix A, we use the 1.5D atmospheric retrievals to interpret the observed data. ¹⁸ This 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 (in particular, 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.

Zoom In Zoom Out Reset image size

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 C22). 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⁵ and 10³ 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₂ thermal dissociation): the temperature in the inversion region of the dayside is indeed hot enough to dissociate most molecules—including water and even more stable volatiles (CO and CO₂) and refractory molecules (FeH, TiO, and VO), along with H₂ (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 subsolar metallicity Z and a supersolar C/O ratio. For the low-resolution case, we estimate the mean chemical characteristics of the planet to be $\mathrm{log}(Z)=-{0.19}_{-0.13}^{+0.11}$, and ${\rm{C}}/{\rm{O}}={0.80}_{-0.05}^{+0.03}$. A more conservative estimate, encompassing the uncertainties from all the reductions and retrievals tested in this work is $-0.77\lt \mathrm{log}(Z)\lt 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 supersolar 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 protoplanetary disk (e.g., Öberg et al. 2011; Mordasini et al. 2016; Brewer et al. 2017; Madhusudhan 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 reuse 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 reintroduced as an observation-driven forcing to enable the physically realistic and case-specific simulations.

As mentioned, C22 has previously reduced the individual transit and eclipse data sets 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 independent ¹⁹ and have better constraints from the phase-curve data, we have decided to reinject 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.

Zoom In Zoom Out Reset image size

For the transits, the three individual observations show the presence of covering hazes or clouds. Specifically, the transit spectra captured during the two phase curves (blue and green) are fully cloudy (i.e., consistent with featureless), while the observation obtained in 2016 (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₂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 ${\rm{\Delta }}\mathrm{ln}(E)=18.3$, and possesses strong water vapor absorption features. Transit 2 (blue) displays an interesting multimodal 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 toward the blue end of the spectrum. Given the unphysical 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, ${\rm{\Delta }}\mathrm{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 ${\rm{\Delta }}\mathrm{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. 2024). However, a 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 nightside, were not recovered by the retrievals on the phase-curve data. Many reasons could explain this result: (i) the emitted signal in the 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⁵ to 10³ Pa region. Importantly, this range is much larger than the average 1σ 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 phase-curve 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.

Zoom In Zoom Out Reset image size

While instrument systematics could remain (see the 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.

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 Section A.5, as well as Skinner & Cho (2021); 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 robustness 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., Figure 2) and prescribed. Profiles at many different times, which have deviated from the prescribed one due to the nonlinear atmospheric motion, are compiled in Figure 5 (right panel): they should be compared with observations (left panel).

Figure 6 shows temperature maps at p = 10⁵ 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³ Pa, and p = 10⁵ Pa, show the vertical distribution, which is strongly barotropic—i.e., vertically aligned (Figure E1); a full movie of this simulation is provided with Figure E1 online. We also show chemical species maps, which are simply postprocessed using the instantaneous density and temperature distributions from the simulation (Figures E2–E4) at p = 1 bar = 10⁵ Pa, and p = 50 mbar = 5 × 10³ 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₂, H, H⁻, e⁻, H₂O, CO, CO₂, CH₄, 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.

Zoom In Zoom Out Reset image size

As can be seen in the figures, the complex motion of the atmosphere—and the organized structures therein—causes 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 point—either sequestered in a long-lived storm or generated near the 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 ≳ 10⁴ Pa vertical region. Skinner et al. (2023) have recently reported a quasiperiodic 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⁵ 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⁵ Pa, for example; see Figure 7 with Figure 4 in Skinner et al. (2023). ²⁰ In both figures, the fluxes exhibit quasiperiodic variations, with excess flux at the substellar point ("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 (see 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. (2023). 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."

Zoom In Zoom Out Reset image size

As expected, the magnitude of the variability in the model domain varies with the p-level: it is greatest at p ∼ 10⁵ Pa. Above p ≈ 2 × 10⁴ 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₂ (Figure E2). On WASP-121 b, while the stellar irradiation can still penetrate as deep as 10⁵ Pa (as in a typical hot Jupiter), the upper layers of the atmosphere likely maintain a much stronger dayside–nightside temperature gradient overall because of the shorter relaxation time (which is ∝T⁻⁴; 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 3σ altitude-averaged temperature range is 680 K between p ∈[10⁵, 10³] 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 versus 1D), the scale of the temperature variations from the eclipse observations and from the theory is compatible.

As explained in Appendix A, for three time frames, t = {49, 62, 137} days, we have postprocessed 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), Ouyang et al. (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. 2023), 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.

Zoom In Zoom Out Reset image size

Note that we obtain a lower level of variability when performing a 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 nonalignment) 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₂ thermally dissociates faster than H₂O; hence, we expect larger variability signals for wavelengths between 1.0 and 1.3 μm.

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, cross-verifiable 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 treatments in the analysis of the observations, 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 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 time ²¹ —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 ground-based 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³, 10⁵] 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 quasiperiodic 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.

This work resulted from an initial discussion at the Exoplanet Symposium 2022, organized at the Flatiron Institute.

Q.C. is funded by the European Space Agency (ESA) under the 2022 ESA Research Fellowship Program. J.W.S. is supported by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). O.V. is funded by the Centre National d'Études Spatiales (CNES), the CNRS/INSU Programme National de Planétologie (PNP), and the Agence Nationale de la Recherche (ANR), EXACT ANR-21-CE49-0008-01. J.N. is funded by the Columbia University/Flatiron Institute joint Research Fellowship program. G.M. is funded by the European Union Horizon 2020 program, Marie Skłodowska-Curie grant agreement. This research received support from the UK Science and Technology Funding Council (STFC; ST/S002634/1 and ST/T001836/1), from the UK Space Agency (UKSA; ST/W00254X/1), and from the European Research Council (ERC; Horizon 2020 758892 ExoAI project). B.E. and A.T. received travel support for this work under the ESA Science Faculty Funds, funding reference ESA-SCI-SC-LE117 and ESA-SCI-SC-LE118.

The Flatiron Institute, DIRAC, and OzSTAR facilities provided the computing resources; this work utilized the Cambridge Service for Data-Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC High Performance Computing Facility (www.dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1. DiRAC is part of the National e-Infrastructure. The OzSTAR program receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government.

This work is based upon observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute (STScI) operated by AURA, Inc. The raw data used in this work are available from the Hubble Archive, which is part of the Mikulski Archive for Space Telescopes. The specific observations analyzed can be accessed via doi:10.17909/7an2-2m33. We are thankful to those who operate these telescopes and their corresponding archives, the public nature of which increases scientific productivity and accessibility (Peek et al. 2019). Our HST reduction pipeline, Iraclis, ²² our light-curve fitting code, PoP, ²³ and our atmospheric retrieval code, TauREx3, ²⁴ are publicly available on Github and can be installed from pypi. The chemical code GGChem ²⁵ in its TauREx3 plugin version is also publicly available. The atmospheric dynamics code BoB is also publicly available from NCAR (Scott et al. 2004).

A.1. Data Reduction with Iraclis

For the reduction and extraction of the spatially scanned spectroscopic images, we use the dedicated and publicly available pipeline Iraclis (Tsiaras et al. 2016a, 2016c, 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-pixel correction, nonlinearity correction, dark current subtraction, gain conversion, sky background subtraction, calibration, flat-field correction, and bad-pixel/cosmic-ray 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 reacquisition events, separating the observations in three distinct segments (see Figure A1). Reacquisitions 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.

Zoom In Zoom Out Reset image size

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; C22; 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.

Zoom In Zoom Out Reset image size

A.2. Light-curve Analysis: the PoP Pipeline

To fit the light curves, we have developed a new pipeline, Pipeline of Pipes (pop), ²⁶ leveraging the flexibility of the TauREx 3.1 framework (Al-Refaie et al. 2022a). While TauREx was originally designed for atmospheric retrievals (Waldmann et al. 2015b, 2015a), 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 the 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. 2016b) 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., ME22), we model the phase-curve variations using a first- or second-order sinusoidal function:

$$\begin{eqnarray}&&\begin{array}{l}{\mathrm{LC}}_{S}={C}_{0}+\displaystyle \frac{{C}_{1}}{2}\left[1-\cos ({\rm{\Phi }}-{C}_{2})\right]+\displaystyle \frac{{C}_{3}}{2}\\ \left[1-\cos (2{\rm{\Phi }}-{C}_{4})\right],\end{array}\end{eqnarray} \tag{ 1 }$$

where Φ is the orbital phase, C₀ is the minimum of the nightside flux, C₁ is the maximum of the dayside flux, and C₂ is the phase-curve offset. C₃ and C₄ correspond to the second-order terms of the sine phase variations. As in Dang et al. (2018), the full phase-curve model M_PC is constructed by the following:

$$\begin{eqnarray}&&{M}_{\mathrm{PC}}={\mathrm{LC}}_{T}+{\mathrm{LC}}_{S}\times {\mathrm{LC}}_{E}.\end{eqnarray} \tag{ 2 }$$

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, 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. 2016c; Kreidberg et al. 2018; Arcangeli et al. 2019). In the case of WASP-121 b, since reacquisition 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 as follows:

$$\begin{eqnarray}&&{{IS}}_{\mathrm{long}}=(1-{A}_{0}^{i}\times (t-{t}_{s}^{i})-{A}_{1}^{i}\times {\left(t-{t}_{s}^{i}\right)}^{2})\times {N}^{i}\end{eqnarray} \tag{ 3 }$$

where A₀ ⁱ is the linear coefficient of the segment i, A₁ ⁱ is the quadratic coefficient of the segment i, and Nⁱ 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. 2016c):

$$\begin{eqnarray}&&{{IS}}_{\mathrm{short}}=1+{B}_{0}\times \exp \left(\displaystyle \frac{t-{t}_{o}^{j}}{{B}_{1}}\right)\end{eqnarray} \tag{ 4 }$$

where B₀ and B₁ are the exponential factors, common to all visits in each observation, and t_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 the following:

$$\begin{eqnarray}&&{{IS}}_{\mathrm{short}}=1+{B}_{0}\times \exp \left(\displaystyle \frac{t-{t}_{o}^{j}}{{B}_{1}R}\right)\end{eqnarray} \tag{ 5 }$$

where R is the correcting function mainly acting on the first orbit, and it is given by

$$\begin{eqnarray}&&R=1+{B}_{2}\times \exp \left(\displaystyle \frac{t-{t}_{v}}{{B}_{3}}\right).\end{eqnarray} \tag{ 6 }$$

B₂ and B₃ 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 independent—i.e., the midtransit time (t_mid), the inclination (i), and the semimajor 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. 2016c). However, we account for additional systematic noise or other unaccounted effects by rescaling 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 the 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 divide-white strategy similar to Kreidberg et al. (2014b), Tsiaras et al. (2016c), and ME22, 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 rescaled 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 the following: Φ ∈ {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 parts per million in the case of WASP-121 b (Morello et al. 2019; ME22; Morello et al. 2021).

A.3. Phase-curve Retrievals

A.4. Transit and Eclipse Retrievals

A.5. Dynamics Modeling

This appendix contains the complementary figures to Section 3, Figures B1–B7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This appendix contains the complementary figures to Section 4, Figures C1–C7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This appendix contains the complementary figures to Section 5, Figures D1–D3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This appendix contains the complementary figures for Section 6, Figures E1–E4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size