A new treatment of telluric and stellar features for medium-resolution spectroscopy and molecular mapping Application to the abundance determination on β Pic b

Molecular mapping is a supervised method exploiting the spectral diversity of integral field spectrographs to detect and characterise resolved exoplanets blurred into the stellar halo. We present an update to the method, aimed at removing the stellar halo and the nui-sance of telluric features in the datacubes and accessing a continuum-subtracted spectra of the planets at R ∼ 4000. We derived the planet atmosphere properties from a direct analysis of the planet telluric-corrected absorption spectrum. We applied our methods to the SINFONI observation of the planet β Pictorisb. We recovered the CO and H 2 O detections in the atmosphere of β Picb by using molecular mapping. We further determined some basic properties of its atmosphere, with T eq =1748 + 3 − 4 K, sub-solar [Fe/H]= − 0 . 235 + 0 . 015 − 0 . 013 dex, and solar C/O=0.551 ± 0.002. These results are in contrast to values measured for the same exoplanet with other infrared instruments. We confirmed a low projected equatorial velocity of 25 + 5 − 6 kms − 1 . We were also able to measure, for the first time and with a medium-resolution spectrograph, the radial velocity of β Picb relative to the central star at MJD=56910.38 with a kms − 1 precision of –11.3 ± 1.1kms − 1 . This result is compatible with the ephemerides, based on the current knowledge of the β Pic system.


Introduction
The system of β Pictoris, with its imaged debris disk of dust, evaporating exocomets and two giant planets, is a stunning window on early stages of planetary systems formation and evolution.At the age of β Pic ∼23±3 Myr (Mamajek & Bell 2014), giant planets have already formed, most of the protoplanetary gas has disappeared from the disk, and Earth-mass planets may be still forming.The discovery of β Pic b in direct high-contrast imaging (Lagrange et al. 2009) has been rapidly recognised as a major finding for several reasons.First, it was, until the discoveries of β Pic c (Lagrange et al. 2019;Nowak et al. 2020), and more recently AF Lep b (Mesa et al, 2023), the shortest period imaged exoplanet, allowing thus "fast" orbit characterisation; second, once its mass is known, it can be used to calibrate brightness-mass models and atmosphere models at young ages; third, it is a precious benchmark for detailed atmosphere and physical characterisation thanks to its proximity to Earth and position with respect to the star; and finally, it is an exquisite laboratory to study disk-planet interactions at a post transition disk stage (Lagrange et al. 2010(Lagrange et al. , 2012)).
Model and age dependent brightness-mass relationships predict β Pic b mass to be within 9-13 M J (Bonnefoy et al. 2013;Morzinski et al. 2015;Chilcote et al. 2017).Its mass is still marginally constrained observationnally, because of significant uncertainties on the amplitude of the radial velocities (RV) vari-⋆ Please send any request to flavien.kiefer@obspm.frations induced by the planets b & c.In particular, the available RV data do not cover the whole β Pic b period, the extrema of the recently discovered β Pic c induced variations are not well constrained with the available data, and the RV variations are strongly dominated by the stellar pulsations (see examples in Lagrange et al. 2019Lagrange et al. , 2020;;Vandal et al. 2020).Gaia was used by several authors to further constrain the planet b mass, <20 M J (Bonnefoy et al. 2013), 13±3 M J (Dupuy et al. 2019), 12.7±2.2M J (GRAVITY Collaboration et al. 2020), 10-11 M J (Lagrange et al. 2020).The most recent determinations combine RV, relative and absolute astrometry, taking into account both planets b & c.They lead to 9.3 +2.6  −2.5 M J using the Hipparcos-GaiaDR2 measurement of astrometric acceleration (Brandt et al. 2021), and 11.7 +2.3  −2.1 M J using the Hipparcos-GaiaDR3 measurement of astrometric acceleration with the same datasets (Feng et al. 2022).We note that the astrometric acceleration measurement, also known as proper motion anomaly (see also Kervella et al. 2019Kervella et al. , 2022)), initially 2.54-σ significant using the DR2 (Kervella et al. 2019), became compatible with zero at 0.86-σ using the DR3 (Kervella et al. 2022).This explains the difference in the derived mass.From dynamical considerations, the mass of β Pic b is thus bounded within 9-15 M J .young imaged planets (Snellen et al. 2014;Brogi et al. 2018;Hoeijmakers et al. 2018; Petit dit de la Roche et al. 2018;Ruffio et al. 2019;GRAVITY Collaboration et al. 2020;Ruffio et al. 2021;Cugno et al. 2021;Petrus et al. 2021;Patapis et al. 2022;Mâlin et al. 2023;Petrus et al. 2023;Miles et al. 2023;Landman et al. 2023) allows us to characterize the atmospheric composition in molecules such as CO, CO 2 , H 2 O, NH 3 and CH 4 .The molecular mapping method was first developed for this objective by Snellen et al. (2014), hereafter S14, for medium or high resolution instruments such as CRIRES (S14, Landman et al. 2023), SINFONI (Hoeijmakers et al. 2018;Cugno et al. 2021;Petrus et al. 2021), Keck/OSIRIS (Petit dit de la Roche et al. 2018;Ruffio et al. 2019Ruffio et al. , 2021) ) or JWST/MRS (Patapis et al. 2022;Mâlin et al. 2023;Miles et al. 2023).This method consists in calculating the cross-correlation function (CCF) of a spectrum emitted from the atmosphere of a planet with a theoretical transmission spectrum, or template, using for instance Exo-REM (Baudino et al. 2015;Charnay et al. 2018).This could reveal the presence of individual molecules.The CCF leads to a similarity score, which if equal to 1 (0) means the spectrum and the template are proportional (totally orthogonal).In general because of noise, systematics, and inaccuracies of models, a CCF never reaches exactly 1.Using the CCF as a template matching score in principle allows us to retrieve simple atmospheric properties such as T eff , log g and relative abundances.
For β Pic b, Hoeijmakers et al. (2018) (H18 hereafter) showed using molecular mapping on the cubes collected by the Spectrograph for INtegral Field Observations in the Near Infrared (or SINFONI) that H 2 O and CO were present in the atmosphere of this young planet, with no evidence of other species.They did a tentative template matching of T eff and log g that led only to large confidence regions of those parameters on their Figure 10.They did not produce any estimation of the planet radial velocity nor rotational broadening v sin i.
With a spectrum-fitting oriented approach, an emitted spectrum of β Pic b was obtained with GRAVITY (GRAVITY Collaboration et al. 2020).Its fit led, in a Bayesian inference framework using Markov chain Monte-Carlo sampling of posteriors, to an effective temperature of 1740±10 K with log g=4.35±0.09, a super-solar metallicity [M/H]∼0.7±0.1 and a sub-solar C/O=0.43±0.041 .This was in good agreement with previous estimations of the planet temperature of 1724 K and a log g=4.2 by Chilcote et al. (2017) and the combined astro-metric+RV planet mass estimation ∼12 M J (Snellen & Brown 2018;GRAVITY Collaboration et al. 2020;Lagrange et al. 2020).However, the metallicity was significantly different if considering the GRAVITY spectrum only (-0.5 dex) or combined with the GPI YJH low-resolution spectra (>0.5 dex).Most recently, Landman et al. (2023) published the analysis of new β Pic b high-resolution spectra taken with the upgraded CRIRES+ instrument that led to similar parameters using atmospheric retrievals, with temperatures slightly higher than in GRAVITY Collaboration et al. ( 2020), a sub-solar metallicity (Fe/H∼-0.4)and a sub-solar C/O=0.41±0.04.Allowed by the high-resolution they obtained a new v sin i measurement of 19.9±1.1 km s −1 .
With an approach similar to (GRAVITY Collaboration et al. 2020), Petrus et al. (2021) used both principal component analysis (PCA) and Halo-subtraction on North-aligned angular differential imaging (nADI) on SINFONI observations of HIP 65426 to extract the emitted spectrum of the planet b keeping the thermal continuum.They then used Bayesian inference with nested sampling (Skilling 2006) to retrieve the basic parameters of planet HIP 65426 b from the spectrum itself, including equilibrium temperature, surface gravity, metallicity ratio [M/H] and C/O.This proved possible to derive a spectrum of a planet observed with SINFONI and that having a spectrum-fitting rather than CCF-optimisation method led to more reliable results.
Here, we perform a new analysis of the β Pic cubes observed with SINFONI.We improve the reduction of the cubes as thoroughly explained in Section 2. We then improve the star removal method used by H18 with a different approach, that corrects for residuals from stellar lines.We discuss H18's method, and explain our improvements in the form of a new method called starem in Section 3.2.Then, in Section 4, we apply molecular mapping and compare to H18 results.We further extract the spectrum of the planet in Section 5. We use a simple grid search as well as a Bayesian framework with an MCMC sampling to fit the observed planet spectrum and measure the parameters of the planet.This is done in Section 6.We discuss the results in Section 7 and give our conclusions in Section 8.

