© The Authors 2025

1 Introduction

The combination of high-contrast imaging (HCI) and spectroscopy offers access to the orbital and atmospheric properties of low-mass companions and super-Jovian exoplanets. Colloquially referred to as direct spectroscopy or high-contrast spectroscopy, this technique relies on the separation of the light coming from the companion from that of its host star using adaptive optics (AO) and post-processing algorithms that disentangle their significantly different spectral signatures (Snellen et al. 2015).

Early results of direct spectroscopy include the detection of carbon monoxide in the atmosphere of the giant exoplanet HR 8799 c (Konopacky et al. 2013) using the OH-Suppressing Infrared Integral Field Spectrograph (OSIRIS; Quirrenbach et al. 2003; Larkin et al. 2006) instrument at the Keck telescope as well as the measurement of the spin velocity of β Pic b (Snellen et al. 2014) using the Cryogenic High-Resolution Infrared Echelle Spectrograph (CRIRES; Kaeufl et al. 2004) from the Very Large Telescope (VLT). Since then, many studies have continued to develop this technique further: OSIRIS has enabled the measurement of the radial velocity (RV) of the HR 8799 b and c planets (Ruffio et al. 2019) and the constraint of the atmospheric parameters of the c planet (Wang et al. 2023). The Spectrograph for INtegral Field Observations in the Near Infrared (SINFONI; Eisenhauer et al. 2003) at the VLT was used to measure the atmospheric ¹²CO/¹³CO isotopologue ratio of the young giant exoplanet TYC 8998-760-1 b (Zhang et al. 2021). Using the NIRSPEC spectrograph (McLean et al. 1998; Lopez et al. 2020; Fitzgerald et al. 2018) from the Keck Planet Imager and Characterizer (KPIC; Delorme et al. 2021), Xuan et al. (2022) measured the carbon-to-oxygen ratio (C/O) and metallicity of the brown dwarf HD 4747 B. More recently, Ruffio et al. (2024) measured the 2.9−5.2 μm spectrum of the T dwarf HD 19467 B using the NIRSpec (Jakobsen et al. 2022; Böker et al. 2023) integral field spectrometer (IFS) on board the James Webb Space Telescope (JWST), leading to the detection of H₂O, CH₄, CO₂, and CO in its atmosphere (Hoch et al. 2024). Beyond the exploration of the atmospheric physics of gas giant exoplanets, the atmospheric composition of low-mass companions might provide insights into their formation history (see, e.g., Öberg et al. 2011; Öberg & Bergin 2021; Turrini et al. 2021; Pacetti et al. 2022). While the details of planet formation are still mostly unclear, these measurements, which are enabled by newer instruments and more powerful analysis tools, might help to put together some pieces of the puzzle.

The latest facility offering high-contrast spectroscopic capabilities on an 8 m class telescope is the Enhanced Resolution Imager and Spectrograph (ERIS; Davies et al. 2023). Installed at the Cassegrain focus of the Unit Telescope 4 (UT4) at the VLT, ERIS is the new AO-assisted instrument replacing the Infrared imaging and spectroscopic capabilities previously offered by NAOS-CONICA (NACO; Lenzen et al. 2003; Rousset et al. 2003) and SINFONI (Eisenhauer et al. 2003). With its refurbished and upgraded IFS SPIFFIER (George et al. 2016) as well as its new AO module (Riccardi et al. 2022), ERIS is expected to improve upon its predecessor SINFONI in terms of contrast, angular separation, and sensitivity. However, its competitiveness in the current landscape of high-contrast spectrographs remains to be proven. In this article, we aim to demonstrate the scientific potential of this new instrument by reporting on the detection of one of the most challenging directly imaged exoplanets discovered so far, AF Lep b, which sits at an angular separation of 0′′.32 and contrast of ΔK = 11.84 mag (De Rosa et al. 2023).

The super-Jovian exoplanet AF Lep b was simultaneously discovered by three teams who were independently following up on the significant astrometric accelerations of its host star measured between the Gaia and Hipparcos catalogues (De Rosa et al. 2023; Mesa et al. 2023; Franson et al. 2023). The exoplanet was later precovered in archival NACO observations taken 11 years prior (Bonse et al. 2025), and since then has been re-observed with the JWST Near Infrared Camera (NIRCam; Franson et al. 2024) and GRAVITY at the VLT Interferometer (VLTI; Balmer et al. 2025). Its dynamical mass was measured at 3.75 ± 0.5 M_jup (Balmer et al. 2025), making it one of the lowest-mass directly imaged exoplanets found to date. By combining the relative astrometry measured by GRAVITY with all archival RV and astrometric data, most of the orbital parameters could be accurately determined: its period is ${24.3}_{-0.4}^{+0.9}$ yr, its semi-major axis is ${8.98}_{-0.08}^{+0.15}$ au and its eccentricity is smaller than 0.02 (Balmer et al. 2025). AF Lep is an F8V star (Gray et al. 2006) at 26.8 pc with a mass of ${1.13}_{-0.04}^{+0.10}$ M_⊙ (De Rosa et al. 2023). It is a likely member of the β Pictoris moving group (Zuckerman et al. 2001; Malo et al. 2013; Ujjwal et al. 2020) with an age of 24 ± 3 Myr (Bell et al. 2015). The latest atmospheric study – based on the currently available photometry and low-resolution spectroscopy – locate the planet at the L–T transition and indicate an effective temperature of T_eff = 800 ± 50 K with evidence of elevated metallicity and chemical disequilibrium (Balmer et al. 2025), which agree with previous atmospheric constraints (Zhang et al. 2023; Palma-Bifani et al. 2024; Franson et al. 2024).

The combined fit on the GRAVITY relative astrometry together with RV and astrometric data was not able to constrain two of the orbital parameters AF Lep b (Balmer et al. 2025). Since direct imaging, and hence interferometry, projects the astrometric path of astronomical objects on the celestial sphere, it is insensitive to the component of the velocity of the planet that is orthogonal to the sky plane, and therefore results in a bimodal distribution for the ascending node and argument of periapsis. This ambiguity prevents the prediction of the orbital phase and the determination of the true obliquity of the planet. The former is relevant for planning follow-up observations in reflected light, while the latter can inform the planet formation history of AF Lep b (Kraus et al. 2020). Orbital motion perpendicular to the sky plane can only be determined by RV measurements. Due to the fast rotation of the star, existing stellar RV measurements of AF Lep are not accurate enough to robustly resolve the reflex motion of the star caused by the planet (De Rosa et al. 2023). Therefore, the only way to uniquely disambiguate the orbital solution is by measuring the RV of the planet directly.

In this article, we present new ERIS/SPIFFIER K-band observations of AF Lep b at moderately high spectral resolution (R ∼ 11 000) obtained within the ERIS Guaranteed Time Observations (GTO) programme allocated to ETH Zurich. We describe the observations and our data reduction (Sects. 2.1 and 2.2), including our custom wavelength calibration (Sect. 2.2.1) as well as our algorithm to remove the stellar point spread function (PSF; Sect. 2.2.3). We compute molecular maps (Hoeijmakers et al. 2018) of the AF Lep system (Sect. 3.1), revealing the planet through the detection of water and carbon monoxide in its atmosphere. We then describe how we measured the RV of AF Lep b relative to its star (Sect. 3.2), which allowed us to resolve the ambiguity in its orbital solution (Sect. 4.1). We discuss the implications of our detection of H₂O and CO (Sect. 4.2), the impact of stellar contamination (Sect. 4.3), and compute the detection limits for our observations (Sect. 4.4). In Sect. 5, we provide our concluding remarks.

Fig. 1 Illustration of the main processing steps applied to one frame of our dataset. Left: raw detector frame. The brick-wall pattern of the slitlets is visible as the telluric absorption lines are shifted in adjacent slitlets. Middle: reduced frame. The missing data in the middle is due to cropping the frames, essentially removing slitlets at the border of the image. Right: residuals after PSF subtraction via spectral PCA. The two horizontal lines at 2.24 and 2.30 μm mark the position of the mask applied to the data to hide the strong stray light features present in slitlets 2 to 6 and 12 to 15 (cf. Fig. C.1). The slitlets are numbered from 1 to 32, and span the image from bottom to top (i.e. from south to north).

2 Methods

2.1 Observations

