© The Authors 2025

1 Introduction

The emission of massive Jovian companions detected by highcontrast imaging can now be characterized in detail from ultraviolet (HST, Zhou et al. 2014; VLT/XSHOOTER, Petrus et al. 2020) to mid-infrared wavelengths (JWST/NIRSpec, 0.70–5.27 μm, Böker et al. 2022; JWST/MIRI, 4.9–27.9 μm, Wells et al. 2015; Argyriou et al. 2023). The mid-infrared is a rich spectral window for exoplanetary atmospheres, as has been hinted at by recent advances in atmospheric modelling: thanks to absorption features of molecules such as CO, CH₄, CO₂, or H₂O, it may encode information about cloud thickness (Charnay et al. 2018), vertical mixing and induced disequilibrium chemistry (Phillips et al. 2020), or the heterogeneity of the cloud cover (Currie et al. 2014). It can help to model the dust environment of young accreting planets, in particular the temperature, radius, and mass of their circumplanetary discs (Wang et al. 2021b; Cugno et al. 2024). Finally, it can possibly hint at auroral signatures in cold planets and brown dwarfs through the detection of CH₄ emission (Faherty et al. 2024). The mid-infrared also has technical advantages: the planet-to-star contrast is lower (e.g. ~4× lower in L′ than in K for β Pic b, Lagrange et al. 2009; Bonnefoy et al. 2011); the atmospheric turbulence has longer coherence lengths and timescales; and the Strehl ratio is higher. The sky emission, however, is higher in the mid-infrared.

One of the first direct emission spectra of an unambiguous¹ exoplanet was obtained in the L band (HR 8799 c with VLT/NACO, Janson et al. 2010). Since then, although mid-infrared photometry has been obtained for many companions, only a few other mid-infrared spectra of directly imaged planets from ground-based telescopes have been published. The list hereafter is exhaustive to the best of our knowledge. High-resolution spectra (R = λ/Δλ > 15 000) were obtained in the L (HR 8799 c, Wang et al. 2018; β Pic b, Janson et al. 2025) and M bands (β Pic b, Parker et al. 2024). High-resolution spectroscopy gives access to many individual spectral lines that can be cross-correlated with molecular template spectra, but not to the spectral energy distribution (SED) of the planet. At low resolution (R = 35), Doelman et al. (2022) presented L-band spectra of HR 8799 c, d, and e with the LBTI/ALES integral field spectrograph, showing their SEDs including broad absorption features. Low- or medium-resolution L-band spectra were also obtained for three isolated or widely separated planetary-mass objects (Miles et al. 2018; Stone et al. 2020). This low number of mid-infrared spectra obtained from the ground may be due to the scarcity of mid-infrared spectrographs assisted by adaptive optics (AO) on 8 m-class telescopes, as high-contrast instruments with extreme AO systems have been developed in the H and K bands so far (VLT/SPHERE, Gemini/GPI, and Subaru/SCExAO, Beuzit et al. 2019; Macintosh et al. 2014; Jovanovic et al. 2015).

More recently, from space, the JWST provided a broad 1–20 μm spectrum of the widely separated (8″) companion VHS 1256-1257 b (Miles et al. 2023; Petrus et al. 2024) and a 5–7 μm spectrum of β Pic b (Worthen et al. 2024). Simulations using molecular mapping or forward modelling techniques applied to the IFU on board the instruments have shown detection limits down to 0.5″ with MIRI (Patapis et al. 2022; Mâlin et al. 2023) and 0.3″ (at constrasts down to 3 × 10⁻⁵) with NIR-Spec (Ruffio et al. 2024). However, characterization below these separations (<3 spaxels and 3λ/D from the host star) remains untested so far, as does the ability to extract the SED of the companion, which is not recovered in molecular mapping.

A push towards spectroscopy at closer angular separations has recently come from ground-based long-baseline interferometry. A dozen planets have been observed with the near-infrared GRAVITY instrument (GRAVITY Collaboration 2017) on the Very Large Telescope Interferometer (see e.g. GRAVITY Collaboration 2019; GRAVITY Collaboration 2020). This has led to precise astrometric measurements enabling well-constrained orbital fits, as well as near-infrared K-band spectra at R = 500, including for companions at separations <0.2″ unreachable by integral field spectrographs (β Pic c, Nowak et al. 2020; HD 206893 c, Hinkley et al. 2023). The 100 m baselines of the VLTI help to efficiently disentangle stellar and planetary photons, providing spectra at close separations and medium spectral resolutions that retain the planet continuum, while resolving some spectral lines.

To extend this capability to the mid-infrared, a new observing strategy, GRA4MAT (Woillez et al. 2024), has been developed to increase the sensitivity of MATISSE, the VLTI’s mid-infrared spectro-interferometer (Lopez et al. 2022). In GRA4MAT, while MATISSE performs science observations, the GRAVITY fringe tracker (Lacour et al. 2019; Nowak et al. 2024) measures the rapidly varying optical path difference (OPD) between telescope light paths introduced by atmospheric turbulence, and corrects the OPD using the VLTI delay lines (instead of GRAVITY’s internal delay lines in a classical GRAVITY observation, GRAVITY Collaboration 2017). As a result, fringes are stabilized for MATISSE in the L, M, and N bands. The narrow off-axis mode adds the possibility of offsetting MATISSE’s optical axis from the fringe tracking source. These implementations lead to two huge improvements in sensitivity: 1) MATISSE can now expose 90× longer than the short exposure times previously required to freeze the turbulence (10 s vs 111 ms), making MATISSE background-limited instead of detector-limited; and 2) by centring MATISSE on the planet, most of the stellar photons fall outside MATISSE’s spatial filter (a 135 mas pinhole, or 1.5λ/D at 3.5 μm), reducing the stellar contamination. A simulation of this effect for the case of β Pic b, with a planet-to-star contrast of 8 × 10⁻⁴ at a separation of 534 mas, is shown in Fig. 1, assuming no turbulence and perfect telescopes. When centred on the planet, the pinhole reduces the speckle flux to the same level as the planetary flux at the planet position.

To demonstrate this new observing window to exoplanets, we observed the iconic β Pictoris b, a giant planet discovered through L′-band imaging at 8 au from its star (Lagrange et al. 2009) inside an edge-on debris disc (Smith & Terrile 1984; Hobbs et al. 1985). β Pictoris is one of the most studied extrasolar planetary systems, whether through direct imaging (Lagrange et al. 2019a; Kammerer et al. 2024), radial velocities hinting at the inner companion β Pic c at 2.7 au (Lagrange et al. 2019b), high-resolution spectroscopy uncovering the fast spin of β Pic b (~20 km/s, Snellen et al. 2014; Landman et al. 2024; Parker et al. 2024), host star astrometry (Snellen & Brown 2018), low- and medium-resolution spectroscopy in the J, H, and M bands (Chilcote et al. 2017; Worthen et al. 2024), molecular mapping with integral field spectroscopy (Hoeijmakers et al. 2018; Kiefer et al. 2024), or K-band spectro-interferometry directly confirming β Pic c (GRAVITY Collaboration 2020; Nowak et al. 2020). The dynamical mass of β Pic b has been constrained from orbital fits of these numerous observations, first using the host star astrometry (Snellen & Brown 2018; Dupuy et al. 2019), then adding radial velocity measurements and planetary astrometry from imaging and interferometry (GRAVITY Collaboration 2020; Vandal et al. 2020; Lagrange et al. 2020; Lacour et al. 2021; Brandt et al. 2021). Dynamical mass estimates from these studies range from 9 to 13 M_Jup. With a flux of 8 mJy (Lagrange et al. 2009), β Pic b is the brightest imaged exoplanet in the L band thanks to its proximity (19.44 ± 0.05 pc, Nielsen et al. 2020) and high temperature of ~1500–1700 K (corresponding to an early L spectral type, Bonnefoy et al. 2013; GRAVITY Collaboration 2020). Although it was discovered in the L′ band, it is still missing an SED in this region. This makes it a prime target for a first demonstration of MATISSE capabilities.

We present in Section 2 the outline of the MATISSE observations of β Pic b. We explain our data reduction in Section 3, including additional steps to the standard MATISSE pipeline. The astrometry and spectrum of the planet is extracted from the reduced data in Section 4. We analyse the spectrum with forward modelling in Section 5, discuss our results in Section 6, and conclude in Section 7.

Fig. 1 Simulation of the gain in contrast achieved by offsetting the pinhole of MATISSE on β Pic b, separated here by 534 mas from its host star. The upper and lower plots depict the stellar and planetary contribution to the coherent flux (the fringe envelope) as seen through the circular MATISSE pinhole (white area, hence the Airy-like shape of the output flux) centred on the star and on the planet, respectively. They are simulated at a wavelength of 3.5 μm and normalized to the maximal intensity of the star-centred stellar PSF. When the pinhole is centred on the planet, the speckle flux is reduced to the same level as the planetary flux. In each simulation, the asymmetric shape of the PSF of the offset object (relative to the pinhole) is due to imperfect spatial filtering of the Airy rings falling in the pinhole.

2 Observations

β Pic b was observed with MATISSE during the commissioning of the GRA4MAT narrow off-axis mode (Woillez et al. 2024), on the nights of November 8, 2022 and February 3, 2023². We obtained observations at R = 500 between 2.75 and 5 μm, covering the full L band (2.8–4.2 μm) and the blue side of the M band (4.5–5 μm). The planet position was predicted with whereistheplanet³ (Wang et al. 2021a), which uses orbital fits of archival planet astrometry (Lacour et al. 2021) with the orbitize! package (Blunt et al. 2020). We followed two different observing sequences in November 2022 and February 2023, summarized in Table 1. In November 2022, we alternated for 2 h between two or three cycles targeting the planet and one targeting the star. A cycle consists of four successive exposures using each of MATISSE’s beam commuting device (BCD) configurations. The total exposure time amounts to 44 min on the planet and 16 min on the star. As is explained in Sect. 3.1, the stellar pointings are used to calibrate the telluric contamination and anchor the stellar contribution to the stellar pointing in the planet observations. In February 2023, we added a so-called ‘anti-planet’ pointing at the opposite position of β Pic b with respect to the star, in order to confirm the astrophysical nature of the signal recorded on the planet. This shorter observation amounted to one cycle on the star (4 min), two on the planet (8 min), and one at the anti-planet position (4 min). We show in Appendix A a comparison of coherent flux ratios at the planet and anti-planet positions. As was expected, no planet-like signal is found at the anti-planet position. Given its short exposure time, higher seeing, higher water vapour, and for data homogeneity, we discard the 2023 dataset and focus on the 2022 dataset in the following sections.

