© The Authors 2025

1 Introduction

With an age of 24±3 Myr (Bell et al. 2015) and a distance of only 19.7 pc (Gaia Collaboration 2016), the A6V Gray et al. (2006) star β Pictoris is one of the best-suited stars in the sky for detailed studies of a young circumstellar environment. Additionally, its circumstellar disk (Aumann 1985) is extraordinarily bright and large, extending out to at least 2000 au (Janson et al. 2021), which led to it becoming the first disk ever outside of the Solar System to be directly imaged (Smith & Terrile 1984). The disk shows an interesting set of sub-structure and features (Kalas & Jewitt 1995), including a warp due to an inclined secondary disk (e.g. Heap et al. 2000; Golimowski et al. 2006); a very large density of comets (Lecavelier des Etangs et al. 2022; Janson et al. 2023) with a range of eccentricities going all the way up to star-grazing orbits (Ferlet et al. 1987; Beust et al. 1990); and the presence of gasses with short expected lifetimes in the strongly irradiated disk, implying recent gas release in large planetesimal disruptions (Brandeker et al. 2016; Cataldi et al. 2018). All these features indicate the presence of planets that dynamically excite various components of the disk. The secondary disk, in particular, has been used for rather precise predictions of a planet orbiting within the disk.

In 2010, the planet β Pic b was indeed confirmed in the system through direct imaging (Lagrange et al. 2010), when common proper motion and orbital motion were confirmed for a candidate that had previously been discovered in a first epoch through reprocessing of archival images (Lagrange et al. 2009). The imaged planet was well consistent with the planet that had been predicted on the basis of the orientation of the secondary disk, which originates from the slight misalignment between the planes of the main disk component and the planetary orbit (Chauvin et al. 2012). Later on, a second planet, β Pic c, was identified first through radial velocity (RV) measurements (Lagrange et al. 2019) and soon after confirmed through direct imaging with the interferometer GRAVITY (Nowak et al. 2020), further demonstrating that β Pic is a rich and dynamic planetary system. The favourable observational properties of the system have led to exquisite orbital determinations of the planets through (primarily) GRAVITY interferometry, to the point that it has been possible to constrain the mass of planet c directly through its impact on the astrometry of planet b (Lacour et al. 2021). The planet masses and semi-major axes are estimated in the ranges of 7–12 M_jup and ~10 au for b, and 7–10 M_jup and ~2.7 au for c (Brandt et al. 2021; Lacour et al. 2021).

Much like in the case of astrometric characterisation, the favourable observational properties of the β Pic system make its planets particularly well suited for spectroscopic characterisation – particularly planet b, due to its larger semi-major axis, leading to a higher angular separation from the star in most epochs. This led to β Pic b being the first exoplanet to be successfully studied with high-resolution spectroscopy¹ (Snellen et al. 2014). The high-resolution (R ~ 50 000–100 000) spectroscopy was acquired in the K band with the CRyogenic high-resolution InfraRed Echelle Spectrograph (CRIRES; Käufl et al. 2004; Dorn et al. 2023) at the Very Large Telescope (VLT) managed by the European Southern Observatories (ESO). The principle behind the technique is based on cross-correlation function (CCF) techniques for effectively summing up the information from whole forests of lines at once, in order to account for the fact that each individual spectral line typically has a very low signal-to-noise ratio (S /N) in the measured spectra. The resulting cross-correlation peak retains many of the features of the individual lines it is based on, and can therefore be used to derive kinematic properties such as the rotational velocity and RV of the planet. Cross-correlation has also been successfully shown to work for medium-resolution (R ~5000) spectroscopy for β Pic b (Hoeijmakers et al. 2018; Kiefer et al. 2024). These techniques are well complemented by spectroscopic studies at a lower resolution that can provide constraints on the absolute flux level as a function of wavelength, which is typically challenging to derive from an analysis based on cross-correlation. Recent results to this end include James Webb Space Telescope (JWST) Mid-InfraRed Imager (MIRI) spectroscopy of β Pic b at 5–7 µm (Worthen et al. 2024), and GRAVITY K band spectroscopy of both planets b and c (GRAVITY Collaboration 2020; Nowak et al. 2020).

Following an upgrade of CRIRES (named CRIRES+) which greatly widens its simultaneous spectral range and enhances its sensitivity, renewed attention has been placed on characterising β Pic b at a high spectral resolution. Recent studies have investigated the planet in the K band (Landman et al. 2024) and in the M band (Parker et al. 2024). Both studies have reported detections of H₂O and CO. The M band study also noted a tentative detection of SiO. In between the K and M bands lies the L band, which as of yet has been less studied spectroscopically. Yet, it is an important wavelength range for the study of young planets, for several reasons. As many direct imaging and spectroscopic studies have noted (e.g. Kasper et al. 2007; Janson et al. 2010; Quanz et al. 2010; Mawet et al. 2013; Janson et al. 2015), the wavelength range around ~4 µm is advantageous in a high-contrast context, due to the favourable contrast between stars and planets for a wide range of planetary temperatures (e.g. Burrows et al. 2006; Linder et al. 2019) in combination with a relatively low thermal background noise in the L band compared to longer wavelengths, as well as the high adaptive optics (AO) stability that can be achieved in the L band. Furthermore, the L band contains strong opacity regions of molecules that are expected to be important in planetary atmospheres, including H₂O, CH₄, and NH₃, and more exotic species such as SiO and SO₂.

In this paper, we present CRIRES+ L band spectroscopy of the β Pic planetary system. While the main target is β Pic b, the CRIRES+ slit alignment automatically also includes β Pic c, and we therefore attempted to extract signal from c as well. The L band range is particularly well suited for studying a contrast-limited planet such as c that would be very challenging to reach at shorter wavelengths, although as we subsequently see, the observational conditions for c were not optimal during these observations. The paper is structured as follows: In Sect. 2, we describe the observations acquired for the purpose of this study. The reduction of the collected data is outlined in Sect. 3, and the subsequent analysis and corresponding results follow in Sect. 4. We discuss the outcomes of the study and consequences for future studies in Sect. 5, and finally briefly summarise the conclusions in Sect. 6.

2 Observations

Observations were carried out with CRIRES+ in Service Mode starting on 29 Nov. 2022. Originally scheduled for ESO observing period P110, only about half of the observations were carried out in that period, the rest being carried over into P111 with the last observation carried out during 1 Nov. 2023. We used four different spectral settings in order to cover all of the L band with minimum gaps and overlaps, and parts of the M band. The observations were divided into 15 observing blocks (OBs), with 3 OBs each dedicated to settings L3262, L3340 and L3426, and 6 OBs dedicated to setting M4318. One incomplete OB executed on 26 Feb. 2023 failed validation and was later re-executed; we have discarded the failed version of the OB. The full set of validated observations is listed in Table 1, where on each of Oct. 12 and 13, two OBs have been merged and are treated as a single observation, since they were acquired consecutively and in the same mode.