AF Lep was observed with ERIS on the 8 November 2023 using SPIFFIER with the smallest image scale (0^′′.8 × 0^′′.8) and the K-long grating (2.19 μm to 2.47 μm at R ∼ 11 000). The observations were taken in field-tracking mode (pupil-tracking is not available at this time) with position angle set to 0^∘ (i.e. north at the top of the image); the star was offset by 0′′.15 to the west to include the exoplanet in the field of view, with a square jitter pattern of width 0′′.1 with offsets every second frame, and sky offsets after each square jitter. We recorded a total of 267 object frames with a detector integration time (DIT) of 40 s and n_DIT = 1. The observing conditions were average, with a median seeing of 0′′.65, a median coherence time of 5.7 ms, and a median integrated water vapour column of 3.32 mm.

2.2 Data reduction

The data were reduced using the ERIS-SPIFFIER pipeline v1.6.0 together with custom tools to improve the quality of the calibration. The standard calibration steps consist of dark subtraction, flat fielding, bad pixel correction, distortion correction, and an initial wavelength solution using daytime arc lamp calibrations. The pipeline normally relies on the OH sky emission lines present in the science data to obtain a good wavelength calibration (Davies 2007); however, our observations are dominated by the stellar PSF, and the OH lines are not visible in the data, which prevented us from using this calibration strategy. Therefore, we had to implement a custom tool to calibrate the wavelength using the telluric absorption lines visible in the stellar PSF. We give a short summary of our custom wavelength calibration in Sect. 2.2.1, while the full description can be found in Appendix A. We describe the additional pre-processing steps necessary to format the output of the ERIS-SPIFFIER pipeline in Sect. 2.2.2, and our PSF subtraction technique in Sect. 2.2.3. The main steps of our pipeline – including the PSF subtraction – are illustrated in Fig. 1. In the following, we use the term spaxel to denote the spectrum extracted from one cube at a given spatial pixel.

2.2.1 Custom wavelength calibration

The ERIS instrument is installed in the Cassegrain focus of the UT4 telescope. Consequentially, internal flexures inevitably impact the position of the spectrum on the detector as the instrument rotates and changes altitude to track its target. This renders daytime wavelength calibrations inaccurate and necessitates leveraging the science frames to produce an accurate wavelength solution. For observations of faint targets with long integration times and large field of views, the ERIS-SPIFFIER pipeline makes use of the OH sky emission lines as references to compute an accurate wavelength solution (Davies 2007) for each science frame. As mentioned above, this is not possible in our case due to the stellar PSF dominating the full frame and hiding the OH emission lines. Instead, we can leverage the telluric absorption lines imprinted on the stellar spectrum to improve the wavelength calibration. As reference, we used a telluric transmission template calculated with SkyCalc (Noll et al. 2012; Jones et al. 2013) and evaluated it with the parameters matching our observing conditions. Our custom wavelength calibration consists of three steps, which are summarised below.

Step 1: We first measure an initial shift in each spaxel of the data by cross-correlating each spaxel with the telluric transmission template. This first step accounts for the spectral curvature, which causes a line of constant wavelength to follow a curve across the detector.

Step 2: We measure higher-order wavelength errors. To do this, we first correct the spectral curvature measured in the first step and combine all spaxels within the same slitlet. We then measure the local error by cross-correlating the data and the telluric transmission template on a moving window along the wavelength axis. This measurement yields the wavelength error as a function of wavelength, which we fit using a spline.

Step 3: We update the wavelength solution. This is done by adding the wavelength error measured in the previous steps to the initial wavelength map computed by the ERIS-SPIFFIER pipeline from arc lamp images taken during daytime calibrations. We then run the eris_ifu_jitter recipe together with the updated wavelength map, which produces the calibrated datacubes.

Figure 2 shows a small region of the mean spectrum contained within the slitlet 11 before and after our wavelength calibration together with the telluric transmission template for comparison. The improvement is clearly visible by eye: the telluric lines present in the calibrated data match the template much more closely. In Appendix A, we give a more exhaustive description of our wavelength calibration procedure (see Fig. A.1 for an illustration of the two main steps) and estimate its accuracy.

Fig. 2 Application of our custom wavelength calibration to the spectrum. Upper panel: spectrum calibrated using daytime arc lamp images as provided by the standard pipeline when OH sky emission lines are hidden by the stellar PSF. Lower panel: spectrum calibrated using our custom pipeline. We provide the telluric transmission template used for the calibration as reference (orange). We only show the region between 2.375 μm to 2.4 μm to easily distinguish the telluric lines by eye.

2.2.2 Pre-processing

After the wavelength calibration, each exposure is a 3D cube consisting of 2155 λ-images of size 64 × 64, covering a wavelength range between 2.187 99 μm and 2.468 14 μm. We implemented a custom pipeline to prepare the data for post-processing. Due to the brick-wall arrangement of the slitlets on the detector (see the left image in Fig. 1), we excluded all λ-images outside the range 2.202 29 μm to 2.451 76 μm to avoid any missing data. Furthermore, we observed the presence of bright stray light extending between 2.24 μm to 2.30 μm and within slitlets 2 to 6 and 12 to 15 (see Fig. 1). The cause for this stray light within ERIS/SPIFFIER is unclear. Since it effectively adds strong correlated noise along the spectral and spatial dimensions, it can hinder the detection of a molecular signal in the atmosphere of a companion or create fiducial features that might be mistaken for one. Therefore, we masked the full spectral region between 2.24 μm to 2.30 μm. We further removed 4 pixels at the edge of each frame to avoid contamination by adjacent slitlets due to inaccuracies of the distortion calibration, leaving λ-images of size 56 × 56. Due to residual instrument flexures between the AO stage and the SPIFFIER detector at the moment of our observations, the target was slowly drifting over time at a rate of approximately 60 mas h⁻¹. This drift forced us to pause the observation after 1.5 h to re-centre the star. Nevertheless, we still had to remove 28 frames for which the star drifted outside the field of view, leaving 239 frames. The position of the star in the remaining frames was measured by fitting the stellar PSF with a Gaussian in the wavelength-averaged frames. At this point, the data had reached their final shape; there were 239 cubes left, each made of 1923 λ-images of size 56 × 56.

After PSF subtraction but before cross-correlating the data with spectral templates (cf. Sections 2.2.3 and 3.1), we convolved each image with an aperture of radius 1.7 px (21.25 mas), i.e., we applied aperture photometry on every pixel. The output of this step is a datacube with the same format as the input cube, but where each spaxel now contains the total flux over its surrounding aperture. This step was carried out to gather all the signal coming from the planet and thereby increasing the signal-to-noise ratio (S/N) of its spectrum. We tried several radii for the apertures, including 1 λ/D (i.e. 58.5 mas); however, they yielded worse results, probably because of remaining systematics from the PSF subtraction that are included in the apertures. The spectrum of the planet is extracted after this convolution by taking the spaxel resulting in the highest S/N in the molecular maps.

2.2.3 Removal of the stellar point spread function with spectral PCA

Since the stellar PSF dominates the signal from the exoplanet, we needed to model and subtract it from the data. As pupil-tracking was not available at the time of our observations, we could not use angular differential imaging (ADI; Marois et al. 2006). Instead, we wanted to leverage the significant differences in spectral features (which can be resolved at high enough spectral resolution) between the companion and the tellurics-imprinted stellar spectrum. This was able to be done by using principal component analysis (PCA; see e.g. Bishop & Nasrabadi 2006).

PCA is widely used by the exoplanet community to suppress systematic noise. Depending on the type of data, different schemes have been developed to apply PCA. In general, PCA models the strongest covariances along what is known as the feature dimension using a set of data points called samples. In the context of HCI, PCA was first introduced to ADI (Amara & Quanz 2012; Soummer et al. 2012). In this framework, PCA explores the spatial covariances of the speckle noise utilising the image dimensions (X × Y) as features and the time dimension (T) as samples. An alternative scheme to apply PCA is PCA-Temporal (Long et al. 2023) or similarly TRAP (Samland et al. 2021). Instead of exploring the spatial covariances, PCA-Temporal models the temporal covariances using the time dimension (T) as features and the image dimensions (X × Y) as samples. For IFS data at low spectral resolution (R < 1000), the wavelength dimension (Λ) can be added to the samples to leverage angular and spectral differential imaging simultaneously (ASDI/CODI; Kiefer et al. 2021). Spectral differential imaging exploits the outward motion of speckles with increasing wavelength due to chromatic PSF broadening to differentiate them from the fixed signal of the planet.

