Unveiling the β Pictoris system, coupling high contrast imaging, interferometric, and radial velocity data

Context. The nearby and young β Pictorissystem hosts a well resolved disk, a directly imaged massive giant planet orbiting at (cid:39) 9 au, as well as an inner planet orbiting at (cid:39) 2.7 au, which was recently detected through radial velocity (RV). As such, it offers several unique opportunities for detailed studies of planetary system formation and early evolution. Aims. We aim to further constrain the orbital and physical properties of β Pictorisb and c using a combination of high contrast imaging, long base-line interferometry, and RV data. We also predict the closest approaches or the transit times of both planets, and we constrain the presence of additional planets in the system. Methods. We obtained six additional epochs of SPHERE data, six additional epochs of GRAVITY data, and ﬁve additional epochs of RV data. We combined these various types of data in a single Markov-chain Monte Carlo analysis to constrain the orbital parameters and masses of the two planets simultaneously. The analysis takes into account the gravitational inﬂuence of both planets on the star and hence their relative astrometry. Secondly, we used the RV and high contrast imaging data to derive the probabilities of presence of additional planets throughout the disk, and we tested the impact of absolute astrometry. Results. The orbital properties of both planets are constrained with a semi-major axis of 9.8 ± 0.4 au and 2.7 ± 0.02 au for b and c, respectively, and eccentricities of 0.09 ± 0.1 and 0.27 ± 0.07, assuming the H IPPARCOS distance. We note that despite these low ﬁtting error bars, the eccentricity of β Pictorisc might still be over-estimated. If no prior is provided on the mass of β Pictorisb, we obtain a very low value that is inconsistent with what is derived from brightness-mass models. When we set an evolutionary model motivated prior to the mass of β Pictorisb, we ﬁnd a solution in the 10–11 M Jup range. Conversely, β Pictorisc’s mass is well constrained, at 7.8 ± 0.4 M Jup , assuming both planets are on coplanar orbits. These values depend on the assumptions on the distance of the β Pictorissystem. The absolute astrometry H IPPARCOS - Gaia data are consistent with the solutions presented here at the 2 σ level, but these solutions are fully driven by the relative astrometry plus RV data. Finally, we derive unprecedented limits on the presence of additional planets in the disk. We can now exclude the presence of planets that are more massive than about 2.5 M Jup closer than 3 au, and more massive than 3.5 M Jup between 3 and 7.5 au. Beyond 7.5 au, we exclude the presence of planets that are more massive than 1–2 M Jup . Conclusions. Combining relative astrometry and RVs allows one to precisely constrain the orbital parameters of both planets and to give lower limits to potential additional planets throughout the disk. The mass of β Pictorisc is also well constrained, while additional RV data with appropriate observing strategies are required to properly constrain the mass of β Pictorisb.


