© The Authors 2025

1 Introduction

Ultra-hot Jupiters (UHJs) are extreme worlds. Their elevated equilibrium temperature (T_eq > 2200 K) introduces atmospheric dynamics unseen on other planets. Among hot Jupiters, they are uniquely characterised by their dayside temperature profile: They have a pressure layer where the temperature increases with increasing altitude (and decreasing pressure), a feature called a temperature inversion (Fortney et al. 2008; Madhusudhan & Seager 2010; Gandhi & Madhusudhan 2019). Temperature inversions on UHJs are not only predicted by theory but have also been observed (e.g. Evans et al. 2017; Mikal-Evans et al. 2019; Coulombe et al. 2023; Mikal-Evans et al. 2022; Evans-Soma et al. 2025). An atmospheric inversion creates emission lines in the dayside spectra instead of absorption lines seen on planets with a non-inverted profile. Furthermore, some molecular species that are present in great quantity in the colder hot Jupiters are expected to be thermally dissociated in UHJs (Parmentier et al. 2018). Notably, H₂O is prone to thermal dissociation, giving rise to significant amounts of hydroxyl (OH) and atomic hydrogen (H) (Parmentier et al. 2018). On cooler planets, the temperature is not high enough for water to dissociate. Therefore, the relative amounts of OH and H₂O present in the atmosphere can be a good tracer of the temperature (Kitzmann et al. 2018).

Hydroxyl was detected for the first time in the emission spectrum of WASP-33b (T_eq ~ 2800 K) (Nugroho et al. 2021; Collier Cameron et al. 2010). Following this, Finnerty et al. (2023) were able to reproduce their detection and performed retrievals to constrain the atmospheric abundance of OH, H₂O, and other prominent species. They concluded that the OH abundance was two orders of magnitude higher than that of H₂O. Brogi et al. (2023) were able to retrieve OH on the colder UHJ WASP-18 b (T_eq ~ 2500 K) (Hellier et al. 2009), finding an abundance marginally higher than that of H₂O. OH was also detected on WASP-76b (T_eq ~ 2200 K) (West et al. 2016) using transit data from CARMENES¹ (Landman et al. 2021; Gandhi et al. 2024). Mansfield et al. (2024) used IGRINS² transit data of WASP-76 b to detect OH and recover its abundance to be comparable but slightly lower than that of H₂O. More recently, OH was found on the dayside of WASP-121 b (Delrez et al. 2016) using IGRINS (Smith et al. 2024).

With its ability to resolve individual spectral lines, high-resolution spectroscopy is an excellent technique to uncover the dynamics in the atmospheres of planets (e.g. Snellen et al. 2010; Louden & Wheatley 2015; Brogi et al. 2016; Ehrenreich et al. 2020; Prinoth et al. 2022; Nortmann et al. 2025). In particular, Doppler shifts of molecular lines in transmission data can be indicative of day-night circulation (e.g. Ehrenreich et al. 2020; Seidel et al. 2021; Wardenier et al. 2021; Bello-Arufe et al. 2022; Seidel et al. 2023). Meanwhile, different velocity shifts between detected species in dayside emission data can be caused by differences in probed pressure levels or an inhomogeneous spatial distribution of the chemical composition across the photosphere (e.g. Cont et al. 2021; Brogi et al. 2023). Since hot Jupiters are expected to be in synchronous rotation, they have a permanent dayside where the irradiation of the star is always present and a nightside that never sees the light of the host star. This day-to-night contrast in irradiation can lead to temperature differences between the day- and nightsides of sometimes more than 1500 K (Dang et al. 2025), which can cause large longitudinal variations in atmospheric chemistry. This important difference in energy between the two hemispheres and the rotation of the planet are responsible for dynamical effects such as the formation of a super-rotating equatorial jet (e.g. Showman et al. 2009).

Global circulation models (GCMs) were developed to better understand three-dimensional dynamic effects present on hot Jupiters (e.g. Showman et al. 2009; Rauscher & Menou 2012). GCMs simulate the 3D climate of a planet using the primitive equations of meteorology and a radiative transfer prescription. External energy sources, such as the heating from the host star, are considered in the models. This naturally gives rise to atmospheric effects such as winds and longitude-dependent temperature profiles (e.g. Harada et al. 2021; Malsky et al. 2021). The effects of H₂ dissociation and recombination (e.g. Bell & Cowan 2018; Tan & Komacek 2019; Roth et al. 2021) and magnetic drag (e.g. Rauscher & Menou 2013; Beltz et al. 2021, 2022) on the heat transport can also be included to obtain a more realistic view of a planet. GCMs have been successfully used to explain dynamical features that cannot be explained using 1D models (e.g. Wardenier et al. 2023, 2024; Kesseli et al. 2024). With the advent of high-resolution spectroscopy at a high signal-to-noise ratio, the dynamics of the atmosphere have been observed in greater detail than ever before, and the use of GCMs is sometimes necessary to properly interpret the results.

WASP-121 b is a well-studied UHJ (T_eq ~ 2350 K) on a 30-hour orbit (Delrez et al. 2016) of a relatively bright (J=9.6) F6V star. With its large radius (R_p = 1.753 ± 0.036 R_Jup, Bourrier et al. 2020), low gravity, and its extreme irradiation, WASP-121 b is an excellent target for both transmission and emission spectroscopy. In transmission, several chemical species have been discovered using high-resolution spectroscopy, including Fe (Gibson et al. 2020); Cr, V, and Fe⁺ (Ben-Yami et al. 2020); Mg, Ca, and Ni (Hoeijmakers et al. 2020); Na, K, Li, and Ca⁺ (Borsa et al. 2021); Sc⁺ (Merritt et al. 2021); and Co, Sr⁺, Ba⁺, and Mn (Azevedo Silva et al. 2022). Furthermore, H₂O and CO have been detected at infrared wavelengths using IGRINS (Wardenier et al. 2024).

Compared to transmission, high-resolution studies of WASP-121 b in emission are sparser. Hoeijmakers et al. (2024) used the ESPRESSO³ spectrograph to uncover the presence of Ca, V, Cr, Mn, Fe, Co, and Ni. Pelletier et al. (2024) combined the ESPRESSO dataset with infrared observations taken with CRIRES+⁴ K-band data to further obtain detections of H₂O and CO on the dayside of WASP-121 b. Finally, Smith et al. (2024) found H₂O, CO and OH using IGRINS.

In this work, we analyse observations of the dayside emission of WASP-121 b taken with the Near InfraRed Planet Searcher (NIRPS) and find evidence of thermal dissociation and dynamical effects in its atmosphere. In Sect. 2, we present the observations used in this analysis and the NIRPS spectrograph. In Sect. 3, an overview of the data processing is provided. The modelling and cross-correlation results are presented in Sect. 4. The retrieval setup and associated results are shown in Sect. 5. The inferred atmospheric composition and dynamics are discussed in Sect. 6. We conclude this work in Sect. 7.

2 Observations

We observed the dayside of WASP-121 b using the NIRPS (Bouchy et al. 2017; Wildi et al. 2022; Artigau et al. 2024; Bouchy et al. 2025) spectrograph. NIRPS is the new infrared spectrograph on the ESO 3.6-m telescope at La Silla Observatory in Chile. With its wavelength coverage between 0.98 and 1.8 μm, NIRPS covers the Y, J and H bands simultaneously. It is a fibre-fed, adaptive-optics-assisted, and ultra-stable spectrograph, making it an excellent tool for high-precision radial velocity measurements and atmospheric characterisation of exoplanets. Two fibres are available for observations: the High-Accuracy (HA) and High-Efficiency (HE) modes. The HA mode uses a 0.4″ fibre with a median spectral resolution of 88 000. The HE mode uses a bigger 0.9″ fibre, at the cost of a lower median spectral resolution of 75 200 (Bouchy et al. 2025). All WASP-121 b observations were taken with the HE mode, with the exception that the 2023 March 5 observations which were taken with the HA fibre.

From radial velocity measurements, NIRPS was evaluated to be extremely stable to below 1.5 m s⁻¹ (Bouchy et al. 2025). This is significantly lower than the required precision for atmospheric characterisation of close-in giant planets, where the planet’s radial velocity generally vary in the order of tens to hundreds of kilometres per second within a given time series of a few to several hours.

NIRPS complements the optical light instrument High Accuracy Radial velocity Planet Searcher (HARPS), installed on the same telescope (Pepe et al. 2002; Mayor et al. 2003). While both instruments can be used simultaneously during an observation, with HARPS adding its 0.38–0.69 μm wavelength range to the NIRPS infrared coverage, we focus the analysis presented in this article on the NIRPS data.

Along with precise radial velocity measurements (Bouchy et al. 2025; Suárez Mascareño et al. 2025), NIRPS has so far been used to characterise the atmospheric physics and chemistry of several hot gas giants. Helium was found on WASP-69 b (Allart et al. 2025) and H₂O was detected in transmission on WASP-127 b (Bouchy et al. 2025). Furthermore, two transits of WASP-189 b were observed with NIRPS and HARPS simultaneously. Although Fe was detected in the HARPS dataset, it was not found in the NIRPS dataset. This non-detection is attributed to the strong opacity of H⁻ in the infrared (Vaulato et al. 2025).

We observed WASP-121 b for a total of 30 hours spread across six nights using NIRPS as part of the guaranteed time observations (Program IDs 60.A-9109 and 112.25P3.001, PI: Bouchy). The dataset is composed of three pre-eclipse observations taken on the nights of 2023 March 5, 2023 December 13 and 2024 February 2. The other three observations covered the post-eclipse phases, taken on the nights of 2023 December 7, 12 and 21 (Fig. 1). The average S/N per pixel across the full wavelength domain was between 18 and 28 for all observations. The exposure times were 300 s with the exception of 2023 March 5 where the exposure times were set to 400 s due to worse seeing conditions. The observation sequences were between 4 and 5 hours per night, with approximately 50 exposures per sequence. We note that the first pre-eclipse observation was done during the last commissioning (Program IDs 60.A-9109) in 2023 March 5. Although this dataset was taken early in the life of NIRPS, its quality is comparable to that of later datasets.

3 Data processing

NIRPS data are pre-processed by two pipelines: the online NIRPS-DRS running in real time at the telescope, adapted from the ESPRESSO pipeline (Pepe et al. 2021) with telluric correction adapted from Allart et al. (2022) and A PipelinE to Reduce Observations (APERO; Cook et al. 2022) running in Montréal, Canada, originally designed for the SPIRou spectrograph (Donati et al. 2020). Both pipelines use similar pre-processing steps including detector effect removal, wavelength calibration, etc. They ultimately produce equivalent data products, namely the spectra per order, hereafter called 2D spectra. One main difference between the two pipelines is the telluric-removal process. We refer the reader to Bouchy et al. (2025) for in-depth analysis of the pipelines.

