© GRAVITY Collaboration 2020

1 Introduction

The ever-increasing number of exoplanet detections (over 4000, at the time of this writing¹) proves that our instrumental capabilities are getting better and better at discovering these other worlds. But even though exoplanets are now routinely being observed, determining their physical properties (temperature, mass, composition), let alone the history of their formation, remains extremely challenging. And yet, these measurements are key to understanding the details of planetary formation processes.

Among all measurable quantities, element abundance ratios are emerging as some of the most promising for understanding planetary formation. The question of the supersolar abundances of heavy elements in the atmosphere of Jupiter is probably what motivated the first attempts to link abundance ratios to planetary formation, and several studies have been carried out to understand how planetesimal accretion can lead to heavy element enrichment (Helled et al. 2006; Helled & Schubert 2009; Owen et al. 1999; Alibert et al. 2005). On the exoplanet front, the work of Öberg et al. (2011) was the first general attempt to show that element ratios in an exoplanet atmosphere can be an imprint of its formation history. This idea has since been investigated further by several authors (e.g., Ali-Dib et al. 2014; Thiabaud et al. 2014; Helling et al. 2014; Marboeuf et al. 2014a,b; Madhusudhan et al. 2014, 2017; Mordasini et al. 2016; Öberg & Bergin 2016; Cridland et al. 2016; Eistrup et al. 2016, 2018). While Öberg et al. (2011) highlighted how gas disk abundances can influence the atmospheric composition, the importance of icy planetesimals for the atmospheric enrichment is stressed in Mordasini et al. (2016), where exoplanet spectra are derived from modeling full formation in the core-accretion paradigm.

Measuring the element ratios is not easy, and requires high-quality data. Madhusudhan et al. (2011) used a free retrieval method on a set of Spitzer and ground-based photometric data in 7 different bands to obtain the first exoplanetary C/O ratio on the hot Jupiter WASP-12b. But the value of C∕O > 1 they obtained has since been ruled out by Kreidberg et al. (2015), showing the difficulty of obtaining reliable abundance ratios. Konopacky et al. (2013) used a different approach in their study of HR 8799 c. They obtained K-band spectroscopic observations of the planet with the spectrograph OSIRIS on the Keck II telescope, and were able to extract an estimate of the C/O ratio using model grid fitting. They found a value of C∕O = 0.65 ± 0.15. Looking at the same planetary system, Lavie et al. (2017) estimated the C/O ratio for four planets (HR 8799 b, c, d, and e), using a retrieval analysis method. In their analysis, they notably emphasized the importance of high-quality K-band spectroscopic data, which they found to be critical for a reliable measurement of the C/O and C/H ratios.

With the recent direct detection of the giant planet HR 8799 e with the GRAVITY instrument (Gravity Collaboration 2019) on the Very Large Telescope Interferometer (VLTI), optical interferometry has become a new arrow in the quiver of exoplanet observers. By taking advantage of the angular-resolution offered by 100+ meter baselines, optical interferometers can separate a dim exoplanet from the overwhelming residual starlight, leading to accurate measurements of the astrometric position (up to 10 μas, Gravity Collaboration 2018), and high signal-to-noise spectroscopic data with absolute calibration of the continuum.

In this paper, we present observations of the giant planet β Pic b obtained with GRAVITY and we investigate the possibility of using this K-band spectro-interferometric data to determine the C/O ratio of the planet. The observations are presented in Sect, 2, together with a brief summary of the data reduction (a complete explanation is given in Appendix A). Section 3 focuses on the orbit and mass of β Pic b. We show in this section how the new GRAVITY astrometric data impacts the best orbital estimate currently available and we provide a new estimate of the dynamical mass of the planet. Section 4 is devoted to the measurement of atmospheric properties and, in particular, to the determination of the C/O ratio, using two different approaches: forward modeling and free retrieval. In Sect. 5, we discuss the C/O ratio obtained in the case of a formation of β Pic b through gravitational accretion and then through core-accretion. Our general conclusions can be found in Sect. 6.

2 Observations and data reduction

2.1 Observations

Observations of β Pictoris b were obtained on September, 22, 2018, using the GRAVITY instrument (Gravity Collaboration 2017), with the four 8 m Unit Telescopes (UTs) of the VLT. The instrument was set up in its medium resolution mode (R = 500), and observationswere conducted in on-axis/dual-field mode.

The observing strategy was similar to the one described in Gravity Collaboration (2019): the fringe-tracker (Lacour et al. 2019)was using the flux from the central star during the observing sequence, while the position of the science fiberwas changed at each exposure, alternating between the central star and the position of the planet. Since the planet was not visible on the acquisition camera, the position used to center the fiber during the planet exposures was a theoretical position, based on predictions from previous monitoring (Wang et al. 2016; Lagrange et al. 2018).

A total of 16 exposures (resp. 17) were acquired on the star (resp. the planet). Each star exposure was made of 50 individual 0.3 s integrations. For the planet, which is ~10 mag fainter thanthe star, the integration time was initially set to 30 s, with 10 integrations per exposure, and reducedto 10 s with 30 integrations at mid-course, since the observing conditions were good (seeing < 0.8″). The complete dataset contains 1.4 h of integration on the planet (and 0.35 h of associated background exposures), and 4 min 30 s of integration on the central star (plus 1 min 15 s of sky background). The observing log is given in Table 1.

2.2 General data reduction

During planet exposures, the science fiber at each telescope is kept at an offset position with respect to the star, reducing significantly the star to fiber coupling ratio. But even though most of the stellar flux is rejected, speckle noise can still couple to the science fiber and dominate the exposures, hence the need for careful data reduction to disentangle the planet signal from the remaining coherent stellar flux.

The general data reduction method used to reduce the VLTI/GRAVITY observations of β Pic b is presented in details in Appendix A. It can be divided into different parts: pipeline reduction (common to all GRAVITY observations), astrometric extraction, and spectrum extraction. These steps are described in Appendices A.2, A.4, and A.5. The end products are an astrometric position for the planet with respect to the star (Δα, Δδ), and a planet-to-star contrast spectrum C(λ) = S_P(λ)∕S_⋆(λ) which is the ratio between the spectra of the planet and of the star.

2.3 K-band spectrum

The contrast spectrum of β Pic b was converted to an absolute spectrum of the planet using a model of the stellar spectrum: S_P (λ) = C(λ) × S_⋆(λ). We used a BT-NextGen model (Hauschildt et al. 1999), with a temperature of 8000 K, a surface gravity of log (g∕g₀) = 4, and a Solar metallicity, as close as possible to the measured value for this star (Lanz et al. 1995; Gray et al. 2006). We scaled this synthetic spectrum to an ESO K-band magnitude of 3.495, taking into account the correct filter (van der Bliek et al. 1996). This strategy, based on the extraction of a contrast spectrum and the use of a model for the star, helps to reduce the impact of Earth’s atmosphere on the final planet spectrum. The result is given in Fig. 1.

2.4 Astrometry

