A Model-Independent Mass and Moderate Eccentricity for $\beta$ Pic b

TL;DR

This study measures the mass and eccentricity of $eta$ Pic b using combined astrometric data, providing a model-independent mass estimate that challenges cold-start formation models and suggests possible orbital migration.

Contribution

First model-independent mass measurement of a directly imaged exoplanet using combined Hipparcos and Gaia data, with implications for planet formation theories.

Findings

Mass of $eta$ Pic b is $13\u00b1 3 M_{Jup}$.

Eccentricity of $eta$ Pic b is $0.24\u00b1 0.06$, higher than previous estimates.

Results are consistent with hot- and warm-start formation models.

Abstract

We use a cross-calibration of Hipparcos and Gaia DR2 astrometry for $\beta$ Pic to measure the mass of the giant planet $\beta$ Pic b $(13\pm3$ $M_{\rm Jup})$ in a comprehensive joint orbit analysis that includes published relative astrometry and radial velocities. Our mass uncertainty is somewhat higher than previous work because our astrometry from the Hipparcos-Gaia Catalog of Accelerations accounts for the error inflation and systematic terms that are required to bring the two data sets onto a common astrometric reference frame, and because we fit freely for the host-star mass $(1.84\pm0.05$ $M_{\odot})$. This first model-independent mass for a directly imaged planet is inconsistent with cold-start models given the age of the $\beta$ Pic moving group $(22\pm6$ Myr) but consistent with hot- and warm-start models, concordant with past work. We find a higher eccentricity $(0.24\pm0.06)$ for $\beta$ Pic b compared to previous orbital fits. If confirmed by future observations, this eccentricity may help explain inner edge, scale height, and brightness asymmetry of $\beta$ Pic's disk. It could also potentially signal that $\beta$ Pic b has migrated inward to its current location, acquiring its eccentricity from interaction with the 3:1 outer Lindblad resonance in the disk.

Tables

Table 1. Table 1 Absolute Stellar Astrometry

Mission	${\mu }_{\alpha \ast }$	$\sigma \left[{\mu }_{\alpha \ast }\right]$	${\mu }_{\delta }$	$\sigma \left[{\mu }_{\delta }\right]$	Corr$[{\mu }_{\alpha \ast },{\mu }_{\delta }]$	$t_{\alpha \ast }$	$t_{\delta }$
	(mas yr-1)		(mas yr-1)			(year)
Hipparcos	4.4	0.4	82.8	0.4	0.002	1991.33	1991.26
Hipparcos–Gaia	4.796	0.027	83.863	0.028	0.025	$\cdots$	$\cdots$
Gaia	2.5	2.5aaGaia DR2 errors have been inflated by a factor of two as recommended by Brandt (2018) for stars like $\beta$ Pic that have large reported proper motion errors in DR2 ($\gtrsim 0.7$ mas yr-1).	82.6	2.5aaGaia DR2 errors have been inflated by a factor of two as recommended by Brandt (2018) for stars like $\beta$ Pic that have large reported proper motion errors in DR2 ($\gtrsim 0.7$ mas yr-1).	0.040	2015.58	2015.67

Table 2. Table 2 MCMC Orbital Posteriors for β 𝛽 \beta Pic b

