The mass of Beta Pictoris c from Beta Pictoris b orbital motion

We aim to demonstrate that the presence and mass of an exoplanet can now be effectively derived from the astrometry of another exoplanet. We combined previous astrometry of $\beta$ Pictoris b with a new set of observations from the GRAVITY interferometer. The orbital motion of $\beta$ Pictoris b is fit using Markov chain Monte Carlo simulations in Jacobi coordinates. The inner planet, $\beta$ Pictoris c, was also reobserved at a separation of 96\,mas, confirming the previous orbital estimations. From the astrometry of planet b only, we can (i) detect the presence of $\beta$ Pictoris c and (ii) constrain its mass to $10.04^{+4.53}_{-3.10}\,M_{\rm Jup}$. If one adds the astrometry of $\beta$ Pictoris c, the mass is narrowed down to $9.15^{+1.08}_{-1.06}\,M_{\rm Jup}$. The inclusion of radial velocity measurements does not affect the orbital parameters significantly, but it does slightly decrease the mass estimate to $8.89^{+0.75}_{-0.75}\,M_{\rm Jup}$. With a semimajor axis of $2.68\pm0.02$\,au, a period of $1221\pm15$ days, and an eccentricity of $0.32\pm0.02$, the orbital parameters of $\beta$ Pictoris c are now constrained as precisely as those of $\beta$ Pictoris b. The orbital configuration is compatible with a high-order mean-motion resonance (7:1). The impact of the resonance on the planets' dynamics would then be negligible with respect to the secular perturbations, which might have played an important role in the eccentricity excitation of the outer planet.


Introduction
The formation and evolution of giant exoplanets is an intense field of research. Several formation scenarios, ranging from gravitational instabilities (e.g., Boss 2011; Nayakshin 2017) to a variety of core-accretion models (e.g., Alibert et al. 2005;Emsenhuber et al. 2020), are still actively debated. The issue is the lack of observables with which to distinguish between the different scenarios. One solution is to analyze the atmospheric composition and search for formation signatures (Öberg et al. 2011;Mordasini et al. 2016;Madhusudhan et al. 2017). Another solution consists in measuring the energy dissipation during formation as a function of the final exoplanetary mass (Mordasini 2013;Marleau & Cumming 2014).
This paper focuses on the latter. Obtaining a dynamical mass for young directly imaged exoplanets is difficult. Only a handful of these objects have published dynamical masses: β Pictoris b (Snellen & Brown 2018;Dupuy et al. 2019;Lagrange et al. 2020;Vandal et al. 2020), β Pictoris c , PDS 51 Pegasi b Fellow. 70 b (Wang et al. 2021), and HR8799 e (Brandt et al. 2021b). One of the main difficulties of a direct measurement is the long orbital period of directly imaged exoplanets. Another is the fact that young stars are often pulsating (for example, β Pictoris is a δ Scuti variable), which makes accurate radial velocity (RV) measurements difficult.
An efficient technique for measuring dynamical masses is observing multi-planetary systems and detecting the mutual influence of planets. The main historical example is the prediction of Neptune from the irregularities in the orbit of Uranus by Le Verrier (1846). A more recent example is the first exoplanetary system detected around PSR B1257+12 (Wolszczan & Frail 1992), where the mutual interactions of the planets were used to confirm their masses (e.g., Konacki & Wolszczan 2003). Last, transit timing variations (Agol et al. 2005) has become a key technique for obtaining the mass of transiting planets (e.g., Demory et al. 2020). In this paper we demonstrate that optical interferometry is now able to detect and measure the mass of an exoplanet solely from the astrometry of another exoplanet at larger separation. We focus on the β Pictoris system, where both the outer planet (b; Lagrange et al. 2010) and the inner planet (c; Lagrange et al. 2019) have well-characterized orbital parameters Nowak et al. 2020). We use the GRAVITY instrument (Gravity Collaboration et al. 2017), a near-infrared interferometer operating at the Very Large Telescope (VLT) at Cerro Paranal. The interferometer has been designed to theoretically reach 10 µas accuracy  and has effectively demonstrated this level of accuracy on SgrA* flares at the center of the Milky Way (Gravity Collaboration et al. 2018). On exoplanets, it has demonstrated a typical 50 µas accuracy (Gravity Collaboration et al. 2019;Lagrange et al. 2020).
In Sect. 2 we present the new data reduction and additional, recent astrometric observations. Section 3 presents the restricted three-body model that we use for the Markov chain Monte Carlo (MCMC) fits. In Sect. 4 β Pictoris c is detected from β Pictoris b astrometry only. In Sect. 5 we perform a fit that includes β Pictoris c astrometry and RVs. In Sect. 6 we discuss dynamical insights from the orbital solutions. Section 7 is the conclusion.