3 Data reduction

In interferometry, the quantity of interest is the coherent flux (also known as complex correlated flux), a complex quantity describing the amplitude and phase of the interferometric fringe pattern. It encodes the Fourier transform of the source spatial brightness distribution, and is sampled by the interferometer at given spatial frequencies defined by the baseline length and orientation relative to the source. A general introduction of the technique can be found in, for example, Buscher (2015). The MATISSE pipeline reduces the fringe pattern images and extracts the coherent flux (Sect. 3.1). We developed custom reduction steps to correct for remaining instrumental and atmospheric effects (Sect. 3.2).

3.1 MATISSE and GRA4MAT standard pipeline

We used the standard MATISSE pipeline (version 2.0.2⁴) in ‘correlated flux’ mode (option ––corrFlux=True) to reduce the observations. The fringes recorded by MATISSE are Fourier-transformed, providing the observer with a complex coherent flux for each frame. The frames are then carefully selected: MATISSE frames during which the GRA4MAT fringe tracker experienced a jump from one fringe to another (i.e. a phase jump of a 2π multiple) are discarded. This is necessary as the fringe tracking and science channels of GRA4MAT are in different bands: a 2π jump in K results in a very different jump in L that is generally not a multiple of 2π. Without knowledge of when jumps happen during a frame, this creates a highly variable fringe contrast that cannot be calibrated. It is therefore preferable to simply flag and remove these frames. In the November 2022 dataset, 60 out of the 360 frames had jumps (~17%). This proportion greatly improved after a later update of the fringe tracker (Nowak et al. 2024; Woillez et al. 2024): no frame was rejected in the February 2023 dataset.

After this selection, the MATISSE pipeline fits and removes in the coherent flux any residual phase created by an achromatic OPD, and coherently averages the frames over the duration of an exposure (1 min). After inspection of the final products, we note that the pipeline does not fully remove residual achromatic OPDs in L and M bands. It also does not correct for chromatic dispersion near the telluric absorption bands. For these reasons, we developed our own phase correction method described in the following subsection. This has the additional advantage of providing individual frame products instead of average exposure products. In this custom phase correction, we use the intermediate pipeline products nrjReal and nrjImag. They contain the real and imaginary parts of the coherent flux extracted from the raw interferogram, after correction of the bias, bad pixels and distortion, but before the achromatic OPD removal.

3.2 Custom steps

3.2.1 Phase correction

In interferometry, an OPD δ(λ) generates a shift of 2πδ(λ)/λ in the fringe phase. The GRA4MAT fringe tracker measures and compensates the K-band OPD originating in the line of sight of GRAVITY, i.e. in the common path of GRAVITY and MATISSE. Additional atmospheric and instrumental phase terms remain, however, as is shown in Fig. 2, an example of phase measurements in two on-star and on-planet frames before the pipeline OPD correction. As is explained in Sect. 4, the phase created by the targets should be zero when pointing at the star (characteristic of an unresolved point source), or be modulated by an oscillatory function of 1/λ when pointing at the planet (characteristic of a resolved binary source, here the stellar halo and the planet). As can be seen here, the actual measured phase contains two additional components: a linear term as function of 1/λ, and high-order terms near the telluric absorption bands defining the edges of the L and M bands. Linear terms should in principle be measured and removed by the fringe tracker. As it appears uncorrected by GRA4MAT, it must originate in a ‘non-common path’ achromatic OPD introduced between GRAVITY and MATISSE, which would be left unseen by the fringe tracker. It is so far only partially fitted and subtracted by the last recipe of the MATISSE pipeline. The high-order terms are created by chromatic variations in the air refractive index (i.e. chromatic dispersion) near the main telluric absorption bands (CO₂ and H₂O bands at ~2.7 μm, CO₂ band at ~4.2 μm). We aim to correct both the non-common path (NCP) achromatic OPD and the common-path (CP) chromatic variations. A detailed demonstration of the method is presented in Appendix B. We present here a summary of the main points.

We first defined each of the OPDs introduced in the successive mediums along the common path: vacuum, atmosphere, and delay lines. We simulated these terms and their correction by the fringe tracker using an air refractive index model (Voronin & Zheltikov 2017) and pessimistic assumptions on the differences in ambient conditions between lines of sight. We found that one of these corrected OPDs, originating in the additional length travelled in space and its compensation by the delay lines in air, is more than ten times larger than the other terms in L and M bands, except at the bluer edge of the bands. This leads us to approximate the corrected common-path OPD as follows: $${\delta }_{\mathrm{C}\mathrm{P}}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}(\lambda )\approx (\frac{n_{\mathrm{D}\mathrm{L}}(\lambda )}{{⟨n_{\mathrm{D}\mathrm{L}}(\lambda )⟩}_K}-1)\mathrm{\Delta }d,$$(1)

in which n_DL(λ) is the air refractive index in the delay lines, ⟨n_DL(λ)⟩_K is its average in K band, and Δd is the additional length travelled in vacuum, compensated in the delay lines. Overall, our phase model after fringe tracking (Φ^corr) is the sum of the target phase (Φ_obj) and the phases originating in the corrected common-path OPD and the non-common path OPD (δ_NCP): $$\begin{array}{rl}{\mathrm{\Phi }}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}(\lambda )& \approx {\mathrm{\Phi }}_{\mathrm{o}\mathrm{b}\mathrm{j}}(\lambda )+\frac{2\pi }{\lambda }[{\delta }_{\mathrm{C}\mathrm{P}}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}(\lambda )+{\delta }_{\mathrm{N}\mathrm{C}\mathrm{P}}]\\ & \approx {\mathrm{\Phi }}_{\mathrm{o}\mathrm{b}\mathrm{j}}(\lambda )+\frac{2\pi }{\lambda }[(\frac{n_{\mathrm{D}\mathrm{L}}(\lambda )}{{⟨n_{\mathrm{D}\mathrm{L}}(\lambda )⟩}_K}-1)\mathrm{\Delta }d+{\delta }_{\mathrm{N}\mathrm{C}\mathrm{P}}].\end{array}$$(2)

Before fitting this model to the data, we removed from both of them their mean L-band phase in order to get the differential phase, a standard step in interferometry when the absolute phase reference is lost during fringe tracking: $${\mathrm{\Phi }}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}(\lambda )={\mathrm{\Phi }}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}(\lambda )-{⟨{\mathrm{\Phi }}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}(\lambda )⟩}_L.$$(3)

During our reduction, only the non-common path term δ_NCP was fitted. We found that an achromatic constant term per frame and baseline is sufficient to reproduce the data well. This likely stems from the differential dispersion between K and L bands, mostly due to water vapour variations; and from a possible slow instrumental drift between GRAVITY and MATISSE during the night, as the instruments are only aligned once before the observations. The common-path term is not fitted but modelled thanks to the path lengths and ambient conditions (temperature, pressure, and humidity) measured in the delay line tunnel and reported in the FITS header of the MATISSE frames. The ambient conditions are fed into air refractive index models, for which we tested three different ones. The model of Voronin & Zheltikov (2017) presents interesting features in comparison to previous widely used models: it takes into account the contribution of the wings of the main infrared absorption lines, which are not considered in Ciddor (1996), and enables the variation in the CO₂ concentration (as well as other gases) contrary to the models of Mathar (2007) tabulated at a fixed value of 370 ppm, which is every year drifting away from the actual CO₂ concentration (~417 ppm in 2022⁵) rising from anthropogenic emissions. This is especially important in the L and M bands that are separated by a strong CO₂ absorption feature.

We show in Fig. 2 the measured differential phase in two frames, one being centred on the star and the other on the planet. We overlay the fits of our model using the three different air refractive index formulae. The one of Voronin & Zheltikov (2017) is clearly favoured, providing a better fit across the L band including very close to the CO₂ line at 4.2 μm. We note that in the planet data, the higher-frequency planet–star modulation does not affect the fit of δ_NCP, if bounds are provided on this term during the fit. We also note some discrepancies between model and data near the red L-band edge and in M band. They are likely caused by the OPD terms we neglected in Eq. (1), which originate in ambient condition differences between lines of sight in the atmosphere and the delay lines. These discrepancies appear almost constant between successive planet and star pointings. They are thus well removed by subtracting from each planet phase the average of the corrected stellar phases from the next observation block (OB). This step is inherently part of the calibration strategy presented in the next section. For this reason, a fit of the ambient condition differences in addition to δ_NCP does not seem to be required at the moment. This could nonetheless be the subject of further study if more precision is required. Our final estimate of the object differential phase is the difference between the measured differential phase and the fit, using the air refractive index of Voronin & Zheltikov (2017): $${\mathrm{\Phi }}_{\mathrm{o}\mathrm{b}\mathrm{j}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}(\lambda )\approx \arg \left[{\mathrm{e}}^{\mathrm{i}{\mathrm{\Phi }}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}(\lambda )}{\mathrm{e}}^{-\mathrm{i}\frac{2\pi }{\lambda }({\delta }_{\mathrm{C}\mathrm{P}}^{\text{corr\hspace{0.17em}}}(\lambda )+{\delta }_{\mathrm{N}\mathrm{C}\mathrm{P}})}\right].$$(4)

Fig. 2 Measured differential phase (grey dots) in two MATISSE frames, one on the star (top) and the other on the planet (bottom). Overlaid are the models presented in this work using the air refractive indices of Voronin & Zheltikov (2017) (blue), Mathar (2007) (green), and Ciddor (1996) (red). In the absence of chromatic or non-common path OPD, the differential phase should be zero on the star or an oscillation on the planet.

3.2.2 Calibration and binning

Frames on the star and the planet are affected by telluric lines. In addition, the on-planet coherent fluxes are contaminated by the stellar halo leaking into the pinhole. To correct it, we used the on-star frames as a calibrator proxy for the on-planet frames. We divided the complex coherent flux of each on-planet frame by the mean coherent flux of the on-star frames in the next observation block (see Table 1). This removed most of the telluric lines and the high-order variations in the stellar contaminating spectrum, only leaving a low-order function created by the chromatic variation in the stellar speckle halo. This can be seen in the data in Fig. 3, in which the amplitudes of the planet-star modulations vary across wavelengths. This contamination will be modelled in the next section.

