Phase-resolving the absorption signatures of water and carbon monoxide in the atmosphere of the ultra-hot Jupiter WASP-121b with GEMINI-S/IGRINS

We observed three transits of WASP-121b with the Immersion GRrating INfrared Spectrometer (IGRINS) on the 8.1-m Gemini-South Telescope, Cerro Pachón, Chile (Yuk et al. 2010; Mace et al. 2018). The observations were taken on 25 February, 10 December, and 24 December 2022 (dates in UTC) as part of programs GS-2022A-Q-242 and GS-2022B-Q-133 (PI: Wardenier). The observations are summarised in Table 2. Additionally, Fig. 1 shows the signal-to-noise ratio (SNR), the airmass, and the relative humidity during each night as a function of the planet’s orbital-phase angle.

IGRINS has a spectral resolution $R\sim$ 45,000. The instrument comprises 54 spectral orders that cover a continuous wavelength range between 1.42 and 2.42 $\mu$m. To perform the observations, we used a standard ABBA nodding pattern and an exposure time of 100 s (nights 1 and 3) or 130 s (night 2) per frame. There is no particular reason for the longer exposure times during night 2 – all observations were planned with a default exposure time of 100 seconds. However, 130 seconds is still well below limit at which Doppler smearing²²2Doppler smearing occurs when the radial velocity of the planet changes by more than one resolution element during a single exposure. For IGRINS, one resolution element amounts to $c/R\approx 6.7$ km/s, with $c$ the speed of light and $R$ the spectral resolution of the instrument. At mid-transit, it takes $\sim$500 s for the radial velocity of WASP-12b to change by this value. becomes an issue. On each night, we observed the full transit of WASP-121b, which takes just under 3 hours, as well as $\sim$30 minutes of stellar baseline before and after the transit. Including ingress and egress, the transit of WASP-121b occurs roughly between phase angles $\pm 17^{\circ}$ (with $0^{\circ}$ corresponding to mid-transit), which is about a tenth of the planet’s full orbit.

An initial data reduction was carried out with the IGRINS pipeline package³³3Available from github.com/igrins/plp. (PLP, Sim et al. 2014; Lee et al. 2017) by the instrument team. We perform further wavelength calibration and a barycentric-velocity correction using routines from the IGRINS_transit package⁴⁴4Available from github.com/meganmansfield/IGRINS$\_$transit. (Weiner Mansfield & Line 2024), which builds on earlier work from Line et al. (2021). For each observation, IGRINS_transit produces a data cube that contains normalised fluxes as a function of wavelength, spectral order, and time (i.e., orbital phase). After wavelength cropping, each order covers 1848 pixels. We discard orders 1, 23$-$27, 53, and 54, which are situated near the edges of the instrument filters, and which coincide with the strongest telluric absorption lines. This results in SNR values that are much lower compared to other parts of the spectrum.

Next, to correct for throughput variations over time, we put all spectra on a “common blaze” (e.g., Gibson et al. 2020) by dividing the fluxes at each orbital phase by a second-order polynomial (we do this for each order separately). This polynomial is a fit to the fluxes in each pixel divided by their median over time. We also perform a baseline subtraction. To this end, we convolve the fluxes with a Gaussian kernel (FWHM = 80 pixels) and we subtract the result from the original spectrum. This removes broadband features from the data that could impact the cross-correlation analysis. Finally, we cut off 100 pixels at the edges of each order to get rid of boundary effects. Because these pixels have inherently lower SNR values (due to the shape of the blaze function), the loss in planet signal is at most a few percent. At the end, each order comprises 1648 pixels.

The fluxes contained in the data cubes are a combination of the stellar spectrum, telluric absorption (especially in the infrared), the planet signal, and photon noise. When it comes to extracting the planet signal from the data, HRS leverages the fact that stellar and telluric lines remain effectively stationary during the observation, while the spectrum of the planet undergoes a Doppler shift induced by its orbital motion (e.g., Snellen et al. 2010; Brogi et al. 2012; Birkby 2018).

To separate the planet signal from the systematics, we perform a singular value decomposition (SVD) on each spectral order, for each night (e.g., De Kok et al. 2013; Giacobbe et al. 2021). The sum of the first $n$ SVD components (which are all rank-1 matrices) is the best rank-$n$ approximation of the data. Therefore, the components associated with the largest singular values will capture the strong, (quasi-)stationary lines caused by telluric/stellar absorption, while the residuals will contain the much weaker planet spectrum buried in noise.

There is no golden rule prescribing how many SVD components to remove from the data – removing too few will cause systematics to persist, while removing too many will affect the planet signal. Moreover, recent work by Smith et al. (2024) on IGRINS data shows that it may be advantageous to remove a different number of SVD components per order and per night, depending on the observing conditions (humidity, airmass, etc.) and the degree of telluric absorption within an order⁵⁵5For further discussion, see Nugroho et al. (2017), Spring et al. (2022), Cheverall et al. (2023), and Klein et al. (2024).. However, designing such an optimisation process is non-trivial, and HRS studies generally opt to remove the same number of SVD components from each order for a given observing night (e.g., Line et al. 2021; Van Sluijs et al. 2023; Boucher et al. 2023; Brogi et al. 2023; Pelletier et al. 2023; Smith et al. 2024). In this work, we adopt a similar approach. Yet, to remain as agnostic as possible with regards to the telluric correction, we carry out separate analyses for $3\leq n\leq 18$, with $n$ the fixed number of SVD components removed from each of the orders, for each night. The motivation for the lower limit of $n$ = 3 is that, from this value onwards, the $K_{\text{p}}$–$V_{\text{sys}}$ maps of the combined transits (see Fig. 3) show clear detections of CO and H${}_{\text{2}}$O near the expected planet position. The upper limit $n$ = 18 was chosen because this is where the SNR associated with the H${}_{\text{2}}$O detection first drops below 4 (Fig. 4).

