© The Authors 2026

1 Introduction

Hot Jupiters are gas giant exoplanets on close-in orbits around their host stars, with orbital periods on the order of a few days. The intense stellar irradiation they receive leads to elevated equilibrium temperatures (T_eq) of about 1500 K, resulting in extended atmospheric scale heights and strong thermal emission fluxes. This makes these exoplanets particularly favorable targets for atmospheric characterization via transmission and emission spectroscopy. Hot Jupiters typically have circular orbits and are expected to exhibit synchronous rotation, both of which are caused by tidal circularization during the early evolution of these planetary systems (Hut 1981).

Despite the intense stellar irradiation they receive, the atmospheres of hot Jupiters are sufficiently cool to allow the presence of significant abundances of molecular species, as most molecules thermally dissociate only at T_eq beyond approximately 2000 K (Parmentier et al. 2018). Both space- and ground-based observations have identified CO and H₂O as dominant contributors to the spectra of hot Jupiters at near-infrared wavelengths (e.g., Snellen et al. 2010; Wakeford et al. 2018). These species contain a substantial fraction of the atmospheric carbon and oxygen and therefore serve as key tracers of a planet’s carbon-to-oxygen (C/O) ratio, which is commonly used to probe exoplanet formation pathways (e.g., Öberg et al. 2011; Mordasini et al. 2016; Williams et al. 2025). Beyond CO and H₂O, numerous other molecular species have been discovered in the atmospheres of hot Jupiters, such as CO₂, CH₄, and SO₂ (e.g., Giacobbe et al. 2021; Carter et al. 2024; Xue et al. 2024).

The opacities of molecular species are generally lower at visible than at infrared wavelengths. As a result, stellar radiation can penetrate efficiently into hot Jupiter atmospheres, while the thermal reemission of the absorbed energy at infrared wavelengths is comparatively less efficient. This effect causes an increase in temperature with atmospheric pressure in hot Jupiter atmospheres, yielding a non-inverted temperature-pressure (T-p) profile (e.g., Fortney et al. 2008; Guillot 2010).

The thermal conditions in hot Jupiter atmospheres allow refractory elements to condense out of the gas phase, leading to the formation of clouds. The pressure level at which these clouds are located depends on a planet’s T-p profile and atmospheric composition (Helling 2019). In addition, photochemical reactions in the upper atmospheres of hot Jupiters result in complex molecular species that form solid particles referred to as hazes (Gao et al. 2021). Clouds and hazes can alter the spectra of hot Jupiter atmospheres. For instance, these aerosols can obscure the spectral features of the underlying atmospheric layers. Moreover, clouds and hazes affect the atmospheric T-p profile through their strong influence on a planet’s albedo as well as by inducing cooling and warming effects above and below the cloud deck, respectively (Heng et al. 2012; Gao et al. 2021).

General circulation models predict the presence of strong dynamical processes in the atmospheres of hot Jupiters. These large-scale atmospheric gas flows are classified into two regimes: dayside to nightside winds and super-rotating equatorial jet streams, both of which can reach velocities of several km s⁻¹ (e.g., Miller-Ricci Kempton & Rauscher 2012; Showman et al. 2013; Tan & Komacek 2019). Detections of Doppler offsets relative to the planetary rest frame and broadening in the spectral line profiles have confirmed the presence of such global-scale wind patterns in observational studies. For instance, Nortmann et al. (2025) revealed the presence of a fast super-rotating jet stream in the atmosphere of the hot Jupiter WASP-127b, while Seidel et al. (2025) identified dayside to nightside winds and jet stream circulation at different pressure levels in the atmosphere of WASP-121b.

High-resolution spectroscopy has proven to be a powerful tool for probing the chemical composition, thermal structure, and atmospheric dynamics of hot Jupiters (e.g., Snellen et al. 2010; Birkby et al. 2013; Line et al. 2021; Rafi et al. 2024; Kanumalla et al. 2024; Finnerty et al. 2025). Beyond providing detailed insights into atmospheric conditions, this technique enables the derivation of key orbital and planetary parameters, such as the planetary mass and orbital inclination. For transiting hot Jupiters, determination of the planetary mass is relatively straightforward by using the radial velocity (RV) method, which traces the periodic Doppler shift of the spectral lines of an exoplanet’s host star. In this specific case, the orbital inclination, which is degenerate with planetary mass in the RV signal, is known from the transit geometry. However, RV measurements of the host star are unable to break the degeneracy between planetary mass and orbital inclination for non-transiting hot Jupiters. This limitation can be overcome by measuring the thermal emission signal from hot Jupiter atmospheres, which allows the planetary orbit to be traced directly. So far, this approach has been successfully applied to the study of the planetary and orbital parameters of a few bright, non-transiting hot Jupiters, including τ Boo b (Brogi et al. 2012), HD 179949 b (Brogi et al. 2014), and 51 Peg b (e.g., Brogi et al. 2013; Birkby et al. 2017).

51 Peg b was the first exoplanet ever detected orbiting a main-sequence star (Mayor & Queloz 1995) and the first hot Jupiter ever discovered. As the planet does not transit in front of the host star from the perspective of the Earth-bound observer, its atmospheric properties cannot be probed via transmission spectroscopy. Atmospheric characterization is therefore only possible via the analysis of reflected stellar light or the planet’s thermal emission spectrum. A number of studies have attempted to detect reflected stellar light from 51 Peg b. Martins et al. (2015) analyzed HARPS-N observations of the planet in the visible wavelength range and reported tentative evidence of a reflection signal. However, this signal could only be explained by assuming an elevated geometric albedo and an inflated planetary radius. Using the same dataset, Borra & Deschatelets (2018) also found evidence of reflected stellar light in the spectrum of 51 Peg b. In contrast, a reanalysis of these observations by Di Marcantonio et al. (2019) yielded a non-detection of the reflected stellar signal. A subsequent study by Scandariato et al. (2021), using the same data along with additional HARPS-N measurements, derived a geometric albedo below 0.15 and concluded that no reflected stellar light could be detected. More recently, Spring et al. (2022) reported the absence of significant reflected light in HARPS-N data. While robust detections of reflected stellar light remain elusive, the thermal emission spectrum of 51 Peg b has been successfully detected at high spectral resolution in the infrared, providing a direct window into its atmospheric composition and thermal structure. Brogi et al. (2013) first identified CO and H₂O spectral lines in the planet’s emission spectrum using K band observations with the old CRIRES. Later, Birkby et al. (2017) reported detection of H₂O in the L band using the same instrument. Additionally, Chiavassa & Brogi (2019) confirmed the presence of CO and H₂O features by reanalyzing the data from Brogi et al. (2013) with improved methods for removing stellar contamination.

This present study provides the first spectroscopic characterization of the benchmark hot Jupiter 51 Peg b using Bayesian inference methods and data from the upgraded CRIRES⁺ instrument, offering a robust quantitative analysis of its atmospheric thermal structure, chemical composition, and dynamical processes. This target was identified as one of the most suitable systems for high-resolution spectroscopic studies of exoplanet atmospheres by Shulyak et al. (2019). All parameters of the 51 Peg system used in this work are summarized in the Table 1. The paper is organized as follows: In Sects. 2 and 3, we describe the observations, data reduction, and preprocessing procedures. In Sect. 4, we update the ephemeris of the 51 Peg system. The methodology used to identify spectral features in the planet’s emission spectrum and the Bayesian framework employed to constrain its atmospheric properties are detailed in Sect. 5. The results are presented and discussed in Sect. 6, and our main conclusions are summarized in Sect. 7.

