© The Authors 2025

1 Introduction

Ultra-hot Jupiters (UHJs) are a class of highly irradiated gas giant planets with equilibrium temperatures exceeding 2000 K. The daysides of these planets are expected to be free of clouds (e.g. Helling et al. 2021), which makes them favourable targets for emission spectroscopy. The extreme temperatures lead to molecular dissociation and ionisation of atoms (Parmentier et al. 2018; Helling et al. 2019), and the stark day-night temperature gradient creates fast winds (e.g. Showman & Guillot 2002; Rauscher & Menou 2010).

Simulations of UHJ atmospheres predict the formation of zonal equatorial jets in the direction of the planet’s rotation, along with day-to-night winds transporting material across the terminator (e.g. Zhang et al. 2017). The intensity of these winds depends on a variety of factors, such as the equilibrium temperature, rotation period, and frictional drag (Tan & Komacek 2019). High-resolution spectroscopy can be used to study the radial velocities produced by these effects, such as velocity shifts due to day-to-night winds (e.g. Flowers et al. 2019; Seidel et al. 2023) and line broadening due to jets (e.g. Cont et al. 2022; Lesjak et al. 2023).

The ultra-hot Jupiter WASP-189 b (Anderson et al. 2018) is thus far one of the hottest planets with a characterised atmosphere (T_eq = 3353 K, Lendl et al. 2020). It orbits a bright A-type star with a period of 2.7 days, and has been previously observed both from space with the CHaracterising ExOPlanets Satellite (CHEOPS, Lendl et al. 2020) and the Transiting Exoplanet Survey Satellite (TESS, Saha 2023), as well as from the ground using the High Accuracy Radial velocity Planet Searcher (HARPS) and HARPS-N (Yan et al. 2020; Prinoth et al. 2022; Stangret et al. 2022) and GIANO (Yan et al. 2022). Near-ultraviolet observations suggest that the upper atmosphere can reach temperatures of 15 000 K (Sreejith et al. 2023). The stellar and planetary parameters can be found in Table 1. Prinoth et al. (2022) investigated the transmission spectrum in the visible wavelength region, detecting the presence of TiO and various metals like chromium, iron, and magnesium. They reported varying velocity offsets from the expected planetary rest frame for different species, suggesting different three-dimensional distributions caused by thermo-chemical stratification. A follow-up study (Prinoth et al. 2023) revealed that the signal strength of individual species varies over the duration of the transit, with most species showing increased absorption closer to egress. However, some species (e.g. Ni and TiO) were instead found to decrease in strength over time, which could be attributed to ionisation and dissociation on the planet’s dayside.

The emission spectrum of WASP-189 b was studied in detail by Yan et al. (2020), detecting Fe and constraining the temperature-pressure (T – p) profile, and Yan et al. (2022) found emission of CO in the near-infrared (NIR) observations with a significant red-shift of 4.5 km s⁻¹. This radial velocity offset could be induced by winds from the dayside to the nightside. However, a general circulation model of WASP-189 b simulated by Lee et al. (2022) resulted in fast zonal equatorial winds with speeds of 5–6 km s⁻¹ in pressure regions from 10⁻² to 10⁻⁴ bar, and only a weak day-to-night wind. The T – p profile showed a temperature inversion at longitudinally dependent heights ranging from 10⁻⁵ bar on the nightside to 10⁻¹ close to the sub-stellar point.

In this work, we present high-resolution dayside observations of WASP-189 b with the CRyogenic InfraRed Echelle Spectrograph⁺ (CRIRES⁺, Dorn et al. 2014, 2023). This instrument is the upgraded version of the original CRIRES (Kaeufl et al. 2004) and enables high signal-to-noise (S/N) observations with sufficient time resolution and a spectral resolution of R ~ 100 000. We retrieve the atmospheric properties of WASP-189 b using a Bayesian retrieval framework, and analyse the spectral line shape to determine whether the atmospheric dynamics are predominantly shaped by day-to-night winds or an equatorial jet.

In Sect. 2, we describe the observations and our data reduction procedure, and in Sect. 3 we detail the model generation and cross-correlation analysis. We present our retrieval analysis in Sect. 4, determine the jet and day-to-night wind speeds by applying a more sophisticated line profile model in Sect. 5, and conclude in Sect. 6.

2 Observations and data reduction

2.1 Observations

We used the CRIRES⁺ instrument at the 8m-class Unit Telescope 3 of the Very Large Telescope to observe WASP-189 b on 01 July 2022 as part of the guaranteed time observations of the CRIRES⁺ consortium. In order to detect the planetary dayside emission, the observation covered a part of the orbit shortly after the secondary eclipse. We used the K2148 wavelength setting, covering the spectral range from 1972 nm to 2452 nm split over six spectral orders. This wavelength range includes spectral lines of the major species expected in UHJ atmospheres (e.g. CO, H₂O, Fe). We chose an exposure time of 300 s to ensure sufficient S/N per exposure while not smearing the planetary spectral lines over many detector pixels, as this effect has been proven to be impactful on high-resolution observations (Boldt-Christmas et al. 2024). CRIRES⁺ is an echelle spectrograph with a science detector array comprised of 3 chips, so that each of the six spectral orders is split into three parts, resulting in 18 segments.

Nodding was employed during the observation, wherein the target position alternates between two different locations along the slit (A and B) with four consecutive exposures per nodding position. Information of both positions are combined during the extraction process to remove the background sky signal and detector artefacts. We treated the reduced spectra of A and B as two distinct sets of data, and applied the data processing and analysis on each time-series individually before combining the results during the cross-correlation.