Additionally, it should be noted that for each of the three observing nights, there is a partial overlap between the telluric absorption lines and the spectrum of WASP-121b in velocity space (see Fig. 2). That is, $V_{\text{bary}}\approx V_{\text{planet}}(\phi_{i})$ for some phase angle $\phi_{i}$, with $V_{\text{bary}}$ the barycentric velocity and $V_{\text{planet}}$ the planet velocity⁶⁶6In each of our observations, the barycentric velocity changes by less than 1 km/s (see Fig. 2). Conversely, the planet velocity changes by more than 120 km/s during the transit.. Such an overlap is nearly unavoidable in transit observations of short-period planets like WASP-121b, as their radial velocity takes on a wide range of values centered around ($V_{\text{bary}}+V_{\text{sys}}$), with $V_{\text{sys}}$ the systemic velocity of the planet. In this light, it is especially important to perform the data analysis for different values of $n$, as choosing a single number of SVD components could lead to remaining systematics being interpreted as planetary signal.

Once the $n$ largest SVD components are subtracted from the data, we are left with residual fluxes containing the planet spectrum buried in noise. To recover the signals of individual molecules in the atmosphere of WASP-121b, we cross-correlate the residual fluxes with a template spectrum at each orbital phase. This gives rise to a 2D cross-correlation function (CCF):

Our template spectra are based on the SPARC/MITgcm model of WASP-121b from Parmentier et al. (2018). The chemical abundances of this model were computed assuming chemical equilibrium and a solar metallicity (the left column in Fig. 8 shows the equatorial plane of the GCM, including the abundances of CO, H₂O, and OH). To obtain the template spectrum for a chemical species $X$, we post-process the model with gCMCRT (Lee et al. 2022b), a 3D Monte-Carlo raditive transfer code. For each template, we only include the opacities of species $X$ and those of the continuum⁷⁷7We include continuum opacities of H₂, He, and H scattering, collision-induced absorption (CIA) by H₂-H₂ and H₂-He, and bound-free and free-free transitions associated with H^-. For references to the opacity data, see Table 2 in Lee et al. (2022a).. We also account for Doppler-broadening of absorption lines caused by planet rotation. As in Wardenier et al. (2023), all templates are computed at mid-transit, at $R\sim$ 135,000. With regards to Eq. 1, this means that the template is only a function of velocity shift, but not of orbital phase. Once the (rotationally broadened) spectra are obtained, we convolve them with a Gaussian kernel whose FWHM corresponds to the spectral resolution of the instrument ($R\sim$ 45,000). Prior to evaluating the CCF, we subtract the continuum from the templates such that they have a horizontal baseline. The cross-correlation is performed for $v$ between $\pm$300 km/s, with steps of 0.5 km/s.

To further assess the robustness of our analysis, we also perform a test in which we apply the SVD to the data cube from IGRINS$\_$transit directly (i.e., without performing the blaze correction and baseline subtraction described in Section 2.2). However, before computing the CCF, we first inject the template spectrum into the main SVD components and then “recover” it through a second SVD, as is routinely done in retrievals on high-resolution datasets (e.g., Brogi & Line 2019; Line et al. 2021; Brogi et al. 2023; Smith et al. 2024). This gives rise to a new template $\tilde{T}_{n}$, with $n$ denoting the number of removed SVD components. $\tilde{T}_{n}$ is a more realistic representation of the planet signal, as it accounts for the subtle changes that may have been imparted on the spectrum during the telluric correction. Reassuringly, as shown in Fig. 13 in Appendix A, the phase-dependence of the Doppler shifts that we observe for CO and H${}_{\text{2}}$O does not change much with this alternative approach.

The CCF maps of CO and H${}_{\text{2}}$O are shown in the first three rows of Fig. 2, which each row corresponding to one observing night. Besides CO and H${}_{\text{2}}$O, we also searched for the presence of OH, TiO, VO, HCN and H${}_{\text{2}}$S, but we did not obtain significant detections⁸⁸8References to line lists: CO (Li et al. 2015), H₂O (Polyansky et al. 2018), OH (Rothman et al. 2010), TiO (Mckemmish et al. 2019), VO (McKemmish et al. 2016), HCN (Barber et al. 2014) and H₂S (Azzam et al. 2016).. We find very tentative evidence for OH (see Section 3.1 and Fig. 3), but the signal is not strong enough to claim a detection.

