Direct confirmation of the radial-velocity planet $\beta$ Pic c

Methods used to detect giant exoplanets can be broadly divided into two categories: indirect and direct. Indirect methods are more sensitive to planets with a small orbital period, whereas direct detection is more sensitive to planets orbiting at a large distance from their host star. %, and thus on long orbital period. This dichotomy makes it difficult to combine the two techniques on a single target at once. Simultaneous measurements made by direct and indirect techniques offer the possibility of determining the mass and luminosity of planets and a method of testing formation models. Here, we aim to show how long-baseline interferometric observations guided by radial-velocity can be used in such a way. We observed the recently-discovered giant planet $\beta$ Pictoris c with GRAVITY, mounted on the Very Large Telescope Interferometer (VLTI). This study constitutes the first direct confirmation of a planet discovered through radial velocity. We find that the planet has a temperature of $T = 1250\pm50$\,K and a dynamical mass of $M = 8.2\pm0.8\,M_{\rm Jup}$. At $18.5\pm2.5$\,Myr, this puts $\beta$ Pic c close to a 'hot start' track, which is usually associated with formation via disk instability. Conversely, the planet orbits at a distance of 2.7\,au, which is too close for disk instability to occur. The low apparent magnitude ($M_{\rm K} = 14.3 \pm 0.1$) favours a core accretion scenario. We suggest that this apparent contradiction is a sign of hot core accretion, for example, due to the mass of the planetary core or the existence of a high-temperature accretion shock during formation.


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

Date
UT Time Nexp / NDIT / DIT airmass tau 0 seeing 2020-02-10 02:32:52 -04:01:17 11 / 32 / 10 s 1.16-1.36 6-18 ms 0.5-0.9" 2020-02-12 00:55:05 -02:05:29 11 / 32 / 10 s 1.12-1.15 12-23 ms 0.4-0.7" 2020-03-08 00:15:44 -01:41:47 12 / 32 / 10 s 1.13-1.28 6-12 ms 0.5-0.9"  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, m sin(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 highcontrast 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 m sin(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 et al. 2019), and atmospheric properties (Mollière et al. 2020) of imaged extrasolar planets with short periods. In this paper, 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.

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 et al. 2017). The first detection was obtained from the allocated time of the Ex-oGRAVITY 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).

Data reduction and detection
The data reduction used for β Pic c follows what has been developed for β Pic b (GRAVITY Collaboration et al. 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

MJD
∆RA  Table 2. Relative astrometry of β Pic c extracted from our VLTI/GRAVITY observations. Due to the interferometric nature of the observations, a correlation coefficient ρ is required to properly describe the confidence intervals, which are not aligned on the sky coordinates. The covariance matrix can be reconstructed using σ∆RA 2 and σ∆DEC 2 on the diagonal, and ρσ∆RAσ∆DEC off-diagonal.
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 Figure 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.
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 Kband 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 Figure 2.

The orbit of β Pic c
The relative astrometry is plotted in Figure 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! 1 (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 & al. (in preparation) 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 et al. 2018). Lastly, for the radial velocity, we used the 12 years of HARPS measurements (Lagrange et al. 2019;Lagrange & al. in preparation).
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 ( 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 • Table 3. Posteriors of the MCMC analysis of the β Pictoris system. Orbital elements are given in Jacobi coordinates. The ascending node is defined to be the node where the motion of the companion is directed away from the Sun. The epoch of periastron is given by orbitize! as a fraction of the period from a given date (here, 58889 MJD). The difference between the β Pic b mass from the prior and posterior indicates a tension between the dynamical data and spectral analysis. The periods are calculated from the posteriors. The K band delta magnitudes are obtained from present and published (GRAVITY Collaboration et al. 2020) GRAVITY spectra.
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 & al. in preparation). 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 ≈ 1 200 µ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 GRAV-ITY 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 & al. (in preparation). 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) 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.

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 fluxcalibrated, the overall flux level also constrains the radius of the planet.  Table 3.
The best fits obtained with Exo-REM and Drift-Phoenix are overplotted to the GRAVITY spectrum in Figure 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 & 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.

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 Marleau & Cumming 2014), or if the shock dissipates sufficient heat into the in-falling gas (Marleau et al. , 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 Figure 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 fluxcalibrated planet spectrum and W as the associated covariance matrix, the total flux in the ESO K-band is: F = T T S. The error bar is given by: ∆F = T W T 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 Kband 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 Figure 4, together with the 18.5 ± 2.5 Myr AMES-Cond isochrones, and two synthetic populations of planets formed through core nucleated accretion 3 . 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 Hot-start AMES-Cond (18.5 ± 2.5 Myr) Core-accretion (nominal) Core-accretion (hot accretion) Pictoris b Pictoris c 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).
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 et al. 2020), and we over-plotted the result in Figure 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 et al. 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 paper 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 coremass effect is taken into account. Interestingly, a similar conclusion has been drawn for β Pic b (GRAVITY Collaboration et al. 2020) based primarily on the estimation of its C/O abundance ratio.

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 veloc-ities. 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. 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 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 Figure 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 & al. in preparation). 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).

Appendix B: Data reduction
The theoretical planetary visibility is given as follows: 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 et al. 2020): 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: 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: 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, Equation (B.3) should be written for t * = t. This means that the factor G from Equation (B.2) and (B.3) do not cancel out perfectly, and an extra factor G(t)/G(t * ) should appear in Equation 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 Equation 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 Equation B.4. The χ 2 sum is then written: with W the covariance matrix affecting the projected visi-bilitiesṼ. This χ 2 , 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 Equation B.4 is flat) and can be taken as a reference: χ 2 no planet = χ 2 planet (0, 0). The values of the χ 2 using a planet model can be compared to this reference, by defining: z(∆α, ∆δ) = χ 2 no planet − χ 2 planet (∆α, ∆δ), (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 Figure 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, Equation (B.4) can be used to extract the contrast spectrum C (GRAVITY Collaboration et al. 2020). The extraction of the astrometry usingṼ model 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 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.
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 REBOUND N-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.  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 (dotted line). The two lower panels show the 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.