SINFONI observations
SINFONI was an infrared instrument, coupling an adaptive optics (AO) module to an integral field spectrograph (IFS) SPIFFI, installed on the Unit Telescope 4 of the Very Large Telescope at Paranal/Chile (Eisenhauer et al. 2003;Bonnet et al. 2004).SINFONI was on-sky from 2004 to 2019.The observations with the SINFONI IFS were performed with different size of fieldof-view (FoV) and spectral resolution (R), and then reduced into data cubes, with two spatial and one spectral dimensions.Here, we focus on observations of the β Pictoris surroundings performed with the 0.8 ′′ ×0.8 ′′ FoV subdivided into 64×64 spaxels2 of size 12.5×12.5-mas 2 at R=4000 along the K-band (2.08-2.45µm).The observations consist in 24 exposures of 60 sec each, recorded on the 10th September 2014 from 08:19:34 UT to 10:05:20 UT.An offset of 0.9-1.1 ′′ from β Pic and a field rotation of -56 • to -19 • was applied, reducing the pollution of the stellar halo upon the planet spaxels with the star decentered outside the FoV, focusing the observations on the surroundings of β Pictoris, and enabling the use of angular differential imaging (Marois et al. 2006).The seeing during the observation varied within 0.8-1.0′′ with an airmass varying from 1.35 to 1.14 between the first and the last exposure.The atmospheric conditions were relatively constant during the observations with fluctuations of pressure and temperature <1%.

From SINFONI raw data to registered cubes
We performed the data reduction of the SINFONI sequence of observations following the scheme described in Petrus et al. (2021) which provides optimally-reduced datacubes for highcontrast science.
The raw data were originally corrected the Toolkit for Exoplanet deTection and chaRacterization with IfS (hereafter TExTRIS; Petrus et al. 2021;Palma-Bifani et al. 2023;Demars et al. 2023, Bonnefoy et al. in prep.) from the so-called "oddeven" effect affecting randomly some pre-amplification channels on SPIFFI's detector (corresponding roughtly to the location of the 25th slitlet).We then used the ESO data handling pipeline version 3.0.0 to reconstruct data cubes from the bi-dimensional science frames.TExTRIS also corrected for the improper registration of the slitlet edges on the detector and for the inaccurate wavelength solution found by the pipeline using synthetic spectra of the telluric absorptions.
Finally, we used TExTRIS to perform a proper registration of the star position outside the field of view.H18 fitted a synthetic function to represent the wings of the star's point spread function (PSF).However, such a method is sensitive to the distribution of flux within the FoV.The later can be affected by (i) the complex evolution of the Strehl ratio which evolved with wavelength and along the sequence, cubes with high Strehl ratio showing strong artefacts due to the telescope spiders while those with low Strehl ratio show a smoother flux distribution, and (ii) the varying part of the halo contained in the FoV due to the field rotation along the sequence.TExTRIS uses instead an initial measurement of the star position in data cubes acquired during short exposures taken at the beginning of the sequence and centered on the star.Then it builds a model of the β Pic centroid positions, which are located outside the FoV in the 24 exposures of the observation sequence, by computing their theoretical wavelength-dependent evolution due to the evolving refraction, the field rotation and the offsets on sky.
A remaining error due to telescope flexure exists (see the ESO user manual) but appears to be below ∼1 pixel in the final registered cubes of β Pictoris.This reduction provides us with 24 data cubes and associated measurements of the offsets and rotation angles that will be used in Section 2.5 to de-rotate and stack the cubes aligned on the position of the planet β Pic b.

Reference stellar spectrum
As it will be used thoroughly in the rest of the study, we define here the method to derive a reference stellar spectrum, free from photons coming from the planet, in the K-band using a SIN-FONI data cube.We found best to use several of the brightest spaxels within a cube and combine them to obtain a stellar spectrum, in order to reduce any pollution from the background and the planet.To find those, we measure for all spaxels of the FoV the flux at the continuum-level of the Br-γ line at ∼2.165 µm by fitting the wings of the line by a 2-degree polynomial and retrieving the level of the continuum at 2.165 µm.From this flux map, we exclude the 10 brighest spaxels to avoid bad spaxels and calculate an average star spectrum from the next 100 brightest ones.Those are the less affected by the background whose level is on the order of ∼20 erg s −1 cm −2 Å −1 while the total flux reaches more than 2000 erg s −1 cm −2 Å −1 .Therefore, its contribution is less than 1% in those spaxels.Fig. 1 shows the absolute total and residual flux distribution among spaxels of the stacked cube obtained in Section 2.5 below.The resulting reference star spectrum is showed on Fig. 2. Note that it includes many telluric lines beyond 2.18 µm, mainly H 2 O, CO 2 and CH 4 lines.