For our analysis, we use reduced data from the APERO pipeline version 0.7.290. We opted for the 2D telluricuncorrected spectra and perform the analysis on an order-by-order basis. As we do not use the telluric-corrected spectra, the difference in pipelines should not significantly affect our results. To verify this, we create cross-correlation maps of the detected species using the NIRPS-DRS spectra and obtain the same detections with similar significances (Fig. A.1). Furthermore, other analyses using NIRPS have shown little differences in the results when using either pipeline (Bouchy et al. 2025; Allart et al. 2025). We also ran the cross-correlation analysis with the APERO telluric-corrected spectra, and we recovered consistent detections, suggesting that our detrending approach is sufficient for telluric line removal.

After pre-processing, we pass the data through the Principal Component Analysis (PCA)-based reduction pipeline described by Pelletier et al. (2021) in order to correct telluric and stellar contamination (see also Pelletier et al. 2023; Bazinet et al. 2024). This method exploits the fact that the spectral lines in the planetary spectrum of WASP-121 b undergo a considerable shift during the observations due to the change in the line-of-sight velocity of the planet (e.g. Birkby 2018). On the other hand, the velocity shifts of telluric and stellar lines do not change significantly during the same observations. Therefore, the features caused by the Earth’s atmosphere or the host star should generally stay on the same pixel on the detector. To keep the planetary signal and remove unwanted contributions, it is thus sufficient to remove features that are constant in wavelength throughout an observation.

In practice, our approach can be summarised in a handful of steps applied to each spectral order of each of the six observing sequences. We started by correcting all pixels that deviate more than 5σ from their temporal mean, calculated separately for each observing sequence. If three or less flagged pixels are next to each other, we corrected them with a linear interpolation using neighbouring pixels. If the group of flagged pixels is bigger than three pixels, they were masked altogether. We also masked regions with high telluric contamination. The wavelengths between 1.3450 and 1.4465 μm were completely masked due to this range being mostly contaminated by atmospheric absorption caused by water vapour. We also masked all lines with depths greater than 70% of the continuum flux. Major OH emission lines from the Earth’s atmosphere were masked using the line list reported by Oliva et al. (2015). This step still leaves weaker OH lines in the data, which can be numerous in the near infrared. This can affect the detections of OH in the atmosphere of a planet. However, in all six datasets, the expected velocity curve of WASP-121 b did not fall close to the velocities of the tellurics such that the spectral lines of WASP-121 b and of Earth’s atmosphere do not overlap (Fig. 2). Furthermore, the observations were taken at different times of the year such that the tellurics do not have the same velocities in the barycentric rest frame between observation. We continued by discarding orders where the masking has removed more than 80% of the data. Afterwards, we applied a high-pass filter to bring each exposure to the same continuum level (Gibson et al. 2020). This filter is a two-step process where we first passed a box filter with length 51 pixels, followed by a Gaussian filter with a standard deviation of 100 pixels. We then divided out a second-order fit of the median spectrum from each exposure. We used PCA to reconstruct the observed spectral time series using the first five principal components and dividing this out from the data. This step removes any remaining unwanted telluric and stellar residuals. Here we opted to remove five principal components, which is similar to previous works (e.g. Line et al. 2021; Mansfield et al. 2024; Wardenier et al. 2024; Smith et al. 2024) and is based on what we find best cleans the data of telluric H₂O and OH residuals without unnecessarily removing too much of the planetary signal (e.g. Holmberg & Madhusudhan 2022; Bazinet et al. 2024; Pelletier et al. 2024). Finally, we masked noisy channels that deviate more than 4σ from the standard deviation of their spectral order.

After these steps, the remaining spectra are composed primarily of photon noise, with the planetary signal buried in that intrinsic noise. The cross-correlation method is then used to stack the signal from many individual spectral lines and recover the unique signature of different chemical species in the planetary atmosphere (Snellen et al. 2010; Brogi et al. 2012).

Fig. 1 Summary of the observations of WASP-121 b. The top panel is the air mass as a function of phase for all observations. The bottom panel shows the S/N per pixel throughout the observations. The grey shaded area centred around phase 0.5 represents the phases where WASP-121 b is eclipsed by its host star. Of the six observations, three of them were obtained during the pre-eclipse phases (phases <0.5), while the other three are of the post-eclipse phases (phases >0.5).

Fig. 2 Cross-correlation function of each observation with an OH model (black-and-white colour map) along with the expected WASP-121 b velocities (black dashed line) compared to where tellurics are in velocity space (blue dotted line). The cross-correlation maps are shown in the BERV-corrected frame. In none of our observation nights do the planet velocities overlap with where the telluric are, lowering the risk that any remaining telluric residuals affect our results.

4 Modelling and cross-correlation analysis

4.1 Modelling steps and cross-correlation setup

As a first step to reveal the atmospheric composition of WASP-121 b, we cross-correlated the cleaned NIRPS spectra with generated synthetic atmospheric models. The atmospheric models used in this analysis were computed using the SCARLET framework (Benneke & Seager 2012, 2013; Benneke 2015; Benneke et al. 2019a,b; Pelletier et al. 2021, 2023, 2024; Bazinet et al. 2024). Given a vertical temperature profile and a chemical composition, SCARLET can calculate a synthetic line-by-line high-resolution (R = 250 000) emission spectrum of a planet. One can also let the composition be in chemical equilibrium. SCARLET uses FastChem (Stock et al. 2018, 2022) to calculate the mixing ratios of species in this scenario. In our models, we considered the opacities of H₂O (Polyansky et al. 2018), OH (Rothman et al. 2010), Fe, Mg (Ryabchikova et al. 2015) (Fig. 3), and H⁻ in both bound-free and free-free form (Gray 2021). SCARLET also considers collision-induced absorption from H₂−H₂ and H₂-He (Borysow 2002). Initially, other molecular and elemental species were considered in the models. These species include Fe⁺, V, Cr, Mn, Ni, Ca, Co, Sc, Si, Al, Ti, Ti⁺, TiO, CH₄, CO, CO₂, NH₃, HCN, H₂S, C₂H₂, and SiO (McKemmish et al. 2019; Ryabchikova et al. 2015; Kurucz 2018; Kramida et al. 2018; Rothman et al. 2010; Li et al. 2015; Hargreaves et al. 2020; Yurchenko et al. 2020, 2022; Harris et al. 2006; Barber et al. 2014; Chubb et al. 2020). However, these were never detected in our observations, and they were ultimately not included in the final models.

The model spectra were first convolved with a Gaussian kernel to decrease their resolution to 80 000 to match the NIRPS instrumental resolution. Broadening caused by the rotation of the planet is also considered, with the equatorial rotational velocity of WASP-121 b being 6.99 km s⁻¹ (assuming synchronous rotation, Bourrier et al. 2020). Finally, a box filter is used to account for the planet’s orbital motion relative to our line of sight over the course of a single 300 s or 400 s exposure, adding a small source of smearing of the planetary lines.

We generated several models with different vertical temperature profiles. These models assume chemical equilibrium with solar elemental composition (Asplund et al. 2021). The temperatures profiles are composed of four points with pressure levels of 10², 10⁻¹, 10⁻⁵ and 10⁻⁸ bar. The bottom two points were set to the same temperature, so that the temperature profile between 10² and 10⁻¹ bar is isothermal. The same was done for the points at 10⁻⁵ and 10⁻⁸ bar. We varied these two temperatures from 1500 to 5000 K with a step size of 500 K, excluding temperature profiles that are entirely isothermal. This ensures that models have a range of temperature gradients, from weak to strong and from inverted to non-inverted.

We cross-correlated each exposure with the planetary models described above at radial velocities that range from −400 to 400 km s⁻¹ with a step size of 2 km s⁻¹. We phase folded the result of the cross-correlations to produce K_p − V_sys maps, where K_p and V_sys are respectively the planet semi-amplitude velocity and the star’s systemic velocity. We considered a grid of Keplerians where the K_p varies from −300 to 300 km s⁻¹ at a resolution of 1 km s⁻¹ and V_sys varies from −100 to 100 km s⁻¹ at a resolution of 1 km s⁻¹. We summed the cross-correlation coefficients following these Keplerians to create the maps. To transform the maps to the S/N scaling, we first subtracted the map by the mean of the map excluding a ±10 km s⁻¹ box around the expected velocities to centre the maps around zero. We then divided by the standard deviation of that same region, outside the ±10 km s⁻¹ box. This was done individually for every map.

Fig. 3 Cross-sections of H2O, OH, Fe, and Mg in the wavelength range of NIRPS. The cross-correlations are calculated at a pressure of 1 mbar and a temperature of 2500 K. H2O contains a forest of lines throughout the entire wavelength range shown here. OH is prominent at wavelengths longer than 1.4 μm but still has strong lines around 1 μm. Fe features are mostly present in the shorter wavelengths, with its most important lines below 1.2 μm. Mg has only a few lines overall.

Fig. 4 Cross-correlation signal-to-noise maps of H2O, CO, Fe, and Mg in the atmosphere of WASP-121 b. The dashed vertical and horizontal white lines represent the stellar Vsys (38.12 km s−1, from Line-by-line Artigau et al. 2022) and Kp (216.8 km s−1, calculated from results in Bourrier et al. 2020) of WASP-121 b, respectively. Tentative signals of H2O, OH, Fe and Mg can be seen near the expected orbital position of WASP-121 b.

4.2 Cross-correlation results

Using the cross-correlation method, we find tentative signatures of H₂O, OH, Fe and Mg in the dayside atmosphere of WASP-121 b (Fig. 4). We find that the model that has temperatures of 2000 K at low altitude and 5000 K at high altitude produces the strongest detections. This model has a temperature inversion, consistent with expectations (Evans et al. 2017; Changeat et al. 2024; Hoeijmakers et al. 2024; Pelletier et al. 2024; Smith et al. 2024). H₂O has a signal with an S/N strength of 4.9, OH with an S/N strength of 3.9, Fe with an S/N strength of 4.0, and Mg with an S/N strength of 4.3. While comparing at the preeclipse phases only map and post-eclipse phases map (Fig. C.1), we see that the pre-eclipse phases have better detections overall. We note that the strength of the detections does not tell anything about the abundance of the species. We only use the maps as a first step to know which chemical species are present in the NIRPS data. We include the volume mixing ratios of the four species as parameters in the retrieval (Sect. 5).

