Precise Dynamical Masses and Orbital Fits for β Pic b and β Pic c

G. Mirek Brandt; Timothy D. Brandt; Trent J. Dupuy; Yiting Li; Daniel Michalik

doi:10.3847/1538-3881/abdc2e

1. Introduction

β Pictoris b (β Pic b) was among the first exoplanets to be directly imaged (Lagrange et al. 2010). Since then, it has been observed dozens of times, resulting in photometry spanning the near-infrared (Quanz et al. 2010; Bonnefoy et al. 2011, 2013; Currie et al. 2011; Males et al. 2014), low-resolution spectroscopy (Chilcote et al. 2015, 2017), and even medium-resolution spectroscopy (Snellen et al. 2014; Gravity Collaboration et al. 2020). Part of the system's importance derives from the well-measured age of β Pic A. The host star is the defining, highest-mass member of the β Pictoris moving group (Navascués et al. 1999; Zuckerman et al. 2001), which has concordant age determinations of ∼20 million years (Myr) from lithium depletion boundary measurements (Shkolnik et al. 2012; Binks & Jeffries 2014), isochrone analysis (Bell et al. 2015), and kinematic traceback (Miret-Roig et al. 2020).

Combining this well-measured age with an independently measured mass and luminosity can constrain β Pic b's initial supply of thermal energy and provide clues to β Pic b's formation mechanism (Marley et al. 2007; Fortney et al. 2008; Spiegel & Burrows 2012; Marleau & Cumming 2014).

Dynamical mass measurements of β Pic b became feasible thanks to absolute astrometry from the Hipparcos (Perryman et al. 1997; van Leeuwen 2007) and Gaia (Gaia Collaboration et al. 2016; Lindegren et al. 2018) missions. Since Gaia's second data release, a number of authors derived masses and orbits of ∼10–15 M_Jup (Snellen & Brown 2018; Dupuy et al. 2019; Nielsen et al. 2020). The picture was recently complicated and enriched by the discovery of a second companion, β Pic c, orbiting roughly 3 au from the host star (Lagrange et al. 2019b; Nowak et al. 2020), interior to the ≈10 au orbit of β Pic b. β Pic c was first discovered using radial velocities (RVs) alone (Lagrange et al. 2019b), but a direct detection with GRAVITY would soon be reported. Prior to the direct detection and publication of relative astrometry of β Pic c, Nielsen et al. (2020) performed a joint orbital fit to the β Pic system. They obtained a mass of 9.4 ± 1 M_Jup for β Pic c and ${8.3}_{-2.6}^{+2.5}\,{M}_{\mathrm{Jup}}$ for β Pic b (both in agreement with cooling models) but a poor constraint (roughly ±13°) on the inclination of β Pic c (Nielsen et al. 2020) due to the lack of relative astrometry.

Nowak et al. (2020) and Lagrange et al. (2020) detected β Pic c with GRAVITY and fit a two-planet Keplerian model to the β Pic system. When these authors adopted an uninformative prior on the mass of β Pic b, their best-fit dynamical mass measurements were 3.2 M_Jup (Lagrange et al. 2020) and 5.6 ± 1.5M_Jup (Nowak et al. 2020). As the authors noted, such low masses are incompatible with cooling models. Cooling models predict much higher masses that are needed to produce the observed flux (e.g., Baraffe et al. 2003; Spiegel & Burrows 2012). Nowak et al. (2020) ultimately adopted a prior of 15 ± 3 M_Jup while Lagrange et al. (2020) used a prior of 14 ± 1 M_Jup. With these priors, the posterior masses are shifted to near ∼10 M_Jup. The necessity to use such an informative prior indicates a tension between the dynamical constraints and model predictions from spectral analyses, as noted by Nowak et al. (2020).

In this paper, we present precise masses and orbits of β Pic b and β Pic c without the need for informative priors on the planets' masses. Our inferred masses are compatible with a range of cooling model predictions and incorporate the new GRAVITY relative astrometry. We structure the paper as follows. In Section 2, we describe the data that we adopt and the method that we use to fit the system's orbit; we include a full validation on a synthetic data set produced using N-body integration. We present our results in Section 3, including the masses and orbital parameters of both planets, an assessment of the relative astrometry, predicted positions of β Pic b and β Pic c, and N-body results. We discuss the details and implications of our work in Section 4. We summarize our results and conclude in Section 5.

2. Data and Fitting

2.1. Data

The available data for the β Pic system comprise more than 15 years of RVs of β Pic A and relative astrometry for β Pic b, three epochs of relative astrometry for β Pic c, a single RV of β Pic b relative to β Pic A, and absolute astrometry of β Pic A from Hipparcos and Gaia. In this section we summarize each of these.