To facilitate the wavelength calibration, metrology was employed prior to the science observations to fine-tune the positions of specified emission lines from Kr- and Ne-lamps. We chose a 0.2″ slit width to ensure a high nominal spectral resolution of R~ 100 000. However, in our observation the adaptive optics system performed very well and the point-spread function (PSF) did not illuminate the whole slit, leading to an even higher resolution (Dorn et al. 2023). We measured the width of the PSF to estimate the true resolution of R~ 140 000, following the methodology of Lesjak et al. (2023).

We reduced the raw spectra with the European Southern Observatory (ESO) CRIRES⁺ data reduction pipeline cr2res (version 1.3.0), which includes dark and flat field corrections, bad pixel removal, and a derivation of the wavelength solution¹. We further applied molecfit (Smette et al. 2015) to the time-averaged spectrum of A and B individually to refine the wavelength solutions for both nodding positions. A further correction of the wavelength solution for each individual exposure was not necessary, as the deviations determined via cross-correlating the telluric lines (see Cont et al. 2024) were sufficiently small.

We discarded 40 pixels at the beginning and end of each segment to avoid detector edge effects, and masked regions that were heavily contaminated by telluric absorption (1997 – 2015 nm, 2162 – 2170 nm, 2314 – 2318 nm, and 2450–2453 nm).

2.2 Normalisation

We removed outliers in the time series of spectral bins with a 5σ-clipping, and corrected for the blaze function by following the three-step process of Gibson et al. (2022) to minimise the distortion of the planetary signal. First, the spectra of each segment were divided by the segment-wide median spectrum. Next, the resulting residual spectra were smoothed by applying a median filter of 501 pixels and then a Gaussian filter with a standard deviation of 100 pixels. Finally, each original spectrum was divided by the smoothed residual spectrum of the corresponding segment. This process removed time-dependent variations and brought all spectra on the same blaze function, which was then removed during the following filtering step. We further excluded deep telluric lines by masking the wavelength regions where the flux falls below 40 % of the continuum level.

2.3 SYSREM

We employed SYSREM (Tamuz et al. 2005; Birkby et al. 2013) to remove the blaze function and the telluric and stellar lines in the spectra, as detailed for example in Gibson et al. (2020) and Lesjak et al. (2023). SYSREM utilises a principal component analysis to construct a model of the linear components in wavelength and time, while incorporating the uncertainty of each data point. The observed spectra contain stellar, telluric, and planetary components with distinct Doppler shifts. While stellar and telluric lines remain nearly static in velocity space throughout the observation, the planetary radial velocity changes significantly. A comprehensive description of SYSREM in the context of high-resolution exoplanet observations can be found in Czesla et al. (2024). We followed the method of Gibson et al. (2022) and first divided both the spectra and uncertainties by the time-averaged spectrum of each segment. Then we used SYSREM to iteratively subtract linear trends in time from each spectral pixel, effectively eliminating the (quasi-)static components of the spectra. This process resulted in residuals mainly comprised of noise and the planetary signal, which shifts in wavelength over time and is not modelled by SYSREM. The linear model of each iteration was refined until the average relative change fell below 0.01. The model was then subtracted from the data before beginning with the subsequent iteration. An example of the resulting residual spectra is shown in Fig. 1. The number of SYSREM iterations impacts the recovered signal strength (Boldt-Christmas et al. 2024), and our methodology for determining this number is outlined in Sect. 3.3.

3 Cross-correlation

3.1 Generating planetary model spectra

We calculated synthetic model spectra for the cross-correlation and retrieval analysis using the radiative transfer code petitRADTRANS (Mollière et al. 2019), and assumed an inverted T – p profile characterised by the parametric model from Guillot (2010). For the models used in the cross-correlation analysis, we chose an irradiation temperature T_irr = 2500 K, an internal temperature T_int = 200 K, an infrared opacity κ_IR = 10^−1.7, and a ratio of visible to infrared opacity γ = 10^1.3. The resulting profile closely matches the one retrieved by Yan et al. (2020) for WASP-189 b.

We used the line lists noted in Table 2 to compute individual models for the following species that could be present under the expected conditions in WASP-189 b’s atmosphere: CO, CO₂, Fe, H₂O, FeH, NH₃, H₂S, HCN, and OH. We set the abundance of each species to a fixed value of log₁₀ VMR = −3, and included continuum absorption of H- and collision-induced absorption (CIA) of H₂-H₂ and H₂–He (Borysow 2002; Richard et al. 2012, and the references therein). The resulting spectra were then converted into planet-to-star flux ratios, convolved to the previously determined resolution of CRIRES⁺ (see Sect. 2) using a Gaussian instrumental profile, and continuum normalised.

Fig. 1 Data reduction steps of a representative wavelength range. Top panel: Unprocessed spectra as they are produced by the CRIRES+ pipeline. Middle panel: After normalisation to bring all spectra onto a common blaze function, and masking of deep telluric lines. Bottom panel: After the removal of stellar and telluric lines with SYSREM.

3.2 Cross-correlation of data and model

We performed the computation of the weighted cross-correlation function (CCF) from the residual spectra and the model as described by (Cont et al. 2022): $$\text{CCF}(v,t)={\displaystyle \sum _{i=0}^N}\frac{R_i(t)\cdot M_i(v)}{{\sigma }_i{(t)}^2},$$(1)

where t represents the time index, i the pixel index, R are the residual spectra with uncertainties σ, and M the model spectrum, which was Doppler shifted with velocities v ranging from −500 km s⁻¹ to +500 km s⁻¹. This resulted in a two-dimensional CCF map for each spectral segment. We combined the information from individual segments by computing the mean of these CCF maps, and merged the CCFs of both nodding positions.

We subsequently derived a K_p–v_offset map from the CCFs by exploring a range of orbital semi-amplitude values (K_p) and shifting the CCF into the corresponding planetary frames based on the Doppler velocities: $$v_p=K_p\sin 2\pi \varphi +v_{\text{sys}}-v_{\text{bary}}+v_{\text{offset}}.$$(2)