While H₂O and OH are expected to be present on the dayside of WASP-121 b and are commonly detected on UHJs with near-infrared spectrographs (Nugroho et al. 2021; Landman et al. 2021; Brogi et al. 2023; Gandhi et al. 2024; Mansfield et al. 2024; Smith et al. 2024), the signals observed for Fe and Mg are more surprising. Fe and other metals are typically detected at optical wavelengths, as their opacities are prevalent in these wavelengths. In this case, however, Fe and Mg lines are observed in the NIRPS bandpass where they still have some spectral lines, detectable with high-resolution spectroscopy (Fig. 3). While the Fe spectral lines in the near-infrared are not as strong as in the optical (Vaulato et al. 2025), the planet-to-star flux contrast in turn is greatly advantageous, vastly increasing the detectability of these species using emission spectroscopy in the near-infrared (as opposed to transmission spectroscopy where the area ratio is relatively similar in the near-infrared compared to the optical). Notably, Mg has relatively few lines in the near-infrared, but they are preferentially on the redder end of the NIRPS wavelength range compared to Fe where the flux contrast is more favourable. The cross-correlation map of Mg has a feature at (K_p, V_sys) = (~258, ~5) km s⁻¹. This feature has a higher significance of 4.8 compared to 4.3 for the feature at (K_p, V_sys) = (~217, ~37) km s⁻¹. These two features are aligned on the diagonal trace created by the cross-correlation of the pre-eclipse phases (Fig. C.1). With few spectral lines, Mg is easily affected by noise and uncertainties, which could explain this abnormality. We cannot exclude that the Mg signal could be an alias of another species, such as Fe. A recent study found that IGRINS is sensitive to refractories (Smith et al. 2024). IGRINS covers the H and K bands and therefore does not cover the short infrared wavelengths such as NIRPS. This is good evidence that the infrared can be harnessed to probe both the volatiles and refractories in UHJ atmospheres. However, we note that the NIRPS wavelength range does not cover the strong CO bands (at 2.3 μm and 4.5 μm), making it difficult to detect CO with NIRPS alone. CO is a major molecular constituent in hot and ultra-hot Jupiter atmospheres. As such, a full view of the volatile content of WASP-121b cannot be accurately inferred from these data alone.

We noticed V_sys shifts between the signals of some species. Namely, H₂O has a visibly redshifted signal compared to all other species (see Sect. 6.2.3).

5 Retrieval analysis

5.1 Retrieval setup

With the tentative detections of molecular and elemental species in the atmosphere of WASP-121 b, we proceed with a retrieval to infer the physical properties of this planet. For this, we used SCARLET coupled with the Bayesian-inference code emcee (Foreman-Mackey et al. 2013).

We used the retrieval to fit the volume mixing ratios of the chemical species found from the cross-correlation analysis. Free constant-with-altitude abundances are fitted for all chemical species. Since WASP-121 b is an UHJ, H₂O is expected to be severely depleted in the upper atmosphere (Parmentier et al. 2018; Pelletier et al. 2024). As a product of H₂O dissociation, OH is also expected to vary considerably as a function of pressure. A constant-with-altitude prescription is therefore not representative of the actual profile of H₂O and OH. However, the objective of this article is to quantify abundances at photospheric level, so that results obtained at similar pressure levels can be compared. In the terminator regions of the slightly colder UHJ WASP-76 b, Gandhi et al. (2024) found that their OH to H₂O ratio did not change when using either the constant-with-altitude prescription or a profile considering dissociation, which supports the idea that a free retrieval is sufficient for our science objective.

We fit the abundance of H⁻ and e⁻, as these are expected to be important sources of continuum opacity in UHJ atmospheres (Arcangeli et al. 2018; Parmentier & Crossfield 2018; Vaulato et al. 2025). These abundances were then used to calculate the H⁻ bound-free and free-free opacity in the models. We also fit the orbital parameters K_p and V_sys.

We introduced a new parameter, ΔV_sys,H2O, which allows water lines to be Doppler shifted compared to other species in the planet spectrum. This is motivated by the fact that the cross-correlation maps show a redshifted H₂O signal, while the other species fall close to the expected velocities (Fig. 4). Indeed, when performing retrievals without the inclusion of this parameter, retrievals would exhibit one of two behaviours. Sometimes, the retrieval would constrain a V_sys near the expected value, with a bounded abundance obtained for OH, but only an upper limit on H₂O. Other times, the retrieval would constrain a V_sys at the redshifted velocity of the H₂O signal in the cross-correlation map, with bounded abundance constraints for H₂O, but only an upper limit for OH. In both cases, the velocity separation between H₂O and OH prevents the retrieval from well-constraining both species simultaneously. While it is possible that using a 2- or 3 D model that takes into account any potential non-uniform distribution of H₂O and OH on the dayside of WASP-121 b could naturally address this issue, this approach would require prohibitively more computational time and is beyond the scope of this paper. The ΔV_sys,H2O parameter allows for H₂O and OH to be simultaneously constrained within our 1D framework, enabling us to better constrain their relative contributions. In particular, relative velocity offsets between different chemical species were observed in UHJ atmospheres (e.g. Brogi et al. 2023), which may lead to biases if assuming a single underlying radial shift for all of them.

For the temperature profile, we used the prescription described in Pelletier et al. (2021), where we fit ten temperature points evenly spaced in log pressure between 10² and 10⁻⁸ bar. The smoothing parameter was set to σ_smooth = 350 K dex⁻² and the prior used was uniform between 1 and 8500 K for each temperature layer. An overview of the fitted parameters and their priors is shown in Table 2. For the likelihood calculation, as in Pelletier et al. (2021), we injected each generated model (projected in time for a given K_p and V_sys) in a PCA-reconstruction of the data and then passing this through the same detrending steps applied to each spectral order in order to mimic altering effect to the model (Brogi & Line 2019; Gibson et al. 2020).

5.2 Retrieval results

Using a free atmospheric retrieval, we constrain the abundances of H₂O, OH, Fe and Mg, as well as the velocities of each species (Fig. 5). The retrieved volume mixing ratios are log₁₀H₂O = $-{3.49}_{-0.93}^{+0.99}$, log₁₀OH = $-{3.59}_{-1.08}^{+1.01}$, log₁₀Fe = $-{2.91}_{-1.36}^{+1.10}$, and log₁₀Mg = $-{2.75}_{-1.37}^{+1.09}$. The uncertainties on the retrieved absolute abundances are relatively large due to the strong correlation between abundance and the temperature-pressure profile. However, relative abundances largely unaffected by this degeneracy and are therefore less dependent on model assumptions (e.g. Brogi & Line 2019; Line et al. 2021; Gibson et al. 2022; Gandhi et al. 2023; Maguire et al. 2023). We retrieve a volume mixing ratio of log₁₀H⁻ < −6.06 (3 σ limit), while log₁₀e⁻ is not well constrained, with the posterior spanning the entire prior space.

The Doppler shifts of all species are retrieved with a K_p value of ${217.21}_{-0.62}^{+0.63}$ km s⁻¹, which is consistent with the known semi-amplitude velocity of WASP-121 b considering the planet’s rotation (216.8 ± 4.5 km s⁻¹, Bourrier et al. 2020). We comment further on this result in Sect. 6.2.2.

The retrieval finds a systemic velocity for the combined OH, Fe and Mg signal of 38.04 ± 0.53 km s⁻¹, consistent with the systemic velocity retrieved by the Line-by-line technique using the NIRPS observations presented in this paper as well as three transmission dataset (38.12 km s⁻¹, Artigau et al. 2022). Our retrieved value is also consistent with previous WASP-121 b results (Brown et al. 2018; Bourrier et al. 2020). We quantify the shift of the H₂O velocity distribution using the retrieval and measure a constrained value of ${4.79}_{-0.97}^{+0.93}$ km s⁻¹. The retrieval therefore finds the V_sys of the H₂O signal to be at ${42.83}_{-0.79}^{+0.80}$ km s⁻¹. We attribute this anomalous shift in the water lines to the 3D nature of the atmosphere of WASP-121 b. We discuss potential mechanisms causing this shift in Sect. 6.2.3.

We retrieve a vertical temperature profile with an inversion layer, as predicted from the cross-correlation detections. This inversion is consistent with modelling and dayside observations of WASP-121 b (Evans et al. 2017; Changeat et al. 2024; Hoeijmakers et al. 2024; Pelletier et al. 2024; Smith et al. 2024).

6 Discussion

6.1 Water dissociation

At the extreme temperatures prevailing on the daysides of UHJs, some of the atmospheric molecular constituents are expected to thermally dissociate. This is the case for H₂O, which dissociates into elemental H and O (Parmentier et al. 2018). These elements recombine to produce OH which, when present in an appreciable amount in the atmosphere, is detectable with high-resolution spectroscopy (Nugroho et al. 2021; Landman et al. 2021; Finnerty et al. 2023; Brogi et al. 2023; Mansfield et al. 2024; Smith et al. 2024). The reaction chain is as follows: $$\begin{array}{rl}& {\mathrm{H}}_2\mathrm{O}\rightleftarrows 2\mathrm{H}+\mathrm{O}\\ & \mathrm{H}+\mathrm{H}+\mathrm{M}\to {\mathrm{H}}_2+\mathrm{M}\\ & \mathrm{O}+{\mathrm{H}}_2\to \mathrm{H}+\mathrm{O}\mathrm{H},\end{array}$$

where M is any chemical species (Parmentier et al. 2018). OH is involved in the recombination of molecular species such as TiO. However, these recombination reactions are significantly slower than the recombination of OH and are not expected to remove a substantial amount of OH in the atmosphere of UHJs (Parmentier et al. 2018).

Another source of OH in an atmosphere is the photolysis of H₂O by UV photons from the host star (Liang et al. 2003). Disequilibrium processes such as photodissociation are believed to be more prominent on colder hot Jupiters and become less important as temperature increases, with UHJ atmospheres being well approximated by thermochemical processes alone (Roudier et al. 2021; Moses et al. 2011). A recent study by Baeyens et al. (2024) looked at the photodissociation of numerous molecular species, including H₂O, in the atmosphere of the UHJ WASP-76 b. They concluded that the production of OH caused by the photolysis of H₂O becomes significant at pressures lower than 10⁻⁵ bar (0.01 mbar). Given that emission and transmission spectroscopy usually probes between the 10⁰ to 10⁻³ bar pressure levels (e.g. Smith et al. 2024; Gandhi et al. 2024), which is significantly higher pressures than the 10⁻⁵ bar balance point, it is safe to assume that the OH and H₂O probed in our observations are mostly affected by thermal dissociation rather than photodissociation. Theoretical work on other hot Jupiters has found similar results, where disequilibrium chemistry is only significant at the top of the atmosphere (e.g. Molaverdikhani et al. 2019; Moses et al. 2011; Moses 2014). WASP-121 b orbits a star similar to that of WASP-76 b; however, it has a shorter orbital period. As such, its equilibrium temperature (T_eq ~ 2350 K) is about 150 K hotter than WASP-76 b (T_eq ~ 2200 K). We should thus expect more thermal dissociation on WASP-121 b compared to WASP-76 b, if probing at the same pressure level. This would make it even more difficult for OH from photodissociation to overcome its thermally dissociated counterpart.

