© ESO 2020

1. Introduction

Giant planets are born in circumstellar disks from the material that remains after the formation of stars. The processes by which they form remain unclear and two schools of thought are currently competing to formulate a valid explanation, based on two scenarios: (1) disk instability, which states that planets form through the collapse and fragmentation of the circumstellar disk (Boss 1997; Cameron 1978), and (2) core accretion, which states that a planetary core forms first through the slow accretion of solid material and later captures a massive gaseous atmosphere (Pollack et al. 1996; Lissauer & Stevenson 2007).

These two scenarios were initially thought to lead to very different planetary masses and luminosities, with planets formed by gravitational instability being much hotter and having higher post-formation luminosity and entropy than planets formed by core accretion (Marley et al. 2007). With such a difference in the post-formation entropy, the two scenarios could be distinguished using mass and luminosity measurements of young giant planets. However, several authors have since shown that so-called ‘hot start’ planets are not incompatible with a formation model of core accretion, provided that the physics of the accretion shock (Marleau et al. 2017, 2019) and the core-mass effect (a self-amplifying process by which a small increase of the planetary core mass weakens the accretion shock during formation and leads to higher post-formation luminosity) is properly accounted for (Mordasini 2013; Mordasini et al. 2017). In this situation, mass and luminosity measurements are of prime interest as they hold important information regarding the physics of the formation processes.

The most efficient method used thus far to determine the masses of giant planets is radial velocity measurements of the star. This method has one significant drawback: it is only sensitive to the product, msin(i), where m is the mass of the planet and i its orbital inclination. In parallel, the most efficient method used to measure the luminosities of giant planets is direct imaging using dedicated high-contrast instruments. Since direct imaging also provides a way to estimate the orbital inclination when the period allows for a significant coverage of the orbit, the combination of radial velocity and direct imaging can, in principle, break the msin(i) degeneracy and enable an accurate measurement of masses and luminosities. But therein lies the rub: while radial-velocity is sensitive to planets orbiting close-in (typically < 1 au) around old (>1 Gyr) stars, direct imaging is sensitive to planets orbiting at much larger separations (≥10 au) around younger stars (<100 Myr). Thus, these techniques are not easily combined for an analysis of a single object.

Significant efforts have been made over the past few years to extend the radial-velocity method to longer periods and younger stars. This is illustrated by the recent detection of β Pic c (Lagrange et al. 2019). The recent characterisation of HR 8799 e also demonstrated the potential of interferometric techniques to determine orbits, luminosities (GRAVITY Collaboration 2019), and atmospheric properties (Mollière et al. 2020) of imaged extra-solar planets with short periods. In this Letter, we attempt to bridge the gap between the two techniques by reporting on the direct confirmation of a planet discovered by radial velocity: β Pic c.

2. Observations and detection of β Pic c

2.1. Observations

We observed β Pic c during the night of the 9^th and 11^th of February 2020, as well as that of the 7^th of March 2020 with the Very Large Telescope Interferometer (VLTI), using the four 8.2 m Unit Telescopes (UTs), and the GRAVITY instrument (GRAVITY Collaboration 2017). The first detection was obtained from the allocated time of the ExoGRAVITY large program (PI Lacour, ID 1104.C-0651). A confirmation was obtained two days later courtesy of the active galactic nucleus (AGN) large programme (PI Sturm, ID 1103.B-0626). The final dataset was obtained with the Director’s Discretionary Time (DDT; ID 2104.C-5046). Each night, the atmospheric conditions ranged from good (0.9″) to excellent (0.4″). The observing log is presented in Table 1, and the observing strategy is described in Appendix A. In total, three hours of integration time were obtained in K-band (2.2 μm) at medium resolution (R = 500).

2.2. Data reduction and detection

The data reduction used for β Pic c follows what has been developed for β Pic b (GRAVITY Collaboration 2020) and is detailed in Appendix B.