2 Observations and data reduction

We observed the thermal emission spectrum of 51 Peg b on 6 October 2024 between 01:17 UT and 04:10 UT with the CRIRES⁺ instrument at the Unit Telescope 3 of the Very Large Telescope (Dorn et al. 2023). The observations cover orbital phases 0.52–0.55, which are close to superior conjunction, when the planetary dayside is most favorably aligned with the observer’s line of sight. A total of 130 consecutive spectra were collected. The adaptive optics (AO) system was disabled during the acquisition of the first 28 spectra due to passing clouds; the remaining 102 spectra were obtained under improved meteorological conditions that allowed for AO usage. The AO system was disabled under cloud coverage because its wavefront sensor has a limited dynamic range that cannot be adjusted during observations. Rapid flux variations from passing clouds can cause the flux to exceed this range, which could damage the wavefront sensor. We verified that including data without AO does not affect our results. Both the cross-correlation and retrieval analyses in Sect. 6 yield consistent outcomes with or without these spectra. Each spectrum was taken with DIT = 15 s and NDIT = 4, resulting in a total integration time of 60 s per frame. Splitting each frame into four sub-integrations was required to prevent saturation of the detector pixels, given the bright host star of our target.

We used an AABBBBAA nodding pattern, which consists of observing the target at two different slit positions (A and B), to facilitate the removal of the sky background during the subsequent data reduction steps. The slit width was set to 0.2″, which corresponds to a nominal resolving power of R ~ 100 000 under the condition of a homogeneous slit illumination. However, during our observations, the seeing conditions were very good, such that the AO system – when active – could achieve a point spread function (PSF) smaller than the slit. The reduced PSF width increased the spectral resolution from R ~ 100 000 (28 spectra with no AO) to R ~ 117 000 (102 spectra with AO). We determined the spectral resolution using molecfit (Smette et al. 2015) via the model kernel full width at half maximum (FWHM) stored in the BEST_FIT_PARAMETERS output file. The FWHM in detector pixels was first converted to wavelength and then to spectral resolution. For the AO assisted spectra affected by super-resolution, the FWHM corresponds to approximately 2.6 pixels. This exceeds the two-pixel minimum required for Nyquist sampling, ensuring adequate sampling of the instrumental profile. The data were acquired in the K2166 wavelength setting, which offers a wavelength coverage between 1921 nm and 2472 nm. This setting consists of seven echelle orders, each divided by two narrow gaps between detectors, yielding a total of 21 wavelength segments. The airmass decreased from 1.5 to a minimum of 1.4 before increasing to 1.6; the mean signal-to-noise (S/N) ratio per detector pixel was 330 and 485 for spectra obtained without and with AO, respectively.

We used the data reduction pipeline cr2res to extract the one-dimensional spectra from the raw frames¹. This procedure includes correction for dark frames and flat fields, removal of bad pixels, and calculation of a wavelength solution. The integration time of the individual frames was shorter than the timescale over which 51 Peg b’s RV changes significantly with respect to the observer (60 s compared to ~400 s to smear the planetary signal over 1 km s⁻¹, which is on the order of one CRIRES⁺ detector pixel). Therefore, we combined adjacent AA and BB frame pairs into one A frame and one B frame, respectively. This limited the subsequent data reduction to the analysis of a 65 frame ABBA nodding sequence, with an integration time of 120 s per frame.

The wavelength solution for the A and B spectra, as calculated by the cr2res pipeline, is based on data collected during daytime calibrations. However, this wavelength solution can differ from that of the nighttime science data due to instrumental drifts or offsets between subsequent frames caused by changes in ambient or observing conditions (e.g., variations in temperature, pressure, or slit illumination). Therefore, we used molecfit to refine the wavelength solution of our spectra. Following Cont et al. (2024), we first derived wavelength solutions for the time averaged A and B spectra. Subsequently, we determined the offsets between the individual spectra and the corresponding A or B wavelength solution. To this end, we computed the cross-correlation between each individual spectrum and the respective A or B time averaged spectrum. Finally, we shifted each A and B spectrum to the rest frame of its time averaged spectrum. Figure A.1 shows how these procedures align the individual spectra to the refined wavelength solution. We treated the A and B spectra as two independent datasets throughout the data processing procedures in Sect. 3, and combined their information in Sect. 5.

3 Preprocessing the spectra

We normalized the extracted spectra of our time series to a common continuum level. To this end, we first computed the median of all spectra in the time series and fitted the resulting master spectrum with a second-order polynomial. Each individual spectrum was then divided by this master and fitted with a first-order polynomial. The final continuum correction was obtained by multiplying these two polynomial fits and dividing each spectrum by the result. The continuum-normalized spectra were then arranged chronologically in a two-dimensional array to form a spectral matrix. We removed outliers by 5σ clipping and masked wavelength regions with deep telluric and stellar lines. In total, 19% of the data points was masked. A more detailed description of these procedures can be found in Cont et al. (2024, 2025).

Telluric and stellar contamination were removed with SYSREM from the normalized spectral matrix (Tamuz et al. 2005). This iterative principle component analysis algorithm has facilitated the study of exoplanet atmospheres at high spectral resolution in recent years (e.g., Birkby et al. 2013; Birkby 2018). We followed the procedure proposed by Gibson et al. (2022), which involves dividing the spectral matrix by the median spectrum before applying the SYSREM algorithm. This approach preserves the relative strengths of the planetary spectral lines and allows one to account for distortions introduced by the algorithm when applying a Bayesian retrieval framework in Sect. 5.2. We ran SYSREM for up to ten consecutive iterations, producing a residual spectral matrix for each iteration and wavelength segment. The residual spectral matrix is expected to contain only the planetary spectral signal, which is dominated by noise.

To determine the optimal number of SYSREM iterations, we assessed the relative change in standard deviation of the residual spectral matrix from one iteration to the next. This approach has been successfully used in previous studies (e.g., Parker et al. 2025; Peláez-Torres et al. 2026). The relative change in standard deviation is defined as $$\mathrm{\Delta }\sigma =\frac{{}^{(i-1)}\sigma -{}^{(i)}\sigma }{{}^{(i-1)}\sigma },$$(1)

where ^(i–1)σ and ⁽ⁱ⁾σ denote the standard deviation of the residual spectral matrix before and after the i-th SYSREM iteration across all wavelength segments. We adopted a conservative approach by determining a common optimal iteration number across all wavelength segments, rather than determining the optimal number for each segment individually. Segment-by-segment optimization has been shown to pose a risk of false-positive detections (e.g., Sánchez-López & Millán 2025). The Δσ value of our data plateaus after six iterations (Fig. A.2). This indicates that SYSREM no longer identifies significant correlated trends beyond this number of iterations.

4 Updated ephemeris from RV analysis

Accurate planetary ephemerides are essential for precise orbital phase assignment of spectroscopic data, which in turn allows for more reliable constraints on atmospheric properties. To update the ephemeris of 51 Peg b, we conducted an RV analysis including publicly available data of 51 Peg from Levay/APF (Vogt et al. 2014) and VLT/ESPRESSO (Pepe et al. 2021), as well as the Lick/Hamilton, Keck/HIRES, OHP/ELODIE, and ESO 3.6-m Telescope/HARPS RV datasets from Birkby et al. (2017). We accessed the APF measurements made available in Rosenthal et al. (2021), which consist of a total of 283 epochs with a typical RV uncertainty of 2.7 m s⁻¹ spanning five years from 16 October 2013 to 14 October 2018. The ESPRESSO data was accessed through the ESO archive and reduced with the instrument data reduction pipeline version 3.3.1². We extracted the RV information via cross-correlating the spectra with a G2 binary mask using default pipeline settings. The 90 ESPRESSO RV data points were obtained on three different nights covering a baseline of 12 days in August 2019. The very high RV precision of 0.15 m s⁻¹ allowed this data to add significant information content to the final fit. Together, the datasets we used span more than 25 years. For the final modeling of the data, we ensured all timestamps were in Barycentric Dynamical Time (BJD_TDB).