Here, ϕ denotes the orbital phase, v_sys and v_bary are the systemic and barycentric velocities, and v_offset accounts for a potential deviation from the planetary rest frame. The shifted CCFs were then collapsed into a one-dimensional vector by averaging along the time axis. Stacking these vectors for every trial K_p value resulted in a two-dimensional K_p-v_offset map, which was converted into units of S/N by dividing by the standard deviation of the map in regions far away from the expected signal position (|v| > 50km s⁻¹).

3.3 Selecting the number of SYSREM iterations

When stellar and telluric lines are removed from the data using SYSREM, the number of iterations influences the resulting K_p-v_offset map. Insufficient iterations leave strong residuals that impede the detection of planetary spectral signatures, while an excessive number of iterations risks removing the planetary signal. To objectively determine the optimal number of iterations, we implemented the method proposed by Cheverall et al. (2023, see their Sect. 3.8).

We first applied SYSREM to the normalised data, generating the cross-correlation function CCF_obs. Then we repeated this process with an artificial signal injected at the expected Doppler velocity to produce CCF_inj, from which we derived the differential cross-correlation function ΔCCF = CCF_inj – CCF_obs and calculated the K_p-v_offset map and S/N for each iteration. We chose the SYSREM iteration with the highest S/N in ΔCCF and applied the same number of iterations to the non-injected data. Figure 2 shows the S/N after each iteration for the example of CO, and the results of this process for all species are summarised in Table 2.

3.4 Cross-correlation signals

We detected a cross-correlation signal for CO (S/N = 14.2) and Fe (S/N = 5.8), while the remaining species were not detectable. While we chose the number of SYSREM iterations using the method detailed above, the signals are stable across a wider range of iterations. The model spectra and K_p-v_offset maps for each individual species are shown in Fig. 3 and A.1. For both detected species we find an offset from the expected planetary rest frame velocity (v_offset = 6 km s⁻¹ for CO and v_offset = 7 km s⁻¹ for Fe), while the K_p agrees with the theoretical value within the error bars. This finding is in agreement with the red-shift of ${4.5}_{-3.0}^{+2.2}\text{\hspace{0.17em}km\hspace{0.17em}\hspace{0.17em}}{\text{s}}^{-1}$ measured in previous emission observations of CO (Yan et al. 2022). As detailed in their work, this could be evidence for global day-to-night winds, which have also been proposed by Prinoth et al. (2022) to explain a blue-shifted transmission signal.

Fig. 2 S/N detection strength as a function of SYSREM iterations for the signal of CO, at nodding positions A (blue colour) and B (orange colour). The top panel shows the signal strength for the real data, and the bottom panel shows the strength of an injected signal using the differential ΔCCF according to the method described in Sect. 3.3. The diamond shapes indicate the SYSREM iteration with the strongest detection of the injected signal.

4 Atmospheric retrieval

4.1 Treatment of the models

We performed a retrieval analysis to quantitatively constrain the physical parameters of the planetary atmosphere. To this end, a Markov chain Monte Carlo (MCMC) algorithm sampled a large number of models with varying parameters, which were then compared to the residual spectra. These synthetic model spectra were generated in a similar way as described in Sect. 3.1, and included all of the species that we probed for in the cross-correlation analysis. However, we did not assume isobaric abundances for the retrieval and instead determined pressure-dependent abundances based on equilibrium chemistry calculations using fastChem (Stock et al. 2018). With this approach, the entire chemistry is determined by the atmospheric carbon-to-oxygen (C/O) ratio and metallicity [M/H], and the effect of thermal dissociation of molecules is included. We assumed that the bulk metallicity corresponds to the abundance of most atmospheric elements ([M/H] = [Fe/H] = [O/H] = [N/H] and so on), while the carbon abundance is set by the C/O ratio and oxygen abundance. As demonstrated by Brogi et al. (2023), isobaric abundances are an inadequate description for the extreme environments found in UHJ atmospheres and lead to conclusions inconsistent with the more appropriate equilibrium chemistry approach. In this work, we did not include disequilibrium effects such as photoionisation, which could influence the atmospheric abundance profiles of the highly irradiated UHJs. For the T – p profile, we fixed T_int = 200 K as the observed spectra are not sensitive to the deeper atmospheric regions that are affected by this parameter.

The models were transformed into the same three-dimensional format as the residual spectra (segment × time × wavelength) and shifted into a planetary rest frame determined by the free parameters K_p and v_offset in the retrieval process. Additionally, we accounted for smearing of the signal due to the planet’s change in velocity during each exposure. We calculated radial velocities at the start and end of each exposure (a difference of ~3 km s⁻¹), shifted the model to ten evenly spaced planetary rest frames between these velocities, and computed the smeared model as the mean over the ten sub-exposures.

To address rotational broadening effects, we computed a rotational profile from the equatorial rotation velocity v_eզ, following the method proposed by Yan et al. (2023) using Eq. (3) in Díaz et al. (2011). We assumed a linear limb darkening law with a coefficient ϵ = 1 and an inclination angle of sin i ≈ 1, given the expected tidal locking of the planet. We convolved the model with this profile in logarithmic wavelength space, and finally brought the spectra to the resolution of CRIRES⁺ using a Gaussian instrumental profile.

The application of SYSREM for the removal of stellar and telluric lines introduces distortions in the planetary signal. This effect should be accounted for to accurately match the models to the observed spectrum. While it is possible to directly apply SYSREM to the models to mimic the distortions, this would be computationally expensive and we instead applied an alternative model filtering following Gibson et al. (2022). This approach involves preserving the individual vectors comprising each SYSREM model, and subsequently calculating a filter matrix from these results and the data uncertainties. Such a filter can be pre-computed and then applied to each individual model during the retrieval process. We analyse the performance of this filter and its ability to preserve subtle differences in the line profile in Appendix B.