An initial data reduction using the GRAVITY pipeline (Lapeyrere et al. 2014) gives the interferometric visibilities obtained on the planet, and on the star, phase-referenced with the metrology system (i.e. in which the atmospheric and instrumental phase distortions are reduced to an unknown constant). The data reduction then proceeds in two steps, namely an initial extraction of the astrometry followed by the extraction of the spectrum. Since the planet and the star are observed successively with the instrument, the data reduction yields a contrast spectrum, defined as the ratio of the planet spectrum to the star spectrum. This observable is more robust to variations of the instrument and atmospheric transmission. The overall process yields a periodogram power map in the field-of-view of the science fiber (see Fig. 1), and a contrast spectrum. The astrometry (see Table 2) is extracted from the periodogram power maps by taking the position of the maximum of the periodogram power, z, and error bars are obtained by breaking each night into individual exposures (11 to 12 per night, see Table 1), so as to estimate the effective standard-deviation from the data themselves.

Fig. 1. Detection of β Pictoris c. The three panels show the periodogram power maps calculated over the fiber field-of-view for each of the three observation epochs and following the subtraction of stellar residuals. The presence of a peak in the power maps indicates a point-like source in the field-of-view of the fiber, with side lobes that are characteristic of the interferometric nature of the observations.

The planet spectrum is obtained by combining the three K-band contrast spectra obtained at each epoch (no evidence of variability was detected) and multiplying them by a NextGen stellar model (Hauschildt et al. 1999) corresponding to the known stellar parameters (temperature T = 8000 K, surface gravity log(g) = 4.0 dex, Zorec & Royer 2012) and a solar metallicity, scaled to an ESO K-band magnitude of 3.495 (van der Bliek et al. 1996). We note that at the resolution of GRAVITY in the K-band, the temperature, surface gravity, and metallicity of the star all have only a marginal impact on the final spectrum. The resulting flux calibrated K-band spectrum is presented in Fig. 2.

Fig. 2. K-band spectrum of β Pic c, with the best fit obtained with the Exo-REM and Drift-Phoenix grids overplotted (see Sect. 4).

3. The orbit of β Pic c

The relative astrometry is plotted in Fig. 3. We used a Markov chain Monte Carlo (MCMC) sampler to determine the orbital parameters of both planets, as well as the mass of each object. The MCMC analysis was done using the orbitize!¹ (Blunt et al. 2020) software. For the relative astrometry, we included all previous direct imaging and GRAVITY data. The NACO, SPHERE, and GRAVITY observations of β Pic b are from Lagrange et al. (2020) and reference therein. The existing GPI data are summarized by Nielsen et al. (2020). Of course, we added the three GRAVITY detections of β Pic c from Table 2. For the absolute astrometry of the star, we used the HIPPARCOS Intermediate Astrometric Data (van Leeuwen 2007) and the DR2 position (Gaia Collaboration 2018). Lastly, for the radial velocity, we used the 12 years of HARPS measurements (Lagrange et al. 2019, 2020).

Fig. 3. Motion of planets b and c around β Pictoris, as observed by GRAVITY (respectively red and green dots) and direct imaging (orange dots). The inset panels are zooms on the GRAVITY observations of β Pic c. Within the insets, the green ellipses correspond to the 1σ error intervals, typically between 100 to 200 μas, elongated because of the shape of the interferometric array. The white dots are the retrieved position at the time of observation from the posteriors in Table 3.

The priors of the MCMC analysis are listed in Table 3. Most of the priors are uniform or pseudo-uniform distributions. We have made sure that all uniform distributions have limits many sigma outside the distributions of the posteriors. We used two Gaussian priors to account for independent knowledge. Firstly, the knowledge of the mass of b that comes from the spectrum of the atmosphere (GRAVITY Collaboration 2020): 15.4 ± 3.0 M_Jup. Secondly, an estimation of the stellar mass of 1.75 ± 0.05 M_⊙ (Kraus et al. 2020).