Property	Median $\pm$1$\sigma$	95.4% c.i.	Prior
Fitted parameters
Companion mass $M_{\mathrm{comp}}$ ($M_{\mathrm{Jup}}$)	${13.1}_{-3.2}^{+2.8}$	7.2, 19.5	$1/M$ (log-flat)
Host-star mass $M_{\mathrm{host}}$ ($M_{☉}$)	$1.84\pm 0.05$	1.74, 1.94	$1/M$ (log-flat)
Parallax (mas)	${51.60}_{-0.39}^{+0.40}$	50.82, 52.37	$\exp \left[-0.5{\left(\left(\varpi -{\varpi }_{\mathrm{DR2}}\right)/\sigma \left[{\varpi }_{\mathrm{DR2}}\right]\right)}^2\right]$
Semimajor axis $a$ (AU)	${11.8}_{-0.9}^{+0.8}$	10.3, 13.7	$1/a$ (log-flat)
Inclination $i$ (∘)	$88.87\pm 0.08$	88.71, 89.04	$\sin \left(i\right)$,
$\sqrt{e}\sin \omega$	$-{0.080}_{-0.029}^{+0.027}$	$-$0.134, $-$0.017	uniform
$\sqrt{e}\cos \omega$	$-0.48\pm 0.05$	$-$0.59, $-$0.36	uniform
Mean longitude at $t_{\mathrm{ref}}=2455197.5$ JD, ${\lambda }_{\mathrm{ref}}$ (∘)	$150\pm 4$	142, 159	uniform
PA of the ascending node $\Omega$ (∘)	$31.65\pm 0.09$	31.48, 31.82	uniform
RV zero point (m s-1)	${73}_{-15}^{+14}$	45, 103	uniform
RV jitter $\sigma$ (m s-1)	$269\pm 6$	257, 281	$1/\sigma$ (log-flat)
Computed properties
Orbital period $P$ (yr)	${29.9}_{-3.2}^{+2.9}$	24.1, 36.8	$\cdots$
Semimajor axis (mas)	${610}_{-50}^{+40}$	530, 700	$\cdots$
Eccentricity $e$	$0.24\pm 0.06$	0.13, 0.35	$\cdots$
Argument of periastron $\omega$ (∘)	${189.3}_{-2.9}^{+3.0}$	182.3, 195.5	$\cdots$
Time of periastron $T_0=t_{\mathrm{ref}}-P\frac{\lambda -\omega }{360{}^{\circ }}$ (JD)	${2456380}_{-60}^{+80}$	2456210, 2456520	$\cdots$
Mass ratio $q=M_{\mathrm{comp}}/M_{\mathrm{host}}$	${0.0068}_{-0.0016}^{+0.0015}$	0.0038, 0.0101	$\cdots$

Table 3. Table 3 Predicted Future Astrometry and Radial Velocities for β 𝛽 \beta Pic b

Epoch		Separation (mas)				PA (∘)				$\Delta \mathrm{RV}$ (km s-1)
		$e=0.10$	$e=0.20$	$e=0.30$		$e=0.10$	$e=0.20$	$e=0.30$		$e=0.10$	$e=0.20$	$e=0.30$
2019 Jan 1		$176\pm 3$	$173\pm 3$	$169\pm 3$		$28.21\pm 0.31$	$27.97\pm 0.33$	$27.79\pm 0.35$		0$2.92\pm 0.23$	$1.24\pm 0.20$	$-0.25\pm 0.19$+
2020 Jan 1		$301\pm 3$	$299\pm 3$	$296\pm 3$		$29.81\pm 0.20$	$29.63\pm 0.21$	$29.49\pm 0.22$		0$5.76\pm 0.29$	$3.68\pm 0.24$	$1.87\pm 0.22$
2021 Jan 1		$407\pm 3$	$412\pm 3$	$413\pm 3$		$30.50\pm 0.15$	$30.32\pm 0.17$	$30.19\pm 0.18$		0$8.07\pm 0.32$	$5.68\pm 0.27$	$3.61\pm 0.25$
2022 Jan 1		$490\pm 3$	$508\pm 3$	$517\pm 3$		$30.92\pm 0.13$	$30.73\pm 0.14$	$30.59\pm 0.15$		0$9.78\pm 0.33$	$7.24\pm 0.28$	$5.00\pm 0.26$
2023 Jan 1		$545\pm 6$	$585\pm 4$	$608\pm 3$		$31.24\pm 0.11$	$31.02\pm 0.12$	$30.86\pm 0.13$		$10.89\pm 0.30$	$8.41\pm 0.28$	$6.10\pm 0.26$
2024 Jan 1		0$571\pm 10$	$642\pm 6$	$684\pm 4$		$31.51\pm 0.11$	$31.25\pm 0.11$	$31.06\pm 0.12$		$11.40\pm 0.24$	$9.22\pm 0.26$	$6.95\pm 0.26$
2025 Jan 1		0$567\pm 16$	$679\pm 9$	$747\pm 6$		$31.78\pm 0.10$	$31.44\pm 0.10$	$31.23\pm 0.11$		$11.29\pm 0.17$	$9.70\pm 0.22$	$7.58\pm 0.25$
2026 Jan 1		0$533\pm 23$	0$695\pm 13$	$794\pm 8$		$32.06\pm 0.11$	$31.62\pm 0.10$	$31.37\pm 0.11$		$10.59\pm 0.18$	$9.86\pm 0.18$	$8.02\pm 0.23$
2027 Jan 1		0$471\pm 30$	0$690\pm 18$	0$827\pm 11$		$32.40\pm 0.14$	$31.80\pm 0.10$	$31.51\pm 0.10$		0$9.30\pm 0.31$	$9.73\pm 0.14$	$8.30\pm 0.20$
2028 Jan 1		0$383\pm 38$	0$664\pm 24$	0$845\pm 15$		$32.88\pm 0.22$	$31.99\pm 0.10$	$31.63\pm 0.10$		0$7.44\pm 0.52$	$9.31\pm 0.13$	$8.42\pm 0.17$