Fig. 3 Model spectra (top), cross-correlation functions (CCFs, middle), and Kp-voffset maps (bottom) of CO, Fe and H2O. The dashed lines in the middle panels indicate the expected course of the planetary trail (which is located between the two lines). The signals of CO and Fe are visible in the CCFs and Kp-voffset maps, while H2O was not detectable.

4.2 Likelihood function

We calculated the logarithmic likelihood function (Hogg et al. 2010; Yan et al. 2020) to compare each model spectrum to the data: $$\ln (L)=-\frac{1}{2}{\displaystyle \sum _{i,j}}\left(\frac{{\left(R_{i,j}-M_{i,j}\right)}^2}{{\left(\beta {\sigma }_{i,j}\right)}^2}+\ln \left(2\pi {\left(\beta {\sigma }_{i,j}\right)}^2\right)\right),$$(3)

where R_i,j denotes the residual spectrum at pixel i and time j, M_i,j represents the 2D matrix of a filtered model spectrum, σ_i,j are the uncertainties of the residuals, and β is a scaling factor for the σ_i,j. The parameter space was sampled using an MCMC algorithm as implemented in emcee (Foreman-Mackey et al. 2013), with 32 walkers and 15 000 steps each. We discarded the first 5000 steps as burn-in, and assessed the convergence by checking the autocorrelation lengths as well as dividing the chains into multiple sets and ensuring that all sub-sets result in similar posteriors. We chose uniform priors for all free parameters, which we summarised in Table 3.

Fig. 4 Posterior distributions from the retrieval of the atmospheric and orbital parameters. The dashed lines indicate the mean value and the 16th and 84th percentile. The included parameters consists of the three parameters defining the T–p profile (Tirr, κIR, γ), the orbital semi-amplitude of the planet (Kp), the offset from the expected planet rest frame (voffset), a scaling factor to the uncertainties (β), the equatorial rotation velocity (veq) and two parameters determining the equilibrium chemistry (C/O ratio and metallicity).

4.3 Retrieval results

The posterior distributions of the free parameters are shown in Fig. 4 and summarised in Table 3. The scaling factor β, which is applied to the data uncertainties, is close to unity, indicating that the uncertainties determined during the data reduction are accurate. The remaining free parameters determining the T – p profile, chemistry, and velocities are discussed in the following.

4.3.1 Temperature–pressure profile

Observations of the dayside emission are mostly sensitive to the properties of the thermal inversion in the atmosphere, because the emission lines originate from that region. We therefore expect to constrain the pressure range at which the inversion is located and the temperature difference between the deeper, colder layer and the upper, warmer layer. The retrieved values of the three free parameters defining the T – p profile (T_irr, κ_IR, and γ) are correlated with each other. This is due to the fact that various combinations of these parameters can lead to similar temperature inversions.

Larger values of the visible to infrared opacity γ are correlated with lower infrared opacities κ_IR and reduced irradiation temperatures T_irr As γ increases, the inversion layer shifts to higher altitudes and the temperature difference increases. This effect can be mitigated by decreasing κ_IR, which results in a lower inversion altitude. Additionally, reducing T_irr lowers the temperature across the entire upper atmosphere down to pressures of ~1 bar. These adjustments thus achieve a T–p profile closely resembling that associated with lower γ values.

The average retrieved T–p profile, derived from 10 000 randomly drawn samples, is shown in Fig. 5. The profile is well constrained between 10⁻¹ bar and 10⁻³ bar, which is the region probed by our high-resolution observation. However, the profile cannot be constrained at lower pressures because the observed spectra remain insensitive to the temperatures at these altitudes. The same is true for the temperature of the deep atmospheric layers, which is artificially constrained in the retrieval because we fixed T_int to a constant value.

Fig. 5 Retrieved temperature and abundance profiles. Left panel: retrieved T–p profile calculated from 10 000 random samples. The shaded region indicates the 1σ uncertainty interval. Middle and right panel: retrieved abundance profiles of the individual species as determined from the C/O and [M/H] values and the T–p profiles of 1000 random samples. The two species in the middle panel have been detected in the cross-correlation analysis, while the species in the right panel were undetectable.

4.3.2 Chemistry

The chemical composition in our models is determined by the atmospheric C/O ratio and metallicity, and we retrieved $\text{C/O}={0.32}_{-0.14}^{+0.41}$ and $[\text{M}/\text{H}]={1.40}_{-0.60}^{+1.39}$ by assuming equilibrium chemistry. However, the posterior distribution of C/O showed a long tail towards higher values and is fully consistent with the solar C/O ratio. A major reason for the inability to tightly constrain the C/O ratio stems from the non-detection of other oxygen-bearing species besides CO, such as H₂O or OH. The information about the abundances of CO and Fe alone was not sufficient to precisely determine the C/O ratio, especially due to the effect the T–p profile has on the emission spectrum. The high C/O values were generally associated with smaller κ_IR and elevated T_irr, resulting in models with similar strengths of CO and Fe emission lines. In our chemistry model, the abundance of oxygen is linked to that of iron, and therefore our retrieved C/O ratio is only accurate if we assume that [O/H] = [Fe/H].

Although we consistently retrieved super-solar metallicities, the precise value remained uncertain due to a degeneracy with C/O. Higher [M/H] correlated with lower C/O, as these combinations yield similar CO abundances. While the abundance measurement of Fe should theoretically help in determining [M/H], it appears insufficient to break the degeneracy in practice. Due to the large uncertainties, the metallicity is consistent with the solar value within 3σ.