As demonstrated in Fig. 2, the planetary trails in the CCF maps are relatively weak. Also, the planet signal on night 2 seems to be strongest when it coincides with the barycentric velocity, which warrants caution. This makes it challenging to assess the phase-dependent behaviour of CO and H${}_{\text{2}}$O on a night-by-night basis. Thus, in order to enhance the signal and average any non-planetary contributions over all nights, we combine the three CCF maps, both in the solar rest frame (fourth row in Fig. 2) and the planet rest frame (bottom row). In the solar rest frame, the planet is visible as a diagonal trail. In the planet rest frame, the planet signal manifests as a vertical trail centered at $v=0$ km/s (modulo any anomalous Doppler shifts). To transform the data to the planet rest frame we use $K_{\text{p}}$ = 217.65 km/s and $V_{\text{sys}}$ = 38.20 km/s, which are the same values as those used in Borsa et al. (2021)⁹⁹9F. Borsa, private communication. Before adding the CCF maps together, we interpolate them onto a common grid with 1,000 points along the phase axis. Upon visual inspection, there appear to be no strong telluric artefacts (i.e., clearly structured features outside the trail) in the combined CCF maps.

Fig. 3 shows the $K_{\text{p}}$–$V_{\text{sys}}$ maps for CO, H${}_{\text{2}}$O, and OH obtained when removing $5$ SVD components from the data. As in Wardenier et al. (2023), we compute the $K_{\text{p}}$–$V_{\text{sys}}$ maps by integrating the values in the CCF map along “orbital trails” of the form $v(\phi)=V_{\text{sys}}+K_{\text{p}}\sin(\phi)$. This gives

We calculate the $K_{\text{p}}$–$V_{\text{sys}}$ maps for $K_{\text{p}}$ between $\pm$350 km/s and $V_{\text{sys}}$ between $\pm$150 km/s. Fig. 3 only shows a subregion of this domain, zooming in on the detection peaks (see Fig. 14 in Appendix A for plots across the whole domain). We use two methods to determine the value of $a$ in Eq. 2. The “box method” (e.g., Line et al. 2021; Nortmann et al. 2024) sets $a$ equal to the standard deviation of all SNR values in a region (a box) far away from the expected planet position. In this work, we use the region $K_{\text{p}}<0$ km/s. The “$\sigma$-clipping method” (e.g., Kasper et al. 2021; Weiner Mansfield et al. 2024) takes the standard deviation of the whole map, but iteratively rejects values that deviate by more than 3$\sigma$ from the median (we use 4 iterations). The $K_{\text{p}}$–$V_{\text{sys}}$ maps in Fig. 3 show the SNR values obtained with the $\sigma$-clipping method.

In each of the individual nights, we detect CO at SNR $\gtrsim$ 4. When combining the observations¹⁰¹⁰10When combining the $K_{\text{p}}$–$V_{\text{sys}}$ maps, we first add the “non-normalised” maps together (i.e., with $a$ = 1 and $b$ = 0), and then compute the values of $a$ and $b$ for the sum (see Eq. 2)., the signal is boosted to SNR $\gtrsim$ 7. Detecting H${}_{\text{2}}$O based on an individual transit appears to be much harder. On night 1, there is a peak at the expected planet position, but its value is not significant. On night 2, we find a relatively strong H${}_{\text{2}}$O peak, but it is offset by $\sim$100 km/s along the $K_{\text{p}}$ axis (possibly due to telluric effects). On night 3, there is a peak with SNR $\sim$ 4 near the expected planet position, but we find a similarly strong feature at a lower $K_{\text{p}}$ value. It is only when combining the observations that an unambiguous signal (SNR $\gtrsim$ 5) emerges near the known $(K_{\text{p}},V_{\text{sys}})$ of the planet. When it comes to OH, we find hints of its presence in the atmosphere of WASP-121b, but the SNR values in Fig. 3 are not significant. We will briefly reflect on this non-detection in Section 6.6.

Fig. 4 shows the detection SNR obtained as a function of the number of SVD components removed from the data. Depending on the method used, we find a maximum SNR of 8–10 for CO (when removing 4 SVD components from all orders for all nights) and a maximum SNR of 5.5–6.5 for H${}_{\text{2}}$O (when removing 8 components). Furthermore, Fig. 5 shows the ($K_{\text{p}}$, $V_{\text{sys}}$) values at which the detection peaks occur for a given number of removed SVD components.

Finally, with regards to the H${}_{\text{2}}$O signal, it should be noted that the $\sim$100 km/s $K_{\text{p}}$ offset on night 2 cannot be physical. This is likely due to some remaining systematics in the data. Therefore, by itself, the H${}_{\text{2}}$O signal from night 2 should be treated with caution. However, we still include night 2 in our combined H${}_{\text{2}}$O observation for several reasons. Firstly, the observation from night 2 still contains valid planet information, as demonstrated by the strong CO detection (for which the peak is not shifted) on the same night. Also, the $K_{\text{p}}$–$V_{\text{sys}}$ map of H${}_{\text{2}}$O does feature a secondary (but less prominent) peak closer near the expected planet position. Furthermore, the observations from night 2 suffered from the lowest humidity and showed the highest SNR (see Table 2 and Fig. 1), so in principle these data should contain the highest-quality signal. In the end, the fact that we combine three independent transits makes the combined observation more robust against noise effects from the individual nights. This is because the complementary planet signals add up, while the impacts of noise and systematics are “diluted” as they are averaged over all nights. Given the currently available data, this is the best we can do to maximise the H${}_{\text{2}}$O signal. Additional transit observations of WASP-121b will be needed to confirm our measurements.

When combining all transits, the detection peaks of CO and H${}_{\text{2}}$O appear offset with respect to the expected planet position in the $K_{\text{p}}$–$V_{\text{sys}}$ maps (Figs. 3, 5). This implies that the Doppler shifts of these species deviate from $v$ = 0 km/s in the planet rest frame (e.g., Wardenier et al. 2021, 2023).