Full text

1. Introduction

Directly imaged planets typically have their masses inferred indirectly from their luminosity and age, using uncalibrated evolutionary models that assume an initial thermal state. Most commonly-used models assume an initially high specific entropy (hot start; e.g., Burrows et al., 1997; Baraffe et al., 2003), but the planet formation process might radiate away a significant amount of energy leading to a much lower initial specific entropy (cold or warm start; e.g., Marley et al., 2007; Spiegel & Burrows, 2012). Furthermore, planet assembly could be slow and only conclude well after the star is formed, in which case young planets could appear even more luminous than hot-start models would predict from the host star’s age. The crucial observations needed to sort out these various possibilities are masses of planets with known age and luminosity.

$\beta$ Pic b was one of the first directly imaged planets to be discovered (Lagrange et al., 2010), and its host star is the namesake of a young moving group of well-determined age ($22\pm 6$ Myr; e.g., Binks & Jeffries, 2014; Shkolnik et al., 2017). We present here a new model-independent dynamical mass for $\beta$ Pic b. We use the methodology of Brandt et al. (2018) to perform a joint orbital analysis of relative astrometry, radial velocities, and host-star astrometry from the cross-calibrated Hipparcos–Gaia Catalog of Accelerations (HGCA; Brandt, 2018). Our new mass is consistent with recent results from Snellen & Brown (2018) but with broader uncertainties owing to our re-assessment of errors reported in Hipparcos and Gaia DR2 catalogs.

2. Data

2.1. Host-Star Astrometry

Brandt (2018) has cross-calibrated Hipparcos and Gaia DR2, placing them on a common reference frame. Figure 1 of Brandt (2018) shows that neither the Hipparcos re-reduction (van Leeuwen, 2007) nor the Gaia DR2 astrometry (Lindegren et al., 2018) are suitable for orbit fitting in their published form: the ensemble of proper motion differences are inconsistent with their formal uncertainties. Moreover, Figure 9 of Brandt (2018) shows that the cross-calibrated HGCA proper motions satisfy the standard assumptions of Gaussianity but that the lowest-precision stars in Gaia (like $\beta$ Pic) have uncertainties that remain underestimated.

HGCA contains three proper motions: a (nearly) instantaneous proper motion near 1991.25, another near 2015.5, and the positional difference between the catalogs scaled by the time between them. The three proper motions are nearly independent. Brandt (2018) also gives the central epoch at which a position was measured; this is the epoch with the minimum positional uncertainty (which differs slightly in right ascension and declination).