In the same table, we also list the expected separation and RV for planets β Pic b and c at each respective epoch, based on orbital parameters in Lacour et al. (2021). Likewise, the barycentric RV of the Earth relative to β Pic is listed for each epoch. Due to the long baseline between first and last epoch, there is significant evolution in these quantities from epoch to epoch, which we take into account when combining the data in Sect. 4.2. The β Pic system RV relative to the barycenter of the Solar System is 20.0 km/s (Gaia Collaboration et al. 2023). For all observations, the position angle of the slit was 31.55 deg, which comfortably encompasses the full orbital planes of planets b and c in the nearly edge-on β Pic system. The slit width was consistently 200 mas, leading to a spectral resolution of ~100 000. Planet b was our main target, but since c is also included at all times, we consider it a secondary target. Unfortunately, the observing epochs were not ideal for c, as we will return to in Sect. 4.4. However, the epochs were (by design) well suited for detection and characterisation of planet b, which resided close to its maximum combined separation and RV relative to the star across the observational baseline.

Observations were acquired using an ABBA nodding pattern in order to enable pairwise subtraction of the high sky level in L band observations. We implemented the metrology procedure offered by ESO at the beginning of each observation in order to ensure reproducibility of the spectral positioning on the detectors between different nights. For all L band settings (L3262, L3340, and L3426), our standard setting was a direct integration time (DIT) of 60 seconds per exposure, with two exposures per nodding position over 15 nodding cycles. For some runs, this was manually adjusted to DIT of 30 s and four exposures per nodding position by ESO staff, in cases where they estimated an elevated risk of saturation. The total on-source integration time per L band observation was therefore 1 hour exactly. For the OBs with the M band setting (M4318), we set DIT to 10 s and the number of integrations per exposure (NDIT) to six with two exposures per nodding position over 13–14 nodding cycles (half with 13 and half with 14) for on-source integration times of 52–56 minutes per OB.

Every spectral setting used for the observations includes six spectral orders, which are each divided over three adjacent detectors in CRIRES+. The data in each epoch can therefore be seen as consisting of 18 spectral regions, which we here refer to simply as ‘chunks’. For all reduction and analysis steps up until the cross-correlation discussed in Sect. 4.2, we operate on each chunk individually. This means that we are able to operate on continuous data sets of ~30 nm width each, in which the point spread function (PSF) is not expected to significantly change from one edge wavelength to the other. For the M band setting M4318, two orders (six chunks) were excluded for each frame: One of the orders overlaps with the fundamental band of telluric CO₂ where the transmission is effectively zero across the whole order, and the second excluded order was the longest-wavelength order in which the thermal background was high enough that the combined star and background flux saturated the detector. The chunks were reduced and analysed as described in the following sections.

3 Data reduction and extraction

3.1 Reduction and calibration

CRIRES+ has a dedicated data reduction pipeline based on the PyReduce software (Piskunov et al. 2021), which under standard circumstances can automatically accomplish the entire data reduction chain from raw data through to fully extracted and calibrated spectra for point source targets. However, in our case, there are two main reasons for why the automatic procedure cannot be directly applied to the data: Firstly, our primary objective is not to extract just the spectrum of the star β Pic A, which dominates the signal in a point-source extraction, but rather to subtract the signal of A and subsequently extract the signals of the fainter b and c companions – this procedure is discussed in Sects. 3.2 and 3.3. Secondly, and of more immediate relevance to the early phases of data reduction, ESO provides no lamp spectra or similar types of calibration files for L or M band observations, which the pipeline needs for its slit tilt, slit curvature, and wavelength calibration steps. We therefore devised a reduction scheme that uses the routines of the standard pipeline whenever possible, with custom procedures mixed in for all steps where this was necessary.

In this way, we produced frames that were corrected for bias, dark, and flatfield effects based on the raw data and standard calibration files. These ‘cleaned’ frames were used for both custom calibration and science extraction purposes, with different processing for each branch. As outlined above, the purpose of the custom calibration in this context was the determination of slit tilt/curvature and wavelength calibration. Both of these calibrations are usually based on lamp spectra, where the full slit is illuminated by emission lines from a known compound. The morphology of the lines on the detector can then be used to determine the orientation and curvature of the slit image as function of wavelength, and the positions of the lines with known emitted wavelengths can be measured to provide the wavelength calibration. In the absence of such lamp spectra, we devised two separate methods for achieving these calibrations based on the science data.

In one method, we used the telluric absorption lines in the β Pic stellar spectra. For each absorption line, the A and B nodding positions provided two data points suitable for determining the slit tilt at that particular wavelength. In principle, the slit image also has a mild curvature in addition to the tilt, which cannot be measured with only two points – however, this curvature is essentially negligible (see Piskunov et al. 2021), so we disregard it in this analysis. A second degree polynomial was then fit to the degree of tilt as function of wavelength, in order to interpolate (or in the case of near-edge pixels, extrapolate) the expected tilt in regions that lacked telluric absorption lines. For wavelength calibration, we then used PySME (Wehrhahn et al. 2023) for identifying the telluric line species and thereby construct a wavelength solution along each spectral trace.

In the second method, we used telluric emission lines instead of absorption lines. For this purpose, we created ‘artificial’ sky frames by dividing every spectral order along the slit direction (near-vertical on the detector) into a lower and upper half. For each A and B nodding frame, we then saved only the half that did not contain the star, put the remaining halves together to create sky frames where both the upper and lower parts were devoid of starlight, and averaged the results into deep frames containing only (to first order) light from the sky. Due to the high image depth of the observations, this provided sky frames with a reasonably rich set of lines in most orders. Slit tilts and wavelength calibrations were then determined in an analogous way as in the absorption-based method. The two methods complemented each other well, in that some regions that were relatively poor in absorption lines still contained a decent number of emission lines, and vice versa. In orders where both methods produced acceptable results, those results were mutually consistent, with the absorption method generally producing higher-quality spectral lines and thus being the preferred choice in such cases.

With this procedure, we were able to produce slit geometry and wavelength calibration, in addition to the standard calibration produced in the default manner. The standard pipeline procedure could then be used to produce a fully reduced and extracted spectrum of the star β Pic, if interpreted as a single point-source. This is useful for normalisation purposes, as we will see in Sect. 3.2. However, the main scientific objectives of the project require us to perform PSF subtraction on the star β Pic A and subsequently extract the signals of the much fainter planets β Pic b and c. That procedure will be discussed in the following section. The customised procedure for determining slit geometry and wavelength calibration described above, for the spectral settings we used (L3262, L3340, L3426, and M4318), is available on github².

3.2 PSF determination and subtraction

The PSF subtraction scheme performed here is a variant of the STARK routine previously implemented on VLT/UVES (Ringqvist et al. 2023) and JWST/NIRCAM (Patel et al. 2024) data, which is available on github³. For the PSF subtraction purposes, we operate on the A – B and B – A pairwise background- subtracted images produced in the procedure described in the previous section. The input pixel values are fluxes in units of detector counts, and each value is also associated with a corresponding error estimation based on the calculated photon count (from non-background subtracted frames) and read noise in each pixel. The main features of the procedure are described below.