There are several sources and numerous measurements of relative astrometry for β Pic b, and three recent measurements for β Pic c. We use all relative astrometry, which is comprised of measurements from Near-Infrared Coronagraphic Imager (NICI) on Gemini-South (Nielsen et al. 2014), NACO on the Very Large Telescope (VLT; Currie et al. 2011; Chauvin et al. 2012), MagAO on Magellan (Nielsen et al. 2014), Gemini Planet Imager (GPI) on Gemini South (Wang et al. 2016; Nielsen et al. 2020), Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) on the VLT (Lagrange et al. 2019a), and GRAVITY on the VLT (7 measurements of β Pic b, 3 of c; Lagrange et al. 2020; Nowak et al. 2020). This corresponds to the Case 6 relative astrometry data set of Nielsen et al. (2020) plus recent observations by GRAVITY.

GRAVITY measurements of β Pic b clustered near 2020 disagree internally by ∼2σ and result in an unacceptable reduced χ² (nearly 3) on the position angle (PA) in the final fit (see Section 3.2). We therefore inflate the errors on the seven GRAVITY measurements of β Pic b by a factor of two to make the PA reduced χ² acceptable and bring the PA of the 2020 measurements into internal agreement.

We use the RVs of β Pic A as presented in Vandal et al. (2020), which are corrected for pulsations via a Gaussian process. We add the five new RVs presented in Lagrange et al. (2020) that are not in the data set of Vandal et al. (2020).⁶ We also use the single measurement of the relative RV of β Pic b and β Pic A from Snellen et al. (2014).

