© The Authors 2025

1 Introduction

In the wake of the large-scale surveys conducted with groundbased extreme adaptive optics (AO) planet imagers such as GPI at Gemini-South (Macintosh et al. 2014) and SPHERE at the VLT (Chauvin et al. 2017), direct imaging has recently achieved a fundamental milestone in the selection and followup of sources using the astrometric acceleration measurements information provided by the Gaia telescope (Kervella et al. 2019; Brandt 2021; Kervella et al. 2022). This astrometric information considerably increases the success rate of surveys based on this prior knowledge, compared to fully blind searches (Nielsen et al. 2019; Vigan et al. 2021), as illustrated by several imaging discoveries of substellar companions over the recent years (Bowler et al. 2021; Bonavita et al. 2022; Currie et al. 2023; Rickman et al. 2024). It also offers the unique possibility of inferring the dynamical masses of the companions, which, in turn, can be used to directly constrain the formation and evolution models of giant planets (Brandt et al. 2021a; Zhang 2024).

An emblematic discovery in that context concerns the accelerating star AF Lep (HD 35850, HR 1817, HIP 25486) around which the discovery of a young giant planet (hereafter, AF Lep b) was announced almost simultaneously by three independent teams (De Rosa et al. 2023; Mesa et al. 2023; Franson et al. 2023). The detection was obtained in direct imaging, but based on a selection of the target using the long time baseline astrometric acceleration, or proper-motion anomaly, measured between the Hipparcos and Gaia missions (Brandt 2021).

AF Lep is a 1.09 ± 0.06 M_⊙ star, of spectral type F8 (Gray et al. 2006), located at a distance of 26.825 ± 0.014 pc (Bailer-Jones et al. 2021), with a super-solar metallicity ([Fe/H] = 0.29 ± 0.03, estimated from a high S/N HARPS spectrum, Perdelwitz et al. 2024). As member of the β Pictoris moving group, the isochronal age of the star has been estimated at 24 ± 3 Myr (Bell et al. 2015). The mass and orbital elements of the young giant planet AF Lep b were recently estimated by combining all available direct imaging data, including an early detection from 2011 obtained from archival VLT/NaCo data (Bonse et al. 2024) and VLTI/GRAVITY data (Balmer et al. 2025): it has an estimated dynamical mass of ${3.68}_{-0.48}^{+0.47}$ M_Jup, and is orbiting at a semi-major axis of ${8.98}_{-0.08}^{+0.15}$ au with an eccentricity $e={0.013}_{-0.024}^{+0.010}$ and an inclination of $i={57.5}_{-0.7}^{+0.6}\circ$.

Recent chemically consistent retrievals for AF Lep b using petitRADTRANS (Mollière et al. 2019) of all previously published emission spectra and photometry spanning 0.9–4.2 µm, with a maximum resolution of R_λ = 30, have confirmed a cool temperature of T_eff ≃ 800 K and low surface gravity of log(g) ≃ 3.7 dex. The retrievals also confirmed the presence of silicate clouds and disequilibrium chemistry in the atmosphere of AF Lep b, as expected for a young, cold, early-T super-Jovian planet (Zhang et al. 2023). This analysis also revealed a metal-enriched atmosphere ([Fe/H] > 1.0 dex), compared to the host star’s metallicity, supported by follow-up study including JWST/NIRCam observations (Franson et al. 2024). These results were independently confirmed by Palma-Bifani et al. (2024), who used a forward-modelling approach with the Exo-REM atmospheric radiative-convective equilibrium model, which includes the effects of non-equilibrium processes and clouds (Charnay et al. 2018). The results of the latter study point towards a lower metallicity of 0.5–0.7 and a C/O ratio solution ranging between 0.4 and 0.8, compatible with a solar to super-solar value (C/O_⊙ = 0.55, Asplund et al. 2009).

Here, we report the results of much higher resolution spectroscopic observations of AF Lep b made with the new HighResolution Imaging and Spectroscopy of Exoplanets (HiRISE, Vigan et al. 2024) instrument at the Very Large Telescope (VLT), which combines the SPHERE exoplanet imager (Beuzit et al. 2019) with the recently upgraded high-resolution spectrograph CRIRES+ (Dorn et al. 2023). HiRISE operates in the H-band at a spectral resolution on the order of R_λ = 140 000.

In Sect. 2, we present the observations and calibration data related to AF Lep. In Sect. 3, we describe the HiRISE data reduction and signal extraction steps applied to the data, in particular, the corrections for the tellurics and instrumental response and the wavelength calibration. The strategy used to model the data is presented in Sect. 4. In Sect. 5, we report the results of our forward modelling atmospheric analysis. In Sect. 6, we present new constraints coming from the radial velocity (RV) measurements of the planet itself and their implications for the orbital solutions derived in combination with previous studies. Finally, we discuss our results in Sect. 7 and present our conclusions in Sect. 8.

2 Observations

The AF Lep system was observed on UT 2023 November 20 and 23 with VLT/HiRISE. These observations are summarised in Table 1. At each epoch, we started the observations by placing the host star on the science fiber to acquire a reference spectrum. We did not use a coronagraph for the observations to maximize the end-to-end transmission of the instrument (Vigan et al. 2024). The centering of the star on the single-mode fiber was optimized using a dedicated procedure performed on the internal source of the SPHERE instrument (El Morsy et al. 2022; Vigan et al. 2024), which provides a typical centering accuracy better than 8 mas (0.2λ/D in the H band). The CRIRES+ spectrograph was set up in the H1567 spectral setting and the detector integration time (DIT) of the science detector was set to 120 sec, resulting in ∼20 000 ADU per exposure, the saturation limit of the CRIRES⁺ detector being 37000ADU (see Table 2 of the CRIRES⁺ user manual). Then, the centering procedure was used again to place the planetary companion AF Lep b on the science fiber. For this, we relied on the astrometric calibration of the tracking camera obtained during the HiRISE commissioning (Vigan et al. 2024) and on the accurate astrometry of the companion measured a few days earlier by VLTI/GRAVITY (Balmer et al. 2025, private communication). For the observations of the companion, the DIT of the science detector was set to 1200 s, resulting in ∼60 ADU per exposure. This DIT was chosen as a compromise between the number of detector readouts over a long exposure, which will increase the overall noise level, and the risk of losing an exposure in case of unstable observing conditions. During all science observations, the internal metrology of CRIRES⁺ is enabled to provide an improved precision in the wavelength solution.

After the science exposures, we also acquired sky backgrounds on an empty region of the sky about 30″ away from AF Lep. The goal of these exposures is mainly to help subtract the leakage term associated to the HiRISE MACAO guide fiber described in Vigan et al. (2024). Although this leakage term has been decreased by a factor ∼10 since commissioning by the addition of an optical attenuator on the guide fiber, it is still visible in long exposures and needs to be properly subtracted before the star and companion’s spectra can be extracted. We note that this strategy is not optimal in very low signal regime because subtracting the same background to all science exposures tends to increase the noise level. Fortunately the leakage term is stable, and the 1200 s DIT was used for all targets during the HiRISE observing run, so we were able to combine the different backgrounds to decrease the noise in the backgrounds. The impact of the number of backgrounds on the derived parameters of the planet is discussed in Appendix B.

The two epochs on AF Lep are mostly identical, except that on the second night we used an offset procedure where the science fiber is moved by ∼10 pixels along the slit axis to help reduce the impact of the bad pixels on the CRIRES⁺ science detectors. This procedure intends to move the science signal with respect to static bad pixels, but the two offsets are not subtracted to each other like in classical nodding. This is why an offset of a few pixels is enough since the full width at half-maximum (FWHM) of the science signal is 2 pixels. On both nights the companion was observed at very low airmass, but the observing conditions during the first night were more stable: the DIMM seeing on the first night remained at 0.6″ with very little variations, while on the second night it varied between 0.7″ and 1.0″. The coherence time, τ₀, was very similar during the two nights, with a value between 4.6 and 5.1 ms. The end-to-end transmission of the system was measured on the star, with 95^th percentile values of 4.0% and 3.8% on the first and second nights, respectively. These values are consistent with instrumental predictions, which confirm that the centering of the stellar PSF on the science fiber was within the specifications of 0.2 λ/D in H-band.

Standard CRIRES⁺ calibrations were acquired automatically the next morning based on the science observations of each night. This includes dark, flat fields and wavelength calibration files, as detailed in the CRIRES⁺ calibration plan. HiRISE observations do not require any specific internal calibrations in addition to standard daily calibrations.

3 Data reduction

3.1 Raw data reduction and signal extraction

A dedicated pipeline was developed to reduce and calibrate HiRISE data (Costes et al. 2024a)¹. This pipeline relies on the official CRIRES⁺ pipeline provided by ESO² to produce the flat-field and wavelength calibrations, and on custom Python routines for data combination, spectral extraction, filtering and wavelength recalibration. We briefly describe below the main reduction steps performed with our Python tools.

The first step is to clean the raw data, independently for the star and the companion. The appropriate backgrounds are median-combined and subtracted to the science frames, and the science frames are divided by the detector flat field generated by the CRIRES⁺ pipeline. The bad pixels, flagged by the CRIRES⁺ pipeline, are all replaced by NaN values in the images. Finally, if multiple science frames are available, they are also mean- combined to increase the signal-to-noise ratio (S/N) of the data.

