© The Authors 2024

1 Introduction

The growing group of directly imaged super-Jupiters at wide separations continue to pose a challenge to planet formation models. It is unclear whether these planets form bottom-up via core accretion or represent the tail end of star formation and form top-down via gravitational instability. Occurrence rates from large surveys have started to provide constraints on the population as a whole (Nielsen et al. 2019; Vigan et al. 2021), offering evidence of multiple formation pathways. However, understanding the formation and evolution history of individual objects remains a challenge. Furthermore, there is still a great deal of uncertainty around the chemical and physical processes in the atmospheres of these planets, namely, the possible presence and composition of (patchy) clouds or disequilibrium chemistry (e.g., Rajan et al. 2017; Mollière et al. 2020). Since we are able to spatially resolve them from their host star, these young, self-luminous planets at wide separations are well suited for atmospheric characterizations, while establishing links between their atmospheric properties and their formation and migration histories.

Multiple tracers of the planet’s formation and migration history have been proposed, which are accessible with present-day observations. First, the elemental and isotopic abundance ratios in the planet’s atmosphere, such as the carbon-to-oxygen (C/O) ratio, have been advocated as useful tracers of where the planet has accreted its gas in the protoplanetary disk (Öberg et al. 2011; Madhusudhan et al. 2017; Zhang et al. 2021a). After measuring such abundance ratios, we can attempt to invert this problem and constrain the planet’s formation and migration history (GRAVITY Collaboration 2020; Mollière et al. 2022). The C/O ratio has been measured for a growing number of directly imaged planets and brown dwarf companions (e.g. Konopacky et al. 2013; Mollière et al. 2020; GRAVITY Collaboration 2020; Wilcomb et al. 2020; Wang et al. 2021b,c, 2023; Petrus et al. 2021; Hoch et al. 2022; Xuan et al. 2022; Palma-Bifani et al. 2023; Brown-Sevilla et al. 2023, see Hoch et al. 2023 for an overview) as well as hot Jupiters (e.g., Madhusudhan et al. 2011; Line et al. 2021; Changeat et al. 2022). While low-resolution spectroscopy generally results in large uncertainties due to degeneracies with, for instance, clouds and surface gravity, medium-to-high-resolution spectroscopy has resulted in accurate and robust estimates of the C/O ratio for a handful of directly imaged planets (e.g., Ruffio et al. 2021; Wang et al. 2021b; Zhang et al. 2021a). Second, the planets’ spin rotation data contain information on the evolution of their angular momentum. This rotation rate may deviate depending on the formation channel or through interactions with the circumplanetary disk. Bryan et al. (2018) found no distinction in the distribution of spin measurements for planetary mass companions and isolated brown dwarfs, indicating that their spin evolution is regulated in the same way. On the other hand, Wang et al. (2021b) found a tentative trend in increasing rotation rate with decreasing companion mass, which could indicate that magnetic braking is less efficient for lower mass companions. While the rotation rates of a handful of planetary mass companions have been measured (Snellen et al. 2014; Schwarz et al. 2016; Bryan et al. 2018; Xuan et al. 2020; Wang et al. 2021b), a larger sample and more accurate measurements are needed to confirm this trend. Finally, the orbit of the planet can reveal its dynamical history and whether it has, for example, experienced a scattering event. High-resolution spec-troscopy allows us to measure the radial velocity of the planet, which can be used to break degeneracies in orbital fits (e.g., Schwarz et al. 2016; Ruffio et al. 2021).

High-resolution spectrographs equipped with adaptive optics (AO) provide the unique capability of measuring robust elemental abundance ratios, spin rotation, and radial velocity for directly imaged planets. One of the best targets for this is the emblematic β Pictoris b (Lagrange et al. 2010) due to its close proximity. Its spin rotation was estimated using the original CRIRES by Snellen et al. (2014), showing a fast rotation rate of v sin i = 25 ± 3 kms⁻¹. Many works have characterized its atmosphere (e.g., Bonnefoy et al. 2013; Currie et al. 2013; Morzinski et al. 2015; Baudino et al. 2015; Hoeijmakers et al. 2018; Stolker et al. 2020), with most studies finding effective temperatures of 1650–1800 K, low surface gravities, and thick clouds. The most detailed study to date was performed by GRAVITY Collaboration (2020), who report a sub-solar C/O ratio of 0.43 ± 0.05 from atmospheric retrievals on VLTI/GRAVITY K-band data in combination with low-spectral resolution YJH data from GPI (Chilcote et al. 2017). Assuming that the host star β Pictoris has a solar C/O ratio, the authors concluded that the planet has most likely formed via core accretion.

In this work, we report on observations with the newly refurbished VLT/CRIRES+ instrument to probe the atmosphere of β Pictoris b and provide an updated measurement of its spin rotation and radial velocity. Section 2 describes the observations and data reduction pipeline. Section 3 details the data modeling of the observations and Sect. 4 describes the atmospheric models of the planet. Section 5 reports the detection of molecules and results from atmospheric retrievals, which are discussed in Sect. 6. Finally, our conclusions are summarized in Sect. 7.