Wavelength calibration correction
The presence of telluric lines is a nuisance to the analysis of stellar and exoplanet spectra.Nonetheless, they can also be used to adjust the wavelength calibration in the SINFONI cubes.Tellurics can be fitted directly in each of the 24 cubes to the normalised star spectrum.Since we are most interested in the planet, the star spectrum is here obtained from the spaxels located at the  position of the planet.The planet position in the derotated cube is determined in Section 4.2 and its PSF of 4-spaxels FWHM in Section 3.2.We calculated for each cube the mean stellar spectrum on a circular area of 6-spaxels radius around the planet location.
We used the ESO code molecfit (Smette et al. 2015;Kausch et al. 2015) v3.13.6 that implements LBLRTM to perform the telluric line fit in this spectrum.Typical site parameters during the observations are taken as inputs, such as MJD, Paranal altitude and coordinates, humidity (∼4%), ambiant pressure (∼740 hPa), ambiant temperature (∼12 • C), mirror temperature (∼10.9 • C) and airmass (sec z∼1.1-1.4).molecfit fits the atmospheric parameters -such as e.g.water column abundance, pressure, temperature etc. -as well as a continuum, a Chebychev polynomial wavelength solution, and a line spread function to the observed telluric lines in the observed spectrum.The error bars are fixed to the square root of the flux divided by the normalisation.The observed reduced χ 2 are consistent within 1.3-1.4all through the cubes time series.We found that the molecules H 2 O, CO 2 , CH 4 and N 2 O dominate the model tellurics spectrum in the K-band, whilst CO, NH 3 , O 2 and O 3 are always either weak and undetermined or fitted to negligible relative abundance values < 10 −4 .To reduce the computation effort, we thus only fit for H 2 O, CO 2 , CH 4 and N 2 O column densities.A polynomial of degree 6 for the fit of the continuum and of degree 1 for the fit of the wavelength solution was adopted.We fixed the LSF to a Gaussian function allowing its width to vary.The LSF is moreover convolved in molecfit by a 1-pixel (0.00025 µm) box to mimic the effect of the slit smearing.Fig. 3 shows a summary of the molecfit solutions along the cubes.Fig. 4 shows the β Pic stellar spectrum compared with the molecfit resulting model.There is a residual timedependent shift of the wavelength solution even after the absolute calibration performed by TexTRIS of about 8 km s −1 from first (cube 0) to last (cube 23) exposure.This agrees with the magnitude of the error on the calibration found in Petrus et al. (2021).Such shift could be due to an effect of flexure of the instrument.We corrected the wavelength solution in all cubes according to this analysis.The figure also shows that the resolving power is varying through the observation series with an LSF FWHM of ∼2.16-2.32pixels at 2.27 µm.This variation is due to the wandering of the planet image on the detector.This leads to an average effective resolving power in our SINFONI K-band spectra of R=4120±90.
The atmospheric molecular column density of H 2 O is relatively constant along the observation sequence.We note that the CO 2 abundance level of ∼500 ppmv retrieved is ∼1.4 times higher than the reference value (∼370 ppmv) in the model Earth atmosphere, while N 2 O abundance ∼0.30 ppmv is about 1.2 its reference level.At the same date and using the whole K-band spectrum (including deep CO 2 and N 2 O lines below 2.07 µm) of another reference star, the CO 2 and N 2 O abundances reach respectively ∼390 ppmv and 0.23 ppmv (Smette, private comm.).Relying on weak lines only, our determination of the CO 2 and N 2 O abundance levels might not be well determined.The absolute values of the species abundances presented here should thus be considered as indicative only.
Applying the same procedure on the star spectrum taken from other locations in the image led to similar molecular abundances.However, it revealed a strongly scattered resolving power within 3500-4900, though the average value of R agreed on an effective resolving power of ∼4000.This is much smaller than the theoretical R=5950, expected for SINFONI in this setup.This might be explained by the degradation due to pixelation when creating the cubes, since the LSF sampling is sub-Nyquist, with a spectral broadening of about 1.5 pixels as noted in the latest SINFONI's documentation3 .
A thorough investigation on how to improve the effective spectral resolution in the reduced data of instruments such as SINFONI is not in the scope of the present study.Nevertheless, in regards of such objective, two leads might be worth mentioning, that are, either i) achieving a finer reconstruction of cubes from the 2D images of the slitlets and the arc lamp calibration, or more simply, but time-consuming, ii) using dithering during the observations, in order to better sample the LSF, at the cost of at least doubling the exposure time.These leads will be explored in further studies.

Science cubes stacking
We align the 24 science cubes taken on the 17th of November upon the star centroid.Then, the cubes are de-rotated in such a way that the planet halo is brought at the same (α * ,δ)-coordinates in every cube and at every wavelength.We use the values obtained with TExTRIS as explained in Section 2.2.This includes a 2D-linear interpolation using the interpolate.griddataroutine from scipy in order for all the shifted-rotated cubes to share a common (α * ,δ)-grid.Then the 24 cubes are stacked together using a simple average.No clipping of flux is applied during this process in order to maximise the signal-to-noise ratio.This gives us the data cube that we use in the rest of the analysis; we name it the "master cube" and it is noted Q obs .

A summary of H18 method
The H18 method used to remove the stellar halo proved to work well for performing molecular mapping of exoplanets.However, it does not allow the extraction of a pure planet atmosphere transmission spectrum where systematic deviations remain.We will summarize the H18 method here to show where is the identified issue.A sketch of the different steps of this method is shown in Fig. 5.
The components of the flux at the planet spaxels i p are as follows: (1) Among which the star flux F i p ,⋆ , the planet flux F i p ,p and a possible background component F i p ,B that we will put aside for the simplicity of the demonstration -it can be added everywhere by duplicating the planet components and changing 'p' by 'B' in the index -and will be discussed in Section 3.3.We will also drop the 'i p ' index in the following to alleviate the equations.
The planet and star spectra are each composed of a continuum, hereafter noted C, multiplied by a 'flat' transmission spectrum whose continuum is normalised to 1 everywhere, hereafter noted η.Both the astrophysical source and telluric lines contribute to this transmission spectrum.It can also be expressed as η=1− A, with A a positive comb of spectral lines with the continuum equal to zero everywhere.We point out that A ⋆ and A p for respectively the star and the planet, are supposed to be spaxelsindependent since intrinsic to the respective sources.Continua C p and C ⋆ are on the other hand spaxel-dependent, because of the point spread function and the wavelength-dependent speckles.
In order to remove the star and single-out the planet, H18 subtract from each spaxel i of the cube Q obs an approximated star spectrum F i,⋆,approx that accounts for wavelength-dependent spaxel-to-spaxel variations.To do so, they divide each spaxel spectrum by the star spectrum determined, as explained in Section 2.3, by averaging some of the brightest spaxels of the master cube.The 'star-free' cube Q star-free thus obtained show low-frequency wavelength-dependent variations that differ in the spectra from one spaxel to the other.They are due to speckle patterns over the detector changing with wavelength.Those speckle patterns are modeled by applying a Gaussian filter G on this star-free spectrum, G * Q star-free .It results in a speckle-proxy cube noted Q spkl in Fig. 5.By multiplying those modelled variations by the star spectrum, they finally obtain the star cube Q ⋆ with at each spaxel i a star spectrum F i,⋆,approx that accounts for wavelength-dependent spaxel-to-spaxel variations.
At any spaxel other than the planet spaxels, subtracting F i,⋆,approx removes the contribution of the star spectrum.However, at the planet spaxel, this is not true.Indeed, the approximated star spectrum contains contributions from both the star and the planet continuum: By slightly overestimating the star contribution, the subtraction of F ⋆,approx rather leads to A supplementary stellar contribution -including telluricsto the residual absorption spectrum remains as emission lines with amplitudes comparable to the planet absorption lines.This is shown in Fig. 6.Note that the persistence of polluting CH 4 lines led H18 to apply several runs of SYSREM (Tamuz et al. 2005) in order to remove them.But this operation does not fix the above issue and stellar lines still remain in the spectrum at the planet spaxels.
The presence of the star spectrum with a non-negligible amplitude is problematic.Moreover, the continuum being subtracted, there are no more reference level to which comparing the strength of molecular lines in the resulting spectrum.As long as the CCF of the star spectrum and the spectra of any of the species found in exoplanet atmosphere is close to zero, this works well for molecular mapping.This is the case here with a 8000 K star and a <2000 K planet.Still, ∆F is not strictly-speaking a pure planet atmosphere spectrum.This issue might also explain the detection of an Br-γ emission line in the PDS-70 b spectrum derived with the same H18 method in Cugno et al. (2021): we now suspect it is an artifact from the stellar spectrum removal.