If the metallicity of WASP-189 b is truly higher than the solar metallicity of the host star (Anderson et al. 2018), this would imply a significant enrichment of heavy elements. In combination with a sub-solar to solar C/O, this would hint at an accretion of icy solids outside of the H₂O snowline (Öberg et al. 2011). An alternative pathway to such a combination of C/O and metallicity is an intensive erosion of the planets core and a distribution of the material across the gaseous envelope (Madhusudhan et al. 2017).

We derived the volume mixing ratios (VMRs) of all included species from the retrieved chemical and T–p parameters by averaging over 1000 random samples (see Fig. 5). The log(VMR) of the detected species CO and Fe was constrained within uncertainties of ~ 1 dex in the altitude range of the inversion layer (10⁻¹ – 10⁻³ bar). However, we could not reliably determine their abundances above these pressure levels, because the observed spectrum is insensitive to changes in these regions. Gandhi et al. (2023) retrieved a Fe abundance of $\log \left({\text{VMR}}_{\text{Fe}}\right)=-{4.94}_{-0.26}^{+0.24}$ from high-resolution transit observations, which is consistent with our abundance profile in the upper atmosphere. It is a reasonable assumption that the signal in their work originates from these higher regions, as they additionally retrieve an opacity deck that blocks contributions from below 10⁻³ bar. This could also be the reason for the low metallicity found in their work, which is based on the Fe abundance and falls significantly below our retrieved metallicity value.

Species without detectable cross-correlation signals were included in the retrieval to discourage combinations of C/O and [M/H] that result in abundances which would be detectable with our observation. According to the retrieved equilibrium chemistry, we expect significant amounts of H₂O and OH in the lower atmosphere, with only traces of FeH, NH₃, CO₂, HCN and H₂S. None of these species persist in the upper atmosphere due to thermal dissociation. The retrieved abundance profiles are consistent with the non-detections in cross-correlation, as confirmed by injection tests.

4.3.3 Velocities

The retrieval resulted in a planetary radial velocity semi-amplitude of $K_p={193.4}_{-2.5}^{+2.4}$ km s⁻¹, falling within the 2σ range of the value derived from literature orbital parameters (K_p,lit = 201 ± 4 km s⁻¹). Additionally, we found a significant velocity offset from the expected planetary rest frame of v_offset = 4.7 ± 0.8 km s⁻¹. Both results agree with the signal position found in the cross-correlation analysis, and with the findings of previous high-resolution emission observations of this planet (Yan et al. 2020, 2022).

We assumed that the planetary atmosphere is rigidly rotating and retrieved an equatorial velocity of $v_{\text{eq\hspace{0.17em}}}={4.35}_{-0.50}^{+0.47}$ km s⁻¹. Including a contribution of 3.0 km s⁻¹ from tidally locked rotation, this indicates an excess speed of ${1.35}_{-0.50}^{+0.47}$ km s⁻¹ for the planetary surface rotation pattern. These velocities offer initial insights into the wind dynamics on the dayside of WASP-189 b: a red-shift of the observed signal can be caused by day-to-night winds, while an increased rotational broadening suggests the presence of zonal equatorial jets, which are expected to form on UHJs (e.g. Lee et al. 2022). To investigate which of these processes is more prevalent, we performed an additional retrieval that includes a more sophisticated modelling of the effects of winds on the line profile in the following section.

Fig. 6 Examples of different radial velocity distributions on the planetary disk and resulting velocity profiles, and the result from the retrieval. The example profiles assume a limb darkening coefficient of ϵ = 1. (a) Tidally locked rotation (according to WASP-189 b’s orbital period), (b) Day-to-night wind of vday–night = 5 km s−1. (c) Zonal jet with a speed of vjet = 5 km s−1 and a width of σjet = 0.3 Rp. (d) Velocity distributions according to the best-fitting parameters from the retrieval, and the average retrieved line profile calculated from 1000 random samples. The shaded region indicates the 1σ uncertainty interval.

5 Retrieval of the wind pattern

5.1 Calculating velocity profiles

To retrieve the atmospheric wind pattern, we determined an empirical line profile by dividing the planetary disc into a grid of quadratic cells, with the planet’s diameter corresponding to 2001 cells. Then we assigned a radial velocity and brightness to each cell and computed a velocity profile from these data. The velocities of different regions of the atmosphere are caused by a combination of multiple effects:

• Tidal locking induces a uniform rotation of the planet with a rotation period equal to the orbital period. In the case of WASP-189 b, this leads to a maximal velocity of 3.0 km s⁻¹. The projected radial velocity grows linearly with the distance from the axis of rotation,

• Day-to-night winds flow from the sub-stellar point to the nightside in all directions, leading to an overall red-shift. We assumed a constant wind speed v_day-night across the entire visible hemisphere, and the resulting radial velocities depend solely on the angle to the line of sight.

• Zonal jets typically cause a super-rotation of the regions close to the equator. We assumed that the velocity decreases towards the poles following a Gaussian profile, and parametrised such a jet with its equatorial velocity v_jet and width σ_jet (describing the standard deviation of the Gaussian distribution in units of R_p).

We did not directly include a hot spot in our simulation in order to reduce the computational complexity. Instead, we approximated its effect by modulating the brightness of each cell following a linear limb-darkening law with a coefficient ϵ as a free parameter. In this way, the central region has a larger contribution compared to the terminator at the outer edge of the planetary disk.

The resulting velocity profile was determined by calculating a histogram of the radial velocities while weighting each cell by its brightness. Examples for the velocity distributions and profiles for tidally locked rotation, a day-to-night wind, and a jet are shown in Figs. 6a, b, and c, respectively. To include this framework in the model spectra calculation, we replaced the rotational profile described in Sect. 4.1 with this more complicated velocity profile, and convolved it with the models in logarithmic wavelength space.