The semi-major axis of β Pic c’s orbit is 2.72 ± 0.02 au, with an orbital period of β Pic c is 3.37 ± 0.04 year. The planet is currently moving toward us, with a passage of the Hill sphere between Earth and the star in mid-August 2021 (±1.5 month). This means that the transit of circumplanetary material could be seen at that time. However, no sign of a photometric event was detected at the time of the previous conjunction, which was between February−June 2018 (Zwintz et al. 2019).

The position angles of the ascending nodes of b and c are 31.82 ± 0.02° and 30.98 ± 0.12°, respectively. Their inclinations are 88.99 ± 0.01° and 89.17 ± 0.50°. This means that the two planets are co-planar to within less than a degree, perpendicular from the angular momentum vector of the star (Kraus et al. 2020). The eccentricity of c from radial velocity alone is 0.29 ± 0.05 (Lagrange et al. 2020). By including the new relative astrometry, we find a higher value, which is consistent within the error bars, of 0.37 ± 0.12. In combination with the eccentricity of b (0.10 ± 0.01), the system has a significant reservoir of eccentricity.

Using the criterion for stability from Petrovich (2015), we establish that the system is likely to be stable for at least 100 Myr. Using an N-body simulation (Beust 2003, SWIFT HJS), we confirm this stability for 10 million years for the peak of the orbital element posteriors. Due to the high mass ratios and eccentricities, the periodic secular variations of the inclinations and eccentricities are not negligible on this timescale. The periods associated with the eccentricities are around 50 000 year and can trigger variations of 0.2 in the eccentricities of both bodies, while the period associated with the relative inclination is around 10 000 year with variations by 20%. Moreover, instabilities that could produce eccentricity may still happen on longer time scales, within a small island of the parameter space or by interaction with smaller planetary bodies. One possibility, for example, is that the planetary system is undergoing dynamical perturbations from planetesimals, as our own solar system may well have undergone in the past (the Nice model, Nesvorný 2018).

With a knowledge of the inclination, the mass of β Pic c can be determined from the radial velocity data. In addition, the mass of c can theoretically be obtained purely from the astrometry of b. This can be explained by the fact that the orbital motion of c affects the position of the central star, which sets the orbital period of c on the distance between b and the star. The amplitude of the effect is twice 2.72 × (M_c/M_star), equivalent to ≈1200 μas. This is many times the GRAVITY error bar, so this effect must be accounted for in order to estimate the orbital parameters of b (see Appendix C). In a few years, once the GRAVITY data have covered a full orbital period of c, this effect will give an independent dynamical mass measurement of c. For now, the combination of RV and astrometry gives a robust constrain of M_c = 8.2 ± 0.8 M_Jup, compatible with the estimations by Lagrange et al. (2020). We observed a mass for β Pic b (M_b = 9.0 ± 1.6 M_Jup) that is 2σ from the mean of the Gaussian prior. If we remove the prior from the analysis, we find a significantly lower mass: M_b = 5.6 ± 1.5 M_Jup. Therefore, the radial velocities bias the data toward a lighter planet, discrepant with other analyses using only absolute astrometry: 11 ± 2 M_Jup in Snellen & Brown (2018), in Kervella et al. (2019) and in Nielsen et al. (2020). These three results, however, should be taken with caution; the analyses were carried out based on a single planet. In our MCMC analysis with two planets, we could not reproduce these high masses. This is because the planet c is responsible for a large quantity of the proper motion observed with HIPPARCOS, converging to a low mass for β Pic b.

4. Physical characteristics of β Pic c

4.1. Atmospheric modelling

The spectrum extracted from our GRAVITY observations can be used in conjunction with atmospheric models to constrain some of the main characteristics of the planet; the overall K-band spectral shape is driven mostly by the temperature and, to a lesser extent, by the surface gravity and metallicity. Since the GRAVITY spectrum is flux-calibrated, the overall flux level also constrains the radius of the planet.