The CRIRES+ spectra are mapped on the detector(s) in such a way that the dispersion (wavelength) direction is primarily horizontal and the spatial direction along the slit is primarily vertical. However, as discussed in Sect. 3.1, the slit image is in fact not completely vertical, but has a slight tilt (see Fig. 1), with a tilt angle that changes slowly across the wavelength axis. In other words, pixels along a detector column (within an order) correspond to similar, but not equal, wavelengths. Thus, as a first step, we consider each pixel within a spectral chunk with pixel indices x_i for the row number and y_i for the column number on the detector. Based on the wavelength and slit tilt calibrations described previously, we calculate the pixel’s wavelength λ_i and its spatial separation along the slit from the stellar trace, ρ_i.

At this point, all relevant pixels are assigned wavelengths and distances from a central coordinate, which as a collective make up a non-uniformly sampled grid. An intuitively natural next step might then be to interpolate a new spectral image, which is sampled uniformly in wavelength-separation space. However, operating with interpolated data at this stage of the process can negatively impact the contrast performance, since imperfections in the interpolation can lead to large errors for any faint planets residing in the bright stellar PSF wings. Instead, we represent the PSF by using a χ² minimisation fit to the data with a model in the form of a 2D spline in wavelength-separation space. A spline does not need to operate on a uniform grid, and can therefore be applied directly to the non-interpolated λ_i and ρ_i coordinates. This procedure also makes it easier to account for hot/bad pixels, which can be identified through sigma clipping and simply ignored in the fitting. To simplify the fitting as much as possible, we renormalise the data by wavelength prior to the fit using the previously extracted stellar 1d-spectrum. This straightens out both tellurics and stellar features, creating a function that always peaks at a value of 1 at the centre of the PSF for all wavelengths. Importantly, when re-scaling the fluxes, we also re-scale their corresponding error estimates by the same factor. The error estimates are included in the χ² fitting, and mean (for example) that data points inside of deep telluric lines acquire a very low weight in the fitting, since the relative error is very large inside such features.

For the spline fitting, we use a uniform grid of knots with ~600 knot points along the spatial (ρ_i) direction, and only 3 knot points along the spectral (λ_i) direction. The specific number of spatial samples is arbitrary, but is set to be high in order to provide an accurate PSF fit along the spatial axis, whilst still maintaining the total number of knot points much smaller than the total number of available data points (~120 000) per chunk. Meanwhile, the reason for the modest spectral sampling is to allow the PSF evolve, but only very smoothly, along the spectral axis. This is sufficient to characterise chromatic PSF effects (which should be small and smooth, given the short wavelength span of the spectral chunks), but simultaneously, it minimises damaging impact on the planetary spectral signal. The continuum flux of any planet around the star inevitably gets included in the PSF fit to the star, and is therefore subtracted out in the PSF subtraction; however, planetary spectral lines are much narrower than the spacing between the spectral spline knots, and thus impossible for the spline to fit. Such lines are therefore efficiently conserved when the fit is subtracted. This principle is confirmed by our injection testing in Sect. 3.4.

In this way, it is possible to produce a 2D spectral PSF model, evaluate it at the exact λ_i and ρ_i points sampled on the detector, and subtract the model from the data, without any need to interpolate the data itself. The result is a residual frame that contains only noise at the location of the star, but retains signal (in the form of spectral lines) from all planets around the star, although such signals are too weak in each individual pixel to be distinguishable by eye. As can be seen in Fig. 2, the procedure works well except in particular spectral regions, which correspond exactly to where there are strong telluric lines in the spectra. Effectively, the PSF as mapped on the detector appears to have a slightly different shape inside and outside of the tel- lurics. The reason for this is either a slight non-linearity in the detector readout, or a consequence of the fact that the PSF is slightly narrower than the slit, leading to a slightly non-uniform slit illumination (see Landman et al. 2024). There are ways to account for these effects that we expect to implement in later versions of the procedure, but for the purpose of this study, the telluric regions are in any case not usable, due to the low transmission and consequential high relative noise in those regions. Instead, we simply let the telluric regions get masked out during the sigma clipping in Sect. 3.3.

Fig. 1 Flux-normalised spectral traces of the star β Pic, for one single spectral chunk. The panels are not plotted as pixeled images (which would not be possible for the bottom panel, since it is not uniformly sampled), but as one coloured dot per sampled data point. Yellow is a normalised flux of 1, blue is 0. Top panel: flux sampled by the xi and yi coordinates on the detector. Bottom panel: the same flux (no interpolation) sampled by wavelength λi and separation ρi.

Fig. 2 Example of residual data after PSF subtraction. The figure is in detector xy-space, since it can then be plotted as a non-interpolated pixeled image. Spatially (vertically), the image is centred on the stellar trace, which after PSF subtraction only remains as a narrow horizontal region of increased noise. Spectrally (horizontally), the image is centred on a strong telluric absorption feature. In this particular region, the PSF subtraction does not perform as well as along the rest of the spectrum, resulting in a vertical net positive feature extending out from the centre. At even larger separations vertically from the centre, the noise in the telluric region is generally enhanced. This feature would remain even if the PSF subtraction had been perfect, since it represents the much lower transmission inside of the telluric line than in the continuum. Planet b is ~9.5 pixels below the stellar trace, but cannot be seen without crosscorrelation techniques.

3.3 Extraction of planetary spectra

From the residual frame after stellar PSF subtraction, we can extract spectra from the circumstellar region in different ways. A straightforward example is to extract the spectra of planets b and c based on a priori knowledge about their locations relative to the star from previous astrometry (Lacour et al. 2021). For doing so, we define synthetic ‘sub-slits’ centred on the known exoplanet location, with a slit length of 5 pixels along the spatial direction (i.e. 2.5 pixels in each direction from the planet centre). The extraction is done epoch by epoch, so the slit location varies with the planetary motion over time. This is not particularly important for planet b, since its motion is not very large (~0.5 pixels) over the observational baseline, but considerably more important for planet c with its much more rapid motion.

The sub-slit now represents an aperture that contains the vast majority of the planetary flux. From this, we extract a 1D spectrum with a wavelength sampling corresponding to one pixel on the detector. Since the λ_i and ρ_i coordinate grid is non-uniform relative to the pixel grid, there is not an equal number of data points in every spectral bin, and the data points are not equally distributed along the planetary PSF, so a straight aperture sum of the flux would not work satisfactorily. Instead, we implement a PSF-normalised weighted aperture sum to account for these effects. For this purpose, we use a shifted version of the already determined stellar PSF from the PSF subtraction procedure, to now represent the planetary PSF. Every data point inside the subslit is normalised by 1/f (ρ_p) where f (ρ_p) is the flux of the PSF (relative to the peak flux) at separation ρ_p from the planetary PSF centre. We normalise the estimated errors by the same factor. An error-weighted mean is then taken inside the spatial aperture for each spectral window to represent the extracted 1D spectral data points.

As references for the planetary spectra and for evaluating any possible systematic noise sources affecting them, we also extract ‘mirror’ spectra from the same separations as the planets, but on the exact opposite side of the star. This yields data sets with similar noise properties as the planetary spectra, but without the planetary signals. Furthermore, we perform a series of extractions in the same way as above, but consecutively centred on integer pixel separations from −15 to 15 pixels from the star. This creates a spectral map which can be used for purposes such as blind searches, or checking systematic noise levels. We subject all extracted spectra to an iterative sigma clipping procedure in which σ is determined from the scatter in the data, after which >4σ deviating data points are removed, followed by a new estimation of σ until convergence is reached. This efficiently removes heavily telluric-affected regions, as well as remaining bad pixels. Due to the normalisations earlier in the procedure, the spectra are conveniently dimensionless, and correspond to the contrast between the (real or hypothetical) planetary flux and the stellar flux at each sampled wavelength. The spectral coverage and the scatter in the contrast at different wavelengths is illustrated in Fig. 3.