Compared to the previously conducted retrieval shown above, we have thus replaced the free parameter v_eզ with v_jet, σ_jet, v_day-night, and ϵ, and assumed the tidally locked rotation to be always present. All other parameters and the included species remained the same.

5.2 Retrieved wind parameters

We found that the observed line profile and velocity offset are best fitted by a model dominated by a fast day-to-night wind. The posterior distributions of the wind parameters are shown in Fig. 7. The velocity of the day-to-night wind was retrieved as $v_{\text{day-night\hspace{0.17em}}}={4.4}_{-2.2}^{+1.8}$ km s⁻¹, and the jet velocity as $v_{\text{jet\hspace{0.17em}}}={1.0}_{-1.8}^{+0.9}\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}$, which is consistent with the over-rotation of the first retrieval. The width of the jet (σ_jet) remained unconstrained, which is likely due to the low velocity of the jet, because the effect of the jet width on the line profile in such a scenario is minor, leading to weak constraints. We were not able to place constrains on the limb darkening coefficient ϵ, likely due to its minute effect on the spectra which was not measurable in our data. The retrieved line profile as calculated from averaging over 1000 random samples is shown in Fig. 6d.

Compared to the retrieval without winds, the velocity offset was significantly reduced to $v_{\text{offset\hspace{0.17em}}}={1.9}_{-1.5}^{+1.6}$ km s⁻¹. This is because a day-to-night wind induces an overall red-shift of the planetary signal, which can account for some of the offset. However, the offset cannot be completely explained by our retrieved wind scenario, because this would require wind speeds upwards of 7 km s⁻¹, which would lead to an extensive broadening that is inconsistent with the observed broadening. Other effects which have not been included in our model could account for the remaining offset from the expected planetary rest frame. The increased flux from a hot spot close to the sub-stellar point might dominate the observed emission signal, and as the hot spot rotates away from the observer following the tidally locked rotation, an overall red-shift of the spectrum could be produced. We assumed that the full dayside is visible during the entire observation, while in reality the planet rotates so that more and more of the nightside comes into view. This effect is not very pronounced in our observation, which is close to the secondary eclipse and only covers a small section of the orbit (0.53 < ϕ < 0.57). According to our simulations, a hot spot alone cannot result in the retrieved line profile, and a day-night wind is required to fit the extended red-shifted wing of the line profile. Nevertheless, the hot spot could introduce a small red-shift as the location of the nightside is associated with negative radial velocities in our viewing geometry. Alternatively, the offset could be due to uncertainties in the ephemerides and system parameters. The propagated error of the orbital period corresponds to a velocity uncertainty of ~1 km s⁻¹, and an error of similar scale could be introduced by the uncertainties in the systemic velocity due to extensive stellar line broadening of the rapidly rotating host star and inaccuracies of the stellar models. Further, a slight eccentricity (0 < e < 0.005) could cause an offset of the correct magnitude while still being consistent with the observed signal within the uncertainties. Additionally, the uncertainty of the planetary radius propagates to the tidally locked rotation velocity, which in turn affects the line broadening. A larger radius would result in stronger broadening from the rotation alone, which in turn would reduce the day-to-night wind speed still consistent with the observed broadening. As the radius uncertainties from the CHEOPS observations of Lendl et al. (2020) are on the order of 1%, we do not expect this effect to significantly influence our results. In contrast to the absolute offset, the relative offset between emission and transmission signal is not affected by these systematic effects. In the future, a joint retrieval of both types of observations could result in a consistent model that can explain all of the observed velocity shifts.

We note that the retrieved posteriors of the day-to-night wind speed are consistent with a wind-free scenario to within ~2σ, and these large uncertainties are mainly caused by the degeneracy between wind speed and velocity offset. In a wind-free scenario, the entire offset of ~5 km s⁻¹ would have to be caused by other effects such as those mentioned above, while in the retrieved best-fitting scenario the remaining offset is reduced to 2 km s⁻¹.

The full corner plot of all free parameters for this retrieval with a dedicated wind line profile is shown in Fig. C.1, and a summary of the priors and retrieved values can be found in Table C.1. We found that the T–p profile parameters, C/O ratio, metallicity, and K_p are consistent within 1 σ with the previous retrieval, which did not include a model of the wind pattern and only accounted for a general line broadening.

While the presence of both jets and day-to-night winds is predicted for atmospheres of hot Jupiters (Zhang et al. 2017), the models of Showman et al. (2013) showed that the equatorial jets can be inhibited due to radiative and frictional damping for planets under extreme stellar irradiation. In this case the day-to-night flows would dominate, which agrees with the best-fitting scenario from our retrieval. The retrieved weak equatorial jet also agrees with the conclusion of Tan & Komacek (2019), who simulated the dependence of the jet strength on the equilibrium temperature and rotation period. They found that higher temperatures lead to weaker jets, because the dissociation of H₂ on the dayside and subsequent recombination on the night-side decreases the day-night temperature gradient. This effect could be one explanation for the apparent slower wind speeds on the ultra-hot WASP-189 b compared to some cooler hot Jupiters, such as WASP-43 b (Lesjak et al. 2023) or WASP-166 b (Seidel et al. 2020). The 3D general circulation models of Lee et al. (2022) predict a fast jet of 5 km s⁻¹ in WASP-189 b’s atmosphere, which is at odds with our result. However, their work did not include the dissociation of H₂, which could lead to an overestimation of the day-night temperature gradient and consequently result in higher jet speeds. Additionally, their model was not able to fully reproduce the lux variation from the CHEOPS phase curve data. While they find that their model agrees with the results of Prinoth et al. (2022), we note that in the transmission spectrum both day-to-night winds and equatorial jets can cause blue-shifts of the detected signals. In the case of a jet, the material at the hotter evening terminator moves towards the observer and contributes more to the signal than the material on the cooler morning terminator. When observing the dayside of a planet, equatorial jets mainly cause a line broadening while the introduced velocity shift is minute. Therefore it is difficult to explain the observed offset in the emission spectrum only with a jet.