To study the precise behaviour of the planet’s absorption lines with orbital phase, we split the CCF signal into four quarters with bin edges $\phi=\{-13.5^{\circ},-6.75^{\circ},0^{\circ},6.75^{\circ},13.5^{\circ}\}$. Note that we exclude the ingress and egress phases here, to make sure that the planet is fully in front of the star across each bin.

Fig. 6 shows the 1D CCF signals that we obtain for CO and H${}_{\text{2}}$O in the planet rest frame, for each quarter of the transit. Although the planet signals are relatively weak, each CCF (except for H${}_{\text{2}}$O in the first quarter) still features a well-defined peak. As illustrated in Fig. 2, there is no H${}_{\text{2}}$O signal in the first quarter of the transit, which is why we do not recover a peak in the corresponding 1D CCF.

To measure the phase-dependent Doppler shifts of CO and H${}_{\text{2}}$O, we use scipy.optimize.curve$\_$fit to fit a Gaussian to each CCF between $\pm$10 km/s (to add more flexibility to the fit, we also allow for the baseline to be negative). This gave better results than performing the fit across a wider range of velocities. The dashed curves in Fig. 6 (right panels) show our Gaussian fits, truncated by a horizontal line at zero. The dotted vertical lines in the plots denote the peak positions inferred from the fits.

As mentioned previously, we perform the above analysis for removing any number of SVD components between 3 and 18. Fig. 7 shows the Doppler shifts that we obtain for CO and H${}_{\text{2}}$O as a function of orbital phase angle. In spite of the scatter in the values that we measure for a given species at a given phase, the qualitative behaviour of the signals is robust under the removal of different numbers of SVD components from the data, which is reassuring. To calculate the “final” Doppler shifts quoted in Table 3, we take the average Doppler shift over all numbers of removed SVD components, weighted by the detection SNR obtained with the $\sigma$-clipping method¹¹¹¹11It should be noted that the weighted average differs by $<0.1$ km/s from the unweighted average, so both methods give essentially the same Doppler shift. (Fig. 4).

The uncertainties quoted in the covariance matrix of Gaussian fits to the 1D CCFs are too small ($<0.1$ km/s) to constitute realistic error bars. Therefore, we use two alternative methods to estimate the error on the measured Doppler shifts from Table 3. Our first error estimate, $\sigma_{\text{SVD}}$, is the square root of the weighted variance of the Doppler shifts across all numbers of removed SVD components. Again, the weights are given by the SNR values obtained from the $\sigma$-clipping method.

For the second error estimate, $\sigma_{\text{jackknife}}$, we assume that 5 SVD components are removed from the data. We then split each quarter of the transit into 5 bins, which span $\sim$1.3^∘ in orbital phase¹²¹²12Individual frames in the dataset span $0.9^{\circ}$ (nights 1 and 3) and $1.1^{\circ}$ (night 2) in orbital phase, respectively.. Subsequently, we measure the Doppler shift in each quarter by fitting a Gaussian to the sum of each combination of 4 bins (leaving out 1 “sample”). Such a jackknife approach gives rise to 5 separate measurements. The error estimate is given by the standard deviation of these values.

As shown in Table 3, both methods produce similar error estimates overall (within a factor 2). The only exception is the Doppler shift of CO in the fourth quarter of the transit (0.5 vs 1.5 km/s). Upon closer inspection, it turns out that one bin contains a large part of the signal, and when this bin is left out, the Doppler shift measurement changes by a few km/s. This is why the error obtained from the jackknife is larger.

We consider four SPARC/MITgcm models of WASP-121b, which are presented in Fig. 8. The SPARC/-MITgcm was first described by Showman et al. (2009). Since, it has been widely used to study the atmospheric physics and chemistry of (ultra-)hot Jupiters (e.g., Fortney et al. 2010; Showman et al. 2013; Kataria et al. 2013; Parmentier et al. 2018; Steinrueck et al. 2021; Tan et al. 2024).

Our first model of WASP-121b is the drag-free atmosphere from Parmentier et al. (2018). As shown in the first column in Fig. 8, the temperature structure of this model is essentially symmetric, such that the morning and evening limbs have similar chemical compositions. The other three models are based on work by Tan et al. (2024). In contrast to the model from Parmentier et al. (2018), these GCMs also account for heat transport due to H${}_{\text{2}}$ dissociation and recombination (e.g., Bell & Cowan 2018; Komacek & Tan 2018; Tan & Komacek 2019; Roth et al. 2021). The idea behind this mechanism is that H${}_{\text{2}}$ thermally dissociates on the dayside, after which atomic hydrogen gets advected to the nightside, where it recombines into H${}_{\text{2}}$ and releases latent heat. When the atmospheric circulation is predominantly eastward, most of this heat is dumped on the evening limb, resulting in a temperature asymmetry between the eastern and western regions of the atmosphere (second column in Fig. 8).