3.4 Injection of artificial planet signals

For the purpose of the sensitivity evaluation analysis that will be discussed in Sect. 4.3, we make tests where we inject artificial planet spectra into the data before PSF subtraction and extract them after PSF subtraction. For example, if we wish to evaluate detectability of CO₂ in planet b, we first generate a theoretical model spectrum with the expected physical properties of planet b (see Sect. 4.1 for a description of how the model spectra are generated), and a CO₂ atmosphere. We transform the flux spectrum into a contrast spectrum by dividing it with a stellar model spectrum for β Pic A. In the PSF determination and subtraction scheme (Sect. 3.2), after the first estimation of the stellar PSF, we then produce an artificial planet by scaling the stellar PSF with the planetary contrast spectrum, and adding the resulting signature into the data at the empty mirror location for β Pic b. We then re-estimate the stellar PSF based on the new data, subtract it, and then extract the signal at the location where the artificial planet was introduced. This results in spectra with noise that is well representative of the noise affecting the real planet, and with a controlled artificial CO₂ abundance for testing whether that abundance should be detectable or not. The process can then be repeated for other molecules for planet b, and in principle also planet c, although such an analysis for c is more complex, as discussed in Sect. 4.4.

The above procedure is useful, but quite time consuming, since every new molecule and every new parameter for each molecule requires a new full PSF subtraction procedure to generate artificial contrast spectra. We therefore also implement a ’pseudo-injection’ scheme, where instead of adding an artificial planet to the 2D spectra prior to PSF subtraction, we calculate a model contrast spectrum for a given molecule and a given set of parameters for that molecule (e.g. temperature, volume mixing ratio, see Sect. 4.1). This model contrast is then added to the extracted spectrum on the opposite side of the star from the planet, in order to include a representative noise level. A normal cross-correlation procedure as discussed above is then performed in order to evaluate the detection strength for the (pseudo)- injected signal. The advantage of this method over actual injection is that it is much faster, and can therefore be applied to a wider range of molecules and molecular parameters. The disadvantage is that it fails to include any self-subtraction effects that take place on real signals during the PSF subtraction step. However, as we will see, these effects are quite minor for the existing data sets, so after running a real injection on water and comparing to a pseudo-injection (see Sect. 4.3) with the same parameters to verify sufficient accuracy, we used the pseudo-injection scheme for most of the injection based analysis of this study.

Fig. 3 Spectral coverage of the observations in this study. Four spectral settings were used: L3262 (upper left), L3340 (upper right), L3426 (lower left) and M4318 (lower right). The coloured regions show the coverage of that particular setting, while the grey regions show the coverage of the other settings. The scatter corresponds to the noise level in each spectral region. Higher scatter can be caused by shorter integration time, worse ambient conditions, or lower instrumental transmission, but here it primarily reflects wavelength ranges of higher thermal background, and/or wavelength ranges of lower atmospheric transmission.

4 Analysis and results

The extracted spectra are analysed using cross-correlation techniques, which are frequently employed in high-resolution exoplanet spectroscopy (e.g. Snellen et al. 2014; Schwarz et al. 2016; Bryan et al. 2018, 2020; Xuan et al. 2020). Here we describe the implementation of the technique, and present the results it has yielded.

4.1 Spectral models

For spectral modelling, we use petitRADTRANS (Mollière et al. 2019, 2020) models based on the relatively well constrained physical properties of planets b and c, with an atmosphere consisting of hydrogen and helium (which contribute atmospheric pressure, but no strong lines in the wavelength region we consider here) plus any individual molecule we intend to search for in the extracted spectra.

Following (Guillot 2010) and the previous recent CRIRES+ spectroscopic studies of β Pic, we use a temperature-pressure profile assuming negligible irradiation contribution, as follows: $$T^4=\left(\frac{1}{2}+\frac{3P{\kappa }_{\text{IR}}}{4g}\right)T_{\text{int}}^4$$(1)

where T is the temperature profile as function of pressure P, κ_IR is the infrared absorption coefficient which we set to κ_IR = 0.005 cm²g⁻¹, 𝑔 the surface gravity of the planet, and T_int is the internal temperature.

For planet b, we use an internal temperature of 1700 K as in Parker et al. (2024), and log 𝑔 = 4.18. For planet c, we set the internal temperature to 1250 K based on Nowak et al. (2020), and log 𝑔 = 4.24. We produce spectra both with and without rotational broadening. The rotational velocity is 20.36 km/s (see Sect. 4.3) for planet b. For planet c, the corresponding velocity is unknown so we set it to a similar value as for b, 20 km/s. For molecule identification in the extracted spectra, we typically use broadened spectra unless stated otherwise. The results based on non-broadened spectra are very similar. We also fit and subtract the continuum from those spectra, in order to distill only the line information. For the synthetic planet injection described in Sect. 3.4, we use broadened spectra. We also include the instrumental broadening, although it is practically negligible relative to the rotation − at R ~ 100 000, the instrumental full width at half maximum (FWHM) is approximately 3 km/s.

4.2 Cross-correlation method

Through cross-correlation of the extracted spectra with the models from Sect. 4.1, we can search for specific molecules by noting which modelled species produce a significant cross-correlation peak; and from the width and velocity shift of any such peaks, we can also determine the rotational velocity and RV of the planet. However, the spectral extraction procedure produces many different spectral chunks from different parts of different spectral orders at different epochs, where each individual chunk is relatively short (~30 nm) and therefore does not necessarily contain sufficient abundances of lines by itself to yield a significant cross-correlation peak. The true wealth of this data set stems from the fact that it contains 204 different spectral chunks (3 chunks per order, 4–6 orders per spectral setting, 4 spectral settings with 2–5 epochs per setting). For the cross-correlations, we combine the various triplets of chunks corresponding to the same order into continuous spectral regions, but that still leaves 68 spectral regions that need to be combined in some way in order to maximise the detection significance and quality of the results. The fundamental procedure to accomplish this is to cross-correlate each spectral chunk with a corresponding model, and average the resulting cross-correlation functions to enhance the signal-to-noise ratio S/N.

However, taking a straight average over all spectral regions would not be a good path towards this aim, because for a given molecular species, some regions may contain many intrinsically strong lines, while other regions may contain essentially nothing. Mixing empty regions in with the rich regions would then only add noise and no signal, decreasing the total S/N. Instead we can, at least in principle, achieve a near-optimal S/N by appropriate weighting during averaging of the different regions. We implement this through a weighted average scheme in which the weights are determined by the expected information content in each region. As proxy for the information content, we use the average spectral gradient in the data. That is to say for each spectral region, we take a derivative of the modelled spectrum in that region, and calculate the mean of its absolute values. This accounts for the variations in the spectrum imposed by the spectral lines: If there are many strong spectral lines within the region, it acquires a correspondingly high weight, while if the there are no lines, the weight is zero and the region is effectively excluded.