STAREM: a new STAr spectrum REMoval method
We propose a different method to subtract the star spectrum, that we named STAREM, that rather makes use of the normalised transmission spectrum η=F/C.In this spectrum, we are going to estimate the contribution of the star in any spaxel spectrum F i by fitting stellar absorption lines, and then subtract it from F i .We show this can lead to a well defined flattened transmission spectrum of the planet atmosphere, from which the star spectrum is fully removed, and in which the line strength is preserved.
First of all, we normalise all spectra of the observed cube, as well as the star spectrum, by fitting a 6th degree polynomial to their continuum, leading to a normalised star spectrum F ⋆,norm and a normalised cube Q obs,norm .Recalling equation (1), the normalised transmission spectrum components of F = C η at one of the planet spaxels in Q obs are Fig. 6.The raw β Pic b spectrum extracted from the planet spaxels after using H18 method of star light removal.SYSREM was not used here.We highlight in red the main artifact due to excess star spectrum lines removal at the Brackett-γ wavelength.Also should be noticed that the continuum being subtracted, the reference level to which compare the strength of planet molecular lines (CO explicited in green) is missing.
We recall that the star and planet pseudo-continua are fixed by the star and planet intrinsic pseudo-continua multiplied by the Strehl ratio and PSF (including speckles) damping of the flux.We introduce the star and planet contribution levels K ⋆ and K p as The transmission spectrum η(λ) can be expressed more compactly as Wishing to subtract the star contribution K ⋆ η ⋆ from this spectrum, we thus need to estimate K ⋆ (λ).This can be achieved by comparing the amplitude of the stellar lines to those of a reference stellar spectrum without contributions from the planet, a ratio of 1 implying a pure stellar spectrum, that is K ⋆ =1.Ideally, with several stellar lines present all through the observed band, the best approach would be to use all the lines so as to obtain a more reliable wavelength-dependent approximation of the K ⋆ (λ) function.In the specific case of the K-band spectrum of β Pic, since the Brackett-γ line at 2.165 µm is the only strong feature in this spectral band, we are only able to derive an approximated constant star contribution ∼K ⋆ (λ 0 =2.165 µm), or K ⋆,2.165 hereafter.The derived values of K ⋆,2.165 and its pendent residual contribution, 1-K ⋆ , through the SINFONI field-of-view around β Pic b are shown in Figure 7.
Removing an approximated contribution K ⋆,2.165 η ⋆ at each spaxel of Q obs,norm leads to a residual cube Q res .In this cube, the residual spectrum obtained on a planet spaxel is Using eq. 5 this could equivalently be written as: This is a flat spectrum whose baseline level is K p,2.165 =1-K ⋆,2.165 .At the planet spaxels, it should contain the normalised spectrum of the planet with an amplitude corresponding to the relative flux of the planet compared to the star at 2.165 µm.Because of the normalisation of the planet spectra by the continuum of the star, this is however only an approximation.Indeed, away from Brackett-γ -more generally from any fitted stellar line -the line amplitudes are impacted by a residual star spectrum component with an amplitude δK(λ)=K p,2.165 -K p (λ).This is a footprint of the planet and the star continua within the flat spectrum.Unfortunately, we cannot assume that δKη ⋆ ≈0, because it generates a warp of the planet spectrum lines deviating from 1 by a few tens of percent, as shown in Figure 8.This will have to be taken into account later on when analysing the spectrum extracted from the planet spaxels.

Tellurics and background removal
In Section 2.4, tellurics were fitted directly in each of the 24 cubes to the normalised reference star spectrum obtained from the brightest spaxels.They were only used for purpose of calibration.We did not wish to remove tellurics at this step because removal might introduce residuals with amplitudes on the order of magnitude, or worse larger, than the planet features.
Having in hand the residual cube Q res that features mostly the planet and background, we now wish to remove the many tellurics lines that are still present in the observed spectrum in the K-band.Most straightfully we could fit a model telluric spectrum to any spaxel spectrum in the cube using molecfit.This is however quite time consuming to run molecfit with the cube featuring about 10,000 spaxels, and moreover the S/N of the resid- ual spectra is poor.Instead, we rather use molecfit on the star spectrum derived from the stacked cube at the position of the planet.This allows us to obtain a reliable telluric model η model , to correct again the wavelength solution, using the same settings as in Section 2.4.Doing so, we determined the effective spectral resolution at the position of the planet in the stacked cube to be 4020±30.
At each spaxel i of the residual cube, we fit out this telluric model η model from η i,res .To do so, we first need to correct for a wavelength-dependent background contribution whose non-zero level is causing an artificial decrease of telluric lines prominence and the appearance of spurious features, such as telluric lines residuals.Those can be seen by comparing a planet-spaxel spectrum and its surrounding background (Fig. 9).The origin of the continuum level of the background is unknown but it could be due to scattered stray light.The features are generated by the differences between the stellar (plus telluric) spectrum actually observed at the given spaxel and the reference stellar (plus telluric) spectrum that was removed in Section 2.34 .
We make an estimation of the background contribution in a spaxel i by taking the median spectrum in a distant ring around this spaxel η i,ring .Since the SINFONI PSF has FWHM∼4 spaxels in radius, we use a ring radius of 6 spaxels with a width of 1 spaxel.This background estimate is then removed from the spaxel i spectrum η i,res,corr = η i,res − η i,ring .Its contribution level at ∼2.165 µm can be viewed in Fig. 7 on the brighter areas away from the planet location.Fig. 1 also shows the 2.165-µm absolute background flux that is ∼1-10% of the absolute total flux at every spaxel.Around the planet location its continuum level is very close to zero, but non-negligible features are seen that have corresponding features in the spectrum at the planet central spaxel.
Adding a background term to eq. 8, the background subtraction can be summarized with the following set of equations: Assuming an equal contribution of background in the central and the neighboring spaxels, i.e.K i,b = K ring,b , this equation simplifies to: We recover Eq. 8, and the previously implicit background is now explicitly suppressed.
The final step is to fit η model the telluric spectrum model to η i,res,corr .We use a simple least-squares optimisation at each spaxel, allowing the telluric model to vary in intensity with η i,tell = a η model + b with a + b = c i if c i is the level of the residuals at spaxel i.Because η model was determined from the reference stellar spectrum that could be a little affected by background, an offset b is added to account for small intensity differences in telluric lines between η model and the telluric spectrum at spaxel i.We divide out η i,tell from η i,res,corr leading to the telluric-free spectrum η i,res, tell-free .
A summary of background subtraction and telluric removal at the central planet spaxel is shown in Fig. 9.

Exo-REM grid of models
To calculate the molecular mapping of diverse species in the atmosphere of β Pic b, as is now most usually performed (see e.g.Snellen et al. 2014;Brogi et al. 2018;Hoeijmakers et al. 2018; Petit dit de la Roche et al. 2018;Ruffio et al. 2019;Cugno et al. 2021;Petrus et al. 2021;Patapis et al. 2022;Mâlin et al. 2023), we calculate a CCF of the observed spectrum at each spaxel with a synthetic spectrum.Here, the templates are taken from an Exo-REM model grid (Baudino et al. 2015;Charnay et al. 2018;Blain et al. 2021).Exo-REM is a 1D radiative-convective model, which self-consistently computes the thermal structure, the atmospheric composition, the cloud distribution and spectra.The model includes the opacities of H 2 O, CH 4 , CO, CO 2 , NH 3 , PH 3 , TiO, VO, FeH, K and Na, and collision-induced absorptions of H 2 -H 2 and H 2 -He.Silicate and iron clouds are included using simple microphysics to determine particle sizes (Charnay et al. 2018).The model grid includes four free parameters totaling 9573 spectra with T eq ranging from 400 to 2000 K with steps of 50 K; log g from 3.0 to 5.0 dex with steps of 0.5 dex; [M/H] from -0.5 to 1.0 dex with steps of 0.5 and [C/O] from 0.1 to 0.8 with steps of 0.05.The synthetic spectra were computed at R=20,000, a compromise between reducing computation speed and reaching the highest resolution possible on atmospheric models to be used as templates of low-to-medium resolution spectra of exoplanets.
We cleansed the model grid from those that did not converge well.To do so, we calculated for each spectrum the integral on the whole wavelength domain from visible to far IR, I s .We compared this integral to the theoretical σ T 4 eff Stefan-Boltzman law.This deviation peaks below ∼5% (Fig. 10).We thus removed 3315 among 9573 spectra (i.e.35%) with a deviation larger than 5%.
Each synthetic spectrum produced by Exo-REM results from individual contributions (or spectrum) of the species out of which it is composed, mainly CO, H 2 O, CH 4 and NH 3 , which we can use as individual molecular templates.