2 Observations and data reduction

2.1 Observations

We observed the β Pictoris system for 40 min of integration time on both 11 November 2021 and 13 November 2021 with the upgraded Cryogenic High-Resolution Infrared Echelle Spectro-graph (CRIRES+; Dorn et al. 2014, 2023) as part of program 108.22HG.001. CRIRES+ is located at the Nasmyth focus B of UT 3 at the VLT and has been upgraded with a cross-disperser and three Hawaii 2RG detectors, increasing the simultaneous wavelength coverage by up to a factor ten. We used wavelength setting K2166 and a slit width of 0.2″, resulting in a spectral resolving power of R ~ 100000. Throughout this work, we do not include the first two of the seven spectral orders, as these are heavily dominated by tellurics. This results in an effective wavelength coverage of 2.06–2.47 µm, but with gaps in between due to the different spectral orders. The slit was placed in such a way that it encompassed both the host star and P Pictoris b. The presence of the star on the slit significantly limited the maximum exposure time in order to avoid saturation. At the time of observation the planet was at a separation of 0.48” and a position angle of 31.25 degrees (Wang et al. 2021a). The observations were taken using ABBA nodding with DIT=15s and NDIT=8, resulting in 10 frames in both nodding positions for both nights. The observations were taken with adaptive optics from MACAO and the conditions were decent with seeing ranging from 0.55″ to 1″ during the observations.

2.2 Data reduction with pycrires

The raw data were reduced using the pycrires¹ package for python (Stolker & Landman 2023). This package contains python wrappers for the CRIRES+ EsoRex pipeline and custom functions specifically developed for direct observations of sub-stellar companions. We used the EsoRex pipeline version 1.2 to perform the basic image processing, including dark and flat field correction and an initial wavelength solution using the uranium-neon lamp and Fabry–Perot etalon calibration files (Seemann et al. 2014). We then used the intermediate data products from the obs_nodding function, which was run on each nodding pair separately, as input to our custom functions. This data product is shown in the top panel of Fig. 1. We subsequently correct for the curvature of the slit on the detector and the slit tilt using interpolation and the information from the trace-wave table produced by the EsoRex pipeline to obtain 2D rectified spectra, as shown in the central panel of Fig. 1. We found that there were slight vertical shifts of about a pixel between different spectral orders in these 2D rectified spectra, likely due to slight errors in the trace determination. To remove this, we fit a Gaussian to the slit illumination, which was calculated by summing the 2D rectified spectra along the spectral dimension. This was done for each spectral order and exposure and we subsequently aligned all the orders using interpolation.

We found the wavelength solution from the EsoRex pipeline to be not sufficiently accurate, so we performed an additional correction on the wavelength solution in a similar way as Holmberg & Madhusudhan (2022). We fitted for a second order polynomial correction to the current best wavelength solution, which is initially obtained from the CRIRES+ pipeline. This was done by maximizing the cross-correlation between the extracted and continuum removed spectrum of β Pic and a template of the telluric absorption generated with SkyCalc (Noll et al. 2012) on a grid of polynomial coefficients. We repeated this process three times and gradually increased the accuracy of the grid around the current best wavelength solution. We found that this resulted in accurate wavelength solutions, even in orders where there are minimal tellurics, which we checked by visual comparison of the resulting spectrum and a telluric model, and through inspection of the resulting cross-correlation functions. An example of the observed spectrum at the planet’s location overlayed with the telluric model is shown in the bottom panel of Fig. 1.

Fig. 1 Visualization of the data extraction pipeline. Top row shows the product of the obs_nodding recipe from the ESO CRIRES+ data reduction pipeline. The middle row shows the rectified 2D spectrum corrected for the slit tilt and curvature using pycrires and the bottom row shows the spectrum at the location of the planet, plotted together with a telluric model generated with SkyCalc.

3 Analysis framework

While the adaptive optics performance was good, the signal at the location of the planet is dominated by the stellar contribution. From the reduced data, we measure the stellar flux to be a factor 100 lower at the location of the planet position compared to the peak stellar flux, as shown in the raw broad-band contrast curve in Fig. 2. Given the star-planet contrast of ΔK = 9.1 (Chilcote et al. 2017), this means that the planet only contributes about ~2% of the flux at its position. It is therefore not straightforward to extract the planet spectrum and derive its atmospheric parameters. In this section, we explain our data modeling, which closely follows the forward modeling approach developed in Ruffio et al. (2019, 2021) for OSIRIS data of the HR8799 planets. Even though the framework from Ruffio et al. (2019) was developed for medium spectral resolution data from an integral field spectrograph, it is almost directly applicable to the high-resolution long-slit data from CRIRES+.

3.1 Data modeling

The reduced ID spectrum, d, at a specific position along the slit can be written as a vector with (as entries) the flux at each wavelength λ. This signal is a combination of the stellar signal, d_s; planet signal, d_p; systematics, d_sys; and noise, η. It can be written as: $$d=d_s+d_{\text{p}}+d_{\text{sys}}+\eta .$$(1)