To explore the effect of atmospheric drag (e.g., Showman et al. 2013; Komacek & Showman 2016; Parmentier & Crossfield 2018), we also consider two models with drag timescales $\tau_{\text{drag}}=10^{6}$ s (weak drag) and $\tau_{\text{drag}}=10^{4}$ s (strong drag), which are shown in the third and fourth column of Fig. 8. The drag timescale encapsulates a variety of physical mechanisms, such as turbulent mixing (Li & Goodman 2010), Lorentz-force braking of ionised winds in the planet’s magnetic field (Perna et al. 2010a), and Ohmic dissipation (Perna et al. 2010b). The drag deposits energy back into the atmosphere. As shown in Fig. 8, increasing the drag strength (and thus lowering $\tau_{\text{drag}}$) slows down winds in the atmosphere¹³¹³13$\tau_{\text{drag}}$ is a tunable parameter in the GCM. It is included as an additional “Rayleigh drag” term $-\vec{v}/\tau_{\text{drag}}$ in the momentum equation solved by the dynamical core. This approach assumes that the drag force is uniform across the whole planet atmosphere. A drag-free scenario corresponds to $\tau_{\text{drag}}\rightarrow\infty$.. This leads to less efficient heat redistribution and a more symmetric temperature structure. Also, the equatorial jet gets suppressed – in the strong-drag model, there is only a day-to-night flow.

Another notable difference between the model from Parmentier et al. (2018) and the models from Tan et al. (2024) is that the latter account for opacities from Fe when evaluating heating and cooling rates, making the thermal inversion on the dayside extend to much lower pressures. We refer to Table 1 in Tan et al. (2024) for the full list of opacities considered in their radiative transfer and to Freedman et al. (2014) for the opacities used by Parmentier et al. (2018). While the extra optical opacity causes the models with H${}_{\text{2}}$ dissociation/recombination to look substantially hotter, they have nearly the same effective temperature ($\sim$2400 K) at the dayside photosphere ($\sim$0.1 bar) as the model from Parmentier et al. (2018).

Table 4 provides a summary of some other important parameters of the four SPARC/MITgcm models. All models were run at a horizontal resolution of C32, which corresponds to roughly 128 cells in longitude and 64 cells in latitude. Before computing phase-dependent spectra of the GCMs with gCMCRT, we bin the outputs down to 32 latitudes and 64 longitudes, as was done in Wardenier et al. (2021, 2023).

Our Doppler-shift measurements of CO and H${}_{\text{2}}$O with IGRINS are not the only phase-resolved transit observations of WASP-121b: Borsa et al. (2021) measured the absorption trail of Fe in the optical with VLT/ESPRESSO (see Fig. 9). We can thus rely on three species to constrain the atmospheric properties of WASP-121b.

We use gCMCRT (Lee et al. 2022b) to compute phase-dependent transmission spectra of the four GCM models across the ESPRESSO (0.38$-$0.79 $\mu$m) and IGRINS (1.42$-$2.42 $\mu$m) bandpasses, accounting for Doppler shifts due to planet rotation and winds. All details regarding the radiative transfer and post-processing can be found in Section 2 in Wardenier et al. (2023), so we just provide a brief summary of the gCMCRT setup below.

Before feeding the GCM outputs into gCMCRT, we map the atmospheric structures onto a 3D grid with altitude (instead of pressure) as a vertical coordinate. We account for the fact that each atmospheric column has a different scale height, which is set by local gravity, temperature, and mean-molecular weight. For each of the four WASP-121b models we simulate 25 spectra (equidistant in orbital phase) between phase angles $\pm$17.6, which also cover the ingress and egress. We assume an edge-on orbit, a semi-major axis of 0.025 AU, a stellar radius of 1.46 $R_{\odot}$, and an orbital period of 1.27 days. At each orbital phase angle, gCMCRT simulates a transmission spectrum by randomly shooting photon packets at the planet limb and evaluating the optical depth encountered by each photon packet. In this calculation, the code accounts for Doppler shifts imparted on the opacities by the radial component of the local wind vector and planet rotation (Wardenier et al. 2021). The transit depth at a given wavelength is then obtained by averaging over all photon packets. We use 10⁵ photon packets per wavelength to accurately model the shapes, depths, and shifts of the spectral lines. Because we do not explicitly treat scattering, the propagation direction of the photon packets does not change throughout the calculation. As in Wardenier et al. (2023), the optical and infrared spectra are computed at resolutions $R=300,000$ and $R=135,000$, respectively. For the radiative transfer, we include the same set of continuum opacities and line species as in Wardenier et al. (2023).

Before computing the CCFs, we convolve the planet spectra and the templates (see Section 3) with a Gaussian kernel corresponding to the resolution of the respective instruments. Because the ESPRESSO observations were performed both in 1-UT ($R\sim$ 138,000) and 4-UT ($R\sim$ 70,000) mode, we convolve the optical spectra and the Fe template to two resolutions, and we calculate the CCF maps for both. Fig. 9 shows the CCF signals of CO, H${}_{\text{2}}$O, and Fe that we obtain for each of the models, with the real data plotted on top. The full CCF maps are shown in Fig. 10.

Indeed, as suggested by the negative $K_{\text{p}}$ offset of CO in the $K_{\text{p}}$–$V_{\text{sys}}$ map (Fig. 3), the absorption lines of CO become increasingly blueshifted during the transit of WASP-121b. This is reminiscent of the CCF signals of Fe and other refractory species that have been reported in the atmosphere of WASP-76b (Ehrenreich et al. 2020; Kesseli & Snellen 2021; Pelletier et al. 2023), whose equilibrium temperature is $\sim$250 K lower than that of WASP-121b. Our observation is also in qualitative agreement with Wardenier et al. (2023), in which we predicted that CO should exhibit similar absorption signatures as refractory species in the atmospheres of UHJs. This is because absorption by CO mainly occurs on the dayside of the planet (see Fig. 7 in Wardenier et al. 2023). At the start of the transit, the CO signal is dominated by the leading (morning) limb, where the redshift associated with planet rotation counteracts the blueshift associated with day-to-night winds. Then, as the dayside of the trailing (evening) limb rotates into view, the Doppler shift should become more negative as planet rotation and day-to-night winds both impart a blueshift to the signal (see also Fig. 2 in Wardenier et al. 2021).