We attempt to quantify the dissociation of H₂O by calculating the ratio of abundance of OH and H₂O. However, we note that while this approach may seem intuitive, there are some caveats to consider. First, the abundances of H₂O and OH in the atmosphere of WASP-121b are expected to decrease at sub-millibar pressures (Figure 6, left panel). Such a gradient in chemical abundance would result in a decrease in line strengths (Parmentier et al. 2018), which our constant-with-altitude abundance profile models cannot reproduce and may fit for lower abundances to compensate. Secondly, the average pressures probed by H₂O and OH lines are not identical (Figure 6, right panel), which means that their ratio may not accurately represent a given pressure level. Nevertheless, with these caveats in mind, our retrieval shows a ratio of log₁₀(OH/H₂O) = −0.15 ± 0.20. If we assume that there is a one-to-one relationship between the dissociated H₂O and the OH production, this result indicates that 41 ± 11% of the H₂O is dissociated⁵.

We started by comparing the retrieved abundances to prediction from chemical equilibrium models. We generated several chemical equilibrium models by assuming a temperature profile sampled from the retrievals. Models have a 5 times solar metallicity, following the results from Smith et al. (2024), who calculated the metallicity using molecular species such as H₂O and OH. With these models, we can find the mean and uncertainties on the abundance of H₂O and OH at every pressure level (Fig. 6). We further calculated the contribution functions for an atmosphere corresponding to the retrieval median fit. These contribution functions show that both H₂O and OH probe close to the mbar level (Fig. 6). The retrieved mixing ratio of H₂O is highly consistent with the chemical equilibrium prediction at 1 mbar. However, the H₂O abundance decreases rapidly at pressures just below the mbar level. The abundance at 1 mbar is around four orders of magnitude greater than at 0.1 mbar. The retrieved OH abundance is consistent within 1 σ with the prediction from chemical equilibrium at 1 mbar. From these chemical models, we calculated the predicted OH/H₂O ratio at each pressure layer (Fig. 6). This ratio varies considerably in the mid-atmosphere, where log₁₀(OH/H₂O) ~−2 at 100 mbar and log₁₀(OH/H₂O) ~4 at 0.01 mbar. The predicted log ratio at 1 mbar is ~0. This is highly consistent with the retrieved ratio of −0.15 ± 0.20. This shows that chemical equilibrium models are good predictors of thermal water dissociation, in agreement with modelling work that predicts UHJs should have atmospheres dominated by thermochemical processes (Roudier et al. 2021; Moses et al. 2011).

However, the above discussions assume that WASP-121 b’s dayside composition is homogeneous throughout the dayside, which might not be realistic. The cross-correlation map shows different dynamics for OH and H₂O. We could therefore be retrieving the H₂O in a region of the atmosphere where there is, for example, a higher abundance than at the location where OH is retrieved. This could also be said for the altitude. However, the contribution functions show that they mostly probe the same pressure levels, with OH probing only slightly lower pressures (Fig. 6).

Gandhi et al. (2024) retrieved the ratio of OH over H₂O in the terminator of the colder UHJ WASP-76b using combined CARMENES and Hubble WFC3⁶ data. Contrary to our results on WASP-121 b, they found that the H₂O is mostly dissociated with an abundance ratio of log₁₀(OH/H₂O) = 0.7 ± 0.3. While it is true that thermal dissociation should be overall more prevalent on hotter planets, the pressure levels probed is also an important factor (Fig. 7). Dissociation is more prominent at lower pressures (or higher altitudes). Since transmission spectroscopy generally probes lower pressures than emission spectroscopy, it is not necessarily surprising that dissociation may appear more pronounced relative to emission spectroscopy. However, the transmission data analysed by Gandhi et al. (2024) probe an average pressure level of 1.5 mbar while the OH detected from our dayside observations originates from around 1 mbar and extends at even lower pressures, contradicting the usual assumption that transmission probes at higher altitude. This discrepancy could be caused by a difference in surface gravity or another such planetary parameter between the two planets. Mansfield et al. (2024) were also able to find the abundance of H₂O and OH in the terminator of WASP-76 b using IGRINS data. They found that the H₂O is dominant by about half to one order of magnitude, depending on the vertical abundance parameterisation used. While this appears to be inconsistent with the results of Gandhi et al. (2024), it is unclear whether the discrepancy between these works may simply be due to their different modelling approaches.

We also compare our results with other UHJ daysides where H₂O and OH were found simultaneously using high-resolution data (Fig. 7). The hottest UHJ on which OH has been detected on its dayside is WASP-33b (Nugroho et al. 2021). There, OH was found to be two orders of magnitude more abundant than H₂O (Finnerty et al. 2023). On the colder UHJ WASP-18b, Brogi et al. (2023) retrieved a constrained H₂O abundance but a lower limit for the OH abundance. They exclude (at 2 σ) scenarios where H₂O is three orders of magnitude more abundant than OH. Their median fit implies that log₁₀(OH/H₂O) ~0.4. The coldest UHJ where the dayside OH and H₂O abundances were retrieved is WASP-121 b, with the results shown in this work. Taken together, these five data points show a possible trend, indicating that warmer planets have more OH, and therefore more thermal dissociation. This is logical since dissociation is directly linked to temperature (Parmentier et al. 2018). However, we note that this trend is expected to fall off at extremely high temperatures. At a certain point, OH is expected to dissociate as well, forming elemental O and H (e.g. Landman et al. 2021; Brogi et al. 2023; Gandhi et al. 2024). The planets follow closely the trend lines from chemical equilibrium models (Fig. 7). The observations follow the trends from equilibrium models at pressures between 0.1 mbar and 0.01 mbar. This does not mean that this is the pressure that is probed by emission spectroscopy. These models assume a uniform temperature profile equal to the equilibrium temperature. On a hot Jupiter, the dayside temperature is higher than the equilibrium temperature. Therefore, the models do not reflect the actual temperature profile in the atmosphere of the planets. Nevertheless, the trend lines can still be useful to compare with observations from a qualitative point of view. Another factor affecting the observed thermal dissociation is the surface gravity of the planet (Parmentier et al. 2018). A difference in gravity will change the pressure level of the photosphere. For a low gravity planet such as WASP-121 b, the photosphere is expected to be at lower pressure, where thermal dissociation of H₂O is more prominent, although this is not what is portrayed in Fig. 7. WASP-121 b has lower gravity than WASP-33b and considerably lower gravity than WASP-18 b. Nevertheless, WASP-33 b, WASP-18 b and WASP-121 b all seem to trace a similar trend line, indicating that the surface gravity might not be such an important factor when it comes to the observations of dissociation in UHJ daysides. However, as all of these measurements were obtained from independent analyses of data sets taken with different instruments, we caution that this mask underlying trends in the observations. This is evidenced by the seemingly different OH/H₂O ratio inferred for WASP-76b (Mansfield et al. 2024; Gandhi et al. 2024). Further uniform analyses with surveys from the same instrument reduced by the same pipeline may shed more insights into the degree of thermal dissociation in different ultra-hot Jupiter atmospheric conditions.

Overall, OH is still a relatively newly discovered molecule in the atmospheres of UHJs, with only a handful of planets in which its abundance relative to H₂O has successfully been retrieved (Finnerty et al. 2023; Brogi et al. 2023; Mansfield et al. 2024; Gandhi et al. 2024). There is still much to learn regarding the trends seen in the thermal dissociation of molecules in UHJ atmospheres. As more OH is detected in other planets with high-resolution infrared spectrograph such as NIRPS, a clearer picture of how dissociation plays a role in shaping the population of UHJ atmospheres will emerge.

Fig. 5 Results from the atmospheric retrieval of WASP-121 b. The left corner plot shows the posteriors of the retrieved parameters, with the marginalised 1D posteriors on the diagonal and the 2D posteriors on the off-diagonals. The contours on the 2D posteriors represent the 1, 2, and 3 σ uncertainties. The first four parameters are the base ten logarithm of volume mixing ratios of the four detected species: H2O, OH, Fe, and Mg. The 1D posteriors for the retrieved abundances show relatively broad distributions. However, 2D posteriors between abundances depict a strong positive correlation, indicating that the abundance ratios are better constrained. logH− and loge− are the volume mixing ratio of the H− radical and of electrons, respectively, used as broadband opacities. Only a weak upper bound is obtained on logH−, while the posterior abundance of loge− spans the full prior. The grey lines on the Kp and Vsys parameters are the expected orbital velocities of the planet (Bourrier et al. 2020; Brown et al. 2018). Both orbital parameters are consistent with expectations. ΔVsys,H2O is the additional Vsys shift added to H2O. The posterior is constrained to a value above 0, indicating that the retrieval also finds the shift seen in the cross-correlation maps (Fig. 4). The top right plot is the retrieved vertical temperature profile, with the shaded areas representing the 1 and 2 σ uncertainties. The profile has a temperature inversion, which is expected for UHJs.

Fig. 6 Comparison between retrieved molecular abundances and equilibrium predictions. We created models in chemical equilibrium using FastChem (Stock et al. 2022, 2018) with a five times solar metallicity composition and where the vertical temperature profile is sampled from the retrieval. In the left panel, the curves and the shaded areas represent the median and 1 σ uncertainties of the mixing ratios of H2O and OH from the equilibrium models. The points are the retrieved abundance for H2O and OH, with their y height at the approximate pressure that they probe. In the middle, the line and the shaded regions are the median and 1 σ uncertainties on the (base ten) log ratio of OH/H2O calculated from the equilibrium models. The point represent the log10(OH/H2O) calculated from the retrieval along with the associated 1 σ uncertainties. The point is at 1 mbar, which is the approximate pressure level where the observations probe. In the right panel, the contribution curves depict where in the atmosphere the signals of H2O (blue) and OH (maroon) are probed. Curves are normalised so that their maximum is equal to one.