To derive parameters of the planet’s atmosphere, we forward-modeled these contributions. We expressed the observed data as a linear combination of stellar, planetary, and systematic components, closely following Ruffio et al. (2019): $$d=M_{\psi }c+\eta ,$$(2)

where M_ψ is the linear model matrix with nonlinear parameters, ψ, and c is a vector with entries corresponding to the amplitude of each linear component, and η is the noise. The different components that make up M_ψ are individually explained in the next subsections.

Fig. 2 Raw white light contrast curve for our CRIRES+ observations. The location of β Pictoris b is indicated in this diagram according to the K-band contrast from Chilcote et al. (2017).

3.1.1 Stellar contribution

We estimated the telluric imprinted stellar master spectrum, f_s, using the on-axis observation of the star. To reduce the influence of noisy datapoints, we estimated this stellar master spectrum by taking the sum over the five rows centered on the host star, which are free from planet signal. However, because speckles move outwards for increasing wavelengths, this stellar contribution is modulated by a low-order function a(λ) off-axis: $$d_{\text{s}}=c_{\text{s}}\alpha (\lambda )f_{\text{s}},$$(3)

where c_s is a unitless linear scaling parameter, part of c in Eq. (2), which we fit for. One can choose how to parameterize the modulation, a(λ). In principle, we can retain the continuum from the planet by jointly fitting the planet signal and this modulation (Ruffio et al. 2023). However, since we are completely dominated by the stellar contribution, which also has almost no spectral lines, this is not possible in our case. Instead, we estimated a directly from the data using: $$\widehat{\alpha }(\lambda )=\frac{ℒ\left[d_{\text{s}}\right]}{ℒ\left[f_{\text{s}}\right]},$$(4)

where ℒ refers to a low-pass filtering operation. However, we do not have access to the isolated contribution of the star, d_s, at the location of the companion, but only the total signal, d. Rewriting

Eq. (1), while ignoring the systematics, and plugging this into Eq. (4) we have: $$\widehat{\alpha }(\lambda )=\frac{ℒ\left[d-d_{\text{p}}\right]}{ℒ\left[f_{\text{s}}\right]}.$$(5)

For a linear low-pass filter we can split these terms up as: $$\widehat{\alpha }(\lambda )=\frac{ℒ\left[d-d_{\text{p}}\right]}{ℒ\left[f_{\text{s}}\right]}.$$(6)

In other words, we have to correct for the planet continuum leaking into the estimate of α, which effectively functions as a highpass filter on the planet model. This correction on the planet signal was included in the forward model of the planet signal, which is discussed in the next subsection and ignored here. We used a second-order Savitsky–Golay filter with a kernel width of 301 pixels as our low-pass filter.

Additionally, because CRIRES is a slit spectrograph, the line spread function (LSF) can change along the slit. For example, in good seeing conditions with adaptive optics, the full width at half maximum (FWHM) of the point spread function (PSF) is smaller than the width of the slit. This means that on-axis the spectral resolution of the stellar contribution is determined by the width of the PSF, while off-axis this is determined by the width of the slit. We account for this effect by including shifted versions of the stellar model in the linear model matrix, which allows for flexibility in fitting the LSF of the stellar contribution.

Considering both effects, we included the following components in M_ψ for the stellar contribution to the observed data: $$m_{\text{s},\text{k}}=\frac{ℒ[d]}{ℒ\left[f_{\text{s}}\right]}{f}_{\text{s}}[\lambda +k\Delta \lambda ],$$(7)

where f_s[λ + kΔλ] refers to the stellar master spectrum shifted by k pixels to account for the changing LSF along the slit. We used integer shifts for k from −3 to 3 wavelength bins throughout this work.

3.1.2 Planet contribution

The observed planet signal can be written as: $$d_{\text{p}}=c_{\text{p}}T(\lambda )F_{\text{p}}(\lambda ,\psi )$$(8)

where T(λ) is the transmission of both the atmosphere and the instrument and includes the unit conversion from flux density to detector counts; F_p(λ, ψ) is a model of the planetary spectrum with nonlinear parameters, ψ and c_p is a unitless linear scaling parameter that is fitted for. The planetary emission model and its parameterization are described in Sect. 4. We estimated the transmission using the stellar master spectrum and a PHOENIX model of β Pic (F_s(λ)): $$T(\lambda )=\frac{f_{\text{s}}}{F_{\text{s}}(\lambda )}$$(9)

Considering that the continuum of the planet model leaks into the estimate of α (Eq. (6)) and, subsequently, the stellar model, we have to correct for this in our model of the planet contribution. Effectively, this applies a high-pass filter to the planet contribution and leads to following corrected planetary contribution, d′_p: $$d_{\text{p}}^{\prime }=d_{\text{p}}-\frac{ℒ\left[d_{\text{p}}\right]}{ℒ\left[f_{\text{s}}\right]}{f}_{\text{s}}=c_{\text{p}}\left(T(\lambda )F_{\text{p}}(\lambda ,\psi )-\frac{ℒ\left[T(\lambda )F_{\text{p}}(\lambda ,\psi )\right]}{ℒ\left[f_{\text{s}}\right]}{f}_{\text{s}}\right).$$(10)