At higher spectral resolution, where the wavelength coverage might not suffice for SDI, the signal from the planet can be disentangled from that of its host star by leveraging their significantly different spectral features. This simple idea is implemented in different ways in the literature; however, most methods follow a similar scheme: (i) first create a stellar PSF model using expert knowledge about HCI and telluric absorption lines, and then (ii) correct any remaining systematics using PCA. Hoeijmakers et al. (2018) implemented this scheme on K-band SINFONI data of β Pic b at R ∼ 4000 and first coined the term molecular maps. Their method first built a PSF model by dividing the 3D IFS data by a master stellar spectrum extracted from the brightest pixels. The result was then low-pass filtered and multiplied by the master stellar spectrum. The role of the low-pass filter is to create a model of the wavelength-dependent stellar PSF; however, it also includes the continuum of the planet, which therefore disappears when the PSF is subtracted. Since telluric absorption lines of CH₄ remained in the residuals, they applied the SYSREM algorithm (Tamuz et al. 2005; Birkby et al. 2013), which reduces to PCA if the data have constant uncertainties, in order to remove these systematics. In this framework, the features are defined along the spectral dimension (Λ) whilst the samples consist of all spaxels of a frame along the two spatial dimensions (X × Y). The same method, where SYSREM was replaced with PCA, was later applied to VLT/MUSE and SINFONI data of the PDS 70 system (Haffert et al. 2019; Cugno et al. 2021), where it was coined high-resolution spectral differential imaging (HRSDI). Ruffio et al. (2019) introduced a forward-model framework to simultaneously fit the planet spectrum together with the stellar PSF, which is constructed as in Hoeijmakers et al. (2018). This approach was later improved (Hoch et al. 2020; Ruffio et al. 2021) to include PCA to model the remaining systematics after an initial stellar PSF fit and subtraction. In this framework, the features are defined as 5 × 5 spaxel patches ($\mathcal{P}$ × Λ, with $\mathcal{P}$ ⊂ X × Y), whilst the samples are defined from the rest of the frame. These two frameworks explicitly enforce the behaviour of the stellar PSF as a function of wavelength: the low-frequency variations come from the chromatic broadening of the stellar PSF and the continuum of the stellar spectrum, whilst the high-frequency variations are due to the telluric absorption lines imprinted on the stellar spectrum. Any remaining systematics due to wavelength calibration errors or variable telluric line strengths are subsequently identified and removed with PCA.

For this work we opted to leave the full modelling of the stellar PSF to PCA (steps (i) and (ii) combined). We defined the features along the spectral dimension (Λ) and the samples consist of all spaxels in the frame (X × Y). We computed a new PCA basis separately for each cube, meaning that the time dimension (T) is not used. Afterwards, each spaxel of each cube was projected onto the first 250 principal components (PCs), thereby creating a model of the stellar PSF for each cube. The number of principal components used to build the PCA basis was determined by a blind search. The planet can be detected with 50 PCs, but the highest S/N was obtained with 250 PCs. Since our approach uses PCA for the full PSF model (i and ii) and not only for the subtraction of the remaining systematics (ii) as in (e.g. Hoeijmakers et al. 2018; Hoch et al. 2020; Ruffio et al. 2021), a significantly larger number of components is needed. The PCA-based noise model was subtracted from the data resulting in a sequence of residual cubes. These residuals were then aligned with the previously derived position of the star in each cube, and median-combined along the time dimension to produce a final residual cube. We refer to this technique as spectral PCA.

Although PCA has been used with features along the spectral axis (Λ) and samples along the two spatial axes (X × Y) before (Hoeijmakers et al. 2018), the novelty of this technique consists in using it to fully model the 3D stellar PSF. To underline the similarities and differences between some of the existing PSF subtraction techniques and Spectral PCA, we report their defining concepts in Table 1.

3 Results

3.1 Molecular maps

Molecular maps are obtained by cross-correlating each spaxel in the residual cube with molecular spectral templates, thereby optimally matching the molecular signature from the exoplanet atmosphere while dampening any remaining systematics. Our spectral templates were calculated using petitRADTRANS (Mollière et al. 2019) with only one molecule as a source of opacity at a time, setting its mass fraction to −2 dex, and adopting the maximum a posteriori pressure–temperature profile constrained by Balmer et al. (2025). We calculated spectral templates for H₂O, CO, CH₄, and CO₂ as they produce significant spectral features in the K band. Additionally, we prepared a full-model template corresponding to the maximum a posteriori parameters derived by Balmer et al. (2025) for AF Lep b (i.e. log g = 4.17 dex, Fe/H = 0.74, C/O = 0.50) containing H₂O, CO, CH₄, CO₂, H₂S, HCN, NH₃, FeH, TiO, PH₃, K, VO, and Na, but ignoring the cloud treatment. We then removed the continuum of the spectral templates using a Gaussian filter with a standard deviation of 2.6 nm. We note that the precise size of the filter does not significantly influence the molecular maps as long as the continuum is effectively removed; however, we did not systematically search for the optimal filter size to create the molecular maps. We cross-correlated the final residual cube with the spectral templates by Doppler-shifting them over a range of radial velocities from −1000 to 1000 km s⁻¹ in steps of 1 km s⁻¹ using the function crosscorrRV from the Python package PyAstronomy (Czesla et al. 2019). Following Petrus et al. (2021), molecular maps were created by evaluating the cross-correlation function (CCF) at the RV of the planet. As an additional step, we standardised the map by subtracting the mean and dividing by the standard deviation over the whole image.

We show the resulting molecular maps in Fig. 3. The planet is clearly revealed in the H₂O, CO, and full-model maps at S/N values of 9.3, 9.7, and 20.7, respectively. The calculations made to quantify the S/N follow Petrus et al. (2021) with the modification to only consider values of the CCF at the same angular separation as the planet as the speckle noise is expected to depend on the angular separation (see e.g. Marois et al. 2008). More details on the computation of the S/N metric are given in Appendix B. The signal is located at a position of Δα ≈ + 0′′.30 and Δδ ≈ + 0′′.06 (i.e. ρ ≈ 0′′.3, θ ≈ 79^∘), as expected from previous direct imaging of the planet (De Rosa et al. 2023; Mesa et al. 2023; Franson et al. 2023; Bonse et al. 2025; Balmer et al. 2025), and shows a (non-calibrated) RV of approximately 20 km s⁻¹.

The extracted spectrum of the planet is shown in Fig. 4 together with our full model template for reference. Although the planet spectrum has a low S/N and might still contain residual stellar light, absorption lines are nevertheless well resolved, for example in the region from 2.30 μm to 2.32 μm where CO dominates the spectrum. The stray light mentioned in Sect. 2.2 can be seen between 2.24 μm and 2.30 μm. The spectral templates used for the cross-correlation are shown in Fig. C.1 together with the spectra at the position of the planet before and after PSF subtraction and aperture photometry. Finally, we show the cross-correlation functions with the different molecules at the position of the planet in Fig. 5.

Fig. 3 PSF intensity (top left) and molecular maps of the AF Lep system. The planet AF Lep b is revealed in the H2O and CO maps, and with the full model from Balmer et al. (2025), but not in the CH4 and CO2 maps. The orange star indicates the position of the star, while the white dashed circle, whose radius is equal to 1 λ/D, indicates the position of the planet.

3.2 Radial velocity