Table 1 lists our HGCA proper motions for $\beta$ Pic, the correlation coefficients between proper motion in RA and Dec, and the central epoch for each measurement. $\beta$ Pic is heavily saturated in Gaia data and thus is among the least-precisely measured stars in the HGCA. Figure 9 of Brandt (2018) indicates that the inflated uncertainties of such stars remain underestimated by as much as factor of two. We have therefore doubled the Gaia DR2 proper motion errors beyond the values in the HGCA. For parallax, we adopt the same 60/40 linear combination of the Hipparcos catalogs as the HGCA and add the same 0.20 mas error inflation in quadrature; this results in a value of $51.61\pm 0.39$ mas.

The uncertainties for the Hipparcos proper motions are much larger in the Brandt (2018) catalog than in the Hipparcos re-reduction, though they are slightly smaller than the uncertainties of the original Hipparcos reduction. This is a generic feature of bright stars in the HGCA. As shown in Figure 1 of Brandt (2018), stars with higher-precision proper motions depart most strongly from the standard normal distribution in their residuals. Even for the most precise 20% of stars, a $\sim$60/40 linear combination of the two Hipparcos reductions gives lower residuals than the van Leeuwen (2007) proper motions alone, and further error inflation is necessary to bring the residuals into agreement with a normal distribution. Given the cross-calibration approach used in the HGCA, it would be infeasible to use Hipparcos epoch astrometry, as in Snellen & Brown (2018), and ensure independence of individual measurements.

2.2. Literature Relative Astrometry & Radial Velocities

We consider all available relative astrometry of $\beta$ Pic b in our orbit analysis, setting aside duplicate measurements when the same data have been analyzed separately in the literature. This includes astrometry from VLT/NaCo (Quanz et al., 2010; Bonnefoy et al., 2011, 2013; Chauvin et al., 2012; Absil et al., 2013; Milli et al., 2014), Gemini-S/NICI (Nielsen et al., 2014), Magellan/MagAO (Nielsen et al., 2014), Gemini-S/GPI (Wang et al., 2016), and VLT/SPHERE (Lagrange et al., 2018a). This comprises 50 measurements spanning sixteen years, with two observations on the northeastern side of the orbit (in 2003 November and 2018 September).

The radial velocity of the host star has been monitored from 2003–2011 with the HARPS spectrograph (Lagrange et al., 2012). We use all 1049 individual published measurements and account for the substantial intrinsic “jitter” that is expected for a young star like $\beta$ Pic. We also use the measurement of the planet’s relative radial velocity ($\Delta{\rm RV}={\rm RV}_{\rm comp}-{\rm RV}_{\rm host}$) from Snellen et al. (2014) in our orbit fit.

3. Orbit Analysis

Relative astrometry from direct imaging has already been shown to constrain many orbital parameters of $\beta$ Pic b given the long time baseline and intensive monitoring (e.g., Wang et al., 2016; Lagrange et al., 2018a). Therefore, as a first step we fit the relative astrometry with a standard seven-parameter Keplerian orbit in order to assess any systematics in combining astrometry from many different instruments and data reduction methods. We found an unreasonably large $\chi^{2}$ of 165 for 93 degrees of freedom (dof), $p(\chi^{2})=6\times 10^{-6}$, when taking all reported astrometric errors at face value. To achieve $p(\chi^{2})=0.5$ we estimated that errors of 4 mas and 0$\fdg$3 would need to be added in quadrature to all separation and PA measurements, respectively. Alternatively, we could exclude a handful of outlier measurements (which have reasons for being suspect) to decrease the $\chi^{2}$ of the maximum likelihood solution to a reasonable value.

Five epochs of VLT/NaCo astrometry from Milli et al. (2014) account for 30% of the $\chi^{2}$ in the relative orbit fit. Chauvin et al. (2012) and Lagrange et al. (2018a) did not use any of these five epochs, even though they each could have used at least some. We therefore exclude all Milli et al. (2014) astrometry, which is contemporaneous with other available measurements. Likewise, Gemini/NICI astrometry from Nielsen et al. (2014) has three highly discrepant measurements (25% of the total $\chi^{2}$) that each were obtained on the same night as another measurement that is more consistent with the orbit fit. We exclude these three measurements as well, using 42 relative astrometry measurements in our final orbital analysis.