We therefore included the following linear component for the corrected planet signal in our model matrix M_ψ: $$m_{\text{p}}(\psi )=T(\lambda )F_{\text{p}}(\lambda ,\psi )-\frac{ℒ\left[T(\lambda )F_{\text{p}}(\lambda ,\psi )\right]}{ℒ\left[f_{\text{s}}\right]}{f}_{\text{s}}.$$(11)

We note that we also have to correct for the planet continuum leaking into the shifted versions of the stellar spectrum used for modeling the LSF (Eq. (7)). However, since the amplitude of the central stellar component dominates over the amplitudes corresponding to the shifted versions, we ignore this effect throughout this work.

3.1.3 Systematics and noise

After an initial subtraction of the stellar contribution from each row, we found from our visual inspection that there still was some correlated structure in the data. This may, for example, originate from the imperfect removal of the stellar and telluric contamination or other systematics. To remove most of these correlated residuals, we included a model of the systematics in our forward model: $$d_{\text{sys }}={\displaystyle \sum _i^N{c}_i}{m}_{\text{sys },i},$$(12)

where c_i are the linear coefficients and m_sys,i the systematics modes. This systematics model was computed the same way as in Wilcomb et al. (2020) and Ruffio et al. (2021). First, we subtracted the estimated stellar contribution from each row by fitting the model without the planet contribution. Subsequently, we scaled the residuals of each row to the same local continuum by dividing by the low-pass filtered data. The systematics were then determined from these residuals using a principal component analysis (PCA), while excluding the four spaxels around the considered position of the companion. This was done to avoid the planet signal leaking into the systematics model. Finally, the PCA components were scaled back to the local continuum at the considered location. We used ten PCA components in the systematics model throughout this work. To remove any residual outliers, we did an initial subtraction of the star and systematics model and masked 4σ outliers. This leaves us with mainly uncorrelated noise, which we checked by inspecting the autocorrelation function of the residuals.

We also found the noise estimates from the pipeline to not be sufficiently accurate, which is crucial for obtaining realistic uncertainties on retrieved atmospheric parameters. Instead, we estimated the noise as follows: For each exposure, the stellar and systematics models were constructed and subtracted from the observed data. The noise was then estimated by taking the standard deviation over time of these residuals. These uncertainties are used as diagonal entries for the covariance matrix Σ₀. Following Ruffio et al. (2019), we then include a scaling parameter, s, to account for a possible underestimation of these uncertainties: Σ = s²Σ₀, which is fitted for using the data.

3.2 Likelihood

Following Ruffio et al. (2019), we can write our forward model in matrix form: $$d=M_{\psi }c+\eta ,$$(13)

where the model matrix M_ψ is given by: $$M_{\psi }=\left[m_{\text{p}}\left(\psi \right),m_{\text{s-3}\Delta \lambda },\dots ,m_{\text{s,3}\Delta \lambda },m_{\text{sys,0,}}\dots ,m_{\text{sys,10}}\right],$$(14)

c are the linear coefficients and η is again the uncorrelated noise. Following the derivation from Ruffio et al. (2019), which uses a prior of P(s) ∝ s^−γ for the noise scaling parameter, s, we can analytically marginalize over the linear parameters c and noise scaling s. This results in the following posterior: $$L(\psi \mid d)=\frac{P(\psi )}{\sqrt{\left|{\Sigma }_0\right|\times \left|M_{\psi }^T{\Sigma }_0^{-1}{M}_{\psi }\right|}}{\left(\frac{1}{{\chi }_0^2}\right)}^{\frac{N_{\text{d}}-N_{\text{c}+\gamma -1}}{2}},$$(15)

with P(ψ) as the prior on the planetary parameters, ψ; N_d is the number of data points; N_c is the number of components in the linear model, and $${\chi }_0^2={\left(d-M_{\psi }\widehat{c}\right)}^T{\Sigma }_0^{-1}\left(d-M_{\psi }\widehat{c}\right),$$(16)

where $\widehat{c}$ is the linear least squares solution such that: $${\widehat{c}}^T{M}_{\psi }^T{\Sigma }_0^{-1}{M}_{\psi }=d^T{\Sigma }_0^{-1}{M}_{\psi },$$(17)

which we find using a numerical least squares solver. For details on the derivation we refer the reader to Appendix D from Ruffio et al. (2019). We used γ = 2 and calculated the likelihood for each nodding position and night separately and added the log-likelihoods. We did this to avoid additional interpolation of the data, as the wavelength solution is slightly different for the different nodding positions and nights. Furthermore, we applied Eq. (15) to each spectral order separately to avoid flux calibration errors between spectral orders having a major impact on the result. This means we are effectively fitting for the continuum of the planet. While this discards information on the continuum, we found that this can lead to more robust results with respect to particular choices in the data reduction. This also means that the noise scaling factor is allowed to deviate per spectral order.