As far as our GCMs of WASP-121b are concerned (Fig. 9), the drag-free and weak-drag models with H${}_{\text{2}}$ dissociation/recombination produce the strongest blueshifts. The fact that both models give rise to essentially the same CCF signal indicates that weak drag has little impact on the wind speeds in the observable part of the atmosphere. Further decreasing $\tau_{\text{drag}}$ does lead to a noticeable change in the line-of-sight velocities, which is why the absorption trail of the strong-drag model is less blueshifted. The jump in the signal from $-$4 km/s to $-$7 km/s after mid-transit is driven by planet rotation (Wardenier et al. 2023).

When H${}_{\text{2}}$ dissociation/recombination is accounted for in our models, the observed CO signal appears to be most commensurate with strong drag in the atmosphere of WASP-121b, although this model still overpredicts the observed blueshift in the second half of the transit by a few km/s. Interestingly, the drag-free model without H${}_{\text{2}}$ dissociation/recombination produces similar Doppler shifts as the strong-drag model with H${}_{\text{2}}$ dissociation/recombination. This result underscores the sensitivity of the wind profile and the resulting Doppler shifts to the treatment of heat transport in the GCM.

Looking at Fig. 9, the Doppler-shift measurement of CO in the second quarter of the transit (just before mid-transit) appears to be somewhat of an outlier compared to the other data points. This is further corroborated by the phase-dependence of the CCF peak strength (see Fig. 6, as well as Fig. 16 in Appendix A). Throughout the transit, the signal strength is essentially constant, except in the second quarter, where the signal is a factor $\sim$2 weaker. This variation is not predicted by any of the models considered in this work (Fig. 10), so this data point should be interpreted with caution.

Contrary to CO and Fe, the H${}_{\text{2}}$O signal of WASP-121b becomes less blueshifted during the transit – in agreement with the positive $K_{\text{p}}$ offset in the combined $K_{\text{p}}$–$V_{\text{sys}}$ map in Fig. 3. As demonstrated in the middle panel of Fig. 9, the only model that can capture the observed phase-dependence of the H${}_{\text{2}}$O trail is the strong-drag model. This is not surprising – in Wardenier et al. (2023), we demonstrated that a decrease in blueshift can only be achieved in a scenario where planet rotation is the dominant contributor to the line-of-sight velocities (e.g., Fig. 8 in Wardenier et al. 2023). That is, winds must be slowed down substantially. Owing to the spatial distribution of H${}_{\text{2}}$O across the atmosphere of an UHJ, most of its signal originates from the nightside, leading to an exact opposite behaviour compared to species that are abundant on the dayside of the planet, such as CO and Fe. At the start of the transit, the trailing (evening) limb dominates the planet signal, while the contribution from the leading (morning limb) is strongest at the end of the transit.

The presence of strong drag in the atmosphere of WASP-121b would be consistent with phase-curve observations from HST, TESS, and JWST in the optical and infrared (Bourrier et al. 2020b; Daylan et al. 2021; Mikal-Evans et al. 2022, 2023), which revealed phase-curves offsets between zero and a few degrees eastward.

The absence of a H${}_{\text{2}}$O signal in the first quarter of the transit remains puzzling. As shown in Fig. 10, none of the GCM models predict zero signal at the start of the transit. However, the models accounting for H${}_{\text{2}}$ dissociation/recombination do show a marginally weaker signal in the first half of the transit than in the second half, so heat transport may form part of the solution. To investigate whether noise or systematics may hamper the detection of H${}_{\text{2}}$O in the first quarter (or whether the lack of signal could be of planetary origin), we perform an injection test whereby we inject a “negative” spectrum¹⁴¹⁴14We multiply each row of the data cube by $1-(-\delta(\lambda))=1+\delta(\lambda)$, with $\delta(\lambda)$ a model of the transit depth of the planet as a function of wavelength $\lambda$. We then perform a regular data analysis as described in Section 2. into the data cube at the planet position. The spectrum is scaled such that it roughly cancels out the CCF peaks of H${}_{\text{2}}$O in the other parts of the transit. This ensures that the injected signal is of the same magnitude as that of the real planet.

Fig. 11 shows the results of the injection test. The CCF of the first quarter features a negative peak at the planet position. At the same time, the other parts of the transit show no clear detections as the (negative) injection cancels out the (positive) planet signal. The fact that there is a negative peak in the first quarter suggests that there is indeed little planet signal to be cancelled out. Furthermore, it demonstrates that the injected signal is detectable given the noise and the systematics of the observation. Thus, had a H${}_{\text{2}}$O signal been present in the first quarter of the transit, we should have been able to measure it.

Although the injection test suggests that the lack of H${}_{\text{2}}$O absorption could be physical, it cannot be fully ruled out that telluric effects play a role (see Section 2). In this light, additional observations with different barycentric-velocity offsets (that do not overlap with the first half of the planet transit) will be beneficial.