Introduction
The nearby (parallax = 51.44 ± 0.12 mas, van Leeuwen 2007) 1 β Pictoris system has been intensively studied since the mid 1980s, when a circumstellar debris disk of dust was first and the infalling comets were attributed to planets in the 1990s, it took more than a decade to detect the first planet in the system (Lagrange et al. 2010). It appeared that β Pictoris b could indeed explain both the disk inner warp (Augereau et al. 2001) and the infalling comets, provided it had a nonzero eccentricity (Beust & Morbidelli 1996, 2000. Many high contrast imaging (hereafter HCI) observations have constrained its orbital properties, and, in particular, excluded highly eccentric (e ≥ 0.25) orbits (Lagrange et al. 2019a). Using just one observation with the GRAVITY long baseline interferometer of the position of β Pictoris b relative to the star combined with previously published GPI (Macintosh et al. 2008) HCI data, GRAVITY Collaboration (2020) found an eccentricity of 0.11-0.20, while Nielsen et al. (2020), using newly re-calibrated GPI data found an eccentricity in the range of 0.09-0.16. However, the lower limits on β Pictoris b's eccentricity remained quite elusive. Finally, using the K-STacker algorithm, which is an approach that is fully different from the Markov-chain Monte Carlo (MCMC) based methods, Le Coroller et al. (2020) derived an eccentricity of 0.07.
β Pictoris b is the only directly imaged exoplanet with a constrained dynamical mass. Radial velocity (RV) data showed that, assuming a zero eccentricity, the planet mass should not exceed 15 M Jup (Lagrange et al. 2012). A first order combined analysis of HCI and RV data (Bonnefoy et al. 2014) led to a similar result. This was also in agreement with brightness-mass relationships, which gave values between 9 and 13 M Jup (Lagrange et al. 2010). Using HIPPARCOS and Gaia data, Snellen & Brown (2018) derived a mass in the range of [10,12] M Jup . Also, based on the GPI relative astrometry and an analysis of HIPPARCOS-Gaia data, Kervella et al. (2019) find a mass of 13.7 +6 −5 M Jup . Nielsen et al. (2020), using mainly NACO, GPI relative astrometry up to November 2018 (using a new calibration of GPI data), as well as HIPPARCOS and Gaia absolute astrometry, find a mass of 12.8 +5.3 −3.8 M Jup . Taking into account HIPPARCOS and Gaia data also increases the semi major axis (sma) of β Pictoris b's orbit to 10.2 au, while they obtain an sma of 9.3 au when they do not take the HIPPARCOS and Gaia data into account. Using NACO, GPI prior to 2017, as well as long base line GRAVITY interferometric data, GRAVITY Collaboration (2020) find a mass of 12.7 ± 2.2 M Jup ; the range of these values is just compatible with the modeling of the Gravity K-band spectrum. All of these estimates that consider only one planet need to be revised as a recent analysis of ESO/HARPS high resolution spectroscopic data obtained since the beginning of the instrument (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) has led to the discovery of a second, inner planet (Lagrange et al. 2019a), with a mass of 9 times the mass of Jupiter, orbiting at 2.7 au on a slightly eccentric (e 0.2) orbit. Taking β Pictoris c into account, Dupuy et al. (2019) derived a mass for β Pictoris b in the range of [10,16] M Jup , while Nielsen et al. (2020) find a mass of 8.35 M Jup , which is much lower than that found without the RV data. Nonetheless, the use of the absolute astrometric data remains tricky at this stage because of the still limited information provided by Gaia DR2 and the fact that this star is saturated in the observations of Gaia. This aspect is addressed in detail in Nowak et al. (2020).
Improving the planets orbit and mass characterization requires both very precise relative astrometry, hopefully on both planets, on a time basis as long as possible, as provided by GRAVITY (see GRAVITY Collaboration 2020, in the case of β Pictoris b), and long time series RV data. In this paper, we present new SPHERE, GRAVITY, and HARPS observations of β Pictoris . We use these new and previous data to combine, in a single MCMC analysis, HCI, RV, and GRAVITY data on β Pictoris b to further constrain both planets.
The system may also host additional planets. The presence of clumps seen at mid-IR (Telesco et al. 2005;Okamoto et al. 2004;Wahhaj et al. 2003) and asymmetries seen in sub-mm (ALMA) data (see Cataldi et al. 2018, and ref. there-in) as well as the presence of the innermost hole (Lagage & Pantin 1994) indicate that planets could be present in the disk. In Lagrange et al. (2019b), we estimated the probability of presence of additional planets from a fraction of an au up to hundreds of au, coupling HCI NACO data with RV data. Companions with masses larger than 15 M Jup were excluded throughout the disk except between 3 and 5 au with a 90% probability. We could, moreover, exclude the presence of planets that are more massive than 3 M Jup closer than 1 au and further than 10 au, with a 90% probability. In the 3-5 au region, we excluded companions with masses larger than 18 (resp. 30) M Jup with probabilities of 60 (resp. 90)%. The present set of data, combined with the results of the properties of β Pictoris b and c, allows a significant improvement in the mass upper limits of possible additional planets.
The various data are described in Sect. 2. The orbital analysis is performed in Sect. 3 and some consequences are discussed in Sect. 4. We combined all HCI (more than 30 epochs) and all RV data to revisit the probability of detection of additional planets in Sect. 5. The conclusion and perspectives are presented in Sect. 6.

SPHERE data
Six new high contrast, IRDIS (Dohlen et al. 2008;Vigan et al. 2010) and IFS (Claudi et al. 2008;Zurlo et al. 2014) coronagraphic SPHERE (Beuzit et al. 2019) data were obtained between October 2018 and February 2020. The first three datasets were taken as part of the SPHERE SHINE GTO (Desidera et al, in prep.), and the last three as part of a classical program (PI Lagrange; 0104.C-0418(D), 0104.C-0418(E), 0104.C-0418(F)). The October 2018 data immediately followed the first observations (in Sept. 2018) of β Pictoris b after its closest approach to the star described in Lagrange et al. (2019a). All but one of the new datasets presented here were obtained in the H2 (central wavelength = 1.59 µm; width = 0.055 µm) and H3 (central wavelength = 1.67 µm; width = 0.056 µm) IRDIS narrow bands, as well as in YJ (0.95â1.35 µm, spectral resolution R 54) with IFS. One data set (March 2019) was taken in the K1 (central wavelength = 2.1025 µm ; width = 0.102 µm) and K2 (central wavelength = 2.255 µm ; width = 0.109 µm) IRDIS narrow bands, and in YH (0.95â1.65 µm, R 33) with IFS. All data were taken with apodized Lyot coronagraphs, including either a 185 mas focal mask (N_ALC_Y JH_S ) or a smaller (145 mas) mask (N_ALC_Y J_S ), when the planet was closer to the star, as well as a pupil mask. Our new data were recorded in stabilized pupil mode so as to perform angular differential imaging (ADI) post-processing (Marois et al. 2006;Galicher et al. 2018), and with the waffle mode on. This mode allows one to apply voltages on the deformable mirror to generate four symmetrical spots in the form of a "cross" or a "plus" shape; these spots are used to measure the star position at the center of the pattern.
Each observing sequence followed the usual following steps: PSF-Coronographic observations-PSF-Sky. The PSF were used as relative photometric references and to obtain an estimate of the PSF just before and just after the observations. This PSF was used to model the planet signal during the characterization step. Sky data were recorded at the end of the coronagraphic sequence to get an estimate of the background level in the science A18, page 2 of 17 A. M. Lagrange et al.: Unveiling the β Pictoris system images. Finally, the pixel scales and true north were measured from IRDIS observations of astrometric fields (see Table 1). We note that in the March 2019 set of data, the IRDIS data were saturated very close to the mask, so only the IFS are considered hereafter for this set. The IRDIS and IFS data were reduced as described in Chauvin et al. (2017) and Delorme et al. (2017), and using the SPECAL tool developed for SPHERE (Galicher et al. 2018). For each epoch, all individual data frames were recentered using the waffles to create a single cube, which was then analyzed to remove the stellar halo and measure the astrometry of the planet.
In the IRDIS data, the β Pictoris b positions relative to the star and associated uncertainties of the PCA images were measured in two ways: (1) fitting a model of the planet to the detected source, and (2) injecting a negative planet into the initial data cube and processing it to minimize the local flux in the ADIprocessed data. These two approaches are described in Sect. 3 of Galicher et al. (2018). They both use the PSF recorded prior to and after the coronagraphic data.
In the IFS data, the positions of β Pictoris b were also measured in two ways: (1) with PCA, as described in Zurlo et al. (2016), and (2) using the PACO software (Flasseur et al. 2020). The astrometric uncertainties were estimated using the errors associated with these fits as well as star centering errors and the errors on the pixel scales and true north.
With these new data, relative astrometric values at 18 epochs are now available. Figure 1 shows the positions obtained with the various methods. The astrometric measurements derived from method 1 and method 2 for IRDIS are in good agreement. In the following, we consider the average of these two values. Table C.1 provides the measurements of β Pictoris b's position in the IRDIS data (average of methods 1 and 2) as well as the previous NACO astrometric data used in the present study. The astrometric measurements on IFS data derived from method 1 and method 2 significantly differ, and the values found with IFS method 1 are closer to those found on IRDIS than those found with IFS method 2. This certainly illustrates that the PCA and PACO methods are intrinsically different, and they do not share the same biases. In the following, we use the values of β Pictoris b's position found with method 1 in the IFS data. They are provided in Table C.1.

GRAVITY data
Six new pieces of high accuracy astrometric data were obtained with the GRAVITY instrument (GRAVITY Collaboration 2017) in the K band. The first astrometric measurement is an updated value of the astrometry published by GRAVITY Collaboration  Table 1. Three of the observations were performed in medium resolution (R = 500), two in high resolution (R = 4000); and during one night (7 Jan 2020), the observations were taken with both modes.
The observation sequence was the same as in GRAVITY Collaboration (2019): the fringe tracker ) uses the central star to correct the atmospheric piston, but also to serve as a reference for the astrometry. The science fiber feeding the spectrometer points, alternatively, between the position of the star and the position of the planet. While the star can be seen on the acquisition camera, the position of the planet was obtained from dedicated software based on GPI data, which is available to the community 2 , or by our estimations when using SPHERE data as input for the orbit estimation.
The reduction starts by using the official GRAVITY pipeline (Lapeyrere 2014) to reduce the images to uncalibrated visibilities. The complex visibilities are then all phased with respect to the fringe tracker phase. From this point on, we used our tools to fit the phase and the spectra of the coherent flux 3 . The interferometric baselines, which are the equivalent of the imaging plate-scale (Woillez & Lacour 2013), use the Earth orientation parameters (EOP) to find the position of the telescopes with respect to the sky. The results of the astrometric fits are presented in Table C.2. The accuracy, presented as σ ∆RA and σ ∆Dec , is on the order of 100 µas. However, the geometry of the interferometer has a profound impact on the error bars. Because the baselines are considerably longer in the north-east direction, that is, 130 m compared to 50 m, the error bars are consistently worse in the north-west direction. This characteristic is accounted for by using the correlation coefficient ρ. This parameter varies between −1 for an anticorrelation, to 1 for a total positive correlation, and zero for no correlation. The correlation coefficients are given in Table C.2, and they are typically between −0.2 and −0.8. They are used in the MCMC analysis.