To measure the RV from the cross-correlation function, we would like to interpret it as the position of the peak in the cross-correlation function. However, each frame might have slightly different systematics. Errors in the wavelength solution, imperfect PSF subtraction, and stray light, for example, might affect the shape of the cross-correlation function as well as the position of its peak. Since the planet is not detected in single frames, we adopted a bootstrapping approach (Beran 2008) to replicate the observation many times by randomly sampling the frames and combining them, thereby quantifying the random effect of the instrument and processing the systematics. Furthermore, we corrected for the relative motion of the Paranal observatory with respect to the AF Lep system. Since we detected the planet with three different spectral templates, we can measure its RV with each model. The exact details of the measurement and calibration of the RV of AF Lep b are given in Appendix C. We report the following results for each model: for H₂O, we obtain $\mathrm{\Delta }{v}_{\mathrm{R},\mathrm{P}\star }^{\left({\mathrm{H}}_2\mathrm{O}\right)}=8.9\pm 2.5$ km s⁻¹; for CO, we obtain $\mathrm{\Delta }{v}_{\mathrm{R},\mathrm{P}\star }^{(\mathrm{C}\mathrm{O})}=6.1\pm 2.8$ km s⁻¹; and for the full model, we obtain $\mathrm{\Delta }{v}_{\mathrm{R},\mathrm{P}\star }^{(\text{Full\hspace{0.17em}})}=7.8\pm 1.7$ km s⁻¹. These three values are consistent with each other and allow the orbital solution to be constrained. We note that since the full model contains H₂O and CO, the three measurements are not independent of each other, which might explain why the RV value obtained with the full model lies between the other two values. This also makes the combination of the three values non-trivial. Instead, we chose the measurement resulting in the highest S/N detection of AF Lep b (i.e. the full model), and take it as our final RV measurement $\mathrm{\Delta }{v}_{\mathrm{R},\mathrm{P}\star }=\mathrm{\Delta }{v}_{\mathrm{R},\mathrm{P}\star }^{(\mathrm{F}\mathrm{u}\mathrm{l}\mathrm{l})}=7.8\pm 1.7$ km s⁻¹.

For a planet that has been discovered simultaneously by three independent teams, it is not surprising that its RV has also been measured by another instrument at a time similar to that of our ERIS observations. Denis et al. (2025) have also observed AF Lep b with the new VLT/HiRISE (Vigan et al. 2024) instrument at high spectral resolution (R ∼ 100 000). Their analysis delivers the value 10.51 ± 1.03 km s⁻¹ (Denis et al. 2025). Our final RV measurement lies within 1.6 σ of their result.

To validate the uncertainty delivered by our bootstrapping approach, we ran a simple simulation using the exposure time calculator of ERIS/SPIFFIER¹. For a field object with an effective temperature of 800 K and an apparent magnitude of 16.77 mag such as AF Lep b, the measured spectrum is estimated to have a S/N ≈ 10 for an identical observing set-up, as described in Sect. 2.1. We simulate 10⁴ noisy spectra at this S/N using our full model of the planet and subsequently measure the RV using the same cross-correlation technique as described in Appendix C (i.e. by measuring the position of the peak of the cross-correlation function evaluated with our full model spectral template). The resulting RV measurements have a standard deviation of 0.1 km s⁻¹. This reflects the uncertainty with which the RV can be measured for this experiment. Since the S/N calculated with the ETC ignores residual systematics from the subtraction of the stellar PSF as well as photon noise from the stellar PSF overlapping with AF Lep b, the S/N is expected to be much smaller. By repeating the experiment with S/N = 1, the resulting uncertainty is 1.1 km s⁻¹, which is already much closer to our uncertainty of 1.7 km s⁻¹. The actual S/N of our extracted spectrum is difficult to estimate, especially considering that the continuum is lost during PSF subtraction. However, since CO absorption lines are visible by eye, the S/N cannot be far from 1. Additionally, as we describe in Sect. 2.2.1, measuring the wavelength error after applying our pipeline showed a residual error of 0.8 ± 0.3 km s⁻¹ over the whole observation (i.e. a bias of 0.8 km s⁻¹ and a random scatter of 0.3 km s⁻¹). Assuming that the bias can be corrected perfectly (see Appendix C for details on the correction of this bias), the random errors arising from the wavelength calibration add another uncertainty of the order of 0.3 km s⁻¹. Our actual uncertainty of 1.7 km s⁻¹ can therefore be quite well explained as a combination of residual wavelength errors and low S/N spectrum.

Fig. 4 Spectrum of AF Lep b after PSF subtraction (black) compared to our full model (orange). The CO-dominated region between 2.30 μm and 2.32 μm matches the model particularly well despite the low S/N of the data. The shaded area represents the mask applied to the data to remove the instrument stray light described in Sect. 2.2.

Fig. 5 Cross-correlation functions with molecular templates at the position of the planet in Fig. 3. The peaks indicate the presence of H2O and CO in the atmosphere of AF Lep b at an (uncalibrated) RV of approximately 20 km s−1 (vertical line). The cross-correlation functions were standardised over the spatial axes in the same way as our molecular maps: the cross-correlation map at each RV was mean-subtracted and normalised by its standard deviation over the whole image.

4 Discussion

4.1 Orbital constraint

Our RV measurement can be used to remove the ambiguity in the previously constrained orbital solution of AF Lep b (De Rosa et al. 2023; Franson et al. 2023; Mesa et al. 2023; Bonse et al. 2025; Zhang et al. 2023; Balmer et al. 2025). The orbital fit obtained by Balmer et al. (2025) on all previous imaging, astrometry, stellar RV, and GRAVITY data can be evaluated to predict the RV of the companion relative to its host star on the night of our ERIS observations. Due to the bimodality of the orbital solution, there are two equiprobable radial velocities: (a) Δ v_R,P⋆ = 10 ± 2 km s⁻¹ associated with the argument of periapsis $\omega ={109}_{-21}^{+13}$ deg and longitude of the ascending node $\mathrm{\Omega }={248.8}_{-0.7}^{+0.4}$ deg, and (b) Δ v_R,P⋆ = −10 ± 2 km s⁻¹ associated with the other solution $\omega ={289}_{-21}^{+13}$ deg, $\mathrm{\Omega }={68.8}_{-0.7}^{+0.4}$ deg (the other orbital parameters are all equal in both cases). Here, the ascending node Ω is defined in the same way as in orvara (Brandt 2021); in other words, it denotes the position angle at which the planet crosses the sky plane moving towards the observer (we note that orbitize! uses the opposite definition; Blunt et al. 2024).

Our RV measurement using the full model template: Δ v_{R,P ⋆} = 7.8 ± 1.7 km s⁻¹. It is consistent with the first orbital solution within 1.3 σ, while excluding the other at 11 σ. Our measurement therefore shows that the orbital solution $\omega ={109}_{-21}^{+13}$ deg, $\mathrm{\Omega }={248.8}_{-0.7}^{+0.4}$ deg is the correct one. This means that (i) the relative RV of AF Lep b was close to its maximum value at the epoch of our observations, (ii) the planet recently crossed the sky plane moving behind the star (i.e. away from the Earth), and (iii) more generally the planet is located behind the star (from the point of view of the Earth) between the position angles ${68.8}_{-0.7}^{+0.4}$ deg and ${248.8}_{-0.7}^{+0.4}$ deg. In particular, it is now possible to compute a phase curve of AF Lep b, which is crucial for the planning of follow-up observations in visible light. Moreover, this delivers one of the two puzzle pieces necessary to determine the true obliquity of AF Lep b (i.e. the angle between its orbital angular momentum and the spin axis of the star). The other piece, the orientation of the spin axis of the star, can be determined by measuring the displacement of the photocentre due to stellar rotation (Le Bouquin et al. 2009; Kraus et al. 2020).

4.2 Implications for the atmospheric chemistry of AF Lep b

In itself, our detection of H₂O and CO in the atmosphere of AF Lep b is not surprising as these molecules have already been found in similar objects. Several atmospheric studies have started to directly constrain the mass fractions of H₂O, CO, and CH₄ in such objects (Zhang et al. 2023; Palma-Bifani et al. 2024; Balmer et al. 2025). Chemical disequilibrium, which is driven by vertical mixing, has been advanced as a potential explanation for the high abundance of CO found by these studies. At an effective temperature of 800 K, CH₄ is expected to be chemically favoured over CO (Mollière et al. 2022), hence the need for an upward stream of CO from hotter deeper layers. Therefore, our detection of H₂O and CO together with our non-detection of CH₄ are compatible with the chemical disequilibrium hypothesis. To investigate whether our new ERIS observations also favour this hypothesis over chemical equilibrium, one would need to constrain the molecular abundances of H₂O and CO. This requires the use of a dedicated atmospheric retrieval framework such as CROCODILE (Hayoz et al. 2023) and will be investigated in a future article.

4.3 Stellar and telluric contamination