Observations
We obtained three additional observations of the β Pictoris system with GRAVITY. Beta Pictoris c was observed during the night of 6 January 2021 and planet b during the nights of 7 January and 27 August 2021. The weather conditions ranged from very good to bellow average (on 27 January). The log of the observations are presented in Table A.1. These observations were part of the ExoGRAVITY large program , taken with a similar observation sequence as for past exoplanet detections (Gravity Collaboration et al. 2019: The fringe tracker  observes the star while the science camera (Straubmeier et al. 2014) alternately observes the star and the exoplanet.
The new observations of β Pictoris b and c were processed with the Public Release 1.5.0 (1 July 2021 1 ) of the ESO GRAV-ITY pipeline (Lapeyrere et al. 2014). This new version has slightly better astrometric capabilities by accounting for the effect of differential astigmatism (GRAVITY Collaboration et al. 2021). We also reprocessed previously published GRAVITY data sets with the same pipeline version. From the seven epochs published by Lagrange et al. (2020), we discarded two data sets: one from 2 November 2019 and one taken on 7 January 2020. In both circumstances, the very low coherence time (τ 0 < 1.5 ms) and limited number of exposures made the calculation of astrometric uncertainties difficult. The updated and new astrometric values are presented in Table 1.

A restricted three-body problem for multi-planetary systems
The Newtonian three-body problem is often addressed in a restricted version, where the mass of one body is negligible with respect to the other two (e.g., a satellite within the gravitational field of Earth and the Sun). Our case is different because the masses of the exoplanets are similar. However, a different restricted version of the three-body problem can be derived for systems with comparable masses. It uses Jacobi coordinates (Plummer 1918). Beust (2003) already developed the formalism for any Nbody system consisting of hierarchically nested orbits. Here we 1 https://www.eso.org/sci/software/pipelines/gravity/ applied it to the specific case of a solar-mass object with two planets. The Lagrangian writes L = T − V, where T is the kinetic energy and V the potential energy. With three bodies interacting by gravitational force, this writes (G is the constant of gravitation): where m and x correspond to the mass and position of the star, and m b , x b , m c , and x c , the mass and position of exoplanets b and c, respectively. To transform these equations into Jacobi coordinates, we set, according to Wisdom & Holman (1991): where q represents the position of planet c relative to the star, Q the position of planet b relative to the center of mass of the system (star, planet c), and R the position of the center of mass of the system. We also set M = m + m b + m c and ν = m + m c .
The vectors x ,b,c can be rewritten in terms of these Jacobi coordinates: In the barycentric reference frame, where the velocity of the center of massṘ is zero, we have: The terms of Eq. (1) are replaced with these new expressions: In the restricted case that matters to us (|x b −x | >> |x c −x |, i.e. |Q| >> |q|), the potential energy becomes, to first order: and the Lagrangian can be approximated by: In this expression, the two Jacobi variables Q and q are decoupled, and the Lagrangian can be written as the sum of two terms, L = L q + L Q , with: The term L q corresponds to the Lagrangian of the classical two-body problem that describes the orbit of planet c around the star, whereas L Q is the Lagrangian of the two-body problem corresponding to planet b orbiting around a virtual particle of mass m c + m located at the center of mass of the system (star, planet c). Both quantities can be solved analytically. This is how we modeled the orbital motion of β Pictoris b and c. The model should be the same as that used in Brandt et al. (2021a). We validate that our orbital model is sufficiently accurate for our GRAVITY astrometry in Appendix B.

Detection of β Pic c "with the point of [a] pen"
Here, instead of a pen 2 , we use the orbitize! code 3 . In this section we fit the relative astrometry of β Pictoris b only to assess whether we can indirectly detect β Pictoris c. We fit both a one-planet model and a two-planet model to only the astrometry of β Pictoris b. The one-planet model is a repeat of the fit done in Gravity Collaboration et al. (2020), using the same eight orbital parameters: semimajor axis (a b ), eccentricity (e b ), inclination (i b ), argument of periastron (ω b ), position angle of the ascending node (Ω b ), epoch of periastron in fractional orbital periods after MJD 59,000 (τ b ), system parallax, and total mass (M tot ). We used all the same priors as Gravity Collaboration et al. (2020). The two-planet model fit adds orbital elements for a second planet (a c , e c , i c , ω c , Ω c , τ c ) as well as replacing total mass with component masses (M * for the star and M b and M c for planets b and c, respectively). The priors on most of the orbital elements of β Pictoris c are the same as for b, except for a log uniform prior from 0.1 to 9 au for a c . We used the same prior on M tot as M * . We used a log uniform prior of 1 to 50 M Jup for M c . We fixed the mass of M b to 10 M Jup since our orbital model described in Sect. 3 cannot particularly constrain M b unless M * is known to < 1% precision.  The one-planet and two-planet model fits that use only the b astrometry are shown in red and cyan, respectively. These models are described in Sect. 4. The two-planet model fit that uses b and c astrometry is also plotted, in blue, and the fit that also uses the RV is plotted in purple (Sect. 5). The top row shows the orbit models and all of the data. The middle row shows the residuals after subtracting a pure Keplerian orbit for planet b based on the orbital parameters from the two-planet model using b and c astrometry. The bottom row shows the residual after accounting for the perturbation of planet c. We note that, although the red one-planet model is a pure Keplerian orbit, it is not a flat line in the middle row because the best-fit one-planet Keplerian model also attempts to fit the perturbations due to the second planet.
In both cases we used the parallel-tempered affine-invariant sampler in ptemcee (Foreman-Mackey et al. 2013;Vousden et al. 2016), using 20 temperatures and 1000 walkers per temperature. We obtained 30,000 samples of the posterior per walker after a "burn-in" of 10,000 steps for each walker in the one-planet model fit. In the two-planet model fit, we obtained 5,000 samples of the posterior from each walker after a burn-in of 55,000 steps for each walker. The posteriors for the parameters are given in Table 2. For the one-planet fit, there are no assumptions -and therefore no constraints -on planet c. The two-planet fit is able to indirectly measure a distinct mass and a c for the second planet.
To assess whether adding a second planet significantly improves the fit to the data, we computed the Bayesian information criterion (BIC) of the maximum likelihood model for both models. The one-planet model gives a BIC of 1247, while the two-planet model fit gives a BIC of 1784. The difference of 537 in the BIC indicates definitively that we have indirectly detected β Pictoris c using only the relative astrometry of β Pictoris b. For comparison, BIC changes between 10 and 100 have been used to show significant detections of outer planets in RV data (Christiansen et al. 2017;Bryan et al. 2019). This is also reflected in the residuals to the orbit fits shown in Fig. 1  The top row shows the data (black error bars) and models (same coloring as Fig. 1). The bottom row shows the same data but after subtraction of the median orbit of the two-planet model fit from b and c astrometry.
The two-planet model that uses only b astrometry (in teal) only weakly constraints the angular separation of planet c, and a direct detection of planet c is needed to obtain precise orbital constraints.
the GPI/SPHERE astrometry and the residual to the GRAVITY astrometry, the one-planet model fit in red fails to fully fit all of the data points and is systematically off. This corroborates the large change in BIC going from the one-planet to the two-planet model.
In addition to the detection of β Pictoris c, we also determined its mass, obtaining 10.04 +4.53 −3.10 M Jup . The estimate is slightly high, although it is still consistent with previous mass derivations (Lagrange et al. 2019;Vandal et al. 2020;Brandt et al. 2021a). The eccentricity was also determined, but with large uncertainties: e c = 0.90 +0.08 −0.26 . The value is high, although the 95% credible interval (in Table 2) is still consistent with previously published estimations. However, we have not included all of the data, and this merely demonstrates the ability of GRAV-ITY to indirectly detect a second planet in the system. The direct detection of β Pictoris c brings much better orbital constraints, as shown in the following section.

Refined masses and orbits of β Pic b and c
More accurate orbital parameters of β Pictoris c can be derived if we include the relative astrometry of c. This, in turn, gives a better mass estimation of c. The two-planet model fit was repeated, but with the addition of the β Pictoris c GRAVITY astrometry. The same MCMC walker parameters were used to obtain the 500,000 samples of the posterior. The posteriors are plotted in blue in both Fig. 1 (β Pictoris b) and Fig. 2 (β Pictoris c). The new β Pictoris c observation now accurately constrains the eccentricity of the orbit (e = 0.32 +0.03 −0.02 ) and the mass of β Pictoris c (9.15 +1.08 −1.06 M Jup ). The mass of β Pictoris c is consistent and comparable in precision to the mass obtained from stellar RVs and absolute astrometry (Lagrange et al. 2019;Vandal et al. 2020;Brandt et al. 2021a). This improves the robustness of the mass estimate since GRAVITY astrometry is not affected by systematic errors in stellar RVs on young stars or in absolute astrometry on bright stars. Last, we performed a two-planet fit that included the RV data of the star from Lagrange et al. (2020) that was reprocessed by Vandal et al. (2020). The inclusion of the RV data allowed us to fit for the mass of β Pictoris b, which we had previously left fixed. We used a log uniform prior between 1 and 50 M Jup for its mass prior. For the MCMC analysis, we used two temperatures and 1000 walkers per temperature. Each walker underwent a 10,000-step burn-in before 1,000 samples of the posterior were obtained from each walker. We find a semimajor axis of 2.68 ± 0.02 au, corresponding to an orbital period of 1221 ± 15 days. The mass estimate of β Pictoris c is still within the previous fit uncertainties, at 8.89 +0.75 −0.75 M Jup , and we estimate a mass of 11.90 +2.93 −3.04 M Jup for β Pictoris b. Both component masses agree well with the latest dynamical mass estimates (Vandal et al. 2020;Brandt et al. 2021a). The two masses are also in very good agreement with a hot-core accretion scenario (Mordasini 2013), as proposed by Nowak et al. (2020).

Stability, resonance, and evaporating bodies
The β Pictoris system has only recently joined the selective group of multi-planet directly imaged systems, along with HR 8799 (Marois et al. 2008(Marois et al. , 2010, PDS 70 Müller et al. 2018;Haffert et al. 2019), and TYC 8998-760-1 (Bohn et al. 2020). Because of the lack of planet detection and the struggle to constrain long-period orbits, the formation and dynamical evolution of cold giant planets remain largely unconstrained.
The eccentricity is thought to trace the formation and dynamical evolution of wide substellar objects. Indeed, wide-orbit giant planets and brown dwarfs appear to have different eccentricity distributions (Bowler et al. 2020). Low eccentricities and nearly planar planetary configurations are nevertheless usually associated with stable systems that formed via core accretion, such as the Solar System, even if post-formation giant impacts, planetplanet interactions, or scattering can excite orbital eccentricities afterward. With their orbits currently separated by ∼ 5 relative Hill radii, the current configuration of the two planets is likely stable, a statement confirmed by Beust et al. (2021, in preparation). However, both planets presumably induce significant secular perturbations on each other due to their high masses. The eccentricity variations caused by these perturbations are detailed in Appendix C. It is interesting to note that, in most solutions from the orbital fit, planet b periodically reduces its eccentricity to a negligible value, while the eccentricity of planet c remains greater than 0.2. This suggests that, contrary to planet c, the eccentricity of planet b may not be primordial, but rather the result of secular perturbations from planet c. Secular interaction also applies to inclination. It has long been noted that the midplane of the inner disk differs from that of the outer (Mouillet et al. 1997;Heap et al. 2000;Augereau et al. 2001) by a few degrees. Secular perturbations arising from both planets themselves having inclination oscillations could indeed cause this phenomenon.
Both HR 8799 and PDS 70 planets are suspected to be in a configuration of mean-motion resonance (MMR; respectively 1 : 2 : 4 : 8 and 2 : 1; Wang et al. 2018Wang et al. , 2021. Although the orbital fits are not yet precise enough to confirm this, commensurable periods are compatible with the astrometric constraints and increase the stability of the systems. In the β Pictoris system, however, the separation between the planets forbids any low-order MMR.  Table 2. Posteriors for the main parameters of the orbital fits. Orbital parameters are in Jacobi coordinates and follow the definitions in Blunt et al. (2020), with the exception of the reference epoch for τ as MJD 59,000 (31 May 2020). For each orbital parameter, the median is reported along with the 68% credible interval in super-and subscript (the 95% credible interval is reported in parentheses).
the high masses of the planets that will introduce non-negligible short-scale variations to the orbital elements. The existence and characteristics of β Pictoris b had notably been postulated long before the first direct-imaging detection as a consequence of the observation of numerous falling evaporating bodies (FEBs) close to the star (Beust & Morbidelli 1996. Those star-grazing planetesimals were thought to come from an asteroid belt in 4 : 1 MMR with a slightly eccentric giant planet with predicted orbital characteristics that approximately match those of β Pictoris b (Thébault & Beust 2001;Beust & Valiron 2007). The presence of another massive planet such as β Pictoris c orbiting well inside β Pictoris b can perturb this picture. As a matter of fact, FEB progenitors trapped in 4:1 MMR with β Pictoris b (i.e., orbiting at ∼4 au from the star) undergo a gradual eccentricity increase with only little changes in semimajor axis before reaching the star-grazing state. This way they inevitably come to cross the orbit of β Pictoris c twice per orbit. Thus, close encounters with that planet would presumably eject most of the comets before they can become observable FEBs. Hopefully, the refined orbital and mass constraints provided in the present paper will help explore this rich dynamical issue. Beust et al. (2021, in preparation) show that the presence of β Pictoris c entirely prevents FEB progenitors from reaching the FEB state. Depending on its eccentricity, their source 4:1 MMR at 4 au may not even be stable. However, various inner MMRs with β Pictoris c (instead of b) could be a valuable alternate source of FEBs, leading to a revision of the FEB scenario.

Conclusions
Our results demonstrate the capabilities of relative astrometry with optical interferometry to characterize multi-planetary systems: 1. We were able to detect β Pictoris c and measure its mass from the astrometry of β Pictoris b alone. Our mass measurement of 10.04 +4.53 −3.10 M Jup is in agreement with the previously published estimations based on RV by Lagrange et al. (2020) and Vandal et al. (2020). This is the first time that the presence of an exoplanet has been derived from the astrometry of another exoplanet. 2. The new detection of β Pictoris c, at 96 mas from its star, is a record in terms of exoplanet detection at small separation. Thanks to this new measurement, the orbital parameters of c, especially its eccentricity at 0.32±0.02, are well determined. Fitting all data sets together, we obtain 8.89 ± 0.75 M Jup for β Pictoris c and 11.90 +2.93