We have additionally tried other methods of weighting and combination of weightings, including weighting for the expected signal-to-noise ratio of each spectral region, and using only binary weighting that either includes a region if it contains spectral lines or excludes the region if it is devoid of lines. However, we find that the different weighting methods make very little difference for the outcomes, hence it appears that the specific weighting method is of little importance – the dominant factor for providing a strong S/N in the cross-correlation output for our current data set is simply the inclusion of data that contains lines, and exclusion of data that contains no lines.

Fig. 4 Cross-correlation S/N map for H2O, showing a clear detection of β Pic b at its expected separation and velocity. The velocity is in the frame of reference of the planet itself. The separation changes between 536 mas and 558 mas during the observational baseline between 29 Nov. 2022 and 1 Nov. 2023.

4.3 Planet b

The by far most prominent spectral feature of planet b, as observed in the L band, is H₂O absorption. In Fig. 4, we show a cross-correlation map of the inner β Pic system, in the velocity frame of β Pic b. While the planet moves in separation by 22 mas (approximately 0.5 pixels) between the first and last epoch of observations, causing some mild vertical smearing in the image, planet b stands out clearly at its expected position with a high level of statistical significance, S/N = 15.0. The S/N is evaluated as the cross-correlation peak value divided by the standard deviation across the cc function at the planetary separation, excluding a range of ±70 km/s centred on the planetary velocity. This is the same definition for the S/N metric as used in Parker et al. (2024), and slightly different from the definition in Landman et al. (2024), which evaluates the noise on the basis of the standard deviation of the cross-correlation function at the symmetrically opposite location of the star, relative to the planet. We assess that the stellar PSF in our data set is sufficiently asymmetric that the noise properties at equal separations in opposite directions from the star are not necessarily fully mutually representative, hence we opt for the former approach. The CCF is quite clean from systematic noise in the continuum around the peak, relative to much of the CCF analysis in the literature. This is most likely a consequence of the fact that the LM band data consists of a wide range of different spectral regions, which help in mitigating such systematic noise, as we discuss in Sect. 5.

In order to evaluate the rotational velocity of the planet, we perform a Markov Chain Monte Carlo (MCMC) procedure over a two-dimensional grid of parameters representing a broadening kernel that we apply to the H₂O model and compare to the observations. The two parameters are the rotational velocity υ_rot and the linear limb darkening parameter c, where c = 0 indicates no limb darkening, and c = 1 indicates very strong limb darkening. This is essentially identical to the corresponding procedure in Landman et al. (2024), except that we use a slightly different broadening procedure from (Carvalho & Johns-Krull 2023). This provides a computation-efficient way of accounting for the variation in broadening as function of wavelength, which is particularly important in our case due to the wide (2800– 4800 nm) wavelength range that we cover with the observations. The results are shown in Fig. 5. We find a best-fit solution of υ_rot= 20.36 ± 0.31 km/s and $c={0.92}_{-0.06}^{+0.07}$. These are very high precision constraints for β Pic b, consistent with the high significance and low systematic noise of the wate feature as discussed above; and fully consistent within error bars with the derived rotational velocities of 19.9±1.0 km/s from Landman et al. (2024) and 22±2 km/s from Parker et al. (2024). Interestingly, our derived limb darkening parameter, $c={0.92}_{-0.06}^{+0.07}$, is very high and pushing on the upper end of the parameter range, just like in Landman et al. (2024) (0.8±0.2). The limb darkening parameter is included in the fitting primarily to avoid potential systematic errors on the rotational velocity estimation that could arise from omitting it. Hence, we do not draw any firm conclusions on the limb darkening itself here, except to note that the values are arguably surprisingly high, given the long wavelengths of both the K band and L band observations, and the fact that limb darkening is generally expected to decrease with increasing wavelength. For reference, Parker et al. (2024) do not fit for a value for the limb darkening parameter, but instead fix it at c = 0.3 for their M band observations, based on models of brown dwarf atmospheres. If we strictly enforce c = 0.3 in the same way, we acquire a lower rotational velocity of 18.7 km/s. This is still marginally consistent with the Landman et al. (2024) and Parker et al. (2024) values at the ∼2σ deviance level, but is a comparatively poor fit to our data, as can be read out from the limb darkening posterior resulting from a uniform prior on c as in Fig. 5.

In addition to H₂O, the only other molecule that has so far been securely detected in β Pic b is CO (Snellen et al. 2014). Unlike for H₂O, which has molecular transitions all across the L band making it a strong region for such detections, CO has essentially no lines at all across the entire L band range (see Fig. 6). In our data, only a single spectral window in our M band setting covers any CO transitions. We therefore expect the data to be unsuitable for finding CO, and indeed, while our cross-correlation with a CO model does show a peak at the velocity of planet b, it is not at a formally significant level, and we therefore cannot put strong additional constraints on CO beyond what has been done in previous studies in wavelength regions with higher CO opacities (e.g. Landman et al. 2024; Parker et al. 2024). However, the presence of CO lines in the M band window covered by our observation is supported by the fact that crosscorrelation with a model that includes both H₂O and CO slightly increases the S/N relative to a cross-correlation with a model containing only H₂O, going from 15.0 without CO to 15.4 with CO included.

A more promising molecule in terms of opacities in the L band range is SiO. A tentative detection of SiO in the atmosphere of β Pic b was reported in Parker et al. (2024), based on cross-correlation with a rotationally broadened SiO model spectrum, giving S/N = 4.3, although the tentative nature of the detection was emphasised by the lack of any clear signal in cross-correlations with a non-broadened model spectrum. The spectral setting used in Parker et al. (2024) had SiO lines only in a single spectral window. In our observations, SiO lines of similar strength are covered by 2–3 different spectral windows, with a similar integration time depth per window. Hence, if the S/N = 4.3 feature in Parker et al. (2024) corresponds to a real SiO detection, we should expect to observe such a signature at a S/N in the range of ~6–7 in our data. We therefore performed cross-correlations with SiO models, where we used both models produced by our standard procedure described in Sect. 4.1, as well as the exact model used in Parker et al. (2024), kindly provided by L. Parker, with and without rotational broadening. Neither cross-correlation however shows any SiO detection. While there is a bump in the rotationally broadened cross-correlation (see Fig. D.1), it is not quite at the right velocity, and it is not statistically significant, with S/N of only 2.2. Given that the previous potential detection was cautioned as being only tentative in Parker et al. (2024), the non-detection at significant levels in this study implies that SiO does not occur at (yet) detectable levels in the photosphere of planet b.