As shown by Brandt et al. (2018), simultaneous measurements of projected relative separation, host-star radial velocity, and host-star astrometric acceleration can provide a direct measurement of companion mass. In practice, observations of directly imaged companions are never truly simultaneous, although for very long long orbital periods on the order of centuries this can be a good approximation. The orbit of $\beta$ Pic b is of order decades (e.g., Chauvin et al., 2012; Nielsen et al., 2014), so more detailed analysis is needed to produce a companion mass from combining these three types of measurements. Our approach is described in detail in Brandt et al. (2018) and briefly here.

Posteriors of orbital parameters were determined using the parallel-tempering Markov chain Monte Carlo (PT-MCMC) ensemble sampler in emcee v2.1.0 (Foreman-Mackey et al., 2013) based on the algorithm described by Earl & Deem (2005). We ran 30 temperatures and 100 walkers fitting for eleven parameters, including the masses of the host star ($M_{\rm host}$) and planet ($M_{\rm comp}$). Eight others define the orbit, including the zero point of the system velocity (RVzero) and the intrinsic RV jitter ($\sigma_{\rm jit}$). We also included parallax ($\varpi$) as a fitted parameter with a Gaussian prior based on the measured value. For the initial step, we drew random values according to our priors across all valid parameter space, where for log-flat priors we used bounds of 0.3–3.0 $M_{\sun}$ in $M_{\rm host}$, 0.001–0.1 $M_{\sun}$ in $M_{\rm comp}$, 1–100 AU in $a$, and 0.3–300 m s*-1* in $\sigma_{\rm jit}$. We used $3\times 10^{5}$ steps in our PT-MCMC analysis, saving every 50th step of our chains. After ensuring that all walkers had stabilized in the mean and standard deviation of the posterior for each of the parameters we discarded all but the last $10^{3}$ samples as the burn-in portion yielding $10^{5}$ PT-MCMC samples across all walkers in the cold chain. Table 2 provides information on all our priors and posteriors, Figure 1 shows our orbit fit compared to the input measurements, and Figure 2 shows posteriors of astrophysically important parameters.

4. Discussion

Previous work has established key aspects of the orbit of $\beta$ Pic b, such as the viewing geometry of the nearly edge-on orbit and total system mass (e.g., Chauvin et al., 2012; Nielsen et al., 2014; Macintosh et al., 2014; Snellen et al., 2014; Wang et al., 2016). The reflex motion induced on $\beta$ Pic by $\beta$ Pic b is this same orbit scaled down by the mass ratio. Because of stellar proper motion, detecting this reflex motion and obtaining a dynamical mass for $\beta$ Pic b requires measuring nonlinear perturbations on the motion of $\beta$ Pic. In principle, RVs could determine a mass for $\beta$ Pic b, but the substantial RV jitter on such a young, active star ($\sigma_{\rm jit}=269\pm 6$ m s*-1*) hampers the measurement of the $\sim$100 m s*-1* expected semiamplitude of the planet.

HGCA reports deviations from constant proper motion of only 1–2$\sigma$ for $\beta$ Pic. Our joint fit of astrometry and RVs yields a mass posterior of $13\pm 3$ $M_{\rm Jup}$ (23% uncertainty) for $\beta$ Pic b. This is not as precise as the value of $11\pm 2$ $M_{\rm Jup}$ (18% uncertainty) from Snellen & Brown (2018) because we adopted cross-calibrated Hipparcos and Gaia DR2 astrometry, which Brandt (2018) found requires error inflation of reported astrometric errors in both catalogs. Moreover, our analysis does not assume a host-star mass, although it broadly supports previous assumptions of 1.75 $M_{\sun}$ with a remarkably precise model-independent mass of $1.84\pm 0.05$ $M_{\sun}$. Our mass determination for $\beta$ Pic b is chiefly driven by the small offset between the Hipparcos proper motion and the Hipparcos-to-Gaia positional difference, as was the case in Snellen & Brown (2018), because the uncertainty in the Gaia DR2 proper motion is very large due to $\beta$ Pic being saturated in Gaia.