HARPS data
More than 1900 new HARPS spectra were obtained at five different epochs since the β Pictoris c discovery paper (Lagrange et al. 2019a) under programs 098.C-0739(A) and 0104.C-0418(A) (PI Lagrange). At each epoch, long sequences of data (several hours per night over 1, 3, 3, 4, and 3 days, respectively) were recorded to allow for the proper extraction of the planet signal. Indeed, the star is a high frequency (about 23-30 min) and high amplitude (up to 2 km s −1 ) pulsator, and the stellar pulsations need to be subtracted from the signal. The RVs were measured as in our previous works, using our SAFIR tool (Galland et al. 2005) and using the same approach as in Lagrange et al. (2019a). Basically, for each spectrum, the RV is computed, in the Fourier plane, with respect to a reference spectrum made of an average of all available and suitable data.
To estimate the stellar pulsations, we used Method 2 described in detail in Lagrange et al. (2019a), which consists in fitting the data with one sine function plus an offset, computing the residuals, repeating the process on the residuals, and so on, nine times. The final offset represents the planets-induced RV plus a constant offset due to the fact that the RVs are relative to a reference spectrum. We refer the readers to Lagrange et al. (2019a) for detailed explanations. The offset values determined for the new epochs as well as previous ones are provided in Table C.3. We note that the offsets presented here, similar to most of those recorded since JDB-2 454 000 = 4000, have low uncertainties. This is because we have significantly increased the duration of the visits (several hours per night and for three or more consecutive nights). The fits are provided in Appendix A. Figure 2 shows the RV time series (7000 data points since 2003) before and after correction from the pulsation signals. Since the last data were published in Lagrange et al. (2019a), the RV have increased again due to β Pictoris c.
The periodograms associated with both the raw and pulsation data show a peak at about 1200 d, which is not present in the window function, due to β Pictoris c. We conclude that the additional data presented in this paper confirm the presence of β Pictoris c.