The stellar PSF is still very strong at the position of the planet: around 10^−3.25 times its maximum value. This means that the planet with its contrast of ΔK = 11.84 mag (De Rosa et al. 2023) is around 30 times fainter than the PSF at its position. There is therefore a risk of contamination of the planet signal by the stellar spectrum. More precisely, since the spectrum of AF Lep only contains atomic lines and at most very weak molecular lines, the risk instead comes from the strong telluric absorption lines imprinted in the stellar spectrum by the atmosphere of the Earth. However, molecular maps remedy this problem by measuring the cross-correlation at the position of the planet and everywhere else in the field of view. Therefore, if the spectral templates correlate with the stellar spectrum, then this contamination can be measured as a cross-correlation signal in a different place than the planet. This is the reason for our chosen S/N metric, which ensures that any contamination by the tellurics-imprinted stellar spectrum is taken into account.

4.4 High-contrast capabilities of ERIS

The high significance of our detection of AF Lep b hints at deeper detection limits. Establishing a robust detection limit for molecular mapping is more complicated than for imaging, due to the high dependence on the atmospheric properties of the companions and the spectral type of the host star. The pressure-temperature profile, the chemical composition of the atmosphere, and the presence and composition of clouds strongly affect the shape and strength of the absorption lines, which in turn will affect the signal picked up by cross-correlation with spectral templates (Bidot et al. 2024). At the same time, the photosphere of the host star can (depending on its spectral type) contain the same molecules of interest as the companion, which might affect the efficacy of the PSF subtraction and thereby the detection limits, hindering the generalisation to another stellar system (Bidot et al. 2024). Accounting for all these factors is beyond the scope of this paper. Nevertheless, we wanted to provide indicative detection limits that can be used as a reference point for the planning of future observations with ERIS/SPIFFIER. To do so, we set up a fake planet injection and retrieval test as is customary in HCI (see e.g. Bonse et al. 2023), but we constrained the analysis to the detection of companions spectrally identical to AF Lep b. Concretely, we used our full model from Sect. 3.1 as planet spectrum.

We created a PSF reference by considering the six frames where the star is closest to the centre of the field of view, cropping them to a size of 0′′.3 × 0′′.3, aligning them, and finally median-combining them. To inject a fake planet at a given contrast, we divided the datacube of the PSF reference by the stellar spectrum and multiplied it by the planet spectrum and the contrast. The stellar and planet spectra were normalised using synthetic photometry with the filter SPHERE/IRDIS K12_2, as it was used to report the contrast of AF Lep b in the K band (De Rosa et al. 2023). For the stellar spectrum, we used a BT-NextGen model (Allard et al. 2012) evaluated at the grid point closest to the stellar parameters reported by the Gaia DR3 (i.e. T_eff = 5900 K, log g = 4 dex, Fe/H = −0.5 dex). Fake planets were injected into each datacube of the reduced data at a fixed relative position with respect to the host star and with no Doppler shift. We injected the fake planets right before the PSF subtraction and subsequently applied the rest of the pipeline to the fiducial datasets.

Figure 3 shows that there is not a lot of room in the 0′′.8 × 0′′.8 field of view of ERIS/SPIFFIER to inject fake planets around the host star. This forced us to inject planets at only one position angle, namely at 105 deg (as measured from north through east), and only up to 0′′.4 angular separation. Concretely, we used a range of angular separation between 1 and 7 λ/D in steps of 1 λ/D 9 (i.e. between 0′′.06 and 0′′.41). For the contrast, we selected a range between 9 and 16 Δmag in steps of 0.5 Δmag. This allowed us to create a contrast grid similar to the ones computed by applefy (Bonse et al. 2023).

The resulting contrast grids for the molecular maps of H₂O, CO, and for the full model are shown in Fig. 6. From these grids, we extracted contrast curves by setting the detection threshold to S/N = 8 to enforce a clear detection. We show the three extracted S/N = 8 detection limits in Fig. 7 together with the K-band contrast and angular separation of other directly imaged exoplanets as reference. To illustrate these contrast grids and our detection limits, we show a few resulting molecular maps for a sample of contrast and angular separation close to our detection limits in Fig. D.1. The injected fake planets are indeed clearly detected at S/N ≥ 8, but not necessarily at lower values, which validates our selected detection threshold.

These detection limits imply that we could have detected AF Lep b even if it had been at an angular separation of 2 λ/D (i.e. 0′′.12). The second row of Fig. D.1 shows that the fake planet is retrieved at S/N of 15.2 and 14.6 with the H₂O template and the full model. This corresponds to the contrast and angular separation close to which β Pic c was detected using interferometry with GRAVITY (Nowak et al. 2020). The close-in exoplanet has not been detected yet with non-interferometric observations. β Pic c certainly has different atmospheric properties than AF Lep b, which is not accounted for in this analysis. To validate the high-contrast potential of molecular mapping with ERIS/SPIFFIER in the 1−2 λ/D range, one would have to actually detect β Pic c. If possible, this would significantly extend the parameter space accessible to characterise exoplanets at moderately high spectral resolution towards smaller angular separation. Molecular mapping with ERIS/SPIFFIER might even allow us to efficiently search for other young and bright super-Jovian planets that are closer in.

We would like to emphasise that these detection limits are only indicative and should not be interpreted as a guarantee to detect a companion with a contrast and angular separation located close to the contrast curves. Many assumptions went into computing them, which is why we would like to emphasise them here. The input spectrum is a simplified model of the spectrum of AF Lep b that does not include any treatment of clouds. The wavelength is assumed to be perfectly calibrated. The S/N metric used here (cf. Petrus et al. 2021) does not take into account effects from small sample statistics and might overestimate the confidence of the detection. For each contrast and angular separation, the fake planets were only injected at one position angle, which can create a bias depending on the position of residual speckles. Finally, the contrast curves are calculated with respect to AF Lep specifically, which means that a brighter star might result in a worse detection limit at close separations due to the presence of brighter speckles, whereas observations of a fainter star might be limited by the sensitivity of ERIS/SPIFFIER. Taking all that into account, our detection limits indicate that a similar observing set-up is likely to result in the detection of an exoplanet with similar atmospheric properties as AF Lep b with a contrast and angular separation that would locate it sufficiently above the detection limits shown in Fig. 7.

Fig. 6 Contrast grids computed by injecting fake planets with our full model spectrum from Sect. 3.1 at different contrast and angular separation from the host star and retrieving them in the molecular maps of H2O (left), CO (centre), and with the full model (right). The dashed red line corresponds to a S/N of 8, which we use as the detection threshold to compute our detection limits.

Fig. 7 Detection limits at S/N = 8 for H2O, CO, and the full model for our ERIS/SPIFFIER observations of the AF Lep system. The detection limits were calculated for the full 178 min total integration time, which required just over 4 h of telescope time, and did not use ADI as it was not yet available at the time of our observations. The contrast and angular separations of the other directly imaged exoplanets indicated in this figure are taken from the following studies: β Pic c from Nowak et al. (2020), β Pic b from Bonnefoy et al. (2011), HR 8799 e from Marois et al. (2010), AF Lep b from De Rosa et al. (2023), and 51 Eri b from Samland et al. (2017).

5 Summary and conclusion

In this article, we demonstrated the first high-contrast spectroscopic application of the new ERIS/SPIFFIER instrument at the VLT. We observed the AF Lep system at moderately high spectral resolution (R ∼ 11 000) in the K band as part of the ERIS GTO programme allocated to ETH Zurich. Since high-contrast spectroscopic applications are not covered by the standard SPIFFIER pipeline, we developed a robust method to calibrate the wavelength solution by cross-correlating the strong telluric absorption lines imprinted on the stellar spectrum with a SkyCalc template (Noll et al. 2012; Jones et al. 2013) and modelling the slight spectral curvature caused by the single diffraction grating within SPIFFIER. Our custom pipeline additionally corrects high-order effects by fitting a spline to the wavelength error. Furthermore, we presented spectral PCA, an innovative application of PCA to model the 3D stellar PSF and subtract it from the data. Our findings are summarised in the following points.

• The new ERIS instrument at the Very Large Telescope can detect molecules in the atmosphere of a planet with a contrast of ΔK = 11.84 mag (De Rosa et al. 2023) and an apparent separation of 0′′.32, and can be used to measure its RV with an uncertainty of 1.7 km s⁻¹. This demonstrates that ERIS/SPIFFIER is a competitive instrument for the spectral and orbital characterisation of directly imaged exoplanets at high contrast and small angular separation.

• We directly detected water and carbon monoxide in the atmosphere of AF Lep b at high S/N. Together with our non-detection of methane, this is compatible with the hypothesis of chemical disequilibrium currently favoured by previous observations. We will further investigate the atmospheric properties of AF Lep b using atmospheric retrievals in a follow-up article together with other targets observed as part of the ERIS GTO programme.