To model the Keplerian orbital parameters, we used PYANETI (Barragán et al. 2019) employing a sinusoid fit and Markov chain Monte Carlo (MCMC) sampler to assess the uncertainties. The datasets of the different instruments were treated independently to allow for offsets in systemic velocity and jitter terms. The Hamilton data were split up into three versions of the instrument to account for instrument interventions, as described in Birkby et al. (2017). For the stellar parameters, we used the most recent literature results from Rosenthal et al. (2021). In a first model run, we treated the orbital eccentricity as a free parameter. After obtaining a two-sigma upper eccentricity limit of 0.0063, which is consistent with a circular orbit, we repeated the RV analysis with the parameter value fixed to zero. The results of this zero-eccentricity fit were then used for the subsequent analyses of 51 Peg b’s atmospheric and planetary properties, based on CRIRES⁺ observations (Sect. 6). Figure 1 shows all data points together with the best-fit Keplerian model, and Table 2 lists the results of the RV analysis.

Fig. 1 Radial velocity curve of 51 Peg with the updated orbital solution. We show the phase-folded RV data points and the best-fit Keplerian model.

5 Atmospheric analysis methods

5.1 Cross-correlation

We used the radiative transfer code petitRADTRANS to model emission spectra of 51 Peg b’s atmosphere (Mollière et al. 2019). In a preliminary cross-correlation analysis using the T-p profile and chemical abundances from Brogi et al. (2013), we achieved a significant detection of the planetary spectral lines in our data. We then applied the retrieval framework described in Sect. 5.2 to derive the atmospheric thermal and chemical properties (see Sect 6.2). The best-fit T-p profile and volume mixing ratios (VMRs) of the considered chemical species were subsequently used to generate the model spectra corresponding to the cross-correlation analysis results presented in Sect. 6.1. The petitRADTRANS code was provided with the spectral line lists of the chemical species under investigation. The CO line data were taken from Li et al. (2015), and the H₂O data were taken from the POKAZATEL database (Polyansky et al. 2018). We generated spectra for CO and H₂O individually, as well as for the combination of the two species.

Each spectral model underwent two normalization steps. First, the models were divided by the blackbody spectrum of the host star. Then, the resulting planet-to-star flux ratio was normalized to the planetary continuum. Convolving the spectra with the instrumental profile yielded the final model spectra used for cross-correlation. The final model spectra used in this study are shown in the left panels of Fig. 2.

We computed the weighted cross-correlation function (CCF) between each model spectrum and the residual spectra. For the i-th residual spectrum of the time series, the CCF is defined as $${\mathrm{CCF}}_i(v)=\sum _j\frac{R_{ij}{m}_j(v)}{{\sigma }_{ij}^2},$$(2)

with R_ij being the residual spectral matrix, σ_ij the corresponding uncertainties, and m_j the spectral model shifted by a Doppler velocity in the range −800 km s⁻¹ < v < +800 km s⁻¹. For each wavelength segment, the CCFs of the spectral time series were arranged in a two-dimensional array. These arrays were then co-added to create a CCF map. Finally, the CCF maps from the two A and B nodding positions were combined into one final CCF map for each model spectrum. The middle panels of Fig. 2 illustrate the final CCF maps obtained using these procedures.

We aligned the CCF map to the rest frame of 51 Peg b over a range of different orbital semi-amplitude velocity (K_p) values. Under the assumption of a circular orbit, the planetary Doppler velocity is described by $$v_{\mathrm{p}}(t)=v_{\text{sys}}+v_{\text{bary}}(t)+K_{\mathrm{p}}\sin 2\pi \varphi (t),$$(3)

where v_sys is the systemic velocity, v_bary (t) is the barycentric velocity correction, and ϕ (t) is the orbital phase. Each shifted CCF map was collapsed into a one-dimensional CCF by averaging along the time axis. The one-dimensional CCF of each alignment was further stacked into a two-dimensional array and normalized by its standard deviation, yielding the so-called S/N map. Following Cont et al. (2024, 2025), the standard deviation was obtained by fitting the distribution of all CCF values with a Gaussian function. The S/N maps obtained in this work are shown in the right panels of Fig. 2.

5.2 Retrieval

The procedures for retrieving the atmospheric properties of 51 Peg b follow a framework similar to that described in Sect. 5.1.1 of Cont et al. (2024). This framework relies on petitRADTRANS to forward model the planetary emission spectrum (Mollière et al. 2019).

To model the atmospheric thermal conditions, we used the T-p profile parameterization introduced by Pelletier et al. (2021). Our choice is motivated by the fact that this parameterization does not make any a priori assumptions about the overall shape of the atmospheric thermal structure. Consequently, it can adapt flexibly to a wide range of possible T-p profiles. We used six free temperature points (T₁–T₆), uniformly distributed in the logarithmic pressure space between 10⁻⁸ bar and 10² bar. Following the methodology of Pelletier et al. (2021), we adopted a smoothing prior on the T-p profile with a standard deviation of σ_s = 150 K dex⁻² to prevent nonphysical temperature oscillations in the retrieved temperature curve. Since the atmosphere of 51 Peg b is predicted to host magnesium-silicate and iron clouds (Goukenleuque et al. 2000), we included a parameterization of an opaque cloud deck in our retrieval framework, which mutes the spectral contribution from atmospheric layers at pressures higher than p_cloud.

The VMRs of the chemical species in 51 Peg b’s atmosphere were modeled with the FastChem chemical equilibrium code as a function of the metallicity ([M/H]) and the C/O ratio (Stock et al. 2022)³. We included the spectral contribution of CO and H₂O into the retrieval. These are the chemical species with the most dominant spectral lines in the K band. We also accounted for Rayleigh scattering from H₂ and He, as well as H₂–H₂ and H₂–He collision-induced absorption effects.

The petitRADTRANS code computes the flux per unit surface area and wavelength emerging from a planet’s atmosphere. Therefore, the overall flux of 51 Peg b is obtained via integration over the entire visible hemisphere, which requires knowledge of the planetary radius (R_p). However, R_p of 51 Peg b is poorly constrained but is expected to lie within the range 1.0–1.3 R_Jup (Birkby et al. 2017). Consequently, we treated R_p as a free parameter and imposed an informed Gaussian prior centered on 1.15 R_Jup with a standard deviation of 0.15 R_Jup.

The surface gravity value (log g) of 51 Peg b is also unknown. This value is a function of R_p and a required input parameter of petitRADTRANS. For this reason, we recalculated log g via $g=G{M}_{\mathrm{p}}/R_{\mathrm{p}}^2$ at each retrieval step with the updated R_p value. The planetary mass (M_p) is calculated using M_p/M_* = K_*/K_p and G denotes the gravitational constant. The K_p value is updated at each retrieval step, and M_* and K_* denote the stellar mass and orbital semi-amplitude. Thus, in addition to the free parameters of our retrieval, we can estimate 51 Peg b’s log g and M_p values.