Properties of β Pictoris b and c
The aim of the present paper is to combine various types of observational data into a single analysis. We used a Markovchain Monte Carlo (MCMC) Bayesian analysis to derive the probabilistic distribution of orbital and physical solutions. The code is an upgrade of the one presented in Grandjean et al. (2019); it uses the affine invariant MCMC ensemble sampler emcee (Foreman-Mackey et al. 2013). The model includes the orbital parameters of β Pictoris b and c as well as the mass and distance of β Pictoris . Walker moves were selected using differential evolution (80% DEMove + 20% DESnookerMove); these move rules ensure that steps taken by walkers are automatically scaled according to the current parameter distributions.
We first used only the astrometry deduced from the NACO + SPHERE HCI data on β Pictoris b, then we added the GRAVITY astrometric measurements on β Pictoris b to determine its orbital parameters. Finally, we took the RV data into account to characterize both β Pictoris b and β Pictoris c. In all cases, the star mass was also estimated. The priors are provided in Table 2.

Analysis of HCI data only
As a first step, we used the NACO + SPHERE IRDIS and IFS PCA astrometric data as input to our MCMC fit. We considered both the average of method 1 and method 2 results for IRDIS data as well as the values derived from method 1 for the IFS data. This set is labeled as set NIRDIFS.   Unless otherwise specified, we adopted a distance to the star of 19.44 pc (HIPPARCOS distance) when dealing with astrometric data only, even though Gaia DR2 provided a more recent estimate. This is to allow for a comparison with other works; this is also because the Gaia DR2 parallax has a larger uncertainty than the HIPPARCOS one due to its saturation, and because the time baseline to model this three-component system is a factor of two lower than the HIPPARCOS one.
The results are displayed in Fig. 3. The sma, eccentricity, and mass distributions are not symmetrical and cannot be fitted by Gaussian profiles. We note this was also the case in previous studies based on SPHERE or GPI results. The sma distribution peaks at 9 au and has a small wing at larger values. We note that this wing is much shallower than in Lagrange et al. (2019a). The eccentricity peaks slightly below 0.01, with a wing toward larger values; this barely rejects a zero eccentricity. The distribution of inclinations can be fitted by a Gaussian profile, with a maximum at 88.94 deg and one sigma of 0.04 deg. The star mass distribution peaks at 1.77 M , with a wing toward lower values. The obtained orbital elements compare well with those described in Lagrange et al. (2019a), as well as with those derived from case 2 using NACO and newly recalibrated GPI data from Nielsen et al. (2020) (who used the HIPPARCOS distance). The mass of the star derived by Nielsen et al. (2020) (case 1) is also comparable. Their eccentricity does not show a peak, but a steady decrease starting from zero.
We now check how these posteriors depend on the various assumptions or inputs by: -Assuming a Gaia distance, in the range of [19.63,19.88] pc, would lead to an sma distribution peaking at 9.17 au and a star mass peaking at 1.83 M . -Removing the two data points closest to the star 4 obtained in Nov 2016 (proj. sep 130 mas) and Sept 2018 (proj. sep 140 mas), that is, at the last observation before the closest approach and the first observation after the closest approach does not change the results. -Considering the IRDIS (average of method 1 and method 2 results) and IFS (method 1) data separately gives posteriors that are compatible within the uncertainties with the NIRDIFS case above; yet, we notice that (1) the median sma and median eccentricity are slightly shifted toward higher values when only using IRDIS data, as opposed to when only using IFS data, and the star mass is slightly lower, and that (2) the posteriors obtained when only using IFS data are closer to those obtained with both IRDIS and IFS (NIRDIFS) than the posteriors obtained using only IRDIS. This is because the IFS measurements are associated with smaller error bars than for the IRDIS measurements. -Finally, using the PACO measurements or the measurements obtained with method 1 as inputs for the IFS data only leads to compatible results, but with a slightly higher median eccentricity, 0.035. In conclusion, the data (IFS versus IRDIS) and the methods used to measure the relative astrometry of β Pictoris b slightly impact the posteriors, but the results are still compatible. Due to the smaller error bars, the IFS measurements (method 1) used in the NIRDIFS case dominate the astrometric information over the IRDIS data. The assumption on the distance to the star also impacts the orbital parameters and the star mass.