4 Atmospheric models

4.1 Nominal model

We used petitRADTRANS (Mollière et al. 2019,2020) to generate model spectra of β Pictoris b. We followed a similar modeling approach as in, among others, GRAVITY Collaboration (2020) and Zhang et al. (2021a,b). Our nominal model uses a free pressure-temperature (P–T) profile, equilibrium chemistry, and the cloud model from Mollière et al. (2020). The P–T profile consists of four free knots located at 1 mbar, 0.1 bar, 1 bar, and 10 bar, respectively. The full P–T profile was then calculated using third-order B-spline interpolation in log-space of pressures and temperature. A prior on the temperature was chosen such that it smoothly decreases towards lower pressures and no inversions are possible. The chemical equilibrium model is detailed in Mollière et al. (2017), which takes as input the

P–T profile, metallicity [Fe/H], C/O ratio, and quenching pressure, P_quench, and returns the abundances of all the species. We include Rayleigh scattering from H2, He, collisional induced absorption of H₂–H₂ and H₂–He and line-by-line opacities of H₂O (Polyansky et al. 2018), CO (Rothman et al. 2010) and CH₄ (Hargreaves et al. 2020). The cloud model is presented in Mollière et al. (2020) and is based on the EddySed model from Ackerman & Marley (2001). We included opacity sources of Fe and MgSiO₃ clouds, assuming crystalline, irregularly-shaped particles, and account for scattering in the radiative transfer. We chose these cloud species as they are likely the most prominent at the expected temperature of β Pictoris b and were retrieved in its atmosphere by GRAVITY Collaboration (2020). Additionally, we fitted for the radial velocity (RV) of the planet and the line broadening caused by the spin of the planet. For the rotation of the planet we used the fastRotBroad function from pyAs-tronomy (Czesla et al. 2019), which has as its free parameters: the projected spin velocity, υ sin(i), and linear limb darkening coefficient, ∈.

Initially, we found the posteriors of the surface gravity, metallicity and P–T profile to be strongly correlated and unstable to changes in the modeling or data reduction, which influenced the constraints on the other parameters. This issue of trying to constrain the surface gravity from just K-band spec-troscopy has been noted in other works (Zhang et al. 2021a; Palma-Bifani et al. 2023). To relieve this issue we use a Gaussian prior on the surface gravity based on the measurement of the dynamical mass from Lacour et al. (2021) of 11.9 ± 3.0 M_Jup and chose a radius prior of 1.4 ± 0.1 R_Jup based on the values of Chilcote et al. (2017); Stolker et al. (2020); GRAVITY Collaboration (2020), which gives log(g) = 4.18 ± 0.13. The priors on the parameters are summarized in Table 1.

Fig. 3 Detection of the planet using model templates of the full planet, H2O and CO using cross-correlation (top row) and the likelihood ratio (bottom row). The radial velocity is shown here in the barycentric restframe. In the top panel we also show the cross-correlation value for a reference spectrum, which was taken at the same separation on the other side of the star. In the bottom panels we indicate the location of the star and the planet and the expected velocity of the planet.

4.2 Alternative models

To test the robustness of the results to model assumptions, we also ran retrievals with alternative models. We tested the following alternative models:

Free composition. Same as the nominal model, but now with free composition, where we fit for the logarithm of the abundance of the species, which is assumed to be constant vertically in the atmosphere. We also include minor species (CO₂ , HCN, NH₃ and ¹³CO) in these retrievals. This is done to search for the presence of these species in a more robust way than the crude cross-correlation analysis (detailed in Sect. 5.1) and allow us to obtain upper limits on the abundance of these species. The retrieved abundances are subsequently converted to a C/O ratio and metallicity.

Forced grey cloud deck. In Sect. 5, we show that the nominal model lead to cloud-free solutions, which is inconsistent with the low spectral resolution data at shorter wavelengths (Chilcote et al. 2017; GRAVITY Collaboration 2020). To test whether this has impact on the values of the other retrieved parameters we include a strong prior on the cloudiness of the planet. Instead of determining the cloud location and opacity from the condensation curves, we follow a similar approach as in Burningham et al. (2017) and fit for the location and opacity of the cloud deck, for which we chose priors such that the cloud is optically thick and in the photosphere of the planet. This alleviates the uncertainty in the composition of the clouds, which can impact the location of the cloud deck.

Strong prior based on low-resolution data. We use the posteriors obtained by GRAVITY Collaboration (2020) as the prior here and fix the P–T profile to the one obtained in GRAVITY Collaboration (2020). This is similar to jointly fitting the GPI, GRAVITY, and our new CRIRES+ data with equal weighting for the combined low-resolution data and the high-resolution data.

5 Results

5.1 Detection of molecules