The cross-correlation function of the planet spectrum
At every spaxel, between 2.08 and 2.43 µm, we calculate the CCF of the observed spectrum and an Exo-REM synthetic spectrum for an assumed T eff =1700 K, log g=4.0 cgs and [Fe/H]=0.0dex, as based on Table 3 in GRAVITY Collaboration et al. (2020).The expected abundances of NH 3 and CH 4 in the atmosphere ∼10 −6 are too low to yield absorption features detectable in the SINFONI spectra.We thus artificially enhance their abundances to 10 −4 , in order to probe possible over-abundance of these species in β Pic b's atmosphere.Details on the contributions of H 2 O, CO, NH 3 and CH 4 are shown in Fig. 11 and compared to the β Pic b's spectrum derived in next Section 5.
We exclude the red edge of the K-band beyond 2.43 µm which displays the strongest telluric lines remnants.Prior to the calculation, we divide the continuum of the observed and synthetic spectra.The continuum is obtained by first fitting a 4th-degree polynomial, and then applying a median filter with a window-width of 0.01 Å, combined to a smoothing Savitzky-Golay filter of order 1.The CCF is calculated by directly crosscorrelating the median-removed observed and synthetic spectra at different shifts.The CCF is finally normalised by the norm of the spectra, thus leading to a zero-normalised CCF.This results in the molecular maps shown in Fig. 12.We then fit the PSF of the planet CCF halo by a 3D Gaussian with respect to (δ, α, v r ).This determines the position and radial velocity of the planet, as well as the broadening of the PSF and of the LSF.The results are summarized in Table 1.They are compared to the values derived in the same fashion but following H18's recipe to suppress the stellar pollution.We only show results for CO, H 2 O and using the full spectrum, confirming, as in H18, the presence of CO and H 2 O but not detecting NH 3 and CH 4 .With some variability from one model to the other, the CCF peak is located at a separation to the central star of ∼351±5 mas, with a PSF broadening standard deviation of ∼2.7 spaxels, that is 34 mas.
A signal-to-noise (S/N) ratio is calculated as the height of the peak divided by the fit residuals on a volume of 200 pix 2 by 2000 km s −1 around the planet (δ,α,v r )-location.Our starem method leads to S/N comparable to those obtained through H18 method.We did not use SYSREM as they did to remove spaxelto-spaxel correlated noise within the cubes, but corrected for the background differently as explained in Section 3.3.We derived slightly better S/N by subtracting the star halo using STAREM than those obtained using the H18 subtraction method, leading to S/N all =19.6,S/N CO =11.8, S/N H 2 O =15.4.This follows the trend found also for HIP 65426 b (Petrus et al. 2021) where taking into account all species leads to the best detection S/N, while for individual contributions of species, H 2 O gives significantly better results than CO.This last properties is best explained by the presence of less prominent but more numerous lines of H 2 O compared to CO all along the K-band, and especially away from regions spanned by telluric lines.
We note that in the H18's paper, the S/N they obtained are larger (all molecules: 22.8 vs. 17.5 here; CO: 13.7 vs. 11.2here; H 2 O: 16.4 vs. 15.7 here).This difference can be explained, because they averaged the CCF of cubes obtained at two different nights, while we used the cubes from only a single night.The smaller difference of S/N for H 2 O can be explained by the presence of traces of telluric lines in the H18 residual cube because they did not recalibrated the wavelength solution before subtracting the reference spectrum (containing the tellurics).We also note that the calculation performed in H18's paper to determine the S/N is slightly different in that they use a distant 3D-ring around the CCF peak to estimate the noise, while we calculate the noise from the residuals of a fit of the CCF peak.

Extracting the planet atmosphere absorption spectrum
The PSF of the planet follows an Airy profile whose central region can be approximated by a Gaussian profile.This can be seen in Fig. 13 where the residual flux of the stellar subtraction at 2.165 µm is shown with respect to the distance of the spaxels to the planet center position.The planet center position is obtained by fitting the residual flux at 2.165 µm by a 2D-Gaussian.
The radial profile of the planet PSF is well modeled by a 2-spaxel wide Gaussian.The planet relative flux rises up to 7% near the center but becomes negligible compared to noise beyond a distance of ∼5 spaxels.
The average spectrum at those spaxels is calculated by integrating the relative flux on all spaxels up to a given radius and applying to each spaxel a weight that depends on its distance to the planet.We model the PSF profile by the Gaussian used above, with F(r)=e −r 2 /2σ 2 PSF where r=∥r i − r c ∥ is the spaxel-distance between a given spaxel i and the planet center position, and σ PSF =2 spaxels at 2.165 µm.We use this profile as the weight function with thus w(r)=e −r 2 /2σ 2 PSF .When integrating the spaxels flux, both planet and noise signals grow, yet limited by the PSF flux dimming, and the signal-to-noise increases like Notes.
Since we cannot assume that the background removal was absolutely perfect in Section 3.3, we again estimate the background pollution within planet spaxels using neighboring spaxels in a ring around the planet.The background level at distance to the planet centroid is close to zero within ±0.01, as can be seen in the PSF profile, Fig. 13.We use the median spectrum within this ring as an estimation of the background spectrum shown in the top panel of Fig. 14.This background spectrum is subtracted from the average planet spectrum, leading to the corrected planetary spectrum, also shown Fig. 14.We found the ring radii optimising the S/N of the corrected planet spectrum to be 4.2 -5.1 spaxels.
As an ultimate step, the spectrum was divided by the continuum, estimated, as in Section 4.2, by applying a median filter with a window-width of 0.01Å combined to a smoothing Savitzky-Golay filter of order 1.The final planet spectrum obtained is compared to a theoretical Exo-REM spectrum of a 1700 K planet with log g=4.0 cgs and [Fe/H]=0.0dex in Fig. 11.

Template-matching of the planet spectrum
We addressed the problem of finding the best matching model of the planet spectrum in two ways: first by exploring the space of available templates and their match with the planet spectrum using a simple grid search, second by running a Markov-Chain Monte Carlo (MCMC) sampling around the optimum spectrum.In both cases we used the forward-modeled Exo-REM (Charnay et al. 2018) spectra models of exoplanet atmosphere.

Grid search
We first explored the template space by a grid search.We tested the χ 2 and zero-normalised CCF scores for optimizing the models.Each model spectrum is convolved by a Gaussian Kernel corresponding to a resolving power of 4020 and by a rotational pro-file with v sin i=25 km s −1 .The rotational profile is the usual belllike profile with limb-darkening coefficient ε=0.6 (Gray 1997).Each model spectrum is then flattened with the same median filtering function, with a window width of 0.01 µm, as used to flatten the observed spectrum.We moreover corrected for the "warp" effect noted in Section 3.2 in the models by adding a supplementary term in (1 − C p /C) η ⋆ , with η ⋆ the normalised star spectrum, C p the model's continuum and C the average continuum in the SINFONI cube at the planet spaxels both normalised to 1 at 2.165 µm.For the ZNCCF the median of both model and data is subtracted before cross-correlation.For the χ 2 , the fit function includes two other parameters applied to the model spectrum, a scale a to allow compensate for the arbitrary level of the spectrum continuum, and a rigid Doppler shift v r .The fit is then performed using the curve_fit procedure from the scipy library.
The error bars of data points are derived from the relation σ data ∝ F ⋆,approx /C, with F ⋆,approx (λ) the reference stellar spectrum defined in Section 2.3 and C(λ) the continuum by which spectra are divided in Section 3.2.The level of noise is normalised to the error directly measured in the planetary spectrum at 2.165 µm from the standard deviation of flux to the median on a width of 0.003 µm ; we take the average of the errors at all spectral channels from 2.15 to 2.18 µm.
The grid search maps with χ 2 -score and CCF-score are shown on respectively Fig. 15 and 16 and the results are summarized in Table 2.For the χ 2 the confidence intervals are bounded by the ∆χ 2 with the 1, 2 & 3-σ regions corresponding to ∆χ 2 =2.3, 6.2 & 11.8.For individual measurements in Table 2 the 1 and 2-σ intervals correspond to ∆χ 2 =1 & 4, interpolated from the 1D ∆χ 2 obtained for each parameter.The small extent of the 1-σ confidence regions around the best fit model is most likely an effect of the discrete grid used to explore the parameter space.Half the 2-σ intervals are certainly more reliable than the 1-σ intervals as error bars.For the CCF grid search we show the ∆CCF regions of 1, 10 and 20% ; they are not translated into confidence intervals and only the optimum model is given in Table 2.
The χ 2 grid search leads to a solution with T eff =1550 K, [Fe/H]=0.0dex, log(g)∼3.5,and C/O∼0.70.The best matching model is compared to the observed spectrum in Fig. 17.The CCF grid search leads to a different solution with much larger temperature saturating at 2000 K and a solar metallicity.But in this case, as shown in Fig. 17 the spectra do not match in line amplitude.This shows that the CCF, due to the removal of the continuum and the normalisation of the spectra, is not well adapted to grid search with a strong degeneracy of models.

