A JWST NIRSpec Phase Curve for WASP-121b: Dayside Emission Strongest Eastward of the Substellar Point and Nightside Conditions Conducive to Cloud Formation

We present the first exoplanet phase curve measurement made with the JWST NIRSpec instrument, highlighting the exceptional stability of this newly-commissioned observatory for exoplanet climate studies. The target, WASP-121b, is an ultrahot Jupiter with an orbital period of 30.6 hr. We analyze two broadband light curves generated for the NRS1 and NRS2 detectors, covering wavelength ranges of 2.70-3.72 micron and 3.82-5.15 micron, respectively. Both light curves exhibit minimal systematics, with approximately linear drifts in the baseline flux level of 30 ppm/hr (NRS1) and 10 ppm/hr (NRS2). Assuming a simple brightness map for the planet described by a low-order spherical harmonic dipole, our light curve fits suggest that the phase curve peaks coincide with orbital phases $3.36 \pm 0.11$ deg (NRS1) and $2.66 \pm 0.12$ deg (NRS2) prior to mid-eclipse. This is consistent with the strongest dayside emission emanating from eastward of the substellar point. We measure planet-to-star emission ratios of $3,924 \pm 7$ ppm (NRS1) and $4,924 \pm 9$ ppm (NRS2) for the dayside hemisphere, and $136 \pm 8$ ppm (NRS1) and $630 \pm 10$ ppm (NRS2) for the nightside hemisphere. The latter nightside emission ratios translate to planetary brightness temperatures of $926 \pm 12$ K (NRS1) and $1,122 \pm 10$ K (NRS2), which are low enough for a wide range of refractory condensates to form, including enstatite and forsterite. A nightside cloud deck may be blocking emission from deeper, hotter layers of the atmosphere, potentially helping to explain why cloud-free 3D general circulation model simulations systematically over-predict the nightside emission for WASP-121b.