First, we aim to detect molecules in the atmosphere of the planet. We do this by calculating the log-likelihood (as defined in Eq. (15)) as a function of position along the slit and radial velocity. The templates used are generated with the nominal model discussed in Sect. 4 with the parameters found in GRAVITY Collaboration (2020) and υ sin(i) = 25 km s⁻¹ (Snellen et al. 2014). In the case of the full model, we compare the likelihood ratio with respect to a model not including any planet signal, so only the stellar model and systematics. For the molecules, we evaluated the likelihood ratio between the full model and the likelihood without the considered species. For example, for water, we have: $$\Delta \ln L_{{\text{H}}_2\text{O}}=\ln L_{\text{full }}-\ln L_{\text{without }}{\text{H}}_2\text{O. }$$(18)

The results are shown in Fig. 3, showing clear detections of H₂O and CO. The planet signal is predominantly detected in two spaxels at a roughly equal strength, which we simply summed for the remainder of the analysis. We also show the cross-correlation function (CCF) for the combined planet signal following Ruffio et al. (2019), given by the linear parameter corresponding to the planet component in the model while including only the specified species in the planet model. We normalized the CCF by the standard deviation of the CCF of the combined two spax- els at the same separation as the planet but on the other side of the star. Using this, we estimate detection signal-to-noise ratios (S/Ns) of 24.5, 21.0, and 8.4 for the full model, water, and carbon monoxide, respectively. We also searched for features of CO2 , CH4 , HCN, and NH3 and we did not find any significant signal, as expected from equilibrium chemistry and the S/N of our observations. The high S/N detections are a major improvement over the results from Snellen et al. (2014), who were unable to detect water using the original CRIRES. In fact, while not shown here, we are able to detect the planet confidently (S /N > 5) in almost all individual exposures of 2 min. This improvement in data quality is mainly due to significantly increased instantaneous wavelength coverage of the instrument and improved detectors. A visual comparison between our cross-correlation map and the one from Snellen et al. (2014) is shown in Fig. A.1.

Fig. 4 Constraints on a subset of the parameters from the atmospheric retrievals using the models discussed in Sect. 4. The RV shown here is in the barycentric restframe.

5.2 Atmospheric retrievals

We ran atmospheric retrievals on the CRIRES+ data for the models specified in Sect. 4. We sampled the posterior using PyMultiNest (Buchner 2016), a Python wrapper for MultiNest (Feroz et al. 2019), with 500 live points and an initial sampling efficiency of 80%, which is the recommended value for parameter estimation. The obtained posterior distribution for a selection of parameters is listed in Table 2 and shown in Fig. 4. The full corner plot for the nominal model is shown in Fig. B.1. The values of the retrieved parameters are discussed in the next section. The retrieved P–T profile for the nominal model is shown in Fig. 5, together with the condensation curves for the best fitting parameters.

To test whether our models provide a good fit to the data, we show the spectral order with the highest detection significance together with the fitted model in Fig. 6. While the data are noisy at the native resolution, the bottom panel in Fig. 6 shows the data smoothed to a spectral resolution of 15 000, which is the intrinsic resolution of the planet signal set by its rotational velocity. This shows good agreement between the data and our model, allowing us to clearly identify individual absorption lines. Finally, we compare the spectral energy distribution of the retrieved models with the low spectral resolution Y-, J-, and H-band data from GPI (Chilcote et al. 2017) and K-band data from GRAVITY Collaboration (2020) in Fig. 7. The data from the different bands are scaled according to the retrieved scaling parameters from GRAVITY Collaboration (2020). The models are generated using the correlated-k mode of petitRADTRANS and include additional opacities from TiO, FeH, VO, Na and K. Since we did not fit for the radius, the models were scaled to match the K-band flux of the GPI data. The figure shows that all models provide a good fit to the data in the K-band, but that there is a large deviation at shorter wavelengths.

Fig. 5 Retrieved P–T profile for the nominal model with associated uncertainties. Also shown: the P–T profile obtained by GRAVITY Collaboration (2020), the condensation curves of Fe and MgSiO3 for the best fitting parameters, and the emission contribution function for the best fitting parameters.

Fig. 6 CRIRES+ data from one of the spectral orders along with the fitted model. Top panel shows the total extracted flux at the location of the planet with the full fitted model including star, planet, and systematics. The second panel shows the data after subtracting the stellar and systematic components, along with the best-fitting planet model. The third panel shows the residuals of this fit. Finally, the bottom panel shows the same as the second panel but then smoothed to a spectral resolution of 15 000 to show the broad planet features more clearly.

Fig. 7 Comparison of the spectral energy distributions of our retrieved models compared to the low-resolution GPI data from Chilcote et al. (2017) and GRAVITY data from GRAVITY Collaboration (2020). Since we do not retrieve radius, the flux of our models is scaled such that it matches the GPI flux in the K-band. The solid line shows the median value and the colored region shows the 16th and 84th percentiles.

6 Discussion

6.1 Abundance ratios and metallicity