The second step is to locate the trace position of the eight orders dispersed on the three science detectors of CRIRES⁺. For this step we use the science data acquired on the star, which is usually at high S/N at the location of the science fiber. For each of the 24 orders segment, we fit a 1d Gaussian function in each of the 2048 spectral channels to provide an accurate position of the trace on the detector and of its FWHM. Some regions of the spectrum are strongly affected by the absorption of telluric lines, which may result in a poor fit. We remove these channels using an iterative sigma-clipping that detects the outliers with respect to a parabolic fit to the position of the trace. The algorithm converges in 3 iterations and the final trace position is defined as the result of the parabolic fit to the trace with all the outliers removed.

The third step is the extraction of the stellar and companion signal in each spectral channel. In each spectral channel, the signal is summed in a 6-pixel window for the companion and the star, centered around the trace position measured at the previous step. We refer to Fig. A.1 for more details. The noise is estimated as the standard deviation in a 20-pixel window located 60 pixels away from the location of the science fiber. In this spectral extraction step, the impact of bad pixels is estimated by computing a weight parameter in each spectral channel. The signal coming through the science fiber is Gaussian, with a standard deviation in the spatial direction of ∼0.85 pix that is stable over all segments of orders. In each spectral channel, the weight is computed by the integrated value of the normalized Gaussian of standard deviation 0.8 pix, centered at the calibrated position of the trace, with a value of zero attributed to the pixels flagged as bad. The weight is equal to 1 when there are no bad pixels in the extracted window, 0 when there are only bad pixels, and intermediate values between 0 and 1 when there are some bad pixels in the extracted window. This method only estimates the impact of bad pixels but cannot be used to compensate for them. For now, in the data analysis, we simply remove all data that result in an extracted window containing at least one bad pixel.

A fourth optional step can be executed when the observations have been performed with an offset procedure, for instance, in the case of our observations of AF Lep on the second night. In this case, the spectra obtained at the two offset position will be affected by different bad pixels and can be combined. The weight vectors computed at the previous step are used to help with the combination: channels where the weights are equal to 1 in both channels will be averaged together, while channels where weight of one spectrum is equal to 1 and the other is not will use the value of the spectrum not affected by bad pixels. Channels where both spectra have a value <1 will be flagged as bad and discarded in the analysis. We found that typically <1% of channels are bad in both spectra in the observation of AF Lep on the second night.

Finally, a filtering step is applied to the extracted signals for the star and the companion to remove strong outliers that may be the result of bad pixels that were not properly flagged in the initial calibrations. This filtering step usually removes less than 0.05% of the data.

Figure 1 illustrates the data extracted on the star and companion from the 2023-11-23 observations. The two offsets have been combined to remove the impact of bad pixels. The spectra are mostly shaped by (i) the overall transmission of the system, which gives an overall parabolic shape, (ii) the blaze function of the spectrograph that adds an additional parabolic shape to each of the eight individual orders, and (iii) the deep telluric lines that affect the data below 1500 nm, above 1750 nm and in some specific regions in between. As the scientific data are located at specific positions of the fibers on the detector, in background exposures without strong starlight diffraction, the main part of the detector contains only the detector’s read-out noise, allowing us to estimate this noise. We estimate a S/N of ∼25 per spectral channel at 1600 nm for the data acquired at the location of the companion. Although the spectrum obtained at the location of the companion is still dominated by the stellar PSF coupling into the science fiber, we refer to the spectrum obtained at that location as the “spectrum of the companion” for practicality.

Fig. 1 Data for the star and companion obtained on 2023-11-23 and extracted using our custom pipeline. The read-out noise of the detector is also plotted to demonstrate that the data obtained on the companion has a mean S/N of ∼25 per spectral channel at the center of the H band. The effect of the telluric lines is clearly visible at the start and end of the band, but also close to the center where many telluric CO bandheads are visible. We remind that data for the companion refers to the spectrum at the location of the companion which contains primarily stellar flux.

3.2 Sky and instrumental response

The sky and instrumental transmissions are an important parameter for the modelling of the data and the search for the planet’s signal (see Sect. 4). At high spectral resolution, the sky imprints many telluric lines over the spectra with varying width and depth. Then, the telescope and instrument also have an impact in shaping the signal, although this instrumental contribution is mostly smooth with wavelength. The main contributors that do not have a flat or almost flat contribution are the SPHERE dichroic filter and the CRIRES⁺ blaze function (see Fig. 10 of Vigan et al. 2022). Finally, the CRIRES⁺ science detectors have numerous bad pixels that have an impact on the extracted signal that needs to be taken into account or modelled.

We use the observations of dedicated early-type stars to compute the telluric and instrumental response, in that case βPicA (A6) on 2023-11-20 andAFLepA itself (F8) on 202311-23. The spectra measured for the calibrators are divided by a PHOENIX stellar model from Husser et al. (2013) at the appropriate T_eff and log g, rotationally broadened using the measured υ sin i for the stars and velocity-shifted by the known RV of the stars. This division effectively removes the stellar effect, leaving mostly the telluric and instrument effects. The response is then normalized to have a median value of one.

3.3 Wavelength recalibration

Another critical aspect of high-resolution spectroscopy is the wavelength calibration of the data. Previous authors have already demonstrated the importance of adopting an additional level of correction to the wavelength solution provided by the CRIRES+ pipeline (e.g., Landman et al. 2024; Nortmann et al. 2025).

For the HiRISE observations of AF Lep, we used a correction based on the telluric lines imprinted into the spectrum of the star. The main limitation of this approach is that telluric lines are not evenly distributed over the H band and there are even small parts of the band where no telluric lines are detectable. This prevents us from using a quadratic wavelength correction in most segments of orders (as in Landman et al. 2024), so we finally settled on a simple constant correction for each of the segments. The resulting corrections are showed in Fig. 2. They are all below 0.02 nm, and even below 0.005 nm for 20 of the 24 segments. The associated RVs of the correction are up to 3 km s⁻¹ as shown in Fig. 2. This demonstrates the importance of the wavelength recalibration. To compute the wavelength correction coefficients for each of the segments, we use a nested sampling algorithm (Skilling 2006). For the loglikelihood function, we use the CCF mapping defined in Zucker (2003). More details on the nested sampling algorithm are given in Sect. 5.1. In conjunction with the wavelength constant corrections, the nested sampling algorithm estimates the uncertainties of these corrections. In the telluric dominated segments, where we have high confidence in our wavelength recalibration, we estimate uncertainties in the RV correction on the order of 100 m/s.

In the future, a more evolved recalibration is foreseen, based on dedicated observations of an M star calibrator. This approach is used for the KPIC instrument (e.g., Wang et al. 2021) and has demonstrated good accuracy (Morris et al. 2020; Ruffio et al. 2023; Horstman et al. 2024). For HiRISE, an M star calibrator was observed every night, but we detected some shifts in wavelength between the calibrator and the science target during the first night, which are not fully understood yet. For the present work, we considered the self-calibration of the data using the telluric lines safer, since the same deep telluric lines are observed at identical positions in both the AF Lep A and AF Lep b data. Finally, the recalibrated wavelength has been corrected for the mean barycentric component computed for the time of the science exposures. The value was computed using the helcorr function of the PyAstronomy package (Czesla et al. 2019).

Fig. 2 Corrections to the wavelength solution for each segments of orders based on the analysis of the telluric lines for the AF Lep data acquired on 2023-11-23. For the two segments centered at 1556 and 1623 nm (light red), there are not enough telluric lines to compute an accurate correction and the default wavelength solution is therefore adopted. The right axis of the main plot shows the RV shift corresponding to the level of correction of the left axis, computed at 1600 nm. The top-left inset shows the effect of the correction at 1666.5 nm (grey vertical line in the main plot), which is a segment that requires one of the largest corrections. The telluric model obtained from SkyCalc is overplotted.

4 Data modeling

As explained in Sect. 3, the extracted 1D spectrum at the location of the companion is dominated by starlight diffraction and speckles coupling into the science fiber. Since we do not have an independent measurement of this stellar contamination (e.g., from a measurement at the same separation but with a different position angle), we must rely on a joint estimation of the stellar contamination and the planetary model. This approach is feasible because the two components have very difference spectral shapes. To achieve this, we adopt a method similar to Landman et al. (2024). The signal extracted from the science fiber at the location of the companion can be decomposed into a planetary signal term, d_p , a starlight contamination term, d_s , and an additional noise term, η: $$d(\lambda )=d_p(\lambda )+d_s(\lambda )+\eta ,$$(1)

where d is the reduced 1D spectrum at the location of the companion. Following Landman et al. (2024), we can modulate the stellar contamination at the location of the companion with the equation $$d_s(\lambda )=c_s\alpha (\lambda )f_s(\lambda ),$$(2)

where c_s is a scaling factor, α is a low-order function and f_s is the stellar master spectrum, that is the observation of the star when the science fiber is centered on the star. Following Landman et al. (2024), we directly estimate this modulation from the data, which corresponds to the following final model for the starlight contamination: $$d_s(\lambda )=c_s\frac{ℒ[d(\lambda )]}{ℒ\left[f_s(\lambda )\right]}{f}_s(\lambda ),$$(3)