Many other molecules that have been discussed in the context of planetary atmospheres have good line regions in the L band spectral range, including CH₄, CO₂, NH₃, HCN, C₂H₂, H₂S, O₃, and PH₃ . We have attempted cross-correlations for all of them, but with no significant detections (see Appendix). In order to assess to which extent each molecule would have been detectable if it were present, we have performed (pseudo-)injections for each molecule as described in Sect. 3.4. Fig. 7 demonstrates that the fast pseudo-injection procedure produces results that are equally well representative of expected signal recovery as the slower full injection scheme. By observing at which volume mixing ratios (VMRs) a given molecule is recoverable, we can exclude the presence of each molecule at those VMR levels, while molecular presence at lower levels than the recoverable limit are still possible. In this context, it is important to note that the detectability of a given molecule only increases with increasing VMR up to a point. Once the VMR becomes high enough that a substantial number of individual lines start saturating and blending, the detectability starts decreasing as VMR increases further. Hence, while very high concentrations (say, 1% VMR) of relatively exotic molecules can probably be seen as physically dubious, they are typically not formally excluded by the spectroscopy. The levels at which each compound is detectable (or not) are summarised in Table 2. Here we use 3σ boundaries since, while a 3σ level would not be sufficient for a definitive detection, it would be sufficient for a tentative detection worthy of further study, similar to the SiO tentative feature from Parker et al. (2024).

The injection procedure demonstrates that the β Pic b spectra would have enabled detection of a wide range of molecules, had they been present in suitable amounts. The sensitivity to CO and CO₂ is however insufficient for detection, which is as expected, given that both molecules have few lines in the spectral windows covered by the observations. An intuitively more surprising feature is the low sensitivity to CH₄, which has lines across the comparatively wide ∼3.0–3.6 µm. The reason for the relatively poor sensitivity to CH₄ is a combination of the facts that its lines are still relatively few and shallow compared to other tested molecules; that it is affected by tellurics to a greater extent than many other molecules; and that the CH₄ lines reside at relatively short wavelengths within the covered span, where the planetary continuum is fainter, which decreases the achievable S/N.

Fig. 5 Results from MCMC fitting of the spin velocity and limb darkening parameter. Upper left: histogram for the spin velocity. Upper right: evolution of the spin velocity (blue) and limb darkening parameter (gold) along the MCMC random walk. Lower left: correlation between spin velocity and limb darkening. Lower right: histogram for the limb darkening parameter.

Fig. 6 Model spectra after continuum subtraction, showing where strong lines occur for four different molecules: H2O in blue, CO in red, CH4 in green, and SiO in orange. The grey regions are the wavelength ranges covered by our observations.

Fig. 7 Comparison of the signal injection scheme compared to the faster ‘pseudo-injection’ scheme described in Sect. 3.4. Top: direct comparison of the injected spectra using regular injection (black line) and pseudo-injection (red line). Also shown as grey dots are extracted spectral points at the location of the injection. Bottom: CCF outputs using injection (black line) and pseudo-injection (red line), compared to the actual water CCF for planet b (dashed line). The injection and pseudo-injection are both well representative of the real signal.

4.4 Planet c

A particularly beneficial aspect of the L band for planet spectroscopy purposes is the relatively low intrinsic contrast between star and planet, for many types of planets. In K band, β Pic c would be very much harder to detect with CRIRES+ than β Pic b, primarily because the noise from the stellar PSF is much higher for c than for b due to the smaller angular separation. This is true also in the L band, but the difference is much smaller, since the balance between background noise versus PSF noise leans much more in the direction of background dominance in the L band. The difference in noise level between the positions of planet c (at maximum separation) and planet b in the CRIRES+ data is only a factor ~3 at the long-wavelength end and a factor ~7 at the shortwavelength end. Hence, we considered it worthwhile to attempt detection of molecular features in the atmosphere of planet c, since it was anyway in the slit that was aligned for including planet b, thanks to the edge-on orientation of the system.

Unfortunately, the timing of the observations was far from ideal for studying planet c. In Nov.–Dec. 2022, its projected separation from the star was <50 mas (and therefore <1 pixel), which is not a workable configuration. In Jan.–Feb. 2023, the separation was 50–80 mas, which is also not ideal. In Sep.–Nov., on the other hand, the separation was essentially perfect at close to 150 mas – however, at this particular orbital configuration (south-western quadrature), it happens to be the case that the systemic RV of the β Pic system and the orbital RV of planet c around its star cancel out almost exactly. The shift between the planet c atmospheric lines and the corresponding telluric lines is then determined almost entirely by the projected barycentric velocity of Earth around the Sun, which is small in the case of β Pic (6–7 km/s in late Autumn 2023) since it is located far from the ecliptic plane. For future efforts, the ideal ephemeris for ground-based high-resolution spectroscopy of β Pic c is when it resides near its north-eastern quadrature, where it is simultaneously at a large angular separation from its star, and at a high velocity shift relative to telluric lines.

For the cross-correlation of a water model with the extracted spectrum of planet c, we use only the epochs from Jan. 2023 onwards when the planet resides at a reasonable separation from its parent star. The result is shown in Fig. 8. Interestingly, a peak is seen at a significance level of S/N = 3.14. The peak is an intriguing feature, although several caveats should be kept in mind for its interpretation. Firstly, the significance level is well below S/N = 5, which we regard as a threshold required for a definitive detection. Secondly, it is not clear that planet c could yield a significantly strong H₂O feature to reach S/N = 3.14 in the data. Injection tests cannot be done in the same way as for planet b, or at least not to the same level of rigor, because for planet c we have much fewer constraints on its absolute brightness level. If planet c had a similar brightness level as planet b, S/N = 3.14 might be quite reasonable. This would be consistent with the fact that the planets have similar masses (Lacour et al. 2021), and presumably, near-equal ages. However, GRAVITY spectra (Nowak et al. 2020) imply that planet c is substantially colder than planet b, as well as smaller. In that context, the predicted line strengths are far weaker, and no significant or near-significant detection would be expected.

For these reasons, we emphasise that the peak in the water cross-correlation for planet c should only be regarded as a tentative feature, requiring further data to confirm its status. We have also tried cross-correlations with the other molecules listed in Sect. 4.3, but without any detections – which is entirely as expected, given the lower signal-to-noise ratios relative to the case of planet b.

Fig. 8 Cross-correlation functions of β Pic b (left) and c (right). The black lines are the cross-correlation functions at the locations of the planets, while the grey lines are the corresponding functions at the exact opposite side of the star relative to the planets. Note the different scaling on the y-axes between the two panels. Planet b is clearly detected, while planet c shows a bump that can only be seen as a tentative feature.

5 Discussion

Thanks to the depth and spectral coverage of the L and M band data in this study, and the strength of the water feature in planet b, it is possible to draw conclusions regarding which observational parameters are the most important for maximising sensitivity in future studies. We illustrate this in Fig. 9. The H₂O line is detected at a significance level of around S/N = 5 in every individual setting at every epoch. When combining different epochs for the same setting, the S/N barely increases at all. By stark contrast, when different spectral settings are combined, the S/N increases rapidly. This implies that, despite the observations being largely thermal background limited for planet b, noise sources that are systematic in the time domain (but variable in the wavelength domain) dominate the cross-correlation S/N. The effect is closely connected to the spectrally non-white oscillations along the baseline of the cross-correlations, which decrease relative to the peak when different spectral settings are combined. We recall that previous spectroscopic analyses based on single spectroscopic settings often show large oscillations of this nature (e.g. Landman et al. 2024), supporting this interpretation.