The MATISSE LM-band detector oversamples the spectral resolution by a factor of five in the medium-resolution mode (R = 500)⁶. In order to gain signal-to-noise (S/N) and prevent correlations between spectral pixels, while conserving the spectral resolution, we binned the coherent flux and their errors into blocks of 5 pixels.

4 Astrometry and spectrum extraction

The extraction of the planet’s astrometry and spectrum follows the outline of the method used by the ExoGRAVITY collaboration (see the appendices of GRAVITY Collaboration 2020; Nowak et al. 2020) with some modifications specific to MATISSE and its lack of absolute metrology. On each baseline, assuming the star and planet are point sources, we modelled the complex coherent fluxes obtained when MATISSE is centred on the planet (F_p) and on the star (F_⋆) as follows: $$F_p(\lambda ,t_p,{\mathit{r}}_p)=T(\lambda ,t_p)[S_p(\lambda )+R(\lambda ,t_p,{\mathit{r}}_p)S_{\star }(\lambda ){\mathrm{e}}^{\mathrm{i}\frac{2\pi \alpha u}{\lambda }}]{\mathrm{e}}^{\mathrm{i}\mathrm{\Phi }(\lambda ,t_p)},$$(5) $$F_{\star }(\lambda ,t_{\star },{\mathit{r}}_{\star })=T(\lambda ,t_{\star })S_{\star }(\lambda ){\mathrm{e}}^{\mathrm{i}\mathrm{\Phi }(\lambda ,t_{\star })},$$(6)

in which the p and ⋆ indices designate the on-planet and on-star frames, respectively. λ, t, and r are the wavelength, the observation time, and the pinhole location on sky; and T, S_p, and S_⋆ are the telluric and instrumental transmission, the planet spectrum, and the stellar spectrum, respectively. R is the star-to-speckle contrast, i.e. the wavelength-dependent stellar point spread function (PSF) at the planet location, which depends on the seeing and AO correction. Together, RS_⋆ thus models the speckle spectrum at the planet location. α is the (Δα, Δδ) angular offset vector of the planet from the star, and u is the (u, v) baseline vector of the observation. Finally, Φ is a phase component not originating in the targets, i.e. uncorrected atmospheric or instrumental phase components. Given the good correction of the chromatic dispersion and non-common path OPD obtained in Sect. 3.2.1, we assume in the following steps that this residual phase is negligible: Φ(λ) ≈ 0.

We use in our study the ratio between the on-planet and on-star coherent fluxes, ℱ: $$\mathcal{F}=\frac{F_p(\lambda ,t_p,{\mathit{r}}_p)}{F_{\star }(\lambda ,t_{\star },{\mathit{r}}_{\star })}\approx \tau (\lambda ,t_p,t_{\star })[C(\lambda )+R(\lambda ,t_p,{\mathit{r}}_p){\mathrm{e}}^{\mathrm{i}\frac{2\pi \alpha u}{\lambda }}],$$(7)

in which C(λ) = S_p(λ)/S_⋆(λ) is the planet-to-star contrast spectrum, and τ(λ, t_p, t_⋆) = T(λ, t_p)/T(λ, t_⋆). To calculate this ratio, the coherent flux, F_p, of each planet frame was divided by the average of the stellar coherent fluxes, F_⋆, in the next stellar OB. We did not calibrate the telluric transmission with dedicated software such as Molecfit (Smette et al. 2015; Kausch et al. 2015). This would be possible on the on-star data but not on the on-planet data, which is too noisy for such an analysis. We instead relied on the small variations in airmass and water vapour between observation blocks (see Table 1), and assumed that the telluric transmission ratio between on-planet and on-star exposures can be modelled by a simple achromatic factor: τ(λ, t_p, t_⋆) = τ(t_p, t_⋆). Following GRAVITY’s two-step extraction procedure, we first extracted the astrometry and the starlight transmission function, and then the contrast spectrum.

Fig. 3 Real (left) and imaginary (right) parts of the measured planet-to-star coherent flux ratio, ℱ, on the six UT baselines, during one 10 s frame. The error range is plotted in the background. The thick lines show the best-fit model as described in Sect. 4.1. It includes the star–planet modulation, the stellar speckle contamination (dashed blue curves), and the planet-to-star contrast assumption (dashed black curve).

4.1 Astrometry and starlight contamination fitting

To first extract the astrometry, we need to assume a model C_model(λ) for the planet-to-star contrast spectrum. We used for this purpose the ratio between a model of stellar spectrum and a model of planetary spectrum. For the planet, we took a BT-Settl (Allard et al. 2013) sub-stellar template close to the known parameters (T_eff = 1700 K, log g = 4.0, [M/H] = 0.0). For the star, we used a BT-NextGen (Allard et al. 2012) template interpolated with species⁷ (Stolker et al. 2020) at T_eff = 7890 K and log g = 3.83 (Swastik et al. 2021), and scaled according to archival photometric fluxes from Tycho (Høg et al. 2000), 2MASS (Skrutskie et al. 2006), and Gaia (Gaia Collaboration 2023), as shown in Fig. 4. The fitted scale factor has an uncertainty of ~1.5% stemming from the photometric flux uncertainties. Each template was convolved at R = 500 and resampled on the MATISSE wavelength grid using spectres (Carnall 2017). The ratio was finally scaled according to the reported L′-band contrast (7.7 mag, Lagrange et al. 2009). Based on Eq. (7), our fit is now described by $${\mathcal{F}}_{\text{model\hspace{0.17em}}}(\lambda )=\mathcal{T}(t_p,t_{\star })C_{\text{model\hspace{0.17em}}}(\lambda )+\mathcal{R}(\lambda ,t_p,{\mathit{r}}_p){\mathrm{e}}^{\mathrm{i}\frac{2\pi \alpha u}{\lambda }},$$(8)

where 𝒯(t_p, t_⋆) = γ τ(t_p, t_⋆) is the product of the telluric and instrumental transmission ratio, τ, with a factor, γ, reflecting possible inaccuracies in the scaling of the contrast template, and ℛ(λ, t_p, r_p) = τ(t_p, t_⋆) R(λ, t_p, r_p) is the product of τ with the stellar contamination R. Following the ExoGRAVITY method, we assume that the stellar contamination varies slowly with wavelength and therefore model ℛ with a low-order polynomial. The fit was performed baseline by baseline (b), frame by frame (t_p), and for each point of a grid of tested planet offsets (Δα, Δδ). We minimized the least-square residuals of the real and imaginary parts of the coherent flux ratios (considered independent), using errors propagated from those provided by the MATISSE pipeline: $$\begin{array}{rl}{\chi }^2(t_p,b,\mathrm{\Delta }\alpha ,\mathrm{\Delta }\delta )=& \sum _{\lambda }\frac{\mathrm{\Re }\mathrm{e}{[\mathcal{F}(\lambda )-{\mathcal{F}}_{\text{model\hspace{0.17em}}}(\lambda )]}^2}{{\sigma }_{\mathrm{\Re }\mathrm{e}[\mathcal{F}(\lambda )]}^2}\\ & +\sum _{\lambda }\frac{\mathfrak{I}\mathrm{m}{[\mathcal{F}(\lambda )-{\mathcal{F}}_{\text{model\hspace{0.17em}}}(\lambda )]}^2}{{\sigma }_{\mathrm{\Im }\mathrm{m}[\mathcal{F}(\lambda )]}^2}.\end{array}$$(9)

For each of these iterations, we fitted the 𝒯 factor and the polynomial coefficients of ℛ. We tried to fit either the same or two different polynomial functions to the real and imaginary parts of the flux, and tried several polynomial orders. Overall, we found that fitting two different first-order polynomial functions for the real and imaginary parts is slightly favoured in terms of reduced chi-squared (${\chi }_r^2$). In total, we thus have five fitted parameters per fit: one for the transmission factor, 𝒯, and four for the two linear functions, ℛ_real and ℛ_imag. We excluded from the fit the regions of strong telluric absorption below 2.87 μm and between 4.15 and 4.57 μm. We also excluded outliers deviating from the median coherent flux ratio by more than 3σ (computed from the median absolute deviation: σ = 1.4826 MAD) in L and M bands separately, amounting to 3–7% of the data. An example of match between the model and the data of a single frame is shown in Fig. 3.

With a ${\chi }_r^2$ value for each tested planet offset, we obtained one χ² map per frame and baseline. We averaged all these maps together to get the total ${\chi }_r^2$ maps presented in Fig. 5. We ran this procedure on two different grids: a large one the size of MATISSE’s pinhole (130 mas) with a resolution of 1 mas, and a narrow one on a 20 × 20 mas square with a resolution of 0.1 mas. The grid is centred on the telescope pointing during the planet observation. We fitted an ellipse on the lowest peak of the narrow map to extract the best-fitting relative astrometry of β Pic $\mathrm{b}(\hat{\mathit{\alpha }})$, its associated errors, and the correlation coefficient, ρ, between Δα and Δδ measurements⁸. The values are reported in Table 2, along with the predictions of the orbital fit of Lacour et al. (2021) based on GRAVITY, GPI, and SPHERE data, computed with whereistheplanet at the same epoch as the MATISSE observations.

The orbital fitting prediction is located 3.38 mas (~2σ) away from our MATISSE astrometry. The GRAVITY measurements it is based on are obtained by averaging the best-fit astrometries per frame or group of frames, which provides precisions of ~0.1 mas, an order of magnitude below the size of individual peaks in the χ² maps. This method requires a consistent peak pattern between maps, which we do not get at the moment with MATISSE. We instead rely on the average of all χ² maps, and are thus limited to a precision of a few milliarcseconds equivalent to the VLTI angular resolution. The stability of the GRAVITY astrometry stems from its internal metrology system (GRAVITY Collaboration 2017), which enables the phases measured on the planet to be anchored to the phases measured on the star, and correcting non-common path aberrations between fringe tracker and science combiner. MATISSE does not have such metrology, and is thus affected by time variations in residual non-common path aberrations between the star and planet pointings. In Figure B.3 we show the OPD drift measured on MATISSE, exhibiting the differential drift between GRA4MAT (fringe tracker) and MATISSE (science instrument). It currently is of the order of one wavelength over one hour. This drift may be reduced to sub-wavelength accuracy by adopting a different observing strategy in the future. To alleviate this issue for the moment, the data reduction is based on self-referenced phases (see Sect. 3.2.1), resulting in an additional astrometric noise of the order of about λ/B, with B the baseline length.