Using the Exo-REM atmospheric models (Charnay et al. 2018) and the species² atmospheric characterization toolkit (Stolker et al. 2020), we find a temperature of T = 1250 ± 50 K, a surface gravity of , and a radius of 1.2 ± 0.1 R_Jup.

The Drift-Phoenix model (Woitke & Helling 2003; Helling & Woitke 2006; Helling et al. 2008) gives a slightly higher temperature (1370 ± 50 K) and a smaller radius (1.05 ± 0.1 R_Jup). The surface-gravity is similar to the one obtained with Exo-REM with larger errors (log(g) = 4.1 ± 0.5 dex).

The best fits obtained with Exo-REM and Drift-Phoenix are overplotted to the GRAVITY spectrum in Fig. 2. We note that at this temperature and surface gravity, both the Exo-REM and Drift-Phoenix models are relatively flat and do not show any strong molecular feature. This is also the case of the GRAVITY spectrum, which remains mostly flat over the entire K-band. The only apparent spectral feature is located at ≃2.075 μm. It is not reproduced by the models, and is most probably of telluric origin.

Interestingly, the temperature and surface gravity derived from the atmospheric models (which are agnostic regarding the formation of the planet) are relatively close to the predictions of the ‘hot start’ AMES-Cond evolutionary tracks (Baraffe et al. 2003). For a mass of 8.2 ± 0.8 M_Jup and an age of 18.5 ± 2.5 Myr (Miret-Roig et al. 2020), the AMES-Cond tracks predict a temperature of T = 1340 ± 160 K and a surface gravity of log(g) = 4.05 ± 0.05.

4.2. Mass/luminosity and formation of β Pic c

Such hot young giant planets ( >  1000 K) have historically been associated with formation through disk instability (Marley et al. 2007). However, disk fragmentation is unlikely to occur within the inner few au (Rafikov 2005) and, in that respect, the existence of a planet as massive as β Pic c at only 2.7 au poses a real challenge.

One possibility is that the planet formed in the outer β Pic system (⪆30 au) and later underwent an inward type II migration that brought it to its current orbit. From a dynamical standpoint, however, the plausibility of such a scenario given the existence of a second giant planet in the system remains to be proven.

An alternative to the formation of β Pic c by disk instability is core-accretion. In recent years, a number of authors have shown that it is possible to obtain hot and young giant planets with such a scenario if the core is sufficiently massive (Mordasini 2013), if the accretion shock is relatively inefficient at cooling the incoming gas (Mordasini et al. 2017; Marleau & Cumming 2014), or if the shock dissipates sufficient heat into the in-falling gas (Marleau et al. 2017, 2019).

A possible argument in favor of a core accretion scenario comes from the marginal discrepancy in luminosity between the AMES-Cond ‘hot-start’ model and our observations. The spectrum given in Fig. 2 is flux-calibrated and can be integrated over the K-band to give the apparent magnitude of the planet. To calculate this magnitude, we first interpolate the ESO K-band transmission curve (van der Bliek et al. 1996) over the GRAVITY wavelength bins to obtain a transmission vector, T. Then, denoting S as the flux-calibrated planet spectrum and W as the associated covariance matrix, the total flux in the ESO K-band is: F = T^TS. The error bar is given by: ΔF = TWT^T. The fluxes are then converted to magnitudes using the proper zero-points (van der Bliek et al. 1996). This gives a K-band magnitude for β Pic c of m_K = 14.3 ± 0.1 (Δm_K = 10.8 ± 0.1). At the distance of the β Pictoris system (19.44 ± 0.05 pc, van Leeuwen 2007), this corresponds to an absolute K-band magnitude of M_K = 12.9 ± 0.1.

For comparison, the AMES-Cond model predicts a K-band magnitude of M_K = 12.3 ± 0.5, which is marginally discrepant with regard to our measurement. This means that AMES-Cond overpredicts either the temperature or the radius of the planet. The fact that the Exo-REM atmospheric modelling agrees with the AMES-Cond tracks on the temperature suggests that the AMES-Cond tracks overpredict the radius. Indeed, our atmospheric modelling yields an estimate of 1.2 ± 0.1 R_Jup for the radius, which is discrepant with regard to the AMES-Cond prediction (R = 1.4 ± 0.05 R_Jup).