Combining our mass for $\beta$ Pic b with the luminosity of \log(\mbox{L_{\rm bol}}/\mbox{L_{\sun}})=-3.78\pm 0.03 dex determined by Morzinski et al. (2015) we calculate upper and lower limits on the substellar cooling age from the hot-start evolutionary models of Saumon & Marley (2008). We use the same method described in Dupuy & Liu (2017), with a uniform prior in age, our orbit posterior as the prior on mass, and rejection-sampling on $L_{\rm bol}$ to select models consistent with $\beta$ Pic b. The posterior on the age is wide, as expected given the low-precision mass. The 3$\sigma$ confidence interval on the cooling age ranges from 7–65 Myr. Combining the 7 Myr lower limit with external age information for $\beta$ Pic and its eponymous young moving group directly determines the amount of time that could have elapsed between the formation of the host star and planet. Adopting the $\beta$ Pic moving group age of $22\pm 6$ Myr from Shkolnik et al. (2017), we thus find an upper limit of $15\pm 6$ Myr in the difference between the times of formation (a.k.a. $t=0$) for $\beta$ Pic b and its host star. This is not particularly constraining on theory, but improved precision in the mass of $\beta$ Pic b in the future will result in stronger tests of the timescale of giant planet formation.

Our dynamical mass of $13\pm 3$ $M_{\rm Jup}$ for $\beta$ Pic b is broadly consistent with hot-start formation models, as these predict a mass of $13.0^{+0.4}_{-0.3}$ $M_{\rm Jup}$ at $22\pm 6$ Myr (Dupuy et al., 2018). The high-mass end of our posterior is also consistent with warm-start models. Our 2$\sigma$ upper limit on the mass of $\beta$ Pic b is 19.5 $M_{\rm Jup}$ (0.019 $M_{\sun}$). Interpolating hot-start evolutionary tracks from Saumon & Marley (2008), an object of this mass should have a luminosity of \log(\mbox{L_{\rm bol}}/\mbox{L_{\sun}})=-3.1 dex at an age of 22 Myr. The actual luminosity of $\beta$ Pic b is 0.7 dex (1.7 mag) fainter than this. This decrement corresponds to the intermediate range of 10-$M_{\rm Jup}$ warm-start models from Spiegel & Burrows (2012, see their Figure 9) and is highly inconsistent with cold-start models that are $\approx$5 mag fainter than hot-start tracks at 22 Myr. However, warm-start models would need to be computed beyond 10 $M_{\rm Jup}$ for a more accurate appraisal.

Our relative orbit for $\beta$ Pic b is consistent with past work within the uncertainties, but our posteriors are notably lacking any near-circular orbits, with $e<0.1$ excluded at $>2\sigma$. Previous work was generally consistent with eccentricities up to 0.1–0.2 but preferred more circular orbits, unlike our orbit fit ($e=0.24\pm 0.06$). Based on tests using various subsets of the relative astrometry, we find that our results are simply the consequence of combining all published measurements in a joint fit. Recent results from Wang et al. (2016) did not have access to the VLT/SPHERE measurements, and Lagrange et al. (2018a) used only VLT astrometry in their analysis. We note that our choice to exclude eight relative astrometry outliers out of 50 measurements does not significantly impact this result. We ran an identical PT-MCMC using all 50 measurements, with errors of 4 mas and 0$\fdg$3 added in quadrature to all separations and PAs in order to make a reasonable $\chi^{2}$. All parameter posteriors were very similar, including a slightly higher eccentricity of $0.28\pm 0.06$.