where L is a lowpass filtering operation. For this filtering operation, we used a Savitzky-Golay filter of order 2 with a kernel size of 301 pixels. This value was empirically determined by maximizing a cross-correlation function (CCF) with template models for several HiRISE targets (see e.g. Sect. 5.3 for AF Lep b). We note that a similar value was also determined and used by Landman et al. (2024) on CRIRES⁺ data. The final results are not very sensitive to the value adopted for this parameter as long as the continuum is removed.

Similarly to Landman et al. (2024), the planetary contribution d_p can be written as $$d_p(\lambda )=c_p\left(M_{p,LSF}(\lambda )T(\lambda )-f_s\frac{ℒ\left[T(\lambda )M_{p,LSF}(\lambda )\right]}{ℒ\left[f_s(\lambda )\right]}\right),$$(4)

where c_p is a linear scaling factor accounting for the brightness of the companion and M_{p,LS F} is the model of the planet convolved at the spectral resolution of the instrument. The terms $f_s\frac{ℒ\left[T(\lambda )M_{p,LSF}(\lambda )\right]}{ℒ\left[f_s(\lambda )\right]}$ comes from the leaking of the planet model continuum into the estimate of α (see Landman et al. 2024). This effectively removes the continuum of the planet model.

The total transmission can be estimated using our stellar master spectrum, f_s . We simply divide the stellar master spectrum by a model of the star: $$T(\lambda )=\frac{f_s(\lambda )}{M_s(\lambda )},$$(5)

where M_s (λ) is a PHOENIX model selected at the known effective temperature and surface gravity, and rotationally broadened (see Sect. 3.2).

For simplicity, we will drop the (λ) notation in the rest of the paper, but it is important to note that all the parameters used in this model depend on the wavelength. Our final model can be written as $$d=c_s\frac{ℒ[d]}{ℒ\left[f_s\right]}{f}_s+c_p\left(M_{p,LSF}T-\frac{f_sℒ\left[T{M}_{p,LSF}\right]}{ℒ\left[f_s\right]}\right)+\eta$$(6)

We can forward model our spectrum in order to jointly estimate the stellar and planetary contributions to the data. Equation (6) can be written in matrix form using the same conventions as in Wang et al. (2021): $$\left(\begin{array}{c}⋮\\ d_i\\ ⋮\end{array}\right)=\left(\begin{array}{cc}⋮& ⋮\\ \frac{ℒ[d]}{ℒ\left[f_s\right]}{f}_s& M_{p,LSF}T-\frac{ℒ\left[M_{p,LSF}T\right]}{ℒ\left[f_s\right]}{f}_s\\ ⋮& ⋮\end{array}\right)\left(\begin{array}{c}{c}_s\\ c_p\end{array}\right)+\eta ,$$(7)

where the left-hand side of the equation corresponds to a column vector with length equal to the number of spectral channels, N_λ. The first matrix in the right-hand side of the equations has an N_λ × 2 dimension. The last term in the right-hand side equations is just a vector of length 2 corresponding to the linear scaling coefficients c_s and c_p that we want to determine. Finally, we can write this model as $$d=M_{\psi }*c+\eta ,$$(8)

where ψ represents the planetary parameters that will determine the shape of the planetary spectrum model, M_{p,LS F} . The parameters include the atmospheric parameters T_eff, log g, [Fe/H] and C/O (see Sect. 5), but also the RV shift and the projected rotational velocity v sin i. For the projected rotational velocity, we used the fastRotBroad function from PyAstronomy (Czesla et al. 2019).

5 Atmospheric characterization

To characterize the atmosphere of AF Lep b, we used forward modeling analysis, consisting in using pre-computed grids of self-consistent models that can cover a range of parameters, including T_eff, log g, [Fe/H], and C/O, or more parameters (see Sect. 5.2). We use the ForMoSA Python package (Petrus et al. 2021)³, which we upgraded to work efficiently with high- spectral resolution data and in which we implemented the model described in Sect. 4. We present our results in Sects. 5.2 and 5.3.

5.1 ForMoSA: A forward modeling analysis tool

ForMoSA has already been extensively described and used in previous works to characterize exoplanets such as HIP 65426 b (Petrus et al. 2021), VHS 1256 AB b (Petrus et al. 2023, 2024), AB Pic b (Palma-Bifani et al. 2023) and AF Lep b (Palma-Bifani et al. 2024). These works were all based on low and medium resolution data.

ForMoSA relies on a nested sampling algorithm involving Bayesian inference (Skilling 2006). The nested sampling method was built to naturally estimate the marginal likelihood, using a multi-surface approach of space parameters exploration during the Bayesian inversion. The marginal likelihood z is defined as $$ɀ={\displaystyle \int \pi }(\theta )L(\theta )d\theta ,$$(9)

with L(θ) the likelihood function and π(θ) the prior distribution. The prior distribution and the likelihood function are defined by the users, given assumptions on the model. Calculating the value of z has the advantage of different model assumptions to be compared through the ratio of their evidence values. This ratio is known as the Bayes factor. Let z₁ and z₂ be the marginal likelihoods of two set of model assumptions. Then the Bayes factor is defined as $$B=\frac{{ɀ}_1}{{ɀ}_2}$$(10)

However, as ForMoSA returns the logarithm of the evidence (log z), we can make use of the logarithm of the Bayes factor. Two model assumptions can then be compared through the difference of the logarithm of the Bayes factor: $$\log B=\log {ɀ}_1-\log {ɀ}_2$$(11)

A positive value of log B denotes a statistical preference of the data to the first model assumption. More quantitatively, Bayes factors can be interpreted against the Jeffrey scale (see Table 2 of Benneke & Seager 2013), in which case a value of log B above 5 is considered a strong statistical preference for the first model assumption compared to the second model assumption.

For this work, we chose the following log-likelihood function: $${\chi }_0^2={\left(d-M_{\psi }\widehat{c}\right)}^T{\text{Σ}}_0^{-1}\left(d-M_{\psi }\widehat{c}\right),$$(12)

where Σ₀ is the covariance matrix of the data and ĉ is the linear least squares solution to Eq. (8) as follows: $$M_{\psi }^T{\text{Σ}}_0^{-1}{M}_{\psi }c=M_{\psi }^T{\text{Σ}}_0^{-1}d.$$(13)

For the covariance matrix, we assumed a simple model where the photon noise, mainly coming from the stellar contamination, is considered as the dominant source of noise in the data. Table 1 shows that the read-out noise is very low, which justifies this assumption. For this paper, we also assumed uncorrelated pixels, which means a diagonal covariance matrix. Practically, this matrix is calculated as the continuum of the extracted 1D spectrum at the location of the companion: $${\text{Σ}}_0=ℒ[d].$$(14)

To account for flux calibration offsets between spectral segments, which can impact the final results, we applied Eq. (6) to each order separately in order to properly rescale the stellar contamination to the level of the flux of the data. This also enables us to better estimate the continuum of the data as it can be problematic to evaluate the continuum with such offsets between segments.

Since telluric lines dominate the spectrum up to 1500 nm and from 1740 nm onward (see Fig. 1), we excluded the spectral segments corresponding to these parts of the spectrum from our analysis. This represents a total of nine spectral segments out of the 24 available. We thus end up using 15 spectral segments, covering a wavelength range approximately between 1500 and 1730 nm.

5.2 Atmospheric models

We ran ForMoSA with a model referred to as Exo-REM/Exo_k. This model is a result of using Exo-REM volume mixing ratio profiles in tandem with Exo_k to produce high-resolution spectra (Radcliffe et al., in prep). On one hand, Exo-REM is a radiative- convective model developed to simulate atmospheres of young giant exoplanets (Baudino et al. 2015; Charnay et al. 2018, 2021). The chemical abundances of each element is defined according to Lodders (2010). It implements disequilibrium chemistry. The flux is solved iteratively assuming 64 pressure levels over the grid, with a minimum pressure of 10⁻⁶ bar and a maximum pressure of 10² bar. This model considers sources of opacities from H₂-H₂ and H₂-H_e collision induced absorption, ro-vibrational bands from nine molecules (H₂O, CH₄ , CO, CO₂, NH₃ , TiO, VO and FeH) and resonant lines from Na and K.

However, the highest resolution Exo-REM models stand at R = 20000. Thus, we use Exo_k, a library constructed to handle radiative opacities from various sources for the subsequent computation emission spectra for 1D planetary atmospheres (Leconte 2021), to recompute the radiative transfer at a higher resolution for a fixed atmospheric structure. We use the cross section data listed in petitRADTRANS (Mollière et al. 2019) and the collision-induced absorption data from HiTRAN (Karman et al. 2019). The resulting model which is used in this work has as free parameters T_eff, log g, [Fe/H] and C/O. T_eff extends from 500 to 1000 K, log g from 3 to 5 dex, [Fe/H] from -0.5 to 1.0 dex, and C/O from 0.1 to 0.8. The models are initially computed from 1 to 5 µm at a spectral resolution of R = 1 000 000. For the analysis in ForMoSA, we downgrade the spectral resolution of the model to 200 000 but keep the bin sampling at 600 000. The Exo- REM/Exo_k model we use does not include the effects of clouds. This is discussed in Sect. 5.4.