The retrieved day-to-night wind speed is in general agreement with the observed blue-shifts in the transmission signal of ~5 km s⁻¹ from previous studies. Prinoth et al. (2022) found different offsets for each of the detected species, and in particular their value for Fe of 3.63 ± 0.43 km s⁻¹ is slightly below our retrieved value but agrees within our uncertainties. We note that our wind properties are mostly determined based on the CO line shape, because the signal of CO was significantly larger than that of Fe.

While our model aims to determine the resulting line profiles of various scenarios of atmospheric dynamics, it is based on some simplifying assumptions that warrant consideration. We assumed uniform day-to-night wind speeds and a fixed viewing geometry, which excluded contributions from the nightside. Additionally, the model did not account for a sub-stellar hot spot, which could increase lux from a localised region of the atmosphere and enhance the ionisation of molecular species. To fully capture the complexities of the planet’s atmospheric circulation, a future study might combine transmission and dayside observations and employ a more complex 3D modelling of the line profile which incorporates these effects.

Fig. 7 Posterior distributions from the retrieval of the atmospheric and orbital parameters. The dashed lines indicate the mean value and the 16th and 84th percentile. The posteriors of the jet width σjet and limb darkening coefficient ϵ fill the entire parameter space and the percentile values of these two parameters are thus not meaningful.

6 Summary

We conducted a comprehensive analysis of the dayside atmosphere of WASP-189 b using high-resolution spectroscopic observations with CRIRES⁺. The cross-correlation analysis confirmed the presence of CO and Fe with signals red-shifted by ~6 km s⁻¹, while a range of other species were not detectable. By performing atmospheric retrieval analyses, we quantitatively constrained the atmospheric and orbital parameters, and we incorporated a detailed modelling of the line profile due to various dynamical effects in the retrieval framework.

Assuming equilibrium chemistry, we found a C/O ratio consistent with the solar value, and a stellar to super-stellar metallicity, consistent with models of late-stage ice accretion or core erosion. The effects of winds on the observed spectra were modelled based on the radial velocity distributions of different wind patterns on a 2D disk representing the planetary dayside. The retrieved line profile fits best with an atmosphere primarily influenced by a fast day-to-night wind, inducing a notable red-shift in the observed spectra. In contrast, the retrieved equatorial jet velocity was significantly lower and the width of the jet remained unconstrained. While some theoretical works predict slower wind speeds in ultra-hot atmospheres, other general circulation models of WASP-189 b result in significantly faster jet speeds that are not consistent with the retrieved line profile. Not the entire offset from the expected planet rest frame can be explained by day-to-night winds, because a further increase of the wind speed not only shifts the signal position but also causes a further broadening, which is not consistent with the observed data. Because of this effect, it is generally difficult to explain large velocity offsets in any dayside spectra purely with day-to-night winds. Possible reasons for the remaining offset in our observation include an inaccurate system velocity and the effect of a sub-stellar hot spot.

Further observations covering a wider range of orbital phases before and after the eclipse will be beneficial to further constrain the dynamical properties. The combination with transmission observations will likely prove useful in explaining the remaining velocity offset and providing a more comprehensive understanding of the atmospheric circulation and thermal structure of WASP-189 b.

Acknowledgements

CRIRES⁺ is an ESO upgrade project carried out by Thüringer Landessternwarte Tautenburg, Georg-August Universität Göttingen, and Uppsala University. The project is funded by the Federal Ministry of Education and Research (Germany) through Grants 05A11MG3, 05A14MG4, 05A17MG2 and the Knut and Alice Wallenberg Foundation. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 109.23HN.002. F.L. acknowledges the support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project number 314665159. D.C. is supported by the LMU-Munich Fraunhofer-Schwarzschild Fellowship and by the Deutsche Forschungs-gemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2094 – 390783311. F.Y. acknowledges the support by the National Natural Science Foundation of China (grant No. 42375118). L.B.-Ch., A.D.R., and N.P. acknowledge support by the Knut and Alice Wallenberg Foundation (grant 2018.0192). S.C. acknowledges the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (CZ 222/5-1). M.R. acknowledges the support by the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (DFG PR 36 24602/41). D.S. acknowledges financial support from the project PID2021-126365NB-C21(MCI/AEI/FEDER, UE) and from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/ 10.13039/501100011033.

Appendix A Additional K_p - v_sys maps

Fig. A.1 Model spectra and Kp-voffset maps of the six undetected species CO2, H2S, HCN, NH3, OH, and FeH.

Appendix B Evaluating the effect of the model filter

Fig. B.1 Comparison of spectral line shapes of synthetic CO models for three scenarios: without any atmospheric dynamics (blue), with a day-to-night wind of 5 km s−1 (red), and with a day-to-night wind of 5 km s−1 and a jet of 3 km s−1 (green). The unprocessed models (dash-dotted lines) were injected into the observed data, to which subsequently SYSREM was applied for 4 iterations. The ratio of the resulting residuals to the residuals of non-injected data show the distortions of the planetary signal introduced by SYSREM (dashed lines). These distortions can be reproduced by applying a filter directly to the original models, resulting in the filtered models (solid lines) which closely resemble the injected models in all three scenarios.

The application of SYSREM can alter the planetary spectrum in the observed data. To account for this effect in our model spectra, we apply a filter based on the approach of Gibson et al. (2022), as outlined in Sect. 4.1. This enables a more accurate comparison between the observed and modelled spectra during the retrieval process. Since our analysis revolves around the line shapes of different planetary features, it is crucial that the filtering process preserves these shapes without introducing major distortions. In the following, we evaluate the filter’s performance by comparing the filtered models to the actual effect of SYSREM on an injected model spectrum.

