11institutetext:  LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne
Université, Université Paris Cité, 5 place Jules Janssen, 92195 Meudon,
France††thanks: Please send any request to<EMAIL_ADDRESS>22institutetext: Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France
33institutetext: Laboratoire Lagrange, Université Cote d’Azur, CNRS,
Observatoire de la Cote d’Azur, 06304 Nice, France

# A new treatment of telluric and stellar features for medium resolution
spectroscopy and molecular mapping

Application to the abundance determination on $\beta$ Pic b

F. Kiefer 11 M. Bonnefoy 22 B. Charnay 11 A. Boccaletti 11 A.-M. Lagrange 1122
G. Chauvin 3322 B. Bézard 11 M. Mâlin 11

(Received – 05/06/2023 ; accepted – 11/02/2024)

Molecular mapping is a supervised method exploiting the spectral diversity of
integral field spectrographs to detect and characterize resolved exoplanets
blurred into the stellar halo. We present an evolution of the method to remove
the stellar halo and the nuisance of telluric features in the datacubes and
access a continuum-subtracted spectra of the planets at R$\sim$4000\. We
derive 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 $\beta$ Pictoris b. We recover the CO and H2O
detections in the atmosphere of $\beta$ Pic b using molecular mapping. We
further determine some basic properties of its atmosphere, with
$T_{\text{eq}}$=$1748^{+3}_{-4}$ K, a sub-solar
[Fe/H]=$-0.235^{+0.015}_{-0.013}$ dex, and a solar C/O=0.551$\pm$0.002 in
contrast with values measured for the same exoplanet with other infrared
instruments. We confirm a low projected equatorial velocity of
25${}^{+5}_{-6}$ km s-1. We are also able to measure, for the first time with
a medium-resolution spectrograph, the radial velocity of $\beta$ Pic b
relative to the central star at MJD=56910.38 with a km/s precision of
-11.3$\pm$1.1 km s-1, compatible with ephemerides based on the current
knowledge of the $\beta$ Pic system.

###### Key Words.:

Exoplanets ; Atmosphere ; Direct imaging ; High contrast imaging ; High
angular resolution ; Spectroscopy

##  1 Introduction

The system of $\beta$ 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 $\beta$ Pic
$\sim$23$\pm$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 $\beta$ 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 $\beta$
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, 2012).