Fig. 4 Scaling of the BT-Nextgen stellar model (Teff = 7890 K, log g = 3.83, [Fe/H] = 0.0) with species using Tycho, Gaia, and 2MASS photometry of β Pictoris. The bottom panel shows the residuals of the fit.

Fig. 5 Detection maps of β Pic b with MATISSE. The maps show the reduced χ2 values after fitting the stellar contamination for each tested planet astrometry (Δα, Δδ) in the grid. The grids are centred on MATISSE’s pointing during the planet observations. The left map is the size of the pinhole (1.5λ/D at 3.5 μm) and has a resolution of 1 mas. The right map was generated at higher resolution (0.1 mas) on a smaller 20 × 20 mas grid. The most likely astrometry from the MATISSE observations is shown as a blue ellipse, found by fitting a 2D Gaussian curve on the lowest ${\chi }_r^2$ peak. The planet location and uncertainties predicted by whereistheplanet from the orbital fitting of GRAVITY observations (Lacour et al. 2021) are shown as a green ellipse. The multiplicity of ${\chi }_r^2$ peaks is due to the interferometric nature of our data. As coherent flux is periodic as a function of α.u (see Eq. (7)), an infinity of astrometric solutions can reproduce the signal. This degeneracy gradually disappears when combining data from different baseline lengths and orientations. The ${\chi }_r^2$ of neighbour peaks can nonetheless still remain close. Using covariances between baselines and wavelengths could help increase the difference between peaks, and will be studied in a future work.

Fig. 6 Spectrum of β Pic b with MATISSE (black, this work) compared to other instruments: NACO (yellow, new data), GPI (blue, Chilcote et al. 2017), GRAVITY (green, new data) and MIRI (red, Worthen et al. 2024). Only the slope of the NACO spectrum can be compared to MATISSE as the absolute level was scaled manually. MATISSE data with a S/N < 5, located in strong telluric bands and not used in this work, are shown in grey. The absorption curves of H2O, CO, and CH4, calculated from the best-fitting Exo-REM model presented in Sect. 5, are plotted at the top to visualize their impact on the spectrum. H2O and CO leaves a visible pattern, but not CH4, likely because it is dominated by H2O absorption or buried in noise.

4.2 Spectrum extraction

Once we had the best-fitting astrometry, $\hat{\mathit{\alpha }}$, we extracted the associated fit of the stellar contamination function, $\hat{\mathcal{R}}$, at this position, for each frame and each baseline. Following Eq. (7), we then estimated the contrast spectrum as follows: $$\hat{C}(\lambda ,t_p,b)=\mathrm{\Re }\mathfrak{e}[\frac{\mathcal{F}(\lambda ,t_p,b)}{\tau (\lambda ,t_p,t_{\star })}-\hat{\mathcal{R}}(\lambda ,t_p,b,\hat{\mathit{\alpha }}){\mathrm{e}}^{\mathrm{i}\frac{2\pi \hat{\alpha }u}{\lambda }}].$$(10)

We neglected the transmission ratio, τ(t_p, t_⋆), as 𝒯(t_p, t_⋆) was consistently fitted close to one in the previous section. With this method, we obtained 1308 estimates of the contrast spectrum (218 frames × 6 baselines). Before averaging them, we first excluded points outlying by more than 3σ (estimated from the median absolute deviation) from the median at each wavelength, which removed ~1% of points inside the L and M bands. We finally computed the mean and the covariance of the mean to get our final estimate of the contrast spectrum and its errors. The covariance matrix was computed by considering the 1308 contrast spectra as samples and the wavelengths as variables.

In order to get the planet spectrum, we have to multiply the contrast by a stellar spectrum. We initially used a β Pic spectrum from the Infrared Space Observatory (Pantin et al. 1999) at R = 2000, but the high-resolution archival data is affected by calibration issues in the M band, and it does not fully cover the K-band wavelengths of GRAVITY which we will use in the modelling. For these reasons, we instead use the same BT-NextGen stellar template as computed in Sect. 4.1. The resulting planet spectrum is shown in Fig. 6, covering the full L band and the blue side of the M band. We added several other spectra of β Pic b from different instruments: Gemini/GPI (Chilcote et al. 2017), JWST/MIRI (Worthen et al. 2024), and two original spectra from VLTI/GRAVITY (S. Lacour, private communication) and VLT/NACO (M. Bonnefoy, private communication). The GRAVITY spectrum is a weighted mean of several epochs since the first publication in GRAVITY Collaboration (2020). We used the contrast spectrum and scaled it with the same stellar model as the MATISSE data. The NACO spectrum was obtained in 2011. The acquisition and reduction of these new spectra are described in Appendix C. We note that the NACO spectrum was scaled manually, so only its slope can be compared to the MATISSE spectrum.

The LM-band MATISSE spectrum of β Pic b seems to be affected by broad absorption features. To verify how the main absorbers, H₂O, CO, and possibly CH₄, impact the SED, we calculated their individual absorption curves as in Palma-Bifani et al. (2024). These curves were obtained from the best-fitting Exo-REM model. This model and the method to generate molecular absorption curves are presented in Sect. 5. The absorption curves are shown at the top of Fig. 6. They show that the L band (2.8–4.2 μm) is shaped by H₂O, while the M-band subregion observed by MATISSE (4.5–5 μm) is shaped primarily by CO, and secondarily by H₂O. The CH₄ band at 3.3 μm does not seem visible, either because H₂O absorption dominates it or because it is buried in noise.

Some small-amplitude periodic features also appear in the spectrum, in particular between 3 and 3.5 μm. They seem associated with off-diagonal periodic correlations in the covariance matrix of the spectrum, presented in Fig. 7. This matrix was built by scaling the contrast covariance with the stellar model: σ_Sp (λ_i, λ_j) = S_⋆(λ_i) S_⋆(λ_j) σ_C(λ_i, λ_j). We believe these periodic correlations come from the instrument and/or the processing rather than the source itself. Their frequency seems indeed related to the star–planet modulations of the coherent flux, i.e. the phasor term in Eq. (7), whose frequency also decreases with wavelength (∝ 1/λ). We note that even if we cannot distinguish these residuals from spectral features, their presence in the covariance matrix should prevent them from biasing the atmospheric parameters fitted through forward modelling (Sect. 5).

These periodic residuals could come from an imperfect fit of the stellar speckle function (ℛ in Eq. (8)). The planet-to-speckle contrast in our observations is close to one, while it was 3–5% in the GRAVITY observations (GRAVITY Collaboration 2020). This is expected from several factors: the 4× brighter intrinsic planet-to-star contrast in L than in K, the larger separation at the time of our observation (534 vs 144 mas), and the larger spatial filter used in MATISSE (135 vs 60 mas). A planet-to-speckle contrast of order unity means that C (the planet-to-star contrast) and ℛ (the star-to-speckle contrast) are on the same order of magnitude, and might thus be more correlated in our fit than they did in GRAVITY observations. To ensure that the assumed contrast model in Sect. 4.1 does not influence the shape of the final contrast spectrum, we ran the same fit using a flat contrast model at 7 × 10⁻⁴. The contrast spectrum extracted with this simplistic assumption has the same slope and absolute level as the one obtained with the realistic assumption, with only the oscillations below 3.5 μm varying slightly. As is seen in Fig. 3, the large bandwidth of MATISSE provides many periods of the star–planet modulations, likely enabling a good fit of the stellar speckle contamination ℛ without large correlation with the chosen contrast model. Nonetheless, a simultaneous fit of C and ℛ is preferable to improve accuracy. Since C(λ) is defined for hundreds of wavelengths, this fit can only be performed by fitting simultaneously on a large number of DITs, as is done for GRAVITY observations. This supposes that the coherent fluxes of all DITs are cophased. We cannot assume this at the moment for MATISSE due to its lack of metrology, which is why we extracted our spectrum DIT by DIT instead of fitting it simultaneously on all the data. Methods to co-phase MATISSE DITs together are being investigated and will be the subject of a future publication.

The relative errors shown on the MATISSE spectrum were extracted from the covariance diagonal. They range from 1% of the flux in the middle of the L band to ~5% in the M band. The absolute MATISSE flux matches well with the MIRI spectrum in their overlapping region, but seems to have a mismatch of ~10% with the GRAVITY spectrum, as well as with previous photometric measurements. We computed an estimated L′ magnitude by convolving this spectrum with the Paranal/NACO L′ filter available on the SVO Filter Profile Service (Rodrigo et al. 2012; Rodrigo & Solano 2020). We find $L_{\text{MATISSE}}^{\mathrm{\prime }}=10.970\pm 0.005$, which is 23.5% higher but within the error bars of the measurement of Lagrange et al. (2009) (11.2 ± 0.3 mag). Flux discrepancies have been noted in many studies combining different instrument spectra, and can be seen here as well between GPI and GRAVITY. Several reasons may explain them in our case.

Firstly, the absolute precision is limited by the 1.5% precision on the scaling of the stellar model (see Sect. 4.1), but this cannot explain the mismatch with GRAVITY as the same model was used. Secondly, the star is partly resolved by MATISSE. Priolet et al. (2025) find a visibility decrease of ~5% on the longest baselines on β Pic, interpreted as the resolved star surrounded by an inner disc. This means that the stellar coherent flux is affected by this visibility and is lower than the stellar spectrum. As a result, the planet spectrum is overestimated if we do not take into account the stellar visibility when we multiply the contrast by the stellar template. The inner disc may also add emission in the stellar spectrum, which is not taken into account in the stellar template. The β Pic interferometric data of Priolet et al. (2025) could be used to estimate and correct these effects. Another factor may be the loss of coherent flux induced by phase jittering during a science integration (Colavita 1999; Tatulli et al. 2007). It is taken into account in the GRAVITY pipeline through the calculation of the so-called ‘vFactor’ based on the high-frequency fringe tracker data⁹. We implemented it in our MATISSE processing by using the simultaneous GRAVITY fringe tracking data, but it results in no substantial difference in the final planet spectrum, possibly because of close vFactors in the on-planet and on-star frames cancelling each other when we take the ratio. Another final possibility is a phase jitter not seen by the fringe tracker when recording on the planet in the narrow VLTI field.

Fig. 7 Covariance (upper right triangle) and correlation (lower left triangle) matrices of the MATISSE spectrum of β Pic b. In order to highlight the off-diagonal correlations in the L band, the covariance values are clipped between 3 × 10−35 and 3 × 10−33 W2 m−4 μm−2 and the correlation values between 0.01 and 1. Two lines of off-diagonal correlations are further highlighted with arrows for readability.

5 Forward modelling