To compare the model to the data in ForMoSA, we first apply the RV and υ sin i correction, then we downgrade the resolution of the model at the spectral resolution of our data, which is estimated at 140000. The resolution of the CRIRES⁺ spectrograph is constant to better than 2.5% in the H band (Dorn et al. 2023) and the value of 140 000 can be estimated from instrumental parameters. This value is also in perfect agreement with the value derived by Nortmann et al. (2025) using on-sky data obtained in very good observing conditions, where the stellar PSF has the same 2-pixel FWHM as the PSF of the HiRISE science fiber. The degradation of the resolution is performed using a Gaussian convolution adapted to the resolution at each wavelength.

For the present work, we used the nested sampling in ForMoSA with 500 living points. We considered two types of priors for the parameters: either uninformative uniform priors, or informative priors for T_eff, log g, [Fe/H] and C/O. In all cases, since we had no prior information on RV and υ sin i, we consistently applied uninformative priors for these parameters, that is u(-100, 100) for RV and u(0, 100) for v sin i. We derived our priors from the final atmospheric parameter values reported in Balmer et al. (2025): T_eff = 800 ± 50 K, log g = 3.7 ± 0.2 dex, [Fe/H] = 0.75 ± 0.25 dex, and C/O = 0.55 ± 0.1. Table 2 presents the priors we use in the forward modeling.

The results are presented in Table 3. In the present section, we present only our results on RV and υ sin i and discuss our results and their implications further in Sect. 7. From the results we can estimate the total RV of the system, which includes both the orbital velocity of the companion and the systemic velocity of the system. For the first night and the first offset of the second night we find RV and υ sin i values that agree well. For the second offset of the second night we find RV values slightly higher, but mostly υ sin i significantly higher than for the other data sets. We note that for the second offset we have a lower integration time for the science and an even lower integration time for the backgrounds, resulting in much noisier data (see Table 1). This could explain the discrepancy between the results of the second offsets and the results of the other datasets. Therefore, we decided to exclude this dataset for the rest of the analysis. We refer the reader to Appendix B for a more detailed analysis on the impact of the number of backgrounds to the derived parameters of the planet. We also note that observing conditions were different between the two nights (see Table 1), with a more variable weather on the second night.

To derive the final values adopted for RV and υ sin i, we combined the datasets of the first night and the first offset of the second night. To do so, we co-added the log-likelihoods associated with the two datasets on ForMoSA. The result of this combination is presented in Fig. 4. We estimate a final RV of ${31.61}_{-0.94}^{+0.96} \text{km} {\text{s}}^{-1}$ and a final υ sin i of ${12.38}_{-2.54}^{+2.30} \text{km} {\text{s}}^{-1}$. This value must be corrected from the systemic velocity, which is estimated to 21.1 ± 0.37 km s⁻¹ by Gaia Collaboration 2023. This value is perfectly consistent with ground-based high-resolution UVES data (Zúñiga-Fernández et al. 2021, 20.90 ± 1.11 km s⁻¹). Using the value from Gaia Collaboration (2023) and propagating our 100 m s⁻¹ uncertainty on the wavelength recalibration (see Sect. 3.3), we were able to finally infer a relative radial velocity between the companion and the star of ${10.51}_{-1.02}^{+1.03} \text{km} {\text{s}}^{-1}$. This value is consistent at 2 σ with the value estimated with the orbital solution derived from previous astrometric measurements (see Sect. 6). Another team of researchers (Hayoz et al. 2025) derived the radial velocity of the planet around the same epoch as our measurement using ERIS/SPIFFIER observations. Their value of 7.8 ± 1.7 km s⁻¹ confirms our measurement.

To assess the quality of the fit, we plot the best Exo- REM/Exo_k ForMoSA fit with the data in Fig. 3. We see that the data are highly dominated by the stellar contribution. We note that for this plot, the continuum of the planet model was removed from the planet model (see Eq. (4)). This is because the continuum of the planet model leaks into the estimate of α (see Eq. (2) and Landman et al. 2024). Using the estimate of c_p defined in Eq. (7), we can give an estimate of the total flux of the companion model we estimate (Table 4). In the H band, De Rosa et al. (2023) derives ∆H = 12.48 ± 0.12 mag. Our estimate is fully consistent with this value for the first night and the first offset of the second night. This confirms that, for these data, we are neither over- nor underestimating the planetary model. For the second offset of the second night, however, the flux of the planet model is slightly overestimated by ∼0.3 mag.

In Fig. 6, we focus on one of the best fitted segments and compare the data after subtraction of the stellar model with the planet model. In this wavelength band, the absorption lines are mainly coming from CH₄ . Some features of the planetary model can be identified in the raw residuals but become obvious when smoothing the data at the intrinsic spectral resolution of the model, which has been rotationally broadened by the fitted value of υ sin i.

5.3 Detection of molecules

Another approach to the analysis high-resolution data is to use CCF analysis; for example, to assess the presence of individual molecules. In this analysis, we cross-correlated template models with the data subtracted by the estimated stellar contribution (see Sect. 4): $$\widehat{d}=d-{\widehat{c}}_s\frac{ℒ[d]}{ℒ\left[f_s\right]}{f}_s,$$(15)

where ĉ_s is the first component of the solution ĉ to the linear least square Eq. (13). We calculate the CCF as $$\mathrm{CCF}(RV)={\displaystyle \sum _{{\lambda }_i}\widehat{d}}\left({\lambda }_i\right)m_{RV}\left({\lambda }_i\right),$$(16)

where m_RV is the template model Doppler-shifted at the radial velocity RV and broadened at the υ sin i given by our results (see Table 3). For this analysis we compute the CCF using a grid of RV ranging from -1000 to + 1000 km s⁻¹ in steps of 0.5 km s⁻¹. We perform the analysis both for best Exo-REM/Exo_k model inferred from the ForMoSA analysis with priors on all the bulk parameters and for model templates of individual molecules. We consider only H₂O and CH₄ as they are the species expected to be dominating in the H band at the T_eff of AFLep b. The templates of individual molecules have been generated using Exo-REM volume mixing ratio profiles in tandem with Exo_k.

The results are presented in Fig. 5 for the night 2023-11-20. For the second night, the results are presented in Appendix C and yield similar conclusions to the ones discussed below. To improve visibility, we plot the RV between -300 and +300 km s⁻¹. The S/N is estimated by normalizing the CCF by its standard deviation computed over two windows 100 km s⁻¹ away from peak of the CCF (30 km s⁻¹), which in this case is approximately between −70 and −1000 km s⁻¹ and between +130 and +1000 km s⁻¹. We overplot the autocorrelation function (ACF) of the model, which we shift to the estimated RV and normalize at the peak S/N of the cross-correlation function.

To accurately estimate the RV associated with the CCF, we fit a Gaussian function to the CCF. We find this method of estimating the RV to be more robust compared to estimating the RV with the maximum of the CCF. With the full model, we estimate a S/N of 6.6 at RV = 31.0 km s⁻¹. For H₂O, we find a detection with a S/N of 6.6 at 31.4 km s⁻¹, consistent with the full model. And finally for CH₄, the detection is at a S/N of 4.4 with a RV of 31.5 km s⁻¹ , again consistent with the full model. We note that these values are within the error bars of the value inferred from our ForMoSA analysis of the first night with the priors on all the bulk parameters. Although the CH₄ CCF is noisier than for the ones for the full model and for H₂O, the location of the peak in RV makes us confident that the detection of CH₄ in the atmosphere of AF Lep b is real, confirming the recent findings of Balmer et al. (2025). This also supports the possibility of inferring a C/O ratio in the H band.

To enhance the confidence in the detections, we overplot in Fig. 5 the CCF between the model and the other three reference fibers of the instrument, which contains signal from the star speckles (see Vigan et al. 2024, for details on these fibers). The CCF with the reference fibers shows similar features to the planet’s CCF. We interpret this by the fact that all fibers see some common signal coming from thermal background noise, which is expected to be at the same level for all fibers, and from star speckles and systematic effects, which are imperfectly removed from the science signal.

Fig. 3 Best-fit model for the 2023-11-23 data with the Exo-REM/Exo_k model. Top panel shows the data, the full model (magenta), the stellar model (blue) and the planet model (red). The last two are scaled to their respective amplitude in the total signal. The stellar component in the data completely dominates the planetary component, so the full model (magenta) is mostly hidden by the star model (blue). Bottom panel shows the residuals, i.e., the data subtracted by the full model. For this panel, the scale is the standard deviation of the residuals. The top of the lower part of the figure depicts the distribution of the residuals, which is well centered around 0.

Fig. 4 Posterior distributions of RV and υ sin i for the two nights combined with the Exo-REM/Exo_k model. These results are the final adopted values for these parameters.

5.4 Considering whether our data are sensitive to clouds