Analysis of HCI and GRAVITY data
We ran new MCMCs taking into account the relative astrometric data provided by GRAVITY on β Pictoris b in addition to the IRDIFS data. This new data set is hereafter referred to as NIRDIFS-GRAV. The obtained sma, eccentricity, and mass distributions are more symmetrical than when only using the NIRDIFS measurements, as seen in Fig. 4. For the sma, we derived a value of 9.64 ± 0.04 au, for the eccentricity, a value of 0.07 ± 0.01, and for the inclination, a value of 88.99 ± 0.01 deg. The star mass is 1.77 ± 0.03 M . Similar values are obtained without taking into account the Nov. 2016 and Sept. 2018 data. When considering the IRDIS and IFS (method 1) measurements separately, or when using the PACO measurements instead of the IFS method 1 ones, we obtain the same tendencies as in the previous section, with even smaller differences. The posteriors are now clearly driven by the GRAVITY data points 5 Finally, we then assumed the Gaia distance instead of the HIPPARCOS one. We got a higher sma (+0.15 au), a similar eccentricity, and a higher star mass (+0.09 M ).
The posteriors found here with NIRDIFS-GRAV are in qualitative-only agreement with those found by GRAVITY Collaboration (2020) from the relative astrometry: sma = 10.6 ± 0.5 au, e = 0.15 ± 0.4, and M = 1.82 ± 0.3 M . This can be explained by a combination of various differences: (1) different HCI data were used (GPI instead of SPHERE), (2) the GPI data suffered from calibration errors (see Nielsen et al. 2020 for the revised values), (3) we use, in the present paper, six 5 As an exercise, we performed the MCMC analysis using GRAVITY data alone as input and using the same priors as before; it turns out that even though they are much less precise (larger uncertainties), the posteriors are compatible with those found in our NIRDIFS-GRAV case, where we used more than 12 yr of HCI data. This illustrates the impact of the GRAVITY data, thanks to their exquisite precision. It has to be noted, though, that GRAVITY needs the HCI data to point and place the fiber at the rough position of the planet before any measurements can be carried out. Both data are then complementary. pieces of GRAVITY data instead of one, and (4) the value of the GRAVITY astrometry used in GRAVITY Collaboration (2020) has been slightly revised for the present paper.
In conclusion, adding the GRAVITY data to the NACO and SPHERE data leads to significantly different orbital parameters and, specifically, higher values of sma and eccentricity for β Pictoris b. The star mass is between that proposed by Crifo et al. (1997), based on photometry and the HIPPARCOS distance and that proposed by Zwintz et al. (2019) based on an asteroseismology study. Using GRAVITY also significantly decreases the fitting uncertainties associated to the orbital parameters, thanks to small error bars. GRAVITY, therefore, dominates the relative astrometric information. As a consequence, the results are also dominated by any biases in GRAVITY. The posteriors are also biased by the assumption of a single planet, and we see below that GRAVITY astrometry is sensitive to the influence of β Pictoris c on the relative astrometry of β Pictoris b. Finally, we note that the posteriors still depend on the assumption made on the distance of the system.

Analysis of HCI, GRAVITY, and RV data
We considered the RV data in addition to the relative astrometric data and constrained the properties of β Pictoris b and c. We used the 32 offsets as input to the MCMC analysis, because using thousands of RV data points in our MCMC analysis would require unrealistically large computing times. In principle, adding the RV data allows one to constrain the masses of both planets in addition to their orbital parameters. We, therefore, added these two additional parameters in our MCMC analysis, and added an RV offset to account for the fact that our RV measurements are not absolute velocities, but they are given with respect to a reference spectrum that is built by averaging a large amount of spectra (see Lagrange et al. 2019a).
As we do not have constraints on the inclination of β Pictoris c, we assumed that its orbit is coplanar with that of β Pictoris b. The distance of the system is considered to be uniform in the [19.39,19.48] pc range, and the mass is uniform within the [1.6,2] M range in the priors. As the GRAVITY precision is typically 0.1 mas, and the total astrometric perturbation introduced by β Pictoris c is about 1.5 mas over 3 yr and that introduced by β Pictoris b is about 4 mas over 20 yr, we took these perturbations into account in our MCMC computations, adding up the barycentric corrections due to each planet.
We first ran an MCMC analysis, using a uniform distribution for the masses of Yet, the posterior on the mass of β Pictoris b is 3.2 M Jup , which is not compatible with any luminosity-mass relation or any previous related estimate. Additionally, this is not compatible with the fact that we did not detect β Pictoris c in single direct images when at quadrature, as expected at about 150 mas, while we did, in fact, detect β Pictoris b at smaller projected separations 6 . With our current knowledge, we attribute this poor constraint on β Pictoris b's mass to the fact that (1) the RV time series does not fully cover β Pictoris b's period as seen in the third row of Fig. 2 7 ) and (2) the data have non-negligible noise. We used a prior of a Gaussian, centered on 14 M Jup , and a one-sigma of 1 M Jup (noted hereafter N(14,1)) for β Pictoris b's mass. This range is representative of the high masses estimated for β Pictoris b in the literature. The results are summarized in Table C.4 and visualized in Figs. 5 and 6. For β Pictoris b, we 6 We note that for this last argument, we assume that both planets formed in the same way. 7 In this figure, the RV induced by the planets correspond to the solutions obtained in the case where the mass of β Pictoris b is constrained in a Gaussian range (14,1), as described below. In the present case, the mass found for β Pictoris b is 3 M Jup , and its impact on the RV would be even smaller.
find an sma of 9.76 ± 0.04 au, an eccentricity of 0.086 ± 0.006, and an inclination of 89.01 ± 0.01 deg. For β Pictoris c, we find an sma of 2.69 ± 0.02 au, and an eccentricity of 0.27 ± 0.07. The star mass is 1.78 ± 0.03 M and the distance is 19.44 ± 0.02 pc. The orbital parameters of β Pictoris c are in good agreement with those found in Lagrange et al. (2019a), where we analyzed only the RV data and made some assumptions about β Pictoris b to correct the RV data from this planet contribution. The masses found for planets b and c are 11.1 ± 0.8 and 7.8 ± 0.4 M Jup , respectively. We see, in Fig. 2 (bottom) that the time baseline is now long enough to possibly see the impact of β Pictoris b on the RV through a modulation of the β Pictoris c induced RV curve.
To investigate the impact of the prior on β Pictoris b's mass, we performed a similar MCMC analysis with a Gaussian N(12,1) M Jup prior. We ended up with a mass of 9.7 ± 0.7 M Jup for β Pictoris b, which is significantly lower than when using the (14,1) M Jup prior and an unchanged mass for β Pictoris c. This confirms that β Pictoris b's mass is not fully constrained yet, while the mass of β Pictoris c is.
To investigate the impact of the RV precision on the results, we set a threshold of 15 m s −1 to the RV uncertainties, that is, all uncertainties lower than 15 m s −1 were set at 15 m s −1 . The  These last tests illustrate how critically the RVs constrain the mass of the planets. Another illustration is provided in Fig. 6 where we show both the RV curves corresponding to our fits of the NIRDIFS-GRAV-RV data, when using a N(14,1) prior (left) and when using a flat prior (right) on β Pictoris b's mass. They give posteriors of 11 and 3.2 M Jup , respectively, for β Pictoris b mass. We see that the RV curve corresponding to the lower mass for β Pictoris b fits the offsets in a slightly better way than the other one, but the difference is small. Longer time baselines for the RV, a better coverage of the RV at both β Pictoris c quadratures, and/or a better precision on the offsets (through a better removal of the stellar pulsations) will help further constrain the mass of β Pictoris b. Gaia full data (IADs) will also help constrain the planet masses.
To investigate the impact of the star distance, we ran an MCMC analysis, assuming the Gaia distance range of [19.63,19.88] pc. Compared to the posteriors obtained using the HIPPARCOS distance, the sma of β Pictoris b and c are slightly increased by, respectively, +0.16 au, which has yet to be compared to the uncertainties of 0.04 au, and +0,04 au, which has yet to be compared to the uncertainties of 0,02 au. The stellar mass is now 1.87 ± 0.04 M , compared to 1.77 ± 0.03 M . The masses of β Pictoris b and c are only marginally increased, and they remain consistent given the uncertainties.