Fig. 7 Ratio of the abundance of OH over the abundance of H2O as a function of equilibrium temperature from high-resolution retrievals. The orange points are from the dayside emission of the planets WASP-121 b (this work), WASP-18 b (left orange point, Brogi et al. 2023), and WASP-33 b (right orange point, Finnerty et al. 2023). The green points are obtained from the transit of WASP-76b (Mansfield et al. 2024; Gandhi et al. 2024). Points with error bars are from papers where the value and uncertainty on log10(OH/H2O) are reported (Gandhi et al. 2024, this work). For points without error bars, we take the median of the posterior of OH and H2O to calculate the ratio. The area of the points is proportional to the surface gravity of the planets. The lines are the prediction from chemical equilibrium models computed using FastChem (Stock et al. 2022, 2018) at different pressure levels. The models used are atmospheric models of WASP-121 b, where we set the metallicity and C/O to a solar value and with a uniform vertical temperature profile at the corresponding equilibrium temperature.

Fig. 8 Contour plot of predicted Kp−Vsys maps by a WASP-121 b GCM. The model used is a strong drag model (τdrag = 104 s) based on the work of Tan et al. (2024). The black dashed lines at (ΔKp, ΔVsys) = (~−5, 0) km s−1 represent the planet position considering the rotation of the planet and assuming that the observed signal is uniform throughout the visible hemisphere. The grey line denotes the expected Kp without the rotation of the planet. The blue and maroon crosses are the retrieved velocities of H2O and OH, respectively, along with their 1 σ uncertainties. To get the contour lines, we first cross-correlate noiseless GCMs with molecular opacity templates, followed by a phase folding in the Kp−Vsys space. The contours represent an elevation of 0.99, 0.9 and 0.7 times the map maximum. The left panel is the resulting map from taking both the pre- and post-eclipse phases. The middle and right plot consider only the pre-eclipse and post-eclipse, respectively.

6.2 Velocity traces

The K_p−V_sys maps show signals that are not at the expected planetary velocities (Fig. 4). Firstly, the retrieved K_p for all species is smaller than the expected planetary K_p. Secondly, the H₂O signal is noticeably redshifted compared to other detected species. In this subsection, we investigate these shifts using GCMs. We conclude that the K_p shift is conceptually predicted by the planet’s rotation, meanwhile, the V_sys relative shift cannot be fully explained by GCMs.

6.2.1 Global circulation models

To further interpret the observed Doppler shifts, we compare our results to the Doppler shifts predicted by four 3D GCMs of WASP-121 b. The first GCM is the drag-free model presented in Parmentier et al. (2018), while the other three models are based on recent work by Tan et al. (2024). The latter includes a dragfree model, a model with weak drag (τ_drag = 10⁶ s) and a model with strong drag (τ_drag = 10⁴ s), where τ_drag is the uniform drag timescale of the atmosphere. The models from Tan et al. (2024) account for additional heat transport between the dayside and the nightside due to hydrogen dissociation/recombination (see also Bell & Cowan 2018; Tan & Komacek 2019; Roth et al. 2021), which impacts the 3D temperature structure of the planet. For a more in-depth discussion of the GCMs we refer to Wardenier et al. (2024), who used the same models to interpret transmission observations of WASP-121 b performed with IGRINS.

To obtain phase-dependent emission spectra for each of the four GCMs, we post-processed the models with gCMCRT (Lee et al. 2022a,b), a 3D Monte Carlo radiative transfer code. gCMCRT is able to account for the Doppler shifts imparted on the opacities in each atmospheric cell by considering the line-of-sight velocities due to the planet’s wind profile and its rotation. For details regarding the inner workings of the code, we refer to Lee et al. (2022b) and Wardenier et al. (2025). We calculated the emission spectra for each of the models during pre-eclipse (phases 0.29–0.46) and post-eclipse (phases 0.54–0.71), with intervals of 10 degrees in orbital phase, at a spectral resolution of 135 000. The radiative transfer was performed in the IGRINS bandpass and considers the same line/continuum species as those used in Wardenier et al. (2023, 2024). This includes H₂O and OH, which are relevant for our study. Although IGRINS (1.42–2.42 μm) and NIRPS (0.98–1.80 μm) have different bandpasses, we expect our models to be equally representative for NIRPS. This is because both bandpasses cover a dense forest of H₂O lines which should, on average, probe similar pressures. Furthermore, the cross-sections of OH are strongest in the wavelength region beyond 1.4 μm, where the instruments overlap (see Fig. 3).

Once all emission spectra are obtained, we computed the K_p−V_sys maps associated with H₂O and OH in pre-/post-eclipse by cross-correlating the spectra with a template spectrum and co-adding the resulting CCFs along the phase axis (similar to Wardenier et al. 2024). The results are shown in Fig. 8. Because gCMCRT ignores the orbital motion of the planet when computing the spectra, the Doppler shifts shown in the plots are inherently due to the 3D structure and dynamics of the atmosphere.

The GCMs predict a decrease in the observed K_p, regardless of the model used. See Sect. 6.2.2 for a detailed discussion of this effect. On the other hand, GCMs are not able to reproduce the magnitude of the V_sys redshift observed in the case of H₂O. However, certain scenarios, such as considering only the preeclipse phases, can lead to a small V_sys shift between H₂O and OH. Possible explanations are discussed in Sect. 6.2.3.

6.2.2 Planetary semi-amplitude K_p

We retrieve a planetary semi-amplitude which is consistent with expectations. The retrieved K_p of ${217.21}_{-0.62}^{+0.63}$ km s⁻¹ is well consistent with 216.8 ± 4.5 km s⁻¹ (Bourrier et al. 2020), which takes into consideration the effect of the rotation of WASP-121 b. To calculate the expected K_p, we used the following equation, which assumes a circular orbit: $$K_p=(V_{\mathrm{o}\mathrm{r}\mathrm{b}}-\frac{2}{3}{V}_{\mathrm{r}\mathrm{o}\mathrm{t}})\sin i=\frac{2\pi }{P}(a-\frac{2}{3}{R}_p)\sin i,$$(1)

where V_orb = 2πa/P is the linear velocity of the planet in its orbit and V_rot = 2πR_p/P is the equatorial rotational velocity assuming synchronous rotation. Traditionally, the rotational velocity, which counters the orbital velocity of the planet, is not considered in the calculation of K_p. However, for large planets on small orbits, the effect of the rotational velocity on dayside observations can be in the order of a percent (Hoeijmakers et al. 2024, see their Fig. 8).

The effect of the rotation of a large planet on its radial-velocity semi-amplitude can be conceptually modelled. The equation dictating the radial velocity of a point of the equator at longitude θ (in radians) on a surface of a zero-eccentricity planet at phase ϕ (in radians) is given by (Hoeijmakers et al. 2024) $$v_r(\varphi ,\theta )=V_{\text{orb}}\sin (\varphi )-V_{\text{rot}}\sin (\varphi -\theta ).$$(2)

This equation also assumes that the planet is tidally locked and has an inclination of i ~ 90°⁷. When not considering the rotation (i.e. V_rot = 0), we get that K_p = V_orb, which is what is mostly used in literature. However, when considering a rotation of the planet, the expected semi-amplitude will vary. Setting θ = 0, we see that the measured K_p should be V_orb − V_rot. One can visualise this decrease in velocity by imagining the velocity of the substellar point throughout the orbit. This point will always be in the interior of the orbit. Therefore, it draws a smaller orbit around the star, and hence it will have a smaller K_p compared to the centre of mass of the planet. However, the signal is never emitted from only one point in the atmosphere. We rather get the signal as a weighted average of the visible hemisphere of the planet. If we suppose that the signal was coming equally throughout the dayside, we expect a decrease of $\frac{2}{3}{V}_{\text{rot}}$ when observing the planet close to the eclipse. This assumption is made in the calculation of the expected planetary semi-amplitude and is the origin of the $\frac{2}{3}$ factor in Eq. (1). However, the weighting depends not only on the flux emitted but also the line depth at each point in the visible atmosphere. As high-resolution spectroscopy probes the depth of chemical lines as opposed to the continuum level, it is therefore important to understand the relationship between the vertical temperature profile and the resulting spectra. Although counterintuitive, the region where the temperature is the highest, the hot spot, is not always where most of the signal is coming from. It was shown that an UHJ with an eastern hot spot offset has an eastern dayside that is usually more isothermal compared to the western dayside (van Sluijs et al. 2023). This means that, even with a colder temperature, the difference between the line and continuum in the western dayside is greater, and the signal is therefore coming mainly from that region (van Sluijs et al. 2023). We are not suggesting that this is the case for WASP-121 b, just that the measured K_p is not trivial to calculate, since it requires the complete temperature profile of the visible hemisphere for the entire duration of the observation.

Since GCMs consider Doppler shifts induced by the rotational velocity of WASP-121 b, they predict the K_p decrease caused by rotation regardless of the treatment of winds (Appendix Fig. D.1). This effect is also seen for pre-eclipse and post-eclipse phases individually (Fig. 8).

Our retrieved planetary semi-amplitude, ${217.21}_{-0.62}^{+0.63}$ km s⁻¹, is consistent within 2 σ with Hoeijmakers et al. (2024), where they found a value of ~216 km s⁻¹ using optical observations of the dayside of WASP-121 b. Our K_p is also consistent with other emission analysis of WASP-121 b (Smith et al. 2024; Pelletier et al. 2024).

A caveat of the discussion above is the uncertainty on the value of the planet semi-amplitude velocity, K_p. Without considering the velocity from the planet’s rotation, one can calculate K_p using one of two equations (Birkby 2018): $$K_p=V_{\text{orb}}\sin i=\frac{2\pi a}{P}\sin i$$(3)

or $$K_p=\frac{M_{\ast }}{M_p}{K}_{\ast }.$$(4)

These two equations should give the same value for K_p. However, when using the median values reported by Bourrier et al. (2020), Eq. (3) gives K_p = 221.4 ± 4.5 km s⁻¹ (used in e.g. Pelletier et al. 2024; Hoeijmakers et al. 2024), while Eq. (4) gives K_p = ${218}_{-16}^{+17}$ km s⁻¹ (used in e.g. Wardenier et al. 2024; Smith et al. 2024). Although these values are consistent, there is a major difference between the uncertainties. This difference stems from the uncertainties of the underlying observables. Equation (3)’s uncertainty is dominated by the semi-major axis a. In Bourrier et al. (2020), a is calculated using (a/R_*)R_*. Therefore, the underlying uncertainty comes from the value of the stellar radius, which is ~2%. As for Eq. (4), the uncertainty on K_p comes from the uncertainties in the radial velocity measurements. M_*, M_p and K_* have uncertainties of ~5%, which leads to an uncertainty on K_p of ~7.5%, which is more than three times larger than using Eq. (3). In literature, the uncertainty on K_p is rarely considered. However, it is crucial to take it into consideration, since a shift in K_p of a few kilometres per second may be incorrectly explained by atmospheric dynamics even if it is still consistent with the expected value.