The position of β Pic c in mass/radius and mass/luminosity is shown in Fig. 4, together with the 18.5 ± 2.5 Myr AMES-Cond isochrones, and two synthetic populations of planets formed through core nucleated accretion³. These two synthetic populations take into account the core-mass effect (Mordasini 2013). They correspond to what is commonly referred to as ‘warm starts’, with their planets much more luminous that in the classical ‘cold start’ scenario (Marley et al. 2007). Compared to the hot-start population, at M = 8 M_Jup, the hot (resp. nominal) core accretion is 30% less luminous (resp. 50%) and has radii that are 5% smaller (resp. 10%). For completeness, we carried out the exact same analysis on β Pic b, using the published GRAVITY K-band spectrum (GRAVITY Collaboration 2020), and we over-plotted the result in Fig. 4. Since the dynamical mass of β Pic b remains poorly constrained, we used a value of 15.4 ± 3 M_Jup, based on atmospheric modelling (GRAVITY Collaboration 2020). The large error bars on the mass of β Pic b makes it equally compatible with AMES-Cond and core-accretion, both in terms of mass/luminosity and mass/radius. This is not the case for β Pic c. Indeed, even if the new GRAVITY measurements presented in this Letter are not yet sufficient to confidently reject a hot-start model such as AMES-Cond, the position of the planet in mass/radius, and, to a lesser extent, in mass/luminosity, seems to be in better agreement with a core-accretion model, particularly when the core-mass effect is taken into account. Interestingly, a similar conclusion has been drawn for β Pic b (GRAVITY Collaboration 2020) based primarily on the estimation of its C/O abundance ratio.

Fig. 4. Position of β Pic b & c in mass/luminosity and mass/radius. The two panels give mass/radius (left-hand side) and mass/radius (right-hand side) diagrams showing the positions of the two known giant planets of the β Pic system and the predictions of planet formation models. The predictions of the AMES-Cond ‘hot start’ model, which extends to ≃14 MJup, are depicted as a red line (with the shaded area corresponding to the uncertainty on the age). Two synthetic populations generated by core accretion are also shown: CB753, corresponding to core accretion with cold nominal accretion and CB752, with hot accretion. These two populations are valid for 20 Myr and take into account the core-mass effect (Mordasini 2013).

5. Conclusion

Thanks to this first detection of β Pic c, we can now constrain the inclination and the luminosity of the planet. The inclination, in combination with the radial velocities, gives a robust estimate of the mass: 8.2 ± 0.8 M_Jup. On the other hand, the mass of β Pic b is not as well-constrained: 9.0 ± 1.6 M_Jup.

To match these masses with the hot start scenario, the hot-start AMES-Cond models show that a bigger and brighter β Pic c would be needed. Our first conclusion is that the formation of β Pic c is not due to gravitational instability but, more likely, to a warm-start core accretion scenario. Our second conclusion is that the mass of β Pic b should be revised in the future, once the radial velocity data covers a full orbital period.

Our final conclusion is that we are now able to directly observe exoplanets that have been detected by radial velocities. This is an important change in exoplanetary observations because it means we can obtain both the flux and dynamical masses of exoplanets. It also means that we will soon be able to apply direct constrains on the exoplanet formation models.

Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme, grant agreement No. 639248 (LITHIUM), 757561 (HiRISE) and 832428 (Origins). A.A. and P.G. were supported by Fundação para a Ciência e a Tecnologia, with grants reference UIDB/00099/2020 and SFRH/BSAB/142940/2018.

References

Appendix A: Observations