For the nominal model, we obtained a sub-solar C/O ratio of 0.41 ± 0.04, which is consistent with the values obtained by GRAVITY Collaboration (2020). The alternative models result in similar constraints, albeit slightly higher for cloudy models. For the free composition model, the measured C/O ratio is the ratio in the gas phase only. In the presence of condensation, a part of the oxygen will be sequestered into clouds and the true C/O ratio will be even lower. While the uncertainty on the C/O ratio is likely underestimated, the retrieved value was found to be relatively robust to changes in the data reduction parameters and modeling choices. This illustrates the power of high-resolution spectroscopy for determining robust elemental abundance ratio’s of exoplanets. The interpretation of this sub-solar C/O ratio is limited by the fact that the oxygen abundance in the host star is unknown. It is therefore hard to make any statements on the formation and migration history of the planet. In the case the C/O ratio of β Pictoris is solar, GRAVITY Collaboration (2020) showed the sub-solar value for β Pictoris b, in combination with its high mass, points towards a formation through core accretion and not gravitational collapse.

We obtained a metallicity of −0.39 ± 0.16 for the nominal model, which is significantly lower than the value of 0.66 ± 0.13 found by GRAVITY Collaboration (2020) on the combined GPI Y, J, and H data with the GRAVITY K-band data, but similar to the value of −0.53 ± 0.3 they obtained using only the GRAVITY K-band data. The derived metallicity is also strongly dependent on the used model, shown by the large deviations in Fig. 4. This is the result of a degeneracy between clouds and metallicity, which is discussed in the next subsection. We note that constraining absolute abundances and thus the planets metallicity using high-resolution spectroscopy is relatively hard, especially when one cannot retain the planets continuum. Combining high-spectral-resolution data with low resolution over a larger wavelength range may help relieve these issues. In our case, this is shown from the retrieval with the priors from

GRAVITY Collaboration (2020), where the metallicity converges to its prior due to the limited information present in the high-resolution spectra.

For the free composition retrieval, we only obtain upper bounds on the minor species. There is a slight preference for a model with ¹³CO, but we consider this marginal and do not discuss it further in this paper.

Fig. 8 Estimation of the radial velocity in the barycentric restframe for individual exposures.

6.2 Clouds and P–T profile

The retrieval for the nominal model results in no constraints on the cloud parameters. This is the result of the retrieved P–T profile not intersecting the condensation curves of the considered cloud species for the retrieved parameters. The retrieved P–T profile is in shape similar to the one obtained by GRAVITY Collaboration (2020) but shifted to slightly higher temperatures. This is because high-resolution spectroscopy, after high-pass filtering, is mostly sensitive to the slope of the P–T profile and not its absolute value. Furthermore, the derived metallicity is much lower, which results in lower condensation temperatures. These effects together lead to cloud-free solution that are inconsistent with the low-resolution data at shorter wavelengths, as shown in Fig. 7.

High-resolution K-band retrieval converging to cloud-free solutions has been seen in other works (e.g., Zhang et al. 2021b; Xuan et al. 2022) and may be the result of the wavelength coverage in combination with the parameterization and priors on the P–T profile. We have tried the P–T parameterization of Mollière et al. (2020) and obtained almost equivalent results to the nominal model. The alternative model with a forced grey cloud deck already agrees better with the low-resolution data. Naturally, the retrieval with the results from GRAVITY Collaboration (2020) as a prior agrees even better with the low- resolution data. Still, the main parameters of interest for this work, C/O, υ sin(i) and RV, seem to be minimally affected by the presence of the clouds, as shown in Fig. 4. The derived metallic- ity is highly correlated with the presence of clouds and changes depending on the modeling choices. This may be caused by the degeneracies between the fitted coefficient of the planet component, clouds, and metallicity, as all these parameters mainly change the depth of the spectral lines. Our observations also do not help constrain the composition and location of the clouds in the atmosphere of β Pic b. Observations over a larger wavelength range with e.g. JWST could resolve this uncertainty, as shown by the recently detected silicate cloud absorption features in VHS 1256–1257 b (Miles et al. 2023).

6.3 Radial velocity and spin rotation

We obtained a radial velocity of 31.9 ±0.3 km s⁻¹ in the barycentric rest frame. Using the systematic velocity, υ_sys = 20 ± 0.7 km s⁻¹ (Gontcharov 2006) for β Pictoris, this gives a radial velocity between the star and the planet of 11.9 ± 0.8 km s⁻¹, where the uncertainty is dominated by the uncertainty on the radial velocity of the star. This is consistent at ≈2σ with the radial velocity, υ_orb = 10.0 ± 0.1 km s⁻¹, predicted from the orbital solution (Wang et al. 2021a). Since we are able to detect the planet in individual frames, we can study the RV evolution over time (shown in Fig. 8), exhibiting no clear variability. We find a rotational velocity υ sin(i) = 19.9 ± 1.0km s⁻¹. This is slightly lower than the value obtained in Snellen et al. (2014), but consistent within 2σ. Using a radius of 1.4 ± 0.1 R_Jup and inclination of 88.95 ± 0.10 degrees, this gives a rotation period of 8.7 ± 0.8 h for β Pic b, assuming that the planet has no obliquity. We note that the limb darkening coefficient is converging to its maximum value in the retrievals. Since there is a correlation between the limb darkening and v sin(i), the rotational velocity could be lower if the limb darkening coefficient turns out to be lower.