We also accounted for spectral line broadening due to the planetary rotation. To this end, we convolved the spectral forward model with a broadening profile. The broadening profile is parametrized by the projected equatorial rotation velocity (v_rot sin i). The analytical description of the broadening profile is given in Eq. (3) of Díaz et al. (2011) and has been used successfully in previous studies (e.g., Yan et al. 2023; Cont et al. 2024, 2025; Lesjak et al. 2025).

As described in Sect. 5.1, the forward modeled spectrum was converted to the planet-to-star flux ratio and convolved with the CRIRES⁺ instrument profile. Subsequently, for each individual spectrum of our time series, we used K_p and v_sys to Doppler shift the spectral forward model to the planetary rest frame. The shifted model spectrum was then interpolated to the wavelength solution of our observational data. These procedures resulted in a two-dimensional forward model matrix, with the same shape as the residual spectral matrix. We applied the filtering method of Gibson et al. (2022) to account for the presence of potential distortion effects introduced into our data by SYSREM.

We calculated the Gaussian log likelihood function $$\ln L=-\frac{1}{2}\sum _{i,j}[\frac{{(R_{ij}-M_{ij})}^2}{{\left(\beta {\sigma }_{ij}\right)}^2}+\ln 2\pi {\left(\beta {\sigma }_{ij}\right)}^2]$$(4)

to compare the model spectrum to the data. In this expression, R_ij and M_ij denote the individual elements of the residual spectral matrix and the filtered forward model matrix, respectively. The uncertainties of the residual spectra are represented by σ_ij and β acts as a noise scaling factor. We used the residual spectral matrix obtained after six consecutive SYSREM iterations, at which Δσ reached a plateau (see Sect. 3). For each wavelength segment and nodding position, the log likelihood function was calculated independently, and the resulting values were added to obtain the combined log likelihood function of all data. We estimated the model parameters by evaluating the combined log likelihood function using MCMC sampling with the emcee package (Foreman-Mackey et al. 2013).

In summary, our high-resolution retrieval framework includes the following free parameters: the T-p profile parameters T₁–T₆; the pressure level of an opaque cloud deck p_cloud; the atmospheric chemical properties represented by [M/H] and the C/O ratio; the planetary radius R_p; the projected equatorial rotation velocity v_rot sin i; the velocity parameters K_p and v_sys; and the noise scaling parameter β. From these free parameters, additional planetary properties can be derived: the surface gravity log g and the planetary mass M_p. For each free parameter, 32 walkers with 15 000 steps were used in the sampling.

Fig. 2 Model spectra, CCF maps, and S/N maps of CO, H2O, and the combination of both species. Left panels: model spectra normalized to the spectral continuum level. Wavelengths covered by the CRIRES+ K2166 setting used in this work correspond to the gray shaded area. Middle panels: cross-correlation function maps. The signal of H2O and the signal obtained from a combination of CO and H2O spectral lines can be identified as diagonal trails in the mappings, indicated by the solid white lines. A residual CO signal, which is consistent with the stellar rest frame and could not be entirely removed with SYSREM, is masked. Right panels: signal-to-noise maps. These maps reveal no clear CO signal but a detection of spectral H2O lines. The position of the planetary signal (Kp value reported by Brogi et al. (2013) and Birkby et al. (2017); vsys determined via the RV analysis in Sect. 4) is indicated by the dashed white lines. We show the CCF maps and S/N maps obtained after six consecutive SYSREM iterations. A non-masked version of the CO and the species-combined CCF and S/N maps is given in Fig. A.3.

6 Results and discussion

6.1 Detection of the planetary spectral lines

We detect line features in the emission spectrum of 51 Peg b’s atmosphere. Figure 2 shows the S/N maps obtained by cross-correlating our CRIRES⁺ spectra with the spectral models of CO, H₂O, and a combination of both species. The white dashed lines indicate the K_p reported by Brogi et al. (2013) and Birkby et al. (2017), and the v_sys from our RV analysis in Sect. 4. The maps obtained under the inclusion of H₂O spectral lines show significant emission signals that coincide in velocity space but differ from the literature values by a few km s⁻¹. We provide a detailed discussion of this offset in Sect. 6.2, together with the results of our Bayesian retrieval analysis.

We identify a clear S/N peak of atmospheric H₂O with a detection strength of S/N = 6.2. In contrast, we do not detect CO at a statistically significant level, with a S/N of only 1.9 at the expected planetary position. This is well below the S/N threshold of 4 commonly adopted in the literature. Nevertheless, we note that a diagonal feature appears at the exact K_p–v_sys position of the robust H₂O detection. While this weak CO feature alone cannot be considered statistically significant given the noise structure in the CO S/N map, its spatial coincidence with the H₂O signal is noteworthy, particularly in light of the previous CO detection reported by Brogi et al. (2013) using the original CRIRES instrument. A direct comparison with their result is not possible, as Brogi et al. (2013) do not provide a CO S/N map but only the combined CO + H₂O S/N map. Given that CRIRES⁺ covers a broader wavelength range containing CO lines than the preupgraded CRIRES, one might expect a stronger CO signal in our data, which makes the absence of a significant detection in our analysis somewhat puzzling. However, the cumulative exposure time of our observations is considerably shorter (~3 h) compared to that of Brogi et al. (2013) (~16 h), which likely explains why the original CRIRES data achieved a significant CO detection despite its more limited spectral coverage. We note that a residual stellar CO signal is present in the data and could, in principle, bias the interpretation of the planetary signal if they overlap in the CCF map. As shown in Fig. 2, this is not the case in our analysis, as the planetary signal is clearly separated from the stellar CO residuals in velocity space. Consequently, the planetary signal and our derived results remain unaffected by these residual features. Furthermore, the inclusion of CO into the retrieval in Sect. 6.2 allowed us to place constraints on the C/O ratio, which would not be possible if the CO signal were purely spurious. Therefore, we consider it plausible that the diagonal CO feature with a S/N of 1.9 indicates a weak contribution of CO to the planetary spectrum. However, we do not claim a CO detection, and additional observations are required.

Cross-correlation with a model spectrum that includes both CO and H₂O spectral lines yields a detection peak with S/N = 7.2. Comparison of this detection strength with that of the individual species reveals that the spectral signature of H₂O dominates the species-combined signal. This is also visible in the model spectra, where the combined spectrum of all species closely matches that of H₂O. The relatively small contribution of CO to the species-combined spectrum can be explained by the different atmospheric layers of 51 Peg b in which CO and H₂O are spectroscopically active. The CO lines mostly form at higher pressures and deeper in the atmosphere, while the H₂O lines form at lower pressures and higher altitudes, as shown in Fig. 3. Consequently, the CO spectral signature is largely obscured by the overlying H₂O spectral lines.

With an effective temperature of 5761 ± 95 K (Rosenthal et al. 2021), 51 Peg is expected to exhibit an elevated number of strong CO lines in its stellar spectrum. We observed that residual stellar CO features persisted in our data throughout the entire range of SYSREM iterations. Therefore, we masked the CO and species-combined CCF maps within a v_sys interval of ±10 km s⁻¹ around the stellar RV. This procedure propagated into the respective S/N maps as a small masked region. However, the affected region does not overlap with the planetary signal of 51 Peg b, and thus it does not affect our conclusion that a robust signal is found in the CO + H₂O combined analysis and that a minor CO feature at the same position as the detected H₂O signal can be identified. An unmasked version of the CCF and S/N maps is provided in Fig. A.3.