To check the physical consistency of the spectrum obtained with MATISSE, we compared it along the new GRAVITY spectrum to the same grid of Exo-REM model spectra (Baudino et al. 2015; Charnay et al. 2018; Blain et al. 2021) that was used in GRAVITY Collaboration (2020). Exo-REM is a self-consistent atmospheric model assuming radiative-convective equilibrium, incorporating simple cloud microphysics that was shown to well reproduce the L-T transition of brown dwarfs and giant planets. For more advanced modelling, we refer to the study of Ravet et al. (2025) which uses all the available archival data of β Pic b (including our MATISSE spectrum) and four different grids of models (Exo-REM, ATMO, SONORA, BT-Settl).

To fit the models, we used the Bayesian inference ForMoSA¹⁰ code (Petrus et al. 2021, 2023; Palma-Bifani et al. 2023), which evaluates posterior density functions on each free parameter of the model grid based on a nested sampling algorithm. The code performs an N-dimensional linear interpolation of the precomputed grids of models between the grid nodes on the fly. The Exo-REM grid has four dimensions, corresponding to the set of free parameters of the self-consistent Exo-REM models: effective temperature (T_eff), surface gravity (log g), carbon-to-oxygen ratio (C/O), and metallicity ([M/H]). In addition to the grid parameters, ForMoSA fits the companion bolometric luminosity (or radius or distance if one of these is provided as fixed prior), radial velocity, and rotational velocity. Finally, it can fit scale factors (α) to different datasets in order to account for flux calibration issues.

β Pic b is part of a limited sample of directly imaged planets with a dynamical mass constrained from orbital fits of astrometric measurements (see Sect. 1). To take advantage of this knowledge, we implemented in ForMoSA the possibility to fit and use the mass as a prior. This replaces in practice the planet radius, which is now set by the sampled mass and surface gravity through the surface gravity definition: $$R=\sqrt{\frac{GM}{{10}^{\log g}}}.$$(11)

This in turn sets the bolometric luminosity of the model via the dilution factor (R/d)² applied to the synthetic spectra to convert them to apparent fluxes. The priors of our fits are listed in Table 3. The mass and distance are set with Gaussian priors based on stellar and planetary astrometry from Nielsen et al. (2020) and GRAVITY Collaboration (2020), respectively. All the other parameters are set with uniform priors. The radial and rotational velocities are not fitted, as their effects are negligible at R = 500.

We fitted models successively on MATISSE only, GRAVITY only, and on GRAVITY and MATISSE spectra together, using 500 live points, and taking into account their covariance matrices. For the joint GRAVITY+MATISSE fits, we perform two fits in which we use a scale factor either on GRAVITY or on MATISSE, to account for the discrepancy between spectra noted in Sect. 4. We show in Fig. 8 the best-fit Exo-REM models together with the MATISSE and GRAVITY data. Our estimates on each parameter and their fitting errors (inferred from the 68% confidence intervals on the posteriors) are reported in Table 3. We caution that the uncertainties are only a projection of the observational errors on the model predictions and do not account for systematic deviations of the models that certainly dominate here. A more systematic study of these deviations is explored in Ravet et al. (2025) based on the full collection of spectrophotometry obtained on β Pic b thus far.

The best models for the joint GRAVITY+MATISSE fits are found for T_eff = 1529 ± 3 K, log g = 3.83 ± 0.04, [M/H] = 0.30 ± 0.03, and C/O = 0.539 ± 0.003. Scaling either the GRAVITY or MATISSE spectra only changes the fitted mass, which sets the radius and thus the absolute flux level in conjunction with log g. We find masses of 10.5 ± 0.7 and 11.8 ± 0.9 M_Jup when scaling MATISSE or GRAVITY, respectively. Figure 9 shows the corner plot of the best-fit model (GRAVITY+MATISSE fixed). We finally note that fitting only the GRAVITY data provides posteriors that are very close to the joint fits.

From the best-fitting model, we generated the absorption curves of individual molecules (H₂O, CO, and CH₄) plotted in Fig. 6. These curves are built in the same way as full Exo-REM spectral models, but only including one molecule at a time. We consider collision-induced absorption (CIA) of H₂–H₂ and H₂-He, Rayleigh scattering, and the cross-section of each molecule. We used the P-T profile and volume-mixing ratios of the best-fitting model, and generated absorption curves for each molecule with the opacity library Exo_k (Leconte 2021). We finally normalized each molecular curve by the CIA-only spectrum for visualization.

Fig. 8 GRAVITY and MATISSE spectra of β Pic b with the best-fit Exo-REM models. The top plot shows the best fits on GRAVITY only, MATISSE only, and GRAVITY+MATISSE, when no scaling is applied on the spectra. The two other plots show the best fits on GRAVITY+MATISSE when a scale factor is applied either on MATISSE (middle) or on GRAVITY (bottom).

6 Discussion

6.1 Solar C/O

We find C/O ratios compatible with or above solar abundances, varying from 0.524 ± 0.004 (GRAVITY data alone) up to 0.775 ± 0.002 (MATISSE data alone). We adopt here the value of 0.539 ± 0.003 found when jointly fitting both datasets. We note that the errors we get from these fits are an order of magnitude smaller than the errors on the C/O of β Pic b obtained by GRAVITY Collaboration (2020). These errors are likely underestimated for the reasons mentioned in Sect. 5.

This new solar C/O is higher than the subsolar values of 0.43 ± 0.05 and 0.41 ± 0.04 found by GRAVITY Collaboration (2020) and Landman et al. (2024), respectively. The former was obtained by forward modelling (with Exo-REM) and atmospheric retrieval (with petitRADTRANS, Mollière et al. 2020) of JHK-band GRAVITY and GPI spectra. The latter was obtained by petitRADTRANS atmospheric retrieval of a K-band high-resolution (R = 100 000) CRIRES spectrum. The solar C/O is, however, in line with the value of 0.551 ± 0.002 found by Kiefer et al. (2024) with molecular mapping of SINFONI data. A high-resolution K-band survey with KPIC also found solar C/O for eight young sub-stellar companions (Xuan et al. 2024), although at larger orbital separations of 50–360 au, which might not be directly comparable with β Pic b at 8 au.

We note that fitting only the new GRAVITY spectrum also favours a solar C/O. This solar value is within ~2σ of the value of GRAVITY Collaboration (2020). This difference could be arising from the higher quality of the new spectrum, based on several more epochs of GRAVITY observations. This higher quality seems to be reflected in the better fit obtained on GRAVITY in this work using the same Exo-REM grid in both studies. The addition of MATISSE data also brings tighter constraints on the C/O thanks to the presence of broad H₂O and CO absorption features.

To examine this issue further, we fitted Exo-REM models with low (0.430 ± 0.001) and high (0.650 ± 0.001) C/O Gaussian priors to our GRAVITY+MATISSE datasets. All the other parameters keep the uniform priors used in Sect. 5, and the scale factor is applied only on GRAVITY. The posteriors of these new fits are listed in Table 4, and the models are plotted in Fig. 10. The solar C/O remains favoured over the low and high C/O both for GRAVITY (${\chi }_r^2$ = 2.38 over 2.75 and 2.98, respectively) and MATISSE (${\chi }_r^2$ = 4.28 over 5.07 and 4.45, respectively). Calculating the likelihood ratio as ΔL = Δχ²/2, this means that the solar C/O model is more likely than the low and high C/O models by factors of ~43 and 67 for GRAVITY, respectively, and ~91 and 20 for MATISSE, respectively.

GRAVITY Collaboration (2020) modelled the variation in C/O as function of the accreted mass of solid material, assuming a solar C/O value for the star and either a gravitational collapse or a core accretion scenario. A low C/O value makes the gravitational collapse scenario unlikely within the effective time available for efficient accretion of solid planetesimals during the pre-collapse stage, and favours rather a core accretion scenario between the CO₂ and H₂O icelines. Following their Fig. 6 and 7, if the planetary C/O is solar as found in our study, this would place gravitational collapse back in the possible scenarios with core accretion. Kiefer et al. (2024) notes that a solar C/O could be reached by gravitational instability anywhere in the disc, or by core accretion close to the H₂O ice line with a moderate planetesimal accretion followed by an outward migration. They consider the latter scenario more likely considering that core accretion is preferred for most compact planetary systems with terrestrial planets, and that the β Pic system has at least two planets within 8 au and small km-sized icy bodies (Lecavelier des Etangs et al. 2022). The planetesimal accretion is corroborated by the enhanced metallicity we retrieve ([M/H] = 0.30 ± 0.03).

Fig. 9 Corner plot of the joint fit of GRAVITY (freely scaled) and MATISSE (fixed) spectra, using ForMoSA and the Exo-REM model grid.

Fig. 10 Fit of Exo-REM models using low (pink) and high (blue) C/O priors, compared to the free uniform C/O prior (green) already presented in Fig. 8.

6.2 Other atmospheric parameters

The T_eff (1529 K) and log g (3.84) of our best model are similar to those obtained by previous forward modelling studies that used Exo-REM, while the radius R (~2 R_Jup) is similar or slightly higher: Baudino et al. (2015) (on photometry: 1550 K, 3.5, 1.76 R_Jup), GRAVITY Collaboration (2020) (on GRAVITY and GPI spectra: 1590 K, 4.0, 1.79 R_Jup), Worthen et al. (2024) (on MIRI, GPI, GRAVITY spectra + photometry: 1471 K, 3.71, 1.97 R_Jup). Studies based on forward modelling of ATMO (Phillips et al. 2020) and DRIFT-PHOENIX (Helling et al. 2008) grids, or on radiative transfer retrieval (petitRADTRANS, Mollière et al. 2019) tend, however, to find temperatures > 1700 K and smaller radii of ~1.4 R_Jup (Chilcote et al. 2017; GRAVITY Collaboration 2020; Worthen et al. 2024). The modelling study including the most data to this day (Ravet et al. 2025) found close results to ours (T_eff ~ 1500–1600 K, log g ~ 4.0–4.5) both using Exo-REM and SONORA, while finding high-T_eff (>1800 K) and low-log g (<3.5) solutions using ATMO and BT-Settl. We refer the reader to their paper for a more comprehensive modelling of the SED of β Pic b, and a discussion of the differences between grids of models.