GRAVITY is a fiber-fed interferometer that uses a set of two single-mode fibers for each telescope of the array. One fiber feeds the fringe tracker (Lacour et al. 2019), which is used to compensate for the atmospheric phase variations, and the second fiber feeds the science channel (Pfuhl et al. 2014). The observations of β Pictoris c presented in this work were obtained using the on-axis mode, meaning that a beamsplitter was in place to separate the light coming from each telescope. The instrument was also used in dual-field mode, in which the two fibers feeding the fringe tracking and science channels can be positioned independently.

The use of the on-axis/dual-field mode is crucial for observing of exoplanets with GRAVITY, as it allows the observer to use a strategy in which the science fiber is alternately set on the central star for phase and amplitude referencing, and at an arbitrary offset position to collect light from the planet (see illustration in Fig. A.1). Since the planet itself cannot be seen on the acquisition camera, its position needs to be known in advance in order to properly position the science fiber. For our observations, we used the most up-to-date predictions for β Pic c (Lagrange et al. 2020). The use of single-mode fibers severely limits the field-of-view of the instrument (to about 60 mas), and any error on the position of the fiber results in a flux loss at the output of the combiner (half of the flux is lost for a pointing error of 30 mas). However, a pointing error does not result in a phase error, and therefore does not affect the astrometry or spectrum extraction other than by decreasing the signal-to-noise ratio (S/N). In practice, any potential flux error is also easily corrected once the exact position of the planet is known (see Appendix B).

Fig. A.1. Schematic representation of the positions of the two fibers when observing in on-axis and dual-field mode with VLTI/GRAVITY. In on-axis mode, a beamsplitter is used to separate the incoming light in amplitude. In dual-field, the fringe-tracking fiber and the science fiber can be moved separately, and can be centered on different targets. For observing exoplanets, two types of configuration are used, with the science fiber either centered on the star, for calibration, or centered at the expected location of the planet. The fringe-tracking fiber is always centered on the star, whose light is used as a reference for the phase of the interferometric visibilities.

Appendix B: Data reduction

The theoretical planetary visibility is given as follows:

(B.1)

where S_planet is the spectrum of the planet, U, V are the u − v coordinates of the interferometric baseline b, (Δα, Δδ) the relative astrometry (in RA/Dec), λ the wavelength, and t the time.

In practice, however, two important factors need to be taken into account: (1) atmospheric and instrumental transmission distort the spectrum, and (2) even when observing on-planet, visibilities are dominated by residual starlight leaking into the fiber. A better representation of the visibility actually observed is given by (GRAVITY Collaboration 2020):

(B.2)

in which G is a transmission function accounting for both the atmosphere and the instrument, and the Q(b, t, λ) are 6th-order polynomial functions of λ (one per baseline and per exposure) accounting for the stellar flux.

In parallel, the on-star observations give a measurement of the stellar visibility, also affected the atmospheric and instrumental transmission:

(B.3)

where the last equality comes from the fact that the star itself is the phase-reference of the observations, and thus V_star = S_star.

Here, if we introduce the contrast spectrum C = S_planet/S_star, we can get rid of the instrumental and atmospheric transmission:

(B.4)

Since the λ → Q(b, t, λ) are polynomial functions, the model is non-linear only in the two parameters, Δα, Δδ.

In reality, since the on-star and on-planet data are not acquired simultaneously, Eq. (B.3) should be written for t^* ≠ t. This means that the factor G from Eqs. (B.2) and (B.3) do not cancel out perfectly, and an extra factor G(t)/G(t^*) should appear in Eq. (B.4). The two sequences are separated by 10 min at most, and the instrument is stable on such durations. Therefore, the main contribution to this factor G(t)/G(t^*) comes from atmospheric variations. In order to limit the impact of these atmospheric variations on the overall contrast level, we use the ratio of the fluxes observed by the fringe-tracking fiber during the on-planet observation and the on-star observation as a proxy for this ratio of G. The fringe-tracking combiner works similarly to the science combiner, and acquires data at the same time (though at a higher rate). The main difference is that the fringe-tracking fiber is always centered on the star. Consequently, the ratio of the fluxes observed by the fringe-tracking fiber is a direct estimate of G(t)/G(t^*). Since the fringe-tracking channel has a much lower resolution that the science channel, this only provides a correction of the overall contrast level (integrated in wavelength).