Using the data reduction method described in Appendix A.4, we found a mean relative planet to star astrometry on all the exposure files of: $$\{\begin{array}{l}\Delta \text{RA}=68.48 \text{mas}\hfill \\ \Delta \text{Dec}=126.31 \text{mas}.\hfill \end{array}$$(1)

The 1 σ confidence interval is given by the covariance matrix of all the 17 exposure files: $$\begin{array}{ll}\text{Covar}\left(\Delta \text{RA},\Delta \text{Dec}\right)\hfill & =\left[\begin{array}{cc}0.0027& -0.0035\\ -0.0035& 0.0045\end{array}\right] {\text{mas}}^2\hfill \end{array}$$

This GRAVITY measurement is shown in the inset plot of Fig. 2.

In its dual-field mode, GRAVITY is limited to observations of planets above the diffraction limit of a single telescope (to separate the planet from the central star), but the relative astrometry derived from these observations still fully benefits from the length of the telescope array.

Fig. 1 Calibrated K-band spectrum of β Pictoris b, atR = 500, extracted from the VLTI/GRAVITY observations (gray points). For comparison, the K-band part of the GPI spectrum from Chilcote et al. (2017) (R ≃ 70) is also overplotted (orange points). The error bars plotted for the GRAVITY spectrum only represent the diagonal part of the full covariance matrix.

3 Orbit and dynamical mass

3.1 Orbital parameters

We fit a Keplerian orbit to the visual astrometry of the planet to characterize its dynamics. As our new GRAVITY point is more than an order of magnitude more precise than any other published astrometric point on the northeastern half of its orbit (c.f., Lagrange et al. 2019a), we expected a better constraint on the eccentricity of the planet’s orbit. We used the published astrometry from Chauvin et al. (2012), Nielsen et al. (2014), and Wang et al. (2016) in this analysis. The orbit was fit using the open-source Python orbit fitting package orbitize! (Blunt et al. 2019). We included a custom likelihood to fit the GRAVITY measurement along the two principal axes of the error ellipse. We fit for the same eight parameters as Wang et al. (2016): semi-major axis (a), eccentricity (e), inclination (i), argument of periastron (ω), position angle of the ascending node (Ω), the first periastron passage after MJD = 55 000 in units of fractional orbital period (τ), system parallax, and total systemmass (M_tot). We generally used relatively unconstrained priors for most of the orbital parameters (see Table 2). For Ω, we constrained it to between π∕10 and π∕2 to account for the fact that Snellen et al. (2014) detected the RV signal of the planet. However, we chose not to explicitly include the RV in the fit as there could be systematics in the reported uncertainties. For the parallax, we used a normal distribution to represent the parallax of 51.44 ± 0.12 mas measured by HIPPARCOS (van Leeuwen 2007). We sampled the posterior using the parallel-temperature affine-invariant sampler in ptemcee (Foreman-Mackey et al. 2013; Vousden et al. 2016) with 20 temperatures, 1000 walkers per temperature. We discarded the first 15 000 steps to allow the walkers to converge. We assessed convergence using the autocorrelation time and by visual inspection of the samples. We then ran each walker for 5000 steps, keeping only every tenth sample to mitigate correlations in the samples produced by any given walker.

Our constraints on the orbit of β Pic b using just astrometry of the planet are collected in Table 2 and plotted in Fig. 2. We find that < 2% of allowed orbits have e < 0.05 and <0.5% of orbits have e < 0.03, although there are still some allowed circular orbits. Dupuy et al. (2019) also proposed an e ≈ 0.25 when including astrometric and radial velocity data on the system. To statistically assess whether eccentric orbits are preferred, we refit the orbit fixing e = 0 and ω = 0 resulting in a fit with two less parameters. Similar to Wang et al. (2018) in assessing the coplanarity of the HR 8799 planets, we compared the Bayesian Information Criterion (BIC) of the fit that allowed eccentric orbits with the fit that fixed the orbit to be circular, and found that the BIC disfavors the circular orbit by 9.9. The reduction in model parameters for a purely circular orbit does not compensate for an increase in fitting residuals, so we disfavor circular orbits for a single planet model. However, additional confusion on this measurement could be due to a second planet in the system (Lagrange et al. 2019b). The second planet β Pic c would induce epicycles in the apparent orbit of β Pic b around the star due to the gravitational influence of the second planet on the orbit of the host star. Using parameters for β Pic c from Lagrange et al. (2019b), the magnitude of these epicycles are several hundred μas, so well detectable by GRAVITY, but hidden beneath the uncertainty of previous astrometry. Thus, they would also bias this single GRAVITY measurement, and continued astrometric monitoring is required to separate out the signal of the separate planet from a possibly eccentric orbit of β Pic b.

However, a moderate eccentricity would fit nicely in the dynamics of the system. An e ≈ 0.15 is consistent with the picture of an eccentric β Pic b launching small bodies towards the star, causing spectroscopic and transiting signatures of exo-comets in observations of the star (Thébault & Beust 2001; Zieba et al. 2019). An interesting question is how such a massive planet acquired a significant eccentricity. The obvious conclusion would point to a second massive planet in the system, such as the radial velocity detected β Pic c (Lagrange et al. 2019b). Otherwise, Dupuy et al. (2019) proposed that if the planet had formed further out and migrated inwards, resonant interactions with the circumstellar disk could pump up its eccentricity to the values we observe today. Characterizing the detailed structure of the circumstellar dust in the system as well as the chemical composition of β Pic b could test this theory.

Generally, the other orbital parameters of β Pic b have already been sufficiently well constrained previously that out results agree with the conclusions drawn in previous works (Millar-Blanchaer et al. 2015; Wang et al. 2016; Lagrange et al. 2019a; Dupuy et al. 2019). We still find that the planet did not transit the star in 2017, and that the Hill sphere of the planet did transit. Assuming a planet mass of 12.9 ± 0.2 M_Jup, we find a Hill sphere ingress at MJD 57852 ± 2 (2017 April 8) and a Hill sphere egress at MJD 58163 ± 2 (2018 February 13). The closest approach, which does not require an assumption on the planet’s mass, is at MJD 58008 ± 1 (2017 September 11), with the planet passing 8.57 ± 0.13 mas from the star (0.166 ± 0.003 au in projection). The precise astrometry of the GRAVITY epoch post conjunction has significantly improved the transit ephemeris from Wang et al. (2016).

Fig. 2 Visualorbit of β Pic b. Plotted in black are possib le orbits randomly drawn from the posterior using only relative astrometry (Sect. 3.1). Previous astrometric measurements used in the orbit fit are in blue. The GRAVITY measurement from this work is in red, with an inset plot that is zoomed in by a factor of ~2000 to display the uncertainties on this measurement.

3.2 Dynamical mass determination

A significant astrometric acceleration for the star β Pic was detected when comparing its average velocity over the course of the HIPPARCOS mission and the average velocity inferred by the change in position of the star between the HIPPARCOS and Gaia missions (Snellen & Brown 2018; Kervella et al. 2019). Assuming this acceleration is due entirely to β Pic b, Snellen & Brown (2018) and Dupuy et al. (2019) used it in conjunction with the visual orbit to measure a dynamical mass for the planet. Snellen & Brown (2018) used the HIPPARCOS intermediate astrometric data (IAD; van Leeuwen 2007) and Gaia DR2 position (Gaia Collaboration 2018) to fit the position, proper motion, and orbital motion of the host star to derive the mass of the planet. Dupuy et al. (2019) used the re-calibrated HIPPARCOS and Gaia proper motions from the HIPPARCOS/Gaia Catalog of Accelerations (HGCA; Brandt 2018) and the stellar radial velocities from Lagrange et al. (2012) to derive the mass of the planet. Being agnostic to which method is more accurate, we repeated both analyses here, now with the new GRAVITY epoch providing strong constraints on a and e, which are otherwise degenerate with the mass of β Pic b, M_b. To repeat the Snellen & Brown (2018) orbit fit, we include five additional parameters in the fit: the position and proper motion of the star in RA and Dec, as well as the mass of the planet. We also switch the prior on parallax to a uniform prior between 50.24 and 52.64 mas, since the HIPPARCOS intermediate astrometric data now constrains this parallax. To repeat the Dupuy et al. (2019) analysis, we only fit for changes in the tangential velocity of the host star, so we do not need to fit for its actual position and proper motion. We only include a RV offset and RV jitter term for the stellar RV data. We modified the orbitize! custom likelihood function to include these measurements of the host star, and repeated the orbit fit.

We list the orbital and mass constraints in Table 2, marginalizing over astrometric parameters of the host star and stellar RV calibration numbers in the two fits. We also plot the posterior probabilties for the mass of β Pic b in Fig. 3. In the fit using the HIPPARCOS IAD, the semi-major axis and eccentricity posteriors now favor slightly higher values by 1σ. We find a dynamical mass of β Pic b of 12.7± 2.2 M_Jup, which is consistent with the values from Snellen & Brown (2018). Conversely, using the recalibrated stellar astrometry from the Brandt (2018) HGCA catalog and the stellar RVs, we find a slightly lower semi-major axis and eccentricity by 1σ than the relative astrometry only fit. Despite these minor differences, all three fits considered in this work favor an eccentricity between 0.1 and 0.2. We do not find orbital solutions with e > 0.25 as has been suggested by Dupuy et al. (2019). In the HGCA fit, we also find a weaker dynamical mass constraint for β Pic b of ${14.2}_{-3.9}^{+3.7}$ M_Jup, which is consistent with Dupuy et al. (2019). The HIPPARCOS IAD method provides more stringent constraints on the planet mass, likely because it has smaller uncertainties. It is unclear whether this better constrain is unbiased, or if the uncertainties are underestimated due to calibration systematics or effects of other planets on the stellar astrometry. However, as seen in Fig. 3, both fits agree with each other, and both dynamical masses are consistent with hot-start derived masses of 12.7± 0.3 M_Jup from Morzinski et al. (2015) and 12.9± 0.2 M_Jup from Chilcote et al. (2017). More accurate stellar astrometry or RVs are necessary to test hot-start evolutionary models more stringently, given that the model-dependent hot-start masses have an order of magnitude better precision than the dynamical masses.

Fig. 3 Dynamical mass estimates of β Pic b using the two different methods described in Sect. 3.2. The shaded grey region is the 2σ uncertainty on the hot-start derived mass from Chilcote et al. (2017).

4 The atmosphere of β Pic b

4.1 Previous work

Physical parameters of β Pictoris b have been reported in a number of previous studies (see Table 15 of Morzinski et al. 2015, Table 2 in Chilcote et al. 2017 for a summary of these results). The temperature of the planet has been estimated by several authors, using atmospheric or evolutionary model grid fitting, on photometric and/or spectroscopic data. The most extensive study to date was performed by Chilcote et al. (2017), who obtained GPI spectroscopic data at R ≃ 50 in Y, J, H, and K-band, as well as photometric points in different bands ranging from 1 to 5 μm. Using different atmospheric models (BT-Settl: Allard et al., 2012; DRIFT-PHOENIX: Woitke & Helling 2003, Helling & Woitke 2006; AMES-DUSTY: Chabrier et al. 2000; Allard et al. 2001), they obtained values ranging from 1650 to 1800 K for the temperature, and 3.0 to 4.5 for log(g∕g₀). These values are similar to what is reported in Bonnefoy et al. (2013, 2014), Chilcote et al. (2015), and Morzinski et al. (2015), with the same models.

The lower limit of the range of temperatures estimated comes from Baudino et al. (2015). Using their Exo-REM model grid, and a set of photometric data only (the GPI spectrum was not available at the time), they derived a temperature of 1550 K, and a surface gravity of log(g∕g₀) = 3.5.

4.2 ExoREM atmospheric grid fitting

Using either the GRAVITY K-band spectrum only, or the GRAVITY K-band and GPI Y JH bands spectra, we performed a grid model fitting using the newest ExoREM grid (Charnay et al. 2018).

We performed a χ²-based grid model fitting on the GRAVITY K-band only data, using the same ExoREM model grid as used to fit the GRAVITY HR 8799 e spectrum in Gravity Collaboration (2019), ranging from 400 to 1800 K in temperature, with a step-size of 50 K, from 3.0 to 5.0 in log (g∕g₀), with a step-size of 0.2, for a metalicity of [Fe∕H] = −0.5, 0, and 0.5, and with a Solar C/O ratio. The best fit was obtained for a Solar metallicity, a temperature of 1750 K, and a log (g∕g₀) of 3.30. However, this best fit also leadsto a mass of 1.3 M_Jup, more than 5 σ away from our estimate given in Table 2. To force the result of the fit to be in agreement with our mass estimate, we added a mass prior in the χ² calculation. We used a weight for the prior similar to the weight of the entire GRAVITY spectrum: $${\chi }^2=n_{\lambda }\frac{{(m-12.7M_{\text{Jup}})}^2}{{(2.2M_{\text{Jup}})}^2}+{\displaystyle \sum \frac{{\left(F_{\text{data}}({\lambda }_k)-F_{\text{model}}({\lambda }_k)\right)}^2}{{\sigma }_F{({\lambda }_k)}^2}}$$(2)

in which F_data and F_model represent the flux from the data and from the model at the different wavelengths, σ_F the error on the data, and m the mass derived from the flux level.

With this new definition of the χ², the same ExoREM grid led to a best fit at T = 1500 K, log (g∕g₀) = 4.0, for a Solar metallicity. The corresponding planet radius is 1.9 R_Jup, and the mass is 14 M_Jup, compatible with the estimate of Sect. 3. However, we find that the fit itself was not very good, with a ${\chi }_{\text{red}}^2$ value of 6.8. The CO region around 2.3 μm was particularly poorly fitted.

The fit was improved by generating a second ExoREM grid, which included the C/O ratio as an additional parameter. The grid was generated on the same range of temperature, surface gravity, and metallicity, for C/O values ranging from 0.3 to 0.8, with a step of 0.05.

Without the mass prior, the new grid yielded a best fit corresponding to a temperature of 1700 ± 50 K, a surface gravity of log(g∕g₀) = 3.5, a metallicity of −0.5 (the lowest value available in our grid), and a C/O ratio < 0.3. The resulting planet mass remained too low, at 2 M_Jup. Adding the mass prior in the definition of the χ², as in Eq. (2) led to a temperature of 1550 ± 20 K, a surface gravity log(g∕g₀) of 4.0, a metallicity of [Fe∕H] = 0.5 (highest value from our grid), and a C/O ratio of 0.41 ± 0.05, for a planet mass of 11.5 M_Jup, in very good agreement with the result of Sect. 3. Contraining the fit to Solar metallicity resulted in a very low C/O ratio of 0.3, with similar temperature and surface gravity.

Including the GPI Y, J, and H band data from Chilcote et al. (2017) and allowing for a multiplicative scaling factor between GPI and GRAVITY resulted in a temperature of T =1590 ± 20 K, with a C/O of 0.43 ± 0.05, for a metallicity of [Fe∕H] = 0.5. For reference, the typical multiplicative factors needed to scale the GPI spectra on the ExoREM grid were ≃ 0.85 for the Y band, and ≃0.9 for J and H bands.

The results of these different fits are summarized in Table 3, and the best fit obtained using GRAVITY+GPI and a mass prior is shown in Fig. 4.

Fig. 4 Best fit obtained with the ExoREM atmospheric model (Charnay et al. 2018) using GPI Y, J, H + GRAVITY K bands, and a mass prior.

4.3 Free retrieval with petitRADTRANS

4.3.1 Retrieval forward model

In addition to fitting a model grid to the β Pic b observation, we carried out a free retrieval. To this end, the spectra were compared to the predictions of a spectral synthesis code, where the atmospheric structure was parametrized. In such an approach more weight is given on atmospheric conditions as constrained by the data, while principles such as radiative-convective equilibrium do not have to be strictly fulfilled. This approach was motivated by the work of Line et al. (2015, 2017); Zalesky et al. (2019) for clear, and Burningham et al. (2017) for cloudy brown dwarfs, in which the power of free retrievals to constrain condensation and cloud processes has been demonstrated.

Our “forward model”, used for predicting the spectra, was constructed using petitRADTRANS (Mollière et al. 2019). Because the atmosphere of β Pic b is expected to be cloudy, we added scattering to petitRADTRANS. We verified the calculations by comparing to spectra of self-consistent models for cloudy, self-luminous planets obtained with petitCODE (Mollière et al. 2015, 2017), which agreed excellently.

One benefit of using a free retrieval is that one of the most uncertain physical processes, namely the formation of clouds, can be parametrized, letting the observations constrain the cloud mass fraction and particle size distribution. A related approach was taken by Burningham et al. (2017), who carried out free retrievals for cloudy brown dwarfs for the first time. Here, we assume that our clouds consist of iron and silicate particles, which fixes the location of the cloud base for a given temperature profile. The cloud parameterization of Burningham et al. (2017) was even more general. One of them retrieved the cloud location (where it becomes optically thick), scale height, the single scattering albedo, as well as the power law slope of the opacity.

For the fits presented here, we parametrized the clouds using the Ackerman & Marley (2001) cloud model. However, in contrast to the usual treatment in grid models (see, e.g. Ackerman & Marley 2001; Marley et al. 2012; Morley et al. 2014; Mollière et al. 2017; Samland et al. 2017; Charnay et al. 2018), we retrieved all of its three parameters. First, the settling parameter f_sed, which is the mass-averaged ratio between the settling and mixing velocity of the cloud particles. This determines the decrease of the cloud mass fraction with altitude, which we set to be ∝ P^f_sed. Second, the atmospheric mixing coefficient K_zz, which sets the average particle size, once f_sed is fixed. In grid models using the Ackerman & Marley (2001) cloud model, this parameter is usually fixed by mixing length theory (with overshooting) or held constant. Third, the width of the log-normal particle size distribution σ_g, which is normally also kept constant. The cloud mass fraction at the bottom of the cloud was a free parameter, whereas the position of the cloud base was found by intersecting the P-T profile with the saturation vapor pressure curves (taken from Ackerman & Marley 2001, in the corrected pressure units) of the cloud species we considered, Fe and MgSiO₃.

In the future we plan to also test the Burningham et al. (2017) models of clouds, as they are more general, and do not assume the prevalence of a certain cloud species. Moreover, retrieving the power law slope and albedo of the cloud opacities may represent a better choice: for us this is encoded in our choice of cloud species, particle sizes and width of the log-normal particle size distribution, in a non-trivial way. Based on their findings, Burningham et al. (2017) suggest that a log-normal particle size distribution may not be the ideal choice, and that a Hansen distribution (Hansen 1971) may be better.

While carrying out verification retrievals of cloudy petitCODE spectra, we found that we had to be very careful with how the temperaturewas parametrized. If the temperature model was too flexible (e.g., independent layers + p-spline interpolation, as used in Line et al. 2015), test retrievals of cloudy synthetic spectra lead to clear, hot atmospheres with shallow temperature gradients, that well matched the synthetic input spectrum, but were inconsistent with the input temperature and cloud structure. This could indicate that the cloud-free solutions occupied a larger prior volume, and were thus favored when using a Markov chain Monte Carlo (MCMC) retrieval.

Specifically, we found it to be necessary to impose a temperature profile in the photospheric region that follows the Eddington approximation, that is $$T_{\text{phot}}^4=\frac{3}{4}{T}_{\text{0}}^4\left(\frac{2}{3}+\tau \right),$$(3)

where T₀ is normally the internal temperature (taken to be a free nuisance parameter here) and τ the optical depth. This shape was used from τ = 0.1 to the radiative-convective boundary, below which we forced the atmosphere onto a moist adiabat. The optical depth was modeled via $$\tau =\delta P^{\alpha },$$(4)

where δ and α are free parameters. A quite strict prior was imposed on α. We rejected all models where $|\alpha -\tilde{\alpha }|>0.1$, where $\tilde{\alpha }$ is the power law index measured from the opacity structure of a given forward model realization. It was obtained from estimating the Rosseland mean opacity using the non-gray opacity of the atmosphere, across the spectral range of the observations. These altitude-dependent values were then used to calculate an optical depth $\tilde{\tau }$, and from this $$\tilde{\alpha }=<\frac{d\text{log}\tilde{\tau }}{d\text{log}P}>.$$(5)

Here $<>$ denotes the average over the photospheric region. This prior ensures that the parametrized, pressure-dependent opacity is consistent with the atmosphere’s non-gray opacity structure. In future applications of the parametrized P-T we will test to not downright reject models with too large $|\alpha -\tilde{\alpha }|$. Instead one could adapt the log-likelihood by adding $$L_{\alpha }=-\frac{{(\alpha -\tilde{\alpha })}^2}{2{\sigma }_{\alpha }^2}-\frac{1}{2}\text{log}\left(2\pi {\sigma }_{\alpha }^2\right)$$(6)

and fitting for σ_α as a free parameter. Moreover, other P-T parametrizations, for example that of Madhusudhan & Seager (2009), should be tested. This parametrization was also used in Burningham et al. (2017).

In order to prevent the location of the Eddington photosphere to be unrealistically deep in the atmosphere, we also rejected models where $$P(\tau =1)>5P(\tilde{\tau }=1).$$(7)

Above the photosphere the temperature was freely variable. We modeled these high altitudes by retrieving the temperature of three locations spaced equidistantly in log(P) space, and spline interpolating between them.

The chemical abundances and moist adiabat of the atmosphere were found by interpolating in a chemical equilibrium table which contained these quantities as a function of T, P, C/O and [Fe/H]. This table was calculated with the equilibrium chemistry code described in Mollière et al. (2017). In addition, we also retrieved a quench pressure P_quench. At pressures smaller than P_quench the abundances of CH₄, H₂O and CO were held constant, so as to model the effect of chemical quenching in regions where the chemical reaction timescales become longer than the mixing timescales (see, e.g., Zahnle & Marley 2014).

The following absorption opacity sources where included: CO, H₂O, CH₄, NH₃, CO₂, H₂ S, Na, K, PH₃, FeH, VO, TiO, H₂-H₂ (CIA), H₂-He (CIA), Fe clouds (crystalline particles, irregularly shaped), MgSiO₃ clouds (crystalline particles, irregularly shaped). The following scattering opacity sources where included: H₂ Rayleigh scattering, He Rayleigh scattering, Fe clouds, MgSiO₃ clouds. The opacity references can be found in Mollière et al. (2019).

Using the setup described above, we were able to successfully retrieve the spectrum and atmospheric parameters for a synthetic observation of a cloudy, self-consistent model obtained with petitCODE. The implementation of the retrieval forward model presented here will be described in detail in an upcoming paper. It will contain a description of how the scattering was added, and the verification thereof, as well as the verification retrieval test.

The parameter estimation was carried out using emcee (Foreman-Mackey et al. 2013). Due to the complex priors resulting from the temperature parametrization, the high dimensionality, and a potentially multimodal posterior, the acceptance fraction is low (of the order of 1–2%) such that one million models were drawn (started around the best-fit position of a pre-burn) to obtain results where the walker positions had converged. As is common for all parameter estimations using an MCMC method, we cannot guarantee that the retrieval results have converged to the true global maximum of the log-probability. While the multi-modality of our model is an inherent property, the acceptance fraction can be improved by setting up the chain closely around the best-fit position of the pre-burn-in run (Foreman-Mackey et al. 2013). In addition, we are currently working on implementing a parameter estimation using nested sampling, which also applies a clustering algorithm for the parameter estimation. This should alleviate thelow acceptance rate problem, and lead to a complete sampling of the parameter space (Feroz & Hobson 2008; Feroz et al. 2009). With these current limitations in mind, we note that we retrieved similar values as in the grid retrieval with Exo-REM in Sect. 4.2, and could successfully retrieve self-consistent cloudy input models when testing our method.

4.3.2 Retrieval parameter results

Our forwardmodel has 17 free parameters: 6 for the temperature model described above, 3 for the abundances (C/O, [Fe/H], P_quench), 5 for the clouds (the cloud mass fraction of the MgSiO₃ and Fe at the cloud base, the settling parameter f_sed, the eddy diffusion coefficient K_zz, the width of the log-normal particle size distribution σ_g), the gravity log(g), the planetary radius R_Pl, and the abundance of FeH (currently not included in the chemical table). We used uniform or log-uniform priors for all parameters. In addition to the parameters above we allowed for an individual scaling of the GPI (Y, J, H) bands by up to ±50%, and by up to ±2.5% for the GRAVITY data. A large value for the scaling of the GPI data was chosen because a similar scaling value was found when comparing the GPI and SPHERE J-band for 51 Eri b in Samland et al. (2017). The maximum scaling we retrieve for β Pic b is 13% in the GPI Y -band, see below.

Similar to the Exo-REM analysis in Sect. 4.2, we ran retrievals for a GRAVITY only and GRAVITY + GPI case. In Fig. 5 we show the results when fitting the GRAVITY and GPI data together, but without imposing a prior on the mass of β Pic b. For producing this plot, we sampled the posterior distribution 100 times, and plot both the spectra and the accordingly scaled data points. We give the median, and 16 and 84 percentile values of some of the free parameters in the figure. The full posterior and resulting temperature confidence envelopes are shown in Appendix B. The effective temperature was obtained from integrating the flux of the sampled spectra from 0.5 to 20 μm. The mass was calculated from the log(g) and R_P values of the posterior sample.

As can be seen in Fig. 5, the GRAVITY data can be well fit. At least two CO bandheads at ~2.3 μm are visible in the data. The GPI data is less well fit, which may be partially due to the high S/N of the GRAVITY data, which dominates the fit. We found that the fit of the GPI data improved when increasing the error bars of the GRAVITY data. Likewise, we found that the slope in the red part of the GRAVITY spectrum can be fit better if theGPI data are neglected. The retrieved parameters are presented in Table 3, together with the ExoRem results for comparison. Most interestingly, petitRADTRANS retrieves a mass which is consistent with the values from the astrometric measurement in Sect. 3, without the need of imposing a prior on the mass. Here we find $M_{\text{P}}={15.43}_{-2.79}^{+2.91} M_{\text{Jup}}$, which is consistent with the values presented in Table 2.

From the retrievals presented here, it appears that the GRAVITY K-band data are useful for constraining the planetary C/O ratio, while adding the GPI Y, J, and H-bands is important for obtaining better constraints on the planet’s gravity, and hence mass (petitRADTRANS retrieves a too low mass if the GPI data are neglected, see Sect. 4.4 below). This is consistent with the sensitivity of the NIR Y JH bands to gravity (atomic and molecular features, as well as the H band shape), while the CO absorption in the K band is not expected to probe the surface gravity very well (see, e.g., Béjar & Martín 2018). As demonstrated here, estimating the planetary mass from a planet’s spectrum alone may become feasible when applying free retrievals over a broad spectral range, such as carried out here with petitRADTRANS.

4.4 Comparison between the grid and free retrieval

In agreement with the Exo-REM fit for GRAVITY+GPI, petitRADTRANS obtains a cloudy atmosphere, and a slightly sub-solar C/O of ${0.43}_{-0.03}^{+0.04}$ (Exo-REM found 0.43 ± 0.05). The free retrieval obtains a metallicity of ${0.68}_{-0.08}^{+0.11}$. This is higher than Exo-REM, where 0.5 was found, but this was at the boundaries of the Exo-REM grid, and could likely be higher. This could also be the reason for the slightly higher log(g∕g₀) value (${4.34}_{-0.09}^{+0.08}$, and 4 for Exo-REM), due to the gravity-metallicity correlation. petitRADTRANS finds an effective temperature which is higher than in the Exo-REM fit by about 150 K (1742 ± 10 K, compared to 1590 K for Exo-REM). The larger radius found by Exo-REM is likely due to the lower temperature it retrieved, so as to conserve the total amount of flux. At the estimated age of β Pic (24 ± 3 Myr; see Bell et al. 2015), a radius of 1.7 R_Jup (the value Exo-REM retrieved) requires masses in excess of 20 M_Jup, and effective temperatures of around 2500 K, when considering hot start models (Spiegel et al. 2011). Core accretion models under the warm² start assumption, which include deuterium burning, require similarly large masses and temperatures, but potentially somewhat younger ages (Mollière & Mordasini 2012; Mordasini et al. 2017), and would put the planet firmly into the mass regime of brown dwarfs currently undergoing deuterium burning. Hence, the large radius retrieved by Exo-REM is likely to be inconsistent with the retrieved mass and temperature. The values of the mass, temperature, and radius (1.36 R_Jup) retrieved bypetitRADTRANS agree with both the cold and hot start predictions, given the age of β Pic.

Also when fitting only the GRAVITY data, the petitRADTRANS results are mostly consistent with Exo-REM. Without a mass prior we find C/O = ${0.35}_{-0.09}^{+0.07}$ (Exo-REM found ≲0.3), [Fe/H] = $-{0.53}_{-0.34}^{+0.28}$ (Exo-REM found −0.5), $\log (g/g_0)={3.3}_{-0.42}^{+0.54}$ (Exo-REM found 3.5), $M={1.4}_{-0.87}^{+3.94}$ M_Jup (Exo-REM found 2 M_Jup). Only the temperatureis larger again, at 1847 ± 55 K (petitRADTRANS), compared to 1700 K (Exo-REM).

In summary, a free retrieval approach gives similar results to a more classical retrieval from a grid of forward models. We note that here a free retrieval appears to lead to physically more consistent results when constraining radii and effective temperatures. Another possible cause for the differences could be how the opacities of gas and clouds are treated. For the gas opacities we note that petitRADTRANS uses the opacity database of petitCODE, the latter of which was successfully benchmarked with Exo-REM in Baudino et al. (2017). Small remaining differences, identified to stem from the use of different line lists in Baudino et al. (2017), have since been removed by updating the opacity database of petitCODE/petitRADTRANS in 2017.

Fig. 5 Results of the combined (GRAVITY+GPI) fit of the β Pic b spectrum with petitRADTRANS. No prior on the mass was used in the fit, and the spectroscopically retrieved mass is consistent with the astrometric value. For producing this plot, 100 samples were drawn from the posterior distribution, for both the model and the data scaling. The 2-d projection of the posterior can be found in Appendix B. Top panel: GPI Y, J and H-band data of Chilcote et al. (2017) are plotted as green, cyan, and orange points with error bars, respectively, the petitRADTRANS models are plotted as purple solid lines. The fit is dominated by the high S/N of the GRAVITY data (shown in the bottom panel), leading to a worse fit in the GPI bands, see text. Bottom panel: GRAVITY data are shown as black points with errorbars, the petitRADTRANS models are plotted as purple solid lines.

4.5 Comparison to Chilcote et al. (2017)

A substantial analysis of the NIR spectrocopy of β Pic b was carried out in Chilcote et al. (2017), using GPI Y JHK band spectra. The data were compared to low gravity and field brown dwarf spectra, the derived bolometric luminosity was compared to evolutionary models, and spectral fits were carried out with four different model grids.

Comparing their bolometric luminosity to evolutionary models (hot start models of Baraffe et al. 2003), Chilcote et al. (2017) found a mass of 12.9 ± 0.2 M_Jup, an effective temperature of 1724 ± 15 K, a surface gravity of log(g∕g₀) = 4.18 ± 0.01 and a radius of 1.46 ± 0.01 R_Jup. Their mass measurement is consistent both with our astrometrically and spectroscopically inferred mass values. Moreover, the other values inferred from the Y JHK fit of petitRADTRANS are close to the evolutionary values of Chilcote et al. (2017), but not within the uncertainties of one another (e.g., 1742 ± 10 K vs. 1724 ± 15 K). As noted in Chilcote et al. (2017), these uncertainties do not contain a contribution of the model uncertainties, and the true uncertainties must be larger. The same holds for our retrievals and fits carried out here. The Exo-REM fits with mass prior lead to a similar agreement in gravity, but the radii and temperatures are further away from the Chilcote et al. (2017) values, with the planet being cooler, and thus larger, in the Exo-REM fits.

The best grid model fit of the combined photometry and spectroscopy in Chilcote et al. (2017) was obtained with Drift-PHOENIX (e.g., Helling et al. 2008), where T_eff = 1651 K, log (g) = 3 and R = 1.58 R_Jup was found, leading to a mass of ~1 M_Jup. The AMES-Dusty (e.g., Allard et al. 2001) fit gave the highest mass (most consistent with our astrometric, spectroscopic and Chilcote et al. 2017’s evolutionary mass), namely 17 M_Jup. The AMES-Dusty best-fit values are also closer to the evolutionary parameters derived in Chilcote et al. (2017), at T_eff = 1706 K, log (g) = 4.5 and R = 1.18 R_Jup, but at an overall worse χ² than the Drift-PHOENIX fit.

In summary, the results of our spectral characterization compare well with the evolutionary values inferred for β Pic b in Chilcote et al. (2017). Especially the parameter values of the free retrieval carried out with petitRADTRANS are close to the evolutionaryvalues. It is also noteworthy that our masses, inferred indepentently with astrometry or spectral retrieval, are consistent with the evolutionary mass of Chilcote et al. (2017).

5 C/O ratio and the formation of β Pic b

5.1 Stellar and planetary C/O ratio

Holweger et al. (1997) have shown that the abundances of several elements (C, Ca, Ti, Cr, Fe, Sr, Ba) on the surface of β Pictoris are Solar. But measuring the abundance of oxygen in stars is a notoriously difficult task, due to line blending, deviations from local thermal equilibrium predictions, or sensitivity to the 3D temperature structure of the star (Asplund 2005). As a consequence, to our knowledge, the abundance of oxygen – and hence the C/O ratio – has not yet been reported in the literature. We note, though, that a subsolar C/O ratio (i.e., ≃0.4) would invalidate most of the following discussion.

In our atmosphere analysis, both the ExoREM grid model fitting and the petitRADTRANS free retrieval point to the same result: the C/O number ratio in β Pictoris b is ≃0.43 ± 0.05, which is subsolar (the solar C/O ratio is 0.55, see Asplund et al., 2009).

A number of studies have been done in the recent years to try to find links between planetary formation processes and element abundances. In particular, Öberg et al. (2011) first attempted to relate the C/O ratio to the position of the different icelines in a protoplanetary system, and to the proportion of gas and solid material accreted by a young planet. They concluded that substellar C/O ratio was a sign of a formation by either gravitational collapse or core-accretion, followed by icy planetesimal enrichment. The objective of this section is to show that the C/O ratio can possibly be used to disentangle the two formation scenarios.

5.2 General model for the evolution of the C/O ratio

In a similar fashion as to Öberg et al. (2011), we assume that the main sources of carbon and oxygen in the protoplanetary disk in which β Pic b formed were CO, CO₂, H₂O, silicates and carbon grains. Assuming a solar C/O ratio for the star, the table of relative abundances given in Öberg et al. (2011) is valid, and we use it as a baseline to set the abundances of each species (see Table 4).

In the framework developed by Öberg et al. (2011), the C/O ratio in the atmosphere of a planet can be calculated from the amount of solid and gaseous material entering its composition. We denote n_X,s (resp. n_X,g) the abundance of element X in the solid phase (resp. gas phase) of the disk, given in number of atoms per unit of disk mass. We also write M_solid (resp. M_gas) the total mass of solid (resp. gas) entering the composition of the atmosphere of the planet, and f_s∕g the dust-to-gas fraction in the disk, which we assume to be equal to 0.01. With these notations, the total number of elements X in the atmosphere of the planet is given by: $$N_X=\frac{n_{X,s}}{f_{\text{s/g}}}\times M_{\text{solid}}+\frac{n_{X,g}}{1-f_{\text{s/g}}}\times M_{\text{gas}}.$$(8)

And the C/O number ratio is then: $$\text{C/O}=\frac{n_{\text{C},s}{f}_{\text{s/g}}{}^{-1}{M}_{\text{solid}}+n_{\text{C},g}{(1-f_{\text{s/g}})}^{-1}{M}_{\text{gas}}}{n_{\text{O},s}{f}_{\text{s/g}}{}^{-1}{M}_{\text{solid}}+n_{\text{O},g}{(1-f_{\text{s/g}})}^{-1}{M}_{\text{gas}}}.$$(9)

Note that both the numerator and the denominator can be given relative to a reference species without affecting the validity of this Eq. (9). In Table 4 and in all the following, we implictly use abundances relative to H₂O.

The exact values of n_C,s, n_C,g, n_O,s, and n_O,g depends on the abundances given in Table 4, and on the state (solid or gaseous) of each species. and hence on the location of the forming planet with respect to the different icelines.

Using ALMA observations, Qi et al. (2015) have shown that the CO iceline in the disk around HD 163296 was likely to be located at ≃ 90 AU from the star. Other observations of the same system, also performed with ALMA, led Notsu et al. (2019) to conclude that the water iceline was located at a distance of ≤20 AU. Since HD 163296 is also an A-type star, these two values give an idea of the possible location of the H₂O and CO icelines in the β Pic system. However, little is known about the relationship between the current orbit of β Pic b and its exact formation location, and about possible variations of the locations of these icelines between systems. Thus, no definitive assumption can be made as to where the planet formed in comparison to the water iceline, and the two options must be considered: a formation within the water iceline, and a formation between the water and the CO₂ icelines.

From there, the terms n_C,s, n_C,g, n_O,s, and n_O,g from Eq. (9) can be determined from the values listed in Table 4. For a planet forming within the water iceline, we have: $$\{\begin{array}{l}{n}_{\text{O},g}=n_{{\text{H}}_2\text{O}}+n_{\text{CO}}+2\times n_{{\text{CO}}_2}=3.33\hfill \\ n_{\text{O},s}=n_{\text{O} (\text{silicates})}=2.12\hfill \\ n_{\text{C},g}=n_{\text{CO}}+n_{{\text{CO}}_2}=2.0\hfill \\ n_{\text{C},s}=n_{\text{C} (\text{grains})}=\mathrm{0.67.}\hfill \end{array}$$(10)

And for a planet forming between the water and CO₂ icelines: $$\{\begin{array}{l}{n}_{\text{O},g}=n_{\text{CO}}+2\times n_{{\text{CO}}_2}=2.33\hfill \\ n_{\text{O},s}=n_{{\text{H}}_2\text{O}}+n_{\text{O} (\text{silicates})}=3.12\hfill \\ n_{\text{C},g}=n_{\text{CO}}+n_{{\text{CO}}_2}=2.0\hfill \\ n_{\text{C},s}=n_{\text{C} (\text{grains})}=\mathrm{0.67.}\hfill \end{array}$$(11)

Fig. 6 Gravitational collapse scenario: evolution of the C/O ratio as a function of the total mass of solid accreted after the initial formation of the protoplanet. The purple curve corresponds to a formation within the H2O iceline, and the brown curve to a formation between the H2O and CO2 icelines. The orange area gives the 68% confidence interval for the value of the C/O ratio. Dashed vertical lines corresponds to different solid accretion limits discussed in the text.

5.3 C/O ratio in the gravitational collapse paradigm

We consider the case of a formation through gravitational collapse (Bodenheimer 1974), a violent mechanism which shares similarities with star formation. In this scenario, an entire region of the circumstellar disk becomes unstable, and rapidly collapses to form a protoplanet, which then slowly contracts and cools down.

The total mass of solid entering in the composition of the atmosphere of a planet formed through gravitational collapse can be separated in two terms: the mass of solid initially contained in the disk fragment which collapsed to create the protoplanet, and the mass of solid planetesimals later accreted by the protoplanet. The solid mass contained in the initial clump is directly related to the dust-to-gas ratio of the disk, and we can write: $$M_{\text{solid}}=f_{\text{s/g}}{M}_{\text{planet}}+M_{\text{accreted}}.$$(12)

This equation assumes that no core has formed in the young protoplanet, which, for a planet as massive as β Pic b is reasonable (Helled & Schubert 2008). For a planet less massive, for which a core could form, sedimentation of a fraction of the initial solid mass on the core should be taken into account.

Injecting the definition of M_solid into Eq. (9), and using f_s∕g = 0.01, M_planet = 12.7 M_Jup, as well as the values for C and O abundances given in Eq. (10) or (11), it is possible to determine the C/O ratio as a function of the mass of accreted planetesimals M_accreted in the gravitational collapse paradigm. The results are given in Fig. 6, for two possible formation locations: within the water iceline, and between the water and carbon dioxide icelines. We have also added the 1 σ confidence intervals of our ExoREM and petitRADTRANS measurements on this graph. This figure shows that a formation bewteen the H₂O and CO₂ icelines is more favorable to a large deviation from the stellar C/O ratio, mainly due to the injection of oxygen coming from solid water ice during planetesimal accretion.

The formation of a planet by gravitational instability can be separated in a few different steps (Bodenheimer 1974): formation of the initial clump in the disk, quasi-equilibrium contraction, hydrodynamic collapse, and a new hydrostatic quasi-equilibrium phase. Accretion of planetesimals is thought to be efficient only during the pre-collapse phase (Helled & Schubert 2009). The duration of this phase decreases with increasing planet mass, and typical values ranges from a few 10⁵ yr for a Jupiter mass planet, to less than 10³ yr for more massive planets (Decampli & Cameron 1979; Bodenheimer et al. 1980). Using the model proposed by Helled & Schubert (2009), the mass of planetesimal accreted during the pre-collapse phase of β Pic b can be estimated using: $$M_{\text{accreted}}={\displaystyle {\int }_0^{t_{\text{collapse}}}\pi }{R}_{\text{capture}}^2(t)\sigma (a,t)\Omega (a)\text{d}t,$$(13)

where t_collapse is the time of collapse, R_capture the protoplanet’s capture radius, σ the surface density of solids in the disk at the location of the protoplanet, and Ω the orbital frequency.

Andrews & Williams (2005) presented a large survey of 153 young stellar objects in the Taurus-Auriga star forming region. Among all these objects, AB Aur and V892 Tau are two A-type stars, for which they give an estimate of the mass: 0.004 M_⊙ and 0.009 M_⊙. Considering all stellar types, the median disk-to-star mass ratio they found is 0.5%. More recent studies of protoplanetray disk demographics based on ALMA observations yielded similar results, with typical dust to star mass ratios of ≃10^−4.5 (Pascucci et al. 2016; Ansdell et al. 2017), i.e. disk-to-star mass ratios of ≃0.3%, assuming a dust-to-gas ratio of 1%.

Considering the upper limit of an extremely massive disk (M_disk = 0.1 M_⊙), and using a power-law for the surface density ($\sigma ={\sigma }_0{(a/5 \text{AU})}^{-\alpha }$, with α = 1.00), the solid density at a = 11 AU is: $$\sigma (11 \text{AU})\simeq 6 \text{g\hspace{0.17em}}{\text{cm}}^{-2}.$$(14)

The orbital period of the planet is ~20 yr (Wang et al. 2016; Lagrange et al. 2018, Sect. 3 of this work). The capture radius decreases with the contraction of the planet, but an optimistic value would be 2 to   3 × 10¹² cm for a 1 M_Jup planet (Helled et al. 2006). For a planet 10 times more massive, the effective radius could be ≃ 5 × 10¹² cm. This yields: $$M_{\text{accreted}}\simeq 4\times M_{\text{Earth}}\times \frac{t_{\text{collapse}}}{1000 \text{yr}}.$$(15)

The corresponding accretion limit has been added to Fig. 6, for a reasonable assumption of t_collapse = 10³ yr.

Taking into account the effective time available for efficient planetesimal accretion during the pre-collapse stage, the low C/O ratio measured with GRAVITY is difficult to explaine, even in the case of a planet forming oustide the H₂O iceline. For the C/O ratio to reach a value of ≃0.43, we need to assume a massive protoplanetary disk and an unusually long time for the pre-collapse phase, or an extremely efficient accretion (with an accretion rate of 4 × 10⁻³ M_Earth yr⁻¹).

Fig. 7 Core-accretion scenario: evolution of the C/O ratio in the atmosphere of β Pic b as a function of the total mass of solid accreted by the protoplanet, for a formation between the H2O and the CO2 icelines, or within the H2O iceline. The orange area gives the 68% confidence interval for the value of the C/O ratio.

5.4 C/O ratio in the core-accretion paradigm

Core-accretion is another formation mechanism, in which an initial solid core forms, and slowly accretes gas from the disk. When the mass of gas is roughly the same as the mass of the core, the protoplanet enters a phase of “runaway gas accretion”, during which it gains significant amount of gas over a short time (Lissauer & Stevenson 2007). In this scenario, the formation of a planet is a much longer process than with gravitational instability, which gives more time to enrich the proto-atmosphere in solid material and to lower its C/O ratio.

Mordasini et al. (2016) explored the effect of planetesimal enrichment coupled with disk composition, in a core-accretion scenario. They focused on the case of Jupiter mass planets migrating to short period orbits (“hot Jupiters”), which is a different archetype than β Pic b. But the general sequence of events they use to form their planets in the core-accretion paradigm can still be applied to β Pic b, only leaving out the inward migration part. First, the core of the planet forms from the accretion of solid material. Then, once the core has formed, the protoplanet starts accreting a gaseous envelope which, during its formation, is enriched by the accretion of disintegrating planetesimals. When the planet reaches a critical mass, runaway accretion occurs, and the mass of the planet significantly increases. This runaway gas accretion clears a gap in the disk, and ends the formation of the planet.

In the gravitational instability scenario, because the formation of the planet happens so quickly compared to typical timescales of disk evolution, the gas and solid making the atmosphere necessarily have a stellar combined C/O. If the solid and gas in the atmosphere are in the same proportion as they are in the disk (M_solid = f_s∕gM_gas), the C/O of the atmosphere is stellar. A deviation of the solid to gas proportion in the atmosphere is required to alter the C/O ratio.

In the case of core-accretion, the situation is different. Without planetismal enrichment before the runaway gas accretion phase, the atmosphere of the planet would not be made of a mixture of gas and solid material, but purely of gas. Thus, without planetesimal enrichment, the atmospheric C/O ratio in the core-accretion paradigm can be expected to be close to the C/O ratio of the gas in the disk, that is, superstellar.

In this core-accretion paradigm, it is still possible to use Eq. (9) to calculate how the final C/O ratio of the atmosphere is impacted by the mass of solid material accreted before the runaway accretion phase. But in this case, all of the solid mass M_solid corresponds to accreted material: M_solid = M_accreted, as opposed to Eq. (12).

In Fig. 7, we show the evolution of the C/O ratio as a function of the mass of accreted planetesimals, for a formation thourgh core-accretion, within the water iceline, or between the water and CO₂ icelines. In this scenario, it is possible to reach C/O values compatible with our GRAVITY measurement with accretion of ≃ 80 M_Earth, if the planet formed between the water and CO₂ icelines. A formation within the water iceline is more diffcult to explain, as it would require at least 150 M_Earth of solid material enrichment to reach the upper limit of the 1 σ interval on the C/O measurement, and up to several 10² M_Earth to reach a value of 0.43.

6 Summary and conclusions

In this work, we presented the first VLTI/GRAVITY spectro-interferometric observation of the giant planet β Pictoris b. Using an adequate data reduction technique detailed in the appendix of this paper, we extracted a high quality K-band spectrum of the planet, at a resolution of R = 500. We also derived the most precise relative astrometry obtained to date on this object, with an error of ≃ 40 μas.

We find that the astrometry disfavors circular orbits for β Pic b, with a value of $e\simeq {0.15}_{-0.04}^{+0.05}$. It remains unclear how a massive planet like β Pic b can acquire such a significant eccentricity. Using this new astrometric datapoint together with previously published visual astrometry and HIPPARCOS/Gaia data, we were able to derive an estimate of the dynamical mass of β Pictoris b, in a similar fashion as to what Snellen & Brown (2018) and Dupuy et al. (2019) did. Our value is compatible with these previous studies, with a best estimate of 12.7 ± 2.2 M_Jup.

We were also able to retrieve a similar mass, albeit with larger error bars, using only the spectral data. Using a free retrieval, including the effect of scattering and clouds, with petitRADTRANS (Mollière & Snellen 2019) to fit the spectrum of β Pic b in Y, J, H, and K bands (Y, J, H from Chilcote et al. 2017, K from this work), we obtained a mass of ${15.43}_{-2.79}^{+2.91}{M}_{\text{Jup}}$. This constitutes a rare case of validation of an atmospheric model with a model-independent measurement.

We performed an in-depth analysis of the K-band spectrum extracted from our GRAVITY observation using two different approaches: forward modeling with the ExoREM code (Charnay et al. 2018), and free-retrieval with petitRADTRANS. We found that both approaches point to a C/O ratio of C∕O = 0.43 ± 0.05.

We showed that, if the C/O ratio of the host star β Pictoris is Solar, it is difficult to explain this C/O ratio with a gravitational collapse formation scenario. This is mainly due to the high mass of β Pictoris b, which hasthe dual consequence of requiring large amount of planetesimal enrichment to lower the initial C/O ratio, while at the same time makingthe whole formation process extremely short. In this case, it appears that a slower formation via core-accretion, somewhere between the H₂O and CO₂ icelines, is more likely. This scenario can potentially explain the subsolar C/O ratio if the planet was enriched in oxygen by icy planetesimal accretion.

The high metal enrichment we retrieve from the spectral fits appears to corroborate this assessment, with the exact value being quite high and at the edge of what is expected from classical core accretion Mordasini et al. (2016).

This model still comes with several important limitations. One of them is that the exact compositition of the initial protoplanetary disk around β Pic remains largely unknown. Another major issue is the efficiency of the planetesimal enrichment, which we have assumed to be of 100% (i.e., all the solid material accreted by the planet is disintegrated in the atmosphere). This is unlikely to be the case, as fraction of this material can be deposited into the planetary core, or can stay at the bottom of the atmosphere. This is particularly true for the core-accretion scenario, in which the solid material is accreted before most of the gas (Mordasini et al. 2016). Strong vertical mixing can potentially mitigate this problem, but further studies are required to be able to take into account these phenomena. Finally, disk chemistry may also play a role. For example, Eistrup et al. (2018) have shown that a large fraction of water molecules can be transformed into dioxygen (O₂) over a few Myr, along a chemical pathway detailed in Walsh et al. (2015). Oustide of the water iceline, such a chemical evolution can potentially deplete the solid material from its oxygen, while enriching the gas.

The observations of β Pictoris b presented in this paper show the potential of long-baseline optical interferometry with VLTI/GRAVITY for exoplanet science. The instrument gives access to medium resolution spectroscopy in K-band and high-precision astrometry, which are both extremely useful to characterise giant exoplanets and to start peering into their formation history.

Acknowledgements

Based on observations collected at the European Southern Observatory under ESO programme 0101.C-0912(A) and 2101.C-5050(A). M.N. acknowledges funding for his PhD from the European Research Council (ERC), under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 639248). P.M. thanks M. Line for insightful discussions. P.M acknowledges support from the ERC under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 694513 and 832428). J.W. thanks R. De Rosa for helpful discussions on Gaia and HIPPARCOS data. J.W. is supported by the Heising-Simons Foundation 51 Pegasi b postdoctoral fellowship. R.G.L. received financial support of Science Foundation Ireland (Grant number 18/SIRG/5597). A-M.L. acknowledges support from the French CNRS and from the Agence Nationale de la Recherche (ANR grant GIPSE ANR-14-CE33-0018). T.H. acknowledges support from the ERC under the Horizon 2020 Framework Program (ERC Advanced Grant Origins 832428).

Appendix A Reduction of the GRAVITY dataset

A.1 Nomenclature and pipeline errors

The data reduction used to extract the β Pictoris b signal from the GRAVITY observations makes heavy use of complex linear algebra, complex error formalism, and maximum likelihood estimation. To avoid confusion and mistakes, complex numbers in this appendix are underlined (e.g., $\underset{\_}{V}$), whereas real numbers arenot (e.g., X).

Most of the GRAVITY data manipulated are quantities which depends on the wavelength λ. These quantities can be represented as vectors of size n_λ (the number of wavelength channels) by concatenating the individual values. These vectors are denoted using a bold font. For example, in the case of the complex visibility obtained on baseline b at time t, we denote: $${\underset{\_}{V}}_{b,t}=\left(\begin{array}{c}\underset{\_}{V}(b,t,{\lambda }_1)\\ \underset{\_}{V}(b,t,{\lambda }_2)\\ ⋮\\ \underset{\_}{V}(b,t,{\lambda }_{n_{\lambda }}).\end{array}\right)$$(A.1)

For a given DIT, it is also possible to concatenate all baselines to create a vector of size n_b × n_λ, where n_b = 6 is the number of baselines. In this case, the subscript b is dropped: $${\underset{\_}{V}}_t=\left(\begin{array}{c}{\underset{\_}{V}}_{b_1,t}\\ ⋮\\ {\underset{\_}{V}}_{b_{n_b},t}.\end{array}\right)$$(A.2)

The complex-conjugate of a complex number $\underset{\_}{V}$ is denoted ${\underset{\_}{V}}^{*}$, and the complex-transpose of a vector or matrix $\underset{\_}{A}$ is denoted ${\underset{\_}{A}}^{†}$. It is defined by: ${\underset{\_}{A}}^{†}={\underset{\_}{A}}^{*}{}^T$ where ^T is the transpose operator.

All the λ-vectors are understood as elements of a n_λ -dimension complex linear space (i.e. a linear space for which the scalar field is the set of complex numbers $ℂ$, rather than the set of real numbers $ℝ$). Adding the natural scalar product operator (i.e. $\langle V_1,V_2\rangle =V_1{}^{†}{V}_2$) makes this linear space an Euclidean space. This mathematical structure allows for several useful concepts: it is possible to compute othogonal projections, to use projector matrices, to define othogonal and/or orthonormal basis, etc.

The data set can be subdivided into two parts: the observations taken with the science fiber on the planet, and the observationstaken on the star (see observing log in Table 1). On-planet and on-star phase-referenced visibilities are calculated from the coherent fluxes measured by GRAVITY, called VISDATA in the FITS files generated by the pipeline. The VISDATA are complex numbers, affected by noise. The GRAVITY pipeline reports these errors in another set of complex numbers, called VISERR. The real part of VISERR contains the uncertainties on the real part of VISDATA, and the imaginary part of VISERR contains the uncertainty on the imaginary part of VISDATA. These errors do not take into account any possible correlation between different spectral channels, or between the real and imaginary parts of the visibility. To take into account such correlations, it is necessary to use the covariance/pseudo-covariance formalism of complex random variables.

In our data reduction algorithm, the GRAVITY pipeline errors are systematically replaced by an empirical estimate of the covariance and pseudo-covariance matrices of the visibilities. We assume that the noise affecting the measurements does not vary significantly over the individual DITs of a single exposure file (~ 5 min), but can vary from file to file. We also allow for correlations between different spectral channels and/or between different baselines. Under these assumptions, the errors on the coherent fluxes are best represented by a set of n_EXP (the number of exposure files) covariance matrices W_k and n_EXP pseudo-covariance matrices Z_k, both of size n_b × n_λ, where n_b = 6 is the number of baselines and n_λ = 235 is the number of wavelength channels. The covariance and pseudo-covariance matrices for each exposure file are estimated directly from the DITs sequence: $$\begin{array}{l}{\underset{\_}{W}}_k=\frac{1}{n_{\text{DIT}}-1}\left({\sum }_{t=1}^{n_{\text{DIT}}}{\underset{\_}{V}}_t{\underset{\_}{V}}_t{}^{†}-\frac{1}{n_{\text{DIT}}}\left({\sum }_{t=1}^{n_{\text{DIT}}}{\underset{\_}{V}}_t\right){\left({\sum }_{t=1}^{n_{\text{DIT}}}{\underset{\_}{V}}_t\right)}^{†}\right)\\ {\underset{\_}{Z}}_k=\frac{1}{n_{\text{DIT}}-1}\left({\sum }_{t=1}^{n_{\text{DIT}}}{\underset{\_}{V}}_t{\underset{\_}{V}}_t{}^T-\frac{1}{n_{\text{DIT}}}\left({\sum }_{t=1}^{n_{\text{DIT}}}{\underset{\_}{V}}_t\right){\left({\sum }_{t=1}^{n_{\text{DIT}}}{\underset{\_}{V}}_t\right)}^T\right)\end{array}$$

where the dummy t runs over the n_DIT DITs of the kth exposure.

The covariance and pseudo-covariance matrices are always related to the covariance of the real and imaginary parts by the following equations: $$\begin{array}{l}\text{cov}\left(\text{Re}\left(\underset{\_}{V}\right),\text{Re}\left(\underset{\_}{V}\right)\right)=\frac{1}{2}\text{\hspace{0.17em}Re}\left(\underset{\_}{W}+\underset{\_}{Z}\right)\\ \text{cov}\left(\text{Im}\left(\underset{\_}{V}\right),\text{Im}\left(\underset{\_}{V}\right)\right)=\frac{1}{2}\text{\hspace{0.17em}Re}\left(\underset{\_}{W}-\underset{\_}{Z}\right)\\ \text{cov}\left(\text{Re}\left(\underset{\_}{V}\right),\text{Im}\left(\underset{\_}{V}\right)\right)=\frac{1}{2}\text{\hspace{0.17em}Im}\left(-\underset{\_}{W}+\underset{\_}{Z}\right)\\ \text{cov}\left(\text{Im}\left(\underset{\_}{V}\right),\text{Re}\left(\underset{\_}{V}\right)\right)=\frac{1}{2}\text{\hspace{0.17em}Im}\left(\underset{\_}{W}+\underset{\_}{Z}\right).\end{array}$$

The covariance and pseudo-covariance matrices can be propagated during the data reduction algorithm by using the complex error propagation equations: $$\begin{array}{l}\text{cov}\left(\underset{\_}{A}\underset{\_}{V}\right)=\underset{\_}{A}\underset{\_}{W}{\underset{\_}{A}}^{†}\\ \text{pcov}\left(\underset{\_}{A}\underset{\_}{V}\right)=\underset{\_}{A}\underset{\_}{Z}{\underset{\_}{A}}^T\end{array}$$

with $\underset{\_}{A}$ any complex matrix of appropriate size.

The $\underset{\_}{W}$ and $\underset{\_}{Z}$ matrices can also be used to resolve linear equations involving complex data. In the case of an unknown real parameter vector X, the solution of the linear problem $\underset{\_}{V}=\underset{\_}{A}X$ (in the sense of maximum likelihood) is: $$\widehat{X}={\left(\text{Re}\left({\underset{\_}{V}}_2{}^{†}{\underset{\_}{W}}_2^{-1}{\underset{\_}{A}}_2\right){\left[\text{Re}\left({\underset{\_}{A}}_2{}^{†}{\underset{\_}{W}}_2^{-1}{\underset{\_}{A}}_2\right)\right]}^{-1}\right)}^T$$(A.9)

where $$\begin{array}{lll}\hfill & \hfill & {\underset{\_}{V}}_2=\left(\begin{array}{c}\underset{\_}{V}\\ {\underset{\_}{V}}^{*}\end{array}\right)\hfill \\ \hfill & \hfill & {\underset{\_}{W}}_2=\left(\begin{array}{cc}\underset{\_}{W}& \underset{\_}{Z}\\ {\underset{\_}{Z}}^{°}& {\underset{\_}{W}}^{*}\end{array}\right)\hfill \\ \hfill & \hfill & {\underset{\_}{A}}_2=\left(\begin{array}{c}\underset{\_}{A}\\ {\underset{\_}{A}}^{*}\end{array}\right).\hfill \end{array}$$

For a complete mathematical derivation of Eq. (A.9), we refer the reader to Appendix B of Nowak (2019).

A.2 Pipeline reduction and phase referencing

The initial step uses the pipeline reduction and is common to all VLTI/GRAVITY observations. It consists in extracting the complex visibilities from the raw data, using the ESO pipeline (Lapeyrere et al. 2014). The pipeline takes care of the background subtraction, flat-field correction, bad-pixel interpolation and P2VM multiplication (Tatulli et al. 2007). It also corrects the phase of the visibilities using the metrology data, and combines all DITs within each exposure. This last step performed by the pipeline (averaging of all DITs within each exposure file) is unwanted for exoplanet observations (see Appendix A.1). Thus, for the β Pic b observations, an intermediate file product generated by the pipeline is used: the “astrored” files, in which all DITs are kept separate. The complex visibilities contained in the “astrored” files are not corrected for the metrology and fringe-tracker zero-point, and thus the correction must be applied manually (see recipe in Nowak, 2019).

For each baseline b, and each time t (i.e. for each DIT), the wavelength-dependent complex visibility VISDATA_onstar and VISDATA_onplanet extracted by the pipeline (and with the above-mentionned corrections) are then “phase-referenced” to the star: $$\begin{array}{l}{\underset{\_}{V}}_{\text{onstar}}=\left|{\text{VISDATA}}_{\text{onstar}}\right|\\ {\underset{\_}{V}}_{\text{onplanet}}={\text{VISDATA}}_{\text{onplanet}}\times {\text{e}}^{-i\arg ({\text{VISDATA}}_{\text{onstar}})}\end{array}$$

where arg (VISDATA_onstar) is the phase of the stellar complex visibility as measured by GRAVITY when the science fiber is positioned on the star.

This phase-referencing step is performed both on the star exposures, and on the planet exposures. When dealing with a star exposure, phase-referencing the visibility is mathematically equivalent to extracting the modulus of the visibility. But when dealing with on-planet exposures, a problem arises: the instrument does not simultaneously observe both the planet and the star. Thus, the quantity arg (VISDATA_onstar) must be estimated from the available star exposures. In the observing strategy used for acquiring the β Pictoris data reported here, a star exposure was performed before and after each on-planet exposure. For each on-planet exposure, the phase reference is then simply estimated by taking the phase of the stellar complex visibily averaged on these two star exposures (before and after the on-planet observation).

A.3 A model for the on-planet visibility

In the absence of stellar flux, the on-planet visibility measured by the instrument and phase-referenced to the star can be written: $${\underset{\_}{V}}_{\text{onplanet}}(b,t,\lambda )=\underset{\_}{G}(b,t,\lambda ){\underset{\_}{V}}_{\text{planet}}(b,t,\lambda )$$(A.15)

in which V_planet is the planet astrophysical visibility phase-referenced to the star, and G is the instrumental response. We note that the visibilities are not calibrated, meaning that the visibility at zero frequency is not 1, but the un-normalized flux. Therefore, as long as the planet remains unresolved by the instrument, its astrophysical visibility is given by: $${\underset{\_}{V}}_{\text{planet}}(b,t,\lambda )=S_{\text{planet}}(\lambda )\times {\text{e}}^{-i\frac{2\pi }{\lambda }\left(\Delta \alpha \times U+\Delta \delta \times V\right)}$$(A.16)

in which (U, V) are the coordinates of baseline b in the UV-plane, (Δα, Δδ) the sky-coordinates of the planet relative to the star, and S_planet(λ) the spectrum of the planet.

Given the typical VLTI baseline lengths (between 45 and 130 m with the UTs), and the expected β Pic b planet-to-star separation at time of observation (≃140 mas), the exponential term in the above equation should produce significant oscillations of the complex visibility over the GRAVITY wavelength range (1.9–2.35 μm). But the phase-referenced on-planet visibilities extracted from our β Pic b observations show no such oscillations. The reason is that the data are dominated by remaining starlight, which needs to be taken into account.

To take into account the coherent starlight leaking into the fiber, Eq. (A.15) must be modified with an additional term, proportional to the stellar phase-referenced visibility V_star. In practice, since the leaking starlight does not originate in the direct coupling of the star to the fiber, but rather in the coupling of speckle noise to the fiber, this term needs to be multiplied by a polynomial in λ to account for its chromaticity. The model is now given by: $${\underset{\_}{V}}_{\text{onplanet}}(b,t,\lambda )=\text{\&} \underset{\_}{Q}(b,t,\lambda )\underset{\_}{G}(b,t,\lambda )V_{\text{star}}(b,t,\lambda )\text{nn\&}+\underset{\_}{G}(b,t,\lambda ){\underset{\_}{V}}_{\text{planet}}(b,t,\lambda )$$(A.17)

with $\lambda \to \underset{\_}{Q}(b,t,\lambda )$ a polynomial function in λ, whose coefficients vary with baseline b and time t.

If the star is not resolved by the instrument, its astrophysical phase-referenced visibility corresponds to its spectrum. If the star is partially resolved by the instrument, the spectrum needs to be multiplied by a term accounting for the resulting drop in visibility, which depends on the angular size of the star, limb-darkening model, etc. Explicitly separating these two terms, the referenced stellar astrophysical visibility can be written using the following equation, in which S_⋆ (λ) is the star spectrum, and J a function accounting for the visibility drop due to the star geometry (typically, J is a bessel function of first order): $$V_{\text{star}}(b,t,\lambda )=S_{\star }(\lambda )J(b,t,\lambda ).$$(A.18)

Going back to Eq. (A.17), the planet term in the right-hand side can be factored by V_⋆ by introducing the planet-to-star contrast spectrum C(λ) = S_planet(λ)∕S_star(λ): $${\underset{\_}{V}}_{\text{onplanet}}=\underset{\_}{Q}\underset{\_}{G}{V}_{\text{star}}+J^{-1}\underset{\_}{G}{V}_{\text{star}}C(\lambda ){\text{e}}^{-i\frac{2\pi }{\lambda }\left(\Delta \alpha U+\Delta \delta V\right)}.$$(A.19)

The on-star equivalent of Eq. (A.17) is simpler, as the reference visibility observed on-star only depends on the stellar referenced visibility and the instrumental response: $${\underset{\_}{V}}_{\text{onstar}}(b,t,\lambda )=\underset{\_}{G}(b,t,\lambda )V_{\text{star}}(b,t,\lambda ).$$(A.20)

This provides a natural way to estimate the term $\underset{\_}{G}{V}_{\text{star}}$ in Eq. (A.19), and thus to calibrate ${\underset{\_}{V}}_{\text{onplanet}}$: $$\underset{\_}{U}(b,t,\lambda )=\frac{{\underset{\_}{V}}_{\text{onplanet}}}{{\underset{\_}{V}}_{\text{onstar}}}(b,t,\lambda )=\underset{\_}{Q}+J^{-1}C(\lambda ){\text{e}}^{-i\frac{2\pi }{\lambda }\left(\Delta \alpha U+\Delta \delta V\right)}.$$(A.21)

Equation (A.21) shows how the physical quantities of interest (i.e. the contrast spectrum C(λ) and the planet separation Δα, Δδ) are encoded in the on-planet data. Even taking into account the filtering of the starlight by the off-axis fiber, as well as the only partly coherent nature of the speckle noise, the polynomial $\underset{\_}{Q}$ in the right-hand side of Eq. (A.21), which model the stellar residuals, is still a factor 20 to 30 superior to the planet signal. It is only because of the phase modulation naturally introduced by the planet separation that this planet signal can be retrieved. To do so, we proceed in two steps: we first extract the star-planet separation vector under some hypothesis on the contrast spectrum, and we then extract the contrast spectrum using the estimated separation vector. The two steps are iterated on time, to check the consistency of the results.

In matrix notations, the multiplications by $\underset{\_}{G}$, J⁻¹, the exponential, or even the polynomial $\underset{\_}{Q}$ can all be represented by diagonal-matrix multiplications. We write: $${\underset{\_}{U}}_{b,t}={\sum }_{k=0}^m{\underset{\_}{a}}_{b,t,k}{\Lambda }^{k}1+J_{b,t}^{-1}{\Phi }_{b,t}^{\Delta \alpha ,\Delta \delta }C,$$(A.22)

where m is the order of the polynomial $\underset{\_}{Q}$, and the ${\underset{\_}{a}}_k$s are complex coefficients used to describe the polynomial. The vector C is defined from C(λ) using the notations introduced in Sect. A.1, 1 is a column vector filled with 1’s, and the matrices Λ, J⁻¹, and $\underset{\_}{\Phi }$ are all diagonal matrices of size n_λ × n_λ defined by: $$\underset{\_}{\Lambda }=\text{diag}\left\{{\lambda }_1,\dots ,{\lambda }_{n_{\lambda }}\right\}{\underset{\_}{J}}_{b,t}^{-1}=\text{diag}\left\{J{(b,t,{\lambda }_1)}^{-1},\dots ,J{(b,t,{\lambda }_{n_{\lambda }})}^{-1}\right\}{\underset{\_}{\Phi }}_{b,t}^{\Delta \alpha ,\Delta \delta }=\text{diag}\left\{e^{-i\frac{2\pi }{{\lambda }_1}\left(\Delta \alpha U(b,t)+\Delta \delta V(b,t)\right)},\dots \right\}.$$(A.23)

A.4 Extracting the astrometry

At the initial iteration, the planet to star contrast spectrum C in Eq. (A.22) can most generally be replaced by a flat spectrum. In the case of β Pictoris b, the temperature and surface gravity of the planet are known from previous work (Chilcote et al. 2017). Thus, the contrast spectrum C can be setto a model value. We use a BT-Settl model (Baraffe et al. 2015), at T = 1700 K and log (g∕g₀) = 4.0 (planet spectrum), divided by a BT-NextGen model (Hauschildt et al. 1999) at T = 8000 K and log (g∕g₀) = 4.0 (the star).

The contrast spectrum being set to a pre-determined value, Eq. (A.22) becomes a model at (m + 1) × n_bn_DIT complex parameters (the ${\underset{\_}{a}}_k$s), and 2 real parameters (Δα and Δδ).

Interestingly, this model is linear in all the ${\underset{\_}{a}}_k$, and nonlinear in Δα, Δδ. To fully benefit from this for the model-fitting, Eq. (A.22) can be re-arranged in a pseudo matrix form, with real parameters.

Introducing: $$x_{b,t}=\left(\begin{array}{c}\text{Re\hspace{0.17em}}({\underset{\_}{a}}_0)\\ \text{Im\hspace{0.17em}}({\underset{\_}{a}}_0)\\ ⋮\\ \text{Re\hspace{0.17em}}({\underset{\_}{a}}_m)\\ \text{Im\hspace{0.17em}}({\underset{\_}{a}}_m).\end{array}\right)$$(A.24)

And ${\underset{\_}{A}}_{b,t}^{\Delta \alpha ,\Delta \delta }$ defined column by column: $${\underset{\_}{A}}_{b,t}^{\Delta \alpha ,\Delta \delta }=\left({\Lambda }^{0}1,i{\Lambda }^{0}1,\dots ,{\Lambda }^{m}1,i{\Lambda }^{m}1,{\underset{\_}{\Phi }}_{b,t}^{\Delta \alpha ,\Delta \delta }C\right).$$(A.25)

We have: $${\underset{\_}{U}}_{b,t}={\underset{\_}{A}}_{b,t}^{\Delta \alpha ,\Delta \delta }{x}_{b,t}.$$(A.26)

For a given Δα and Δδ, the corresponding best estimate of x_b,t (in the sense of the maximum likelihood) is given by: $${\widehat{x}}_{b,t}={\left(\text{Re}\left({\underset{\_}{U}}_{b,t,2}{}^{†}{\underset{\_}{W}}_2^{-1}{\underset{\_}{A}}_{b,t,2}\right){\left[\text{Re}\left({\underset{\_}{A}}_{b,t,2}{}^{†}{\underset{\_}{W}}_2^{-1}{\underset{\_}{A}}_{b,t,2}\right)\right]}^{-1}\right)}^T$$(A.27)

where ${\underset{\_}{W}}_2$ is a matrix composed of the covariance and pseudo covariance matrices of ${\underset{\_}{U}}_{b,t}$ as defined in Eq. (A.11).

The log-likelihood $\log {\mathcal{L}}_{b,t}(\Delta \alpha ,\Delta \beta )$, restricted to baseline b, DIT t, and to the nonlinear parameter Δα and Δδ is then given by the following equation, in which the dependance in Δα, Δδ of the right-hand side is implicit in the definition of $\underset{\_}{A}$ and $\widehat{x}$. $$-\log {\mathcal{L}}_{b,t}={\left[{\underset{\_}{U}}_{b,t}-{\underset{\_}{A}}_{b,t}{\widehat{x}}_{b,t}\right]}_2{}^{†}{\underset{\_}{W}}_{b,t,2}^{-1}{\left[{\underset{\_}{U}}_{b,t}-{\underset{\_}{A}}_{b,t}{\widehat{x}}_{b,t}\right]}_2.$$(A.28)

The total log-likelihood can be obtained by summing over all baselines b and DITs t: $$-\log \mathcal{L}(\Delta \alpha ,\Delta \delta )=-{\sum }_{b,t}\log {\mathcal{L}}_{b,t}(\Delta \alpha ,\Delta \delta ).$$(A.29)

With the expression of ${\underset{\_}{A}}_{b,t}$ from Eq. (A.25) and $\widehat{x}$ from Eq. (A.27), this gives a closed-form expression of the log-likelihood in Δα, Δδ, from which a map can be calculated, in order to extract the best estimate with the associated error bars.

A.5 Extracting the spectrum

Extracting the spectrum from the GRAVITY observations is more difficult than extracting the astrometry for two reasons: first, due to the dimensionality of the problem (n_λ > 200), a $\log \mathcal{L}$ map approach is impractical; second, the stellar residuals affecting the on-planet visibility can lead to a degenerated solution for the contrast spectrum C.

The impact of the stellar residuals on the contrast spectrum can be quantified by using Eq. (A.22) again. The calculations in the rest of this section are simpler when considering all visibilities shifted to the planet position. To do so, we multiply both sides of Eq. (A.22) by the inverse of ${\underset{\_}{\Phi }}_{b,t}^{\Delta \alpha ,\Delta \delta }$. From now on, we will denote ${\widetilde{\underset{\_}{V}}}_{b,t}$ any “shifted” visibility ${\underset{\_}{V}}_{b,t}$. Equation (A.22) becomes: $${\widetilde{\underset{\_}{U}}}_{b,t}={\sum }_{k=0}^m{\underset{\_}{a}}_{b,t,k}{\Lambda }^k\underset{\_}{\tilde{1}}+J_{b,t}^{-1}C.$$(A.30)

From there, we can introduce the subspace ${ℂ}^m[\Lambda ]\underset{\_}{\tilde{1}}$ of ${ℂ}^{n_{\lambda }}$, defined as the subspace generated by the family of m + 1 vectors ${\Lambda }^0\underset{\_}{\tilde{1}},\dots ,{\Lambda }^m\underset{\_}{\tilde{1}}$ (i.e., the subspace of vectors which are linear combinations of these m + 1 vectors). Introducing this subspace if of course motivated by the fact that the stellar residual term is part of it. We can then introduce the projector matrix orthogonal to this subspace, which we denote³ ${\underset{\_}{P}}_{{ℂ}^m[\Lambda ]{\widetilde{\underset{\_}{1}}}_{b,t}}$. ProjectingEq. (A.22) then gives: $${\underset{\_}{P}}_{{ℂ}^m[\Lambda ]{\underset{\_}{\tilde{1}}}_{b,t}}{\underset{\_}{\tilde{U}}}_{b,t}={\underset{\_}{P}}_{{ℂ}^m[\Lambda ]{\tilde{1}}_{b,t}}{J}_{b,t}^{-1}C.$$(A.31)

The new Eq. (A.31) is a representation of the exact information content ofEq. (A.22) regarding the contrast spectrum. Since ${\underset{\_}{P}}_{{ℂ}^m[\Lambda ]\underset{\_}{\tilde{1}}}$ is a projector matrix, it is necessarily of rank < n_λ, with an exact value which depends on the dimension of the subspace generated by the ${\Lambda }^k\underset{\_}{\tilde{1}}$. Thus, Eq. (A.31) is not invertible, and the contrast spectrum cannot be fully recovered from it.

Fortunately, the dataset acquired on β Pictoris b contains several baselines and DITs. For each baseline and each DIT, the projector matrix is different, since it depends on the matrix ${\underset{\_}{\Phi }}_{b,t}$ through the vector ${\underset{\_}{\tilde{1}}}_{b,t}$. These variations can be leveraged to unambiguously recover the complete contrast spectrum C.

The proper way to proceed is to start by extracting a set of linearly independant equations from Eq. (A.31). This can be done by using a diagonal representation of the projector matrix, for example using a singular value decomposition. We can introduce an hermitian matrix ${\underset{\_}{H}}_{b,t}$ and a diagonal matrix D_b,t such that: $$\begin{array}{l}{\underset{\_}{H}}_{b,t}{\underset{\_}{P}}_{{ℂ}^m[\Lambda ]{\underset{\_}{\tilde{1}}}_{b,t}}{\underset{\_}{H}}_{b,t}{}^{†}=D_{b,t}\\ {\underset{\_}{H}}_{b,t}{\underset{\_}{H}}_{b,t}{}^{†}=I.\end{array}$$

Since ${\underset{\_}{P}}_{{ℂ}^m[\Lambda ]{\underset{\_}{\tilde{1}}}_{b,t}}$ is a projector matrix, its eigenvalues are either 1 or 0. We can assume D to be of the form: $$D_{b,t}=\left(\begin{array}{cc}{I}_{r(b,t)}& 0\\ 0& 0\end{array}\right)$$(A.34)

where I_r(b, t) is the identity matrix of size given by the rank of the projector: $r(b,t)=\text{rank}({\underset{\_}{P}}_{{ℂ}^m[\Lambda ]{\underset{\_}{\tilde{1}}}_{b,t}})$.

The matrix ${\underset{\_}{H}}_{b,t}$ can also be written in blocks: $${\underset{\_}{H}}_{b,t}=\left(\begin{array}{cc}{\underset{\_}{H}}_{b,t}^{11}& {\underset{\_}{H}}_{b,t}^{12}\\ {\underset{\_}{H}}_{b,t}^{21}& {\underset{\_}{H}}_{b,t}^{22}\end{array}\right).$$(A.35)

From there, multiplying both sides of Eq. (A.31) by ${\underset{\_}{H}}_{b,t}$, noting that ${\underset{\_}{H}}_{b,t}{\underset{\_}{P}}_{{ℂ}^m[\Lambda ]{\underset{\_}{\tilde{1}}}_{b,t}}=D_{b,t}{\underset{\_}{H}}_{b,t}$, and using a block calculation gives: $$\left[\begin{array}{cc}{\underset{\_}{H}}_{b,t}^{11}& {\underset{\_}{H}}_{b,t}^{12}\end{array}\right]{\widetilde{\underset{\_}{U}}}_{b,t}=\left[\begin{array}{cc}{\underset{\_}{H}}_{b,t}^{11}& {\underset{\_}{H}}_{b,t}^{12}\end{array}\right]J_{b,t}C$$(A.36)

and a dummy equation 0 = 0.

Since r(b, t) < n_λ, the linearly independent system defined by Eq. (A.36) is underdetermined (the matrix $\left[\begin{array}{cc}{\underset{\_}{H}}_{b,t}^{11}& {\underset{\_}{H}}_{b,t}^{12}\end{array}\right]$ has more columns than rows).

To invert the problem, it is necessary to combine all the equations obtained for the different baselines b and t. In matrix notation, this is just a matter of concatenating all the sub-matrices: $$\underset{\_}{H}=\left(\begin{array}{ccc}\left[\begin{array}{cc}{\underset{\_}{H}}_{1,1}^{11}& {\underset{\_}{H}}_{1,1}^{12}\end{array}\right]& & 0\\ & \ddots & \\ 0& & \left[\begin{array}{cc}{\underset{\_}{H}}_{n_b,n_{\text{DIT}}}^{11}& {\underset{\_}{H}}_{n_b,n_{\text{DIT}}}^{12}\end{array}\right]\end{array}\right)$$(A.37) $$\widetilde{\underset{\_}{\mathbb{U}}}=\left(\begin{array}{c}{\widetilde{\underset{\_}{U}}}_{1,1}\\ {\widetilde{\underset{\_}{U}}}_{2,1}\\ ⋮\\ {\widetilde{\underset{\_}{U}}}_{n_b,n_{\lambda }}\end{array}\right)$$(A.38) $$\mathbb{J}=\left(\begin{array}{c}{J}_{1,1}{}^{-1}\\ J_{2,1}{}^{-1}\\ ⋮\\ J_{n_b,n_{\lambda }}{}^{-1}\end{array}\right)$$(A.39)

This gives: $$\underset{\_}{H}\widetilde{\underset{\_}{\mathbb{U}}}=\underset{\_}{H}\mathbb{J}C$$(A.40)

Which hasthe form of a linear problem: $$\underset{\_}{\mathbb{Y}}=\underset{\_}{ℍ}C,$$(A.41)

where $\underset{\_}{\mathbb{Y}}$ is a linear transformation of the calibrated visibilities defied by: $$\mathbb{Y}=\underset{\_}{H}\widetilde{\underset{\_}{\mathbb{U}}}.$$(A.42)

And $\underset{\_}{ℍ}$ is the collapsed matrix: $$\underset{\_}{ℍ}=\left(\begin{array}{c}\left[\begin{array}{cc}{\underset{\_}{H}}_{1,1}^{11}& {\underset{\_}{H}}_{1,1}^{12}\end{array}\right]J_{1,1}^{-1}\\ \left[\begin{array}{cc}{\underset{\_}{H}}_{2,1}^{11}& {\underset{\_}{H}}_{2,1}^{12}\end{array}\right]J_{2,1}^{-1}\\ ⋮\\ \left[\begin{array}{cc}{\underset{\_}{H}}_{n_b,n_{\text{DIT}}}^{11}& {\underset{\_}{H}}_{n_b,n_{\text{DIT}}}^{12}\end{array}\right]J_{n_b,n_{\text{DIT}}}^{-1}\end{array}\right)$$(A.43)

The problem defined by Eq. (A.41) can be solved using the maximum likelihood formalism, adapted to complex random variables. The expression of the maximum likelihood solution is given by Eq. (A.9). For the contrast spectrum, we have: $$\widehat{C}={\left(\text{Re}\left({\tilde{\mathbb{Y}}}_2^{†}{\underset{\_}{W}}_2^{-1}{\underset{\_}{ℍ}}_2\right){\left[\text{Re}\left({\underset{\_}{ℍ}}_2^{†}{\underset{\_}{W}}_2^{-1}{\underset{\_}{ℍ}}_2\right)\right]}^{-1}\right)}^T$$(A.44)

where the extended vectors and matrices ${\widetilde{\underset{\_}{Y}}}_2$, ${\underset{\_}{ℍ}}_2$, and ${\underset{\_}{W}}_2$, are defined by: $$\widetilde{{\underset{\_}{Y}}_2}=\left(\begin{array}{c}\widetilde{\underset{\_}{Y}}\\ {\widetilde{\underset{\_}{Y}}}^{*}\end{array}\right){\underset{\_}{ℍ}}_2=\left(\begin{array}{c}\underset{\_}{ℍ}\\ {\underset{\_}{ℍ}}^{*}\end{array}\right){\underset{\_}{W}}_2=\left(\begin{array}{cc}\underset{\_}{W}& \underset{\_}{Z}\\ {\underset{\_}{Z}}^{†}& {\underset{\_}{W}}^{*}\end{array}\right).$$(A.45)

The uncertainty on this best estimate of the contrast spectrum can be obtained with a direct error propagation all the way to the real covariance matrix on $\widehat{C}$. The covariance and pseudo covariance matrices of ${\underset{\_}{\mathbb{Y}}}_2{}^{†}{\underset{\_}{W}}_2^{-1}{\underset{\_}{ℍ}}_2$ are given by: $$\begin{array}{l}\text{cov}\left({\underset{\_}{\mathbb{Y}}}_2{}^{†}{\underset{\_}{W}}_2^{-1}{\underset{\_}{ℍ}}_2\right)={\underset{\_}{ℍ}}_2{}^{†}{\underset{\_}{W}}_2^{-1}{}^{†}\text{cov}\left({\underset{\_}{\mathbb{Y}}}_2\right){\underset{\_}{W}}_2^{-1}{\underset{\_}{ℍ}}_2\\ \text{pcov}\left({\underset{\_}{\mathbb{Y}}}_2{}^{†}{\underset{\_}{W}}_2^{-1}{\underset{\_}{ℍ}}_2\right)={\underset{\_}{ℍ}}_2{}^T{\underset{\_}{W}}_2^{-1}{}^T\text{pcov}\left({\underset{\_}{\mathbb{Y}}}_2\right){\underset{\_}{W}}_2^{-1}{\underset{\_}{ℍ}}_2.\end{array}$$

A proper combination of these covariance and pseudo-covariance matrices gives the covariance matrix of the real part of ${\underset{\_}{\mathbb{Y}}}_2{\underset{\_}{W}}_2^{-1}{\underset{\_}{ℍ}}_2$: $$\left\{\text{Re}\left({\underset{\_}{\mathbb{Y}}}_2{}^{†}{\underset{\_}{W}}_2^{-1}{\underset{\_}{ℍ}}_2\right)\right\}=\frac{1}{2}\text{Re}\left(\text{cov}\left({\underset{\_}{\mathbb{Y}}}_2{}^{†}{\underset{\_}{W}}_2^{-1}{\underset{\_}{ℍ}}_2\right)+\text{pcov}\left({\underset{\_}{\mathbb{Y}}}_2{}^{†}{\underset{\_}{W}}_2^{-1}{\underset{\_}{ℍ}}_2\right)\right),$$(A.48)

which canthen be propagated to give the final full covariance matrix on $\widehat{C}$: $$\left(\widehat{C}\right)=\text{Re}{\left({\underset{\_}{ℍ}}_2{}^{†}{\underset{\_}{W}}_2^{-1}{\underset{\_}{ℍ}}_2\right)}^{-1}{}^T\text{cov}\left(\text{Re}\left({\tilde{\mathbb{Y}}}_2^{†}{\underset{\_}{W}}_2^{-1}{\underset{\_}{ℍ}}_2\right)\right)\text{Re}{\left({\underset{\_}{ℍ}}_2{}^{†}{\underset{\_}{W}}_2^{-1}{\underset{\_}{ℍ}}_2\right)}^{-1}.$$(A.49)

Appendix B Posterior of the petitRADTRANS fit

In Fig. B.1 we show the posteriors of the GRAVITY+GPI fit with petitRADTRANS, described in Sect. 4.3. Panel a shows the corner plot for all but the temperature nuisance parameters. Panel b shows the retrieved temperature uncertainty envelopes.

The parameters shown are the following: the C/O, adjusted by changing the oxygen abundance at a given [Fe/H]. The metallicity [Fe/H], which was used to scale the number fraction of all atomic elements (except H and He) by 10^[Fe/H]. The quench pressure P_quench of the atmosphere, here converged to a low enough value such that non-equilibrium chemistry effects are negligible. The mass fractions of Fe and MgSiO₃ at the cloud base. Here they are expressed in units of a log-ed decrease factor, which gets multiplied with the maximally allowed mass fraction, based on the elemental composition of the atmosphere. The stoichiometric factors of MgSiO₃ are used for finding this upper limit. The cloud settling parameter f_sed, as described in Sect. 4.3. The log-ed eddy diffusion coefficient K_zz, in units of cm² s⁻¹. This is used for calculating the cloud particle size, as described in Sect. 4.3. The planet’s surface gravity log (g). The planetary radius R_P, in units of Jupiter radii. The planetary mass M_P, in units of Jupiter masses. This is calculated using the sampled log (g) and R_P values. The width of the log-normal cloud particle size distribution σ_g, as described in Sect. 4.3. The vertically constant mass fraction of FeH, expressed in units of a log-ed decrease factor. This factor gets multiplied with the maximally allowed mass fraction of FeH, based on the elemental composition of the atmosphere and the Fe atoms not yet incorporated into the Fe clouds. Finally the f_GPI-Y, f_GPI-J, f_GPI-H and f_GRAV factors describe the multiplicative scaling of the individual bands, which were allowed to vary by 50% in the case of GPI, and by 5% in the case of GRAVITY.

Fig. B.1 Panel a: projected 2-d posterior of the GRAVITY+GPI fit with petitRADTRANS (spectrum shown in Fig. 5), described in Sect. 4.3. See the text in Appendix B for a description of the parameters. Panel b: pressure-temperature envelopes obtained for the same retrieval. At every pressure, we plot the 16 to 84-percentile envelopes in dark blue, and the 2.5 to 97.5 percentile envelopes in light blue. If the temperature values were following a Gauss distribution, this would correspond to the 1 and 2 σ envelopes, respectively.

The temperature envelopes in panel b of Fig. B.1 are obtained by plotting, at every pressure, the 16 to 84-percentile envelopes, and the 2.5 to 97.5 percentile envelopes. If the temperature values at a given pressure layer were following a normal distribution, this would correspond to the 1 and 2 σ envelopes, respectively.

References