The detection of absorption lines in the dayside emission spectrum of 51 Peg b provides strong evidence of the presence of a non-inverted T-p profile. This finding is in line with previous investigations of other hot Jupiter atmospheres (e.g., Bazinet et al. 2024; Finnerty et al. 2024).

Fig. 3 Contribution functions of the CO, H2O, and species-combined model spectra. The CO signal originates from deeper regions of the atmosphere, whereas the H2O lines form at higher altitudes. Consequently, in the species-combined model, the CO spectral lines are largely muted by the overlying H2O signal. The contribution of CO to the combined signal becomes most evident in the region around the CO bandhead at approximately 2.3 μm.

6.2 Retrieval of atmospheric and planetary properties

Our retrieval is able to constrain the atmospheric thermal and chemical structure, and it allowed us to derive additional planetary properties from the Doppler shift of 51 Peg b’s orbital motion. Figure B.1 shows the posterior distributions, together with the correlations between the model parameters. A summary of the best-fit parameters of the retrieval is provided in Table 3.

6.2.1 Free parameters

Figure 4 shows the retrieved T-p profile, which confirms the presence of a non-inverted thermal structure already detected in Sect. 6.1. This finding is consistent with the thermal properties derived in previous observational studies of 51 Peg b (Brogi et al. 2013; Birkby et al. 2017; Chiavassa & Brogi 2019). Our retrieval yields a T-p profile spanning between approximately 1400 K to 3000 K from the upper to the lower planetary atmosphere. The retrieved thermal conditions of 51 Peg b’s upper atmosphere are consistent with the planetary T_eq, estimated via Eq. (4) of Cowan & Agol (2011) as approximately 1330 and 1700 K in the extreme cases of full and no heat redistribution from the planet’s dayside to nightside. Toward the bottom of 51 Peg b’s atmosphere the retrieved T-p profile is hotter than the calculated T_eq values. This is in line with the expectation that hot Jupiters cannot efficiently reradiate significant amounts of energy received from their host stars back to space from their lower atmospheres (e.g., Fortney et al. 2008; Guillot 2010).

Our retrieval constrains the pressure level below which an opaque cloud deck may exist to p_cloud > 10^−1.6 bar. In Figure 4, we overplot the condensation curves of Fe, Mg₂SiO₄, and MgSiO₃ from Visscher et al. (2010) on the retrieved T-p profile. In the deeper layers of 51 Peg b’s atmosphere temperatures are too high for these species to condense. The condensation curves intersect the retrieved T-p profile at approximately the same pressure as our derived p_cloud limit. At low atmospheric pressures, the condensation curves remain close to the T-p profile, suggesting that condensation could in principle also occur higher up in the atmosphere. However, the retrieval excludes the presence of a cloud deck in this low-pressure region. Consequently, the most likely cloud formation region lies near a pressure of ~10^−1.6 bar, where the condensation curves first intersect the T-p profile and where the retrieval allows for the presence of a cloud deck. The log p_cloud posterior distribution is very similar to those obtained from high-resolution spectroscopy studies of other hot Jupiters with comparable equilibrium temperatures and surface gravities, such as WASP-43b (Lesjak et al. 2023) and HD 189733b (Boucher et al. 2021; Finnerty et al. 2024), noting that we report a two-sigma lower limit rather than median values. In contrast, the cloud deck pressure retrieved for τ Boo b by Pelletier et al. (2021) is higher, consistent with the higher log g of that planet.

The C/O ratio in 51 Peg b’s atmosphere is constrained to ${0.54}_{-0.18}^{+0.23}$, while we derived [M/H] = ${2.63}_{-1.00}^{+0.93}$ dex relative to the solar value. The derived C/O ratio of 51 Peg b is consistent with that of its host star, estimated as approximately 0.5 and 0.7 when using the abundance values of Brewer et al. (2016) and Maldonado et al. (2018), respectively. The planetary [M/H] value is significantly higher that that of the host star, reported as 0.2 dex in both studies. The retrieved value for [M/H] is high even when compared to the Saturn mass planet HD 149026b which is considered the most metal-rich giant planet atmosphere detected to date (Bean et al. 2023). However, 51 Peg b’s mass lies at the border to the hot Saturn mass regime, where planet formation and interior models predict substantially enhanced heavy-element enrichment. This is also reflected in the observed relation between planet mass and atmospheric (and bulk) metallicity (Kreidberg et al. 2014; Welbanks et al. 2019; Sun et al. 2024). It is possible that the high atmospheric metallicity measured for 51 Peg b in this study could be a reflection of a quenched atmosphere dominated by heavy elements from the planet interior. In addition, we note that our retrieved [M/H] uncertainties of approximately ±1 dex are relatively large. Thus, lower metallicity values at the 1.5 dex level, which have previously been found in other gas giant exoplanets, are consistent with our measurement.

Figure 4 shows how the retrieved [M/H] and C/O values translate into the VMRs of CO and H₂O as a function of atmospheric pressure. The VMRs are approximately constant throughout the investigated pressure range, except for H₂O at the lowest pressures, where a slight decrease can be observed. At a representative atmospheric pressure of 10⁻⁴ bar, the retrieved logarithmic VMRs of CO and H₂O are $-{1.10}_{-1.02}^{+0.53}$ and $-{1.45}_{-0.78}^{+0.49}$, respectively, which is higher than previous estimates via cross-correlation (Brogi et al. 2013; Birkby et al. 2017). We note that the relatively constant VMR profiles of CO and H₂O hint – as mentioned in the previous paragraph – toward significant chemical quenching, i.e., atmospheric mixing timescales are shorter than chemical reaction timescales, causing the molecular abundances to remain fixed at the values set in the deeper atmosphere (Moses 2014; Heng & Tsai 2016). However, our retrieval framework assumes thermochemical equilibrium and does not explicitly include quenching. We acknowledge that the assumption of thermochemical equilibrium represents a limitation of our retrieval, as the exclusion of quenching could potentially bias the inferred chemical abundances. This may explain the high [M/H] value retrieved, as a higher metallicity is required under the assumption of equilibrium chemistry to match the elevated molecular abundances resulting from quenching.

To ensure that the high [M/H] value is not the result of including CO emission lines, which contribute negligibly to the planetary detection in Sect. 5.1, we performed an additional retrieval considering only H₂O and excluding CO (the C/O ratio is fixed to the stellar value in this retrieval to exclude any dependence on CO). This retrieval converged on the same parameter values as the retrieval that included CO. Therefore, we conclude that our results are driven by the strong H₂O signal and that the elevated [M/H] value is not a consequence of including CO in our analysis.

For 51 Peg b’s planetary radius, we applied a Gaussian prior based on the previous estimate of 1.01–1.32 R_Jup from Birkby et al. (2013). Our retrieval yields R_p = 1.07 ± 0.14 R_Jup, which is in agreement with this earlier estimate.