An MCMC sampling of models
We use emcee (Foreman-Mackey et al. 2013) to apply 1,000,000 iterations of Markov-Chain Monte-Carlo (MCMC) sampling around the best solution found by gridsearch in Section 6.1 to assess the true posterior distribution and correlations among the varied parameters.The fitted physical parameters are T eff , log g, [Fe/H], C/O, v r and v sin i.We also fit for the spectral resolving power R spec , in order to account for possible deviations from R=4020±30.We also sample a rescaling a to adjust the continuum of model and observed spectra together.The preprocessings of data and models were the same as those used adopted in the gridsearch.The error bars of the data points σ data (as defined in Section 6.1) might be over or under estimated by a certain multiplicative amount f err and under-estimated by a con-  stant jitter σ.We thus use the error model σ 2 err = ( f err σ data ) 2 + σ 2 with σ data depending on wavelength.We add f err and σ as hyperparameters in the MCMC.The prior distributions, either normal N or uniform U, assumed for all parameters are listed in Table 3.All parameters, except R spec , follow uniform priors.
The Exo-REM model spectra are interpolated through the initial 4D-grid at the specific values of parameters chosen by the Table 2. χ 2 and CCF grid search results (see Section 6.1.The reduced optimum χ 2 has been normalised to 1 and is therefore not given here.Uncertainties are obtained by interpolation of the 1D ∆χ 2 maps.Notes. ( †) The 1-σ error bars for v r , v sin i, and a are obtained for the optimal model using scipy's curve_fit procedure.The v sin i and a are not optimized for the CCF, assuming v sin i=25 km s −1 and a=1.MCMC walkers at each new step.For a quick calculation, we apply the LinearNDInterpolator from scipy on the smallest Hull simplex, where each vertex is a model of the Exo-REM grid, surrounding each MCMC sampled points.This accounts for missing models at several grid nodes (see Section 4.1 above).
The MCMC runs for at most 1,000,000 steps with 20 walkers.It stops whenever the number of steps is smaller than 50× the largest auto-correlation length of the samples.Table 4 sum 4, in this case, the MCMC leads to a larger temperature of 1746 +4 −3 K, and a smaller C/O of 0.551±0.002compatible with a solar value.
We also tried fixing T eff 's prior to Chilcote et al. (2017)'s value 1724±15 K and both log g and T eff .These trials are also shown in Table 4 and are compatible with the fixed-log g trial.We calculated the Akaike information criterion AIC=2k − 2 ln L for all models.The maximum likelihood estimator (MLE) model obtained with the T eff fixed minimizes the Akaike information criterion (AIC) and is thus preferred over all other models.The large difference in AIC of ∆AIC∼36 means that the T efffixed model is 2.5×10 −8 times more likely than the model with uninformed priors, and 1.4× more likely than the model with both log g and T eff fixed.This preferred solution has a T eff of 1748 +3 −4 K, a log g=4.22±0.03slightly larger than the Chilcote's value of 4.18, a sub-solar Fe/H=-0.235+0.013  −0.015 dex, a solar C/O ∼0.551±0.002,and a v sin i=25 +5 −6 km s −1 .In all cases, we found a jitter σ that is compatible with zero, with a correcting factor of error bars f err ∼1 +0.05 −0.11 very close to 1.Both imply that our determination of flux error bars σ data introduced in 6.1 is reasonable.Fig. 19 shows our preferred model compared to the observed planet spectrum, and the corner plot of posterior distributions is shown in Fig. 18.
The radial velocity of the planet is found to be around 0.6±0.9km s −1 in the Earth reference frame.With a barycentric correction of 8.1 km s −1 this leads to a radial velocity of planet β Pic b of 8.7±0.9 km s −1 .Compared to β Pic systemic RV of ∼20±0.7 km s −1 (Gontcharov 2006), it implies an RV of the planet relative to β Pic's central star of -11.3±1.1 km s −1 at MJD=56910.38.This value agrees at 0.3-σ with the predicted RV of the planet at this MJD extrapolating the ephemerides of β Pic b's orbital motion from the RV+astro+imaging solution of Lacour et al. (2021) of -11.6 km s −1 .It also compares well to the RV measured at high resolution using CRIRES a year earlier, -15.4±1.7 km s −1 at MJD=56,644.5 (Snellen et al. 2014).

A new estimation of v sin i
In this work, we found a v sin i of 25 +5 −6 km s −1 in good agreement with the S14's result v sin i=25±3 km s −1 and the most recent Landman et al. (2023)'s v sin i estimation of 19.9±1.1 km s −1 .The larger confidence region of our measurement is explained by the much lower resolving power of SIN-FONI (R=4030±30 as determined in Section 2.4) compared to that of CRIRES (R=75,000).
S14 noticed that given the mass of the planet, the spin velocity of 25 km s −1 was too low to comply with the log-linear massspin law followed by Solar System objects from Jupiter down to asteroids (e.g.Hughes 2003), that would imply a spin velocity of ∼50 km s −1 .Figure 20 shows this relationship for the Solar system planets, with a log-linear law well-fitted to Earth+Moon5 , Mars, Jupiter, Saturn, Uranus and Neptune mass and equatorial speed, with values taken from Hughes (2003).Mercury and Venus are recognized as deviating from this law due to a loss of momentum during their lifetime through tidal interactions with the Sun (Fish 1967;Burns 1975;Hughes 2003).
The difference in v sin i for β Pic b compared to the expected equatorial speed at a mass of ∼11 M J is mostly explained by the young age of the planet (∼23 Myr; Mamajek & Bell 2014).Indeed, the planet is currently contracting from ∼1.5 R J down to ∼1 R J (S14, Schwarz et al. 2016) and its spin should thus be accelerating.Fig. 21 shows the effect of dilatation/contraction through time on the equatorial radius from two different evolution models, ATMO2020 (Phillips et al. 2020) and Baraffe-Chabrier-Barman (Baraffe et al. 2008).According to the conservation of momentum, its spin velocity is expected to increase up to 40±20 km s −1 at 4.5 Gyr in better agreement with the solar system law.Fig. 22 shows the predicted v sin i at 23 Myr evolved backward using the ATMO2020 models, down from an equatorial velocity determined by the spin-mass law at 4.5 Gyr shown in Fig. 20.It can be seen on Fig. 22 that even taking into account contraction there is still a tension between the v sin i measurements of β Pic b and its mass, especially if the mass is contained within 10-14 M J .The v sin i and mass overlap only over regions for v sin i>22 km s −1 , with a mass either <10 M J or >14 M J .These masses are marginally supported by the combined astrometric and RV measurements that rather favour a planet mass within 9-15 M J (Dupuy et al. 2019;GRAVITY Collaboration et al. 2020;Lagrange et al. 2020;Vandal et al. 2020;Brandt et al. 2021;Feng et al. 2022).
This remaining discrepancy between mass and v sin i could be explained by the random aspect of moment exchange during planet formation.However, it might also be the hint of a tilt of the planet's equator compared to its orbital plane.In such case, the projected spin velocity v sin i is smaller than the true equatorial velocity.A direct compatibility at 1-σ of the predicted v sin i at 23 Myr for a mass at 11 M J and our measurement 25±6 km s −1 leads to a tilt compared to edge-on >15 • .

A solar C/O
We derive here a value of the C/O for β Pic b that is solar, 0.551±0.002,while GRAVITY Collaboration et al. ( 2020) and Landman et al. (2023) found a sub-solar C/O of respectively 0.43±0.05and 0.41±0.04.Forming β Pic b in situ along the core accretion scenario (Pollack et al. 1996) would imply a C/O ratio in the atmosphere of the planet largely super-solar >0.8, because of the expected abundances of the different gases in the disk from one ice line to the other when assuming a disk with a static composition all through the main phase of planet formation (Öberg et al. 2011).In this framework, a proposed scenario to reach solar and sub-solar C/O is to consider accretion of icy planetesimals from beyond the ice lines, where most of the H 2 O and CO 2 of the disk is in condensed phase.Alternatively, a solar C/O is more naturally reached by forming either by gravitational instability (Boss 1997) anywhere in the disk, or by core accretion close to the H 2 O ice line with a moderate planetesimal accretion followed by an outward migration.Given that core accretion is preferred for most compact planetary systems with also terrestrial planets -as e.g. the Solar system -and that the β Pic system has at least two planets, with one within 5 au (Lagrange et al. 2019) plus small km-sized icy bodies (Ferlet et al. 1987;Beust et al. 1990;Kiefer et al. 2014;Lecavelier des Etangs et al. 2022), we further consider this scenario as the most likely for forming the β Pic's planets.
In this framework, we can estimate the location of the ice lines compared to β Pic b location in the disk.Following Öberg et al. (2011), the typical temperature profile in a protoplanetary disk is given by (Andrews & Williams 2005, 2007): Here, T 0 is the average temperature in the disk as if it were located at 1 R ⋆ from the center of the star.Given the effective temperature of β Pic (Saffe et al. 2021), we have that T 0 =8000 K. Considering the typical evaporation temperatures of H 2 O, CO 2 and CO summarized in Table 5 with R ⋆ =1.7 R ⊙ (Kervella et al. 2004), we derived the radii of the different ice lines around β Pic.They are shown in Table 5 as well.It results that with a b ∼9.8 au (Lagrange et al. 2020), planet b is located between the H 2 O (6±1 au) and the CO 2 (31±5 au) ice lines.Fig. 23 shows the variation of the ice lines through time, as derived using Dart-mouth12 evolutionary tracks for pre main sequence stars.Notes.
(⋆) For σ, given the shape of the posterior distribution peaking close to 0 and compatible with 0 at less than 2-σ, we only give the upper-limit at the 84th percentile.
( †) max(τ λ ) is the maximum auto-correlation time among all varied parameters λ.  (Holweger et al. 1997).Our C/O of 0.551 +0.003 −0.002 is in good agreement with the expected solar C/O for a solar metallicity.The low metallicity of the planet that we found, -0.235 +0.015 −0.013 , is also in agreement with new-generation planetary population synthesis (NGPPS, Schlecker et al. 2021) simulations through core-accretion, on the same range of semi-major axis ∼10 au, especially if the mass of β Pic b is larger than 10 M J .The correlation between bulk metallicity and planet mass, as obtained through the NGPPS simulations performed around stellar hosts with metallicities ranging from -0.5 to 0.5, can be seen in Fig. 10 of Petrus et al. (2021).
To explain the C/O value found when assuming a static disk composition during planet formation, as in Öberg et al. (2011), we propose that i) planets b & c underwent an inner migration during the first Myr of their formation, allowing them to gather large amounts of gas with solar composition; followed by ii) an outward migration, with planet b & c in 7:1 mean motion resonance leading them to reach their actual location (Beust et al., priv. comm.).This scenario would allow avoiding a fine-tuned icy planetesimal accretion within the planet to reach an almost perfect solar C/O.
Alternatively, there exists a scenario that do not require invoking outward migration of β Pic b.Considering non-static disk composition, Mollière et al. (2022) found that pebble evaporation and the dilution of water and CO in-between the H 2 O and CO iceline (Fig. 6 in Mollière et al. 2022) can lead to a nearly stellar C/O ratio in the circumstellar gas in about 1 Myr.That could enable in-situ formation of β Pic b provided most of its gas was not accreted before 1 Myr.
To conclude, a solar C/O for β Pic b is challenging for planet formation models.Interestingly, it fits in the mass-C/O relation obtained by Hoch et al. (2022).They found a clear threshold at 4 M J beyond which imaged planet mostly have C/O consistent with solar ∼0.55, while below 4 M J , transiting planet have C/O with any values from 0.2 to 2.0.This threshold has been interpreted as distinguishing two main formation pathways for planets either less massive than 4 M J and or either more massive than 4 M J .This 4-M J threshold was already discovered in the distribution of stellar metallicities among massive Giant exoplanets with a similar interpretation (Santos et al. 2017).Therefore, the solar C/O of β Pic b might be the sign of a formation scenario that differ from the usual core accretion scenario invoked for close-in less massive planets, either still by core accretion but considering perhaps a distinct pathway for planets formed at large separation that did not undergo inward migration, or either more simply by gravitational instability.