7 Summary and outlook

We report on observations of β Pictoris b using the recently upgraded CRIRES+. We present dedicated features in pycrires for the reduction of observations of directly imaged planets with this instrument. By forward-modeling the data in combination with free atmospheric retrievals, we attempted to characterize the atmosphere of the planet. We find a slightly sub-solar C/O ratio, υ sin(i) = 19.9 ± 1.0km s⁻¹, which gives a rotation period of 8.7 ± 0.8 h and a radial velocity of 31.9 ± 0.3 km s⁻¹ in the barycentric restframe. However, the results from the nominal retrieval on solely the CRIRES+ data are inconsistent with the low-resolution data at shorter wavelengths. By either forcing the presence of clouds or using the results from GRAVITY Collaboration (2020) as a prior we obtain solutions that are more consistent with the low-resolution data. Fortunately, the

main parameters of interest are minimally affected by this. This shows the power of high-resolution spectroscopy in breaking degeneracies in atmospheric retrievals and deriving robust elemental abundance ratios. Furthermore, we are able to confidently detect the planet in individual exposures of 2 min, showing the significantly improved capabilities of CRIRES+. Searching for exomoons using radial velocity or monitoring the spectral variability of some of the most favourable objects may thus already be within reach, especially with the arrival of HiRISE (Vigan et al. 2024), which couples the high-contrast imager SPHERE to CRIRES+.

Acknowledgements

R.L., I.S., S.d.R. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program under grant agreement no. 694513. T.S. acknowledges the support from the Netherlands Organisation for Scientific Research (NWO) through grant VI. Veni.202.230. A.S.L. acknowledges financial support from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/ 10.13039/501100011033. This work used the Dutch national e-infrastructure with the support of the SURF Cooperative using grant no. EINF-4556.

Appendix A Visual comparison with Snellen et al. 2014

Fig. A.1 Visual comparison between the cross-correlation signal-to-noise ratio map of β Pic b using the original CRIRES (Snellen et al. 2014) and our new CRIRES+ observations. The relative position along the slit was matched to the convention used in Snellen et al. (2014).

Appendix B Full corner plot for the nominal model

Fig. B.1 Full corner plot of the atmospheric retrieval for the nominal model.

References

All Tables

All Figures

	Fig. 1 Visualization of the data extraction pipeline. Top row shows the product of the obs_nodding recipe from the ESO CRIRES+ data reduction pipeline. The middle row shows the rectified 2D spectrum corrected for the slit tilt and curvature using pycrires and the bottom row shows the spectrum at the location of the planet, plotted together with a telluric model generated with SkyCalc.
In the text

	Fig. 2 Raw white light contrast curve for our CRIRES+ observations. The location of β Pictoris b is indicated in this diagram according to the K-band contrast from Chilcote et al. (2017).
In the text

Fig. 3 Detection of the planet using model templates of the full planet, H2O and CO using cross-correlation (top row) and the likelihood ratio (bottom row). The radial velocity is shown here in the barycentric restframe. In the top panel we also show the cross-correlation value for a reference spectrum, which was taken at the same separation on the other side of the star. In the bottom panels we indicate the location of the star and the planet and the expected velocity of the planet.

In the text

	Fig. 4 Constraints on a subset of the parameters from the atmospheric retrievals using the models discussed in Sect. 4. The RV shown here is in the barycentric restframe.
In the text

	Fig. 5 Retrieved P–T profile for the nominal model with associated uncertainties. Also shown: the P–T profile obtained by GRAVITY Collaboration (2020), the condensation curves of Fe and MgSiO3 for the best fitting parameters, and the emission contribution function for the best fitting parameters.
In the text

Fig. 6 CRIRES+ data from one of the spectral orders along with the fitted model. Top panel shows the total extracted flux at the location of the planet with the full fitted model including star, planet, and systematics. The second panel shows the data after subtracting the stellar and systematic components, along with the best-fitting planet model. The third panel shows the residuals of this fit. Finally, the bottom panel shows the same as the second panel but then smoothed to a spectral resolution of 15 000 to show the broad planet features more clearly.

In the text

Fig. 7 Comparison of the spectral energy distributions of our retrieved models compared to the low-resolution GPI data from Chilcote et al. (2017) and GRAVITY data from GRAVITY Collaboration (2020). Since we do not retrieve radius, the flux of our models is scaled such that it matches the GPI flux in the K-band. The solid line shows the median value and the colored region shows the 16th and 84th percentiles.

In the text

	Fig. 8 Estimation of the radial velocity in the barycentric restframe for individual exposures.
In the text

	Fig. A.1 Visual comparison between the cross-correlation signal-to-noise ratio map of β Pic b using the original CRIRES (Snellen et al. 2014) and our new CRIRES+ observations. The relative position along the slit was matched to the convention used in Snellen et al. (2014).
In the text

	Fig. B.1 Full corner plot of the atmospheric retrieval for the nominal model.
In the text