At the low temperature of AF Lep b, condensate clouds should sink below the photosphere, reducing their effect in the emission spectrum (Lodders & Fegley 2006). However, Zhang et al. (2023) noted the possible presence of silicate clouds in the atmosphere of AF Lep b (see Fig. 15 of Zhang et al. 2023). As illustrated in Fig. 6 of Xuan et al. (2022), high-resolution data are sensitive to lower pressures compared to low-resolution data. This is because high-resolution data resolve the cores of absorption lines (from H₂O and CH₄ molecules in this case) which, compared to the continuum, comes from regions of high molecular opacity. The optical depth in a layer of the atmosphere is proportional to the integral of opacity and abundance, integrated over the path length. Since the opacity at line cores is high, the path length must be small, so these regions come from lower pressures. As a consequence, our very high-resolution data may be sensitive to pressures lower than the silicate clouds base pressure inferred in Zhang et al. (2023).

To verify this hypothesis, we plot the best P-T profile inferred from our results of the combination of the two nights in Fig. 7. We can derive the condensation curve of MgSiO₃ from Eq. (20) of Visscher et al. (2010), which depends on the metallicity. To plot the 1σ condensation curve of MgSiO₃, we use the 1σ uncertainties on metallicity derived in Balmer et al. (2025) ([Fe/H] = 0.75 ± 0.25 dex). We also estimate the pressure levels to which our data are most sensitive by calculating the brightness temperature of our best model at each wavelength and match it to the pressure in the PT profile. The relatively narrow region to which our data are most sensitive is a result of the relatively small wavelength range covered by our data (1.43–1.77 µm). As we see in Fig. 7, our data are sensitive to layers in the atmosphere significantly above the silicate cloud base pressure level, justifying our assumption that our data should not be sensitive to clouds and that we can safely use cloudless models.

6 Orbit fitting

We used Orvara (Brandt et al. 2021b) to constrain the orbit of AF Lep b. We ran the parallel-tempering Markov Chain Monte Carlo (MCMC) sampler (ptemcee) of Orvara (Foreman Mackey et al. 2013; Vousden et al. 2016) with 10 temperatures, 100 walkers, 500 000 steps per walker. The chains were saved every 50 steps and in each chain we discarded the first 3000 saved steps as burn-in. We fit for the mass of the star (M_⋆), the mass of the companion (M_p), the orbit semi-major axis (a), the inclination angle (i), eccentricity (e) and argument of periastron (ω) in parametrized forms ($\sqrt{e}$ sinω and $\sqrt{e}$ cos ω), the longitude of the ascending node (Ω), and the mean longitude at reference epoch J2010.0 (λ).

In Sect. 6.1, we describe the datasets used in the analysis. To assess the contribution of our RV measurements to the orbit of AF Lep b, we first analyze the orbit excluding our RV measurement (Sect. 6.2) and then we compare with the results obtained using our HiRISE measurement (Sect. 6.3).

Fig. 5 Cross-correlation functions for the first night with the Exo-REM/Exo_k model. In grey, the cross-correlation between the signal of the 3 reference fibers and the best model. In blue, the auto-correlation between the best model and itself. The model has been broadened at the v sin i given by the fit.

Fig. 6 Comparison between the 2023-11-23 data after subtraction of the stellar contribution and the best fit Exo-REM/Exo_k model in a segment of order dominated by CH4 absorption lines. The top panel shows the raw data subtracted by the estimated stellar contribution at the full spectral resolution, while in the bottom panel, the data has been smoothed to a spectral resolution of ∼30 000 that takes into account the decrease of resolution induced by v sin i.

6.1 Datasets

In our orbit analysis we used relative astrometry spanning a baseline of 11 years. This includes the SPHERE astrometry (Mesa et al. 2023; De Rosa et al. 2023), the Keck/NIRC2 astrometry (Franson et al. 2023), the archival VLT/NaCo astrometry (Bonse et al. 2024), and the recent VLTI/GRAVITY observations (Balmer et al. 2025). We also included 20 RV measurements of the star obtained with Keck/HIRES (Butler et al. 2017), spanning a baseline of 12 years.

A recent analysis of these datasets lead by Balmer et al. (2025) measured a mass of ${3.68}_{-0.48}^{+0.47}{M}_{\text{Jup }}$ for the companion, a semi-major axis of ${8.98}_{-0.08}^{+0.15}$ au and a quasi-circular orbit $\left(\text{e}={0.013}_{-0.010}^{+0.024}\right)$. It also concluded in a spin-orbit alignment with an inclination angle $i={57.12}_{-0.71}^{+0.67}°$. This value is consistent with the stellar inclination angle $i_{\star }={54}_{-9}^{+11}°$ derived in Franson et al. (2023).

Fig. 7 P-T profile with 1σ and 2σ uncertainties derived from the results of combining the 2 nights. In green we show the 1σ condensation curve of MgSiO3 . The shaded blue region corresponds to the pressure levels to which our data are most sensitive.

6.2 Orbit analysis

These datasets provide good constraints on the main orbital parameters, but they leave some ambiguity regarding the direction of motion of the companion in the line of sight. Two families of orbit co-exist depending on the sign of the RV of the companion relative to the host star at epoch J2024.0. This was first demonstrated by Zhang et al. (2023). This leaves two degenerate solutions for the mean longitude of the star at epoch J2010.0 (λ) and the longitude of ascending node (Ω) (see Fig. D.1).

From the orbit, we can derive the RV of the star (RV_⋆) at any reference epoch t_ref using the following equation: $${\text{RV}}_{\star ,ref}=K_{\star }\left(e\cos (\omega )+\cos \left(v_{ref}+\omega \right)\right)$$(17) $$K_{\star }={\left(\frac{2\pi G}{P}\right)}^{1/3}\frac{M_p}{{\left(M_p+M_{\star }\right)}^{2/3}}\frac{sin(i)}{\sqrt{1-e^2}},$$(18)

where v_ref is the true anomaly of the companion at epoch t_ref. From this equation, we can infer the RV of the companion (RV_p) with the following relation: $$M_p{\text{RV}}_p=-M_{\star }{\text{RV}}_{\star }.$$(19)