Conclusion
In this study we have derived a new infrared spectrum of the young giant planet β Pic b observed with SINFONI.Doing so, we have shown that the actual spectral resolving power of SIN- FONI in the K-band at a spaxel-resolution of 12.5×25 mas 2 is ∼4000.Then, using a novel method of stellar halo removal, we have been able to directly extract the spectrum and the molecular lines of the planet without the need of using molecular mapping techniques.We have fitted the spectrum to models from the forward-modeling Exo-REM library.This led to different results, depending on assumptions on the planet mass and radius: without any prior constraints, we have obtained T eff =1555 +22 −29 K at a log g=3.12 +0.12 −0.09 , with sub-solar metallicity of -0.325 +0.065 −0.045 dex and a super-solar C/O=0.79 +0.01 −0.11 ; assuming a prior on the T eff =1724±15 K based on Chilcote et al. (2017) independent photometric characterisation, we have rather found higher T eff =1748 +3 −4 K, with again sub-solar metallicity of -0.235 +0.015 −0.013 dex and a now solar C/O=0.551±0.002.
Our preferred parameters are those derived imposing T eff =1724±15 K, as it better reflects the gravitational mass and the geometric radius derived independently using photometry.We find a projected rotation speed of β Pic b's equator of 25 +5 −6 km s −1 agreeing with the 25±3 km s −1 found by Snellen et al. (2014)     This tends to confirm the current estimation of the orbits and mass of the β Pic's planets, and add a new RV point for future dynamical characterisation of the system.
The presented work shows that ground-based infrared medium-resolution spectroscopy with even a modest resolving power of 4000 -close to that of the JWST/MIRI-MRS -without the need of doing molecular mapping and with a careful treatment of wavelength calibration, as well as star halo and telluric lines subtraction at any wavelength, can allow deriving key properties of imaged exoplanets, including equatorial rotation velocity.Proper determination of the atmospheric parameters though still requires independent priors, such as on T eff or log g, due to fit degeneracy of the unresolved spectral lines.

Fig. 1 .
Fig. 1.Absolute total (blue) and residual (orange) flux per spaxel.The spaxels are ordered from the brightest to the faintest in absolute total flux per spaxel.The absolute total flux is the raw flux obtained as output of the stack phase in Section 2.5.The residual flux is the remaining absolute flux once the star spectrum is removed (see Section 3.2 for more details).