The fitted mass of our best model shifted from an initial prior of 12.7 ± 2.2 M_Jup (from the dynamical mass of GRAVITY Collaboration 2020) to posteriors of 10.5 ± 0.7 (GRAVITY fixed) or 11.8 ± 0.9 M_Jup (MATISSE fixed). These posteriors are well in agreement with the latest dynamical mass estimates using GRAVITY, radial velocity and imaging: ${9.3}_{-2.6}^{+2.5}$ (Brandt et al. 2021) and ${11.90}_{-3.04}^{+2.93}$ M_Jup (Lacour et al. 2021). If we can be relatively confident about the mass, the parameters correlated to it, log g and R (the latest being derived from the two others) are off compared to what is expected from evolutionary models. Based on an estimated age of 24 ± 3 Myr and a bolometric luminosity log L/L_⊙ = −3.76 ± 0.02, and comparing them to Baraffe et al. (2003) hot-start evolutionary models, Chilcote et al. (2017) found expected log g and R of 4.18 ± 0.01 dex and 1.46 ± 0.01 R_Jup, respectively. The reasons for this discrepancy could be both on the model and data sides. On the model side, inconsistencies have already been noted between evolutionary and atmospheric models, as explained in Carter et al. (2023). Additionally, model uncertainties are not estimated (both for evolutionary models and spectral grids), which would result in higher uncertainties to the fitted parameters if provided. On the data side, as explained in Sect. 4.2, the flux calibration of our GRAVITY and MATISSE spectra could be too high, either from overestimated stellar flux or interferometric stellar visibility (if the star is actually resolved). The stellar visibilities measured by Priolet et al. (2025) could help us to derive a more accurate flux calibration.

Fig. 11 Signal-to-noise ratio calculated on the MATISSE contrast spectrum of β Pic b.

6.3 The potential of LM-band exoplanet interferometry

We obtain a high median S/N of ~80 in L band, with maximal values peaking at ~120 around 3.5 μm, as shown in Fig. 11. With this high quality, we expect to be able to characterize fainter and closer-in sub-stellar companions in the future. Separation ranges and L-band fluxes of a sample of planets and brown dwarfs are reported in Fig. 12.

The spectrum of β Pic b in L and M bands presents broad absorption bands of CO and H₂O that can be used to constrain its atmosphere, but it does not show narrow absorption bands or lines, due to its relatively high temperature of ~1500 K. Colder giant planets and brown dwarfs at the L-T transition or below are expected to have deeper absorption bands, including narrow bands of CH₄ at 3.3 μm (seen in early-L to T-type brown dwarfs and young planet analogues, Sorahana & Yamamura 2012; Miles et al. 2023) and CO₂ at 4.2 μm (not accessible from the ground). At medium and high resolution (R > 500), these absorption bands and some spectral lines are resolved. Separating them from each other helps them to be matched to line lists and makes the estimation of molecular abundances easier. These abundances are tracers of the planet birthplace and formation scenario, through C/O and other elemental ratios (Öberg et al. 2011). These abundances can also constrain the disequilibrium chemistry at play in these atmospheres, which becomes particularly important for giant planets below the L–T transition (<1200 K, Charnay et al. 2018; Phillips et al. 2020).

Interferometry is so far the only technique that has obtained medium-resolution spectra of directly detected exoplanets at separations shorter than 0.2″. Like GRAVITY, MATISSE has the potential to observe companions at these separations and provide a complementarity to the JWST, which is not equipped with coronagraphy on its spectroscopic modes and may thus have difficulty getting spectra at these close separations. Ruffio et al. (2024) demonstrated a sensitivity of 3 × 10⁻⁵ at 0.3″ (only 3 spaxels away from the star) for JWST/NIRSpec, but contrast limits below this separation are still unknown, as well as the ability to retain the continuum of the planetary spectrum. For now, companions characterized by interferometry have come mainly from the sample of a few dozen sub-stellar companions discovered by direct imaging. A few interferometric targets have been detected first by radial velocity (Nowak et al. 2020; Hinkley et al. 2023) and astrometric surveys (Gaia DR3, Pourré et al. 2024; Winterhalder et al. 2024). With the fourth Gaia Data Release in 2026, based on more than five years of data, potentially several tens of thousands of long-period Jupiter-mass planets are expected to be discovered (Perryman et al. 2014). Among them, several dozens could be accessible to interferometry (GRAVITY+ Collaboration 2022). Obtaining dynamical masses and KLM fluxes on an order of magnitude more planets with GRAVITY and MATISSE could break the current degeneracy between mass, age, and luminosity in planetary evolutionary models, which results from poor constraints on the post-formation luminosity of giant planets (Mordasini et al. 2017). This will be extremely useful for exoplanet imaging as it relies on a planet’s luminosity to estimate its mass through evolutionary models. This should thus improve the precision on the masses of directly imaged planets, which are often poorly constrained at the moment. It will also better constrain the accretion mechanism in forming planets, which is responsible for the heat accumulated by a planet at the end of its formation.

Finally, the new AO system of the VLTI provided by the GRAVITY+ project, GPAO (Millour et al. 2024), will improve the performances of MATISSE even further in the near future. Switching from MACAO to GPAO, the Strehl ratio is expected to increase from 70% to 90% in L band, which will increase the planet flux injection into the spatial filter of MATISSE by 30%. This better AO correction will in addition reduce the stellar speckle background by a factor of 2 to 5 within the AO correction radius. These two improvements should lead to an increase in planetary S/N by a factor ranging from 1.1× (for planet-to-speckle contrasts ≫ 0.5) to 1.8–2.9x (for planet-to-speckle contrasts ≪0.5, depending on the separation) with GPAO compared to MACAO. These two improvements should lead to a 1.2–2.9× increase in the S/N on the planet with GPAO compared to MACAO, depending on the contrast and separation. Furthermore, the frequency of fringe jumps, which was already well reduced by the latest fringe tracker update (Nowak et al. 2024; Woillez et al. 2024), will decrease even more as result of the higher AO stability, increasing the quality of the MATISSE data. In addition, a dark hole observing strategy (an active speckle suppression technique, Malbet et al. 1995) is being developed on GRAVITY (Pourré et al. 2022). Using the new GPAO system, it could bring the GRAVITY detection limits from a current contrast limit of 4 × 10⁻⁵ at 75 mas down to 3 × 10⁻⁶ at 60 mas (Pourré et al. 2024). A similar strategy could be implemented on MATISSE to reach planets at higher contrasts and closer separations.

Fig. 12 Angular separations and L-band fluxes of a sample of exoplanets and brown dwarfs. Fluxes are either directly taken from the literature (blue points) or computed with ATMO and/or ExoREM models using atmospheric parameters fitted in the literature (red points). Minimal and maximal separations are represented as segments, while the current position of the planet is marked with an asterisk (except for β Pic c and HD 206893 c that are currently well below the MATISSE inner working angle).

7 Conclusions

We observed an exoplanet for the first time with MATISSE using the newly offered GRA4MAT narrow off-axis mode. We developed a new method of correcting chromatic dispersion and non-common path features in the fringe phase. After extraction, we obtained the spectrum of β Pictoris b in the L and M bands at a spectral resolution of 500, showing broad absorption features of H₂O and CO. We used the ForMoSA nested sampling tool and the Exo-REM grid to model the MATISSE spectrum jointly with a new GRAVITY spectrum stacking several years of observations. By imposing a mass prior based on the dynamical mass of β Pic b, we found a best model at T_eff = 1529 ± 3 K, log g = 3.83 ± 0.04, [M/H] = 0.30 ± 0.03, and C/O = 0.539 ± 0.003. This solar C/O value was found both on the fits using only the new GRAVITY spectrum, and the ones using GRAVITY and MATISSE jointly. It is higher than the value found by GRAVITY Collaboration (2020) but in line with Kiefer et al. (2024). This solar C/O indicates that gravitational collapse is not excluded as a formation scenario for β Pic b, although core accretion might still remain favoured by other characteristics of the β Pic system.

The high S/N observed in our spectrum (median of 80 per spectral channel in the L band, with values as high as 120 at 3.5 μm) with only 36 min of integration on the planet indicates that fainter and closer-in companions should be accessible to MATISSE, which opens exciting perspectives for this observing technique and for mid-infrared spectroscopy. MATISSE should be able to complement JWST at short separations at which it cannot obtain spectra due to the absence of coronagraphs on its spectroscopic modes. This new window onto exoplanets comes at an exciting time when both the instrument capability and the exoplanet sample should extend in the near future. The new VLTI AO system, GPAO, has been commissioned since the end of 2024 and is already providing much higher performances than MACAO, thus injecting more planetary flux and less stellar contamination into the VLTI instruments. Finally, in 2026, Gaia DR4 is expected to bring potentially thousands of new exoplanet detections, of which dozens should be accessible to GRAVITY and MATISSE. By providing dynamical mass and KLM-band luminosity measurements on these Gaia companions, future large programs with GRAVITY and MATISSE have the potential to break the current mass-age-luminosity degeneracy in planetary evolutionary models. This will prove extremely useful for future observations with Extremely Large Telescope (ELT) instruments (HARMONI, Thatte et al. 2021; Houllé et al. 2021; METIS, Brandl et al. 2021; Bowens et al. 2021), which will require precise mass estimates to properly interpret their exquisite high-resolution spectroscopic data.

Acknowledgements

We thank the referee and the editor for their helpful comments. This work was supported by the French Agence Nationale de la Recherche (ANR), under grants ANR-21-CE31-0017 (EXOVLTI) and ANR-21-CE31-0018 (MASSIF). M. Bonnefoy acknowledges support in France from the French National Research Agency (ANR) through grants ANR-20-CE31-0012 and ANR-23-CE31-0006. 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 885593). J. Varga is funded from the Hungarian NKFIH OTKA project nos. K-132406 and K-147380, and acknowledges support from the Fizeau exchange visitors program. J. J. Wang is supported by NASA XRP Grant 80NSSC23K0280. The research leading to these results has received funding from the European Union’s Horizon 2020 research and innovation programme under Grant Agreement 101004719 (ORP). This research has made use of the Jean-Marie Mariotti Center Aspro service¹¹, and of the SVO Filter Profile Service “Carlos Rodrigo”, funded by MCIN/AEI/10.13039/501100011033/ through grant PID2020-112949GB-I00.

References

Appendix A Planet versus anti-planet coherent flux ratios