Residual features of telluric absorption could cause oscillations in the cross-correlation functions, although since the different epochs are taken over a wide time span, the telluric features move around in velocity space relative to the planet. They should therefore average down when combining different epochs of the same setting, which is not observed. Probable contributing effects to the oscillations are inaccuracies in the water models (for example, the line lists) or regularities in the spacings between individual water lines, causing partial degeneracies in the cross-correlation output. Whatever the underlying cause, these results have implications for cross-correlation studies of planetary atmospheres in general, implying that an optimal observational approach for maximising S/N is to spread the wavelength coverage as much as practically possible. In other words, if several spectral windows are available in which a planet is expected to be observable and have interesting molecular lines, then multiple separate shallow observations covering all of those spectral windows is likely to be a more fruitful strategy than spending the corresponding observing time on deep observations of a single window. These results also imply that in a future study, there could be potentially much to gain from combining all of the existing high-resolution data sets of β Pic b, covering K band, L band and M band, for maximising S/N and potentially detect a wider range of molecules in the atmosphere.

The main advantages of the L band range for high- resolution spectroscopy, aside from sensitivity to a wide range of molecules, are the favourable planet-to-star contrast for a wide range of planetary temperatures, and the relative ease to acquire high adaptive optics performance and PSF stability. The main disadvantage, meanwhile, is the high level of thermal background. This limits the utility of 8m-class telescopes and their instruments (such as VLT/CRIRES+) for L band exoplanet spectroscopy. VLT/CRIRES+ is excellently suited for studying β Pic since the β Pic planets are the brightest planets known, thanks to their proximity, youth, and high masses. Fainter planets will increasingly become dominated by the thermal noise, generally necessitating impractically long integration times with CRIRES+. Instead, the major application of L band high-resolution spectroscopy will come with the next generation of large-aperture (20–40 m diameter) telescopes. Such telescopes feature an enormously improved background-limited performance, due to the strong dependence of telescope sensitivity on aperture size. The largest, and therefore least background limited, of such telescopes will be ESO’s 39 m Extremely Large Telescope (ELT). One of the first-generation instruments on ELT will be METIS (Brandl et al. 2021), a mid-infrared imager and spectrograph with a spectral resolution in the L and M bands of up to ∼100 000. Thanks to the beneficial properties of the L band range once the background limit is improved, METIS will be an ideal instrument for detailed characterisation of a large number of directly imaged planets. Similar instruments with L band high-resolution capability are foreseen also for the other next-generation facilities: The Giant Magellan Telescope Near-Infrared Spectrograph (GMTNIRS; Jaffe et al. 2016) on the 20 m Giant Magellan Telescope, and the Mid-IR Camera, High-disperser & IFU spectrograph (MICHI; Packham et al. 2018) on the 30 m Thirty Meter Telescope (TMT).

Fig. 9 Cross-correlation results (normalised to units of S/N) for different ways of combining the individual data sets. In each case, the black dashed line represents the result when combining all the data together. Top left: orange lines are the CCFs from each individual observation taken in the L3262 setting. The red line is the combination of all L3262 observations. No dramatic improvement in S/N is acquired from this combination. Top right: same, but for the case of L3340. Middle left: same, for the case of L3426. Middle right: same, for the case of M4318. Bottom: here, a set of single observations from each spectral setting is plotted as orange lines. Their combination is plotted as a red line. In this case, a significant improvement is gained from combining the data sets, almost reaching the S/N level from the full data set.

6 Conclusions

In this study, we have analysed the results of CRIRES+ LM band spectroscopy of the β Pic system, primarily for the purpose of characterising β Pic b, and secondarily for searching for detectable features in the atmosphere of β Pic c. We outline our main conclusions below:

1. We detected H₂O in the atmosphere of β Pic b at a significance level of S/N = 15.0. If CO is included in the model used for the cross-correlation, the S/N increases to 15.4. The rotational velocity was measured with a high precision as υ_rot = 20.36 ± 0.31 km/s. These results are well consistent with previous studies (Landman et al. 2024; Parker et al. 2024);

2. We do not confirm a previous tentative indication of SiO (Parker et al. 2024) in the atmosphere of β Pic b. A mild peak appears in the CCF, but it is not at the expected velocity, and is not statistically significant at S/N = 2.2. Based on the theoretical SiO spectrum, our sensitivity should be the highest acquired so far. Still, injection results imply that detection of actual SiO should be challenging even in our data set;

3. We detected a peak in the CCF for H₂O at the expected separation and velocity of β Pic c. However, it is not a statistically significant detection (S/N = 3.14). It is also unclear if planet c is intrinsically bright enough to cause the feature. We therefore consider the feature as tentative. Further studies would be necessary in order to establish whether or not the feature is real;

4. Combining data sets with different spectral settings rapidly increases S/N for the H₂O feature of planet b, while combining data sets with equal spectral settings results only in small improvements at best. This implies that the widest possible spectral coverage will be a more important factor than the deepest possible integrations when planning CCF-based exoplanet studies.

Altogether, the observations demonstrate the power of high- resolution spectroscopy in the LM band range. This wavelength range is currently limited for exoplanet purposes by high thermal background, but will greatly expand in applicability with the advent of next-generation telescopes and instrumentation, including ELT/METIS (Brandl et al. 2021), GMT/GMTNIRS (Jaffe et al. 2016), and TMT/MICHI (Packham et al. 2018).

Acknowledgements

M.J. gratefully acknowledges funding from the Knut and Alice Wallenberg Foundation and the Swedish Research Council. The authors thank Thomas Marquart, Luke Parker, Rico Landman and Jayne Birkby for useful discussions regarding the analysis of CRIRES+ data. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The study has made use of the CDS, NASA ADS, and NASA exoplanet archive services, and python packages matplotlib (Hunter 2007) and numpy (Harris et al. 2020).

Appendix A Additional molecular lines

Fig. A.1 shows the expected lines in the L band region (and its surroundings) for 7 molecules that have been included in the CCF analysis but not discussed at depth in the main text.

Fig. A.1 Theoretical spectra after continuum subtraction of molecules that we have searched for in the cross-correlation analysis. A similar figure for H2O, CO, CH4 and SiO is shown in the main text, Fig. 6. From top to bottom, the additional molecules shown here are CO2 in cyan, NH3 in magenta, HCN in light blue, C2H2 in dark red, H2S in pink, O3 in purple, and PH3 in orange. Grey areas correspond to the wavelength windows covered by our observations.

Appendix B H₂O injection

Fig. B.1 shows the results of the injection analysis for H₂O. The injection procedure predicts strong detectability of H₂O for a wide range of VMRs, consistent with the actual detection in the CCF at the location of the planet.

Fig. B.1 Left: CCFs in units of S /N for H2O in the atmosphere of β Pic b. Thick red line: Actual CCF for the location of the planet. Thick black line: CCF for the location at the opposite side of the star. Thin lines: CCF for injected H2O, colour coded by different volume mixing ratios as shown in the colour bar. Right: S /N of injected signals as function of VMR. The green shaded area shows at which range of VMRs the molecule would have been marginally detectable.

Appendix C CO injection

Fig. C.1 shows the results of the injection analysis for CO.

Fig. C.1 Same as Fig. B.1, for CO.

Appendix D SiO injection