Discussion
The present analysis provides significantly improved orbital parameters for β Pictoris b, compared to previous studies, as well as a good estimate of the orbital properties of β Pictoris c. We note that (1) the eccentricity obtained for β Pictoris c (0.27) might be overestimated due to the limited signal-to-noise of the RV data and their temporal sampling (see Shen & Turner 2008), and (2) even with this eccentricity, its orbit would be dynamically stable.
As seen before, the mass of β Pictoris b is poorly constrained and additional long observing RV sequences (typically 3-6 h long each night on three continuous nights) over a few more years are needed to measure this dynamical mass. An improved correction of the pulsations would also allow for a better mass estimate. Previous studies have used HIPPARCOS-Gaia data to constrain the mass of β Pictoris b, but only using one planet in their fit, which is not adapted. We, therefore, attempted to use HIPPARCOS-Gaia to further constrain the mass of this planet. We conclude that the RV data fully drives the solution on the planet masses, as detailed in Appendix B.
On the contrary, β Pictoris c M · sin i is much better constrained. This is because the period of planet c is much smaller than that of planet b, and it has a much better temporal coverage over its period. The present data do not allow for one to constrain the inclination of the planet orbit, but a significant inclination would be very surprising. This work allowed us to estimate the previous and next transit times for both planets. Using the fit of NIRDIFS-GRAV-RV, and assuming a Gaussian distribution N(14,1) for β Pictoris b, we find that the last transit of β Pictoris b occurred on Sept 2017, and that of β Pictoris c occurred on Nov 2017. The next transit of β Pictoris b will occur in July-August 2040 and that of β Pictoris c in March 2021. The dates vary within a few days for the next transit of β Pictoris c and months for β Pictoris b if we change the hypothesis on β Pictoris b's mass, on the star distance, and even more if we change the "thresholded" RV (see above).
Given the mass of β Pictoris c and using COND03 models, we find an absolute magnitude of 12.7 at H band and, hence, we expect a 10.7 mag contrast with the star. Given the measured magnitude of β Pictoris b at H band, we expect a contrast of 0.8-0.9 mag between both planets (hence a flux ratio of 6-8). We examined the possibility of direct detection of β Pictoris c in our SPHERE data, using the positions computed with the orbital parameters found in the previous section. Even if this is not an independent confirmation of its presence, detecting a signal at the planet position would be useful because it would set the relative inclination of the two orbits and provide an estimate of the relative contrast.
It appears that at most epochs, the projected position of β Pictoris c is behind the coronagraphic mask. The projected separation is greater than 100 mas in the following months: Dec. 2014;Feb. 2015;Sept. Oct, and Nov. 2016;Sept., Oct., and Dec. 2018;Dec. 2019;and Jan. 2020. Yet, in Sept., Oct., and Nov. 2016and Sept., Oct., and Dec. 2018, the two planets are within 3 pixels (with a separation down to 9 mas in Nov. 2016) in the plane of the sky. We are then left with four epochs in which the planet is expected to be out of the coronagraph and sufficiently far from b: Feb. 2015 and Dec. 2018, when it should be NE of the star, as well as Dec. 2019 and Feb. 2020, when it is expected to be SW of the star. An examination of the corresponding images does not clearly reveal β Pictoris c. Further work will include the combination of all data set, for example, as done in Le Coroller et al. (2020) to try and directly detect β Pictoris c in the SPHERE data.