• We measured the radial velocity of AF Lep b relative to its host star: Δ v_{R,P ⋆} = 7.8 ± 1.7 km s⁻¹. This value indicates that the orbital solution with argument of periapsis $\omega ={109}_{-21}^{+13}$ deg and longitude of ascending node $\mathrm{\Omega }={248.8}_{-0.7}^{+0.4}$ deg is the correct one, removing the bimodal distribution found by Balmer et al. (2025). As explained in Sect. 4.1, this now allows the computation of the phase curve of AF Lep b necessary for the planning of follow-up observations in visible light, and delivers one piece of the puzzle necessary for the measurement of its true obliquity.

• We computed the detection limits for our observations by injecting fake planets at different contrast and angular separations from the host star. These detection limits were derived with many assumptions listed in Sect. 4.4, the most critical being that we used a simple model of the spectrum of AF Lep b, whilst the sensitivity of molecular mapping might strongly depend on the atmospheric parameters (i.e. pressure–temperature profile, chemistry, clouds). Taken at face value, our detection limits seem to indicate that ERIS/SPIFFIER might be able to reach a contrast of ΔK = 12 mag at an angular separation of 2 λ/D (i.e. 0′′.12). If confirmed, this would significantly extend the parameter space available for the search and investigation of the atmospheric and orbital properties of low-mass companions and young exoplanets at moderately-high spectral resolution, as to date the only known exoplanet in this regime is β Pic c, which has only been directly imaged using interferometry with the VLTI/GRAVITY instrument (Nowak et al. 2020).

• Finally, ADI (Marois et al. 2006) is not strictly necessary for the detection and characterisation of exoplanets when observed at high enough spectral resolution. Pupil-tracking had not yet been commissioned for ERIS/SPIFFIER at the time of our observations. Our results demonstrate that stellar speckles can be effectively disentangled from the signal of the companion using PSF subtraction via spectral PCA. Future work should investigate whether molecular mapping can be combined with ADI to reach even deeper detection limits. However, if ADI is not necessary, it might relax some functionality requirements under investigation for the design of the next generation of high-contrast facilities aiming at characterising terrestrial exoplanets, such as the Planetary Camera and Spectrograph (PCS, Kasper et al. 2021) at the upcoming Extremely Large Telescope.

Acknowledgements

We thank the anonymous referee for their constructive feedback which helped us to greatly improve our article. JH and SPQ gratefully acknowledge the financial support from the Swiss National Science Foundation (SNSF) under project grant number 200020_200399. Part of this work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. GC thanks the Swiss National Science Foundation for financial support under grant numbers P500PT_206785 and P5R5PT_225479. This project received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 819155). Contributions: JH proposed the project, wrote the observing proposal, prepared and analysed the observations, and wrote the article. MJB provided insightful discussions on PSF subtraction techniques and helped writing the corresponding section. AAB, YC, and RD provided help with the implementation of our custom wavelength calibration. MJB, FD, and HMS provided insightful discussions about orbital fitting. The first 16 co-authors listed after JH all provided useful feedback during the project or while writing the manuscript. All other co-authors are included for their past contributions to designing, building, and commissioning the ERIS instrument, and are ordered alphabetically.

Appendix A Wavelength calibration

In this appendix, we give a full description of our wavelength calibration procedure as summarised in Sect. 2.2.1. We then estimate its accuracy and indicate the average wavelength correction applied as a function of time during the whole observation. As already stated in the main text, our custom wavelength calibration consists of three steps, which we extensively describe in the following.

Step 1: We first measure an initial shift in each spaxel of the data. This is done by cross-correlating each spaxel with the telluric template using the function phase_cross_correlation from skimage (van der Walt et al. 2014) with an upsampling factor of 20. We found that the error varied approximately linearly over the extent of each slitlet (see the upper panel in Fig. A.1). Moreover, the slope of the error also varies from the slitlets on one side of the detector to the slitlets on the other side. This effect – which is hardly noticeable by eye on raw detector frames (cf. Fig. 1) – is due to spectral curvature, which causes a line of constant wavelength to follow a curve from one side to the other. We modelled this spectral curvature by fitting the slope of the wavelength error over all slitlets using a third-order polynomial (see the lower left panel in Fig. A.1). We then integrated this third-order polynomial and fit the y-intercept of the resulting function to model the wavelength error within each slitlet (see the orange line in the upper panel in Fig. A.1).

Fig. A.1 Illustration of the steps involved in our custom wavelength calibration applied to an arbitrary frame. Panel (a): Measurement of the wavelength error across all slitlets by cross-correlating each spaxel with a telluric transmission template (black). We only show every fourth measurement for readability. Bad measurements are indicated by red dots, while the red line indicates the median error over the whole detector. The blue line indicates an initial linear fit of the error over each slitlet. The orange line represents the model created in the first step of our wavelength calibration as described in Sect. 2.2.1 and shown in Fig. A. Panel (b): Slope of the linear fit (blue line in Fig. A) to the wavelength error in each slitlet (black) as described in step 1 of Sect. 2.2.1. The orange line shows the third-order polynomial fit to the data. Panel (c): Measurement of the wavelength error along the wavelength dimension within slitlet 11 (black) as described in step 2 of Sect. 2.2.1. The orange line shows the spline fit to the data. The shaded areas correspond to the slitlets at the top and bottom of the field of view which are cropped out after the wavelength calibration. The figures at the top and on the bottom left illustrate step 1 whilst the figure on the bottom right shows step 2.

Step 2: We then measure higher-order wavelength errors. Since the telluric lines cover the whole extent of the K-long grating, we can locally measure the wavelength error by cross-correlating the telluric transmission template on subintervals of the spectrum. To do so, we first combine all spaxels within the same slitlet after correcting them by the initial error measured previously. We then define subintervals of size 200 wavelength channels that we shift by 50 channels along the wavelength dimension such that they overlap. The size of the intervals is chosen to ensure that enough telluric lines are contained for a robust measurement, while the overlaps increase the sampling along the wavelength dimension. We then cross-correlate each subinterval with the telluric transmission template. This yields a list of measurements of the wavelength error as a function of wavelength (see lower right panel in Fig. A.1 for the calibration of the slitlet 11 taken as an example). We fit the error with a third-order spline and a smoothing factor of 0.4 using the function UnivariateSpline from SciPy (Virtanen et al. 2020). The smoothing factor was chosen by inspecting the resulting fits by eye and ensuring that noisy error measurements did not impact the spline. For the 100 channels at either side of the spectrum, we set the wavelength error as equal to the value of the spline at the last measurement point. This is done for the sake of robustness because the spline can grow polynomially on either side and can result in extreme values.

Step 3: Finally, we update the wavelength solution. The standard pipeline provides an initial wavelength map computed from arc lamp images taken during daytime calibrations. This map relates each detector pixel to its corresponding wavelength. We update it by adding the errors from the two previous steps according to the column number and wavelength bin associated to each pixel. The updated wavelength map is then used together with the normal recipe eris_ifu_jitter.

To validate our calibration, we repeat the measurement of the wavelength error described in steps 1 and 2 above after applying our custom calibration. Figures A.2 and A.3 present the results for the same frame as in Fig. A.1: most of the structure present in the uncalibrated data (cf. Fig. A.1) has been corrected, including both the linear trend across the slitlets and the wavy pattern along the wavelength axis. Unfortunately, some slitlets – more precisely slitlets 12, 17, and 18 – still show some trends of the order of 0.25 to 0.5 px; however, the trend is significantly smaller than before the calibration. From a median error of 3.85 px (i.e. 0.5 nm or 64 km s⁻¹) initially measured for the frame considered in Fig. A.1, our pipeline delivered a calibration with a median residual error of 0.10 px (i.e. 0.01 nm or 1.7 km s⁻¹). To obtain an overview of the quality of the calibration, we repeated this exercise by measuring the median residual error of each frame after having applied our custom calibration. Over the whole observation, the median of the wavelength error is 0.05 px and the standard deviation 0.02 px (i.e. 6.5 ± 2.4 pm or 0.8 ± 0.3 km s⁻¹). These results validate our custom wavelength calibration routine.