To this end, we injected a model spectrum into the observed data, applied SYSREM, and divided the result by the non-injected residuals. The ratio of injected to non-injected residuals represents the planetary signal post-SYSREM, which showed clear distortions compared to the original model. Figure B.1 illustrates the comparison between a single CO emission line from the original model, the post-SYSREM model, and the result after applying the filter. Shown are three different scenarios: without any atmospheric dynamics, with a day-to-night wind of 5 km s⁻¹, and with a day-to-night wind as well as a jet of 3 km s⁻¹. Although this example used four iterations of SYSREM, a similar behaviour is observed across other iterations.

SYSREM primarily reduces the amplitude of spectral lines and introduces troughs on either side of the lines. The filtered model accurately replicates these effects, with only minor deviations from the injected model. Importantly, the differences in line shape between the models – such as the broader line width and extended redshifted wing in the wind-included models – are preserved. Although a wind speed of 5 km s⁻¹ was used in this demonstration, the line profile differences are retained also at lower wind speeds. Thus, we conclude that applying this filter to the models before comparing them with SYSREM-processed data enables reliable determination of wind speeds based on the spectral line shapes.

Appendix C Full corner plot and parameters of the wind pattern retrieval

Fig. C.1 Full corner plot of the posterior distributions from the retrieval with a dedicated wind line profile. The dashed lines indicate the mean value and the 16th and 84th percentile. The posteriors of the jet width σjet and limb darkening coefficient ϵ fill the entire parameter space and the percentile values of these two parameters are thus not meaningful.

References

All Tables

All Figures

	Fig. 1 Data reduction steps of a representative wavelength range. Top panel: Unprocessed spectra as they are produced by the CRIRES+ pipeline. Middle panel: After normalisation to bring all spectra onto a common blaze function, and masking of deep telluric lines. Bottom panel: After the removal of stellar and telluric lines with SYSREM.
In the text

Fig. 2 S/N detection strength as a function of SYSREM iterations for the signal of CO, at nodding positions A (blue colour) and B (orange colour). The top panel shows the signal strength for the real data, and the bottom panel shows the strength of an injected signal using the differential ΔCCF according to the method described in Sect. 3.3. The diamond shapes indicate the SYSREM iteration with the strongest detection of the injected signal.

In the text

	Fig. 3 Model spectra (top), cross-correlation functions (CCFs, middle), and Kp-voffset maps (bottom) of CO, Fe and H2O. The dashed lines in the middle panels indicate the expected course of the planetary trail (which is located between the two lines). The signals of CO and Fe are visible in the CCFs and Kp-voffset maps, while H2O was not detectable.
In the text

Fig. 4 Posterior distributions from the retrieval of the atmospheric and orbital parameters. The dashed lines indicate the mean value and the 16th and 84th percentile. The included parameters consists of the three parameters defining the T–p profile (Tirr, κIR, γ), the orbital semi-amplitude of the planet (Kp), the offset from the expected planet rest frame (voffset), a scaling factor to the uncertainties (β), the equatorial rotation velocity (veq) and two parameters determining the equilibrium chemistry (C/O ratio and metallicity).

In the text

Fig. 5 Retrieved temperature and abundance profiles. Left panel: retrieved T–p profile calculated from 10 000 random samples. The shaded region indicates the 1σ uncertainty interval. Middle and right panel: retrieved abundance profiles of the individual species as determined from the C/O and [M/H] values and the T–p profiles of 1000 random samples. The two species in the middle panel have been detected in the cross-correlation analysis, while the species in the right panel were undetectable.

In the text

Fig. 6 Examples of different radial velocity distributions on the planetary disk and resulting velocity profiles, and the result from the retrieval. The example profiles assume a limb darkening coefficient of ϵ = 1. (a) Tidally locked rotation (according to WASP-189 b’s orbital period), (b) Day-to-night wind of vday–night = 5 km s−1. (c) Zonal jet with a speed of vjet = 5 km s−1 and a width of σjet = 0.3 Rp. (d) Velocity distributions according to the best-fitting parameters from the retrieval, and the average retrieved line profile calculated from 1000 random samples. The shaded region indicates the 1σ uncertainty interval.

In the text

	Fig. 7 Posterior distributions from the retrieval of the atmospheric and orbital parameters. The dashed lines indicate the mean value and the 16th and 84th percentile. The posteriors of the jet width σjet and limb darkening coefficient ϵ fill the entire parameter space and the percentile values of these two parameters are thus not meaningful.
In the text

	Fig. A.1 Model spectra and Kp-voffset maps of the six undetected species CO2, H2S, HCN, NH3, OH, and FeH.
In the text

Fig. B.1 Comparison of spectral line shapes of synthetic CO models for three scenarios: without any atmospheric dynamics (blue), with a day-to-night wind of 5 km s−1 (red), and with a day-to-night wind of 5 km s−1 and a jet of 3 km s−1 (green). The unprocessed models (dash-dotted lines) were injected into the observed data, to which subsequently SYSREM was applied for 4 iterations. The ratio of the resulting residuals to the residuals of non-injected data show the distortions of the planetary signal introduced by SYSREM (dashed lines). These distortions can be reproduced by applying a filter directly to the original models, resulting in the filtered models (solid lines) which closely resemble the injected models in all three scenarios.

In the text

	Fig. C.1 Full corner plot of the posterior distributions from the retrieval with a dedicated wind line profile. The dashed lines indicate the mean value and the 16th and 84th percentile. The posteriors of the jet width σjet and limb darkening coefficient ϵ fill the entire parameter space and the percentile values of these two parameters are thus not meaningful.
In the text