Similarly, a difference in G can occur if the fiber is misaligned with the planet (the fiber is always properly centered on the star, which is observable on the acquisition camera of the instrument). This effect is easily corrected once the correct astrometry of the planet is known, by taking into account the theoretical fiber injection in Eqs. (B.2) and (B.3). In practice, the initial pointing of the fiber is good and this effect remains small (<3%).

To assess the presence of a planet in the data, we first look for a point-like source in the field-of-view of the fiber, assuming a constant contrast over the K-band: C(λ) = c. To do so, we switch to a vector representation, and write V_onplanet = (V_onplanet(b, t, λ))_b, t, λ the vector obtained by concatenating all the datapoints (for all the six baselines, all the exposures, and all the wavelength grid points). Similarly, for a given set of polynomials Q, and for given values of c, Δα, and Δδ, we derive V_Δα, Δδ[Q, c], the model for the visibilities described in Eq. (B.4). The χ² sum is then written:

(B.5)

with W the covariance matrix affecting the projected visibilities .

This χ², minimized over the nuisance parameters, Q and c, is related to the likelihood of obtaining the data given the presence of a planet at Δα, Δδ. The value at Δα = 0 and Δδ = 0 corresponds to the case where no planet is actually injected in the model (the exponential in Eq. (B.4) is flat) and can be taken as a reference: . The values of the χ² using a planet model can be compared to this reference, by defining:

(B.6)

The quantity z(Δα, Δδ) can be understood as a Bayes factor, comparing the likelihood of a model that includes a planet at (Δα, Δδ) to the likelihood of a model without any planet. It is also a direct analog of the periodogram power used, for example, in the analysis of radial-velocity data (Scargle 1982; Cumming 2004).

The resulting power maps obtained for each of the three nights of β Pic c observation are given in Fig. 1. The astrometry is extracted from these maps by taking the position of the maximum of z, and the error bars are obtained by breaking each night into individual exposures so as to estimate the effective standard-deviation from the data themselves.

Once the astrometry is known, Eq. (B.4) can be used to extract the contrast spectrum C (GRAVITY Collaboration 2020). The extraction of the astrometry using relies on the implicit assumption that the contrast is constant over the wavelength range. To mitigate this, the whole procedure (astrometry and spectrum extraction) is iterated once more starting with the new contrast spectrum, to check for consistency of results.

Appendix C: Note on multi-planet orbit fitting

For the MCMC analysis, we used a total of 100 walkers performing 4000 steps each. A random sample of 300 posteriors are plotted in Fig. C.1. The two lower panels show the trajectory of both planets with respect to the star. Both trajectories are not perfectly Keplerian because each planet influence the position of the star. Therefore, the analysis must be global.

Fig. C.1. β Pic planetary system: absolute astrometry, radial velocities, and relative astrometry observations. Absolute astrometry is showed in the upper panel as deviation from position RA = 5:47:17.08346 Dec = −51:04:00.1699 at year 1991.25, with a proper motion as in Table 3. Because individual measurements of HIPPARCOS data are 1 dimension only, they are projected in this display along the axis of the β Pictoris system, ie. 31.8 deg. Radial velocity in the second panel includes the v0 term of −23.09 m s−1 (dotted line). Two lower panels: orbital motion of the planets with respect to the star, also projected on an axis inclined by 31.8 deg to the east of the north.