Fig. D.1 shows the results of the injection analysis for SiO.

Fig. D.1 Same as Fig. B.1, for SiO.

Appendix E CH₄ injection

Fig. E.1 shows the results of the injection analysis for CH₄.

Fig. E.1 Same as Fig. B.1, for CH4.

Appendix F CO₂ injection

Fig. F.1 shows the results of the injection analysis for CO₂.

Fig. F.1 Same as Fig. B.1, for CO2.

Appendix G NH₃ injection

Fig. G.1 shows the results of the injection analysis for NH₃.

Fig. G.1 Same as Fig. B.1, for NH3.

Appendix H HCN injection

Fig. H.1 shows the results of the injection analysis for HCN.

Fig. H.1 Same as Fig. B.1, for HCN.

Appendix I C₂H₂ injection

Fig. I.1 shows the results of the injection analysis for C₂H₂.

Fig. I.1 Same as Fig. B.1, for C2H2.

Appendix J H₂S injection

Fig. J.1 shows the results of the injection analysis for H₂S.

Fig. J.1 Same as Fig. B.1, for H2S.

Appendix K O₃ injection

Fig. K.1 shows the results of the injection analysis for O₃.

Fig. K.1 Same as Fig. B.1, for O3.

Appendix L PH₃ injection

Fig. L.1 shows the results of the injection analysis for PH₃.

Fig. L.1 Same as Fig. B.1, for PH3.

References

All Tables

All Figures

Fig. 1 Flux-normalised spectral traces of the star β Pic, for one single spectral chunk. The panels are not plotted as pixeled images (which would not be possible for the bottom panel, since it is not uniformly sampled), but as one coloured dot per sampled data point. Yellow is a normalised flux of 1, blue is 0. Top panel: flux sampled by the xi and yi coordinates on the detector. Bottom panel: the same flux (no interpolation) sampled by wavelength λi and separation ρi.

In the text

Fig. 2 Example of residual data after PSF subtraction. The figure is in detector xy-space, since it can then be plotted as a non-interpolated pixeled image. Spatially (vertically), the image is centred on the stellar trace, which after PSF subtraction only remains as a narrow horizontal region of increased noise. Spectrally (horizontally), the image is centred on a strong telluric absorption feature. In this particular region, the PSF subtraction does not perform as well as along the rest of the spectrum, resulting in a vertical net positive feature extending out from the centre. At even larger separations vertically from the centre, the noise in the telluric region is generally enhanced. This feature would remain even if the PSF subtraction had been perfect, since it represents the much lower transmission inside of the telluric line than in the continuum. Planet b is ~9.5 pixels below the stellar trace, but cannot be seen without crosscorrelation techniques.

In the text

Fig. 3 Spectral coverage of the observations in this study. Four spectral settings were used: L3262 (upper left), L3340 (upper right), L3426 (lower left) and M4318 (lower right). The coloured regions show the coverage of that particular setting, while the grey regions show the coverage of the other settings. The scatter corresponds to the noise level in each spectral region. Higher scatter can be caused by shorter integration time, worse ambient conditions, or lower instrumental transmission, but here it primarily reflects wavelength ranges of higher thermal background, and/or wavelength ranges of lower atmospheric transmission.

In the text

	Fig. 4 Cross-correlation S/N map for H2O, showing a clear detection of β Pic b at its expected separation and velocity. The velocity is in the frame of reference of the planet itself. The separation changes between 536 mas and 558 mas during the observational baseline between 29 Nov. 2022 and 1 Nov. 2023.
In the text

	Fig. 5 Results from MCMC fitting of the spin velocity and limb darkening parameter. Upper left: histogram for the spin velocity. Upper right: evolution of the spin velocity (blue) and limb darkening parameter (gold) along the MCMC random walk. Lower left: correlation between spin velocity and limb darkening. Lower right: histogram for the limb darkening parameter.
In the text

	Fig. 6 Model spectra after continuum subtraction, showing where strong lines occur for four different molecules: H2O in blue, CO in red, CH4 in green, and SiO in orange. The grey regions are the wavelength ranges covered by our observations.
In the text

Fig. 7 Comparison of the signal injection scheme compared to the faster ‘pseudo-injection’ scheme described in Sect. 3.4. Top: direct comparison of the injected spectra using regular injection (black line) and pseudo-injection (red line). Also shown as grey dots are extracted spectral points at the location of the injection. Bottom: CCF outputs using injection (black line) and pseudo-injection (red line), compared to the actual water CCF for planet b (dashed line). The injection and pseudo-injection are both well representative of the real signal.

In the text

Fig. 8 Cross-correlation functions of β Pic b (left) and c (right). The black lines are the cross-correlation functions at the locations of the planets, while the grey lines are the corresponding functions at the exact opposite side of the star relative to the planets. Note the different scaling on the y-axes between the two panels. Planet b is clearly detected, while planet c shows a bump that can only be seen as a tentative feature.

In the text

Fig. 9 Cross-correlation results (normalised to units of S/N) for different ways of combining the individual data sets. In each case, the black dashed line represents the result when combining all the data together. Top left: orange lines are the CCFs from each individual observation taken in the L3262 setting. The red line is the combination of all L3262 observations. No dramatic improvement in S/N is acquired from this combination. Top right: same, but for the case of L3340. Middle left: same, for the case of L3426. Middle right: same, for the case of M4318. Bottom: here, a set of single observations from each spectral setting is plotted as orange lines. Their combination is plotted as a red line. In this case, a significant improvement is gained from combining the data sets, almost reaching the S/N level from the full data set.

In the text

Fig. A.1 Theoretical spectra after continuum subtraction of molecules that we have searched for in the cross-correlation analysis. A similar figure for H2O, CO, CH4 and SiO is shown in the main text, Fig. 6. From top to bottom, the additional molecules shown here are CO2 in cyan, NH3 in magenta, HCN in light blue, C2H2 in dark red, H2S in pink, O3 in purple, and PH3 in orange. Grey areas correspond to the wavelength windows covered by our observations.

In the text

Fig. B.1 Left: CCFs in units of S /N for H2O in the atmosphere of β Pic b. Thick red line: Actual CCF for the location of the planet. Thick black line: CCF for the location at the opposite side of the star. Thin lines: CCF for injected H2O, colour coded by different volume mixing ratios as shown in the colour bar. Right: S /N of injected signals as function of VMR. The green shaded area shows at which range of VMRs the molecule would have been marginally detectable.

In the text

	Fig. C.1 Same as Fig. B.1, for CO.
In the text

	Fig. D.1 Same as Fig. B.1, for SiO.
In the text

	Fig. E.1 Same as Fig. B.1, for CH4.
In the text

	Fig. F.1 Same as Fig. B.1, for CO2.
In the text

	Fig. G.1 Same as Fig. B.1, for NH3.
In the text

	Fig. H.1 Same as Fig. B.1, for HCN.
In the text

	Fig. I.1 Same as Fig. B.1, for C2H2.
In the text

	Fig. J.1 Same as Fig. B.1, for H2S.
In the text

	Fig. K.1 Same as Fig. B.1, for O3.
In the text

	Fig. L.1 Same as Fig. B.1, for PH3.
In the text