In the February 2023 dataset, we acquired 4 min of data on the so-called ‘anti-planet’ position, i.e. the antipodal point at the same angular separation but with a position angle 180° away from the planet with respect to the star. The objective is to compare the planet signal to the background flux level. In Figure A.1, we show for each baseline the modulus of the ratio between the coherent flux on the planet F_p and on the star F_⋆ (blue curve), and the modulus of the ratio between the coherent flux on the anti-planet point F_antiplanet and on the star F_⋆ (red curve). The complex coherent flux has been averaged over 1 min of exposure time. We expect to see coherent signal from the planet only in the blue curve, and none in the red one. As was expected, coherent flux oscillations typical of a binary target are present at the expected planet position, and absent at its antipodal position.

Fig. A.1 Planet-to-star (blue line) vs anti-planet-to-star (red line) coherent flux ratios in the L and M bands on the six VLTI-UT baselines. Coherent flux oscillations typical of a binary target are present at the expected planet position, and absent at its antipodal position.

Appendix B Phase correction

The measured fringe phase, Φ, is composed of the object fringe phase, Φ_obj, and additional OPD terms produced in the successive mediums along the line of sight: $$\mathrm{\Phi }(\lambda )={\mathrm{\Phi }}_{\mathrm{o}\mathrm{b}\mathrm{j}}(\lambda )+\frac{2\pi }{\lambda }[{\delta }_{\mathrm{v}\mathrm{a}\mathrm{c}}+{\delta }_{\mathrm{a}\mathrm{t}\mathrm{m}}(\lambda )+{\delta }_{\mathrm{D}\mathrm{L}}(\lambda )+{\delta }_{\mathrm{N}\mathrm{C}\mathrm{P}}],$$(B.1)

in which δ_vac, δ_atm, δ_DL, and δ_NCP are the OPDs between telescopes in vacuum, the atmosphere, the delay lines, and the non-common path between GRAVITY and MATISSE, respectively. A schematic view is offered in Fig. B.1. Due to the line of sight inclination relative to the telescope baseline, the light along one of the paths travels, in vacuum, an additional distance, Δd, resulting in the following achromatic vacuum OPD: $${\delta }_{\mathrm{v}\mathrm{a}\mathrm{c}}=d_{\mathrm{v}\mathrm{a}\mathrm{c}}-(d_{\mathrm{v}\mathrm{a}\mathrm{c}}+\mathrm{\Delta }d)=-\mathrm{\Delta }d.$$(B.2)

This additional travel length is then compensated in air in the VLTI delay lines located in the common path, generating a chromatic common-path OPD: $$\begin{array}{rl}{\delta }_{\mathrm{D}\mathrm{L}}& =n_{\mathrm{D}\mathrm{L},1}(\lambda )(d_{\mathrm{D}\mathrm{L}}+\mathrm{\Delta }d)-n_{\mathrm{D}\mathrm{L},2}(\lambda )d_{\mathrm{D}\mathrm{L}}\\ & =\mathrm{\Delta }{n}_{\mathrm{D}\mathrm{L}}(\lambda )d_{\mathrm{D}\mathrm{L}}+n_{\mathrm{D}\mathrm{L},1}(\lambda )\mathrm{\Delta }d.\end{array}$$(B.3)

In these equations and the next, the n_m,i and d_m,i terms designate, respectively, the refractive index and the geometrical length along path i in medium m. Δn_m is the difference of refractive indices between two paths in medium, m. In the approximation of a locally plane-parallel atmosphere, we can consider that light travels the same distance along each path in the atmosphere. The chromatic OPD introduced by the atmosphere is therefore: $${\delta }_{\mathrm{a}\mathrm{t}\mathrm{m}}=[n_{\mathrm{a}\mathrm{t}\mathrm{m},1}(\lambda )-n_{\mathrm{a}\mathrm{t}\mathrm{m},2}(\lambda )]d_{\mathrm{a}\mathrm{t}\mathrm{m}}=\mathrm{\Delta }{n}_{\mathrm{a}\mathrm{t}\mathrm{m}}(\lambda )d_{\mathrm{a}\mathrm{t}\mathrm{m}}.$$(B.4)

The sum of these three OPDs forms the common-path OPD that is seen both by GRAVITY and MATISSE: $$\begin{array}{rl}{\delta }_{\mathrm{C}\mathrm{P}}(\lambda )& ={\delta }_{\mathrm{v}\mathrm{a}\mathrm{c}}+{\delta }_{\mathrm{a}\mathrm{t}\mathrm{m}}(\lambda )+{\delta }_{\mathrm{D}\mathrm{L}}(\lambda )\\ & =(n_{\mathrm{D}\mathrm{L},1}(\lambda )-1)\mathrm{\Delta }d+\mathrm{\Delta }{n}_{\mathrm{D}\mathrm{L}}(\lambda )d_{\mathrm{D}\mathrm{L}}+\mathrm{\Delta }{n}_{\mathrm{a}\mathrm{t}\mathrm{m}}(\lambda )d_{\mathrm{a}\mathrm{t}\mathrm{m}}.\end{array}$$(B.5)

The fringe tracker, GRA4MAT, measures the average of δ_CP(λ) in K band: ⟨δ_CP(λ)⟩_K. It then calculates and sends the according shift to the delay lines to correct the K-band OPD to zero: d_FT = ⟨δ_CP(λ)⟩_K/⟨n_DL(λ)⟩_K. The common-path OPD corrected by the fringe tracker is therefore $$\begin{array}{rl}{\delta }_{\mathrm{C}\mathrm{P}}^{\text{corr\hspace{0.17em}}}(\lambda )& ={\delta }_{\mathrm{C}\mathrm{P}}(\lambda )-n_{\mathrm{D}\mathrm{L}}(\lambda )d_{\mathrm{F}\mathrm{T}}\\ & ={\delta }_{\mathrm{C}\mathrm{P}}(\lambda )-\frac{n_{\mathrm{D}\mathrm{L}}(\lambda )}{{⟨n_{\mathrm{D}\mathrm{L}}(\lambda )⟩}_K}{⟨{\delta }_{\mathrm{C}\mathrm{P}}(\lambda )⟩}_K\\ & =\sum _j({\delta }_j(\lambda )-\frac{n_{\mathrm{D}\mathrm{L}}(\lambda )}{{⟨n_{\mathrm{D}\mathrm{L}}(\lambda )⟩}_K}{⟨{\delta }_j(\lambda )⟩}_K),\end{array}$$(B.6)

in which δ_j is each of the three terms constituting δ_CP in Eq. (B.5).

In order to evaluate the higher limits of these three terms, we simulated them using an air refractive index model and some pessimistic assumptions on the ambient conditions in the lines of sight. We use the air refractive index description of Voronin & Zheltikov (2017) that provides a compact generalization of the Sellmeier equation from the ultraviolet to the far-infrared. We assumed a pessimistic difference of 2° C in temperature and 5% in humidity between two delay lines; and 5° C and 40% between two telescopes’ lines of sight in the atmosphere. We also assumed a maximal length difference between delay lines. The whole set of parameters is listed in Table B.1. We simulated the three common-path OPD terms of Eq. (B.5), corrected by the fringe tracker as in Eq. (B.6). The results are shown in Fig. B.2. We find that (n_DL (λ) − 1) Δd is in general more than 10× larger than the other terms in L and M, except at the beginning of the bands. This term is in addition the easiest and more accurate to model, as sensors are available in the delay line tunnel to get the average delay-line temperature and humidity. The other terms are difficult to model as there is no differential sensing between delay lines and between paths in the atmosphere. We thus neglect them in our corrected common-path OPD model: $$\begin{array}{rl}{\delta }_{\mathrm{C}\mathrm{P}}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}(\lambda )& \approx (n_{\mathrm{D}\mathrm{L}}(\lambda )-1)\mathrm{\Delta }d-\frac{n_{\mathrm{D}\mathrm{L}}(\lambda )}{{⟨n_{\mathrm{D}\mathrm{L}}(\lambda )⟩}_K}{⟨(n_{\mathrm{D}\mathrm{L}}(\lambda )-1)\mathrm{\Delta }d⟩}_K\\ & \approx (\frac{n_{\mathrm{D}\mathrm{L}}(\lambda )}{{⟨n_{\mathrm{D}\mathrm{L}}(\lambda )⟩}_K}-1)\mathrm{\Delta }d.\end{array}$$(B.7)

Our model for the fringe phase is therefore finally $${\mathrm{\Phi }}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}(\lambda )\approx {\mathrm{\Phi }}_{\mathrm{o}\mathrm{b}\mathrm{j}}(\lambda )+\frac{2\pi }{\lambda }[(\frac{n_{\mathrm{D}\mathrm{L}}(\lambda )}{{⟨n_{\mathrm{D}\mathrm{L}}(\lambda )⟩}_K}-1)\mathrm{\Delta }d+{\delta }_{\mathrm{N}\mathrm{C}\mathrm{P}}].$$(B.8)

In our reduction, we only fitted δ_NCP. The ${\delta }_{\text{CP}}^{\text{corr}}(\lambda )$ was modelled thanks to the path lengths, temperature, pressure, and humidity measured in the delay line tunnel and reported in the FITS header of the MATISSE frames. Before fitting, we additionally remove from the data and model the mean L-band phase to get the differential phase, a standard step in interferometry as the absolute phase reference is lost during fringe tracking: $${\mathrm{\Phi }}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}(\lambda )={\mathrm{\Phi }}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}(\lambda )-{⟨{\mathrm{\Phi }}^{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}(\lambda )⟩}_L.$$(B.9)

The final estimate of the object differential phase is the difference between the differential phase data and the fitted OPD model: $${\mathrm{\Phi }}_{\mathrm{o}\mathrm{b}\mathrm{j}}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}(\lambda )\approx \arg \left[{\mathrm{e}}^{\mathrm{i}{\mathrm{\Phi }}^{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}(\lambda )}{\mathrm{e}}^{-\mathrm{i}\frac{2\pi }{\lambda }({\delta }_{\mathrm{C}\mathrm{P}}^{\text{corr\hspace{0.17em}}}(\lambda )+{\delta }_{\mathrm{N}\mathrm{C}\mathrm{P}})}\right].$$(B.10)

We show in Fig. B.3 the non-common path OPDs fitted in each frame and each baseline. We can see systematic differences between BCD configurations, and a slow linear drift during the two hours, highlighted by the linear fits on one of the BCD configurations. The fitted OPD drifts range from 0.8 to 7.0 μm/h depending on the baseline.

Fig. B.1 Schematics of a MATISSE + GRA4MAT observation, with the terms described in equations 1 to 5.

Fig. B.2 Simulation of the common-path OPD terms. The vertical dashed line shows the fringe tracking wavelength (2.2 μm), where the OPD is corrected to zero (horizontal dashed line).