The orbitize! (Blunt et al. 2020) software does not permit N-body simulations that account for planet-planet interactions. However, it allows for the simultaneous modelling of multiple two-body Keplerian orbits. Each two-body Keplerian orbit is solved in the standard way by reducing it to a one-body problem where we solve for the time evolution of the relative offset of the planet from the star instead of each body’s orbit about the system barycenter. As with all imaging astrometry techniques, the GRAVITY measurements are relative separations, so it is convenient to solve orbits in this coordinate system. However, as each planet orbits, it perturbs the star which itself orbits around the system barycenter. The magnitude of the effect is proportional to the offset of each planet times M_planet/M_tot where M_planet is the mass of the planet and M_tot is the total mass of all bodies with separation less than or equal to the separation of the planet. Essentially, the planets cause the star to wobble in its orbit and introduce epicycles in the astrometry of all planets relative to their host star. We modelled these mutual perturbations on the relative astrometry in orbitize! and verified it against the REBOUNDN-body package that simulates mutual gravitational effects (Rein & Liu 2012; Rein & Spiegel 2015). We found a maximum disagreement between the packages of 5 μas over a 15 year span for β Pic b and c. This is an order of magnitude smaller than our error bars, so our model of Keplerian orbits with perturbations is a suitable approximation. We are not yet at the point we need to model planet-planet gravitational perturbations that cause the planets’ orbital parameters to change in time.

All Tables

All Figures

Fig. 1. Detection of β Pictoris c. The three panels show the periodogram power maps calculated over the fiber field-of-view for each of the three observation epochs and following the subtraction of stellar residuals. The presence of a peak in the power maps indicates a point-like source in the field-of-view of the fiber, with side lobes that are characteristic of the interferometric nature of the observations.

In the text

	Fig. 2. K-band spectrum of β Pic c, with the best fit obtained with the Exo-REM and Drift-Phoenix grids overplotted (see Sect. 4).
In the text

Fig. 3. Motion of planets b and c around β Pictoris, as observed by GRAVITY (respectively red and green dots) and direct imaging (orange dots). The inset panels are zooms on the GRAVITY observations of β Pic c. Within the insets, the green ellipses correspond to the 1σ error intervals, typically between 100 to 200 μas, elongated because of the shape of the interferometric array. The white dots are the retrieved position at the time of observation from the posteriors in Table 3.

In the text

Fig. 4. Position of β Pic b & c in mass/luminosity and mass/radius. The two panels give mass/radius (left-hand side) and mass/radius (right-hand side) diagrams showing the positions of the two known giant planets of the β Pic system and the predictions of planet formation models. The predictions of the AMES-Cond ‘hot start’ model, which extends to ≃14 MJup, are depicted as a red line (with the shaded area corresponding to the uncertainty on the age). Two synthetic populations generated by core accretion are also shown: CB753, corresponding to core accretion with cold nominal accretion and CB752, with hot accretion. These two populations are valid for 20 Myr and take into account the core-mass effect (Mordasini 2013).

In the text

Fig. A.1. Schematic representation of the positions of the two fibers when observing in on-axis and dual-field mode with VLTI/GRAVITY. In on-axis mode, a beamsplitter is used to separate the incoming light in amplitude. In dual-field, the fringe-tracking fiber and the science fiber can be moved separately, and can be centered on different targets. For observing exoplanets, two types of configuration are used, with the science fiber either centered on the star, for calibration, or centered at the expected location of the planet. The fringe-tracking fiber is always centered on the star, whose light is used as a reference for the phase of the interferometric visibilities.

In the text

Fig. C.1. β Pic planetary system: absolute astrometry, radial velocities, and relative astrometry observations. Absolute astrometry is showed in the upper panel as deviation from position RA = 5:47:17.08346 Dec = −51:04:00.1699 at year 1991.25, with a proper motion as in Table 3. Because individual measurements of HIPPARCOS data are 1 dimension only, they are projected in this display along the axis of the β Pictoris system, ie. 31.8 deg. Radial velocity in the second panel includes the v0 term of −23.09 m s−1 (dotted line). Two lower panels: orbital motion of the planets with respect to the star, also projected on an axis inclined by 31.8 deg to the east of the north.

In the text