6.2.3 The V_sys shifts

The visible redshift of the H₂O signal is detected with a value of ${4.79}_{-0.97}^{+0.93}$ km s⁻¹, which is ~5 σ away from zero. We also observed small V_sys shifts in the Fe and Mg signals (~ ± 1 km s⁻¹). However, we focus our efforts on H₂O and OH, as they are intrinsically linked through thermal dissociation and their relative shift is more significant. We analyse the relative velocity shifts of H₂O and OH rather than their absolute shifts, as the absolute radial velocity of the system is affected by effects such as gravitational redshift, convective blueshift, and systematic instrumental offsets. Such uncertainties can affect the absolute inferred value of V_sys but would not affect the relative shift between H₂O and OH.

When considering both pre- and post-eclipse phases, GCMs cannot reproduce the shift between H₂O and OH (Fig. 8). However, the pre-eclipse phases map does predict a redshift of the H₂O signature. This effect rises from the fact that the distribution of species is not constant throughout the dayside. At the substellar point, H₂O is depleted and OH is enriched (e.g. Wardenier et al. 2023). Meanwhile, the regions closer to the terminators are dominated by H₂O. In the pre-eclipse phases, the eastern terminator is visible, and it has a positive relative speed (moving away from us), explaining the redshift of the H₂O feature. This same reasoning can be applied to the post-eclipse phases, where the western terminator is in view. This will, however, create a blueshifted H₂O signal. We find that, out of the four GCMs considered, the strong drag model has the greatest V_sys shift between molecules in the pre-eclipse phases. The winds that transport heat to the eastern dayside increases the OH abundance in the eastern dayside, which in turn reduces the relative velocity shift between OH and H₂O. Therefore, a scenario with restricted wind circulation (strong drag) should have a stronger abundance contrast in the eastern dayside creating a stronger velocity shift when observing at the pre-eclipse phases. We note that the relative shift seen in the pre-eclipse phases is in the order of 2 km s⁻¹, while we observe a shift of about 5 km s⁻¹. The complex interactions in GCMs are difficult to model accurately. These hard-to-model effects include wind speed. Therefore, the absolute shifts seen in GCMs should be taken with caution, and they should be looked at from a qualitative point of view.

The fact that heterogeneous distributions of species can create velocity shifts can be explained using Eq. (2). A difference in longitude of the signal will not directly create a shift in V_sys, but instead a shift in the phase of the radial velocity curve. When varying θ in Eq. (2), the resulting sine curve will not have the same phase as the phase of the planet around its host star. This makes sense since when the planet is at phase 0, for example, the point of (non-zero) longitude θ will either be trailing or leading the centre of the planet. However, our dataset only covers a part of the full orbit, where the dayside is in view. It happens that the pre- and post-eclipse phases is when the phase shift creates a similar wavelength shift, which can be approximated by a V_sys shift. Therefore, if we assume that the radial velocity curve of the chemical species has the same phase as the orbital movement of the planet and all observations are taken on the dayside, we should expect to see a shift in V_sys, even if, in reality, this shift arises from a difference in the phase of the radial velocity curve.

If the H₂O signal was only coming from the eastern dayside, we would expect a V_sys shift, similar to the observations. This would also imply that the OH signal contribution is symmetrically distributed around the substellar point. However, if the H₂O contribution is too much to the east, we should not be able to detect the H₂O in the post-eclipse phases. We propose that the H₂O signal on the dayside of the planet is a gradient, with the strongest signal on the eastern part of the dayside which decreases as it approaches the western dayside. As such, the observed spectra from all angles would come predominantly from the eastern part of the visible hemisphere. This is consistent with the observation that the pre-eclipse cross-correlation map seem to have a stronger signal than the post-eclipse map (Fig. C.1). Overall, our observations suggest that WASP-121 b is in a strong drag regime where the H₂O signature is weaker on the western dayside. This weakening of the H₂O signal might be caused by several effects. Firstly, it could simply be that the H₂O is less abundant in the western dayside. Secondly, the vertical temperature profile might be close to isothermal in the western dayside. Thirdly, clouds could hide the H₂O features. However, clouds are not expected to be on the daysides of UHJs due to the extreme temperatures (e.g. Helling et al. 2021; Gao et al. 2020), and therefore we favour the two scenarios mentioned above. Because of the correlation between temperature, abundance, and cloud-top pressure, it is impossible to distinguish between these scenarios with high resolution alone.

The strong drag scenario was used to explain observations of other UHJs. WASP-18 b showed a symmetrical temperature profile around the substellar point (Coulombe et al. 2023), which the authors attribute this strong drag to magnetic effects (Rauscher & Menou 2013; Beltz et al. 2021, 2022). UHJs atmospheres contain ions due to the thermal stripping of electrons around metallic elements. The interaction of these ions with the magnetic field of the planet causes drag that is reflected by low wind speeds and poor heat distribution in the atmosphere.

Day-to-night winds could also induce V_sys shifts. These strong winds generate from the substellar point and travel towards the nigthside such that species that are carried by that wind are expected to be redshifted when observing the dayside. The H₂O might be one of these molecular species. We expect H₂O to be close to the terminator regions, where the day-to-night wind is strong. It is predicted that this wind could be as fast as 10 km s⁻¹ at the mbar level, which is consistent with our observations (Kempton & Rauscher 2012). The OH, on the other hand, is expected to be present closer to the substellar point, due to the strong thermal dissociation of H₂O in this region. The air in the substellar point should rise in altitude; therefore, a species in that wind should not be redshifted but blueshifted instead (Costa Silva et al. 2024). This might explain the small blue shift seen in the OH signal; however, this blue shift is not significant enough to draw conclusive results. GCMs naturally include day-to-night winds. However, if our observed redshift of H2O is caused by day-to-night winds, the model is under-predicting its strength. One way to increase the day-to-night wind speeds is to decrease the drag (i.e. increase the drag timescale) in the GCM, but such a model would be at odds with previous observations of WASP-121 b (e.g. Mikal-Evans et al. 2023; Wardenier et al. 2024).

Other works observed shifts between species in the K_p−V_sys maps (e.g. Brogi et al. 2023; Prinoth et al. 2022; Cont et al. 2021). Notably, Brogi et al. (2023) found velocity shifts between H₂O, OH and CO in a single dayside observation of WASP-18 b. They attributed these shifts mainly to atmospheric dynamics. They explained the difference in shift between H₂O and OH by the different distribution of the molecular species caused by the increased dissociation of H₂O into OH in the hot spot. This explanation is somewhat suitable for observations of only post-eclipse phases or only pre-eclipse phases, but not both. As mentioned above, we require a mechanism to weaken the H₂O signal in the western dayside to explain our relative shift. Interestingly, even with the velocity shifts between molecular detections, Brogi et al. (2023) retrieved constrained abundances for all of their detected species. This is possibly because their dataset only had one observation night. With our six nights, covering both the pre-eclipse and post-eclipse phases, the cross-correlation signal is more local and does not have an extended trail like we would expect from having only one night. Therefore, our retrieval with six nights cannot find a common point in velocity space where it can get the contribution of all chemical species with a relative shift present.

In WASP-121 b, Pelletier et al. (2024) found a blueshifted H₂O signal using CRIRES+ pre- and post-eclipse dayside observations. Although less significant than our observed shift, their H₂O signal is offset by ~−2.5 km s⁻¹ relative to other detected species, in disagreement with our redshifted feature. They point out that a possible explanation is that the western dayside is cooler than the eastern dayside. In this case, the H₂O should be underabundant in the eastern dayside due to thermal dissociation. The western dayside would contain a greater concentration of H₂O, which would contribute more to the signal. Meanwhile, Smith et al. (2024) found a redshifted H₂O feature, along with a marginally blueshifted OH signal using pre- and post-eclipse IGRINS observations. Their relative redshift between OH and H₂O is ~2 km s⁻¹, significantly lower than our observed shift. Their hypotheses to explain this shift include day-to-night winds and a signal coming from a region away from the substellar point, similar to our hypotheses.

The contradicting results between Pelletier et al. (2024), Smith et al. (2024), and this work support the hypothesis that noise might play an important role in our cross-correlation. The relatively low S/N in our observations could introduce some additional noise source. The H₂O spectral lines could preferentially be redshifted due to physical, systematic, or instrumental effects.

Fig. 9 Posterior of the ratio of Fe over Mg in the atmosphere of WASP-121 b from retrieval results. The blue point is the median of the blue distribution, with the error bars being the 1 σ uncertainties of that same posterior. The yellow and red areas respectively represent the 1 σ uncertainties of the Fe/Mg for the Sun (Asplund et al. 2021) and WASP-121 (Polanski et al. 2022). The retrieved ratio is consistent with both the solar and stellar value.

6.3 Metallic ratio

Using a free retrieval, we constrain the abundance of Fe and Mg (Fig. 5) to log₁₀Fe = $-{2.91}_{-1.36}^{+1.10}$ and log₁₀Mg = $-{2.75}_{-1.37}^{+1.09}$. Albeit with large uncertainties, these abundances are marginally elevated compared to a solar atmospheric composition. Observations and numerical simulations predict that WASP-121 b’s metallicity is super-solar with a metal enrichment of 10–30 times the solar value (Daylan et al. 2021; Mikal-Evans et al. 2019; Evans et al. 2018; Thorngren & Fortney 2019). We find that an atmospheric model with a metallicity of ~30 has Fe and Mg mixing ratios that are close to the retrieved abundances. Here, the model used has a vertical temperature profile corresponding to the median retrieved profile. However, we note the high uncertainties on the retrieved abundances indicate that a large range of metallicities is consistent with our results.

The calculated abundance ratio between Fe and Mg using the retrieval results gives a value of log₁₀(Fe/Mg) = $-{0.14}_{-0.39}^{+0.36}$ (Fig. 9). This is consistent with both the solar value of −0.09 ± 0.05 (Asplund et al. 2021) and the stellar value of 0.00 ± 0.07 (Polanski et al. 2022). Our retrieval value is also consistent with the Fe/Mg retrieved by Maguire et al. (2023) from three different ESPRESSO transits of WASP-121 b and by Smith et al. (2024) from two IGRINS dayside observations, who both found a log ratio close to 0. Our result is in agreement with planet formation theories that hot Jupiters, and planets in general, are expected to have metallic ratios that are similar to their host star (e.g. Thiabaud et al. 2015; Adibekyan et al. 2021; Chachan et al. 2023). This trend has observed in other planets as well (e.g. Dorn et al. 2017; Pelletier et al. 2023).