Model and age dependent brightness–mass relationships predict $\beta$ Pic b
mass to be within 9-13 M${}_{\text{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) variations induced by the planets b & c. In particular,
the available RV data do not cover the whole $\beta$ Pic b period, the extrema
of the recently discovered $\beta$ 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. 2019,
2020; Vandal et al. 2020). Gaia was used by several authors to further
constrain the planet $b$ mass, $<$20 M${}_{\text{J}}$ (Bonnefoy et al. 2013),
13$\pm$3 M${}_{\text{J}}$ (Dupuy et al. 2019), 12.7$\pm$2.2 M${}_{\text{J}}$
(GRAVITY Collaboration et al. 2020), 10-11 M${}_{\text{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.5}^{+2.6}$ M${}_{\text{J}}$ using the Hipparcos–GaiaDR2 measurement
of astrometric acceleration (Brandt et al. 2021), and 11.7${}^{+2.3}_{-2.1}$
M${}_{\text{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. 2019, 2022), initially 2.54–$\sigma$ significant using
the DR2 (Kervella et al. 2019), became compatible with zero at 0.86–$\sigma$
using the DR3 (Kervella et al. 2022). This explains the difference in the
derived mass. From dynamical considerations, the mass of $\beta$ Pic b is thus
bounded within 9–15 M${}_{\text{J}}$.

The study of infrared spectra emitted by transiting and non-transiting hot
Jupiter (Brogi et al. 2012, 2013; de Kok et al. 2013; Birkby et al. 2013;
Lockwood et al. 2014; Brogi et al. 2014; Piskorz et al. 2016, 2017; Birkby et
al. 2017; Guilluy et al. 2019; Cont et al. 2021, 2022b; Yan et al. 2022; Cont
et al. 2022a) and 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, CO2, H2O, NH3 and CH4.
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. 2019, 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_{\text{eff}}$, $\log g$ and relative abundances.

For $\beta$ 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 H2O 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_{\text{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 $\beta$ 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$\pm$10 K with $\log
g$=4.35$\pm$0.09, a super-solar metallicity [M/H]$\sim$0.7$\pm$0.1 and a sub-
solar C/O=0.43$\pm$0.04111The commonly adopted solar abundance values are
$\log N_{M,\odot}/N_{H,\odot}$=$\log Z_{\odot}/X_{\odot}$=-1.74 and C/O=0.55
(Asplund et al. 2009).. 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 astrometric+RV planet mass estimation $\sim$12
M${}_{\text{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 $\beta$ 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$\sim$-0.4) and a
sub-solar C/O=0.41$\pm$0.04. Allowed by the high-resolution they obtained a
new $v\sin i$ measurement of 19.9$\pm$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 $\beta$ 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.

##  2 SINFONI data pre-processing

###  2.1 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 field-of-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 $\beta$
Pictoris surroundings performed with the 0.8″$\times$0.8″FoV subdivided into
64$\times$64 spaxels222Initially there are 64$\times$32 rectangular spaxels of
size 12.5$\times$25 mas2, but are then subdivided into 64$\times$64 square
spaxels by splitting one spaxel into two equal flux spaxels. of size
12.5$\times$12.5-mas2 at $R$=4000 along the $K$-band (2.08-2.45 $\mu$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
$\beta$ 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 $\beta$
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%.

###  2.2 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 high-contrast 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
”odd-even” 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 $\beta$ 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 $\sim$1 pixel in the final registered cubes of $\beta$
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
$\beta$ Pic b.

###  2.3 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 SINFONI 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-$\gamma$ line at $\sim$2.165 $\mu$m by fitting the
wings of the line by a 2-degree polynomial and retrieving the level of the
continuum at 2.165 $\mu$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 $\sim$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 $\mu$m, mainly H2O, CO2 and CH4
lines.

Figure 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). Figure
2: Star spectrum calculated from the brightest spaxels. The flux is normalised
to the pseudo-continuum flux at the top of the Br-$\gamma$ line at 2.165
$\mu$m.

###  2.4 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 ($\sim$4%), ambiant
pressure ($\sim$740 hPa), ambiant temperature ($\sim$12∘C), mirror temperature
($\sim$10.9∘C) and airmass ($\sec z$$\sim$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 $\chi^{2}$ are consistent
within 1.3–1.4 all through the cubes time series. We found that the molecules
H2O, CO2, CH4 and N2O dominate the model tellurics spectrum in the K-band,
whilst CO, NH3, O2 and O3 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 H2O, CO2, CH4 and N2O 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 $\mu$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 $\beta$ Pic stellar spectrum compared with the molecfit resulting model.
There is a residual time-dependent 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 $\sim$2.16-2.32 pixels at 2.27 $\mu$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$\pm$90.

The atmospheric molecular column density of H2O is relatively constant along
the observation sequence. We note that the CO2 abundance level of $\sim$500
ppmv retrieved is $\sim$1.4 times higher than the reference value ($\sim$370
ppmv) in the model Earth atmosphere, while N2O abundance $\sim$0.30 ppmv is
about 1.2 its reference level. At the same date and using the whole K-band
spectrum (including deep CO2 and N2O lines below 2.07 $\mu$m) of another
reference star, the CO2 and N2O abundances reach respectively $\sim$390 ppmv
and 0.23 ppmv (Smette, private comm.). Relying on weak lines only, our
determination of the CO2 and N2O 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 $\sim$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
documentation333https://www.eso.org/sci/facilities/paranal/
decommissioned/sinfoni/doc/VLT-MAN-ESO-14700-3517
-P103.1.pdf.

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.

Figure 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 H2O, CO2, CH4 and N2O. Figure 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.

###  2.5 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 ($\alpha^{*}$,$\delta$)-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.griddata` routine from `scipy` in order for all the
shifted–rotated cubes to share a common ($\alpha^{*}$,$\delta$)-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_{\text{obs}}$.

##  3 Star spectrum removal with STAREM

###  3.1 A summary of H18 method

Figure 5: Steps of star spectrum removal methods compared. $Q_{\text{obs}}$
stands for the stacked cube observed, while $\\{Q^{(i)}_{\text{obs}}\\}$ refer
to the full series of cubes collected during the night.

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:

$F_{i_{p}}(\lambda)=F_{i_{p},\star}(\lambda)+F_{i_{p},p}(\lambda)+F_{i_{p},B}(\lambda)$
(1)

Among which the star flux $F_{i_{p},\star}$, 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 $\eta$. Both the astrophysical
source and telluric lines contribute to this transmission spectrum. It can
also be expressed as $\eta$=$1-A$, with $A$ a positive comb of spectral lines
with the continuum equal to zero everywhere. We point out that $A_{\star}$ and
$A_{p}$ for respectively the star and the planet, are supposed to be spaxels-
independent since intrinsic to the respective sources. Continua $C_{p}$ and
$C_{\star}$ 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_{\text{obs}}$ an approximated star spectrum
$F_{i,\star,\text{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_{\text{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_{\text{star-free}}$. It results in a speckle-proxy cube noted
$Q_{\text{spkl}}$ in Fig. 5. By multiplying those modelled variations by the
star spectrum, they finally obtain the star cube $Q_{\star}$ with at each
spaxel $i$ a star spectrum $F_{i,\star,\text{approx}}$ that accounts for
wavelength-dependent spaxel-to-spaxel variations.

At any spaxel other than the planet spaxels, subtracting
$F_{i,\star,\text{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:

$\displaystyle F_{\star,\text{approx}}=$
$\displaystyle\left(C_{p}(\lambda)+C_{\star}(\lambda)\right)\eta_{\star}(\lambda)$
(2)

By slightly overestimating the star contribution, the subtraction of
$F_{\star,\text{approx}}$ rather leads to

$\displaystyle\Delta F(\lambda)=$ $\displaystyle
C_{p}(\lambda)\left(A_{\star}(\lambda)-A_{p}(\lambda)\right)$ (3)

A supplementary stellar contribution – including tellurics – to 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 CH4 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.

Figure 6: The raw $\beta$ 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-$\gamma$ 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.

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, $\Delta F$ is not strictly-speaking a pure planet
atmosphere spectrum. This issue might also explain the detection of an
Br-$\gamma$ 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.

###  3.2 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
$\eta$=$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_{\star,\text{norm}}$ and a normalised cube
$Q_{\text{obs,norm}}$. Recalling equation (1), the normalised transmission
spectrum components of $F=C\,\eta$ at one of the planet spaxels in
$Q_{\text{obs}}$ are

$\eta(\lambda)=\frac{C_{\star}(\lambda)\,\eta_{\star}(\lambda)+C_{p}(\lambda)\,\eta_{p}(\lambda)}{C(\lambda)}$
(4)

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_{\star}$ and $K_{p}$ as

$\displaystyle K_{\star}(\lambda)$ $\displaystyle=1-K_{p}(\lambda)$ (5)
$\displaystyle K_{p}(\lambda)$
$\displaystyle=\frac{C_{p}(\lambda)}{C(\lambda)}$

.

The transmission spectrum $\eta(\lambda)$ can be expressed more compactly as

$\displaystyle\eta(\lambda)$
$\displaystyle=K_{\star}(\lambda)\,\eta_{\star}(\lambda)+K_{p}(\lambda)\,\eta_{p}(\lambda).$
(6)

Wishing to subtract the star contribution $K_{\star}\,\eta_{\star}$ from this
spectrum, we thus need to estimate $K_{\star}(\lambda)$. 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_{\star}$=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_{\star}(\lambda)$ function. In the specific case of
the K-band spectrum of $\beta$ Pic, since the Brackett-$\gamma$ line at 2.165
$\mu$m is the only strong feature in this spectral band, we are only able to
derive an approximated constant star contribution $\sim$
$K_{\star}(\lambda_{0}$=$2.165\,\mu\text{m})$, or $K_{\star,2.165}$
hereafter. The derived values of $K_{\star,2.165}$ and its pendent residual
contribution, 1-$K_{\star}$, through the SINFONI field-of-view around $\beta$
Pic b are shown in Figure 7.

Figure 7: Star flux fraction $K_{i_{p},\star}$ at
$\lambda$=$2.165\,\mu\text{m}$ (top panel) and the corresponding planet flux
fraction 1-$K_{i_{p},\star}$ at the same wavelength (bottom panel). It can be
seen that some background contributions participate to the residuals.

Removing an approximated contribution $K_{\star,2.165}\,\eta_{\star}$ at each
spaxel of $Q_{\text{obs,norm}}$ leads to a residual cube $Q_{\text{res}}$. In
this cube, the residual spectrum obtained on a planet spaxel is

$\displaystyle\eta_{\text{res}}(\lambda)=$
$\displaystyle\left(K_{\star}(\lambda)-K_{\star,2.165}\right)\,\eta_{\star}(\lambda)+K_{p}(\lambda)\,\eta_{p}(\lambda)$
(7)

Using eq. 5 this could equivalently be written as:

$\displaystyle\eta_{\text{res}}(\lambda)=$
$\displaystyle\left(K_{p,2.165}-K_{p}(\lambda)\right)\,\eta_{\star}(\lambda)+K_{p}(\lambda)\,\eta_{p}(\lambda)$
(8)

This is a flat spectrum whose baseline level is
$K_{p,2.165}$=1-$K_{\star,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 $\mu$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-$\gamma$ – more generally from any fitted stellar line – the line
amplitudes are impacted by a residual star spectrum component with an
amplitude $\delta K(\lambda)$=$K_{p,2.165}$-$K_{p}(\lambda)$. This is a
footprint of the planet and the star continua within the flat spectrum.
Unfortunately, we cannot assume that $\delta K\eta_{\star}$$\approx$$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.

Figure 8: Top panel: the warp function $\delta K(\lambda)$ (orange) normalised
to 1 at the Brackett-$\gamma$ wavelength, is compared to a normalised Exo-REM
model with $T$=1750K, $\log g$=4.0, $[$Fe/H$]$=0.0 and C/O=0.55 (blue) and the
normalised stellar spectrum $\eta_{i_{p},\star}$ (green) of $\beta$ 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).

###  3.3 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_{\text{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 residual
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 $\eta_{\text{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$\pm$30.

At each spaxel $i$ of the residual cube, we fit out this telluric model
$\eta_{\text{model}}$ from $\eta_{i,\text{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.3444Spaxel-to-spaxel variations in the telluric spectrum
arise from small differences in the wavelength calibration from one spaxel to
the others, and from adding a non-zero spaxel-dependent background continuum
when fitting and subtracting the reference stellar (plus telluric) spectrum in
Section 3.2.

We make an estimation of the background contribution in a spaxel $i$ by taking
the median spectrum in a distant ring around this spaxel
$\eta_{i,\text{ring}}$. Since the SINFONI PSF has FWHM$\sim$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
$\eta_{i,\text{res,corr}}=\eta_{i,\text{res}}-\eta_{i,\text{ring}}$. Its
contribution level at $\sim$2.165 $\mu$m can be viewed in Fig. 7 on the
brighter areas away from the planet location. Fig. 1 also shows the
2.165-$\mu$m absolute background flux that is $\sim$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:

$\displaystyle\eta_{i,\text{res,corr}}=$
$\displaystyle\left(K_{i,p,2.165}+K_{i,b,2.165}-K_{i,p}-K_{i,b}\right)\eta_{\star}+K_{i,p}\eta_{p}+K_{i,b}\eta_{b}$
$\displaystyle-\left(K_{\text{ring},b,2.165}-K_{\text{ring},b}\right)\eta_{\star}-K_{\text{ring},b}\eta_{b}$
(9)

Assuming an equal contribution of background in the central and the
neighboring spaxels, i.e. $K_{i,b}=K_{\text{ring},b}$, this equation
simplifies to:

$\displaystyle\eta_{i,\text{res,corr}}=$
$\displaystyle\left(K_{i,p,2.165}-K_{i,p}\right)\eta_{\star}+K_{i,p}\eta_{p}$
(10)

We recover Eq. 8, and the previously implicit background is now explicitly
suppressed.

The final step is to fit $\eta_{\text{model}}$ the telluric spectrum model to
$\eta_{i,\text{res,corr}}$. We use a simple least-squares optimisation at each
spaxel, allowing the telluric model to vary in intensity with
$\eta_{i,\text{tell}}=a\,\eta_{\text{model}}+b$ with $a+b=c_{i}$ if $c_{i}$ is
the level of the residuals at spaxel $i$. Because $\eta_{\text{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 $\eta_{\text{model}}$ and the telluric
spectrum at spaxel $i$. We divide out $\eta_{i,\text{tell}}$ from
$\eta_{i,\text{res,corr}}$ leading to the telluric-free spectrum
$\eta_{i,\text{res, tell-free}}$.

A summary of background subtraction and telluric removal at the central planet
spaxel is shown in Fig. 9.

Figure 9: Background subtraction and telluric removal from a planet-spaxel
spectrum. Top: planet spaxel spectrum and average background on a distant
ring; middle: $\eta_{\text{res,corr}}$ compared to its fitted telluric model
$\eta_{\text{tell}}$; bottom: planet-spaxel spectrum corrected from background
and tellurics.

##  4 Molecular mapping dedicated to SINFONI spectra

###  4.1 Exo-REM grid of models

To calculate the molecular mapping of diverse species in the atmosphere of
$\beta$ 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 H2O, CH4, CO, CO2, NH3, PH3, TiO,
VO, FeH, K and Na, and collision-induced absorptions of H2–H2 and H2–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${}_{\text{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
$\sigma\,T_{\text{eff}}^{4}$ Stefan-Boltzman law. This deviation peaks below
$\sim$5% (Fig. 10). We thus removed 3315 among 9573 spectra (i.e. 35%) with a
deviation larger than 5%.

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

Each synthetic spectrum produced by Exo-REM results from individual
contributions (or spectrum) of the species out of which it is composed, mainly
CO, H2O, CH4 and NH3, which we can use as individual molecular templates.

###  4.2 The cross-correlation function of the planet spectrum

At every spaxel, between 2.08 and 2.43 $\mu$m, we calculate the CCF of the
observed spectrum and an Exo-REM synthetic spectrum for an assumed
$T_{\text{eff}}$=1700 K, $\log g$=4.0 cgs and $[$Fe/H$]$=0.0 dex, as based on
Table 3 in GRAVITY Collaboration et al. (2020). The expected abundances of NH3
and CH4 in the atmosphere $\sim$ $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 $\beta$ Pic b’s atmosphere. Details on the contributions of H2O, CO, NH3
and CH4 are shown in Fig. 11 and compared to the $\beta$ Pic b’s spectrum
derived in next Section 5.

Figure 11: Top: the SINFONI planet absorption spectrum of $\beta$ 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.0 and
$[$Fe/H$]$=0.0 with $R$=6000. Individual contributions of species CO, H2O, CH4
and NH3 are represented from top to bottom with diverse colors. For all the
spectra, their continuum were divided out as explained in Section 4.2.

We exclude the red edge of the K-band beyond 2.43 $\mu$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 cross-correlating 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.

COH2OCH4$10^{-4}$NH3$10^{-4}$ Figure 12: Molecular CCF maps of $\beta$ Pic b
for species CO, H2O, CH4 (10-4) and NH3 (10-4). At each spaxel, the maps show
the CCF height at 0 km s-1, a value close to the expected radial velocity of
the planet $b$, $\left|\text{VR}\right|$$<$10 km s-1given the instrument LSF
of $\sim$75 km s-1.

We then fit the PSF of the planet CCF halo by a 3D Gaussian with respect to
$(\delta,\alpha,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, H2O and using the full spectrum, confirming, as in H18,
the presence of CO and H2O but not detecting NH3 and CH4. With some
variability from one model to the other, the CCF peak is located at a
separation to the central star of $\sim$351$\pm$5 mas, with a PSF broadening
standard deviation of $\sim$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 pix2 by 2000 km s-1 around the planet
($\delta$,$\alpha$,$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 spaxel-to-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${}_{\text{all}}$=19.6, S/N${}_{\text{CO}}$=11.8,
S/N${}_{\text{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, H2O
gives significantly better results than CO. This last properties is best
explained by the presence of less prominent but more numerous lines of H2O
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.2 here; H2O: 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 H2O 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.

Table 1:  Comparing results of the different star spectrum subtraction schemes
on the normalised CCF of the residual map with a 1700 K, $\log g$=4.0,
$[$Fe/H$]$=0.0 Exo-REM model, detailing the CO and H2O individual
contributions. We also tried the molecular mapping on CH4, NH3, FeH, and CO2,
none leading to a significant detection within the area around planet spaxels.

Method | CCF strength | Radial velocity | Planet PSF
---|---|---|---
| peak | noise | S/N | $v_{r}$ | $\sigma_{\text{broad}}$ |  $X_{c}$††footnotemark: † |  $Y_{c}$††footnotemark: † | $\sigma_{\text{PSF}}$
| height | level |  | (km/s) | (km/s) | (mas) | (mas) | (pixel)
Full-model: 1700 K, $\log g$=4.0, $[$Fe/H$]$=0.0, $R_{\text{spec}}$=6000
H18 | 0.8952$\pm$0.0076 | 0.0511 | 17.52$\pm$0.15 | 1.03$\pm$0.35 | 46.80$\pm$0.38 | 218.0$\pm$0.25 | -277.13$\pm$0.25 | 2.34$\pm$0.01
STAREM | 0.8452$\pm$0.0060 | 0.0431 | 19.60$\pm$0.14 | -0.25$\pm$0.25 | 41.87$\pm$0.30 | 216.63$\pm$0.25 | -276.75$\pm$0.25 | 2.71$\pm$0.01
only CO
H18 | 0.4837$\pm$0.0049 | 0.0432 | 11.20$\pm$0.11 | -4.29$\pm$0.41 | 45.91$\pm$0.49 | 215.38$\pm$0.25 | -274.25$\pm$0.25 | 2.47$\pm$0.02
STAREM | 0.5094$\pm$0.0048 | 0.0431 | 11.81$\pm$0.11 | -3.24$\pm$0.35 | 41.85$\pm$0.42 | 214.00$\pm$0.25 | -273.50$\pm$0.25 | 2.66$\pm$0.02
only H2O
H18 | 0.7887$\pm$0.0076 | 0.0523 | 15.07$\pm$0.15 | 2.28$\pm$0.38 | 41.64$\pm$0.33 | 218.62$\pm$0.25 | -278.63$\pm$0.25 | 2.28$\pm$0.02
STAREM | 0.6947$\pm$0.0056 | 0.0451 | 15.41$\pm$0.12 | 0.76$\pm$0.29 | 41.64$\pm$0.33 | 217.13$\pm$0.25 | -279.00$\pm$0.25 | 2.75$\pm$0.02

555
$${}^{\dagger}$$$${}^{\dagger}$$footnotetext: The star is located at coords
$(X_{\star},Y_{\star})$=(0,0) mas.

##  5 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 $\mu$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 $\mu$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 $\sim$5
spaxels.

Figure 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.

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\sigma_{\text{PSF}}^{2}}$
where $r$=$\|{\bf r_{i}}-{\bf r_{c}}\|$ is the spaxel-distance between a given
spaxel $i$ and the planet center position, and $\sigma_{\text{PSF}}$=2 spaxels
at 2.165 $\mu$m. We use this profile as the weight function with thus
$w(r)$=$e^{-r^{2}/2\sigma_{\text{PSF}}^{2}}$. 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 $\sqrt{N_{\text{spaxels}}}$. We found
reasonable to integrate the flux up to 4–$\sigma$ from the PSF centroid, with
$\sigma(\lambda)$=$\sigma_{\text{PSF}}\times 2.165\,\mu$m/$\lambda$.

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 $\pm$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.

Figure 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
$\beta$ Pic b’s atmosphere absorption lines, while the main stellar feature
(Brackett-$\gamma$) has been well removed.

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.0 dex in
Fig. 11.

##  6 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.

###  6.1 Grid search

We first explored the template space by a grid search. We tested the
$\chi^{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 profile with $v\sin i$=25 km s-1. The
rotational profile is the usual bell-like profile with limb-darkening
coefficient $\varepsilon$=0.6 (Gray 1997). Each model spectrum is then
flattened with the same median filtering function, with a window width of 0.01
$\mu$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)\,\eta_{\star}$, with $\eta_{\star}$ 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 $\mu$m. For
the ZNCCF the median of both model and data is subtracted before cross-
correlation. For the $\chi^{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
$\sigma_{\text{data}}$$\propto$$\sqrt{F_{\star,\text{approx}}}/C$, with
$F_{\star,\text{approx}}(\lambda)$ the reference stellar spectrum defined in
Section 2.3 and $C(\lambda)$ 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 $\mu$m from the standard deviation of flux
to the median on a width of 0.003 $\mu$m ; we take the average of the errors
at all spectral channels from 2.15 to 2.18 $\mu$m.

The grid search maps with $\chi^{2}$-score and CCF-score are shown on
respectively Fig. 15 and 16 and the results are summarized in Table 2. For the
$\chi^{2}$ the confidence intervals are bounded by the $\Delta\chi^{2}$ with
the 1, 2 & 3–$\sigma$ regions corresponding to $\Delta\chi^{2}$=2.3, 6.2 &
11.8. For individual measurements in Table 2 the 1 and 2–$\sigma$ intervals
correspond to $\Delta\chi^{2}$=1 & 4, interpolated from the 1D
$\Delta\chi^{2}$ obtained for each parameter. The small extent of the
1-$\sigma$ 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-$\sigma$ intervals are certainly more reliable than the 1-$\sigma$ intervals
as error bars. For the CCF grid search we show the $\Delta$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 $\chi^{2}$ grid search leads to a solution with $T_{\text{eff}}$=1550 K,
$[$Fe/H$]$=0.0 dex, $\log(g)$$\sim$$3.5$, and C/O$\sim$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.

Figure 15: Grid search $\Delta\chi^{2}$ maps for T${}_{\text{eff}}$-$\log g$
(top panel) T${}_{\text{eff}}$-C/O (middle panel) and
T${}_{\text{eff}}$-$]$Fe/H$]$. The solid, dashed and dotted lines indicate 1,
2 and 3-$\sigma$ confidence regions. The red cross locates the optimum. Figure
16: Grid search CCF maps for T${}_{\text{eff}}$-$\log g$ (top panel)
T${}_{\text{eff}}$-C/O (middle panel) and T${}_{\text{eff}}$-$]$Fe/H$]$. The
solid, dashed and dotted lines indicate 1, 10 and 20% $\Delta$CCF levels
compared to the optimal CCF. The red cross locates the optimum.

Figure 17: Best matching models, using the $\chi^{2}$ optimisation (top) and
the CCF optimisation (bottom). Table 2: $\chi^{2}$ and CCF grid search results
(see Section 6.1. The reduced optimum $\chi^{2}$ has been normalised to 1 and
is therefore not given here. Uncertainties are obtained by interpolation of
the 1D $\Delta\chi^{2}$ maps.

Parameter | Unit |  $\chi^{2}$ grid | CCF grid
---|---|---|---
|  | optimum | 1-$\sigma$ | 2-$\sigma$ | optimum
$T_{\text{eff}}$ | (K) | 1500 | (+3)(-1) | (+100)(-5) | 2000
$\log g$ | (cgs) | 3.5 | (+0.0)(-0.3) | (+0.0)(-0.5) | 3.5
$[$Fe/H$]$ | (dex) | 0.0 | (+0.0)(-0.3) | (+0.0)(-0.5) | 0.0
C/O |  | 0.70 | (+0.01)(-0.00) | (+0.10)(-0.01) | 0.65
††footnotemark: †$v_{r}$ | (km s-1) | 0.78 |  $\pm$0.79 |  | -1.10
††footnotemark: †$a$ |  | 0.999 |  $\pm$0.025 |  | fixed to 1
††footnotemark: †$v\sin i$ | (km s-1) | fixed to 25 |  |  | fixed to 25
R |  | fixed to 4020 |  |  | fixed to 4020
$\chi^{2}_{r}$ |  | 1 |  |  |
max(CCF) |  |  |  |  | 0.69
666
$$\dagger$$$$\dagger$$footnotetext: The 1-$\sigma$ 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.

###  6.2 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_{\text{eff}}$, $\log g$, $[$Fe/H$]$, C/O, $v_{r}$ and $v\sin i$. We also
fit for the spectral resolving power $R_{\text{spec}}$, in order to account
for possible deviations from $R$=4020$\pm$30\. We also sample a rescaling $a$
to adjust the continuum of model and observed spectra together. The pre-
processings of data and models were the same as those used adopted in the
gridsearch. The error bars of the data points $\sigma_{\text{data}}$ (as
defined in Section 6.1) might be over or under estimated by a certain
multiplicative amount $f_{\text{err}}$ and under-estimated by a constant
jitter $\sigma$. We thus use the error model
$\sigma_{\text{err}}^{2}=(f_{\text{err}}\sigma_{\text{data}})^{2}+\sigma^{2}$
with $\sigma_{\text{data}}$ depending on wavelength. We add $f_{\text{err}}$
and $\sigma$ as hyper-parameters in the MCMC. The prior distributions, either
normal ${\mathcal{N}}$ or uniform ${\mathcal{U}}$, assumed for all parameters
are listed in Table 3. All parameters, except $R_{\text{spec}}$, follow
uniform priors.

Table 3: List of prior probability density functions (PDF) and parameters used for the MCMC runs. For $T_{\text{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(\tau_{\lambda})$, is larger than 50 times the actual number of iterations. Parameter | Prior PDFs or value
---|---
$T_{\text{eff}}$ |  ${\mathcal{U}}(400;2000)$ or ${\mathcal{N}}(1724;15)$
$\log g$ |  ${\mathcal{U}}(3.0;5.0)$ or ${\mathcal{N}}(4.18;0.01)$
$[$Fe/H$]$ | ${\mathcal{U}}(-0.5;1.0)$
C/O | ${\mathcal{U}}(0.1;0.8)$
$v\sin i$ | ${\mathcal{U}}(0,+\infty)$
$R_{\text{spec}}$ | ${\mathcal{N}}(4020;30)$
$v_{r}$ | ${\mathcal{U}}(-\infty,+\infty)$
$f_{\text{err}}$ | ${\mathcal{U}}(0.1,10)$
$\sigma$ | ${\mathcal{U}}(0,0.5)$
$a$ | ${\mathcal{U}}(0.5,1.5)$
$N_{\text{steps, max}}$ | 1,000,000
$N_{\text{walkers}}$ | 20

The Exo-REM model spectra are interpolated through the initial 4D-grid at the
specific values of parameters chosen by the 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$\times$ the largest auto-correlation
length of the samples. Table 4 summarizes 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. 24. The synthetic spectrum at the median of
these posterior distributions is compared to the observed data on Fig. 25. 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–$\sigma$ 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 $\sim$0.79${}^{+0.01}_{-0.11}$ and a $v\sin i$=31$\pm$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$\pm$20 K). Also our derived metallicity and $\log g$ 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$$\sim$3.1 found above, as well as by the GRAVITY collaboration
(GRAVITY Collaboration et al. 2020), implies a planet mass close to 2
M${}_{\text{J}}$. It disagrees with the independent dynamical constraints on
the mass of the planet $\sim$12 M${}_{\text{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$$\sim$4.18$\pm$0.01. As shown in Table 4, in this case, the MCMC leads to a
larger temperature of 1746${}_{-3}^{+4}$ K, and a smaller C/O of
0.551$\pm$0.002 compatible with a solar value.

We also tried fixing $T_{\rm eff}$’s prior to Chilcote et al. (2017)’s value
1724$\pm$15 K and both $\log g$ and $T_{\rm eff}$. These trials are also shown
in Table 4 and are compatible with the fixed-$\log g$ trial. We calculated the
Akaike information criterion ${\rm AIC}$=$2k-2\ln{\mathcal{L}}$ for all
models. The maximum likelihood estimator (MLE) model obtained with the $T_{\rm
eff}$ fixed minimizes the Akaike information criterion (AIC) and is thus
preferred over all other models. The large difference in AIC of
$\Delta$AIC$\sim$36 means that the $T_{\rm eff}$–fixed model is
2.5$\times$10-8 times more likely than the model with uninformed priors, and
1.4$\times$ more likely than the model with both $\log g$ and $T_{\rm eff}$
fixed. This preferred solution has a $T_{\rm eff}$ of 1748${}_{-4}^{+3}$ K, a
$\log g$=4.22$\pm$0.03 slightly 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
$\sim$0.551$\pm$0.002, and a $v\sin i$=25${}^{+5}_{-6}$ km s-1. In all cases,
we found a jitter $\sigma$ that is compatible with zero, with a correcting
factor of error bars $f_{\rm err}$$\sim$1${}^{+0.05}_{-0.11}$ very close to
$1$. Both imply that our determination of flux error bars
$\sigma_{\text{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$\pm$0.9 km 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 $\beta$ Pic b of 8.7$\pm$0.9 km s-1.
Compared to $\beta$ Pic systemic RV of $\sim$20$\pm$0.7 km s-1 (Gontcharov
2006), it implies an RV of the planet relative to $\beta$ Pic’s central star
of -11.3$\pm$1.1 km s-1 at MJD=56910.38. This value agrees at 0.3–$\sigma$
with the predicted RV of the planet at this MJD extrapolating the ephemerides
of $\beta$ 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$\pm$1.7 km s-1 at
MJD=56,644.5 (Snellen et al. 2014).

Table 4:  Table summarizing the MCMC results for the run. Unless stated
otherwise, all parameters, except $R_{\text{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.

Parameter | Unit | All priors uniform | Prior $\log g$=4.18$\pm$0.01 | Prior $T_{\rm eff}$=1724$\pm$15 K | Prior $\log g$ & $T_{\rm eff}$
---|---|---|---|---|---
|  | median$\pm$$\sigma_{68.3\%}$ | MLE | median$\pm$$\sigma_{68.3\%}$ | MLE | median$\pm$$\sigma_{68.3\%}$ | MLE | median$\pm$$\sigma_{68.3\%}$ | MLE
$T_{\text{eff}}$ | K | 1555${}_{-29}^{+22}$ | 1525 | 1746${}_{-3}^{+4}$ | 1746 | 1748${}_{-4}^{+3}$ | 1749 | 1745${}_{-2}^{+3}$ | 1746
$\log g$ | cgs | 3.12${}_{-0.09}^{+0.12}$ | 3.24 | 4.185${}_{-0.010}^{+0.010}$ | 4.182 | 4.216${}_{-0.031}^{+0.027}$ | 4.210 | 4.183${}_{-0.010}^{+0.010}$ | 4.190
$[$Fe/H$]$ | dex | -0.325${}_{-0.045}^{+0.065}$ | -0.246 | -0.235${}_{-0.011}^{+0.014}$ | -0.226 | -0.235${}_{-0.013}^{+0.015}$ | -0.223 | -0.237${}_{-0.010}^{+0.013}$ | -0.233
C/O | — | 0.79${}_{-0.11}^{+0.01}$ | 0.68 | 0.551${}_{-0.002}^{+0.002}$ | 0.550 | 0.551${}_{-0.002}^{+0.002}$ | 0.550 | 0.551${}_{-0.002}^{+0.002}$ | 0.550
$v\sin i$ | km s-1 | 31.4${}_{-5.0}^{+4.6}$ | 33.0 | 27.1${}_{-5.4}^{+4.5}$ | 24.9 | 24.9${}_{-6.2}^{+5.2}$ | 25.4 | 26.6${}_{-5.2}^{+4.4}$ | 28.6
$R_{\text{spec}}$ | — | 4018${}_{-30}^{+30}$ | 4042 | 4019${}_{-30}^{+30}$ | 4024 | 4018${}_{-30}^{+30}$ | 4032 | 4019${}_{-30}^{+30}$ | 4020
$v_{r}$ | km s-1 | 0.77${}_{-0.97}^{+1.05}$ | 1.20 | 0.58${}_{-0.89}^{+0.86}$ | 0.30 | 0.62${}_{-0.89}^{+0.87}$ | 1.03 | 0.63${}_{-0.88}^{+0.84}$ | 0.28
$f_{\text{err}}$ | — | 1.02${}_{-0.13}^{+0.05}$ | 1.08 | 0.98${}_{-0.12}^{+0.05}$ | 0.98 | 0.99${}_{-0.11}^{+0.05}$ | 0.99 | 0.98${}_{-0.12}^{+0.05}$ | 1.00
⋆⋆footnotemark: ⋆$\sigma$ | — | $<0.013$ | 0.002 | $<0.013$ | 0.007 | $<0.012$ | 0.006 | $<0.013$ | 0.005
$a$ | — | 1.000${}_{-0.001}^{+0.001}$ | 1.000 | 1.000${}_{-0.001}^{+0.001}$ | 1.000 | 1.000${}_{-0.001}^{+0.001}$ | 1.000 | 1.000${}_{-0.001}^{+0.001}$ | 1.000
Goodness-of-fit at MLE |  |  |  |  |  |  |  |  |
$\chi^{2}$ |  |  | 1345 |  | 1375.1 |  | 1366.0 |  | 1340.5
$\chi^{2}_{\text{red}}$ |  |  | 1.002 |  | 1.029 |  | 1.022 |  | 1.003
$\log{\mathcal{L}}/N_{\text{DoF}}$ |  |  | 3.274 |  | 3.284 |  | 3.285 |  | 3.281
AIC |  |  | -8721.6 |  | -8755.6 |  | -8758.3 |  | -8757.6
Diagnostics |  |  |  |  |  |  |  |  |
††footnotemark: †$N_{\text{step}}$/max($\tau_{\lambda}$) |  | 17 | 52 | 50 | 51
Acceptance rate |  | 0.21 | 0.31 | 0.31 | 0.31
$N_{\text{DoF}}$ |  | 1335 | 1336 | 1336 | 1337

777
$$\star$$$$\star$$footnotetext: For $\sigma$, given the shape of the posterior
distribution peaking close to 0 and compatible with 0 at less than 2–$\sigma$,
we only give the upper-limit at the 84th percentile.
$$\dagger$$$$\dagger$$footnotetext: max($\tau_{\lambda}$) is the maximum auto-
correlation time among all varied parameters $\lambda$.

Figure 18:  Corner plot summarizing the MCMC results for the fit of the
$\beta$ Pic b IR SINFONI spectrum by Exo-REM models with $T_{\rm eff}$ fixed
at the Chilcote et al. (2017)’s value (see text). Figure 19:  Plot comparing
in the top-panel the $\beta$ Pic b SINFONI spectrum (blue) and the median Exo-
REM model (orange) from the MCMC posteriors with $T_{\rm 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.

##  7 Discussion

###  7.1 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$\pm$3 km s-1 and the most recent
Landman et al. (2023)’s $v\sin i$ estimation of 19.9$\pm$1.1 km s-1. The
larger confidence region of our measurement is explained by the much lower
resolving power of SINFONI ($R$=4030$\pm$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 mass-spin law followed by Solar
System objects from Jupiter down to asteroids (e.g. Hughes 2003), that would
imply a spin velocity of $\sim$50 km s-1. Figure 20 shows this relationship
for the Solar system planets, with a log-linear law well-fitted to
Earth+Moon888Due to major exchange of momentum between Earth and Moon, the
Earth+Moon couple is considered as if the Moon was put back into the Earth and
the total momentum conserved (Hughes 2003)., 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 $\beta$ Pic b compared to the expected
equatorial speed at a mass of $\sim$11 M${}_{\text{J}}$ is mostly explained by
the young age of the planet ($\sim$23 Myr; Mamajek & Bell 2014). Indeed, the
planet is currently contracting from $\sim$1.5 R${}_{\text{J}}$ down to
$\sim$1 R${}_{\text{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$\pm$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 $\beta$ Pic b and its mass, especially if the mass is
contained within 10-14 M${}_{\text{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${}_{\text{J}}$ or $>$14 M${}_{\text{J}}$. These masses are marginally
supported by the combined astrometric and RV measurements that rather favour a
planet mass within 9-15 M${}_{\text{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-$\sigma$ of the
predicted $v\sin i$ at 23 Myr for a mass at 11 M${}_{\text{J}}$ and our
measurement 25$\pm$6 km s-1 leads to a tilt compared to edge-on $>$15∘.

Figure 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 $\beta$ Pic b derived in this work (pink) and in S14 (cyan) and
Landman et al. (2023) (red). Figure 21: Illustration of evolutionary tracks of
equatorial radius through time with ATMO2020 models (Phillips et al. 2020) and
Baraffe-Chabrier-Barman (Baraffe et al. 2008) or BCB08 models at different
planet mass.

Figure 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^{\circ}$ compared
to edge-on. Green and blue area show the $v\sin i$ confidence regions on
$\beta$ Pic b based on this work’s, S14’s (blue) and Landman et al. (2023)’s
(red)

results.

###  7.2 A solar C/O

We derive here a value of the C/O for $\beta$ Pic b that is solar,
0.551$\pm$0.002, while GRAVITY Collaboration et al. (2020) and Landman et al.
(2023) found a sub-solar C/O of respectively 0.43$\pm$0.05 and 0.41$\pm$0.04.
Forming $\beta$ 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 H2O and CO2 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 H2O 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 $\beta$ 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 $\beta$
Pic’s planets.

In this framework, we can estimate the location of the ice lines compared to
$\beta$ 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):

$T=T_{0}\left(\frac{r}{1\,\text{R}_{\star}}\right)^{-0.62}$ (11)

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 $\beta$
Pic (Saffe et al. 2021), we have that $T_{0}$=8000 K. Considering the typical
evaporation temperatures of H2O, CO2 and CO summarized in Table 5 with
$R_{\star}$=1.7 R⊙ (Kervella et al. 2004), we derived the radii of the
different ice lines around $\beta$ Pic. They are shown in Table 5 as well. It
results that with $a_{b}$$\sim$9.8 au (Lagrange et al. 2020), planet $b$ is
located between the H2O (6$\pm$1 au) and the CO2 (31$\pm$5 au) ice lines. Fig.
23 shows the variation of the ice lines through time, as derived using
Dartmouth12 evolutionary tracks for pre main sequence stars.

Figure 23: Ice lines locations through time of H2O, CO2 and CO derived for $\beta$ Pictoris. The location of planet $b$ is marked as a red circle at 23$\pm$3 Myr and prolonged down to 10 kyr with a black dotted line. Table 5: Ice lines of CO, CO2 and H2O. The values of $T_{\text{evap}}$ are taken from Öberg et al. (2011) and the corresponding ice lines radii are calculated as explained in the text.  Molecule |  $T_{\text{evap}}$ (K) |  $r_{\text{ice}}$ (au)
---|---|---
CO | 20$\pm$2 | 125$\pm$20
CO2 | 47$\pm$5 | 31$\pm$5
H2O | 135$\pm$15 | 5.7$\pm$1.0

Following Petrus et al. (2021), our C/O measurement fits well to the Nissen
(2013) planet-C/O to star-Fe/H linear relation. They indeed report a positive
correlation between the C/O ratio and the [Fe/H] of stars with an average
C/O=0.58$\pm$0.06 at [Fe/H]=0. Elemental abundances and spectral matching
agree on a metallicity of $\beta$ Pic that is solar up to a factor $\sim$2
(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 $\sim$10 au, especially if the mass of $\beta$ Pic b is
larger than 10 M${}_{\text{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 $\beta$ 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 H2O 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 $\beta$ Pic b provided most of its gas was not
accreted before 1 Myr.

To conclude, a solar C/O for $\beta$ 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${}_{\text{J}}$ beyond which
imaged planet mostly have C/O consistent with solar $\sim$0.55, while below 4
M${}_{\text{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 MJ and or either more massive
than 4 MJ. This 4-M${}_{\text{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
$\beta$ 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.

##  8 Conclusion

In this study we have derived a new infrared spectrum of the young giant
planet $\beta$ Pic b observed with SINFONI. Doing so, we have shown that the
actual spectral resolving power of SINFONI in the K-band at a spaxel-
resolution of 12.5$\times$25 mas2 is $\sim$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_{\text{eff}}$=1555${}^{+22}_{-29}$ K at a $\log
g$=3.12${}_{-0.09}^{+0.12}$, 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_{\rm eff}$=1724$\pm$15 K based on Chilcote et al.
(2017) independent photometric characterisation, we have rather found higher
$T_{\text{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$\pm$0.002.

Our preferred parameters are those derived imposing $T_{\rm eff}$=1724$\pm$15
K, as it better reflects the gravitational mass and the geometric radius
derived independently using photometry. We find a projected rotation speed of
$\beta$ Pic b’s equator of 25${}^{+5}_{-6}$ km s-1 agreeing with the 25$\pm$3
km s-1 found by Snellen et al. (2014) at high spectral resolution with CRIRES
and the most recent CRIRES+ Landman et al. (2023)’s $v\sin i$=19.9$\pm$1.1 km
s-1. However, our measurement of a solar C/O is in stark contrast with the
sub-solar C/O$\sim$0.45 obtained from GRAVITY spectra (GRAVITY Collaboration
et al. 2020) and by atmospheric retrieval of CRIRES spectra in Landman et al.
(2023).

Conversely to the conclusions of GRAVITY Collaboration et al. (2020), with a
C/O for the star also close to 0.55, the stellar C/O disfavours the scenario
of icy-planetesimals injection. It would indeed imply fine-tuning of the
planetesimals accretion flux and of the duration of the phenomenon so that it
reaches a value close to the stellar C/O. Such value of $\beta$ Pic b’s C/O is
on the other hand in better agreement with in situ gravitational instability
scenario, or formation next to the H2O ice line followed by migration. This
latter scenario could be supported by the proximity of the planets b & c
orbital periods to a 7:1 resonance, whose orbital parameters may have
significantly evolved during the first Myrs of the system.

Finally, we measure a radial velocity of $\beta$ Pic b relative to the central
star of -11.3$\pm$1.1 km s-1 at MJD=56,910.38. It agrees at 0.3–$\sigma$ with
the ephemerides for the orbit of this planet based on the current knowledge of
the system from (Lacour et al. 2021). This tends to confirm the current
estimation of the orbits and mass of the $\beta$ 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_{\rm eff}$ or $\log g$, due to fit degeneracy of the unresolved spectral
lines.

###### Acknowledgements.

We are thankful to Alain Smette for his careful and constructive reading of
the manuscript that helped improving the analysis significantly. This work was
based on observations collected at the European Southern Observatory under ESO
programme 093.C-0626(A). It was granted access to the HPC resources of MesoPSL
financed by the Region Ile de France and the project EquipMeso (reference
ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the
Agence Nationale pour la Recherche. This project has received funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (COBREX; grant agreement n° 885593). This
work was also funded by the initiative de recherches interdisciplinaires et
stratégiques (IRIS) of Université PSL ”Origines et Conditions d’Apparition de
la Vie (OCAV)”. FK also acknowledge funding from the Action Incitative
Exoplanètes de l’Observatoire de Paris.

## References

  * Andrews & Williams (2005) Andrews, S. M. & Williams, J. P. 2005, ApJ, 631, 1134
  * Andrews & Williams (2007) Andrews, S. M. & Williams, J. P. 2007, ApJ, 659, 705
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  * Baraffe et al. (2008) Baraffe, I., Chabrier, G., & Barman, T. 2008, A&A, 482, 315
  * Baudino et al. (2015) Baudino, J. L., Bézard, B., Boccaletti, A., et al. 2015, A&A, 582, A83
  * Beust et al. (1990) Beust, H., Lagrange-Henri, A. M., Vidal-Madjar, A., & Ferlet, R. 1990, A&A, 236, 202
  * Birkby et al. (2013) Birkby, J. L., de Kok, R. J., Brogi, M., et al. 2013, MNRAS, 436, L35
  * Birkby et al. (2017) Birkby, J. L., de Kok, R. J., Brogi, M., Schwarz, H., & Snellen, I. A. G. 2017, AJ, 153, 138
  * Blain et al. (2021) Blain, D., Charnay, B., & Bézard, B. 2021, A&A, 646, A15
  * Bonnefoy et al. (2013) Bonnefoy, M., Boccaletti, A., Lagrange, A. M., et al. 2013, A&A, 555, A107
  * Bonnet et al. (2004) Bonnet, H., Abuter, R., Baker, A., et al. 2004, The Messenger, 117, 17
  * Boss (1997) Boss, A. P. 1997, Science, 276, 1836
  * Brandt et al. (2021) Brandt, G. M., Brandt, T. D., Dupuy, T. J., Li, Y., & Michalik, D. 2021, AJ, 161, 179
  * Brogi et al. (2014) Brogi, M., de Kok, R. J., Birkby, J. L., Schwarz, H., & Snellen, I. A. G. 2014, A&A, 565, A124
  * Brogi et al. (2018) Brogi, M., Giacobbe, P., Guilluy, G., et al. 2018, A&A, 615, A16
  * Brogi et al. (2012) Brogi, M., Snellen, I. A. G., de Kok, R. J., et al. 2012, Nature, 486, 502
  * Brogi et al. (2013) Brogi, M., Snellen, I. A. G., de Kok, R. J., et al. 2013, ApJ, 767, 27
  * Burns (1975) Burns, J. A. 1975, Icarus, 25, 545
  * Charnay et al. (2018) Charnay, B., Bézard, B., Baudino, J. L., et al. 2018, ApJ, 854, 172
  * Chilcote et al. (2017) Chilcote, J., Pueyo, L., De Rosa, R. J., et al. 2017, AJ, 153, 182
  * Cont et al. (2021) Cont, D., Yan, F., Reiners, A., et al. 2021, A&A, 651, A33
  * Cont et al. (2022a) Cont, D., Yan, F., Reiners, A., et al. 2022a, A&A, 668, A53
  * Cont et al. (2022b) Cont, D., Yan, F., Reiners, A., et al. 2022b, A&A, 657, L2
  * Cugno et al. (2021) Cugno, G., Patapis, P., Stolker, T., et al. 2021, A&A, 653, A12
  * de Kok et al. (2013) de Kok, R. J., Brogi, M., Snellen, I. A. G., et al. 2013, A&A, 554, A82
  * Demars et al. (2023) Demars, D., Bonnefoy, M., Dougados, C., et al. 2023, arXiv e-prints, arXiv:2305.09460
  * Dupuy et al. (2019) Dupuy, T. J., Brandt, T. D., Kratter, K. M., & Bowler, B. P. 2019, ApJ, 871, L4
  * Eisenhauer et al. (2003) Eisenhauer, F., Abuter, R., Bickert, K., et al. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, ed. M. Iye & A. F. M. Moorwood, 1548–1561
  * Feng et al. (2022) Feng, F., Butler, R. P., Vogt, S. S., et al. 2022, ApJS, 262, 21
  * Ferlet et al. (1987) Ferlet, R., Hobbs, L. M., & Vidal-Madjar, A. 1987, A&A, 185, 267
  * Fish (1967) Fish, F. F. 1967, Icarus, 7, 251
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Gontcharov (2006) Gontcharov, G. A. 2006, Astronomical and Astrophysical Transactions, 25, 145
  * GRAVITY Collaboration et al. (2020) GRAVITY Collaboration, Nowak, M., Lacour, S., et al. 2020, A&A, 633, A110
  * Guilluy et al. (2019) Guilluy, G., Sozzetti, A., Brogi, M., et al. 2019, A&A, 625, A107
  * Hoch et al. (2022) Hoch, K. K. W., Konopacky, Q. M., Theissen, C. A., et al. 2022, arXiv e-prints, arXiv:2212.04557
  * Hoeijmakers et al. (2018) Hoeijmakers, H. J., Schwarz, H., Snellen, I. A. G., et al. 2018, A&A, 617, A144
  * Holweger et al. (1997) Holweger, H., Hempel, M., van Thiel, T., & Kaufer, A. 1997, A&A, 320, L49
  * Hughes (2003) Hughes, D. W. 2003, Planet. Space Sci., 51, 517
  * Kausch et al. (2015) Kausch, W., Noll, S., Smette, A., et al. 2015, A&A, 576, A78
  * Kervella et al. (2019) Kervella, P., Arenou, F., Mignard, F., & Thévenin, F. 2019, A&A, 623, A72
  * Kervella et al. (2022) Kervella, P., Arenou, F., & Thévenin, F. 2022, A&A, 657, A7
  * Kervella et al. (2004) Kervella, P., Thévenin, F., Morel, P., et al. 2004, in Stars as Suns : Activity, Evolution and Planets, ed. A. K. Dupree & A. O. Benz, Vol. 219, 80
  * Kiefer et al. (2014) Kiefer, F., Lecavelier des Etangs, A., Boissier, J., et al. 2014, Nature, 514, 462
  * Lacour et al. (2021) Lacour, S., Wang, J. J., Rodet, L., et al. 2021, A&A, 654, L2
  * Lagrange et al. (2012) Lagrange, A. M., Boccaletti, A., Milli, J., et al. 2012, A&A, 542, A40
  * Lagrange et al. (2010) Lagrange, A. M., Bonnefoy, M., Chauvin, G., et al. 2010, Science, 329, 57
  * Lagrange et al. (2009) Lagrange, A. M., Gratadour, D., Chauvin, G., et al. 2009, A&A, 493, L21
  * Lagrange et al. (2019) Lagrange, A. M., Meunier, N., Rubini, P., et al. 2019, Nature Astronomy, 3, 1135
  * Lagrange et al. (2020) Lagrange, A. M., Rubini, P., Nowak, M., et al. 2020, A&A, 642, A18
  * Landman et al. (2023) Landman, R., Stolker, T., Snellen, I., et al. 2023, arXiv e-prints, arXiv:2311.13527
  * Lecavelier des Etangs et al. (2022) Lecavelier des Etangs, A., Cros, L., Hébrard, G., et al. 2022, Scientific Reports, 12, 5855
  * Lockwood et al. (2014) Lockwood, A. C., Johnson, J. A., Bender, C. F., et al. 2014, ApJ, 783, L29
  * Mâlin et al. (2023) Mâlin, M., Boccaletti, A., Charnay, B., Kiefer, F., & Bézard, B. 2023, arXiv e-prints, arXiv:2301.02116
  * Mamajek & Bell (2014) Mamajek, E. E. & Bell, C. P. M. 2014, MNRAS, 445, 2169
  * Marois et al. (2006) Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556
  * Miles et al. (2023) Miles, B. E., Biller, B. A., Patapis, P., et al. 2023, ApJ, 946, L6
  * Mollière et al. (2022) Mollière, P., Molyarova, T., Bitsch, B., et al. 2022, ApJ, 934, 74
  * Morzinski et al. (2015) Morzinski, K. M., Males, J. R., Skemer, A. J., et al. 2015, ApJ, 815, 108
  * Nissen (2013) Nissen, P. E. 2013, A&A, 552, A73
  * Nowak et al. (2020) Nowak, M., Lacour, S., Lagrange, A. M., et al. 2020, A&A, 642, L2
  * Öberg et al. (2011) Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16
  * Palma-Bifani et al. (2023) Palma-Bifani, P., Chauvin, G., Bonnefoy, M., et al. 2023, A&A, 670, A90
  * Patapis et al. (2022) Patapis, P., Nasedkin, E., Cugno, G., et al. 2022, A&A, 658, A72
  * Petit dit de la Roche et al. (2018) Petit dit de la Roche, D. J. M., Hoeijmakers, H. J., & Snellen, I. A. G. 2018, A&A, 616, A146
  * Petrus et al. (2021) Petrus, S., Bonnefoy, M., Chauvin, G., et al. 2021, A&A, 648, A59
  * Petrus et al. (2023) Petrus, S., Chauvin, G., Bonnefoy, M., et al. 2023, A&A, 670, L9
  * Phillips et al. (2020) Phillips, M. W., Tremblin, P., Baraffe, I., et al. 2020, A&A, 637, A38
  * Piskorz et al. (2016) Piskorz, D., Benneke, B., Crockett, N. R., et al. 2016, ApJ, 832, 131
  * Piskorz et al. (2017) Piskorz, D., Benneke, B., Crockett, N. R., et al. 2017, AJ, 154, 78
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62
  * Ruffio et al. (2021) Ruffio, J.-B., Konopacky, Q. M., Barman, T., et al. 2021, AJ, 162, 290
  * Ruffio et al. (2019) Ruffio, J.-B., Macintosh, B., Konopacky, Q. M., et al. 2019, AJ, 158, 200
  * Saffe et al. (2021) Saffe, C., Miquelarena, P., Alacoria, J., et al. 2021, A&A, 647, A49
  * Santos et al. (2017) Santos, N. C., Adibekyan, V., Figueira, P., et al. 2017, A&A, 603, A30
  * Schlecker et al. (2021) Schlecker, M., Mordasini, C., Emsenhuber, A., et al. 2021, A&A, 656, A71
  * Schwarz et al. (2016) Schwarz, H., Ginski, C., de Kok, R. J., et al. 2016, A&A, 593, A74
  * Skilling (2006) Skilling, J. 2006, Bayesian Analysis, 1, 833
  * Smette et al. (2015) Smette, A., Sana, H., Noll, S., et al. 2015, A&A, 576, A77
  * Snellen et al. (2014) Snellen, I. A. G., Brandl, B. R., de Kok, R. J., et al. 2014, Nature, 509, 63
  * Snellen & Brown (2018) Snellen, I. A. G. & Brown, A. G. A. 2018, Nature Astronomy, 2, 883
  * Tamuz et al. (2005) Tamuz, O., Mazeh, T., & Zucker, S. 2005, MNRAS, 356, 1466
  * Vandal et al. (2020) Vandal, T., Rameau, J., & Doyon, R. 2020, AJ, 160, 243
  * Yan et al. (2022) Yan, F., Pallé, E., Reiners, A., et al. 2022, A&A, 661, L6

##  Appendix A MCMC run with all free parameters

Figure 24:  Corner plot summarizing the MCMC results for the fit of the
$\beta$ Pic b IR SINFONI spectrum by Exo-REM models with all free parameters.
The blue lines and dot show the median parameters ; the orange lines and dot
are those found for the preferred solution (see Section 6.2). Figure 25:  Plot
comparing the $\beta$ Pic b SINFONI spectrum (blue) and the median Exo-REM
model (orange) from the MCMC posteriors (Fig. 24). The lower-panel shows the
residuals. The uncertainties assumed for the observed flux are plotted as
light-blue vertical lines.


continues to pose a challenge for galaxy formation models.

##  Appendix B Computing Power Spectra

We quantify the kinetic state of the ISM in Fourier space by first calculating
the kinetic energy power spectra and then divide it by the mass within the box
to get the specific kinetic energy power spectra (see Section 4.3). This
improves the comparison between simulations. The kinetic energy power spectra
of the weight $w$ is defined as

$P\left(\mathbf{k}\right)=\tilde{w}\left(\mathbf{k}\right)\cdot
w\left(\mathbf{k}\right)^{*}\quad,$ (11)

where $\tilde{w}\left(\mathbf{k}\right)$ is the Fourier transform of the real
array $w\left(\mathbf{k}\right)$, $\mathbf{k}$ is the wave vector, and ‘*’
refers to the complex-conjugate. To get the Fourier transform of the kinetic
energy field we used $w=\sqrt{\rho}\mathbf{v}\left(\mathbf{r}\right)$. Here
$\mathbf{v}\left(\mathbf{r}\right)$ is the three dimensional velocity vector
at cell center $\mathbf{r}$ (measured in our case from the center of the
galaxy) and $\rho$ the total density within each simulation cell. With this
choice of weights, the power spectrum directly measures the kinetic energy
(see also Grisdale et al., 2017). Another weight that is often used
($w=\rho^{1/3}\mathbf{v}$) corresponds to a kinetic energy flux and has been
shown to reproduce energy spectra with Kolmogorov scaling
($e_{\mathrm{\mathbf{k}}}\propto\mathbf{k}^{-5/3}$) in super-sonic flows see
Kritsuk et al., 2007; Kowal & Lazarian, 2007; Schmidt et al., 2008; Federrath
et al., 2009, for this choice.

We calculate the power spectrum within a $\sim 10\times 10\times 10$ kpc3
region centred on the galaxy999More precisely, we calculate the size of the
box for the power spectrum such that a power of two number of cells covers the
region using the given cell resolution.. For all of the spectra we use a
uniform grid of resolution $\sim$ 40 pc and account for the non-periodic
boundary by zero-padding the calculated cubes. The length of the void domain
is in each dimension set to be twice as big as the chosen box size. We
experimented with more padding as well as different box sizes and confirmed
that the main results and findings are not altered.

To obtain the power spectrum $P(\mathbf{k})$ of the specific kinetic energy we
first Fourier transform the cube, bin $P(\mathbf{k})$ in wave vectors
$k=\left|\mathbf{k}\right|$, normalise by the area of the spherical surface
with radius k + 0.5, and finally divide it by the total mass within the box.
This finally gives the ‘angle-averaged’ specific kinetic energy spectrum Joung
& Mac Low, 2006b, e.g.. $\left<P(k)\right>$, where the physical scale $l$ is
connected to the wave vector via $k=2\pi/l$.