We observe a spectral line width that exceeds the broadening expected from 51 Peg b’s planetary rotation under the assumption of tidal locking. The derived projected equatorial rotation velocity is v_rot sin i = ${4.07}_{-0.60}^{+0.53}$ km s⁻¹, whereas the expected tidally locked value lies between 1.2 km s⁻¹ and 1.6 km s⁻¹ for R_p values between 1.0 R_Jup and 1.3 R_Jup. The presence of a velocity excess in the retrieved v_rot sin i value suggests strong atmospheric dynamics, as phenomena such as super-rotating winds or turbulent flow can produce additional line broadening beyond what is expected from planetary solid body rotation. Strong jet stream winds have already been observed in the atmospheres of other gas giant exoplanets, such as WASP-127b and WASP-121b (Nortmann et al. 2025; Seidel et al. 2025). Our interpretation that jet stream winds or turbulent mixing are responsible for the excess line broadening in 51 Peg b’s spectrum rather than a dayside to nightside flow regime is motivated by the fact that the observed signal is not affected by an offset in the retrieved v_sys (e.g., Alonso-Floriano et al. 2019; Sánchez-López et al. 2019; Yan et al. 2022). Indeed, Lesjak et al. (2025) have shown that a dayside to nightside wind flow would result in an asymmetric spectral line shape, producing a net offset from the expected position in the v_sys value.

We emphasize that care must be taken to ensure that the measured spectral line broadening is not biased by instrumental effects, data processing, or overly simplistic modeling. Several effects, such as overestimation of the instrumental spectral resolution, wavelength instabilities during the spectral time series, erosion of line shapes by SYSREM, or neglecting atmospheric pressure broadening, can introduce biases into the inferred rotational broadening if not properly accounted for. In our analysis, we have addressed these effects as follows: the spectral resolution was estimated via molecfit fitting, which we consider a reliable method. Spectral drifts were corrected during data preprocessing steps. The potential deformation of line shapes by SYSREM was accounted for by filtering our forward model spectra according to Gibson et al. (2022). Finally, atmospheric pressure broadening is included in the spectral forward models via the opacities of the radiative transfer calculation. We therefore conclude that the observed line broadening is likely of atmospheric origin, consistent with findings from other hot gas giant exoplanets that also exhibit excess spectral line broadening (e.g., Cont et al. 2022; Lesjak et al. 2025; Zhang et al. 2026).

The noise scaling term β is close to one. This indicates a proper estimate of the CRIRES⁺ uncertainties.

Fig. 4 Atmospheric temperature and molecular abundance profiles. Left panel: median T-p profile (solid line) with its one-sigma uncertainty range (shaded area). The vertical dotted and dashed lines indicate the expected Teq of 51 Peg b, calculated for the two limiting cases of full and no dayside to nightside heat redistribution, respectively. The horizontal gray lines at 10−1.6 bar mark the retrieved lower limit on pcloud. In addition, we plot the condensation curves of Fe, Mg2SiO4, and MgSiO3, the most important cloud constituents expected in 51 Peg b’s atmosphere. Middle panel: retrieved molecular VMRs (solid lines) together with their one-sigma uncertainties (shaded areas). Both CO and H2O remain approximately constant throughout the entire pressure range. Right panels: histograms of the VMRs of the two molecular species, sampled at a representative pressure of 10−4 bar. The median and one-sigma uncertainty intervals are indicated by the dotted black lines.

6.2.2 Deviation from previously reported K_p values

The systemic velocity is derived as v_sys = $-{35.21}_{-2.35}^{+2.11}$ km s⁻¹, which is consistent with the value inferred from the stellar spectral lines in Sect. 4. The retrieved K_p is equal to ${102.8}_{-9.1}^{+8.3}$ km s⁻¹. Previous cross-correlation studies of 51 Peg b have reported higher K_p values than those obtained from our retrieval. Brogi et al. (2013) and Birkby et al. (2017) both found values around 133 km s⁻¹, which is approximately 30 km s⁻¹ higher than our results.

We tested if adopting an orbital phase shift (Δϕ) could align the signal of 51 Peg b to the K_p values reported by these previous studies. To this end, we conducted an additional retrieval run in which we fixed all parameters to their best-fit values and allowed only for an additional Δϕ as a free parameter. This retrieval resulted in a best-fit phase shift of Δϕ = −0.007 for aligning our signal with a K_p of 133 km s⁻¹. We verified whether the uncertainty on the orbital phase, resulting from the propagation of the uncertainties in the planetary time of conjunction (T_C) and P_orb, can allow for such a phase shift. The orbital phase uncertainty is ±0.0007, which is an order of magnitude smaller than the phase shift value required to match the K_p from the previous studies. Therefore, we rule out the possibility that uncertainties in T_C and P_orb are the cause of the observed K_p discrepancy.

Alternatively, a nonzero eccentricity could cause an orbital phase shift. The observations of Brogi et al. (2013), Birkby et al. (2017), and the present study probe different orbital phase intervals. Consequently, an eccentric orbit would produce distinct phase shifts in each dataset, leading to different offsets in the derived K_p values. However, we consider it unlikely that a mature hot Jupiter such as 51 Peg b has a significantly eccentric orbit. Indeed, in Sect. 4, we derived a two-sigma upper eccentricity limit at 0.0063, which is consistent with a circular orbit. However, we note that this does not rule out the very small eccentricity of 0.001 and the offset in the argument of periastron of 6 deg that would align the signal of Birkby et al. (2017) with the expected v_sys value.

Direct comparison with the results of Brogi et al. (2013) and Birkby et al. (2017) is complicated by the phase shifts they applied to align their planetary signals with the expected v_sys. We note that the phase shift of 0.1 in Birkby et al. (2017) is erroneously expressed in radians rather than in orbital phase units (confirmed by the authors; priv. comm.), corresponding to an orbital phase shift of 0.016, which is close to the 0.0095 phase shift applied by Brogi et al. (2013). Both previous studies find higher K_p values than retrieved in our work. Considering the results for individual nights in Brogi et al. (2013) (see their Table 1; K_p values of ${122.3}_{-24.8}^{+15.9}$ km s⁻¹ and ${134.7}_{-21.3}^{+13.7}$ km s⁻¹) and the one-sigma uncertainty interval of Birkby et al. (2017) (see their Fig. 8, with the one-sigma region extending down to a K_p of approximately 107 km s⁻¹), our result of K_p = ${102.8}_{-9.1}^{+8.3}$ km s⁻¹ is close to agreement at the upper uncertainty margin. However, our K_p cannot be fully reconciled with the combined result of all nights from Brogi et al. (2013). Part of this discrepancy may stem from the different orbital phases sampled: Brogi et al. (2013) observed closer to the extrema of 51 Peg b’s RV curve, where the sensitivity to K_p is highest, whereas our observations near superior conjunction are less sensitive to the parameter. Additionally, discrepancies between the studies may potentially arise from a different treatment of wavelength calibration and telluric removal methods, both of which are known caveats in high-resolution spectroscopy of exoplanet atmospheres. We therefore present our measurement as an alternative orbital solution, and note that resolving this discrepancy will require future observations.

To assess the robustness of our results, we performed an additional retrieval with K_p fixed to 133 km s⁻¹, close to the value reported by Brogi et al. (2013) and Birkby et al. (2017). The posterior distributions resulting from this additional test are shown in Fig. B.2. All retrieved parameters are consistent with those from our free K_p retrieval, except for v_sys, which is red-shifted by approximately 6 km s⁻¹ relative to the value derived in Sect. 4 from the stellar spectral lines. Such an offset in v_sys could in principle result from a strong dayside to nightside wind. However, we deem this scenario unlikely, as dayside to nightside circulation is expected at T_eq values higher than that of 51 Peg b, for which jet stream circulation is predicted (e.g., Komacek & Showman 2016; Showman et al. 2020). Furthermore, the retrieved spectral line broadening is slightly larger in the fixed K_p retrieval than in the free K_p retrieval, suggesting that K_p ~ 103 km s⁻¹ achieves a better alignment of the planetary spectral lines to a common rest frame.