The strong preference of our fit for non-circular orbits has implications for the origin of $\beta$ Pic b and its history of dynamical interactions with the disk. While the focus of the literature has been on the planet–disk interaction as a means to explain the disk warp (e.g., Dawson et al., 2011), the eccentricity may help explain other observations (Apai et al., 2015; Millar-Blanchaer et al., 2015; Wang et al., 2016). The higher eccentricity for the planet ($e\sim 0.2$) is also consistent with the exo-comet hypothesis put forth to explain the occasional absorption features in the host star’s spectrum. Thébault & Beust (2001) note that the frequency of observed events is well explained by the excitation of cometary bodies in a 3:1 resonance with a massive perturber at roughly 10 AU. However, the location from which these comets would be launched lies somewhat inside the inner edge of the disk as fit by Millar-Blanchaer et al. (2015). Simulations of the planet–disk interaction at lower values of $e$ have not successfully explained the observed inner edge and disk scale height (Millar-Blanchaer et al., 2015; Nesvold & Kuchner, 2015). These authors also note that a low eccentricity ($e<0.1$) is unable to account for the observed northeast–southwest brightness asymmetry in the disk. Although several authors have proposed a possible unseen second planet to explain these features (Apai et al., 2015; Millar-Blanchaer et al., 2015; Nesvold & Kuchner, 2015; Wang et al., 2016), further modeling with $\beta$ Pic b alone at a higher eccentricity may be warranted.

While the higher eccentricity for $\beta$ Pic b may help explain some present-day disk observations, it is consistent with a wide range of past formation scenarios. Giant planets are thought to form on relatively circular orbits due to efficient damping in the natal protoplanetary disk (Armitage, 2011), but many mechanisms can subsequently pump their eccentricities. High eccentricities can easily be generated by secular, resonant, or scattering interactions with a massive perturber in the form of another planet or nearby stars (Holman et al., 1997; Laughlin & Adams, 1998; Rasio & Ford, 1996; Ford & Rasio, 2008). The perturber responsible for $\beta$ Pic b’s eccentricity need not remain in the system and be observable today; it could have been ejected by a strong scattering event. As noted above, the presence of a second planet is favored in some models to explain disk structures. Detailed analysis of long-term RV monitoring excludes much of the parameter space for additional planets (Lagrange et al., 2018b), and while absolute astrometry can potentially rule out more planets (Kervella et al., 2018), the fact that $\beta$ Pic b is only marginally detected in current observations complicates the interpretation of additional astrometric signals due to more planets.

A second massive planet is not, however, required to generate an eccentricity as high as $e=0.2$–0.3. In principle, migration of a very massive planet like $\beta$ Pic b through a massive gas disk (e.g., Papaloizou et al., 2001; Ragusa et al., 2018) or even a planetesimal disk (Murray et al., 1998) can generate substantial eccentricity growth. While more modest-mass planets have their eccentricities damped by the disk, planets with masses \gtrsim 10\,\mbox{M_{\rm Jup}} can have their eccentricities pumped through interaction with the 3:1 outer Lindblad resonance: the planet excites eccentricity in the disk, which back-reacts to excite eccentricity in the planet (Kley & Nelson, 2012). Papaloizou et al. (2001) explicitly predict that a massive eccentric planet inside a disk cavity is a natural outcome of this process. A key implication of this mechanism is that $\beta$ Pic b formed exterior to its current orbit.

Given the still limited observational coverage of the $\approx$30-year orbit of $\beta$ Pic b, and a particular lack of data on the northeastern side, its eccentricity is still relatively uncertain. In Table 3, we provide predicted astrometry and RVs for three representative eccentricities from our PT-MCMC posterior ($e=0.1$, 0.2, and 0.3). In the near term, the RV of $\beta$ Pic b is the most discriminating between different eccentricities, but after a few years separation measurements will cleanly define the orbit. Lower eccentricity orbits predict smaller separations in the next decade and a more imminent turnaround toward decreasing separation.

We thank the referee for a thoughtful and timely review. This work has made use of data from the European Space Agency mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. T.J.D. acknowledges research support from Gemini Observatory. T.D.B. gratefully acknowledges support from the Heising-Simons foundation and from NASA under grant #80NSSC18K0439. Our research has employed NASA ADS; SIMBAD; VizieR; and J. R. A. Davenport’s IDL implementation of the cubehelix color scheme (Green, 2011).