INTRODUCTION
Corresponding author: Thomas Mikal-Evans tmevans@mpia.de * Flatiron Research Fellow It is crucial to explore the coupling between the dayside and nightside hemispheres of tidally locked planets if we are to develop a global understanding of their atmospheres. One of the most effective means of constraining the chemical, thermal, and dynamical properties of both hemispheres is to measure planetary emission spectra arXiv:2301.03209v2 [astro-ph.EP] 16 Feb 2023 over the course of a full planetary orbit (e.g. Knutson et al. 2007;Cowan et al. 2012;Stevenson et al. 2014;Kreidberg et al. 2018;Arcangeli et al. 2019;Feng et al. 2020;Irwin et al. 2020;Mikal-Evans et al. 2022). To this end, our team observed a spectroscopic phase curve for the ultrahot Jupiter WASP-121b using JWST NIRSpec with the G395H grating. In this Letter, we present an analysis of two broadband light curves generated from this dataset, one for each of the NRS1 and NRS2 detectors. As well as providing some preliminary scientific results, our purpose is to provide a brief report on the stability of JWST and NIRSpec when performing longstare time-series measurements lasting tens of hours.
The target, WASP-121b, is an inflated gas giant with a mass of 1.183 +0.064 −0.062 M J , a radius of 1.753 ± 0.036 R J measured at red-optical wavelengths, and an orbital period of 30.59820 ± 0.00001 hr (Delrez et al. 2016;Bourrier et al. 2020). Secondary eclipse measurements made with HST have shown that the dayside atmosphere has a thermal inversion, with a near-infrared brightness temperature close to 2,700 K (Evans et al. 2017;Mikal-Evans et al. 2020). Evidence has been uncovered for numerous UV-optical absorbers in the atmosphere that are likely responsible for the dayside thermal inversion, such as VO, SiO, V, Fe, Mg, and Ca (e.g. Evans et al. 2018;Hoeijmakers et al. 2020;Lothringer et al. 2021).
Phase curve measurements with HST have shown that temperatures on the nightside of WASP-121b are likely low enough for refractory condensates containing V, Fe, Mg, and Ca to form (Mikal-Evans et al. 2022). Since these latter elements are observed in the gas phase at the day-night terminator, if such clouds do form, they are presumably recirculated back to the dayside hemisphere where they are vaporized before settling to the deeper layers of the atmosphere. Meanwhile, nondetections of Ti and Al at the day-night terminator reported by Hoeijmakers et al. (2020Hoeijmakers et al. ( , 2022 are consistent with these elements being cold-trapped in deeper layers of the nightside hemisphere as perovskite (CaTiO 3 ) and corundum (Al 2 O 3 ).
From these previous studies, we are already observing the important connection between the dayside and nightside hemispheres of WASP-121b. One of the ultimate goals of our JWST program will be to provide significantly tighter constraints on the dayside and nightside properties. This should be possible thanks to the JWST data having a much higher signal-to-noise and spectral resolution than the existing HST data, as well as the more favorable planet-to-star emission ratio at the longer wavelengths covered by NIRSpec. Our team observed the WASP-121 system over a full planetary orbit on 2022 October 14-15 using the NIRSpec instrument as part of program GO-1729 (P.I. Mikal-Evans, co-P.I. Kataria). Observations were made without interruption for 37.8 hr, commencing 145 min prior to secondary eclipse ingress and continuing until 105 min after egress of the following secondary eclipse. We used the G395H grating, providing wavelength coverage between 2.70-5.15 µm at R ∼ 3000. For reading the detector, we employed the SUB2048 subarray option with the NRSRAPID readout pattern. Each integration consisted of 42 groups, translating to integration times of 38.8 s and an overall duty cycle of 99%.
To extract the target spectra from the data frames, we used the (Fast InfraRed Exoplanet Fitting for Lightcurves) FIREFLy code (Rustamkulov et al. 2022a,b). The data reduction begins with a customized reduction of stage 1 from the JWST Python pipeline (Bushouse et al. 2022). Detector "1/f " noise Rustamkulov et al. 2022a;Schlawin et al. 2020) is removed at the group level, masking out the spectral trace pixels and removing the median value from each column of the remaining background pixels. Bad pixels and cosmic rays are then flagged as outliers, both in the time-series and spatially on the 2D images. Movement of the spectra along the x-and y-axes of the detector are measured using cross-correlation, exhibiting standard deviations at the level of 0.002 pixels in both directions for NRS1 and in the y-direction for NRS2, while a higher standard deviation of 0.005 pixels is measured in the x-direction for NRS2 ( Figure 1). This level of jitter is similar to previously reported on-sky performance (Rigby et al. 2022;Espinoza et al. 2022), though demonstrated here on a much longer timescale. The 2D data frames are then adjusted using interpolation to account for the small pointing shifts. We find that this alignment procedure can dampen x and y correlations in the subsequent spectrophotometry, as it effectively takes into account the redistribution of flux across adjacent pixels. However, residual x and y trends may still remain due to intrapixel sentivities ). Next, a fourth-order polynomial was used to trace the center of the point-spread-function (PSF) along the detector. Counts were then summed within an aperture centered on the trace with a width of 4.9 pixels for NRS1 and 5.5 pixels for NRS2, linearly interpolating the counts in those pixels where only a fraction was contained within the aperture. The aperture widths were selected after trialing a number of values and selecting those that minimized the scatter of the in-eclipse photometric time-series. To obtain the wavelength calibration, we extrapolated JWST pipeline data products across the detector edge pixels that did not have an assigned wavelength.