If the absence of a H${}_{\text{2}}$O signal in the first quarter of the transit is indeed physical, it could the the result of clouds on the nightside of WASP-121b (e.g., Helling et al. 2021). Because CO and Fe mainly probe the dayside of the atmosphere, their signals should remain largely unaffected (Wardenier et al. 2023). Crucially, however, clouds should only mute the H${}_{\text{2}}$O signal on the trailing limb, but not on the leading limb. To explore this scenario, we recompute the transit spectra of the strong-drag model, but this time we include an optically thick cloud deck between longitudes 0${}^{\circ}<\varphi<$ 180^∘. Our cloud presciption is identical to that in Savel et al. (2022) and Wardenier et al. (2023). That is, we place the cloud deck at all temperatures that lie below the Al₂O₃ condensation curve, while restricting its vertical extent to 10 scale heights (see top panel in Fig. 12). We note that this is a very crude way to account for the impact of clouds, so the model should be seen as a limiting case.

The bottom panel in Fig. 12 shows the resulting H${}_{\text{2}}$O signal. In the first quarter of the transit, H${}_{\text{2}}$O absorption is indeed strongly suppressed (also see Fig. 16 in Appendix A), which is in better agreement with our observations. On the other hand, there is a larger discrepancy between the measured and modelled Doppler shifts in the second quarter of the transit. This is a consequence of the muted trailing-limb contribution. It should be noted, though, that our cloud model is by no means an quantitative fit to the data. The cloud parameters (opacity, vertical extent, etc.) could be optimized to better match the Doppler shifts, but such a retrieval is beyond the scope of this work.

Besides considering clouds, we also investigate whether chemical transport can give rise to a strong asymmetry in the H${}_{\text{2}}$O abundance. The idea is that it takes a nonzero time for H and OH to recombine into H${}_{\text{2}}$O as material is advected from the dayside to the nightside of the planet. Hence, the true H${}_{\text{2}}$O abundances in the limb region do not just depend on pressure and temperature, but also on atmospheric circulation and recombination timescales (e.g., Parmentier et al. 2018). To simulate this scenario, we feed the equatorial abundances and wind speeds from the GCMs into the 2D VULCAN photochemical model (Tsai et al. 2024). For both the weak-drag and the strong-drag model, we find that the nightside H${}_{\text{2}}$O abundances only change substantially (by up to 3 orders of magnitude) in small regions of the atmosphere around 10^-5 bar (see Fig. 15 in Appendix A). Muting the H${}_{\text{2}}$O signal, however, would require such changes across multiple dex in pressure (just like the cloud deck in Fig. 12 has a large vertical extent). Therefore, chemical transport alone cannot explain the muted H${}_{\text{2}}$O signal in the first quarter of the transit.

If “evening clouds” are indeed the most plausible mechanism behind the strongly varying H${}_{\text{2}}$O signal of WASP-121b, one would simultaneously have to explain why clouds do not prevail on the leading (morning) limb of the planet. One scenario could be that condensates form on the trailing limb, and settle gravitationally as they are advected across the nightside. Another driver could be cloud-patchiness caused by global-scale cloud transport (e.g., Komacek et al. 2022).

When it comes to confirming the cloud hypothesis, the JWST/NIRSpec phase curve of WASP-121b (Mikal-Evans et al. 2023) will be a valuable dataset. Cloudy regions on the nightside should emit less flux than predicted by (cloud-free) GCMs, as clouds block radiation from deeper, hotter layers of the atmosphere (e.g., Parmentier & Crossfield 2018; Parmentier et al. 2021; Bell et al. 2024). The JWST dataset also includes a transit spectrum. At low resolution, evening clouds (if present) should also mute H${}_{\text{2}}$O features at the start of the transit.

The right panel in Fig. 9 shows the Fe trails of our GCM models, along with the Doppler shifts of Fe measured by Borsa et al. (2021) with ESPRESSO. Especially for the strong-drag model, there is a considerable difference between the absorption trails at both ESPRESSO resolutions. This is because the trails only show the location of the maximum CCF value as a function of phase. Fig. 10 illustrates what is going on. At $R\sim$ 138,000, the CCF map consists of two modes, each associated with one of the planet limbs – see Nortmann et al. (2024) for an extreme example of this effect. At such high resolution, the two modes are still “spectrally” separated. However, when the spectra are convolved to $R\sim$ 70,000, the modes blend into one, changing the location of the CCF maximum.

For $R\sim$ 70,000 (4-UT mode) there is impressive agreement between the ESPRESSO data and the absorption trail of the strong-drag model, especially since the GCM is by no means an optimised fit to the data. For $R\sim$ 138,000 (1-UT mode), the situation is more complicated. Firstly, the strong-drag model does not manage to reproduce the observed Doppler shifts of Fe in the first half of the transit. Rather, the data show much better agreement with the weak-drag and drag-free models that include H${}_{\text{2}}$ dissociation/recombination. On the contrary, the drag-free model without H${}_{\text{2}}$ dissociation/recombination does not provide a good match – it underpredicts the blueshift of Fe across the entire transit. In the second half of the transit, none of the GCM models produce a perfect fit to the data at $R\sim$ 138,000. While the strong-drag model matches the observed trail better in terms of shape (note the Doppler shifts during the egress phase of the transit), the weak-drag model shows better agreement in terms of overall Doppler shift.