Fig. B.3 Temporal evolution of the non-common path OPD fitted on the data. The colours represent different BCD positions. The dashed line shows a linear fit on the frames using the IN-IN BCD position (blue crosses). The fit parameters and coefficient of determination are printed in each subplot.

Appendix C Additional original datasets

Appendix C.1 GRAVITY

The new GRAVITY spectrum used in this paper has been built by taking the covariance-weighted mean of the contrast spectra of five β Pic b epochs from 2018 to 2022. The data was reduced with the exogravity pipeline (GRAVITY Collaboration 2020). The planet spectra of all epochs are shown in Fig. C.1, as well as the covariance-weighted average spectrum. We scaled the GRAVITY average contrast spectrum with the same stellar model as the one used for MATISSE.

Appendix C.2 NACO

β Pic b was observed on December 5, 2011 (programme 088.C-0196(A), PI S. Quanz) with VLT/NACO (Lenzen et al. 2003; Rousset et al. 2003), using the 172 mas slit, the L27_1_SL spectroscopic mode, the detector cube mode, and the HighDynamic / Double_RDRSTRD readout. These settings provide a wavelength coverage from 2.6 to 4.2 μm at a spectral resolution of ~350.

The slit was oriented so that both the central star and the planet (located at a separation of 450 mas and a position angle of 211.6° at the time of observation) were included. Nodding with an ABBA pattern was used to subtract background, but no extra jitter to ensure a maximum relative positioning stability. Star recentring was employed roughly every 40 minutes. A 180° camera rotation was also employed to reduce further instrumental artifacts such as ghosts. In all exposures the peak flux of the central star was kept below the linearity limit, roughly two thirds of the full dynamic range, including at the longest wavelengths (~4.1 μm) where the background emission dominates.

Each sequence of data was reduced separately using custom-made IDL scripts. A careful frame subtraction and spectrum extraction was then applied to the data cubes, as summarized hereafter. First, cubes for which the AO loop was opened and/or the Strehl ratio was low were rejected using statistical analysis. Then, sub-pixel shifts of the star position due to AO jitter were corrected using cross-correlation and sub-pixel shifts with a tanh interpolation kernel. Finally, positions of 17 telluric absorption features from ESO’s online DR05 Sky model (v1.2.0) were measured on the β Pic spectra and compared to the models to derive a wavelength calibration of the spectra. The result is a data cube containing the star and planet spectra along the slit in a 2D spatial-spectral map.

In order to subtract stellar contamination and extract the planet spectrum, we then applied a custom version of the so-called ‘spectral deconvolution’ algorithm (Sparks & Ford 2002; Vigan et al. 2008, 2012) to suppress most of the flux arising from the star. The algorithm uses the chromatic dependency of the stellar PSF and all its substructures (speckles and quasi-static aberrations) through the slit with respect to the fixed position of a companion to deblend the spectra of the two objects: each spectral channel i is rescaled by a factor λ_m/λ_i with λ_m the longest wavelength covered by the spectrum and λ_i the processed channel one. The rescaled slit is then composed of a series of nearly identical spectra of the star that can be median-combined to model the star spectrum while the companion spectrum is moving diagonally through the slit. This model spectrum is flux-normalized and subtracted from each stellar spectrum, leaving only the companion spectrum. As a final step, the spectral channels of the residual companion spectrum are scaled back to their original spatial resolution.

The level of residuals is strongly correlated to the position of the star within the slit, as the third Airy ring contaminates the planet spectrum at wavelengths longer than 3.5 μm for some stellar positions. We selected all non-contaminated spectra and extracted separately the planet spectrum in each. Spectra of β Pic b were divided by the stellar spectrum and multiplied by a blackbody function at the temperature of the star (8052 K, Gray et al. 2006) to remove telluric features. The averaged spectrum is shown in Figure 6. We masked zones with strong telluric residuals or where the extraction process left strong residuals.

Fig. C.1 Compilation of GRAVITY spectra obtained at different epochs on β Pic b. The spectrum used in this paper is the weighted mean of all these spectra.

All Tables

All Figures

Fig. 1 Simulation of the gain in contrast achieved by offsetting the pinhole of MATISSE on β Pic b, separated here by 534 mas from its host star. The upper and lower plots depict the stellar and planetary contribution to the coherent flux (the fringe envelope) as seen through the circular MATISSE pinhole (white area, hence the Airy-like shape of the output flux) centred on the star and on the planet, respectively. They are simulated at a wavelength of 3.5 μm and normalized to the maximal intensity of the star-centred stellar PSF. When the pinhole is centred on the planet, the speckle flux is reduced to the same level as the planetary flux. In each simulation, the asymmetric shape of the PSF of the offset object (relative to the pinhole) is due to imperfect spatial filtering of the Airy rings falling in the pinhole.

In the text

Fig. 2 Measured differential phase (grey dots) in two MATISSE frames, one on the star (top) and the other on the planet (bottom). Overlaid are the models presented in this work using the air refractive indices of Voronin & Zheltikov (2017) (blue), Mathar (2007) (green), and Ciddor (1996) (red). In the absence of chromatic or non-common path OPD, the differential phase should be zero on the star or an oscillation on the planet.

In the text

Fig. 3 Real (left) and imaginary (right) parts of the measured planet-to-star coherent flux ratio, ℱ, on the six UT baselines, during one 10 s frame. The error range is plotted in the background. The thick lines show the best-fit model as described in Sect. 4.1. It includes the star–planet modulation, the stellar speckle contamination (dashed blue curves), and the planet-to-star contrast assumption (dashed black curve).

In the text

	Fig. 4 Scaling of the BT-Nextgen stellar model (Teff = 7890 K, log g = 3.83, [Fe/H] = 0.0) with species using Tycho, Gaia, and 2MASS photometry of β Pictoris. The bottom panel shows the residuals of the fit.
In the text

Fig. 5 Detection maps of β Pic b with MATISSE. The maps show the reduced χ2 values after fitting the stellar contamination for each tested planet astrometry (Δα, Δδ) in the grid. The grids are centred on MATISSE’s pointing during the planet observations. The left map is the size of the pinhole (1.5λ/D at 3.5 μm) and has a resolution of 1 mas. The right map was generated at higher resolution (0.1 mas) on a smaller 20 × 20 mas grid. The most likely astrometry from the MATISSE observations is shown as a blue ellipse, found by fitting a 2D Gaussian curve on the lowest ${\chi }_r^2$ peak. The planet location and uncertainties predicted by whereistheplanet from the orbital fitting of GRAVITY observations (Lacour et al. 2021) are shown as a green ellipse. The multiplicity of ${\chi }_r^2$ peaks is due to the interferometric nature of our data. As coherent flux is periodic as a function of α.u (see Eq. (7)), an infinity of astrometric solutions can reproduce the signal. This degeneracy gradually disappears when combining data from different baseline lengths and orientations. The ${\chi }_r^2$ of neighbour peaks can nonetheless still remain close. Using covariances between baselines and wavelengths could help increase the difference between peaks, and will be studied in a future work.

In the text

Fig. 6 Spectrum of β Pic b with MATISSE (black, this work) compared to other instruments: NACO (yellow, new data), GPI (blue, Chilcote et al. 2017), GRAVITY (green, new data) and MIRI (red, Worthen et al. 2024). Only the slope of the NACO spectrum can be compared to MATISSE as the absolute level was scaled manually. MATISSE data with a S/N < 5, located in strong telluric bands and not used in this work, are shown in grey. The absorption curves of H2O, CO, and CH4, calculated from the best-fitting Exo-REM model presented in Sect. 5, are plotted at the top to visualize their impact on the spectrum. H2O and CO leaves a visible pattern, but not CH4, likely because it is dominated by H2O absorption or buried in noise.

In the text

Fig. 7 Covariance (upper right triangle) and correlation (lower left triangle) matrices of the MATISSE spectrum of β Pic b. In order to highlight the off-diagonal correlations in the L band, the covariance values are clipped between 3 × 10−35 and 3 × 10−33 W2 m−4 μm−2 and the correlation values between 0.01 and 1. Two lines of off-diagonal correlations are further highlighted with arrows for readability.

In the text

	Fig. 8 GRAVITY and MATISSE spectra of β Pic b with the best-fit Exo-REM models. The top plot shows the best fits on GRAVITY only, MATISSE only, and GRAVITY+MATISSE, when no scaling is applied on the spectra. The two other plots show the best fits on GRAVITY+MATISSE when a scale factor is applied either on MATISSE (middle) or on GRAVITY (bottom).
In the text

	Fig. 9 Corner plot of the joint fit of GRAVITY (freely scaled) and MATISSE (fixed) spectra, using ForMoSA and the Exo-REM model grid.
In the text

	Fig. 10 Fit of Exo-REM models using low (pink) and high (blue) C/O priors, compared to the free uniform C/O prior (green) already presented in Fig. 8.
In the text

	Fig. 11 Signal-to-noise ratio calculated on the MATISSE contrast spectrum of β Pic b.
In the text

Fig. 12 Angular separations and L-band fluxes of a sample of exoplanets and brown dwarfs. Fluxes are either directly taken from the literature (blue points) or computed with ATMO and/or ExoREM models using atmospheric parameters fitted in the literature (red points). Minimal and maximal separations are represented as segments, while the current position of the planet is marked with an asterisk (except for β Pic c and HD 206893 c that are currently well below the MATISSE inner working angle).

In the text

	Fig. A.1 Planet-to-star (blue line) vs anti-planet-to-star (red line) coherent flux ratios in the L and M bands on the six VLTI-UT baselines. Coherent flux oscillations typical of a binary target are present at the expected planet position, and absent at its antipodal position.
In the text

	Fig. B.1 Schematics of a MATISSE + GRA4MAT observation, with the terms described in equations 1 to 5.
In the text

	Fig. B.2 Simulation of the common-path OPD terms. The vertical dashed line shows the fringe tracking wavelength (2.2 μm), where the OPD is corrected to zero (horizontal dashed line).
In the text

	Fig. B.3 Temporal evolution of the non-common path OPD fitted on the data. The colours represent different BCD positions. The dashed line shows a linear fit on the frames using the IN-IN BCD position (blue crosses). The fit parameters and coefficient of determination are printed in each subplot.
In the text

	Fig. C.1 Compilation of GRAVITY spectra obtained at different epochs on β Pic b. The spectrum used in this paper is the weighted mean of all these spectra.
In the text