LIGHT CURVE ANALYSIS
We generated two light curves by summing each extracted 1D spectrum between dispersion columns 300-2042 for NRS1 and 5-2010 for NRS2, conservatively encompassing the full passbands of both detectors. The resulting raw light curves are shown in the top two panels of Figure 2. We modeled these light curves using the starry package (Luger et al. 2019) for the star-planet signal, multiplied by a linear trend in time to account for instrumental drift. For the planet brightness distribution, a simple dipole map was adopted, comprising the Y 0,0 and Y 1,0 spherical harmonics. For the stellar brightness distribution, a quadratic limb-darkening profile was assumed. In defining the data log-likelihood, we assumed that the data were normally distributed with standard deviations given by the pipeline Poisson un-certainties and an additional high-frequency systematics noise term combined in quadrature.
We fitted the NRS1 and NRS2 light curves jointly, with the following parameters shared across both light curves: the stellar mass (M ) and logarithmic radius (log 10 R ); the differences in the orbital period (∆P ) and transit mid-time (∆T c ) from the values reported by Bourrier et al. (2020); the orbital inclination (i); and the planetary mass (M p ). Separate sets of the remaining model parameters were fitted for each light curve in the joint fit, namely: the stellar limb-darkening coefficients (u 1 , u 2 ); the logarithmic planetary radius (log 10 R p ); the overall amplitude of the planetary brightness map (A) and the coefficient of the Y 1,0 spherical harmonic (y 1,0 ); the rotational offset of the planetary brightness map relative to a map with a brightness peak centered at the substellar point (∆φ); the coefficients of the linear instrumental drift trend (f 0 , f 1 ); and the high-frequency systematics noise term (σ syst ). The priors we adopted for each parameter are shown by the dashed lines in Figure 3. For M , log 10 R , i, and M p , these were normal priors based on the posteriors reported in Bourrier et al. (2020). Broad normal priors were also adopted for log 10 R p , ∆T c , u 1 , u 2 , A, y 1,0 , ∆φ, f 0 , and f 1 . Halfnormal priors only allowing positive values were adopted for σ syst . We marginalized the resulting posterior distribution using the No-U-Turn-Sampling (NUTS; Hoffman & Gelman 2014) method implemented by PyMC3 (Salvatier et al. 2016).