Taken together, the Fe observations are consistent with at least some degree of drag ($10^{4}\lesssim\tau_{\text{drag}}\lesssim\mbox{1}0^{6}$ s) in the atmosphere of the planet, but this drag may not necessarily be uniform. Moreover, as opposed to WASP-76b, the Fe trail of WASP-121b does not provide clear evidence of a strong thermochemical asymmetry between the morning and evening terminator (e.g., Wardenier et al. 2021; Savel et al. 2022). As illustrated in Fig. 8, the limbs of the weak-drag and strong-drag models have very similar temperatures and compositions.

While our GCM models are able to broadly mimic the observed Doppler-shift trends for CO, H${}_{\text{2}}$O, and Fe, there still exists some tension between the measured and modelled absorption signals of WASP-121b.

One issue is that, for each of the four GCM models, the predicted Doppler shifts of Fe and CO are very similar, within $\sim$2 km/s at every phase (see also Wardenier et al. 2023). However, in the first and final quarters of the transit, the observed Doppler shifts of both species differ by 4–5 km/s. These measurements indicate that CO and Fe must (on average) probe different pressures on the dayside, and that the vertical wind shear between these layers is probably stronger than expected (e.g., Miller-Ricci Kempton & Rauscher 2012; Seidel et al. 2021). One limitation of the GCM in this regard, is that its upper boundary lies at 2 $\mu$bar. When performing the radiative transfer at lower pressures with gCMCRT, we simply assume that the wind speeds are the same as in the upper layers of the GCM output. While the bulk of the iron lines must probe pressures $>$2 $\mu$bar (Wardenier et al. 2021, 2023), this could explain part of the discrepancy.

There also exists tension between the signals of H${}_{\text{2}}$O and Fe. To produce the decreasing blueshift observed for H${}_{\text{2}}$O, the GCM requires strong drag (i.e., low wind speeds). Even in the strong-drag scenario, our models struggle to match the $\sim$2 km/s redshift measured in the second half of the transit. At the same time, the strong-drag model fails to reproduce the Fe signal from Borsa et al. (2021) at both VLT resolutions. The fact that the measured Doppler shifts of Fe do not change between both resolutions (see Fig. 9) favours a scenario with weaker drag.

In the light of these observations, it should be noted that our GCMs assume uniform drag across the atmosphere of WASP-121b. In reality, however, $\tau_{\text{drag}}$ will depend on local conditions. For example, in the case of magnetic drag, $\tau_{\text{drag}}$ is a function of the local temperature, magnetic field strength, number density, and ionisation fraction (Rauscher & Menou 2013; Beltz et al. 2023; Soriano-Guerrero et al. 2023). Beltz et al. (2023) showed that accounting for these dependencies with a more sophisticated magneto-hydrodynamical (MHD) prescription can have significant impact on the Doppler shifts of species measured in transmission. Thus, it may not be surprising that a model with uniform drag fails to match the signals of all species at all orbital phases, at different spectral resolutions.

Although we find a hint of OH in the atmosphere of WASP-121b (see Figs. 3, 14), the associated signal (SNR $\sim$ 2.5) is not significant. This is remarkable given that OH was confidently detected (SNR $\sim$ 6.1) in the atmosphere of WASP-76b based on one transit visit with CARMENES ($R\sim 80,000$; Landman et al. 2021), which is on a smaller, 3.5-m telescope. Since the equilibrium temperature of WASP-121b is $\sim$250 K higher than that of WASP-76b, one explanation could be that OH has further dissociated into atomic O and H on the hotter dayside of WASP-121b, such that its features in the transmission spectrum are weaker. A retrieval would allow to shed light on the upper limits of the OH abundance in the region probed by our transit data.

Another aspect to consider is the instrument. Based on two transit visits with IGRINS, Weiner Mansfield et al. (2024) also reported OH in the atmosphere of WASP-76b. However, with SNR values of $\sim$3 and $\sim$4, respectively, their detections are weaker compared to the result from Landman et al. (2021). In this light, observations with other high-resolution spectrographs like ESO-3.6/NIRPS (Wildi et al. 2022) or VLT/CRIRES+ (Dorn et al. 2023) may help to get a better handle on the OH signal of WASP-121b.

Our observations of WASP-121b highlight some challenges associated with performing high-resolution transmission retrievals in the infrared (see also the discussion in Wardenier et al. 2023). In the infrared, the main absorbing species on UHJs are CO and H${}_{\text{2}}$O, which have distinct 3D spatial distributions across the atmosphere. Therefore, the absorption lines of both species will probe different temperatures (CO is more sensitive to the dayside, while H${}_{\text{2}}$O is more sensitive to the nightside), but they will also be subject to different phase-dependent Doppler shifts. This leads to species-dependent peak offsets in the $K_{\text{p}}$–$V_{\text{sys}}$ map (e.g., Fig. 3).

The fact that CO and H${}_{\text{2}}$O probe different parts of the atmosphere suggests that a single temperature (profile) and a single pair of ($\Delta K_{\text{p}}$, $\Delta V_{\text{sys}}$) parameters may not be sufficient to accurately represent the physics underlying the observation in a forward model. To our knowledge, it has not been studied how abundance measurements from high-resolution transmission retrievals are affected when assuming the same Doppler shift and/or temperature for CO and H${}_{\text{2}}$O.

—