Fig. A.2 Validation of our custom wavelength calibration. The measurement described in step 1 of Sect. 2.2.1 was repeated after applying the wavelength calibration (black) to the same frame as in Fig. A.1. We only show every fourth measurement for readability. The shaded areas correspond to the slitlets at the edges of the field of view and are cropped out after the wavelength calibration.

Over the whole observation, the wavelength solution varied by a total of 5 px, i.e. 0.65 nm or 83 km s⁻¹, as shown in Fig. A.4. This large drift of the wavelength solution over time is due to internal flexures within SPIFFIER as the telescope – and thus the Cassegrain focus where ERIS is installed – changes altitude to track its target. The light-weighted diffraction gratings of SINFONI suffered from cryogenic deformation (George et al. 2017), which significantly impacted its spectral profile. Therefore, heavier gratings were used for SPIFFIER which helped to greatly improve the spectral profile at the cost of significant flexures with changes in altitude and rotation. The impact of the flexures in SPIFFIER can be mitigated with a good wavelength calibration strategy – either by using OH sky emission lines (Davies 2007) for faint targets and long integration times as is done by default by the ERIS-SPIFFIER pipeline, or by leveraging telluric absorption lines for observations of bright targets such as presented in this work.

Fig. A.3 Validation of our custom wavelength calibration. The measurement described in step 2 of Sect. 2.2.1 was repeated after applying the wavelength calibration (black).

Fig. A.4 Median correction applied to the wavelength solution shown as a function of time since the start of the observation. The star was recentred after 100 min, hence the break. The significant drift is due to flexures of the diffraction gratings of ERIS/SPIFFIER with changes in altitude and rotation.

Appendix B Statistical test for the detection of molecules

In this appendix, we describe the statistical test that we used to interpret the values of the cross-correlation function (CCF) as a measure for the confidence in the presence of molecular signal in each spaxel. This is done by setting a detection threshold on the values of the CCF based on a desired confidence level. Several such metrics have been used in the literature, as described in Garvin et al. (2024). We choose the test introduced by Petrus et al. (2021), as it has been used a few times already (Patapis et al. 2022; Mâlin et al. 2023), with a small modification to evaluate the CCF only at the angular separation of the planet. Concretely, we consider only the pixels of the standardised CCF (i.e. mean-subtracted and normalised by the standard deviation) at the RV of the planet (i.e. 20 km s⁻¹) and in an annulus of width 1 λ/D (i.e. 58.5 mas) and radius equal to the angular separation of the planet (0.30). We also mask out the 1 λ/D region containing the planet. This results in a total of 214 px considered for the noise estimate. We fit a Gaussian function to the histogram of their values. The S/N is then given by S/N = (C_p − μ)/σ, where μ and σ are the mean and standard deviation of the fitted Gaussian, and C_p is the value of the CCF at the position of the planet. Figure B.1 illustrates our measurement for the maps of H₂O, CO, and for the full model.

Fig. B.1 Top: molecular maps of H2O, CO, and full model. The dashed annulus highlights the region considered to quantify the S/N of the detection. The dashed circle represents the position of AF Lep b and is masked for the purpose of quantifying the S/N. Bottom: histograms of the values of the CCF contained in the annular region. The dark line shows the Gaussian fit (with parameters μ and σ shown in the legend), while the orange line shows the value Cp of the CCF at the position of the planet.

Appendix C Measurement of the radial velocity

In this appendix, we describe the measurement of the RV of the planet.

Preparation of the data: we extracted two sets of spectra. We used an aperture of radius 1.7 px (i.e. 21.25 mas) centred on the planet in the data both before and after PSF subtraction. The reason for extracting the stellar spectrum at the position of the planet before PSF subtraction is to investigate remaining calibration errors of the wavelength solution. We removed the continuum of the stellar spectrum using a Gaussian filter of size 40 (i.e. 5.2 nm).

Preparation of the models: as a template for the telluric, we calculated a transmission spectrum of the Earth atmosphere with the same weather conditions as in the middle of our observations using SkyCalc-iPy,² a Python wrapper to access ESO’s SkyCalc tool (Noll et al. 2012; Jones et al. 2013). To obtain precise measurements of the RV, we used spectral templates at a higher bin sampling rate than the data (we chose a factor of 10) and convolved them down to the spectral resolution of the instrument using a Gaussian kernel with size $\sigma =(2\sqrt{2\ln 2}{)}^{-1}(\lambda /R)=0.055$ nm, where λ = 2.327 μm is the average wavelength of the K-long grating and R = 11 000 is the spectral resolution of ERIS/SPIFFIER. The continuum of the planet model was then removed using a Gaussian filter with the same physical size as the one used for the data, i.e. 5.2 nm. We calculated the CCFs separately for each frame on the range −80 to 80 km s⁻¹ and in steps of 0.1 km s⁻¹. The spectral templates used to measure the RV of AF Lep b are shown in Fig. C.1 together with the extracted spectrum of the planet before and after PSF subtraction.

Fig. C.1 Spaxel at the position of the planet before and after PSF subtraction with spectral PCA and aperture convolution (red and black), and spectral templates used to compute the molecular maps presented in Fig. 3. All spectra are shown after continuum subtraction. Additionally, they are normalised and offset vertically for readability. The shaded area between 2.24 μm and 2.3 μm represents the masked wavelength range where stray light is present.

Propagation of the uncertainty: as described in the main text (see Sect. 3.2), we adopted a bootstrapping approach (Beran 2008) to estimate the uncertainty arising from the systematics affecting the CCF. We describe the exact process in the following. We started by extracting the 239 planet spectra after PSF subtraction (i.e. before frame combination). We then cross-correlated these spectra with the spectral templates of H₂O, CO, and the full model. We then picked 10 000 bootstrap samples for each spectral template, each consisting of a random choice with replacement of the 239 CCFs. For each sample, we mean-combined the CCFs and extracted the RV as the peak of the combined CCF. We noticed that, after bootstrapping, the measured RVs of the planet and the tellurics depended on the size of the filter and the method (Gaussian or median) used to remove the continuum. To remove this source of bias, we repeated the bootstrapping framework for both filter types and with filter sizes between 10 and 100 in steps of 2 (i.e. 1.3 nm and 13 nm in steps of 0.26 nm), and pooled the resulting 112 sets of samples together, essentially weighting each filter equally.

Intermediate results: we obtained the following values for the RV of the planet and for the spectral templates of H₂O, CO, and the full model: $v_{\mathrm{R},\mathrm{P}}^{\left({\mathrm{H}}_2\mathrm{O}\right)}$ = 20.6 ± 2.4 km s⁻¹, $v_{\mathrm{R},\mathrm{P}}^{(\mathrm{C}\mathrm{O})}$ = 17.8 ± 2.7 km s⁻¹, $v_{\mathrm{R},\mathrm{P}}^{(\mathrm{F}\mathrm{u}\mathrm{l}\mathrm{l})}$ = 19.5 ± 1.6 km s⁻¹ and for the RV of the tellurics v_{R, Tell} = 5.0 ± 0.2 km s⁻¹. At first glance, the fact that the tellurics are not at rest in the data seems to invalidate the extensive wavelength calibration discussed in Sect. 2.2.1. We report a systematic wavelength error of the order of 0.8 km s⁻¹ over the whole observation, which is six times smaller than the measurement of the RV of the tellurics reported here. We can resolve this contradiction by comparing the position of the planet in the molecular maps (Fig. 3) with the residual wavelength error presented in Fig. A.2. The planet is unfortunately located within the slitlet 18, which still shows some trend in the wavelength error from one side to the other side of the slitlet. Along slitlet 18, the remaining wavelength error after calibration varies by 0.5 px, i.e. 8.5 km s⁻¹. The planet being on the side of the slitlet, this causes our wavelength calibration to systematically leave a residual error in that location in the image. Fortunately, we can easily correct for this effect by subtracting the RV of the tellurics measured at the position of the planet from the RV of the planet. Future works could explore how to more accurately model the spectral curvature to better stabilise the wavelength calibration across the field of view of SPIFFIER in the cases not covered by the standard pipeline, i.e. when the stellar PSF dominates the image.