RESULTS
The maximum a posteriori (MAP) light curve models are shown in Figure 2. Posterior distributions for each of the parameters that were varied in the fitting are shown as dark blue histograms in Figure 3, while the light blue histograms show the corresponding distributions for other parameters that were calculated from the latter, namely: the stellar radius R ; the planetto-star radius ratio R p /R ; the planet-to-star flux ratio F p /F for the dayside and nightside hemispheres of the planet; and the corresponding planetary brightness temperatures T b . The brightness temperatures were calculated by assuming F p /F = (R p /R ) 2 B p (T b )/F , where B p denotes a blackbody function for the planetary flux, integrated over the instrument passband. In deriving distributions for T b , uncertainties in the host star properties were accounted for by first randomly drawing values for the stellar effective temperature (T ), surface gravity (log g ), and metallicity ([M/H]) from the posterior distributions reported in Delrez et al. (2016), which were assumed to be normal. For each of these sets of stellar properties, we used pysynphot (STScI Develop- . Panels (e) and (f ): Black points show the residuals between the data and best-fit models for the NRS1 and NRS2 detectors, respectively. Orange lines show the residuals after applying a Gaussian filter to smooth the random noise. Green and purple shading show the times of transit and eclipse, respectively, with darker ranges corresponding to ingress and egress times.   −39 K for the NRS2 passband. We note that the brightness temperature uncertainties are dominated by the stellar model uncertainties. When we only account for the uncertainties in F p /F and R p /R derived from the light curve fits while holding the stellar properties fixed to the best-fit values reported by Delrez et al. (2016), the uncertainties on the dayside brightness temperatures reduce by an order-of-magnitude from around 30-40 K (Figure 3) to 3 K. By contrast, the uncertainties associated with dayside brightness temperatures measured previously with the Spitzer Space Telescope in the 3-5 µm wavelength range have been dominated by the F p /F uncertainty. For example, Garhart et al. (2020) reported dayside brightness temperatures for WASP-121b of 2, 490 ± 77 K and 2, 562 ± 66 K in the 3.6 µm and 4.5 µm Infrared Array Camera (IRAC) passbands, which cover similar wavelength ranges to NRS1 and NRS2.
The long-term instrumental drift is found to be mild for both detectors. For NRS1, we obtain a negative drift of 1,172 ppm over the course of the observation, corresponding to a rate of approximately 30 ppm/hr. NRS2 is even better behaved, exhibiting a positive drift of only 416 ppm from start to finish, at a rate of approximately 10 ppm/hr. The standard deviations of the light curve residuals are 127 ppm and 161 ppm for the NRS1 and NRS2 detectors, respectively. This is 39% and 17% higher than the Poisson uncertainties derived from the instrument pipeline, respectively. As such, we derive high-frequency systematics noise values (σ syst ) of 90 ± 2 ppm for NRS1 and 80 ± 4 ppm for NRS2.

DISCUSSION
The small positive values that we measure for ∆φ in both the NRS1 (3.36 ± 0.11 • ) and NRS2 (2.66 ± 0.12 • ) passbands translate to phase curves with maxima shifted prior to mid-eclipse. This in turn suggests that the eastern half of WASP-121b's dayside hemisphere is hotter on average than the western hemisphere, which is qualitatively consistent with the predictions of general circulation model (GCM) simulations of hot Jupiters (e.g. Showman & Guillot 2002;Showman et al. 2009;Komacek et al. 2022). However, GCMs typically predict eastward phase curve shifts that are significantly higher than what we have measured here for WAP-121b; closer to 10 • or more (e.g. Parmentier & Crossfield 2018). Small phase curve shifts relative to GCM predictions have also been measured for a number of other hot Jupiters with properties similar to WASP-121b (e.g. Zhang et al. 2018;Kreidberg et al. 2018). For example, Kreidberg et al. (2018) measured an eastward phase shift of 2.0 ± 0.7 • for WASP-103b in the Spitzer IRAC 3.6 µm channel, compared to a shift of 9.2 • predicted by a GCM presented in the same study. For ultrahot planets like WASP-121b and WASP-103b -which both have dayside brightness temperatures well above 2, 500 K -Lorenz drag resulting from the motion of thermally ionized gas through the planetary magnetic field may inhibit the advection of heat throughout the atmosphere, reducing the observed phase curve offsets (Perna et al. 2010). Alternatively, even if the hottest region of the dayside hemisphere is located significantly eastward of the substellar point, the presence of nightside clouds can reduce the observed phase shift in the resulting light curve (Parmentier et al. 2021). We plan to investigate the latter possibility further in future analyses by including higher-order spherical harmonics terms for the planetary brightness map, which will allow for a sharper delineation in the brightness of the dayside and nightside hemispheres.
Our dayside and nightside emission measurements are plotted in Figure 4 against the 3D GCM predictions for WASP-121b from , which assume 1× and 5× solar metallicity, chemical equilibrium, and do not include clouds. The models exhibit spectral features of H 2 O and CO, which are in emission on the dayside due to a thermal inversion and absorption on the nightside due to a cooling temperature profile (Figure 4). The agreement between the measured dayside emission (3,924±7 ppm) and the 5× solar GCM prediction (3,916 ppm) across the NRS1 passband is very close to 1σ. However, the dayside emission measured across the NRS2 passband (4,924±9 ppm) is 19σ higher than the 5× solar GCM prediction (3,916 ppm). For the nightside emission, the measurements across the NRS1 (136 ± 8 ppm) and NRS2 (630 ± 10 ppm) passbands are > 10σ lower than both the 1× solar metallicity (588 ppm for NRS1 and 1,080 ppm for NRS2) and 5× solar metallicity (413 ppm for NRS1 and 789 ppm for NRS2) GCM predictions. Such large discrepancies are unsurprising, given the small measurement uncertainties and the fact that the GCMs were not tuned to match the data. Still, it is notable that the measured nightside brightness temperatures are both substantially lower than predicted by the GCM simulations. Furthermore, the  GCMs did not include the dissociation and recombination of hydrogen, which if anything would result in the predicted nightside temperatures being higher than those shown in Figure 4 due to the release of latent heat (Bell & Cowan 2018;Tan & Komacek 2019;Komacek et al. 2022), exacerbating the discrepancy.
Clouds -which were also not included in the Parmentier et al. (2018) GCMs -might help explain our lower-than-predicted nightside brightness temperature measurements. Previous HST observations have shown that the nightside temperature profile of WASP-121b does not have a thermal inversion at the near-infrared photosphere, instead cooling with decreasing pressure (Mikal-Evans et al. 2022). If an optically thick cloud deck blankets the nightside hemisphere, it could block the emission from deeper, hotter layers of the atmosphere, lowering the observed brightness temperature (Mendonça et al. 2018;Keating et al. 2019;Beatty et al. 2019;Parmentier et al. 2021). Even if the nightside cloud distribution is patchy rather than uniform, the brightness temperature can be reduced significantly relative to that of a cloud-free nightside (Komacek et al. 2022).
The nightside brightness temperatures that we measure for WASP-121b in the NRS1 (926 +12 −12 K) and NRS2 (1, 122 +10 −10 K) passbands fall below the condensation temperatures of silicates, such as enstatite and forsterite (Visscher et al. 2010;Wakeford et al. 2017), which are expected to be abundant in the nightside atmospheres of hot Jupiters (e.g. Gao & Powell 2021). Therefore, our measurements reveal nightside conditions for WASP-121b that do appear to be conducive to the formation of silicate clouds, strengthening the possibility that clouds may prove key to understanding why our measured nightside brightness temperatures are significantly lower than predicted by the  GCMs shown in Figure 4. Further work will be required to investigate if the difference in the nightside brightness temperatures obtained for the two passbands could be caused by the wavelength-dependent opacity of atmospheric layers above a cloud deck. For these future  Note that the measurement uncertainties are smaller than the diamond symbols. Predictions from the cloud-free 3D GCM simulations of  are also shown for heavy element enrichments of 1× solar (blue lines) and 5× solar (orange lines). Circle and square symbols show, respectively, the 1× and 5× solar GCM predictions binned to the light curve passbands, which are shown as green lines in the top panel.
analyses, the dayside and nightside emission spectra will be extracted and interpreted, rather than the broad passbands considered here.
It should also be stressed that the analysis we have presented here is preliminary. In particular, we have adopted a very simple light curve model and it is evident by eye that low-amplitude correlations still remain in the time-series of residuals ( Figure 2). We did experiment with including the x and y pointing coordinates shown in Figure 1 as additional linear decorrelation variables in the light curve fits, but found that this did not appreciably reduce the scatter in the residuals, nor did it affect the derived parameter distributions shown in Figure 3. Tidal deformation of WASP-121b could also subtly affect the observed phase curve (Cowan et al. 2012;Wahl et al. 2021). We did use starry to perform some preliminary light curve fits treating the planetary oblateness as a free parameter, but found that our basic conclusions were unaffected. Further work is required to determine if an oblate shape for the planet is supported by the data. Additional possibilities yet to be considered include adding higher-order spherical har-monics terms for the planetary brightness map and allowing for nonlinear instrumental baseline trends. These investigations are ongoing and could conceivably affect the constraints that we ultimately obtain for parameters of interest, such as the planetary nightside emission and brightness phase offsets.
In the meantime, our phase curve measurement for WASP-121b demonstrates the overall high level of stability that JWST NIRSpec is capable of maintaining for a single-stare observation lasting 37.8 hr. The approximately linear drift observed in the baseline flux level is extremely mild in comparison to the instrumental systematics that have affected past HST and Spitzer datasets (e.g. Cowan et al. 2012;Stevenson et al. 2014;Kreidberg et al. 2018;Mikal-Evans et al. 2022). We are hopeful that the additional 80-90 ppm high-frequency systematics noise observed here in our preliminary analysis can be further reduced as more on-sky calibrations become available, along with continued refinement of the data analysis, such as improved treatment of the 1/f detector noise. Based on other recent analyses of NIRSpec data (Alderson et al. 2022;Rustamkulov et al. 2022b), we also anticipate that simple models fitted to light curves generated over narrower wavelength ranges than the broad passbands we have considered in the present study will achieve residual scatters closer to Poisson predictions.

CONCLUSION
We have measured a full-orbit phase curve for the ultrahot Jupiter WASP-121b using JWST NIRSpec. The resulting light curves generated across broad passbands for the NRS1 and NRS2 detectors exhibit minimal systematics over the 37.8 hr observation. We find that the phase curve peaks are shifted prior to mid-eclipse by 3.36 ± 0.11 • (NRS1) and 2.66 ± 0.12 • (NRS2), suggesting that the eastern region of the dayside hemisphere is hotter on average than the western region. The mea-sured dayside emission in the NRS1 passband is in good agreement with a cloud-free GCM assuming 5× solar metallicity; however, the same GCM underpredicts the dayside emission in the NRS2 passband by 19σ. For the nightside emission, cloud-free GCM simulations assuming 1× and 5× solar metallicity significantly overpredict the data. This observation could possibly be explained by nightside clouds blocking the emission from deeper, hotter layers of the atmosphere. The corresponding nightside brightness temperatures are < 1200 K in both passbands, which is cool enough for various condensates to form, including silicates such as enstatite and forsterite.