Although the retrieval does not constrain the absolute abundances of Fe and Mg well due to degeneracies with both the temperature structure and the continuum set here by the H⁻ bound-free opacity (Gray 2021; Arcangeli et al. 2018; Parmentier & Crossfield 2018; Vaulato et al. 2025), their abundance ratio can still be precisely determined. It has been shown that retrievals give mostly better constraints on elemental ratios rather than absolute abundance (e.g. Brogi & Line 2019; Line et al. 2021; Gandhi et al. 2023; Maguire et al. 2023). Indeed, the 2D distributions (Corner plot of Fig. 5) reveal positive correlations between species. Since the log of a ratio is the subtraction of the log abundances, we can imagine collapsing the 2D distribution from bottom left to top right. When there is a strong positive correlation between parameters, this will remove the axis with the highest variance and will give a better constraint than the individual 1D distribution.

We note that the Fe and Mg signal comes primarily from one night (2024 February 2). Several factors can explain differences in the detection strengths between observing nights. Such factors include atmospheric variability (Ouyang et al. 2023), but also non-planetary effects. Given the relatively low S/N of our dataset, it is virtually impossible to distinguish which effect is primarily responsible for the observed differences. To determine whether this single night substantially biases the posteriors, we ran a retrieval without the night of 2024 February 2. The resulting retrieval is still able to constrain the Fe and Mg. However, the constraints on the abundances are wider, with a metallic ratio of log₁₀(Fe/Mg) = $-{1.02}_{-1.56}^{+0.59}$ being retrieved. Notably, this retrieval prefers less Fe relative to Mg than the retrieval where all six nights are considered, resulting in a lower Fe/Mg (but still consistent within 2 σ).

7 Conclusions

In this work, we have presented a dayside analysis of the UHJ WASP-121 b using near-infrared high-resolution spectra taken with NIRPS. We report detections of H₂O, OH, Fe, and Mg. The simultaneous detections of H₂O and OH is evidence of thermal dissociation in the dayside photosphere of WASP-121 b. This makes WASP-121 b one of the few UHJs where OH was detected. The atmospheric retrieval analysis returns a ratio of log₁₀(OH/H₂O) = −0.15 ± 0.20, indicating that the abundance of OH is similar to that of H₂O in the photosphere. This ratio is in agreement with the expectation from chemical equilibrium models. However, since water dissociation is highly correlated with the pressure, the ratio varies considerably throughout the atmosphere and an uncertainty on the pressure level probed can lead to drastically different conclusions. We also find a preliminary trend that hotter UHJs exhibit more dissociation in their atmosphere than cooler UHJs. This trend closely follows the expectations from chemical equilibrium models. These results indicate that UHJ atmospheres can be well approximated by models using only thermochemical processes.

All the detections have a measured K_p that is fully consistent with the expected planetary semi-amplitude considering the planet’s rotational velocity. We explored two ways of calculating the expected K_p and noticed that their uncertainties are significantly different. We urge caution when using the value of K_p without taking its uncertainty into account.

The cross-correlation maps show a significant redshift of planetary H₂O signature compared to other detected species. Using a retrieval approach, we quantified the shift to ${4.79}_{-0.97}^{+0.93}$ km s⁻¹. We compared our observations to GCM predictions and find that the models cannot fully explain the observed relative shift between molecular species. A possible explanation is that the H₂O signal on the western dayside is weakened by some atmospheric process, such as an under abundance of H₂O or an isothermal temperature profile. Day-to-night winds where different species are trapped in winds of different regimes could be responsible for a strong shift in systemic velocity. Alternatively, systematic or instrumental noise might be responsible for the shift.

We included the abundance of Fe and Mg as free parameters in the retrieval. The posteriors on the abundance of the individual species are broad. However, the ratio of these species is well constrained with a value of log₁₀(Fe/Mg) = $-{0.14}_{-0.39}^{+0.36}$. This is consistent with the solar value and the value of the host star WASP-121 within 1 σ.

With its wide wavelength range and high resolving power, NIRPS is able to detect OH alongside H₂O in the atmosphere of WASP-121 b. Furthermore, it reveals an intriguing velocity shift between molecular species. These results represent one of the first demonstrations of the power of NIRPS for atmospheric characterisation.

Acknowledgements

We express our gratitude to the referee, whose comments significantly helped improve the quality of the paper. LB, RA, BB, JPW, NJC, ÉA, FBa, RD, LMa, AB, CC, L-PC, AD-B, LD, PL, OL, JS-A, PV & TV acknowledge the financial support of the FRQ-NT through the Centre de recherche en astrophysique du Québec as well as the support from the Trottier Family Foundation and the Trottier Institute for Research on Exoplanets. LB acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC). RA acknowledges the Swiss National Science Foundation (SNSF) support under the Post-Doc Mobility grant P500PT_222212 and the support of the Institut Trottier de Recherche sur les Exoplanètes (IREx). This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. TF, XB, XDe, ACar & VY acknowledge funding from the French ANR under contract number ANR24CE493397 (ORVET), and the French National Research Agency in the framework of the Investissements d’Avenir program (ANR-15-IDEX-02), through the funding of the “Origin of Life” project of the Grenoble-Alpes University. KAM acknowledges support from the Swiss National Science Foundation (SNSF) under the Postdoc Mobility grant P500PT_230225. ÉA, FBa, RD, LMa, TA, J-SM, MO, JS-A & PV acknowledges support from Canada Foundation for Innovation (CFI) program, the Université de Montréal and Université Laval, the Canada Economic Development (CED) program and the Ministere of Economy, Innovation and Energy (MEIE). SCB, NCS, ARCS & EC acknowledge the support from FCT – Fundação para a Ciência e a Tecnologia through national funds by these grants: UIDB/04434/2020, UIDP/04434/2020. SCB acknowledges the support from Fundação para a Ciência e Tecnologia (FCT) in the form of a work contract through the Scientific Employment Incentive program with reference 2023.06687.CEECIND and DOI https://doi.org/10.54499/2023.06687.CEECIND/CP2839/CT0002. The Board of Observational and Instrumental Astronomy (NAOS) at the Federal University of Rio Grande do Norte’s research activities are supported by continuous grants from the Brazilian funding agency CNPq. This study was partially funded by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior—Brasil (CAPES) — Finance Code 001 and the CAPES-Print program. BLCM & AMM acknowledge CAPES postdoctoral fellowships. BLCM acknowledges CNPq research fellowships (Grant No. 305804/2022-7). NBC acknowledges support from an NSERC Discovery Grant, a Canada Research Chair, and an Arthur B. McDonald Fellowship, and thanks the Trottier Space Institute for its financial support and dynamic intellectual environment. DBF acknowledges financial support from the Brazilian agency CNPq-PQ (Grant No. 305566/2021-0). Continuous grants from the Brazilian agency CNPq support the STELLAR TEAM of the Federal University of Ceara’s research activities. JRM acknowledges CNPq research fellowships (Grant No. 308928/2019-9). XDu acknowledges the support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement SCORE No 851555) and from the Swiss National Science Foundation under the grant SPECTRE (No 200021_215200). DE acknowledge support from the Swiss National Science Foundation for project 200021_200726. The authors acknowledge the financial support of the SNSF. JIGH, RR, ASM, FGT, JLR, NN & AKS acknowledge financial support from the Spanish Ministry of Science, Innovation and Universities (MICIU) projects PID2020-117493GB-I00 and PID2023-149982NB-I00. ICL acknowledges CNPq research fellowships (Grant No. 313103/2022-4). CMo acknowledges the funding from the Swiss National Science Foundation under grant 200021_204847 “PlanetsInTime”. Co-funded by the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. GAW is supported by a Discovery Grant from the Natural Sciences and Engineering Research Council (NSERC) of Canada. 0 This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project SPICE DUNE, grant agreement No. 947634). This material reflects only the authors’ views and the Commission is not liable for any use that may be made of the information contained therein. ARCS acknowledges the support from Fundaçao para a Ciência e a Tecnologia (FCT) through the fellowship 2021.07856.BD.LD acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) and from the Fonds de recherche du Québec (FRQ) – Secteur Nature et technologies. FG acknowledges support from the Fonds de recherche du Québec (FRQ) – Secteur Nature et technologies under file #350366. H.J.H. acknowledges funding from eSSENCE (grant number eSSENCE@LU 9:3), the Swedish National Research Council (project number 2023-05307), The Crafoord foundation and the Royal Physiographic Society of Lund, through The Fund of the Walter Gyllenberg Foundation. NN acknowledges financial support by Light Bridges S.L, Las Palmas de Gran Canaria. NN acknowledges funding from Light Bridges for the Doctoral Thesis “Habitable Earth-like planets with ESPRESSO and NIRPS”, in cooperation with the Instituto de Astrofísica de Canarias, and the use of Indefeasible Computer Rights (ICR) being commissioned at the ASTRO POC project in the Island of Tenerife, Canary Islands (Spain). The ICR-ASTRONOMY used for his research was provided by Light Bridges in cooperation with Hewlett Packard Enterprise (HPE). CPi acknowledges support from the NSERC Vanier scholarship, and the Trottier Family Foundation. CPi also acknowledges support from the E. Margaret Burbidge Prize Postdoctoral Fellowship from the Brinson Foundation. AKS acknowledges financial support from La Caixa Foundation (ID 100010434) under the grant LCF/BQ/DI23/11990071. TV acknowledges support from the Fonds de recherche du Québec (FRQ) – Secteur Nature et technologies under file #320056.

References

Appendix A Comparison between pre-processing pipelines

NIRPS has two data reduction pipelines: APERO and NIRPS-DRS. Figure A.1 shows the cross-correlation maps of the four detected species, where the only difference in the data reduction steps is the choice of the pre-processing pipeline. The maps are similar. All features, such as the relative shifts, are present regardless of the pipeline used. Minor differences between maps can be attributed by the different approaches to the telluric correction.

Fig. A.1 Comparison of the cross-correlation maps of the two pre-processing pipelines. The top panels use data pre-processed by APERO (Same as Fig. 4). The bottom panels use the NIRPS-DRS pipeline.

Appendix B Principal components removal