This is coming from the conservation of angular momentum of the system and the fact that the companion moves in opposite direction of the star. The RV of the companion itself is $${\text{RV}}_{p,ref}=K_p(e\cos (\omega )+\cos \left(v_{ref}+\omega \right),$$(20) $$K_p=-{\left(\frac{2\pi G}{P}\right)}^{1/3}\frac{M_{\star }}{{\left(M_p+M_{\star }\right)}^{2/3}}\frac{sin(i)}{\sqrt{1-e^2}}.$$(21)

Finally, the relative RV between the planet and the host star can be expressed as $$\begin{array}{c}{\text{ΔRV}}_{ref}=-{\left(\frac{2\pi G\left(M_{\star }+M_p\right)}{P}\right)}^{1/3}\frac{\sin (i)}{\sqrt{1-e^2}}\\ \times \left(e\cos (\omega )+\cos \left(v_{ref}+\omega \right)\right).\end{array}$$(22)

At epoch J2023.88, corresponding to 2023-11-20, we find, using the dataset described in Sect. 6.1 and Eq. (22), one population of orbit with $\text{RV}={8.76}_{-0.35}^{+0.22}$ km s⁻¹ and another population of orbit with $\text{RV}={8.78}_{-0.21}^{+0.33}$ km s⁻¹. The first population is consistent within 2σ with our HiRISE measurement of ${10.51}_{-1.01}^{+1.03}$km s⁻¹ for AFL ep b (see Sect. 5.2).

Fig. 8 Astrometric orbit of AF Lep b with 1000 orbits randomly drawn from the posterior. The left panel shows the solutions obtained without the HiRISE RV measurement and the right panel shows the solutions including our measurement. The inset plot in each panel shows the distributions of the relative RV between the planet and the star. As illustrated on the plot, the determination of the RV of the planet offers useful information on its phase.

Fig. 9 Predictions of the relative RV between the companion and the star obtained with 1000 orbits randomly drawn from the posterior. As in Fig. 8, the left panel shows the solutions without the HiRISE RV measurement and the right panel shows the solutions including our measurement. The bottom left of the right panel, shows a zoom around the HiRISE measurement.

6.3 HiRISE RV measurement

In practice, Orvara only uses the relative RV measurements between the secondary and the primary. However, the ratio of the RVs between the primary and the secondary is scaled by a factor of $\frac{M_p}{M_{\star }}$, which in most cases is significantly lower than 1. It is the case here as $\frac{M_p}{M_{\star }}={0.00289}_{-0.00036}^{+0.00036}$ (see Table 5). As a result, the RV of the star is entirely dominated by the uncertainty in the RV of the companion. We thus directly included the ${10.51}_{-1.01}^{1.03}$ km s⁻¹ RV measurement of the companion on Orvara.

The results of adding our RV measurement are presented in Fig. 8 (right) and Table 5. The astrometric orbit with our RV measurement appears visually similar to the one without it. But, with the RV measurement included, we rule out an entire family of orbital solutions. The bottom left of each panel in Fig. 8 presents the distribution of the relative RV between the companion and the star. We clearly see two widely separated populations of orbits when our RV measurement is not included (left panel), which vanish when we include our measurement (right panel). Also, our measurement sets the value of the argument of periastron ω, albeit still with large uncertainties, the ascending node Ω and the mean longitude at reference epoch J2010.0 λ (see Fig. 10). In addition, the very precise astrometric orbit, together with the RV information, offers precise information on the phase of the planet at a given epoch. This will be of particular interest for the next generation of instruments probing the reflected light of exoplanets.

Figure 9 presents the prediction of the relative RV between the companion and the star as a function of time. Similarly, our RV measurement resolves the ambiguity in the sign of the RV. However, our measurement does not provide additional constraints on its exact value. This is illustrated in the bottom left inset of the right panel, which presents a zoom around our HiRISE measurement.

Fig. 10 Posterior distribution for the ascending node (Ω), the argument of periastron (ω) and the mean longitude at reference epoch J2010.0 (λ). Grey: without our RV measurement. Red: with our RV measurement.

7 Discussion

7.1 Effective temperature and surface gravity

Previous studies have suggested a lower T_eff, on the order of 800 K (Palma-Bifani et al. 2024; Zhang et al. 2023). This parameter, as well as log g, impact mostly the continuum of the spectra as well as the depth of molecular lines. Our HiRISE data tend to favour low log g values and high T_eff values, compared to what is expected (see Table 3.) We also find strong correlations between T_eff and log g in our results (see Figs. E.1 and E.2), indicative that a more consistent lower T_eff would also decrease the estimated value of log g. We note that low values for log g are physically inconsistent with the estimated radius of AF Lep b (1.3 ± 0.15 R_Jup, Balmer et al. 2025) as the value of log g strongly affects the radius estimated from Newton’s law (see Table 3). For a log g of 3.22 dex, corresponding to the value obtained for the first night with uninformative priors on all parameters, the self consistent radius using Newton’s law would be 2.37 R_Jup, which is highly inconsistent with the results of Zhang et al. (2023), Palma-Bifani et al. (2024) and Balmer et al. (2025). This demonstrates the importance of using an accurate value on log g when fitting a spectra. To address this issue, we used Gaussian priors on T_eff and log g (T_eff ~ N(800,50), log g ~ N(3.7,0.2)). However, even with these Gaussian priors, our results tend to favour T_eff values slightly higher and log g slightly lower than the expected values.

As a result, it is difficult to constrain these parameters using only our HiRISE data. These two parameters are generally best estimated with high S/N low-resolution data over a large wavelength range that offers accurate information on the depth of the spectral lines as well as on the continuum; however, this is not the case here since we have removed the continuum (see Sect. 4). It would therefore be interesting to combine low-resolution data with HiRISE data in order to have a robust estimate of these parameters without the need to apply strong priors to these parameters.

7.2 Carbon-to-oxygen ratio and metallicity

The literature reports enriched metallicities for AF Lep b, higher than 0.4 dex for Palma-Bifani et al. (2024) and around 0.75 dex for Balmer et al. (2025), but based on low-resolution data. When we leave all parameters free, we find values of [Fe/H] inconsistent with these results at 3 σ for the first night and 2 σ for the second night. The difficulty to constrain [Fe/H] with high- resolution data alone has been previously noted by Landman et al. (2024). This is mainly the results of the low S/N of our data. Indeed, Xuan et al. (2024) and Costes et al. (2024b) show that at high S/N, high-resolution data can give very accurate [Fe/H] estimates. Again, combining HiRISE data with other datasets at lower resolution with high S/N and over a larger wavelength range would certainly help provide better constraints on this parameter.

For the C/O ratio, the detection of CH₄ in the atmosphere of AFLepb (see Sect. 5.3), as also supported by Balmer et al. (2025), reinforces the possibility of deriving a C/O ratio in the H band. The C/O ratio estimated using HiRISE data alone for the second night appears inconsistent with the near-solar values reported in the literature when all parameters are free. However, this parameter tends to vary when we apply a prior on [Fe/H] (see Table 3). For both nights, the C/O ratio converges to near-solar or slightly super-solar values without the need for a constrained prior on this parameter as depicted in Fig. 11. This demonstrates that combining high-resolution with low-resolution data may help break degeneracies that can exist between these parameters.

The C/O ratio is often assumed as a tracer of the formation location of the planet in the protoplanetary disk (Öberg et al. 2011; Mordasini et al. 2016; Molyarova et al. 2017; Madhusudhan et al. 2017; Booth et al. 2017). This is because gaseous H₂O, CO₂ and CO molecules condensate into solid icy grains at various snowlines (Öberg et al. 2011), therefore affecting the C/O ratio in the gas and solid phases. However, caution should be taken when interpreting the C/O ratio as it is particularly sensitive to planet formation and evolution assumptions (see e.g., Mollière et al. 2022; Hoch et al. 2023). As shown in Balmer et al. (2025), the C/O ratio and [Fe/H] for AF Lep b derived from GRAVITY data could either lead to a core formed beyond the CO iceline, which would have migrated inwards, or to an in situ formation, depending on the assumptions made on the chemical evolution of the disk.

Fig. 11 First offset of the second night posterior distribution of [Fe/H] and C/O for the Exo-REM/Exo_k model with priors on Teff and log g. Three cases are considered: no prior (blue), Prior on [Fe/H] (orange), Priors on [Fe/H] and C/O (red).

7.3 Orbit, radial velocity, and projected rotational velocity

Previous orbital analyses of AF Lep b have been conducted by De Rosa et al. (2023), Franson et al. (2023), and Mesa et al. (2023), who independently concluded in very different mass estimates: ${4.3}_{-1.2}^{+2.9}$ M_Jup, ${3.2}_{-0.6}^{+0.7}$ M_Jup, and ${5.24}_{-0.10}^{+0.09}$ M_Jup, respectively. These discrepencies come from the fact that their relative astrometry was measured at different epochs and over different baselines. By combining all the astrometry measurements used in these analyses, Zhang et al. (2023) gave an update on the orbital parameters of AF Lep b, but they were still unable to capture the circular nature of the planet’s orbit because of the short time-baseline covered by the relative astrometry data. The addition of the recovered archival VLT/NaCo astrometry of 2011 relieves this issue and confirms the quasi-circular nature of the planet’s orbit (Bonse et al. 2024). This is because the NaCo astrometric data point is almost opposite to the more recent astrometry data points with respect to the star (Fig. 8). However, for the same reason, it does not help to better constrain the inclination of the orbit. Balmer et al. (2025) has further confirmed the circular nature of the orbit with the addition of the GRAVITY measurements, making a significant improvement in the constraints on the planet’s orbital parameters. With this dataset, the amplitudes of RVs of the star and the companion are well constrained, but still leave some ambiguity regarding the sign of the RV. Our HiRISE RV measurement on the companion itself allows us to resolve this ambiguity.

The estimated radial velocity and projected rotational velocity do not vary when we change the priors configuration, as shown in Table 3 (see also Figs. E.3 and E.4). This comforts us in the derived value for these parameters. The measured projected rotational velocity gives us an idea of the period of rotation of the companion. The rotation rates of exoplanets are an important parameter since they determine the magnitude of their Coriolis force, which is a key parameter in understanding the climate of exoplanets. It is also linked to the formation of exoplanets since it is related to the accretion of angular momentum during the formation phase in the proto-planetary disk (Lissauer & Kary 1991; Dones & Tremaine 1993; Johansen & Lacerda 2010; Batygin 2018). From the combination of the orbit fitting and forward modeling, we infer a rotational velocity of $\frac{v\sin i}{sin(i)}={14.63}_{-2.87}^{+2.63}$ km s⁻¹ for AFLepb, assuming the orbit and rotation of the planet are coplanar $\left(i={57.78}_{-0.58}^{+0.57°}\right)$. From this rotational velocity, we can compute the rotation period with $$P_{\text{rot}}=\frac{2\pi R}{v_{\text{rot}}},$$(23)

with R the radius. We find a period of rotation of ${10.86}_{-2.08}^{+3.00}$ hours, given a radius of 1.3 ± 0.15 R_Jup (Balmer et al. 2025). The rotation period of Jupiter is (~ 10 hours), corresponding to a rotational velocity of 12.20 km s⁻¹. This result is in line with the tentative trend in spin velocity with planet mass (Hughes 2003; Wang et al. 2021; Morris et al. 2024; Xuan et al. 2024; Hsu et al. 2024). However, caution should be made when inferring the rotation period from v sin i and sin i, or equivalently, when inferring sin i from the rotation period and v sin i. Indeed, because v sin i and sin i are correlated, combining these estimates is somewhat more complex (Masuda & Winn 2020; Bryan et al. 2021; Bowler et al. 2023; Zhang et al. 2024).

7.4 Limitations and possible improvements

Despite the analysis conducted above, it is important to keep in mind that high-resolution models still perform poorly in predicting observations. Discrepancies such as continuum mismatches (Lim et al. 2023; Jahandar et al. 2024), unidentified spectral features in observations (Jahandar et al. 2024), and shifts in spectral features (Tannock et al. 2022; Jahandar et al. 2024) can strongly bias parameter estimations. Systematics between models (Ravet et al., in prep.; Petrus et al. 2024) also play a significant role. These issues could be mitigated by fitting only the spectral lines or parts of the spectrum we can trust (Jahandar et al. 2024; Petrus et al. 2024). Additionally, the effects of systematics could be incorporated into a refined modeling of the covariance matrix. At present, our covariance matrix modeling assumes uncorrelated noise. However, correlated residuals across adjacent pixels may arise due to rotational broadening, imperfections in the template model, or systematics. This issue is generally addressed through the use of Gaussian processes (Czekala et al. 2015; Kawahara et al. 2022; Iyer et al. 2023; de Regt et al. 2024).

A second limitation lies in the extraction of the signal itself. Currently, we extract the spectrum by simply summing the signal within a window centered around the trace, extending over 6 pixels (see Sect. 3). This approach has been shown to be suboptimal and could be improved by applying weighting to the sum based on the shape of the fitted trace (i.e., using an “optimal extraction” algorithm; Horne 1986; Marsh 1989; Zechmeister et al. 2014; Piskunov et al. 2021; Holmberg & Madhusudhan 2022), typically offering a 70% improvement in effective exposure time. This would probably result in more robust estimates (lower error bars) of the parameters of the planet.

A final limitation lies in the observation strategy. The current wavelength calibration correction is likely not fully optimal. Future observations could leverage M-dwarf or M-giant stars, which have deep and stable spectral lines, to achieve robust wavelength calibration across all order segments (see Sect. 3.3).

8 Conclusions

We obtained new H-band high-resolution (R ≈ 140000, 1.4−1.8 µm) spectroscopic data of AFLepb with the new VLT/HiRISE instrument. These observations result in a detection of AF Lep b with the first direct measurement of its projected rotational velocity and RV. We also detected methane absorption in the planet’s atmosphere at a 3σ significance level. This detection, a major finding of this paper, is consistent with the observations of Balmer et al. (2025). It also supports our derivation of the C/O ratio, as the only significant carbon-bearing species in the H-band is CH₄ .

We give new constraints on the orbit of AF Lep b, ruling out a whole family of orbital solutions. Combining our results with the low-resolution results from Balmer et al. (2025), we can give an update on the atmospheric and orbital parameters of the planet: T_eff = 800 ± 50 K, log g = 3.7 ± 0.2 dex, R = 1.3 ± 0.15 R_Jup, [Fe/H] = 0.75 ± 0.25 dex, C/O = 0.55 ± 0.10, and $\text{RV}={8.76}_{-0.35}^{+0.22}$ km s⁻¹ at epoch J2023.88 (2023-11-20).

Our data uniquely provide an estimation of the planet’s projected rotational velocity, $v\sin i={12.58}_{-2.54}^{+2.30}$ km s⁻¹. We predict the rotation period of AF Lep b to be 10.9 ± 3 hours, assuming that the orbit and rotation of the planet are coplanar. This issue could be resolved by acquiring time-resolved high-resolution photometric measurements of AF Lep b to measure variations the in brightness of its atmosphere. This was done for the first time on the 4 M_Jup directly imaged companion of 2M1207b (Zhou et al. 2016), but also for VHS 1256b (Zhou et al. 2020). However, for high-contrast planets, variability monitoring from the ground is limited by varying speckle noise (Biller et al. 2021; Wang et al. 2022). From space, the stability and high-contrast imaging capabilities of the James Webb Space Telescope could be promising for the purposes of variability monitoring on high- contrast planets such as AF Lep b. Franson et al. (2024) presents a detection along with the first attempt at variability monitoring of AF Lep b with the James Webb Space Telescope, but shows no evidence of variability. However, the authors used five epochs of observations spread over two years. Given the v sin i estimated for the planet, it would be interesting to obtain photometric data points distributed over a few hours. Still, the value derived for the projected rotational velocity of AF Lep b, in line with the known trend in spin velocity with planet mass, suggests that the planetary orbit and rotation axes should be aligned as is the case for Jupiter, the main gas giant in our system.

With the characterization of AF Lep b with HiRISE, we are getting closer to characterizing Jupiter-mass companions in direct imaging at high spectral resolution. HiRISE demonstrates a promising instrument for the characterization of such faint, low-mass companions; moreover, it excels at low-separation detection in direct imaging standards. This is particularly interesting since Santos et al. (2017) and Schlaufman (2018) have reported the existence of two populations of exoplanets which split at 4 M_Jup. Planets with masses above 4 M_Jup are found to be orbiting less metal-rich stars, compared to planets with masses below 4 M_Jup, where the host stars are more metal-rich. The latter is associated with the accretion formation pathway. Similarly, Hoch et al. (2023) determined there were two populations of companion when looking at C/O ratio and mass. Interestingly, they also concluded that the split is around 4 M_Jup. Hopefully, the ESA/Gaia DR4 results expected will provide additional Jupiter- mass planetary candidates in the coming years, which could help confirm such a trend.

The orbit of AF Lep b is now well constrained, both in astrometry and in RV. Balmer et al. (2025) proposed that astrometric monitoring of AF Lep b could detect deviations from the well-established astrometry, which might be caused by an inner planet or a massive exomoon in the system. This method has proven to be effective, even allowing constraints to be placed on the orbit and mass of the inner planet in the β Pictoris system (Lacour et al. 2021). Another way to do this would be through the acquisition of a high S/N spectrum of the star to study the variations of RV of the star. The instrument VLT/ESPRESSO would be a good candidate in such a study (Pepe et al. 2021).

Compared to transiting planets, directly imaged planets orbit farther away, which results in larger Hill spheres and makes these planets more favorable for moon formation and retention. In this context, it has been proposed that studying RV variations could be a promising way to detect exomoons around directly imaged planets (Vanderburg et al. 2018; Vanderburg & Rodriguez 2021; Ruffio et al. 2023). HiRISE, as well as future generations of instruments, could be promising, not necessarily for carrying out detections, but for setting initial constraints on the parameters of binary planets and exomoons around AF Lep b. Altogether, it is clear that AF Lep b will remain a benchmark for exoplanet characterization, providing valuable insights into planet formation and evolution both within and beyond our Solar System.

Acknowledgements

This project has received funding from Agence Nationale de la Recherche (ANR) under grant ANR-23-CE31-0006-01 (MIRAGES). The HiRISE instrument has been developed with funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme, grant agreements No. 757561 (HiRISE) and 678777 (ICARUS), from the Commission Spécialisée Astronomie-Astrophysique (CSAA) of CNRS/INSU, and from the Action Spécifique Haute Résolution Angulaire (ASHRA) of CNRS/INSU co-funded by CNES, and from Région Provence-Alpes-Côte d’Azur under grant agreement 2014-0276 (ASOREX). This research has made use of computing facilities operated by CeSAM data center at LAM, Marseille, France. S.P. is supported by the ANID FONDECYT postdoctoral program No. 3240145. N.G. acknowledges funding by the European Union (ERC, ESCAPE, projet No 101044152). Views and opinions expressed are however those of the author(s) only and do not reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

Appendix A Detector traces

Figure A.1 presents the tracing of the companion data for a single order. The top panel depicts a 2D image of the data when the fiber is centered at the companion. In the 2D image of the companion, we can observe the residual MACAO internal source fiber signal around the y-pixel number 90. Secondary signals are also visible around y-pixel number 115 and y-pixel number 85, originating from 2 reference fibers. The signal from the last reference fiber is too faint to be seen on this image. The bottom panel shows the shape of the signal between y-pixels 90 and 110 at different wavelengths (x-pixels 340,480 and 620). The dashed line is the profile of the companion data, whereas the solid line is the profile pf the star data. We can see that the profiles of the companion and star star are very similar across different wavelengths. Extracting the signal over 6 pixels is a reasonable choice to capture the entire signal at each wavelength without adding too much noise. By analyzing the trace at different positions on the x-axis, we see that the profiles of the companion and star data are globally very consistent.

Fig. A.1 Trace of the CRIRES+ detector for a single order. Top panel: 2D detector image of the companion. The hatched area corresponds to the noise estimation region. The spectrum of the companion is located at the center of the image, around the y-pixel number 100. Around y-pixel 90 is the residual MACAO internal source fiber signal after background subtraction. Around y-pixel number 115 and y-pixel number 85, we can see the signal of 2 of the 3 reference fibers. The bad pixels are depicted in white. Bottom panel: normalized profile between y-pixels 90 and 110 of the signal at x-pixels 340, 480 and 620 for the star (solid lines) and the companion (dashed lines). The extracted signal area is represented by the vertical red dashed lines.

Appendix B Impact of S/N ratio on RV and v sin i estimations

In Sect. 5.2, we decided to exclude the data of the second offset of the second night from the analysis, mainly due to the lower number of available backgrounds for this second offset. To verify this hypothesis, we ran ForMoSA on the data of the first night with an increasing numbers of backgrounds used in the data reduction, going from one to five backgrounds. The backgrounds are mainly used to subtract the MACAO guide fiber leakage (Vigan et al. 2024). The removal of this signal leaves some residual noise which decreases as we increase the number of backgrounds used. The results of this analysis are shown in Fig. B.1. We see that the more backgrounds we use to reduce the data, the more constrained the final RV and v sin i are. This effect is particularly strong when using a small number of backgrounds (less than three), where some distributions of RV and v sin i are bi- modal, and diminishes as the number of backgrounds increases. This is expected since we theoretically expect a reduction in the S/N of the reduced data of $\sqrt{N}$ with N being the number of backgrounds used, up to a certain point. The differences between the results with four and five backgrounds become very small, which indicates that we are probably reaching a fundamental limit in our data, such as photon noise or a mixture of photon noise and instrumental systematics. We recall that the photon noise in our data comes mainly from the stellar contribution at the location of the companion, which is the dominant source in our data.

Fig. B.1 First night posterior distribution of RV and v sin i for different number of backgrounds used in the data reduction. The number of backgrounds impacts the final S/N of the reduced data, which has a direct impact on the valued derived for the RV and v sin i of the companion.

Appendix C CCF analysis for the second night

Figures C.1 and C.2 present the cross-correlation analysis for the first and second offset of the second night. Compared to the first night, the detections have lower S/N, confirming the better data quality of the first night. The estimated RV values are 31.9 km s⁻¹, 31.7 km s⁻¹ and 31.8 km s⁻¹ for the full model, H₂O and CH₄ with the first offset. These values are consistent with the values estimated by ForMoSA (see Table 3). For the second offset, we estimate RVs of 33.1 km s⁻¹ with the full model and 34.0 km s⁻¹ with H₂O. The CCF with CF₄ does not show a convincing detection of CH₄, albeit it exhibits a small peak in the CCF around 34 km s⁻¹.

Fig. C.1 Cross-correlation functions for the first offset of the second night.

Fig. C.2 Cross-correlation functions for the second offset of the second night.

Appendix D Orvara results

Figure D.1 presents the posterior distribution obtained with Orvara. 2 cases are depicted: without including our RV measurement (grey) and with our RV measurement (red). The inclusion of our RV measurement does not impact the estimation of the orbital parameters except for the argument of periastron ω, the mean longitude at epoch 2010 λ, and the longitude of ascending node Ω. The latter two present bimodal distributions when we do not include our radial velocity measurement.

Fig. D.1 Posterior distribution for the orbit of AF Lep b with the astrometry from SPHERE, KECK/NIRC2 and GRAVITY (grey) and with our RV points included (red).

Appendix E ForMoSA results

This section presents some of the results obtained with ForMoSA. Figure E.1 presents the posterior distribution for the first night without any prior. Figure E.3 compares the posterior distributions for the first night with fixed T_eff and log g under the different configurations of priors adopted (See Table 3). Similarly, Figs. E.2 and E.4 present the posterior distribution for the first offset of the second night under the different configurations adopted for the priors.

Fig. E.1 First night posterior distribution of atmospheric parameters for the Exo-REM/Exo_k model in the case where all parameters are free

Fig. E.2 First offset second night posterior distribution of atmospheric parameters for the Exo-REM/Exo_k model in the case where all parameters are free

Fig. E.3 First night posterior distribution of atmospheric parameters for the Exo-REM/Exo_k model with priors on Teff and log g. Three cases are considered : No prior (blue), Prior on [Fe/H] (orange), Prior on [Fe/H] and C/O (red).

Fig. E.4 First offset of the second night posterior distributions of atmospheric parameters for the Exo-REM/Exo_k model with priors on Teff and log g. Three cases are considered: No prior (blue), Prior on [Fe/H] (orange), Prior on [Fe/H] and C/O (red).

References

All Tables

All Figures

Fig. 1 Data for the star and companion obtained on 2023-11-23 and extracted using our custom pipeline. The read-out noise of the detector is also plotted to demonstrate that the data obtained on the companion has a mean S/N of ∼25 per spectral channel at the center of the H band. The effect of the telluric lines is clearly visible at the start and end of the band, but also close to the center where many telluric CO bandheads are visible. We remind that data for the companion refers to the spectrum at the location of the companion which contains primarily stellar flux.

In the text

Fig. 2 Corrections to the wavelength solution for each segments of orders based on the analysis of the telluric lines for the AF Lep data acquired on 2023-11-23. For the two segments centered at 1556 and 1623 nm (light red), there are not enough telluric lines to compute an accurate correction and the default wavelength solution is therefore adopted. The right axis of the main plot shows the RV shift corresponding to the level of correction of the left axis, computed at 1600 nm. The top-left inset shows the effect of the correction at 1666.5 nm (grey vertical line in the main plot), which is a segment that requires one of the largest corrections. The telluric model obtained from SkyCalc is overplotted.

In the text

Fig. 3 Best-fit model for the 2023-11-23 data with the Exo-REM/Exo_k model. Top panel shows the data, the full model (magenta), the stellar model (blue) and the planet model (red). The last two are scaled to their respective amplitude in the total signal. The stellar component in the data completely dominates the planetary component, so the full model (magenta) is mostly hidden by the star model (blue). Bottom panel shows the residuals, i.e., the data subtracted by the full model. For this panel, the scale is the standard deviation of the residuals. The top of the lower part of the figure depicts the distribution of the residuals, which is well centered around 0.

In the text

	Fig. 4 Posterior distributions of RV and υ sin i for the two nights combined with the Exo-REM/Exo_k model. These results are the final adopted values for these parameters.
In the text

	Fig. 5 Cross-correlation functions for the first night with the Exo-REM/Exo_k model. In grey, the cross-correlation between the signal of the 3 reference fibers and the best model. In blue, the auto-correlation between the best model and itself. The model has been broadened at the v sin i given by the fit.
In the text

Fig. 6 Comparison between the 2023-11-23 data after subtraction of the stellar contribution and the best fit Exo-REM/Exo_k model in a segment of order dominated by CH4 absorption lines. The top panel shows the raw data subtracted by the estimated stellar contribution at the full spectral resolution, while in the bottom panel, the data has been smoothed to a spectral resolution of ∼30 000 that takes into account the decrease of resolution induced by v sin i.

In the text

	Fig. 7 P-T profile with 1σ and 2σ uncertainties derived from the results of combining the 2 nights. In green we show the 1σ condensation curve of MgSiO3 . The shaded blue region corresponds to the pressure levels to which our data are most sensitive.
In the text

Fig. 8 Astrometric orbit of AF Lep b with 1000 orbits randomly drawn from the posterior. The left panel shows the solutions obtained without the HiRISE RV measurement and the right panel shows the solutions including our measurement. The inset plot in each panel shows the distributions of the relative RV between the planet and the star. As illustrated on the plot, the determination of the RV of the planet offers useful information on its phase.

In the text

	Fig. 9 Predictions of the relative RV between the companion and the star obtained with 1000 orbits randomly drawn from the posterior. As in Fig. 8, the left panel shows the solutions without the HiRISE RV measurement and the right panel shows the solutions including our measurement. The bottom left of the right panel, shows a zoom around the HiRISE measurement.
In the text

	Fig. 10 Posterior distribution for the ascending node (Ω), the argument of periastron (ω) and the mean longitude at reference epoch J2010.0 (λ). Grey: without our RV measurement. Red: with our RV measurement.
In the text

	Fig. 11 First offset of the second night posterior distribution of [Fe/H] and C/O for the Exo-REM/Exo_k model with priors on Teff and log g. Three cases are considered: no prior (blue), Prior on [Fe/H] (orange), Priors on [Fe/H] and C/O (red).
In the text

Fig. A.1 Trace of the CRIRES+ detector for a single order. Top panel: 2D detector image of the companion. The hatched area corresponds to the noise estimation region. The spectrum of the companion is located at the center of the image, around the y-pixel number 100. Around y-pixel 90 is the residual MACAO internal source fiber signal after background subtraction. Around y-pixel number 115 and y-pixel number 85, we can see the signal of 2 of the 3 reference fibers. The bad pixels are depicted in white. Bottom panel: normalized profile between y-pixels 90 and 110 of the signal at x-pixels 340, 480 and 620 for the star (solid lines) and the companion (dashed lines). The extracted signal area is represented by the vertical red dashed lines.

In the text

	Fig. B.1 First night posterior distribution of RV and v sin i for different number of backgrounds used in the data reduction. The number of backgrounds impacts the final S/N of the reduced data, which has a direct impact on the valued derived for the RV and v sin i of the companion.
In the text

	Fig. C.1 Cross-correlation functions for the first offset of the second night.
In the text

	Fig. C.2 Cross-correlation functions for the second offset of the second night.
In the text

	Fig. D.1 Posterior distribution for the orbit of AF Lep b with the astrometry from SPHERE, KECK/NIRC2 and GRAVITY (grey) and with our RV points included (red).
In the text

	Fig. E.1 First night posterior distribution of atmospheric parameters for the Exo-REM/Exo_k model in the case where all parameters are free
In the text

	Fig. E.2 First offset second night posterior distribution of atmospheric parameters for the Exo-REM/Exo_k model in the case where all parameters are free
In the text

	Fig. E.3 First night posterior distribution of atmospheric parameters for the Exo-REM/Exo_k model with priors on Teff and log g. Three cases are considered : No prior (blue), Prior on [Fe/H] (orange), Prior on [Fe/H] and C/O (red).
In the text

	Fig. E.4 First offset of the second night posterior distributions of atmospheric parameters for the Exo-REM/Exo_k model with priors on Teff and log g. Three cases are considered: No prior (blue), Prior on [Fe/H] (orange), Prior on [Fe/H] and C/O (red).
In the text