6.2.3 Derived parameters

As described in Sect. 5.2, the retrieval parameters R_p and K_p enable the derivation of quantitative constraints on additional planetary properties. These include the surface gravity log g and planetary mass M_p. In addition, assuming the planet’s rotation axis is aligned with the inclination angle of its orbit, the orbital inclination can be derived from the retrieved K_p value via sin i = P_orbK_p/2πa, where P_orb is the orbital period and a the orbital semi-major axis. Consequently, the orbital inclination i represents a third parameter of 51 Peg b that can be constrained from the retrieval, in addition to log g and M_p. Figure B.3 shows the posterior distributions of the three derived parameters, all of which result in bounded constraints, yielding the following values: log g = 3.15 ± 0.12 log cm s⁻², M_p = ${0.61}_{-0.05}^{+0.06}$ M_Jup, and i = ${49.8}_{-5.7}^{+5.8}$ deg. For 51 Peg b, this is the first direct measurement of log g, with a value that is consistent with those typically found for other hot Jupiters (e.g., Southworth et al. 2007; Torres et al. 2008). The two parameters M_p and i were previously estimated by Brogi et al. (2013) and Birkby et al. (2017). The values of 0.46 M_Jup and 0.48 M_Jup in these studies differ from M_p determined in this work, due to the direct dependence of the parameter on K_p. Our derived orbital inclination is lower than the previously reported values of approximately 80 deg. Again, this discrepancy can be attributed to the dependence on K_p. Our inferred value of i = 50 deg is supported by the recent measurement of the stellar inclination of 51 Peg of ${53}_{-11}^{+25}$ deg (Metcalfe et al. 2024), under the assumption that the planetary orbit is aligned with the stellar spin axis. We note, however, that due to the large uncertainties, this measurement is also consistent with the inclination value reported by Birkby et al. (2017) and only marginally inconsistent with that of Brogi et al. (2013).

The determination of the orbital inclination is needed for understanding which regions in the planetary atmosphere are being probed and provides context to the retrieved atmospheric thermal, chemical, and dynamical properties. While the results of Brogi et al. (2013) and Birkby et al. (2017) suggest that around superior conjunction the planetary dayside is almost perfectly aligned with the observer’s line of sight, the lower inclination value derived in the present work indicates that a substantial fraction of the planetary nightside may still be visible during these orbital phases. Therefore, the atmospheric conditions obtained from our retrieval are likely not corresponding to pure dayside conditions, but rather to a combination of dayside and nightside conditions present on the hemisphere aligned with the observer’s line of sight. To resolve the discrepancy between previous studies and the present work in the derived M_p and i values, and to enable a more comprehensive understanding of 51 Peg b’s atmospheric properties, future observations covering a broader range of orbital phases will be required.

The derivation of K_p provides an interesting additional constraint on the stellar mass that is implicitly included in the calculation of the orbital parameters. We can write it in an explicit way if we use Kepler’s third law to derive mass constraints that depend only on the orbital period P_orb, eccentricity e, and the semi-amplitude velocities K_p and K_*: $$M_{\mathrm{p}}(\sin i{)}^3=\frac{P_{\mathrm{o}\mathrm{r}\mathrm{b}}}{2\pi G}{(1-e^2)}^{\frac{3}{2}}\frac{{(K_{\ast }+K_{\mathrm{p}})}^3}{1+\frac{K_{\mathrm{p}}}{K_{\ast }}},$$(5) $$M_{\ast }(\sin i{)}^3=\frac{P_{\mathrm{o}\mathrm{r}\mathrm{b}}}{2\pi G}{(1-e^2)}^{\frac{3}{2}}\frac{{(K_{\ast }+K_{\mathrm{p}})}^3}{1+\frac{K_{\ast }}{K_{\mathrm{p}}}}.$$(6)

From our results, we obtained M_p(sin i)³ = 0.27 M_Jup and M_*(sin i)³ = 0.48 M_⊙. This is consistent with the mass determined from spectroscopic analysis and isochrone models, M_* = 1.07 M_⊙, which also yields i = 50 deg and M_p = 0.61 M_Jup.

7 Conclusions

We have characterized the atmospheric properties and derived planetary parameters of the non-transiting hot Jupiter 51 Peg busing high-resolution spectroscopy in the K band with CRIRES⁺. In addition, we have used archival data from multiple high-stability, high-resolution spectrographs to refine the planetary ephemeris. By combining RV analysis methods, the cross-correlation technique, and a Bayesian retrieval framework, we were able to obtain the following key results:

1. Our RV analysis provides the most up-to-date estimates of the time of conjunction and orbital period of 51 Peg b. Moreover, we derived an upper limit on the orbital eccentricity that is consistent with a circular orbit (Sect. 4). The retrieved orbital semi-amplitude of 51 Peg b is lower than reported in previous cross-correlation studies (Sect. 6.2.2);

2. We confirmed the spectral signature of H₂O and found a possible indication for CO spectral lines in the planetary emission spectrum. The CO signal appears weaker because it originates from deeper atmospheric layers than the H₂O spectral lines and is therefore largely muted (Sect. 6.1);

3. We inferred a non-inverted T-p profile, consistent with previous studies of 51 Peg b and with the expected equilibrium temperature. This indicates the absence of strong optical absorbers capable of producing thermal inversions. The retrieval constrains the pressure level of an opaque cloud deck at the bottom of 51 Peg b’s atmosphere, which is consistent with the condensation curves of magnesium-silicate and iron clouds (Sect. 6.2.1);

4. Our retrieval yields a C/O ratio consistent with the solar value. The atmospheric metallicity retrieved under the assumption of chemical equilibrium is higher than typically expected for hot Jupiters. This elevated value may result from a quenched chemical structure of 51 Peg b’s atmosphere, which is caused by the uplift of heavy elements from the planet’s interior (Sect. 6.2.1);

5. We inferred the presence of spectral line broadening that cannot be explained by 51 Pegb’s tidally locked rotation alone. This excess broadening is most likely caused by super-rotation in the planetary atmosphere. The absence of an offset in the retrieved systemic velocity from the expected value disfavors a dayside to nightside wind flow as the main cause of the observed line broadening (Sect. 6.2.1);

6. Finally, we derived additional properties of 51 Peg b from the retrieved RV semi-amplitude, including the planetary surface gravity, mass, and orbital inclination. Given that our RV semi-amplitude differs from previous estimates, the derived mass is higher and the orbital inclination lower than reported in earlier studies (Sect. 6.2.3).

Beyond analyzing the atmosphere of 51 Peg b, the present study demonstrates that high-resolution spectroscopy is a powerful tool for investigating fundamental parameters – such as planetary mass or orbital inclination – that are otherwise difficult to access in non-transiting exoplanet systems. Several massive, short-period exoplanets that are not observable in transit still lack measurements of these parameters (e.g., Cullen & Bayliss 2024). Therefore, future high-resolution emission spectroscopy of non-transiting exoplanets will enable a significant expansion of the sample of planets with well constrained fundamental properties, as well as provide detailed atmospheric characterization.

Data availability