Fig. 2 .
Fig. 2. Star spectrum calculated from the brightest spaxels.The flux is normalised to the pseudo-continuum flux at the top of the Br-γ line at 2.165 µm.

Fig. 3 .
Fig. 3. Summary of molecfit results.First panel: wavelength solutions showing the Doppler shift against the wavelength through the K-band with respect to cubes series number.Second panel: the measured resolving power from FWHM of the fitted LSF.Third and fourth panels: fitted ppmv abundances of H 2 O, CO 2 , CH 4 and N 2 O.

Fig. 4 .
Fig. 4. A stellar spectrum observed with SINFONI (blue) compared to the model Earth telluric spectrum calculated with molecfit (orange) and the continuum model (green).Upper panel: the two spectra directly compared.Lower panel: the stellar spectrum then divided by the telluric spectrum model.

Fig. 5 .
Fig.5.Steps of star spectrum removal methods compared.Q obs stands for the stacked cube observed, while {Q(i)  obs } refer to the full series of cubes collected during the night.

Fig. 7 .
Fig. 7. Star flux fraction K ip,⋆ at λ=2.165 µm (top panel) and the corresponding planet flux fraction 1-K ip,⋆ at the same wavelength (bottom panel).It can be seen that some background contributions participate to the residuals.

Fig. 8 .
Fig. 8. Top panel: the warp function δK(λ) (orange) normalised to 1 at the Brackett-γ wavelength, is compared to a normalised Exo-REM model with T =1750K, log g=4.0, [Fe/H]=0.0and C/O=0.55 (blue) and the normalised stellar spectrum η ip,⋆ (green) of β Pic on the K-band.Residual telluric absorptions can be seen in the stellar spectrum.Bottom panel: the naked Exo-REM model (blue) compared to the warped model (orange).

Fig. 9 .
Fig. 9. Background subtraction and telluric removal from a planetspaxel spectrum.Top: planet spaxel spectrum and average background on a distant ring; middle: η res,corr compared to its fitted telluric model η tell ; bottom: planet-spaxel spectrum corrected from background and tellurics.

Fig. 10 .
Fig. 10.Histogram of relative difference between the exo-REM models emittance and Stefan-Boltzman law.

Fig. 11 .
Fig. 11.Top: the SINFONI planet absorption spectrum of β Pic b along the K-band once corrected from star, background and telluric pollutions.Bottom: Exo-REM spectrum of simulated atmosphere of 1700 K, log(g)=4.0and [Fe/H]=0.0with R=6000.Individual contributions of species CO, H 2 O, CH 4 and NH 3 are represented from top to bottom with diverse colors.For all the spectra, their continuum were divided out as explained in Section 4.2.

Fig. 13 .
Fig.13.The planet PSF radial profile at a varying distance from 0 to 8 spaxels.The relative flux is obtained by dividing out the continuum and subtracting the star spectrum, as explained in Section 3.2.

Fig. 14 .
Fig. 14. Background removal.Top panel: average spectrum about the planet centroid (blue) compared to the average background spectrum on a ring around the planet (orange).Common features can be well distinguished.Bottom panel: planet spectrum corrected from the background.The features remaining are β Pic b's atmosphere absorption lines, while the main stellar feature (Brackett-γ) has been well removed.

Fig. 16 .
Fig. 16.Grid search CCF maps for T eff -log g (top panel) T eff -C/O (middle panel) and T eff -]Fe/H].The solid, dashed and dotted lines indicate 1, 10 and 20% ∆CCF levels compared to the optimal CCF.The red cross locates the optimum.
marizes the results of two different runs, one with all parameters freely varying within uninformed priors (except for R).The corner-plot showing the posterior distributions with all parameters freely varying is plotted on Fig. A.1.The synthetic spectrum at the median of these posterior distributions is compared to the observed data on Fig. A.2.The derived parameters are in good agreement with the initial guess from the grid-search in Section 6.1, and lead to more reliable confidence regions, with an effective temperature ranging at 1-σ within 1555 +29 −22 K, log g=3.12 +0.12 −0.09 , Fe/H=-0.325 +0.065 −0.045 dex, a super-solar C/O ∼0.79 +0.01 −0.11 and a v sin i=31±5 km s −1 .Our temperature estimate agrees with the GRAVITY Collaboration et al. (2020) results for the GRAVITY + GPI YJH band data fitted with exo-REM models (1590±20 K).Also our derived metallicity and Table 3. List of prior probability density functions (PDF) and parameters used for the MCMC runs.For T eff two types of priors are used, either uniform along the range of possible values.The MCMC run stops whenever the maximum autocorrelation length among all parameters, max(τ λ ), is larger than 50 times the actual number of iterations.Parameter Prior PDFs or value T eff U(400; 2000) or N(1724; 15) log g U(3.0; 5.0) or N(4.18; 0g agree with the fit of the GRAVITY-only spectrum.But we found a super-solar C/O where they found it sub-solar.The log g∼3.1 found above, as well as by the GRAVITY collaboration (GRAVITY Collaboration et al. 2020), implies a planet mass close to 2 M J .It disagrees with the independent dynamical constraints on the mass of the planet ∼12 M J (Snellen & Brown 2018; Lagrange et al. 2020; GRAVITY Collaboration et al. 2020).We thus try and fix the log g at a higher value, around 4.18, as suggested by Chilcote et al. (2017) work, defining a Gaussian prior probability distribution for the log g∼4.18±0.01.As shown in Table

Fig. 18 .
Fig. 18.Corner plot summarizing the MCMC results for the fit of the β Pic b IR SINFONI spectrum by Exo-REM models with T eff fixed at the Chilcote et al. (2017)'s value (see text).
at high spectral resolution with CRIRES and the most recent CRIRES+ Landman et al. (2023)'s v sin i=19.9±1.1 km s −1 .However, our measurement of a solar C/O is in stark contrast with the sub-solar C/O∼0.45 obtained from GRAVITY spectra (GRAVITY Collaboration et al. 2020)

Fig. 19 .
Fig. 19.Plot comparing in the top-panel the β Pic b SINFONI spectrum (blue) and the median Exo-REM model (orange) from the MCMC posteriors with T eff fixed at the Chilcote et al. (2017)'s value (see text).The bottom-panel shows the residuals.The uncertainties assumed for the observed flux are plotted as light-blue vertical lines.

Fig. 20 .
Fig. 20.The log-linear relationship between the equatorial velocity and mass of Solar System planets.In blue, possible laws compatible with planets Mars, Jupiter, Saturn, Uranus and Neptune.The ellipses show the v sin i and mass of β Pic b derived in this work (pink) and in S14 (cyan) and Landman et al. (2023) (red).

Fig. 22 .
Fig. 22.The v sin i prediction (gray region and black solid line) of the planet for different masses at 23 Myr, assuming a tilt of 0 • compared to edge-on.Green and blue area show the v sin i confidence regions on β Pic b based on this work's, S14's (blue) and Landman et al. (2023)'s (red) results.

Fig. 23 .
Fig. 23.Ice lines locations through time of H 2 O, CO 2 and CO derived for β Pictoris.The location of planet b is marked as a red circle at 23±3 Myr and prolonged down to 10 kyr with a black dotted line.

Table 4 .
Tablesummarizingthe MCMC results for the run.Unless stated otherwise, all parameters, except R spec , follow uniform priors.We describe the posterior distribution of any parameter by the median and the deviation of the 16th and 84th percentiles to the median.The maximum likelihood estimator (MLE) is given as well, at which point the goodness-of-fit diagnostics are determined.

Table 5 .
Ice lines of CO, CO 2 and H 2 O.The values of T evap are taken fromÖberg et al. (2011)and the corresponding ice lines radii are calculated as explained in the text.Pic that is solar up to a factor ∼2