A crucial step in our data detrending procedure is the removal of the first few principal components (PCs). The choice of how many PCs should be removed from a given data set has been the subject of several studies (e.g. Cheverall et al. 2023). Typically, a few PCs are needed to remove the dominant contribution from tellurics; however, removing too many components can risk removing the planetary signal of interest. Here, similar as in Holmberg & Madhusudhan (2022), we explore the evolution of the observed planetary signal as a function of the number of PCs removed (Fig. B.1). All models exhibit the same general trend. When the number of components removed is low, the strength is weak. The strength increases as the number of PCs removed increases, up to a certain point where the strength remain somewhat constant. This point varies for one model to another; however, all of them become constant after ~4 PCs.

To determine whether the H₂O shift is dependent on the reduction steps, we also generate K_p-V_sys signal-to-noise maps of H₂O varying the number of PCs removed (Fig. B.2). The position of the H₂O signal stays between V_sys of 42 and 43 km s⁻¹, regardless of the number of PCs removed. For our main analysis, we choose to remove 5 PCs. However, we also run a retrieval where the only difference is that we set the number of components removed to 10. The results are consistent in both cases (Table B.1).

Fig. B.1 Evolution of signals strengths as a function of the number of principal components. Signal strengths are defined as the maximum in the cross-correlation signal-to-noise Kp-Vsys map in a 12 km s−1 by 12 km s−1 region centred on the expected velocity. This is done for an atmospheric model that includes all four species: H2O, OH, Fe and Mg (black solid line) as well as separately for each species (coloured dashed lines).

Fig. B.2 Cross-correlation signal-to-noise maps of H2O where we vary the number of principal components removed (top left number of each panel). The H2O signal position does not change significantly based on the number of components removed.

Appendix C Pre- and post-phases cross-correlation maps

Figure C.1 shows the cross-correlation maps where we only consider the pre-eclipse phases observation or only consider the post-eclipse phases observations. Overall, the pre-eclipse phases seem to depict better detections.

Fig. C.1 Same as Fig. 4 but where the top row is the map only considering the pre-eclipse phases observations and the bottom map is only considering the post-eclipse phases.

Appendix D GCM cross-correlation maps

Figure D.1 shows the cross-correlation maps of the four considered GCMs, along with the detection locations. All models predict the decrease of the expected K_p by the synchronous rotation of the planet.

Fig. D.1 Same as Fig. 8 but for four different GMCs. The top row has models that has no drag and does not consider latent heat release from dissociation and recombination of H2 (Parmentier et al. 2018). The next three rows are maps of models that consider the dissociation and recombination of H2 (Tan et al. 2024). The difference comes in the drag timescale: τdrag = ∞ (no drag), τdrag = 106 s (weak drag), and τdrag = 104 s (strong drag). The left column maps consider both pre- and post-eclipse phases. The middle column includes the maps when considering only the pre-eclipse phases. In the right column are the maps of the post-eclipse phases only.

All Tables

All Figures

Fig. 1 Summary of the observations of WASP-121 b. The top panel is the air mass as a function of phase for all observations. The bottom panel shows the S/N per pixel throughout the observations. The grey shaded area centred around phase 0.5 represents the phases where WASP-121 b is eclipsed by its host star. Of the six observations, three of them were obtained during the pre-eclipse phases (phases <0.5), while the other three are of the post-eclipse phases (phases >0.5).

In the text

Fig. 2 Cross-correlation function of each observation with an OH model (black-and-white colour map) along with the expected WASP-121 b velocities (black dashed line) compared to where tellurics are in velocity space (blue dotted line). The cross-correlation maps are shown in the BERV-corrected frame. In none of our observation nights do the planet velocities overlap with where the telluric are, lowering the risk that any remaining telluric residuals affect our results.

In the text

Fig. 3 Cross-sections of H2O, OH, Fe, and Mg in the wavelength range of NIRPS. The cross-correlations are calculated at a pressure of 1 mbar and a temperature of 2500 K. H2O contains a forest of lines throughout the entire wavelength range shown here. OH is prominent at wavelengths longer than 1.4 μm but still has strong lines around 1 μm. Fe features are mostly present in the shorter wavelengths, with its most important lines below 1.2 μm. Mg has only a few lines overall.

In the text

Fig. 4 Cross-correlation signal-to-noise maps of H2O, CO, Fe, and Mg in the atmosphere of WASP-121 b. The dashed vertical and horizontal white lines represent the stellar Vsys (38.12 km s−1, from Line-by-line Artigau et al. 2022) and Kp (216.8 km s−1, calculated from results in Bourrier et al. 2020) of WASP-121 b, respectively. Tentative signals of H2O, OH, Fe and Mg can be seen near the expected orbital position of WASP-121 b.

In the text

Fig. 5 Results from the atmospheric retrieval of WASP-121 b. The left corner plot shows the posteriors of the retrieved parameters, with the marginalised 1D posteriors on the diagonal and the 2D posteriors on the off-diagonals. The contours on the 2D posteriors represent the 1, 2, and 3 σ uncertainties. The first four parameters are the base ten logarithm of volume mixing ratios of the four detected species: H2O, OH, Fe, and Mg. The 1D posteriors for the retrieved abundances show relatively broad distributions. However, 2D posteriors between abundances depict a strong positive correlation, indicating that the abundance ratios are better constrained. logH− and loge− are the volume mixing ratio of the H− radical and of electrons, respectively, used as broadband opacities. Only a weak upper bound is obtained on logH−, while the posterior abundance of loge− spans the full prior. The grey lines on the Kp and Vsys parameters are the expected orbital velocities of the planet (Bourrier et al. 2020; Brown et al. 2018). Both orbital parameters are consistent with expectations. ΔVsys,H2O is the additional Vsys shift added to H2O. The posterior is constrained to a value above 0, indicating that the retrieval also finds the shift seen in the cross-correlation maps (Fig. 4). The top right plot is the retrieved vertical temperature profile, with the shaded areas representing the 1 and 2 σ uncertainties. The profile has a temperature inversion, which is expected for UHJs.

In the text

Fig. 6 Comparison between retrieved molecular abundances and equilibrium predictions. We created models in chemical equilibrium using FastChem (Stock et al. 2022, 2018) with a five times solar metallicity composition and where the vertical temperature profile is sampled from the retrieval. In the left panel, the curves and the shaded areas represent the median and 1 σ uncertainties of the mixing ratios of H2O and OH from the equilibrium models. The points are the retrieved abundance for H2O and OH, with their y height at the approximate pressure that they probe. In the middle, the line and the shaded regions are the median and 1 σ uncertainties on the (base ten) log ratio of OH/H2O calculated from the equilibrium models. The point represent the log10(OH/H2O) calculated from the retrieval along with the associated 1 σ uncertainties. The point is at 1 mbar, which is the approximate pressure level where the observations probe. In the right panel, the contribution curves depict where in the atmosphere the signals of H2O (blue) and OH (maroon) are probed. Curves are normalised so that their maximum is equal to one.

In the text

Fig. 7 Ratio of the abundance of OH over the abundance of H2O as a function of equilibrium temperature from high-resolution retrievals. The orange points are from the dayside emission of the planets WASP-121 b (this work), WASP-18 b (left orange point, Brogi et al. 2023), and WASP-33 b (right orange point, Finnerty et al. 2023). The green points are obtained from the transit of WASP-76b (Mansfield et al. 2024; Gandhi et al. 2024). Points with error bars are from papers where the value and uncertainty on log10(OH/H2O) are reported (Gandhi et al. 2024, this work). For points without error bars, we take the median of the posterior of OH and H2O to calculate the ratio. The area of the points is proportional to the surface gravity of the planets. The lines are the prediction from chemical equilibrium models computed using FastChem (Stock et al. 2022, 2018) at different pressure levels. The models used are atmospheric models of WASP-121 b, where we set the metallicity and C/O to a solar value and with a uniform vertical temperature profile at the corresponding equilibrium temperature.

In the text

Fig. 8 Contour plot of predicted Kp−Vsys maps by a WASP-121 b GCM. The model used is a strong drag model (τdrag = 104 s) based on the work of Tan et al. (2024). The black dashed lines at (ΔKp, ΔVsys) = (~−5, 0) km s−1 represent the planet position considering the rotation of the planet and assuming that the observed signal is uniform throughout the visible hemisphere. The grey line denotes the expected Kp without the rotation of the planet. The blue and maroon crosses are the retrieved velocities of H2O and OH, respectively, along with their 1 σ uncertainties. To get the contour lines, we first cross-correlate noiseless GCMs with molecular opacity templates, followed by a phase folding in the Kp−Vsys space. The contours represent an elevation of 0.99, 0.9 and 0.7 times the map maximum. The left panel is the resulting map from taking both the pre- and post-eclipse phases. The middle and right plot consider only the pre-eclipse and post-eclipse, respectively.

In the text

Fig. 9 Posterior of the ratio of Fe over Mg in the atmosphere of WASP-121 b from retrieval results. The blue point is the median of the blue distribution, with the error bars being the 1 σ uncertainties of that same posterior. The yellow and red areas respectively represent the 1 σ uncertainties of the Fe/Mg for the Sun (Asplund et al. 2021) and WASP-121 (Polanski et al. 2022). The retrieved ratio is consistent with both the solar and stellar value.

In the text

	Fig. A.1 Comparison of the cross-correlation maps of the two pre-processing pipelines. The top panels use data pre-processed by APERO (Same as Fig. 4). The bottom panels use the NIRPS-DRS pipeline.
In the text

Fig. B.1 Evolution of signals strengths as a function of the number of principal components. Signal strengths are defined as the maximum in the cross-correlation signal-to-noise Kp-Vsys map in a 12 km s−1 by 12 km s−1 region centred on the expected velocity. This is done for an atmospheric model that includes all four species: H2O, OH, Fe and Mg (black solid line) as well as separately for each species (coloured dashed lines).

In the text

	Fig. B.2 Cross-correlation signal-to-noise maps of H2O where we vary the number of principal components removed (top left number of each panel). The H2O signal position does not change significantly based on the number of components removed.
In the text

	Fig. C.1 Same as Fig. 4 but where the top row is the map only considering the pre-eclipse phases observations and the bottom map is only considering the post-eclipse phases.
In the text

Fig. D.1 Same as Fig. 8 but for four different GMCs. The top row has models that has no drag and does not consider latent heat release from dissociation and recombination of H2 (Parmentier et al. 2018). The next three rows are maps of models that consider the dissociation and recombination of H2 (Tan et al. 2024). The difference comes in the drag timescale: τdrag = ∞ (no drag), τdrag = 106 s (weak drag), and τdrag = 104 s (strong drag). The left column maps consider both pre- and post-eclipse phases. The middle column includes the maps when considering only the pre-eclipse phases. In the right column are the maps of the post-eclipse phases only.

In the text