The extracted data are available at https://zenodo.org/records/17565920 (Cont & Lavail 2025).

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. This project is based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under the ESO programme 114.279X. D.C. is supported by the LMU-Munich Fraunhofer-Schwarzschild Fellowship and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2094 – 390783311. F.L. acknowledges support from the European Union (ERC-CoG, EVAPORATOR, Grant agreement No. 101170037). 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. 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. N.P., L.B.-Ch., and A.D.R. acknowledge support by the Knut and Alice Wallenberg Foundation (grant 2018.0192). O.K. acknowledges support by the Swedish Research Council (grant agreements no. 2023-03667) and by the Swedish National Space Agency. E.N. acknowledges the support from the Deutsches Zentrum für Luft- und Raumfahrt (DLR, German Aerospace Center) – project number 50OP2502. M.R. acknowledges the support by the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (DFG PR 36 24602/41). This work has made use of the following Python packages: Astropy (Astropy Collaboration 2013), CMasher (van der Velden 2020), corner (Foreman-Mackey et al. 2013), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), PyAstronomy (Czesla et al. 2019), and SciPy (Virtanen et al. 2020).

References

Appendix A Additional figures

Fig. A.1 Radial velocity offsets before (top panel) and after (bottom panel) the alignment to the molecfit wavelength solution. Each data point represents an individual spectrum and wavelength segment. Points corresponding to the A and B nodding positions are shown in red and blue, respectively. The ±0.1 km s−1 interval around the average A and B wavelength solutions is indicated by the gray dashed lines.

Fig. A.2 Relative change in the standard deviation of the residual spectral matrix between consecutive SYSREM iterations. After six iterations, a plateau is reached, indicating that SYSREM no longer identifies significant correlated trends.

Fig. A.3 Same as in Fig. 2 but without masking the stellar CO residuals of the CCF maps when computing the S/N maps. For CO, the stellar residuals introduce strong artifacts into the S/N map. For the combined CO and H2O signal, the artifacts are weaker, but still present, given the reduced S/N value compared to Fig. 2.

Appendix B Posterior distributions

Fig. B.1 Retrieval results for 51 Peg b. Corner plot: Posterior distributions and correlations between the retrieval parameters. The dashed vertical lines in each posterior distribution indicate the median and one-sigma intervals (16th, 50th, and 84th percentiles) for bounded parameters. For parameters with lower or upper limits, the two-sigma interval (5th or 95th percentile) is reported. Top-right panel: Median temperature curve with uncertainty intervals. The vertical dotted and dashed lines indicate the expected Teq value of the planet, for the cases of full and no dayside to nightside heat redistribution, respectively. The horizontal gray lines indicate the retrieved lower pressure limit pcloud of an optically thick cloud deck.

Fig. B.2 Same as Fig. B.1, but with Kp fixed to 133 km s−1. All parameters except the systemic velocity are consistent with the free Kp retrieval. We note that that the spectral line broadening is slightly higher than in the free Kp retrieval. This indicates that a better alignment of the planetary spectral lines is achieved with the Kp value of approximately 103 km s−1 obtained by the free retrieval.

Fig. B.3 Same as Figs. B.1, but for the derived planetary parameters. To illustrate the correlations between the free parameters Rp and Kp, from which the derived parameters are obtained, we also include these free parameters in the figure.

All Tables

All Figures

	Fig. 1 Radial velocity curve of 51 Peg with the updated orbital solution. We show the phase-folded RV data points and the best-fit Keplerian model.
In the text

Fig. 2 Model spectra, CCF maps, and S/N maps of CO, H2O, and the combination of both species. Left panels: model spectra normalized to the spectral continuum level. Wavelengths covered by the CRIRES+ K2166 setting used in this work correspond to the gray shaded area. Middle panels: cross-correlation function maps. The signal of H2O and the signal obtained from a combination of CO and H2O spectral lines can be identified as diagonal trails in the mappings, indicated by the solid white lines. A residual CO signal, which is consistent with the stellar rest frame and could not be entirely removed with SYSREM, is masked. Right panels: signal-to-noise maps. These maps reveal no clear CO signal but a detection of spectral H2O lines. The position of the planetary signal (Kp value reported by Brogi et al. (2013) and Birkby et al. (2017); vsys determined via the RV analysis in Sect. 4) is indicated by the dashed white lines. We show the CCF maps and S/N maps obtained after six consecutive SYSREM iterations. A non-masked version of the CO and the species-combined CCF and S/N maps is given in Fig. A.3.

In the text

Fig. 3 Contribution functions of the CO, H2O, and species-combined model spectra. The CO signal originates from deeper regions of the atmosphere, whereas the H2O lines form at higher altitudes. Consequently, in the species-combined model, the CO spectral lines are largely muted by the overlying H2O signal. The contribution of CO to the combined signal becomes most evident in the region around the CO bandhead at approximately 2.3 μm.

In the text

Fig. 4 Atmospheric temperature and molecular abundance profiles. Left panel: median T-p profile (solid line) with its one-sigma uncertainty range (shaded area). The vertical dotted and dashed lines indicate the expected Teq of 51 Peg b, calculated for the two limiting cases of full and no dayside to nightside heat redistribution, respectively. The horizontal gray lines at 10−1.6 bar mark the retrieved lower limit on pcloud. In addition, we plot the condensation curves of Fe, Mg2SiO4, and MgSiO3, the most important cloud constituents expected in 51 Peg b’s atmosphere. Middle panel: retrieved molecular VMRs (solid lines) together with their one-sigma uncertainties (shaded areas). Both CO and H2O remain approximately constant throughout the entire pressure range. Right panels: histograms of the VMRs of the two molecular species, sampled at a representative pressure of 10−4 bar. The median and one-sigma uncertainty intervals are indicated by the dotted black lines.

In the text

Fig. A.1 Radial velocity offsets before (top panel) and after (bottom panel) the alignment to the molecfit wavelength solution. Each data point represents an individual spectrum and wavelength segment. Points corresponding to the A and B nodding positions are shown in red and blue, respectively. The ±0.1 km s−1 interval around the average A and B wavelength solutions is indicated by the gray dashed lines.

In the text

	Fig. A.2 Relative change in the standard deviation of the residual spectral matrix between consecutive SYSREM iterations. After six iterations, a plateau is reached, indicating that SYSREM no longer identifies significant correlated trends.
In the text

	Fig. A.3 Same as in Fig. 2 but without masking the stellar CO residuals of the CCF maps when computing the S/N maps. For CO, the stellar residuals introduce strong artifacts into the S/N map. For the combined CO and H2O signal, the artifacts are weaker, but still present, given the reduced S/N value compared to Fig. 2.
In the text

Fig. B.1 Retrieval results for 51 Peg b. Corner plot: Posterior distributions and correlations between the retrieval parameters. The dashed vertical lines in each posterior distribution indicate the median and one-sigma intervals (16th, 50th, and 84th percentiles) for bounded parameters. For parameters with lower or upper limits, the two-sigma interval (5th or 95th percentile) is reported. Top-right panel: Median temperature curve with uncertainty intervals. The vertical dotted and dashed lines indicate the expected Teq value of the planet, for the cases of full and no dayside to nightside heat redistribution, respectively. The horizontal gray lines indicate the retrieved lower pressure limit pcloud of an optically thick cloud deck.

In the text

Fig. B.2 Same as Fig. B.1, but with Kp fixed to 133 km s−1. All parameters except the systemic velocity are consistent with the free Kp retrieval. We note that that the spectral line broadening is slightly higher than in the free Kp retrieval. This indicates that a better alignment of the planetary spectral lines is achieved with the Kp value of approximately 103 km s−1 obtained by the free retrieval.

In the text

	Fig. B.3 Same as Figs. B.1, but for the derived planetary parameters. To illustrate the correlations between the free parameters Rp and Kp, from which the derived parameters are obtained, we also include these free parameters in the figure.
In the text