Calibration of the radial velocity: to obtain the RV of AF Lep b relative to its host star, we first subtracted the RV of the tellurics to express our measurement in the rest frame of the observatory. We then transformed to the rest frame of the barycentre of the Solar System by correcting for the rotation of the Earth, the motion of the Earth around the Earth–Moon barycentre, and the motion of the Earth–Moon barycentre around the Sun. To do so, we used the function helcorr from PyAstronomy (Czesla et al. 2019) at the Julian date of the middle of our observation 2 460 256.729 24 d. This yielded a correction of v_R,BC = −13.96 km s⁻¹ (as measured towards AF Lep). We note that the correction varied between −14.18 and −13.70 km s⁻¹ between the start and end of our observations. We did not track the mean Julian date within each sampled frame used for our bootstrapping procedure; however, we expect the mean Julian date within each sample to vary by less than 3% due to the large number of frames (239) from which the samples are selected randomly, thus producing an additional uncertainty of the order of 0.015 km s⁻¹ (i.e. two orders of magnitude below our measurement uncertainty). For the RV of AF Lep, we adopted the value measured by GAIA: v_{R, ⋆} = 20.61 ± 0.51 km s⁻¹ (Brown et al. 2018). Finally, the RV of AF Lep b relative to its host star is given by $\mathrm{\Delta }{v}_{\mathrm{R},\mathrm{P}\star }^{(i)}=v_{\mathrm{R},\mathrm{P}}^{(i)}$ − v_{R, Tell} − v_R,BC − v_R,⋆, where i denotes any of the three spectral templates used (i.e. H₂O, CO, and the full model).

Appendix D Selected sample of fake planet injection tests

Fig. D.1 Sample of fake planet injection tests to support the detection limits extracted from the contrast grids computed in Sect. 4.4. Each row shows a different experiment, and each column shows the corresponding molecular map computed with a different spectral template. The contrast and angular separation of the injected planets are annotated in white at the top of each figure. The dotted white circles show the position of the injected planets as well as the annulus considered to compute the S/N. The molecular maps show the cross-correlation values at the RV of the injected planet, i.e. 0 km s−1, which is why they look slightly different than Fig. 3. The actual planet, AF Lep b, is still clearly visible due to the broad peak of the CCF (cf. Fig. 5).

References

All Tables

All Figures

Fig. 1 Illustration of the main processing steps applied to one frame of our dataset. Left: raw detector frame. The brick-wall pattern of the slitlets is visible as the telluric absorption lines are shifted in adjacent slitlets. Middle: reduced frame. The missing data in the middle is due to cropping the frames, essentially removing slitlets at the border of the image. Right: residuals after PSF subtraction via spectral PCA. The two horizontal lines at 2.24 and 2.30 μm mark the position of the mask applied to the data to hide the strong stray light features present in slitlets 2 to 6 and 12 to 15 (cf. Fig. C.1). The slitlets are numbered from 1 to 32, and span the image from bottom to top (i.e. from south to north).

In the text

Fig. 2 Application of our custom wavelength calibration to the spectrum. Upper panel: spectrum calibrated using daytime arc lamp images as provided by the standard pipeline when OH sky emission lines are hidden by the stellar PSF. Lower panel: spectrum calibrated using our custom pipeline. We provide the telluric transmission template used for the calibration as reference (orange). We only show the region between 2.375 μm to 2.4 μm to easily distinguish the telluric lines by eye.

In the text

	Fig. 3 PSF intensity (top left) and molecular maps of the AF Lep system. The planet AF Lep b is revealed in the H2O and CO maps, and with the full model from Balmer et al. (2025), but not in the CH4 and CO2 maps. The orange star indicates the position of the star, while the white dashed circle, whose radius is equal to 1 λ/D, indicates the position of the planet.
In the text

	Fig. 4 Spectrum of AF Lep b after PSF subtraction (black) compared to our full model (orange). The CO-dominated region between 2.30 μm and 2.32 μm matches the model particularly well despite the low S/N of the data. The shaded area represents the mask applied to the data to remove the instrument stray light described in Sect. 2.2.
In the text

Fig. 5 Cross-correlation functions with molecular templates at the position of the planet in Fig. 3. The peaks indicate the presence of H2O and CO in the atmosphere of AF Lep b at an (uncalibrated) RV of approximately 20 km s−1 (vertical line). The cross-correlation functions were standardised over the spatial axes in the same way as our molecular maps: the cross-correlation map at each RV was mean-subtracted and normalised by its standard deviation over the whole image.

In the text

	Fig. 6 Contrast grids computed by injecting fake planets with our full model spectrum from Sect. 3.1 at different contrast and angular separation from the host star and retrieving them in the molecular maps of H2O (left), CO (centre), and with the full model (right). The dashed red line corresponds to a S/N of 8, which we use as the detection threshold to compute our detection limits.
In the text

Fig. 7 Detection limits at S/N = 8 for H2O, CO, and the full model for our ERIS/SPIFFIER observations of the AF Lep system. The detection limits were calculated for the full 178 min total integration time, which required just over 4 h of telescope time, and did not use ADI as it was not yet available at the time of our observations. The contrast and angular separations of the other directly imaged exoplanets indicated in this figure are taken from the following studies: β Pic c from Nowak et al. (2020), β Pic b from Bonnefoy et al. (2011), HR 8799 e from Marois et al. (2010), AF Lep b from De Rosa et al. (2023), and 51 Eri b from Samland et al. (2017).

In the text

Fig. A.1 Illustration of the steps involved in our custom wavelength calibration applied to an arbitrary frame. Panel (a): Measurement of the wavelength error across all slitlets by cross-correlating each spaxel with a telluric transmission template (black). We only show every fourth measurement for readability. Bad measurements are indicated by red dots, while the red line indicates the median error over the whole detector. The blue line indicates an initial linear fit of the error over each slitlet. The orange line represents the model created in the first step of our wavelength calibration as described in Sect. 2.2.1 and shown in Fig. A. Panel (b): Slope of the linear fit (blue line in Fig. A) to the wavelength error in each slitlet (black) as described in step 1 of Sect. 2.2.1. The orange line shows the third-order polynomial fit to the data. Panel (c): Measurement of the wavelength error along the wavelength dimension within slitlet 11 (black) as described in step 2 of Sect. 2.2.1. The orange line shows the spline fit to the data. The shaded areas correspond to the slitlets at the top and bottom of the field of view which are cropped out after the wavelength calibration. The figures at the top and on the bottom left illustrate step 1 whilst the figure on the bottom right shows step 2.

In the text

	Fig. A.2 Validation of our custom wavelength calibration. The measurement described in step 1 of Sect. 2.2.1 was repeated after applying the wavelength calibration (black) to the same frame as in Fig. A.1. We only show every fourth measurement for readability. The shaded areas correspond to the slitlets at the edges of the field of view and are cropped out after the wavelength calibration.
In the text

	Fig. A.3 Validation of our custom wavelength calibration. The measurement described in step 2 of Sect. 2.2.1 was repeated after applying the wavelength calibration (black).
In the text

	Fig. A.4 Median correction applied to the wavelength solution shown as a function of time since the start of the observation. The star was recentred after 100 min, hence the break. The significant drift is due to flexures of the diffraction gratings of ERIS/SPIFFIER with changes in altitude and rotation.
In the text

Fig. B.1 Top: molecular maps of H2O, CO, and full model. The dashed annulus highlights the region considered to quantify the S/N of the detection. The dashed circle represents the position of AF Lep b and is masked for the purpose of quantifying the S/N. Bottom: histograms of the values of the CCF contained in the annular region. The dark line shows the Gaussian fit (with parameters μ and σ shown in the legend), while the orange line shows the value Cp of the CCF at the position of the planet.

In the text

Fig. C.1 Spaxel at the position of the planet before and after PSF subtraction with spectral PCA and aperture convolution (red and black), and spectral templates used to compute the molecular maps presented in Fig. 3. All spectra are shown after continuum subtraction. Additionally, they are normalised and offset vertically for readability. The shaded area between 2.24 μm and 2.3 μm represents the masked wavelength range where stray light is present.

In the text

Fig. D.1 Sample of fake planet injection tests to support the detection limits extracted from the contrast grids computed in Sect. 4.4. Each row shows a different experiment, and each column shows the corresponding molecular map computed with a different spectral template. The contrast and angular separation of the injected planets are annotated in white at the top of each figure. The dotted white circles show the position of the injected planets as well as the annulus considered to compute the S/N. The molecular maps show the cross-correlation values at the RV of the injected planet, i.e. 0 km s−1, which is why they look slightly different than Fig. 3. The actual planet, AF Lep b, is still clearly visible due to the broad peak of the CCF (cf. Fig. 5).

In the text