Additional planets in the β Pictoris system
Here, we used all available (25) high contrast IRDIS images obtained with NACO and SPHERE and combined them with our 7000+ pulsation-corrected RVs to estimate the probability of the presence of other planets in the system. The approach is the same as the one described in Lagrange et al. (2018) where we used 7 NACO and 2100 pieces of RV data available up to 2016. In a nutshell, we generated (Monte Carlo) planets with various orbital properties and masses and tested whether these simulated planets are detected at least for one epoch in direct imaging or if they are detected in RV data. In direct imaging, the planet is detected if its flux (derived from COND03 brightness-mass relationships, Baraffe et al. 2003) is higher than the three-sigma above noise in a given image. For the RV, we computed the RV series associated to each simulated planet at the time of observations and we tested whether the planet is or is not detectable by comparison with the observed RV time-series (see details in Lannier et al. 2017). In Lagrange et al. (2018), we used NACO L' data only. We upgraded the code to be able to combine images obtained with different instrumental set-ups. The results are given in Fig. 7. The comparison with Fig. 5 of Lagrange et al. (2018) shows the tremendous improvement in detection limits, with detection limits reduced by a factor of 2-4 at almost all separations. To obtain tighter limits, we subtracted the β Pictoris b and c RV signals from the pulsation-corrected data, using the parameters found in the previous section. The result is shown in Fig. 7. From this figure, we can exclude, with a probability greater than 90%, planets that are more massive than about 2.5 M Jup closer than 3 au, and those that are more massive than 3.5 M Jup between 3 and 7.5 au. The intermediate region will be further constrained by the Gaia final release.
The present study shows that there are no massive planets in the system, other than β Pictoris b and c, and that the overall dynamics of the β Pictoris system is controlled by β Pictoris b and c. This property could help us to achieve a better constraint on their orbits. The β Pictoris b and c system can indeed be considered as dynamically isolated from the rest of the disk. A stability analysis that eliminates the unstable configurations could help us refine our knowledge of the orbits.
Using the estimates of Lazzoni et al. (2018) of chaotic zones of planets (Mustill & Wyatt 2012), we find that the presence of β Pic b and c with masses of 11 and 7.8 M Jup , sma of 9.8 and 2.7 au, and eccentricities of 0.08 and 0.28, respectively, is enough to clear out the material in the inner 10 au region. We now focus on the region beyond 10 au.
The planetesimal belt in the β Pic disk is located between 50-130 au (Augereau et al. 2001;Dent et al. 2014). Assuming that the belt initially extended down to the inner region, we could set a lower limit to the masses of planets that could be hidden in the inner region and that could have cleared a gap in about A18, page 11 of 17 A&A 642, A18 (2020) 20 Myr. We based our calculations on the numerical simulations by Shannon et al. (2016), which show that two planets of a minimum of 0.3 M Jup would be needed between β Pic b and the inner location of the outer planetesimal belt to explain the cavity. There can, therefore, be two additional planets with masses between 0.3 and 2.5 M Jup between 10 and 50 au that could be explored further observationally using more SPHERE data and combining them as in this present work, or more sensitive instruments on the JWST or on the ELT. We note that gas observations in β Pic Cataldi et al. 2018) show that the atomic gas does not extend in the inner regions, as expected from previous models (Kral et al. , 2017, and it is mainly colocated with the outer planetesimal belt, which may be a sign that viscous spreading is halted by a planet inward of 50 au, as demonstrated in Kral et al. (2020).

Concluding remarks and perspectives
We obtained new astrometric data on β Pictoris b with SPHERE in high contrast direct imaging and with GRAVITY in long base interferometry, as well as new RV data. The new RV data confirm the detection of β Pictoris c, and the modulation of the RV data by β Pictoris b is becoming visible.
We combined these HCI, GRAVITY, and RV data, as well as previously published data, in a single MCMC analysis to constrain the orbital properties of both planets. The orbital properties of β Pictoris b are driven by the GRAVITY measurements, which have much smaller uncertainties than the HCI ones. This illustrates the exquisite complementarity between both techniques. We confirm that β Pictoris b has a small but nonzero eccentricity (0.09 ± 0.01) and it does not transit the star. The orbital properties of β Pictoris c are well constrained, provided its orbital plane indeed coincides with that of β Pictoris b. The mass of β Pictoris c is well constrained at 7.8 ± 0.4 M Jup , assuming both orbits of the planets are coplanar.
While the orbital properties of β Pictoris b are now well constrained, its mass is not. This is because the baseline of the RV data is still shorter than the planet period, and the data are noisy. Additional data recorded with our updated observing strategy will allow constraining β Pictoris b's mass.
The absolute astrometry HIPPARCOS-Gaia data are consistent with the solutions presented here at the 2σ level, but these solutions are fully driven by the relative astrometry plus RV data. HIPPARCOS-Gaia data should help to further constrain β Pictoris b's mass once more detailed Gaia data are available.
Using the NACO and IRDIS HCI images as well as the RV data, once the pulsations and the signals of β Pictoris b and c were corrected for, we computed the detection limits of additional planets in the system. We can now exclude the presence of planets that are more massive than about 2.5 M Jup closer than 3 au, and more massive than 3.5 M Jup between 3 and 7.5 au. This represents an improvement of more than 5-10 in detection limits compared to our previous study in the 1-10 au region. Between 10 and 100 au, we exclude planets with masses of 1-2 M Jup . We conclude that the dynamics of the system as a whole is controlled by β Pictoris b and c, and that additional planets with masses equal to or less than Jupiter mass have to be present to explain the inner void of material in the 10-70 au region. They certainly represent interesting targets for forthcoming JWST and ELT observations. Finally, given the properties found in this paper, β Pictoris b and c can explain the relative void of material in the inner 10 au disk.

Appendix B: On the use of HIPPARCOS-Gaia data
We used the Template Model Builder (TMB 8 , Kristensen et al. 2016), which implements automatic differentiation (AD), to robustly and quickly solve nonlinear models. This fast method (obtained using the DR2 nominal scanning law provided within the Gaia auxiliary data 9 ), we simulated the reflex motion of the primary due to the presence of the two companions. Knowing, through the scanning law, the scan angle and the parallax factor of each observatoin, we simulated fake intermediate astrometric data (IAD) in the HIPPARCOS way and the five parameter solution was simply derived through a linear least squares solution (e.g., van Leeuwen & Evans 1998). Generating such a solution for each MCMC trial can be cumbersome, but it is quick for the AD method.
To be able to combine the HIPPARCOS IAD of van Leeuwen (2007) with Gaia DR2 data, we needed to first put both systems in the same reference frame, in terms of proper motions (Lindegren et al. 2018;Brandt 2018;Kervella et al. 2019) and of the parallax zero point (Arenou et al. 2018). For simplicity in the handling of the HIPPARCOS IAD, we converted the Gaia DR2 data to the HIPPARCOS reference system. To derive the rotations, we used the 1 arcsec cone search HIPPARCOS-Gaia cross-match described in Marrese et al. (2019) and selected only the HIPPARCOS stars not already known to be binaries (same sample as used in Arenou et al. 2018). We tested that using this full sample and we obtain rotations consistent with Brandt (2018) and Kervella et al. (2019). We found a significant variation of the rotation with magnitude and a smaller one with color. As we do not have enough bright stars in common to select on both magnitude and color, we concentrated only on the main magnitude effect and selected all stars with G < 5. We derived the following rotations: the first one is from the proper motion derived from the Gaia-HIPPARCOS position difference divided by the ∼24-yr baseline (HG) to the HIPPARCOS proper motion system: ω HG2Hip = (0.071, −0.142, −0.057) with an uncertainty on each of those coefficients of around 0.02; the second one is from the Gaia DR2 proper motions to the HIPPARCOS ones: ω Gaia2Hip = (−0.168, −0.255, −0.126) with an uncertainty of 0.03. Similarly, while the global parallax zero point difference between Gaia and HIPPARCOS is −0.12 mas (Arenou et al. 2018), it increases for bright stars. We derived for our sample with G < 5: G − H = 0.086 mas. We also know that the errors are underestimated in a magnitude-dependant way (Lindegren et al. 2018;Arenou et al. 2018) with an additive term for HIPPARCOS and an inflation factor for Gaia (Brandt 2018). Using our bright sample, we derived an additive term of σ = 0.36 mas and an inflation factor of the Gaia errors of 2.1. Those were applied to the covariance matrix of the Gaia DR2 to derive the log likelihood of the tested parameters. For TMB code simplicity, we neglected the impact of the Rv of Beta Pic on the position propagation from HIPPARCOS to the Gaia epoch (which is of 0.05 mas).
All the parameters were let free without priors. When we did not fit the HIPPARCOS IADs, we fixed the parallax to the HIPPARCOS value. Else the five astrometric parameters derived from our code correspond to the system barycenter at J1991.25. As a default starting point for system parameters, we used those of Table C.3, which are unconstrained and a mutual inclination of Nielsen et al. (2020), the HIPPARCOS values for the five astrometric parameters, and zero for the RV offset V0 10 . We obtained the same results as with the MCMC method within the quoted uncertainties. We also successfully tested the convergence robustness of the TMB method by using the results of the MCMC chains as starting points. From those tests, we conclude that the RVs are fully driving the solution concerning the mass of the components. Adding the HIPPARCOS-Gaia constrain does not change the results significantly. The HIPPARCOS-Gaia data are consistent with the solutions presented here at the 2σ level. 10 We also implemented the possibility to add an additional offset after a given epoch for test purposes.
A18, page 15 of 17 A&A 642, A18 (2020)  Notes. The coefficient ρ is the correlation between the error in RA and in Dec. The covariance matrix can be reconstructed using σ ∆RA 2 and σ ∆Dec 2 on the diagonal, and ρσ ∆RA σ ∆Dec off-diagonal. The * indicates new, unpublished data thus far.