We use the absolute astrometry of the Hipparcos–Gaia Catalog of Accelerations (HGCA; Brandt 2018). These astrometric measurements adopt the Gaia DR2 parallax values as priors to all Hipparcos data. Gaia is usually much more precise than Hipparcos, but β Pic A is at the saturation limit of Gaia (G-band magnitude of 3.7). This strongly impacts the astrometric performance of Gaia (Lindegren et al. 2018). Thus, the formal parallax uncertainties of the two missions are comparable (assuming a substantial error inflation to the parallax of the Hipparcos re-reduction in line with the HGCA's inflation of proper-motion errors). Because the HGCA adopts Gaia parallaxes as a prior, we take the Gaia DR2 parallax value of 50.62 ± 0.33 mas (Lindegren et al. 2018) as our prior for the orbital fit. This value is consistent to within 1% with the Hipparcos values (Perryman et al. 1997; van Leeuwen 2007). Regardless, the precise distance to β Pic does not drive our results.

The HGCA argues for a factor of ∼2 inflation for all Gaia DR2 proper-motion errors. We further inflate the HGCA Gaia DR2 proper-motion errors on β Pic by another factor of 2 (a net factor of ∼4 over the Gaia DR2 errors). This is due to systematics in the astrometric fit for very bright stars and is justified by the black histogram (worst 5% of stars) in Figure 9 of Brandt (2018). The Gaia DR2 proper motion has a negligible impact on our results with this large error inflation.

2.2. Orbit Code

We use orvara (Brandt et al. 2020b) along with htof (Brandt et al. 2020a; Brandt & Michalik 2020) to fit for the motion of the β Pic system. orvara fits one or more Keplerian orbits to an arbitrary combination of RV, relative, and absolute astrometry. For the present analysis, we added the ability to fit the single relative RV measurement by Snellen et al. (2014). orvara treats the full motion of the system as a linear combination of Keplerian orbits: an orbit between β Pic b and the combined β Pic A/c system, and a second Keplerian orbit between β Pic A and β Pic c. When computing relative astrometry between β Pic A and β Pic c, orvara neglects interactions with β Pic b. For relative astrometry between β Pic A and β Pic b, orvara computes the displacement of β Pic A from its center of mass with β Pic c and adds this to the displacement of β Pic b from the center of mass of the β Pic A/c system. In other words, orvara only adds astrometric perturbations due to inner companions, not due to outer companions. For the RVs and absolute astrometry of β Pic A, orvara adds the perturbations from the two Keplerian orbits. The perturbation from planet c on the relative RV measurement is negligible.

orvara uses htof to derive positions and proper motions from synthetic epoch astrometry relative to the system's barycenter. htof uses the known Hipparcos observation times and scan angles and the predicted Gaia observation times and scan angles⁷ (with dead times removed) and solves for the best-fit position and proper motion relative to the barycenter. orvara then compares these positions and proper motions to the equivalent values in the HGCA.

orvara marginalizes out the RV zero-point, the parallax, and the barycenter proper motion. We fit a total of 16 parameters to the system using Markov Chain Monte Carlo (MCMC) with ptemcee (Foreman-Mackey et al. 2013; Vousden et al. 2016). These are the six Keplerian orbital elements for each of planets b and c, the mass of each companion, the mass of β Pic A, and a RV jitter to be added in quadrature with the RV uncertainties. We adopt uninformative priors on all parameters: uniform priors on all parameters except for RV jitter (a log-uniform prior) and inclination (a geometric prior).

2.3. Validation

Given the approximate treatment of the three-body problem in orvara, we test its fidelity on data integrated forward using REBOUND (Rein & Liu 2012). We initialize a 1.8 M_⊙ star with two planets of 9 and 8 M_Jup; we give these planets the best-fit orbital elements of β Pic b and c, respectively, found by Lagrange et al. (2020). We then integrate the system forward to produce synthetic RVs and relative astrometry for both companions with REBOUND. We take 52 measurements of relative astrometry for each planet distributed over 17 years, each of which has the 100 μas precision typical of GRAVITY (Gravity Collaboration et al. 2020; Lagrange et al. 2020). We fit 52 RV points, each with an uncertainty of 1 m s⁻¹. We add a single relative RV between β Pic b and A (the synthetic analog to the relative RV of Snellen et al. 2014) with an uncertainty of 1 km s⁻¹. We then fit these synthetic data with orvara.

orvara is able to fit all data satisfactorily. Figure 1 shows that the unmodeled three-body effects are a factor of ∼5 below the level detectable by GRAVITY (see the middle panel). The superposition of the two Keplerian orbits shows up in the relative astrometry of β Pic b, where synthetic GRAVITY observations clearly detect the orbit of β Pic A about its center of mass with β Pic c (bottom panel of Figure 1). Unmodeled RV residuals are well below 1 m s⁻¹ (the reduced χ² of the RV fit is 0.05). We derive masses that agree well, but not perfectly, with the input masses: the derived masses of β Pic A and β Pic b are each ∼3% larger than their true values. These systematics are a factor of ∼5 lower than the uncertainties we derive for β Pic c in the following section and are negligible for β Pic b.

orvara returns two-body elements for each planet about the star. In the three-body system that we initialized in REBOUND, the two-body input orbital elements (semimajor axis, eccentricity, etc.) cease to have a strict meaning unless a primary is specified (e.g., the barycenter or β Pic A). However, we still expect the recovered orbital elements to roughly be equal to those that were used as inputs. We expect the semimajor axes to be close but not exactly equal to the inputs, because, e.g., the input semimajor axis of β Pic b was defined relative to the barycenter of β Pic A and c—yet β Pic b will orbit the total system barycenter during integration. Likewise, we expect the argument and time of periastron to be biased slightly. Elements like the inclination i and PA of the ascending node Ω should be returned exactly—the three-body interactions should not rotate the orientation of either orbit over a ∼20 yr integration.

We find that orvara recovers i and Ω exactly, with a residual less than 10⁻³ of a degree (nearly equal to the formal error) on both. Although unexpected, we recover the eccentricity exactly: the residual is less than 10⁻⁴ and the formal error is 2 × 10⁻⁴. The three elements recovered with biases follow. The argument of periastron and mean longitude at the reference epoch are recovered to within 02. The semimajor axes of both planets are recovered to within 0.1 au.

We conclude that our approximation to the three-body dynamics is more than sufficiently accurate for the β Pic system: the biases induced in the parameters inferred from the test data are much smaller than the formal errors on the measured parameters. Our accounting of only inner companions when perturbing relative astrometry recovers the masses to within a few percent. Figure 1 shows that a full N-body integration of the β Pic system will remain unnecessary even with future GRAVITY relative astrometry.

3. Results

We infer masses and orbital parameters using a parallel-tempered MCMC with 15 temperatures; for each temperature we use 100 walkers with one million steps per walker.⁸ Our MCMC chains converged after 40,000 steps; we conservatively discard the first 250,000 as burn in and use the remainder for inference.⁹

We check convergence informally by confirming that we obtain the same posterior distributions, for every parameter, from any several percent portion of our chains. Next, the acceptance fraction of the coldest chain is satisfactory (∼0.15). Lastly, multiple MCMC analyses starting with different, and in many cases poor, initial guesses converge to the same posterior distributions. We quantitatively confirm convergence with the Gelman–Rubin Diagnostic (GRD; Gelman & Rubin 1992; Roy 2019). Perfect convergence for a parameter is suggested if the GRD is 1, and a common threshold adopted for convergence is 1.1 (Roy 2019). Our chains have GRD values better than 1.0001 for all parameters, although one should keep in mind that the GRD was designed for chains with independent walkers.

3.1. Orbital Analysis of the β Pic System

Table 1 lists the six Keplerian orbital elements for both β Pic b and β Pic c, along with the other five fitted parameters.

Table 1. Posteriors of the β Pic System from an orvara MCMC Analysis

Parameter	Prior Distribution	Posteriors ±1σ
Stellar mass	Uniform	1.83 ± 0.04 M_☉
Parallax (ϖ)	50.62 ± 0.33 mas (Gaia DR2)	50.61 ± 0.47 mas
Barycenter proper motions^a	Uniform	μ_α = 4.80 ± 0.03 mas yr⁻¹ and μ_δ = 83.87 ± 0.03 mas yr⁻¹
RV zero-point	Uniform	33 ± 13 m s⁻¹
RV jitter	Log uniform over (0, 300 m s⁻¹)	50 ± 8 m s⁻¹

Parameter	Prior Distribution	Posterior on β Pic b ± 1σ	Posterior on β Pic c ± 1σ

Semimajor axis (a)	Uniform	10.26 ± 0.10 au	${2.738}_{-0.032}^{+0.034}$ au
Eccentricity (e)	Uniform	0.119 ± 0.008	^b ${0.21}_{-0.09}^{+0.16}$
Inclination (i)	$\sin i$ (geometric)	88.94 ± 0.02 degrees	89.1 ± 0.66 degrees
PA of ascending node (Ω)	Uniform	211.93 ± 0.03 degrees	${211.1}_{-0.2}^{+0.3}$ degrees
Mean longitude at t_ref (λ_ref)	Uniform	−36.7 ± 0.9 degrees	$-{50}_{-14}^{+13}$ degrees
Planet mass (M)	Uniform	${9.3}_{-2.5}^{+2.6}{M}_{\mathrm{Jup}}$	8.3 ± 1.0M_Jup

Argument of Periastron (ω)	(derived quantity)	${22.6}_{-2.9}^{+2.8}$ degrees	${119}_{-7.0}^{+30}$ degrees
Periastron time (T₀)	(derived quantity)	${2456656}_{-64}^{+61}$ BJD	${2455789}_{-63}^{+95}$ BJD
Period	(derived quantity)	${8864}_{-113}^{+118}$ days	${1222}_{-17}^{+18}$ days
		${24.27}_{-0.31}^{+0.32}$ years	${3.346}_{-0.045}^{+0.050}$ years

`orvarar` eference epoch (t_ref)	2455197.50 BJD	⋯	⋯

Notes. Orbital elements all refer to the orbit of the companion about the barycenter. The orbital parameters for β Pic A about each companion are identical except ω_A = ω + π. We use ± when the posteriors are Gaussian. In the case of non-Gaussian posteriors we denote the value by median ${}_{-l}^{+u}$ where u and l denote the 68.3% confidence interval about the median. The reference epoch t_ref is not a fitted parameter and has no significance within the fit itself, it is the epoch at which the mean longitude (λ_ref) is evaluated.

^a μ_α and μ_δ refer to the proper motions in right-ascension and decl., respectively.^bThe posterior on the eccentricity of β Pic c is not Gaussian. However, eccentricities below 0.1 and above 0.7 are strongly disfavored (see Figure 3).

Download table as: ASCII Typeset image

Every fitted element of β Pic b results in a nearly Gaussian posterior (see Figure 2). The elements of β Pic c are also well constrained except for eccentricity and the mean longitude at the reference epoch λ_ref. The mean longitude at the reference epoch is poorly constrained because of the poor constraint on the eccentricity, which results from having only three relative astrometric measurements closely spaced in time. We show the variances and covariances between the fitted parameters in Figure 3 for β Pic c as a corner plot. There is a modest covariance between the semimajor axis and eccentricity resulting from the short time baseline of relative astrometry on β Pic c.

**Figure 2.** Best-fit orbital elements for β Pic b from the `orvara` MCMC chain. Orbital elements are with respect to the star. The elements, in the same order as plotted, are: the primary mass in solar masses, M_pri; the planet mass in Jupiter masses, ${M}_{\sec };$ the semimajor axis in astronomical units, a; the eccentricity, e; the inclination in degrees, i; the PA of the ascending node in degrees, Ω; and the mean longitude at the reference epoch (2455197.50 BJD) in degrees, λ_ref. In the 1D histograms, the vertical dashed lines about the center dashed lines give the 16% and 84% quantiles around the median. In the 2D histograms, the contours give the 1σ, 2σ, and 3σ levels.
Download figure:
Standard image High-resolution image

**Figure 3.** Best-fit orbital elements for β Pic c. See Figure 2 for the description.
Download figure:
Standard image High-resolution image

The best-fit orbit and nearby (in parameter space) orbits agree well with all data: the pulsation-corrected RVs, the Snellen et al. (2014) relative RV, the relative astrometry from VLT/NACO, Gemini-South/NICI, Magellan/MagAO, Gemini-South/GPI, and GRAVITY, and absolute astrometry from the HGCA.

Figure 4 shows the agreement between the calibrated Hipparcos and Gaia proper motions from Brandt (2018) and the best-fit orbit. The sum of the χ² of the fits to both proper motions is very good (nearly 1, see Table 4). There are six measurements, but the unknown barycenter proper motion removes two degrees of freedom. The reflex motion of β Pic c with a period of ∼three years is clearly seen, as well as the long-term oscillation from the ∼24 yr orbit of β Pic b. Here the constraining power of the Hipparcos proper motion is visible: the Hipparcos proper motion is much more precise than that of Gaia DR2 for β Pic b and can exert a sizable tug on the mass and mass uncertainty of β Pic b.

Figures 5 and 6 show the agreement in relative separation and PA from our set of relative astrometry (Case 3 from Nielsen et al. (2020) plus the seven GRAVITY measurements on β Pic b and three GRAVITY measurements on β Pic c). Figure 7 shows the agreement between the RVs from Vandal et al. (2020) and Lagrange et al. (2020) and the best-fit orbit. The jitter parameter found by the MCMC analysis is 50 ± 8 m s⁻¹. Lower masses for β Pic c slightly favor lower eccentricities. The Snellen et al. (2014) relative RV χ² is 1.7 (indicating a ∼1.3σ residual). However, our posteriors are completely identical within rounding if we exclude the single Snellen et al. (2014) measurement.

**Figure 5.** Left: relative separation of β Pic c. Right: PA of β Pic c. All three data points are from GRAVITY (Nowak et al. 2020) and are not error inflated. The best-fit orbit is shown in black. A random sampling of other orbits from the MCMC chain are shown and are color coded by the mass of β Pic c.
Download figure:
Standard image High-resolution image

**Figure 6.** Left: relative separation of β Pic b. Right: PA of β Pic b. The GRAVITY errors have been inflated by a factor of 2 to make the reduced χ² of the fit acceptable. A random sampling of orbits from other MCMC steps are shown and are color coded by the mass of β Pic b. The best-fit orbit is shown in black.
Download figure:
Standard image High-resolution image

**Figure 7.** The best-fit orbit (black) agrees well with the observed β Pic pulsation-corrected RVs. β Pic c has an eccentricity of e = 0.30 in the best-fit orbit while b has e = 0.120. Top panel: the observed RVs overplot with the best-fit orbit and a random sampling of other orbits from the MCMC chain. Bottom panel: the RV residuals with respect to the best-fit orbit. Both panels: the random sampling of other orbits from the MCMC chain are color coded by the mass of β Pic c. RVs are from Vandal et al. (2020) with the most recent five points from Lagrange et al. (2020). The black error bars give the observed errors reported by Vandal et al. (2020) and Lagrange et al. (2020). The red error bars include the best-fit jitter of ∼50 m s⁻¹ added in quadrature to the observed errors.
Download figure:
Standard image High-resolution image

We display an additional corner plot in Figure 8 that showcases select covariances between the orbital parameters of β Pic b and β Pic c. The inferred mass of each planet is relatively insensitive to the orbital parameters of the other (see the two appropriate covariances in the left hand columns of Figure 8). In particular, the mass of β Pic b is nearly independent of the mass of c. However, owing to the three-body interaction between the planets, the inferred eccentricity of β Pic b varies slightly with the eccentricity of the inner planet, β Pic c. Improved relative astrometry on β Pic b mildly improves constraints on the eccentricity of β Pic c; an identical orbital fit excluding the SPHERE relative astrometry on β Pic b results in a slightly worse eccentricity constraint on β Pic c. The inferred semimajor axis of β Pic c covaries modestly with β Pic b's eccentricity. Despite uncertainties in the eccentricity of β Pic c, we find that β Pic b and β Pic c are coplanar to within a half-degree at 68% confidence and coplanar to within one degree at 95% confidence.

We use our new constraints on the orbital parameters of β Pic b and c to predict their on-sky positions over the next 5 years at 15 day intervals. Tables 2 and 3 give a truncated version of the predicted positions of β Pic b and c. The supplementary data contain the full tables. β Pic c will be less than ∼50 mas from the star by 2021 March. β Pic c will re-emerge (once again being further than ∼50 mas from the star) in 2021 October. Our predicted positions from our orbit analysis localize both β Pic b and c to within ±40 mas, which is well within the fiber field of GRAVITY (Nowak et al. 2020), at any point over the next 5 years.

Table 2. Predicted Positions of β Pic b

Date	δ	σ_δ	α	σ_α	ρ_{α δ}	Sep	σ_Sep
	(mas)	(mas)	(mas)	(mas)		(mas)	(mas)
2020 Dec 30	351.9	0.5	211.6	0.3	0.951	410.6	0.1
2021 Jan 14	355.2	0.6	213.7	0.4	0.954	414.5	0.1
2021 Jan 29	358.4	0.6	215.8	0.4	0.957	418.4	0.2
2021 Feb 13	361.6	0.6	217.9	0.4	0.959	422.2	0.2
2021 Feb 28	364.7	0.7	219.9	0.4	0.961	425.9	0.2
2021 Mar 15	367.9	0.7	222.0	0.4	0.963	429.6	0.2
2021 Mar 30	370.9	0.8	224.0	0.5	0.965	433.3	0.2
2021 Apr 14	374.0	0.8	226.0	0.5	0.967	437.0	0.3
2021 Apr 29	377.0	0.8	227.9	0.5	0.969	440.6	0.3
⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
2025 Dec 19	460	10	290	7	0.999	540	10

Note. The offsets, and their errors (σ), from the star in right-ascension (α), decl. (δ), and separation (Sep), are given in milli-arcseconds (mas). ρ_{α
δ} is the correlation coefficient between right-ascension and decl. A machine-readable version of this table with all 122 epochs is available with the supplementary data online. This portion is shown here for guidance regarding its form and content.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

Table 3. Predicted Positions of β Pic c

Date	δ	σ_δ	α	σ_α	ρ_{α δ}	Sep	σ_Sep
	(mas)	(mas)	(mas)	(mas)		(mas)	(mas)
2020-12-30	−82	8	−52	5	0.911	97	5
2021-01-14	−76	9	−49	6	0.927	90	5
2021-01-29	−70	10	−45	6	0.939	83	4
2021-02-13	−60	10	−41	7	0.950	75	2
2021-02-28	−60	10	−37	8	0.958	68	3
2021-03-15	−50	10	−33	9	0.965	60	8
2021-03-30	−40	10	−29	9	0.970	50	10
2021-04-14	−40	20	−20	10	0.975	40	20
2021-04-29	−30	20	−20	10	0.978	30	10
⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
2025-12-19	60	20	40	10	0.992	80	20

Note. See the table note of Table 2 for a description of the columns.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

3.2. Assessing Consistency of Relative Astrometry

Table 4 shows quantitatively the goodness of the orbital fit in terms of χ² = ∑(data − model)²/σ² for PA, separation, RV, and the three proper motions. The reduced χ² for all the astrometry (which takes into account the GRAVITY covariances between separation and PA) is 1.06. A good fit should have χ²/N ≈ 1 where N is the number of degrees of freedom.

Table 4. The Goodness of the orvara Orbital Fit to the Various Data in the β Pic System

Data	Points (N)	χ²
Separation	56	65.6
PA	56	59.0
All astrometry	112	^a118.8
RV	41	40.7
β Pic b–A relative RV	1	1.67
Hipparcos μ (HGCA)	2	0.33
Gaia μ (HGCA)	2	0.75
HGCA long baseline μ	2	0.001

Notes. The χ² quoted here include both companions and are for the maximum likelihood orbits. The χ² for μ includes both μ_δ and μ_α. N is the number of data points in the corresponding data set.

^aThis χ² is slightly less than the sum of the χ² in PA and separation because of the covariance between PA and separation in the GRAVITY observations.

Download table as: ASCII Typeset image

Nielsen et al. (2020) argued for a systematic offset between the SPHERE relative astrometry from Lagrange et al. (2019a) and the relative astrometry from Gemini-South/GPI. Nielsen et al. (2020) investigated fitting for an offset in both separation and PA within the SPHERE data. The SPHERE data do appear to be systematically offset in PA relative to the best-fit orbit (see the bottom right panel of Figure 6). However, a fit without the 12 SPHERE observations reduces the χ² in PA and separation by roughly the expected 12 points, suggesting that the data are consistent with the astrometric record. Moreover, the reduced χ² including SPHERE is acceptable (∼65 points of χ² for 56 data points) and so we include SPHERE in our final analysis.

We find evidence for either an underestimate in the PA uncertainties from GPI or an offset in PA between GPI and one or more of the other astrometric data sets (see the Wang et al. 2016 GPI data in the right panel of Figure 6). Removing the 15 GPI relative astrometry measurements decreases the χ² in PA by roughly 40. Using the χ² survival function, a change of that magnitude corresponds to roughly 2.5σ evidence in favor of a PA offset. However, including GPI still results in an acceptable overall χ² (See Table 4), and so we include GPI in our final fit.

Whether or not we include one, both or neither of GPI and SPHERE, our results are nearly identical. The best-fit masses on both β Pic b and c shift by less than 0.5 M_Jup between all three cases, and the confidence intervals on their masses are identical to within 5%. This speaks to the constraining power of the GRAVITY measurements, and to the robustness of our results with respect to the details of how the relative astrometry is analyzed.

In Figure 5, the χ² of the β Pic c fit to the relative astrometry is much less than 1 because the relative astrometry is effectively overfit: the RVs primarily constrain the mass, phase, and semimajor axis of β Pic c while the four remaining orbital parameters have substantial freedom to fit the three relative astrometry points (six coordinates). By contrast, β Pic b is overconstrained by the data and the reduced χ² of the fit is near 1. The right-hand side of the bottom-most panel for both separation and PA in Figure 6 shows the GRAVITY points for β Pic b. GRAVITY points near the same epoch (in both PA and separation) disagree by ≲1σ after error inflation. Without error inflation, GRAVITY observations near the same epoch disagree by ∼2σ and the reduced χ² in PA of the best-fit jumps to nearly 3 for β Pic b.

The three GRAVITY measurements of β Pic c do not have χ² or agreement issues. We leave the errors on β Pic c as they are in Nowak et al. (2020). However, inflating the errors by a factor of 2 on β Pic c does not significantly change our results: the resulting posteriors and errors are identical except for the errors on β Pic c's inclination, which are doubled.

3.3. N-body Simulations

We expect the orbital parameters of β Pic b and c to vary over time due to the mutual influence between these two massive planets. The evolution of the eccentricity and orbit of β Pic b depends heavily on the eccentricity of β Pic c, which is poorly constrained. In Figure 9, we show the evolution of the β Pic system over 0.1 million years (Myr), integrated forward using the ias15 integrator of REBOUND (Rein & Liu 2012; Rein & Spiegel 2015), assuming the median orbital parameters presented in Table 1 for each planet. We vary the eccentricity of β Pic c within the posterior constraints. The gray shaded region shows how the eccentricity of β Pic c and β Pic b could evolve over the next 10⁵ yr. The two planets exchange eccentricity with a period of ∼50,000 yr. We found numerically that the system is stable and the periodic variability in Figure 9 repeats for at least the next 10 Myr.

**Figure 9.** Top panel: the eccentricity evolution of β Pic b computed using `REBOUND`'s ias15 integrator (Rein & Spiegel 2015) for the median (black, 0.21) and 68.3% confident bounds on the eccentricity of β Pic c from Table 1 (0.12 is red and 0.37 is blue). Bottom panel: the eccentricity evolution of β Pic c for its median (black) and 68.3% confident eccentricities. The parameter space spanned by the 68.3% confident range of eccentricities is shaded gray. The two planets exchange eccentricity over a ∼50,000 year cycle.
Download figure:
Standard image High-resolution image

4. Discussion

Our mass measurements for β Pic b and c agree within 1σ compared to previous work by Snellen & Brown (2018), Dupuy et al. (2019), Nielsen et al. (2020), and Vandal et al. (2020). Our analysis is the first to incorporate the new GRAVITY measurements with uninformative priors while obtaining masses in the expected expected range. Our error bars on the mass of β Pic b are larger than all but Dupuy et al. (2019) because, like that work, we adopt the inflated errors on the Hipparcos proper motions as recommended by Brandt (2018). Our mass posteriors do not change if we exclude the Snellen et al. (2014) relative RV measurement. We were unable to reproduce the ≈3 M_Jup and ≈5 M_Jup (when using a uniform prior) mass estimates from Nowak et al. (2020) and Lagrange et al. (2020). Using their slightly different data set, we find ${9.5}_{-1.8}^{+2.0}\,{M}_{\mathrm{Jup}}$ for β Pic b and ${9.2}_{-0.8}^{+1.0}\,{M}_{\mathrm{Jup}}$ for β Pic c.

We corroborate the findings by Nielsen et al. (2020) and Nowak et al. (2020) that β Pic c and β Pic b are coplanar. Nowak et al. (2020) found inclinations for β Pic b and c of 8899 ± 001 and 8917 ± 050, respectively, with a Gaussian prior on the mass of β Pic b. We confirm these inclinations without an informative prior. We find 8894 ± 002 and 891 ± 07.

β Pic is surrounded by an extended debris disk and an inner disk that is slightly misaligned with respect to the primary (Smith & Terrile 1984; Heap et al. 2000). The extended debris disk around β Pic is inclined at 90.0 ± 0.1 degrees (Ahmic et al. 2009; Kraus et al. 2020). β Pic b is thus misaligned by 1.06 ± 0.11 degrees with respect to the debris disk. Our inferred inclination for β Pic c slightly favors misalignment but does not exclude alignment.

Nowak et al. (2020) found that the β Pic system exhibited an oscillating eccentricity for both bodies over a timescale of ≈5 × 10⁴ yr using their orbital parameter posteriors. We find variations in eccentricity over a similar timescale and confirmed numerically with REBOUND that the system is stable for at least the next 10 Myr. We find that it is moderately likely to observe the current eccentricity of the system amidst all the possible eccentricities over a 10 million year time span.

The first observational evidence that β Pic b has a significant, nonzero eccentricity was presented by Dupuy et al. (2019). They discussed the implications of an eccentricity as high as ≈ 0.2 in the context of both single- and multiplanet scenarios; at the time β Pic c was not known. The scenario in which β Pic b formed on a circular orbit but gained eccentricity from interactions with the disk and migrated inward to its current location, with no influence from β Pic c, is still plausible. Such a pathway is available to any sufficiently massive planet. Given that we find that β Pic c is also massive (8.3 ± 1 M_Jup), it may have also opened a gap in the disk, migrating inward and acquiring eccentricity from gravitational interactions with β Pic b and the disk. Indeed, with two such massive planets in close proximity it is natural to expect that both should have significantly nonzero eccentricities by a system age of ≈20 Myr.

Our masses follow from uniform priors, allowing us to independently assess the agreement of the dynamical masses with model predictions. To simplify our model comparisons, we assume an age of 20 Myr for the system, compatible with all available age determinations for the β Pic moving group (Binks & Jeffries 2014; Mamajek & Bell 2014; Miret-Roig et al. 2020). We examine the hot-start Cond (Baraffe et al. 2003) models, the Saumon & Marley (2008) models with a hybrid cloud treatment (which we denote as SM08), and the warm-start Spiegel & Burrows (2012) models (SB12) with hybrid clouds and solar metallicity but a range of initial entropies. We perform our comparisons in the K band, as this is the only measurement available for β Pic c (Nowak et al. 2020). We convert luminosities to K-band magnitudes for the SM08 models using Cond colors at the SM08 effective temperatures.

Figure 10 shows our results. We find that our dynamical mass measurement for β Pic c is consistent with all models except those with low initial entropies (≲10 k_B/baryon). Our dynamical mass for β Pic b is roughly 1σ below the predictions of hot-start models, and rules out cold starts. Similarly to previous work (Dupuy et al. 2019; Vandal et al. 2020), none of the disagreements with models are significant beyond ∼1σ, and the precisions of the dynamical masses are insufficient to distinguish between most of the models shown. Stronger tests of models will require significantly better precision, especially for β Pic b.

As Figure 10 shows, reaching 0.1–0.5 M_Jup levels of precision on the mass of β Pic b is crucial to accurately discern between evolutionary models. The best prospect for improving the mass of β Pic b is long-term RV monitoring over the next decade. Even drastically improved absolute astrometry (e.g., Gaia DR4) will only provide a modest improvement to the mass measurement of β Pic b. If we assume optimistically that Gaia at the end of its mission will achieve the same precision on the G = 3.7 mag β Pic A as it has on G ≈ 6 mag stars (the brightest for which the mission was originally designed), then it would achieve a factor of ∼100 improvement on the proper motion of β Pic A.¹⁰ The uncertainty on the mass of β Pic b would shrink by 35%, to ±1.7 M_Jup, if the proper-motion precision is improved by a factor of 100—using the same MCMC analysis as presented here with otherwise the same data. Assuming more conservatively that Gaia reaches only a factor of 10 better precision on the proper motion of β Pic, the uncertainty on the mass of β Pic b improves by 25%.

5. Conclusions

In this paper we have derived masses and orbits of both planets in the β Pictoris system with uninformative priors. We validated our approach against synthetic data from a full N-body integration. Our masses and orbital parameters are derived from two decades of observational data. The GRAVITY data show clear evidence of the gravitational perturbations of β Pic c ( $P={3.346}_{-0.045}^{+0.050}$ yr) on the orbit of β Pic b relative to A (P = 24.27 ± 0.32 yr). The resulting model-independent masses allow us to compare the observed properties of β Pic b and c with predictions from models of the formation and evolution of giant planets. We summarize our main results below.

1.
We find a mass of ${9.3}_{-2.5}^{+2.6}\,{M}_{\mathrm{Jup}}$ for β Pic b and 8.3 ± 1.0 M_Jup for β Pic c with uninformative priors for all orbital parameters. The mass constraint on β Pic c is superior due to the RVs covering many orbital periods and due to the impact of β Pic c on the relative astrometry of β Pic b.
2.
β Pic b and β Pic c are both consistent with the Spiegel & Burrows (2012) warm-start models with initial entropies of at least 10 k_B/baryon. They are also both consistent with a 20 Myr age under the hot-start COND evolutionary tracks (Baraffe et al. 2003) and the Saumon & Marley (2008) models using a hybrid cloud model. In all cases, consistency with models would favor a mass for β Pic b that is ∼1σ higher than our dynamical measurement.
3.
We find an eccentricity of 0.119 ± 0.008 for β Pic b and ${0.21}_{-0.09}^{+0.16}$ for c. These modest eccentricities could have been generated by interactions with the disk, or via the mutual interactions between b and c. The eccentricity and mean longitude of β Pic c are poorly constrained because there are only three relative astrometric observations, and these are closely spaced in time. There is a mild covariance between the eccentricity of β Pic b and c owing to the three-body dynamics in the system. Additional GRAVITY relative astrometry on β Pic c will help constrain the eccentricity of β Pic b and especially β Pic c.
4.
The mass constraint on β Pic b needs to be improved by a factor of ∼3–5 in order to more reliably constrain its age or formation conditions. Long-term RV monitoring over the coming years or decade is needed for better mass constraints on β Pic b. An improved proper motion from a future Gaia data release will offer up to a 35% better constraint on the mass of β Pic b (assuming Gaia reaches a better precision on the brightest stars).

The new GRAVITY relative astrometry (Lagrange et al. 2020; Nowak et al. 2020) appeared to create tension between dynamical and spectral mass constraints on β Pic b. Our analysis dissolves this tension and results in masses for β Pic b and β Pic c that are consistent with warm and hot-start evolutionary models. Additionally, the system is dynamically interesting—with eccentricities of both planets varying by ∼50% over 10⁴–10⁵ yr timescales. The precision on the masses and eccentricities of β Pic b and c will improve with continued astrometric and RV monitoring. The planets around β Pic A will continue to provide some of the best tests of super-Jovian planet formation and evolution.

G.M.B. is supported by the National Science Foundation (NSF) Graduate Research Fellowship under grant No. 1650114.

We thank Kaitlin Kratter for fruitful discussions and comments. We thank the anonymous referee for the constructive comments that improved the quality of our work.

This work made use of the REBOUND code which is freely available at http://github.com/hannorein/rebound.

This work made use of the orvara code. The exact version we used is available via the githash 51757c58dc78db406a39503b509a929253b9.¹¹

This work made use of the htof code. We used version 0.3.1¹² (Brandt & Michalik 2020).

Software: astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), scipy (Virtanen et al. 2020), numpy (Oliphant 2006; van der Walt et al. 2011), pandas (McKinney 2010; pandas development team 2020), orvara (Brandt et al. 2020b), htof (Brandt & Michalik 2020; Brandt et al. 2020a), REBOUND (Rein & Liu 2012), corner (Foreman-Mackey 2016), Jupyter.

Precise Dynamical Masses and Orbital Fits for β Pic b and β Pic c

Article metrics

Permissions

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction