<figure><img src="image_1.png" /></figure>

HAL
open science

# The Gemini Planet Imager Exoplanet Survey: Dynamical Mass of the Exoplanet β Pictoris b from Combined Direct Imaging and Astrometry

Eric L. Nielsen, Robert J. de Rosa, Jason J. Wang, Johannes Sahlmann, Paul Kalas, Gaspard Duchêne, Julien Rameau, Mark S. Marley, Didier Saumon, Bruce Macintosh, et al.

► To cite this version:

Eric L. Nielsen, Robert J. de Rosa, Jason J. Wang, Johannes Sahlmann, Paul Kalas, et al.. The Gemini Planet Imager Exoplanet Survey: Dynamical Mass of the Exoplanet β Pictoris b from Combined Direct Imaging and Astrometry. The Astronomical Journal, 2020, 159, 10.3847/1538-3881/ab5b92. insu-03705224

HAL Id: insu-03705224

https://insu.hal.science/insu-03705224v1

Submitted on 4 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# The Gemini Planet Imager Exoplanet Survey: Dynamical Mass of the Exoplanet *β* Pictoris b from Combined Direct Imaging and Astrometry

Eric L. Nielsen¹, Robert J. De Rosa¹, Jason J. Wang²,<fnref n="35" />, Johannes Sahlmann³, Paul Kalas⁴,⁵,⁶, Gaspard Duchêne⁴,⁷, Julien Rameau⁷,⁸, Mark S. Marley⁹, Didier Saumon¹⁰, Bruce Macintosh¹, Maxwell A. Millar-Blanchaer¹¹,<fnref n="36" />, Meiji M. Nguyen⁴, S. Mark Ammons¹², Vanessa P. Bailey¹¹, Travis Barman¹³, Joanna Bulger¹⁴,¹⁵, Jeffrey Chilcote¹⁶, Tara Cotten¹⁷, Rene Doyon⁸, Thomas M. Esposito⁴, Michael P. Fitzgerald¹⁸, Katherine B. Follette¹⁹, Benjamin L. Gerard²⁰,²¹, Stephen J. Goodsell²², James R. Graham⁴, Alexandra Z. Greenbaum²³, Pascale Hibon²⁴, Li-Wei Hung²⁵, Patrick Ingraham²⁶, Quinn Konopacky²⁷, James E. Larkin¹⁸, Jérôme Maire²⁷, Franck Marchis⁵, Christian Marois²⁰,²¹, Stanimir Metchev²⁸,²⁹, Rebecca Oppenheimer³⁰, David Palmer¹², Jennifer Patience³¹, Marshall Perrin³, Lisa Poyneer¹², Laurent Pueyo³, Abhijith Rajan³, Fredrik T. Rantakyrö²⁴, Jean-Baptiste Ruffio¹, Dmitry Savransky³², Adam C. Schneider³¹, Anand Sivaramakrishnan³, Inseok Song¹⁷, Remi Soummer³, Sandrine Thomas²⁶, J. Kent Wallace¹¹, Kimberly Ward-Duong¹⁹, Sloane Wiktorowicz³³, and Schuyler Wolff³⁴

¹ Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, Stanford, CA 94305, USA; enielsen@stanford.edu

² Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA

³ Space Telescope Science Institute, Baltimore, MD 21218, USA

⁴ Department of Astronomy, University of California, Berkeley, CA 94720, USA

⁵ SETI Institute, Carl Sagan Center, 189 Bernardo Ave., Mountain View, CA 94043, USA

⁶ Institute of Astrophysics, FORTH, GR-71110 Heraklion, Greece

⁷ Univ. Grenoble Alpes/CNRS, IPAG, F-38000 Grenoble, France

⁸ Institut de Recherche sur les Exoplanètes, Département de Physique, Université de Montréal, Montréal QC, H3C 3J7, Canada

⁹ NASA Ames Research Center, MS 245-3, Mountain View, CA 94035, USA

¹⁰ Los Alamos National Laboratory, P.O. Box 1663, Los Alamos, NM 87545, USA

¹¹ Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA

¹² Lawrence Livermore National Laboratory, Livermore, CA 94551, USA

¹³ Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721, USA

¹⁴ Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA

¹⁵ Subaru Telescope, NAOJ, 650 North A'ohoku Place, Hilo, HI 96720, USA

¹⁶ Department of Physics, University of Notre Dame, 225 Nieuwland Science Hall, Notre Dame, IN 46556, USA

¹⁷ Department of Physics and Astronomy, University of Georgia, Athens, GA 30602, USA

¹⁸ Department of Physics & Astronomy, University of California, Los Angeles, CA 90095, USA

¹⁹ Physics and Astronomy Department, Amherst College, 21 Merrill Science Drive, Amherst, MA 01002, USA

²⁰ University of Victoria, 3800 Finnerty Rd, Victoria, BC, V8P 5C2, Canada

²¹ National Research Council of Canada Herzberg, 5071 West Saanich Rd, Victoria, BC, V9E 2E7, Canada

²² Gemini Observatory, 670 N. A'ohoku Place, Hilo, HI 96720, USA

²³ Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA

²⁴ Gemini Observatory, Casilla 603, La Serena, Chile

²⁵ Natural Sounds and Night Skies Division, National Park Service, Fort Collins, CO 80525, USA

²⁶ Large Synoptic Survey Telescope, 950N Cherry Ave., Tucson, AZ 85719, USA

²⁷ Center for Astrophysics and Space Science, University of California San Diego, La Jolla, CA 92093, USA

²⁸ Department of Physics and Astronomy, Centre for Planetary Science and Exploration, The University of Western Ontario, London, ON N6A 3K7, Canada

²⁹ Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794-3800, USA

³⁰ Department of Astrophysics, American Museum of Natural History, New York, NY 10024, USA

³¹ School of Earth and Space Exploration, Arizona State University, P.O. Box 871404, Tempe, AZ 85287, USA

³² Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, NY 14853, USA

³³ Department of Astronomy, UC Santa Cruz, 1156 High St., Santa Cruz, CA 95064, USA

³⁴ Leiden Observatory, Leiden University, 2300 RA Leiden, The Netherlands

Received 2019 April 12; revised 2019 November 20; accepted 2019 November 23; published 2020 January 28

## Abstract

We present new observations of the planet *β* Pictoris b from 2018 with the Gemini Planet Imager (GPI), the first GPI observations following conjunction. Based on these new measurements, we perform a joint orbit fit to the available relative astrometry from ground-based imaging, the *Hipparcos* Intermediate Astrometric Data (IAD), and the *Gaia* DR2 position, and demonstrate how to incorporate the IAD into direct imaging orbit fits. We find a mass consistent with predictions of hot-start evolutionary models and previous works following similar methods, though with larger uncertainties: $12.8^{+5.3}_{-3.2} M_{Jup}$. Our eccentricity determination of $0.12^{+0.04}_{-0.03}$ disfavors circular orbits. We consider orbit fits to several different imaging data sets, and find generally similar posteriors on the mass for each combination of imaging data. Our analysis underscores the importance of performing joint fits to the absolute and relative astrometry simultaneously, given the strong covariance between orbital elements. Time of conjunction is well-constrained within 2.8 days of 2017 September 13, with the star behind the planet's Hill sphere between 2017 April 11 and 2018 February 16 (±18 days). Following the recent radial velocity detection of a second planet in the system, *β* Pic c, we perform additional two-planet fits combining relative astrometry, absolute astrometry, and

<footnote n="35">³⁵ 51 Pegasi b Fellow.</footnote>

<footnote n="36">³⁶ NASA Hubble Fellow.</footnote>
stellar radial velocities. These joint fits find a significantly smaller mass ($8.0 \pm 2.6 M_{Jup}$) for the imaged planet $\beta$ Pic b, in a somewhat more circular orbit. We expect future ground-based observations to further constrain the visual orbit and mass of the planet in advance of the release of Gaia DR4.

*Unified Astronomy Thesaurus concepts: Astrometry (80); Direct imaging (387); Exoplanets (498); Coronographic imaging (313); Orbit determination (1175)*

# 1. Introduction

Masses of exoplanets detected by the radial velocity (RV) method can be directly measured to within $sin(i)$, as can the mass ratio between microlensing planets and their parent star, and masses can be inferred for transiting planet systems by modeling transit timing variations. The masses of directly imaged planets, however, must be inferred from evolutionary models if only imaging data are available. These models predict the mass of the planet as a function of age of the system and luminosity of the planet. While the COND models (Baraffe et al. 2003) have been consistent with upper limits on directly imaged planet masses (Lagrange et al. 2012; Wang et al. 2018), direct measurements of the mass allow for a more robust testing of the models. Giant planets are most easily imaged around young stars (~100 Myr), which tend to be too active for precise RV measurements; e.g., Lagrange et al. (2012) describe searching for a ~10 $m s^{-1}$ signal in RV data with a ~3 $km s^{-1}$ peak-to-peak RV variation. Astrometry, however, is less affected by stellar activity, and represents a way forward to determining the dynamical mass of these planets from stellar reflex motion.

In particular, the second Gaia data release (DR2) gives independent measurements of ~2016 position and proper motions for one billion stars. Recently, Snellen & Brown (2018) combined an orbit fit to direct imaging data by Wang et al. (2016) with Hipparcos Intermediate Astrometric Data (IAD) and Gaia positions for the planet $\beta$ Pictoris b. This combination of the orbital period from imaging, with absolute positions in ~1991 and ~2016 resulted in a measurement of the planet mass of $11 \pm 2 M_{Jup}$. A similar analysis was undertaken by Dupuy et al. (2019) earlier this year.

$\beta$ Pic is a young, nearby ($d = 19.44$ pc), intermediate-mass (~1.8 $M_{\odot}$) star that hosts a bright edge-on debris disk (Smith & Terrile 1984; Kalas & Jewitt 1995; Wahhaj et al. 2003; Weinberger et al. 2003; Golimowski et al. 2006; Nielsen et al. 2014). It is part of the $\beta$ Pic moving group (Barrado y Navascués et al. 1999; Zuckerman et al. 2001; Binks & Jeffries 2014; Bell et al. 2015), which sets the age of the star to 26 ± 3 Myr (Nielsen et al. 2016). $\beta$ Pic b was one of the first directly imaged exoplanets, first observed on the northeast side of the star in 2003 (Lagrange et al. 2009), before being confirmed after it passed behind the star to the southwest side (Lagrange et al. 2010). Subsequent observations allowed the orbit of the planet to be determined to increasing accuracy (Currie et al. 2011; Chauvin et al. 2012; Macintosh et al. 2014; Nielsen et al. 2014; Millar-Blanchaer et al. 2015; Wang et al. 2016). The planet's orbital plane has been found to be very similar to the plane of the disk, and though a transit-like event was observed in 1981 (Lecavelier Des Etangs & Vidal-Madjar 2009), additional relative astrometry has ruled out the possibility of the planet itself transiting, although the planet's Hill sphere does pass in front of the star (e.g., Millar-Blanchaer et al. 2015; Wang et al. 2016).

The Gemini Planet Imager (GPI; Macintosh et al. 2014) is an extreme adaptive optics system on the Gemini South 8 m

telescope optimized for detecting self-luminous giant exoplanets. $\beta$ Pic was observed multiple times by GPI since 2013, tracking the orbit of the planet as it moved closer to the star (Wang et al. 2016). Here, we present new observations from GPI in 2018, following conjunction, and a joint fit of the imaging data and Hipparcos and Gaia astrometry, along with an estimate of the mass of the planet.

# 2. Observations

## 2.1. New GPI data

$\beta$ Pic was observed in 2018 after a hiatus in which the planet passed too close to the star in angular projection (~0.15). Due to the close angular separation of the planet and the star, we chose to observe $\beta$ Pic b in the J band in order to maximize sensitivity at ~150 mas radius while maintaining a favorable flux ratio of the planet. In this paper, we present two epochs of GPI J-band integral field spectroscopy observations of the planet. The first epoch was taken on 2018 September 21 between 8:42 and 10:02 UT. After discarding frames in which the AO loops opened, we obtained a total of 59 exposures with integration times of 60 s. A total of 36.8 of field rotation was obtained for angular differential imaging (ADI; Marois et al. 2006a). The second epoch was taken on 2018 November 18 between 5:51 and 9:13 UT with a total of 145 exposures, each of which is comprised of four coadded 14.5 s frames. These observations were better timed, and a total of 96.9 of field rotation was obtained for ADI.

The data were first reduced using the automated GPIES data reduction pipeline (Wang et al. 2018), with one notable exception. During the night of 2018 September 21, GPI was not able to access the Gemini Facility Calibration Unit and could not obtain an Argon arc lamp snapshot before each observation sequence to correct for instrument flexure (Wolff et al. 2014). For the $\beta$ Pic observations, we corrected instrument flexure manually through visual inspection. This did not significantly impact the spatial image reconstruction of the three pixel box extraction algorithm used in the GPI Data Reduction Pipeline (DRP; Perrin et al. 2014, 2016), but it likely affected our spectral accuracy. However, for the purpose of astrometry, we collapse the spectral data cubes into a broadband image, so the impact on astrometry is minimal. In both epochs, we used the satellite spots, four fiducial diffraction spots centered on the location of the star (Marois et al. 2006b; Sivaramakrishnan & Oppenheimer 2006), to locate the star behind the coronagraph in each wavelength slice of each spectral data cube (Wang et al. 2014). The stellar point-spread function (PSF) was then subtracted out using *pyKLIP* (Wang et al. 2015), which uses principal component analysis (Sommer et al. 2012; Pueyo et al. 2015) constructed from images taken at other times (ADI) and wavelengths (spectral differential imaging (SDI); Sparks & Ford 2002). The reductions of the two epochs (shown in Figure 1) use 20 principal components to model and subtract out the stellar PSF, and are averaged over time and wavelength. We estimated
<figure><img src="image_4.png" /><figcaption>**Figure 1.** GPI images of β Pic b processed with the automatic GPIES pipeline. Images are rotated north-up/east-left, and have not been flux-calibrated. Colors are presented on a linear scale, and the white arrows point to the location of the planet.</figcaption></figure>

signal-to-noise ratios of 4.5 and 11.7 for the September and November data sets, respectively.

To measure the position of β Pic b in each data set, we follow the same technique that was outlined in Wang et al. (2016), where the signal of the planet is forward modeled through the data reduction process and the forward model is then fit to the data. In these reductions, we found it was optimal to discard frames from the sequences, due to varying image quality. For both data sets, we ordered data cubes by the contrast in each single data cube at 250 mas. For the September 21st epoch, we only used the best 40 data cubes, resulting in a total integration time of 40 minutes. For the November 18th epoch, we used the best 120 frames, resulting in a total integration time of 116 minutes.

We then used the astrometry modules in *pyKLIP* to run a stellar PSF subtraction that simultaneously forward models the PSF of the planet. For both epochs, we built 15 principal components from the 150 more correlated reference PSFs, where the reference PSFs are drawn from frames at other wavelengths and times where β Pic b moved at least one pixel in the image, due to a combination of ADI and SDI. For the September 21st epoch, we broke up the image between 6.5 and 25.6 pixels from the star into three concentric annuli of 4.0, 6.7, and 8.4 pixels in width, respectively. We then broke each annuli into four sectors, and ran our stellar PSF subtraction and forward modeling on each sector. For the November 18th epoch, we only used one annulus, centered on the star, with an inner radius of 6.5 pixels and an outer radius of 19.2 pixels. We did not split up this annulus. The annuli geometries were defined by the focal plane mask and the edge of the field of view. The planet's position in the data was then fit over a box nine pixels wide, centered on the estimated location of the planet. In this box were pixels that fell inside the focal plane mask and were not included in our reduction. We did not consider these pixels in the fit, thus reducing the number of data points by a few. The fit was done using the Bayesian framework described in Wang et al. (2016), where we used a Gaussian process to model the correlated speckle noise present in the data. Due to the close separation of the planet in these two epochs, we did not trust the assumption of Gaussian noise used in our Bayesian framework when estimating uncertainties on the planet's location. To empirically quantify this and any residual biases in the forward model, we injected simulated planets into the data sets with a spectrum from a model fit to β Pic b's spectrum at the same separation as β Pic b, but at position angles that are at least three full widths at half maximum apart from the measured position of the planet. We

injected one simulated planet at a time, measured its astrometry, and compared it to the true position at which we injected it. We found a scatter in the position of 0.3 pixels for the September 21st epoch, and a scatter in the position of 0.13 pixels for the position of the November 18th epoch. We found the average measured astrometry of the simulated planets was biased by <0.02 pixels, so we conclude that fitting biases are negligible. We use the scatter in the simulated planet positions as the uncertainty in the position of β Pic b. To obtain relative astrometry of the host star, we assumed a star centering precision of 0.05 pixels (Wang et al. 2014), a plate scale value of 14.161 ± 0.021 mas/pixel, and a residual north angle correction of 0°.45 ± 0°.11 (De Rosa et al. 2019). The relative astrometry is reported in Table 1.

## 2.2. Previously Published Data Sets

As in Nielsen et al. (2014) and Wang et al. (2016), we compile relative astrometry of β Pic b from the literature to extend the time baseline. Chauvin et al. (2012) presented nine epochs of data from VLT/NACO, including the two initial discovery epochs of 2003 (Lagrange et al. 2009) and 2009 (Lagrange et al. 2010), up until 2011. An additional seven epochs of data from 2009 to 2012 were reported from Gemini-South/NICI by Nielsen et al. (2014), as well as two 2012 epochs from *Magellan/MagAO* (Morzinski et al. 2015). Twelve epochs of Gemini-South/GPI data were presented by Wang et al. (2016), running from 2013 to 2016. An additional attempt was made to observe β Pic b with GPI on UT 2016 November 18; however, due to its proximity to the host star and the poor seeing that night, the planet was not detected in this data set. Recently, Lagrange et al. (2019a) published eleven epochs of relative astrometry from VLT/SPHERE between 2014 and 2016, as well as an epoch from 2018 September 17 when the planet reappeared on the northeast side of the star.

Due to the timing issue and a change in the astrometric calibration (De Rosa et al. 2019), we also recomputed the astrometry of the epochs published in Wang et al. (2016) using the same reduction parameters as the previous work. The parallactic angles in each data set were recomputed with the correct time in the header following the procedure outlined in De Rosa et al. (2019). We also used the new plate scale value of 14.161 ± 0.021 mas/pixel and varying residual north angle correction from De Rosa et al. (2019). We used a residual north angle of 0°.23 ± 0°.11 for the 2013 epochs, 0°.17 ± 0°.14 for the 2014 November 8 and 2015 April 2 epochs, and 0°.21 ± 0°.23 for the remaining 2015 and 2016 epochs. The recomputed astrometry is listed in Table 1. The most significant change to the astrometry presented here, compared to Wang et al. (2016), is the change in assumed north angle, from -0°.2 to approximately +0°.2, shifting all position angles to larger values by ~0°.4. Additionally, we include an additional epoch from 2015 January 24, which had been initially rejected in Wang et al. (2016) due to artifacts at the location of the planet. With the rereduction, the artifacts are no longer visible, so we include this epoch in our final data set.

# 3. Orbit Fitting

## 3.1. Hipparcos Intermediate Astrometric Data

The *Hipparcos* mission performed detailed astrometric monitoring of bright stars, with the majority of stars (including
<table><caption>**Table 1**
Relative Astrometry of β Pic b</caption><thead><tr><th>Epoch</th><th>Sep (")</th><th>PA (deg)</th><th>Instrument</th><th>References</th></tr></thead><tbody><tr><td>2008 Nov 11</td><td>0.210 ± 0.027</td><td>211.49 ± 1.9</td><td>VLT/NACO</td><td>Currie et al. (2011)</td></tr><tr><td>2003 Nov 10</td><td>0.413 ± 0.022</td><td>34 ± 4</td><td>VLT/NACO</td><td>Chauvin et al. (2012)</td></tr><tr><td>2009 Oct 25</td><td>0.299 ± 0.014</td><td>211 ± 3</td><td>VLT/NACO</td><td>Chauvin et al. (2012)</td></tr><tr><td>2009 Dec 29</td><td>0.306 ± 0.009</td><td>212.1 ± 1.7</td><td>VLT/NACO</td><td>Chauvin et al. (2012)</td></tr><tr><td>2010 Apr 10</td><td>0.346 ± 0.007</td><td>209.9 ± 1.2</td><td>VLT/NACO</td><td>Chauvin et al. (2012)</td></tr><tr><td>2010 Sep 28</td><td>0.383 ± 0.011</td><td>210.3 ± 1.7</td><td>VLT/NACO</td><td>Chauvin et al. (2012)</td></tr><tr><td>2010 Nov 16</td><td>0.387 ± 0.008</td><td>212.4 ± 1.4</td><td>VLT/NACO</td><td>Chauvin et al. (2012)</td></tr><tr><td>2010 Nov 17</td><td>0.390 ± 0.013</td><td>212 ± 2</td><td>VLT/NACO</td><td>Chauvin et al. (2012)</td></tr><tr><td>2011 Feb 1</td><td>0.408 ± 0.009</td><td>211.1 ± 1.5</td><td>VLT/NACO</td><td>Chauvin et al. (2012)</td></tr><tr><td>2011 Mar 26</td><td>0.426 ± 0.013</td><td>210.1 ± 1.8</td><td>VLT/NACO</td><td>Chauvin et al. (2012)</td></tr><tr><td>2009 Dec 3</td><td>0.339 ± 0.010</td><td>209.2 ± 1.7</td><td>Gemini-South/NICI</td><td>Nielsen et al. (2014)</td></tr><tr><td>2009 Dec 3</td><td>0.323 ± 0.010</td><td>209.3 ± 1.8</td><td>Gemini-South/NICI</td><td>Nielsen et al. (2014)</td></tr><tr><td>2010 Dec 25</td><td>0.407 ± 0.005</td><td>212.8 ± 1.4</td><td>Gemini-South/NICI</td><td>Nielsen et al. (2014)</td></tr><tr><td>2011 Oct 20</td><td>0.452 ± 0.003</td><td>211.6 ± 0.4</td><td>Gemini-South/NICI</td><td>Nielsen et al. (2014)</td></tr><tr><td>2011 Oct 20</td><td>0.455 ± 0.005</td><td>211.9 ± 0.6</td><td>Gemini-South/NICI</td><td>Nielsen et al. (2014)</td></tr><tr><td>2012 Mar 29</td><td>0.447 ± 0.003</td><td>210.8 ± 0.4</td><td>Gemini-South/NICI</td><td>Nielsen et al. (2014)</td></tr><tr><td>2012 Mar 29</td><td>0.448 ± 0.005</td><td>211.8 ± 0.6</td><td>Gemini-South/NICI</td><td>Nielsen et al. (2014)</td></tr><tr><td>2012 Dec 2</td><td>0.461 ± 0.014</td><td>211.9 ± 1.2</td><td>Magellan/MagAO</td><td>Nielsen et al. (2014)</td></tr><tr><td>2012 Dec 4</td><td>0.470 ± 0.010</td><td>212.0 ± 1.2</td><td>Magellan/MagAO</td><td>Nielsen et al. (2014)</td></tr><tr><td>2013 Nov 16</td><td>0.4308 ± 0.0015</td><td>212.43 ± 0.17</td><td>Gemini-South/GPI</td><td>This Work<sup>a</sup></td></tr><tr><td>2013 Nov 16</td><td>0.4291 ± 0.0010</td><td>212.58 ± 0.15</td><td>Gemini-South/GPI</td><td>This Work<sup>a</sup></td></tr><tr><td>2013 Nov 18</td><td>0.4302 ± 0.0010</td><td>212.46 ± 0.15</td><td>Gemini-South/GPI</td><td>This Work<sup>a</sup></td></tr><tr><td>2013 Dec 10</td><td>0.4255 ± 0.0010</td><td>212.51 ± 0.15</td><td>Gemini-South/GPI</td><td>This Work<sup>a</sup></td></tr><tr><td>2013 Dec 10</td><td>0.4244 ± 0.0010</td><td>212.85 ± 0.15</td><td>Gemini-South/GPI</td><td>This Work<sup>a</sup></td></tr><tr><td>2013 Dec 11</td><td>0.4253 ± 0.0010</td><td>212.47 ± 0.16</td><td>Gemini-South/GPI</td><td>This Work<sup>a</sup></td></tr><tr><td>2014 Nov 8</td><td>0.3562 ± 0.0010</td><td>213.02 ± 0.19</td><td>Gemini-South/GPI</td><td>This Work<sup>a</sup></td></tr><tr><td>2015 Apr 2</td><td>0.3173 ± 0.0009</td><td>213.13 ± 0.20</td><td>Gemini-South/GPI</td><td>This Work<sup>a</sup></td></tr><tr><td>2015 Nov 6</td><td>0.2505 ± 0.0015</td><td>214.14 ± 0.34</td><td>Gemini-South/GPI</td><td>This Work<sup>a</sup></td></tr><tr><td>2015 Dec 5</td><td>0.2402 ± 0.0011</td><td>213.58 ± 0.34</td><td>Gemini-South/GPI</td><td>This Work<sup>a</sup></td></tr><tr><td>2015 Dec 22</td><td>0.2345 ± 0.0010</td><td>213.81 ± 0.30</td><td>Gemini-South/GPI</td><td>This Work<sup>a</sup></td></tr><tr><td>2016 Jan 21</td><td>0.2226 ± 0.0021</td><td>214.84 ± 0.44</td><td>Gemini-South/GPI</td><td>This Work<sup>a</sup></td></tr><tr><td>2014 Dec 8</td><td>0.35051 ± 0.00320</td><td>212.60 ± 0.66</td><td>VLT/SPHERE</td><td>Lagrange et al. (2019a)</td></tr><tr><td>2015 May 5</td><td>0.33242 ± 0.00170</td><td>212.58 ± 0.35</td><td>VLT/SPHERE</td><td>Lagrange et al. (2019a)</td></tr><tr><td>2015 Oct 1</td><td>0.26202 ± 0.00178</td><td>213.02 ± 0.48</td><td>VLT/SPHERE</td><td>Lagrange et al. (2019a)</td></tr><tr><td>2015 Nov 30</td><td>0.24205 ± 0.00251</td><td>213.30 ± 0.74</td><td>VLT/SPHERE</td><td>Lagrange et al. (2019a)</td></tr><tr><td>2015 Dec 26</td><td>0.23484 ± 0.00180</td><td>213.79 ± 0.51</td><td>VLT/SPHERE</td><td>Lagrange et al. (2019a)</td></tr><tr><td>2016 Jan 20</td><td>0.22723 ± 0.00155</td><td>213.15 ± 0.46</td><td>VLT/SPHERE</td><td>Lagrange et al. (2019a)</td></tr><tr><td>2016 Mar 26</td><td>0.20366 ± 0.00142</td><td>213.90 ± 0.46</td><td>VLT/SPHERE</td><td>Lagrange et al. (2019a)</td></tr><tr><td>2016 Apr 16</td><td>0.19749 ± 0.00236</td><td>213.88 ± 0.83</td><td>VLT/SPHERE</td><td>Lagrange et al. (2019a)</td></tr><tr><td>2016 Sep 16</td><td>0.14236 ± 0.00234</td><td>214.62 ± 1.10</td><td>VLT/SPHERE</td><td>Lagrange et al. (2019a)</td></tr><tr><td>2016 Oct 14</td><td>0.13450 ± 0.00246</td><td>215.50 ± 1.22</td><td>VLT/SPHERE</td><td>Lagrange et al. (2019a)</td></tr><tr><td>2016 Nov 18</td><td>0.12712 ± 0.00644</td><td>215.80 ± 3.37</td><td>VLT/SPHERE</td><td>Lagrange et al. (2019a)</td></tr><tr><td>2018 Sep 17</td><td>0.14046 ± 0.00312</td><td>29.71 ± 1.67</td><td>VLT/SPHERE</td><td>Lagrange et al. (2019a)</td></tr><tr><td>2015 Jan 24</td><td>0.3355 ± 0.0009</td><td>212.88 ± 0.20</td><td>Gemini-South/GPI</td><td>This Work</td></tr><tr><td>2018 Sep 21</td><td>0.1419 ± 0.0053</td><td>28.16 ± 1.82</td><td>Gemini-South/GPI</td><td>This Work</td></tr><tr><td>2018 Nov 18</td><td>0.1645 ± 0.0018</td><td>28.64 ± 0.70</td><td>Gemini-South/GPI</td><td>This Work</td></tr></tbody></table>

**Note.**

β Pic) being fit by a five-parameter solution, R.A. and decl. of the star (as would be observed from the solar system barycenter) at a reference epoch of 1991.25, parallax, and proper motion in R.A. and decl. Individual measurements were made of each star along a one-dimensional scan referred to as the abscissa, with no published constraints in the direction perpendicular to the scan direction. The direction of the scan changed from orbit to orbit as the satellite surveyed the sky, allowing a two-dimensional motion to be reconstructed from a series of one-dimensional measurements. van Leeuwen (2007a)

provides IAD from the rereduction of the *Hipparcos* data in the form of a DVD-ROM attached to the book, including scan directions, residuals from the fit, and errors on the measurement, for each epoch of data.

While the IAD do not contain the abscissa measurements themselves, the measurements can be reconstructed from these values. We extract from the van Leeuwen (2007b) IAD the epoch of the orbit in decimal years (*t*), scan direction (sin(*φ*) and cos(*φ*)), residual to the best fit (*R*), and error on the original measurement (*ε*). This is combined with the best-fitting

<footnote><sup>a</sup> These epochs originally appeared in Wang et al. (2016), but have been recomputed here as a result of changes in the GPI pipeline—most noticeably the assumed north angle.</footnote>
solution from the van Leeuwen (2007a) catalog for the star, which provides the five astrometric parameters, $\alpha_0$, $\delta_0$, $\pi$, $\mu_{\alpha^*}$, $\mu_\delta$: the R.A. and decl. at the *Hipparcos* reference epoch of 1991.25 in degrees, the parallax in mas, and the proper motion in R.A. and decl. in mas yr$^{-1}$. The notation $\mu_{\alpha^*}$ indicates offsets and velocities in R.A. are multiplied by $\cos \delta_0$, to prevent a constant factor between the magnitude of offsets in R. A. and decl.

We first find the ephemeris for the star over the epochs of *Hipparcos* measurements ($t$) from the best-fit astrometric parameters:

$$ \Delta\alpha^*(t) = \pi(X(t)\sin(\alpha_0) - Y(t)\cos(\alpha_0)) \\ + (t - 1991.25)\mu_{\alpha}^* \quad (1) $$

and

$$ \Delta\delta(t) = \pi(X(t)\cos(\alpha_0)\sin(\delta_0) \\ + Y(t)\sin(\alpha_0)\sin(\delta_0) - Z(t)\cos(\delta_0)) \\ + (t - 1991.25)\mu_{\delta}. \quad (2) $$

Here, $\Delta\alpha^*(t)$ and $\Delta\delta(t)$ represent the offset from the catalog position ($\alpha_0$, $\delta_0$) at the solar system barycenter of the photocenter from proper motion and parallax only. The values of $X$, $Y$, and $Z$ in au represent the location of the Earth in barycentric coordinates. With this ephemeris, we can then reconstruct the abscissa measurement for each *Hipparcos* epoch. The residual gives the difference between this ephemeris and the *Hipparcos* measurement at a time $t$, along the scan direction $\phi$. The abscissa measurement, then, is a line that passes through the point:

$$ \alpha_a^*(t) = R(t)\cos(\phi(t)) + \Delta\alpha^*(t) \\ \delta_a(t) = R(t)\sin(\phi(t)) + \Delta\delta(t), \quad (3) $$

where for convenience, $\alpha_a^*$ and $\delta_a$ are offsets from (0,0), taken to be the *Hipparcos* catalog values of $\alpha_0$ and $\delta_0$. The *Hipparcos* measurement is one-dimensional, and so consists of a line through the point ($\alpha_a^*(t)$, $\delta_a(t)$), but perpendicular to the scan direction. We define such a line by two points each separated by 1 mas from ($\alpha_a^*(t)$, $\delta_a(t)$),

$$ \alpha_M^*(t) = [-1, 1] \times \sin(\phi(t)) + \alpha_a^*(t) \\ \delta_M(t) = [1, -1] \times \cos(\phi(t)) + \delta_a(t). \quad (4) $$

Thus, the *Hipparcos* measurement at epoch $t$ is given by a line passing through the points defined by $\alpha_M^*(t)$ and $\delta_M(t)$. The error from van Leeuwen (2007b) ($\epsilon$) is the distance in mas from this line in the perpendicular direction (along the scan direction). These measurements and errors can then be fit with any astrometric model, either the five-parameter fit performed by van Leeuwen (2007a) or a more complicated combination of these parameters and orbital parameters. For arbitrary functions that give calculated values of position as a function of time $\alpha_C^*(t)$ and $\delta_C(t)$, $\chi^2$ can be calculated by first finding the residual separation ($d$) from the measurement in the perpendicular direction (along the *Hipparcos* scan direction), using the equation for the distance from a point ($x_0$, $y_0$) to a line defined

by the points ($x_1$, $y_1$) and ($x_2$, $y_2$):

$$ d = \frac{|(y_2 - y_1)x_0 - (x_2 - x_1)y_0 + x_2y_1 - y_2x_1|}{\sqrt{(y_2 - y_1)^2 + (x_2 - x_1)^2}}. \quad (5) $$

When we substitute ($x_1$, $x_2$) = $\alpha_M^*(t)$, ($y_1$, $y_2$) = $\delta_M(t)$, $x_0 = \alpha_C^*(t)$, and $y_0 = \delta_C(t)$, the expression for $d$ simplifies to:

$$ d(t) = |(\alpha_a^*(t) - \alpha_C^*(t))\cos(\phi(t)) \\ + (\delta_a(t) - \delta_C(t))\sin(\phi(t))|, \quad (6) $$

which allows us to calculate the $\chi^2$ of a given model from

$$ \chi^2 = \sum_t \left( \frac{d(t)}{\epsilon(t)} \right)^2. \quad (7) $$

To test the consistency of this method, we extract the abscissa measurements of $\beta$ Pic from van Leeuwen (2007a) and van Leeuwen (2007b), which consist of 111 epochs between 1990.005 and 1993.096, and we then refit them with the same five-parameter model. We use a Metropolis-Hastings MCMC procedure (Nielsen et al. 2014) to sample the posterior of the five parameters $\alpha_{H0}^*$, $\delta_{H0}$, $\pi$, $\mu_{\alpha^*}$, $\mu_\delta$, and compare to the values and errors given by van Leeuwen (2007a). We define $\alpha_{H0}^*$ and $\delta_{H0}$ as the offsets in mas of the photocenter in 1991.25 from the *Hipparcos* catalog positions $\alpha_0$ and $\delta_0$ as measured from the solar system barycenter. Thus, in this five-parameter fit, our model has values for $\alpha_C^*(t)$ and $\delta_C(t)$ of

$$ \alpha_C^*(t) = \alpha_{H0}^* + \pi(X(t)\sin(\alpha_0) - Y(t)\cos(\alpha_0)) \\ + (t - 1991.25)\mu_{\alpha}^* \quad (8) $$

and

$$ \delta_C(t) = \delta_{H0} + \pi(X(t)\cos(\alpha_0)\sin(\delta_0) \\ + Y(t)\sin(\alpha_0)\sin(\delta_0) - Z(t)\cos(\delta_0)) \\ + (t - 1991.25)\mu_{\delta}. \quad (9) $$

A simple fit of the extracted abscissa values and errors produces posteriors with median values that match the catalog values, but with standard deviations that are ~10% too large. This discrepancy arises because the catalog errors are renormalized to achieve $\chi^2 = 1$; to reproduce this renormalization, we multiply the individual errors on each abscissa measurement ($\epsilon(t)$) by a factor $f$:

$$ f = D \left( G \sqrt{\frac{2}{9D}} + 1 - \left( \frac{2}{9D} \right)^3 \right), \quad (10) $$

where $D$ is the number of degrees of freedom ($N_{\text{data}} - N_{\text{parameters}} - 1 = N_{\text{epochs}} - 6$) and $G$ is the goodness of fit (Michalik et al. 2014). The value for $G$ for $\beta$ Pic is -1.63, as given by van Leeuwen (2007a). Figure 2 shows the comparison after performing this renormalization of the errors, with our fit in the filled red histogram. The van Leeuwen (2007a) *Hipparcos* catalog values are represented as the black curve, taken to be a Gaussian with mean equal to the catalog measurement, and standard deviation the catalog error. The two match to within the numerical precision of the catalog values. We conclude that the abscissa measurements we extract from the *Hipparcos* IAD are suitable for including in our orbit fits of the system.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Refit of extracted *Hipparcos* IAD abscissa measurements, with the 1D posterior on each parameter from our MCMC fit to the IAD shown in the red filled histogram, and a Gaussian probability distribution using the *Hipparcos* catalog values and errors for each parameter shown as an overplotted black curve. We find excellent agreement between our MCMC fit and the van Leeuwen (2007a) *Hipparcos* catalog values.</figcaption></figure>

## 3.2. Gaia DR2

The Gaia DR2 magnitude of β Pic is $G = 3.72$, and it is therefore a star that lies outside the nominal magnitude range of the Gaia mission (Gaia Collaboration et al. 2016). It is being observed because small improvements to the onboard detection parameters were made before routine operations began (Martín-Fleitas et al. 2014; Sahlmann et al. 2016b). However, it can be expected that the degraded astrometric performance for bright stars in the range $G = 5–6$ observed in DR2 (e.g., Lindegren et al. 2018, Figure 9) is even more pronounced for brighter stars like β Pic. The data for β Pic in Gaia DR2 therefore must be treated with additional caution.

To establish a notion of the quality of the DR2 data, we compared several quality indicators for a comparison sample of stars, chosen to have magnitudes within ±1 of β Pic. We used pygacs<fnref n="37" /> to query the Gaia archive, and retrieved 1494 very bright stars with $G = 2.72–4.72$. Figure 3 shows a small selection of DR2 catalog parameters, and we inspected many more. From this comparison, β Pic appears to be a “typical” very bright star in terms of excess noise, parameter uncertainties, and number of Gaia observations, with no indication of being particularly problematic.

In particular, the `astrometric_excess.noise` of 2.14 mas is large when compared to stars in the nominal Gaia magnitude range, but not outstanding when compared to other very bright stars. If the excess noise were normally distributed, we would expect it to average out with $1/\sqrt{\text{astrometric_matchdObservations}} = 30$,

yielding 0.39 mas, which is comparable to the DR2 errors in positions and parallax (0.32–0.34 mas).

As in Snellen & Brown (2018), we also make use of the Gaia DR2 data (Gaia Collaboration et al. 2018) to further constrain the orbit. The IAD from Gaia are not yet publicly available, and so we can only utilize the catalog values from the five-parameter fit. As Snellen & Brown (2018) note, $\alpha_G$ and $\delta_G$, the solar system barycentric coordinates of the star at a Gaia reference epoch of 2015.5, strongly constrain the proper motion, given the long time baseline to the 1989–1993 Hipparcos data. As both measurements are in the solar system barycentric frame (ICRS J2000), the offset between $(\alpha_G, \delta_G)$ and the Hipparcos values $(\alpha_0, \delta_0)$ at the reference epoch of 1991.25 should be a combination of proper motion of the system and orbital motion.

## 3.3. Orbit Fitting Results

### 3.3.1. Relative Astrometry Only

Before including the *Hipparcos* and *Gaia* data, we begin by fitting an orbit to the direct imaging data alone. We again utilize the MCMC Metropolis Hastings orbit fitting procedure described previously in Nielsen et al. (2014, 2016) and Nielsen et al. (2017). We perform a fit in seven parameters, with the typical priors for visual orbits, semimajor axis ($a$) uniform in $\log(a)$ ($\frac{dN}{d\log a} \propto C$, which is equivalent to $\frac{dN}{da} \propto a^{-1}$), uniform eccentricity ($e$), inclination angle ($i$) uniform in $\cos(i)$, and uniform in argument of periastron ($\omega$), position angle of nodes ($\Omega$), epoch of periastron passage ($T_0$), and total mass ($M_{\text{tot}}$). Period ($P$) is then derived from $a$ and $M_{\text{tot}}$ using Kepler's third law. The distance for this fit is set to be fixed at the *Hipparcos*

<footnote n="37">³⁷ https://github.com/Johannes-Sahlmann/pygacs</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Gaia DR2 parameters of β Pic (large green circle), compared with ~1500 stars with similar magnitudes. The x-axis is the star sequence number.</figcaption></figure>

value of 19.44 pc (van Leeuwen 2007a). To avoid systematic offsets between different instruments as much as possible, we do not fit the SPHERE data from Lagrange et al. (2019a) and we limit our fit to the data set of Wang et al. (2016). Fits to imaging data sets have a well-known degeneracy in orbital parameters between [$\omega, \Omega$] and [$\omega+180^\circ, \Omega+180^\circ$], a degeneracy that is classically broken with RV observations. In the case of β Pic b, an RV measurement has been made for the planet itself by Snellen et al. (2014), who find the RV of the planet, with respect to the host star, to be $-15.4 \pm 1.7 \text{ km s}^{-1}$, at 2013 December 17. We include that RV data point in this and subsequent fits. We refer to this data set as “Case 1.” Given the changes to the GPI astrometry, we find an orbit fit that is shifted toward lower periods and more circular orbits. Wang et al. (2016) reported [$a, e, i, \omega, \Omega, \tau, P, M_{\text{tot}}$] of [$9.66_{-0.64}^{+1.12}$ au, $0.080_{-0.053}^{+0.091}$, $88_{-0.11}^{+0.12}$, $205_{-13.0}^{+52.6}$, $31_{-0.09}^{+0.80}$, $0.73_{-0.41}^{+0.14}$, $22.47_{-2.26}^{+3.77}$ yr, $1.80_{-0.04}^{+0.03} M_{\odot}$], compared to our values for “Case 1” of [$8.95_{-0.32}^{+0.30}$ au, $0.0360_{-0.022}^{+0.029}$, $88_{-0.97}^{+80}$, $290_{-73.8}^{+80}$, $32_{-0.97}^{+0.09}$, $1.14_{-0.26}^{+0.22}$, $20.18_{-0.97}^{+1.05}$ yr, $1.75 \pm 0.03 M_{\odot}$]. Wang

et al. (2016) define epoch of periastron passage, $\tau$, as the number of orbital periods from MJD = 50,000 (1995.7726), and we converted our value of $T_0$ to this convention for this comparison.

Next, we repeat the orbit fit, but including the two additional epochs of GPI data from 2018 described in Section 2.1, which we refer to as “Case 2.” We display the posteriors for this fit in Figure 4. The orbits themselves are shown in Figure 5, with posteriors for this and all other orbits given in Table 3.

Generally, low-eccentricity orbits are preferred, with a peak at $e = 0$, and with a strong correlation between eccentricity and semimajor axis. Periastron is preferred to be near 2014 (2013.5$_{-0.7}^{+3.4}$), with nonzero probability across 20 yr, corresponding to circular orbits where periastron is undefined.

Figure 6 compares posteriors on five parameters for the Case 1 and Case 2 fits. Including GPI data after conjunction results in higher probability of more eccentric orbits, larger periods, and larger total mass for the system.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Triangle plot for the orbit fit to β Pic b using only the imaging data from NaCo, NICI, Magellan, and GPI (Case 2). A strong degeneracy exists between eccentricity and semimajor axis, with more eccentric orbits having longer periods.</figcaption></figure>

### 3.3.2. Relative and Absolute Astrometry

We next include the *Hipparcos* and *Gaia* data in our fit. Here, we add six more parameters to the previous seven (*a*, *e*, *i*, *ω*, *Ω*, *T*₀, *M*<sub>tot</sub>), for a total of thirteen. The additional parameters are: mass of the planet, in *M*<sub>Jup</sub> (*M*<sub>P</sub>); the location of the star from the solar system barycenter at the *Hipparcos* reference epoch of 1991.25, (α<sub>H0</sub><sup>*</sup>, δ<sub>H0</sub><sup>*</sup>), both values being expressed as an offset from the van Leeuwen (2007a) catalog positions, in mas; parallax (π), in mas; and proper motion (μ<sub>α</sub><sup>*</sup>, μ<sub>δ</sub><sup>*</sup>), in mas yr<sup>-1</sup>. As before, α<sub>H0</sub><sup>*</sup> and μ<sub>α</sub><sup>*</sup> indicate α<sub>H0</sub> cos(δ<sub>0</sub>) and μ<sub>α</sub> cos(δ<sub>0</sub>), in order to correct for the nonrectilinear nature of the

coordinate system. Uniform priors are assumed for all six additional parameters.

Our data set includes the imaging data and planet RV used in the previous fit, as well as our extracted abscissa measurements and errors from the *Hipparcos* IAD, and the *Gaia* DR2 values of α<sub>G</sub> and δ<sub>G</sub> and associated errors. Thus, χ<sup>2</sup> has four components. The first component consists of the standard separation and position angles for the imaging data and errors, with calculated values taken from the seven imaging data orbital parameters, and the distance taken from the parallax parameter. The second component is the CRIRES RV of the planet from Snellen et al. (2014), with reported errors.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Orbit tracks for the orbit fit using only the imaging data (Case 2). The black line shows the lowest $\chi^2$ orbit, while the blue curves are 100 sets of orbital parameters drawn from the posterior.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Posterior probability distributions for semimajor axis (au), eccentricity, inclination angle (degrees), total mass ($M_{\odot}$), and period (yr), for the imaging-only fit with the new 2018 GPI data (Case 2, green) and without (Case 1, pink). The new data, following conjunction, result in more probability at orbits with larger eccentricity and semimajor axis.</figcaption></figure>

The third component, the IAD contribution, comes from Equation (7) and all thirteen parameters, with the position as a function of *Hipparcos* epoch calculated from the standard five astrometric parameters, and additional displacement given by the motion of β Pic around the center of mass of the star/planet system. We approximate β Pic b as having zero flux in the *Hipparcos* and *Gaia* bandpasses. From the BT-Settl models (Baraffe et al. 2015), at 26 Myr and 20 $M_{Jup}$ (well above the expected mass of ~12 $M_{Jup}$), β Pic b would have an apparent magnitude in the *Gaia* G bandpass of 16.9 mags, which is 13.1 mags fainter than β Pic. From our MCMC fit to the visual data alone, the maximum value of apastron reached was 0.8; even at this value, the offset between the photocenter and the star itself in the *Gaia* data is 0.005 mas, well below the precision of any of the measurements. The parameters for the visual orbit give the motion of the planet around the star ($\Delta\alpha_{\text{V}}^*$, $\Delta\delta_{\text{V}}$), and so the motion of the star around the barycenter is then $\Delta\alpha_{\text{s}}^* = -\Delta\alpha_{\text{V}}^* M_{\text{P}} / M_{\text{tot}}$, and similar for $\Delta\delta_{\text{s}}$. The values of $\Delta\alpha_{\text{s}}^*$ and

$\Delta\delta_{\text{s}}$ are calculated at 1991.25 and subtracted from each *Hipparcos* epoch to give the relative motion since the reference values of $\alpha_0$ and $\delta_0$.

The final components of the $\chi^2$ come from fits to the *Gaia* values of $\alpha_G$ and $\delta_G$. We fit the offset between these two values and our fit parameters, $(\alpha_G - \alpha_0 - \alpha_{H0}^*) \times \cos\delta_0$ and $\delta_G - \delta_0 - \delta_{H0}$, with errors given by the stated errors in the *Gaia* DR2 catalog. We then fit this offset from the combination of the astrometric motion from $\mu_{\alpha^*}$ and $\mu_{\delta}$, as well as the orbital motion using the same method as for the *Hipparcos* IAD. We do not incorporate corrections to the nonrectilinear coordinate system or relativistic effects described by Butkevich & Lindegren (2014), given that the *Gaia* error bars are significantly larger than the magnitude of these effects.

We refer to this orbit fit, to the *Hipparcos* and *Gaia* absolute astrometry, the CRIRES RV, and the relative astrometry from NACO, NICI, Magellan, and GPI, as “Case 3.” These results are presented in Figures 7 and 8, Table 2 and are the boldface
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Triangle plot for the orbit fit to the imaging data set of NaCo, NICI, MagAO, and GPI, the CRIRES RV, as well as the astrometric data from *Hipparcos* and *Gaia* (Case 3). With the addition of the astrometry, slightly larger eccentricities are preferred, and thus slightly larger orbital periods.</figcaption></figure>

adopted values in Table 3. In this combined fit, the eccentricity has shifted upward slightly, with eccentricity $\le 0.05$ no longer allowed. The other imaging parameters are similar to our previous imaging-only fit. Astrometric parameters are similar to the van Leeuwen (2007a) *Hipparcos* catalog values as well. Offsets from the van Leeuwen (2007a) reference location ($\alpha_{H0}^*$ and $\delta_{H0}$) are $0.06 \pm 0.11$ mas and $0.03 \pm 0.13$ mas, respectively. Parallax of $51.44 \pm 0.13$ mas is essentially the same as the *Hipparcos* catalog value of $51.44 \pm 0.12$ mas. Meanwhile, as expected for significant reflex motion, we infer the proper motion of the system ($\mu_{\alpha^*}$, $\mu_{\delta}$) to be ($+4.94 \pm 0.02, +83.93_{-0.04}^{+0.03}$) mas yr$^{-1}$, different from their

catalog values of ($+4.65 \pm 0.11, +83.10 \pm 0.15$) mas yr$^{-1}$, by $2.2\sigma$ and $4.4\sigma$, respectively.

In Figure 9, we plot the predicted proper motion from our Case 3 orbit fit of $\beta$ Pic as a function of time. The proper motion is well-constrained by the *Hipparcos* IAD measurement between 1990 and 1993, and matches our accuracy on the system proper motion for this orbit fit ($+4.94 \pm 0.02, +83.93_{-0.04}^{+0.03}$) mas yr$^{-1}$. Though we do not include the *Gaia* DR2 proper motion measurement in this fit, we mark its location as points with error bars at 2015.5. We note that the *Gaia* DR2 proper motion errors ($\pm 0.68$ mas yr$^{-1}$) are considerably larger than the *Hipparcos* values of van
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Orbit tracks for the orbit fit to both the imaging and astrometric data sets (Case 3). While the uncertainty in eccentricity remains, with zero eccentricity orbits no longer allowed, longer orbital periods are preferred.</figcaption></figure>

<table><caption>**Table 2**
Properties of the β Pic System</caption><thead><tr><th></th><th>β Pic</th><th>β Pic b</th><th>References</th></tr></thead><tbody><tr><td>α (deg)</td><td>86.82123366090</td><td></td><td>Gaia Collaboration et al. (2018)</td></tr><tr><td>δ (deg)</td><td>-51.06614803159</td><td></td><td>Gaia Collaboration et al. (2018)</td></tr><tr><td>μ<sub>α*</sub> (mas yr<sup>-1</sup>)</td><td>4.94±0.02</td><td></td><td>This Work</td></tr><tr><td>μ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td>83.93<sup>+0.03</sup><sub>-0.04</sub></td><td></td><td>This Work</td></tr><tr><td>π (mas)</td><td>51.44 ± 0.13</td><td></td><td>This Work</td></tr><tr><td>d (pc)</td><td>19.44 ± 0.05</td><td></td><td>This Work</td></tr><tr><td>M</td><td>1.77 ± 0.03 M<sub>☉</sub></td><td>12.8<sup>+5.5</sup><sub>-3.2</sub> M<sub>☉</sub></td><td>This Work</td></tr><tr><td>log L/L<sub>☉</sub></td><td></td><td>-3.76 ± 0.02</td><td>Chilcote et al. (2017)</td></tr><tr><td>a (au)</td><td></td><td>10.2<sup>+0.4</sup><sub>-0.3</sub></td><td>This Work</td></tr><tr><td>e</td><td></td><td>0.12<sup>+0.04</sup><sub>-0.03</sub></td><td>This Work</td></tr><tr><td>i (deg)</td><td></td><td>88.88 ± 0.09</td><td>This Work</td></tr><tr><td>ω (deg)</td><td></td><td>198±4</td><td>This Work</td></tr><tr><td>Ω (deg)</td><td></td><td>32.05 ± 0.07</td><td>This Work</td></tr><tr><td>T<sub>0</sub></td><td></td><td>2013.7±0.2</td><td>This Work</td></tr><tr><td>P (yr)</td><td></td><td>24.3<sup>+1.5</sup><sub>-1.0</sub></td><td>This Work</td></tr></tbody></table>

Leeuwen (2007a). While the proper motion in decl. is a good match to the tracks, the R.A. proper motion is significantly off from the tracks. It is unclear whether this is a result of systematics in extracting astrometry from bright stars or this offset is the effect of attempting to fit an acceleration in proper motion over a 1.5 yr time baseline with a five-parameter fit. If future Gaia data releases are able to reach <0.1 mas yr⁻¹ proper motion precision, it should greatly reduce the errors in the measurement of the mass of the planet.

### 3.3.3. Independent Analysis

To probe the robustness of our results against different methods and algorithms, we performed a second, independent analysis of the same data set—in this case, the data set discussed above, as well as the SPHERE relative astrometry (referred to as “Case 5” below). We reconstructed the HIP2 (van Leeuwen 2007b) IAD abscissa using the method described in Sahlmann et al. (2011), Section 3.1. When fitting the standard linear five-parameter model, we recovered the HIP2 catalog parameters and obtained a residual rms of 0.79 mas. When adding the Gaia DR2 position of β Pic (Gaia DR2 4792774797545105664) the rms in the HIP2 residuals increases to 0.89 mas. To correctly include the parallax-free Gaia DR2 catalog position in the fit, we set the corresponding parallax factors to zero.

In combination with the ground-based relative astrometry of β Pic b, the Hipparcos and Gaia absolute astrometry allows us to determine model-independent dynamical masses of β Pic and its planetary companion, under the assumption that the space-based astrometry is unbiased (see next section). We performed an MCMC analysis similar to Sahlmann et al. (2016a, 2013). The 13 free parameters are *P*, *e*, *i*, *ω*, *T*<sub>P</sub>, *M*<sub>⋆</sub>, *M*<sub>b</sub>, Ω (eight parameters for the orbital motion) and α<sub>2012</sub><sup>*</sup>, δ<sub>2012</sub><sup>*</sup>, ϖ, μ<sub>α*</sub><sup>*</sup>, μ<sub>δ</sub><sup>*</sup> (five parameters for the standard astrometric model), where we defined ω as the argument of periastron for the barycentric orbit of the primary (in the previous sections, ω referred to the relative orbit). In the MCMC, we adjusted the pair √e cos ω and √e sin ω instead of *e* and ω, to mitigate the effect of correlations that naturally exists between those parameters. We also chose the reference epoch at year 2002 (between the Hipparcos and Gaia epochs), to mitigate correlations between positional offsets and proper motions. Additionally, values of α<sub>2012</sub><sup>*</sup> and δ<sub>2012</sub><sup>*</sup> here correspond to the
<table><caption>**Table 3**
Orbit Fitting Posterkors</caption><thead><tr><th></th><th>a (au)</th><th>e</th><th>i (deg)</th><th>ω (deg)</th><th>Ω (deg)</th><th>T<sub>0</sub></th><th>P (yr)</th><th>M<sub>tot</sub></th><th>M<sub>P</sub> (M<sub>Jup</sub>)</th><th>α<sup>*</sup><sub>H0</sub> (mas)</th><th>δ<sup>*</sup><sub>H0</sub> (mas)</th><th>π (mas)</th><th>μ<sub>α</sub><sup>*</sup></th><th>μ<sub>δ</sub></th><th>ρ<sub>S</sub>/ρ<sub>G</sub></th><th>θ<sub>S</sub> - θ<sub>G</sub> (deg)</th></tr></thead><tbody><tr><td colspan="17">Case 1: NACO, NICI, MagAO, GPI, prior to 2018</td></tr><tr><td>Median</td><td>8.95</td><td>0.036</td><td>88.803</td><td>297.33</td><td>32.025</td><td>2018.55</td><td>20.2</td><td>1.75</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>68% CL min.</td><td>8.63</td><td>0.014</td><td>88.684</td><td>222.43</td><td>31.939</td><td>2014.59</td><td>19.2</td><td>1.72</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>68% CL max.</td><td>9.25</td><td>0.065</td><td>88.921</td><td>355.43</td><td>32.110</td><td>2021.47</td><td>21.2</td><td>1.78</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>95% CL min.</td><td>8.30</td><td>0.002</td><td>88.564</td><td>194.72</td><td>31.853</td><td>2009.06</td><td>18.2</td><td>1.69</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>95% CL max.</td><td>9.87</td><td>0.106</td><td>89.038</td><td>478.25</td><td>32.193</td><td>2023.33</td><td>23.3</td><td>1.81</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td colspan="17">Case 2: NACO, NICI, MagAO, GPI</td></tr><tr><td>Median</td><td>9.32</td><td>0.038</td><td>88.835</td><td>203.63</td><td>32.005</td><td>2013.49</td><td>21.3</td><td>1.78</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>68% CL min.</td><td>9.01</td><td>0.009</td><td>88.734</td><td>192.91</td><td>31.924</td><td>2010.05</td><td>20.3</td><td>1.76</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>68% CL max.</td><td>9.96</td><td>0.101</td><td>88.935</td><td>475.83</td><td>32.084</td><td>2014.23</td><td>23.5</td><td>1.81</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>95% CL min.</td><td>8.81</td><td>0.001</td><td>88.633</td><td>183.19</td><td>31.843</td><td>2005.52</td><td>19.5</td><td>1.74</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>95% CL max.</td><td>10.59</td><td>0.158</td><td>89.034</td><td>535.92</td><td>32.163</td><td>2024.12</td><td>25.8</td><td>1.83</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td colspan="17">Case 3: NACO, NICI, MagAO, GPI, Hipparcos, Gaia (Adopted Parameters)</td></tr><tr><td>Median</td><td><b>10.18</b></td><td><b>0.122</b></td><td><b>88.877</b></td><td><b>18.12</b></td><td>-147.954</td><td><b>2013.72</b></td><td><b>24.3</b></td><td><b>1.79</b></td><td><b>12.82</b></td><td><b>0.057</b></td><td><b>0.034</b></td><td><b>51.439</b></td><td><b>4.936</b></td><td><b>83.927</b></td><td></td><td></td></tr><tr><td>68% CL min.</td><td><b>9.89</b></td><td><b>0.094</b></td><td><b>88.782</b></td><td><b>14.09</b></td><td>-148.029</td><td><b>2013.47</b></td><td><b>23.3</b></td><td><b>1.76</b></td><td><b>9.58</b></td><td>-0.058</td><td>-0.096</td><td><b>51.309</b></td><td><b>4.911</b></td><td><b>83.890</b></td><td></td><td></td></tr><tr><td>68% CL max.</td><td><b>10.60</b></td><td><b>0.159</b></td><td><b>88.971</b></td><td><b>22.04</b></td><td>-147.880</td><td><b>2013.95</b></td><td><b>25.8</b></td><td><b>1.81</b></td><td><b>18.29</b></td><td><b>0.173</b></td><td><b>0.166</b></td><td><b>51.570</b></td><td><b>4.960</b></td><td><b>83.962</b></td><td></td><td></td></tr><tr><td>95% CL min.</td><td><b>9.69</b></td><td><b>0.074</b></td><td><b>88.687</b></td><td><b>9.16</b></td><td>-148.103</td><td><b>2013.17</b></td><td><b>22.6</b></td><td><b>1.73</b></td><td><b>7.07</b></td><td>-0.173</td><td>-0.226</td><td><b>51.177</b></td><td><b>4.880</b></td><td><b>83.841</b></td><td></td><td></td></tr><tr><td>95% CL max.</td><td><b>11.10</b></td><td><b>0.199</b></td><td><b>89.064</b></td><td><b>26.32</b></td><td>-147.806</td><td><b>2014.19</b></td><td><b>27.6</b></td><td><b>1.84</b></td><td><b>27.07</b></td><td><b>0.290</b></td><td><b>0.298</b></td><td><b>51.700</b></td><td><b>4.984</b></td><td><b>84.000</b></td><td></td><td></td></tr><tr><td colspan="17">Case 4: NACO, NICI, MagAO, GPI, Hipparcos, Gaia, prior to 2018</td></tr><tr><td>Median</td><td>10.07</td><td>0.116</td><td>88.871</td><td>22.28</td><td>-147.943</td><td>2013.96</td><td>24.0</td><td>1.77</td><td>13.61</td><td>0.065</td><td>0.045</td><td>51.437</td><td>4.941</td><td>83.936</td><td></td><td></td></tr><tr><td>68% CL min.</td><td>9.78</td><td>0.091</td><td>88.760</td><td>14.72</td><td>-148.026</td><td>2013.52</td><td>23.0</td><td>1.74</td><td>9.91</td><td>-0.051</td><td>-0.087</td><td>51.306</td><td>4.915</td><td>83.896</td><td></td><td></td></tr><tr><td>68% CL max.</td><td>10.52</td><td>0.154</td><td>88.981</td><td>31.35</td><td>-147.862</td><td>2014.48</td><td>25.6</td><td>1.81</td><td>20.02</td><td>0.182</td><td>0.178</td><td>51.567</td><td>4.965</td><td>83.972</td><td></td><td></td></tr><tr><td>95% CL min.</td><td>9.60</td><td>0.074</td><td>88.649</td><td>8.05</td><td>-148.108</td><td>2013.10</td><td>22.5</td><td>1.71</td><td>7.24</td><td>-0.167</td><td>-0.218</td><td>51.176</td><td>4.878</td><td>83.837</td><td></td><td></td></tr><tr><td>95% CL max.</td><td>11.17</td><td>0.202</td><td>89.090</td><td>41.66</td><td>-147.781</td><td>2015.06</td><td>27.8</td><td>1.84</td><td>29.28</td><td>0.300</td><td>0.312</td><td>51.699</td><td>4.989</td><td>84.010</td><td></td><td></td></tr><tr><td colspan="17">Case 5: NACO, NICI, MagAO, GPI, SPHERE, Hipparcos, Gaia</td></tr><tr><td>Median</td><td>10.03</td><td>0.108</td><td>89.003</td><td>17.81</td><td>-147.921</td><td>2013.69</td><td>23.9</td><td>1.77</td><td>14.25</td><td>0.069</td><td>0.050</td><td>51.436</td><td>4.944</td><td>83.940</td><td></td><td></td></tr><tr><td>68% CL min.</td><td>9.79</td><td>0.085</td><td>88.925</td><td>13.23</td><td>-147.991</td><td>2013.41</td><td>23.0</td><td>1.75</td><td>10.32</td><td>-0.047</td><td>-0.081</td><td>51.305</td><td>4.921</td><td>83.907</td><td></td><td></td></tr><tr><td>68% CL max.</td><td>10.36</td><td>0.139</td><td>89.081</td><td>21.94</td><td>-147.850</td><td>2013.94</td><td>25.0</td><td>1.80</td><td>20.60</td><td>0.185</td><td>0.181</td><td>51.567</td><td>4.966</td><td>83.974</td><td></td><td></td></tr><tr><td>95% CL min.</td><td>9.63</td><td>0.069</td><td>88.846</td><td>8.12</td><td>-148.062</td><td>2013.10</td><td>22.5</td><td>1.72</td><td>7.51</td><td>-0.162</td><td>-0.211</td><td>51.174</td><td>4.896</td><td>83.869</td><td></td><td></td></tr><tr><td>95% CL max.</td><td>10.77</td><td>0.174</td><td>89.159</td><td>27.09</td><td>-147.780</td><td>2014.23</td><td>26.5</td><td>1.83</td><td>29.69</td><td>0.302</td><td>0.314</td><td>51.698</td><td>4.990</td><td>84.011</td><td></td><td></td></tr><tr><td colspan="17">Case 6: NACO, NICI, MagAO, GPI, SPHERE, Hipparcos, Gaia, with calibration offsets</td></tr><tr><td>Median</td><td>10.10</td><td>0.114</td><td>88.986</td><td>17.82</td><td>-148.337</td><td>2013.70</td><td>24.0</td><td>1.78</td><td>13.71</td><td>0.063</td><td>0.044</td><td>51.436</td><td>4.941</td><td>83.935</td><td>1.001</td><td>-0.471</td></tr><tr><td>68% CL min.</td><td>9.82</td><td>0.087</td><td>88.906</td><td>13.78</td><td>-148.483</td><td>2013.45</td><td>23.1</td><td>1.75</td><td>10.04</td><td>-0.053</td><td>-0.086</td><td>51.304</td><td>4.917</td><td>83.899</td><td>0.998</td><td>-0.613</td></tr><tr><td>68% CL max.</td><td>10.49</td><td>0.148</td><td>89.064</td><td>21.84</td><td>-148.192</td><td>2013.94</td><td>25.4</td><td>1.82</td><td>19.90</td><td>0.179</td><td>0.176</td><td>51.568</td><td>4.964</td><td>83.970</td><td>1.004</td><td>-0.329</td></tr><tr><td>95% CL min.</td><td>9.62</td><td>0.069</td><td>88.827</td><td>8.84</td><td>-148.629</td><td>2013.15</td><td>22.5</td><td>1.72</td><td>7.34</td><td>-0.168</td><td>-0.217</td><td>51.173</td><td>4.888</td><td>83.852</td><td>0.996</td><td>-0.757</td></tr><tr><td>95% CL max.</td><td>11.01</td><td>0.190</td><td>89.143</td><td>26.77</td><td>-148.046</td><td>2014.21</td><td>27.2</td><td>1.85</td><td>29.06</td><td>0.297</td><td>0.309</td><td>51.699</td><td>4.988</td><td>84.009</td><td>1.007</td><td>-0.186</td></tr><tr><td colspan="17">Case 7: NACO, SPHERE, Hipparcos, Gaia</td></tr><tr><td>Median</td><td>9.82</td><td>0.101</td><td>89.115</td><td>23.34</td><td>-148.111</td><td>2014.03</td><td>23.5</td><td>1.72</td><td>15.77</td><td>0.082</td><td>0.068</td><td>51.431</td><td>4.951</td><td>83.952</td><td></td><td></td></tr><tr><td>68% CL min.</td><td>9.57</td><td>0.081</td><td>88.953</td><td>13.63</td><td>-148.376</td><td>2013.44</td><td>22.8</td><td>1.65</td><td>11.06</td><td>-0.035</td><td>-0.064</td><td>51.300</td><td>4.928</td><td>83.918</td><td></td><td></td></tr><tr><td>68% CL max.</td><td>10.18</td><td>0.130</td><td>89.277</td><td>34.22</td><td>-147.847</td><td>2014.65</td><td>24.6</td><td>1.79</td><td>22.89</td><td>0.199</td><td>0.201</td><td>51.561</td><td>4.973</td><td>83.985</td><td></td><td></td></tr></tbody></table>
<table><caption>**Table 3**
(Continued)</caption><thead><tr><th></th><th>a (au)</th><th>e</th><th>i (deg)</th><th>ω (deg)</th><th>Ω (deg)</th><th>T<sub>0</sub></th><th>P (yr)</th><th>M<sub>tot</sub></th><th>M<sub>P</sub> (M<sub>Jup</sub>)</th><th>α<sup>*</sup><sub>H0</sub> (mas)</th><th>δ<sup>*</sup><sub>H0</sub> (mas)</th><th>π (mas)</th><th>μ<sub>α</sub><sup>*</sup></th><th>μ<sub>δ</sub></th><th>ρ<sub>S</sub>/ρ<sub>G</sub></th><th>θ<sub>S</sub> - θ<sub>G</sub> (deg)</th></tr></thead><tbody><tr><td>95% CL min.</td><td>9.38</td><td>0.066</td><td>88.788</td><td>4.96</td><td>-148.641</td><td>2012.82</td><td>22.3</td><td>1.59</td><td>7.86</td><td>-0.151</td><td>-0.196</td><td>51.169</td><td>4.902</td><td>83.877</td><td></td></tr><tr><td>95% CL max.</td><td>10.71</td><td>0.170</td><td>89.436</td><td>48.57</td><td>-147.583</td><td>2015.31</td><td>26.2</td><td>1.86</td><td>32.21</td><td>0.316</td><td>0.335</td><td>51.692</td><td>4.996</td><td>84.022</td><td></td></tr></tbody></table>

Note. The boldface type refers to the adopted parameters of Case 3.
<figure><img src="image_13.png" /><figcaption>**Figure 9.** The observed proper motion of β Pic, including the system proper motion and the reflex motion due to the orbit of β Pic b, with the tracks (color-coded by planet mass) drawn from the posterior, again using the orbit fit with all imaging data except SPHERE, as well as *Hipparcos* and *Gaia* (Case 3). Dark gray bars mark the timeframe of the *Hipparcos* and *Gaia* observations, with the light gray bar representing the expected remaining extent of the full 7 yr *Gaia* mission. The *Hipparcos* IAD constrains the proper motion well between 1990 and 1993, and a more precise *Gaia* proper motion measurement can greatly reduce the error bars on the planet mass.</figcaption></figure>

location of the β Pic barycenter at the reference epoch, while in the previous fit, $\alpha_{H0}^*$ and $\delta_{H0}$ referred to the location of the system photocenter at the reference epoch. All priors are flat and seed values and their uncertainties for the MCMC chains were set based on either the five-parameter fit above or previous orbital solutions. We used 160 walkers with 44,000 steps each and discarded the first 25% of samples, which yields more than five million samples per parameter.

The MCMC chains exhibit stable convergence, and the posterior distributions show clearly peaked shapes. The residual rms in the absolute astrometry (*Hipparcos* and *Gaia*) with the median orbital model is 0.80 mas, which is significantly smaller than the 0.89 mas obtained with the linear model. This confirms that orbital motion is detected in the absolute astrometry.

In terms of system parameters, the results of the two independent analyses are in excellent agreement, as illustrated by Figure 10, with 1D posteriors overlapping. This gives us confidence in the accuracy of both fitting algorithms.

### 3.3.4. Bias in the Gaia DR2 Catalog Parameters Due to Orbital Motion

The source parameters in *Gaia* DR2 were obtained by fitting either a five-parameter model or a two-parameter model to the astrometric data collected by the satellite (Lindegren et al. 2018). For the five-parameter solution of β Pic, this means that any orbital motion present in the *Gaia* astrometry was not accounted for specifically. Orbital motion may instead manifest itself as an increased excess noise or a bias in the DR2 parameters, which is worse for our purposes.

In an attempt to quantify the bias in the DR2 position caused by orbital motion, we simulated the individual *Gaia* observations. We used the *Gaia* Observation Forecast Tool (gost; https://gaia.esac.esa.int/gost/) to predict the *Gaia* focal plane crossings of β Pic in the time range considered in DR2 after the ecliptic pole scanning (2014 August 22T21:00:00–2016 May 23T11:35:00; Lindegren et al. 2018).<fnref n="38" /> Unfortunately, the earliest date accepted by gost is 2014 September 26T00:00:00, but we corrected for the missing month as described below. In the queried time range, gost predicted 26 *Gaia* focal plane crossings in 16 visibility periods. The *Gaia* DR2 catalog reports 30 *astrometric_match-ed_observations* in 15 visibility periods over the slightly longer time range included in DR2. This validates that the gost predictions are a reasonable approximation of the actual *Gaia* observations. To account for the missing first month in the gost prediction, we duplicated the last two gost predictions and prepended them to the list of predictions with timestamps that correspond to the start of the DR2 time range. Our simulated *Gaia* observation setup thus includes 28 focal plane crossings in 17 visibility periods.

We used a set of the 13 parameters fitted in the previous section to compute noiseless *Gaia* along-scan measurements (equivalent to the *Hipparcos* abscissa) that include the orbital motion, setting the reference epoch to 2015.5. Observation times, parallax factors, and scan angles were specified according to the gost predictions. We also computed the model position of the star at epoch 2015 including barycentric orbital motion and proper motion (zero by definition of the reference epoch), but not parallax (by setting the parallax factor to zero), to replicate the parallax-free DR2 catalog position.

We then fitted the standard five-parameter linear model to the simulated *Gaia* data of β Pic and compared the 2015.5 model position to the best-fit position offsets. The difference between the two corresponds to our estimate of the DR2 position bias. When no significant orbital motion is present (e.g., the planet mass is set to zero), both the model position at 2015.5 and the fitted coordinate offsets of the five-parameter fit are zero, and the input proper motions and parallax are recovered. When orbital motion is present, the actual and the linear-fit position are different.

Because we cannot be certain about the fidelity of the gost predicted DR2 epochs, we estimated the uncertainty in the bias estimation by repeating random draws of 28, 26, and 24 out of 28 predicted epochs. We also incorporated the varying fit parameters by using samples from the MCMC chains in the previous section.

We considered three cases. Case (a) has nominal *Gaia* DR2 positions and uncertainties, with no bias correction. In case (b), when drawing 28 or 26 epochs for the solution in the previous section and 10000 random draws and parameter sets, the DR2 coordinate bias due to orbital motion is estimated to $\epsilon_{R.A.} = 0.017 \pm 0.003$ mas and $\epsilon_{Decl.} = 0.013 \pm 0.003$ mas, which is negligible given the DR2 position uncertainties of ~0.3 mas. In case (c), if we draw 24 epochs, these estimates increase to $\epsilon_{R.A.} = 0.085 \pm 0.142$ mas and $\epsilon_{Decl.} = -0.040 \pm 0.142$ mas, so the bias essentially increases the uncertainty in the DR2 positions and introduces a minor shift.

<footnote n="38">³⁸ Snellen & Brown (2018) mention a time range between 2014 October 1 and 2016 April 19 for β Pic measurements.</footnote>
<figure><img src="image_14.png" /><figcaption>Figure 10. Comparison of the two orbit-fitting techniques for Case 5 shows excellent agreement of the two sets of posteriors.</figcaption></figure>

We repeated the MCMC analysis in all three cases. When debiasing the DR2 position, we subtracted $\epsilon$ from the catalog coordinate before including it in the fit, and we added the bias uncertainty in quadrature to the DR2 position uncertainty. We found that the effect of the position bias as estimated above on the solution parameters is negligible. We illustrate this in Figure 11, where we show posteriors for several fit parameters that are essentially indistinguishable. The same applies to all other parameters.

Although the bias of the DR2 parameters due to orbital motion is negligible for our purposes, this is certainly not the general case. For instance, we found that the bias in β Pic's DR2 proper motion is significant: $\epsilon_{\mu_{\alpha^*}} = 0.37 \pm 0.08$ mas yr$^{-1}$ and $\epsilon_{\mu_{\delta}} = 0.62 \pm 0.13$ mas yr$^{-1}$ (the corresponding parallax bias is smaller than 3 $\mu$as). A bias of $\sim 0.4$ mas yr$^{-1}$ in the R.A. direction is not enough to explain the offset seen in Figure 9, where the Gaia value of $\mu_{\alpha^*}$ is $\sim 2$ mas yr$^{-1}$ from the orbit tracks, so the full cause of this offset is still unclear. Likewise, the $\mu_{\delta}$ bias moves the Gaia proper motion even further from the orbit tracks. Caution is therefore necessary when using the DR2 parameters of systems exhibiting orbital motion, and in

particular when determining orbital parameters from the Gaia DR2 catalog in combination with other surveys (e.g., Brandt et al. 2019; Kervella et al. 2019).

### 3.3.5. The Effect of Using Different Data Sets

We consider different combinations of relative astrometry to investigate how different combinations influence the derived mass. In addition to the fit to Hipparcos, Gaia, CRIRES, NACO, NICI, Magellan, and GPI discussed above (“Case 3”), we also consider the effect of the 2018 GPI data by performing a second fit, but without these two GPI data points in 2018 (“Case 4”). We perform three additional fits as well, all using the Hipparcos, Gaia, and CRIRES data: one including the SPHERE data of Lagrange et al. (2019a) (“Case 5”) as presented; another including this SPHERE data, but fitting for additional offset terms for the GPI separation and position angle (“Case 6”); and a third using relative astrometry only from ESO instruments, NACO, and SPHERE (“Case 7”). We present the full set of posteriors in Table 3 for each of these orbit fits.
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Posteriors on star mass, planet mass, eccentricity, and period (yr) with (dashed line: case (b), dashed-dotted line: case (c)) and without (solid line: case (a)) DR2 position bias correction.</figcaption></figure>

Given the small errors on the imaging data (~1 mas for the early GPI data), the importance of the astrometric calibration becomes key, especially when combining data sets from different instruments. An uncorrected offset in calibration (either in plate scale or true north) results in a large acceleration between data points, which the orbit fitter will attempt to compensate for with a more eccentric orbit, where orbital speed can be varied near the problematic epochs. It has been speculated that a calibration offset is the likely cause for different predictions for the Hill Sphere crossing and closest approach of β Pic b in the previous two years (Wang et al. 2016). We examine the influence of this effect by performing multiple orbit fits combining imaging and the *Hipparcos* and *Gaia* data, with different combinations of instruments.

Figure 12 shows posteriors for planet mass, eccentricity, and period for these multiple orbit fit cases. Fits that do not utilize SPHERE data (Cases 3 and 4) give generally lower masses, with slightly larger eccentricity and period, compared to fits that include SPHERE data (Cases 5, 6, 7). The largest difference is from Case 3 (GPI but not SPHERE) with a planetary mass of $12.8_{-3.2}^{+5.5} M_{\text{Jup}}$ to Case 7 (SPHERE but not GPI), where the mass measurement is $15.8_{-4.7}^{+7.1} M_{\text{Jup}}$, with combinations of the two instruments falling in between.

Following the updates to the north angle in the GPI pipeline, we find evidence for a systematic position angle offset between GPI and SPHERE. The Case 6 fit introduced two additional parameters into the fit, a multiplicative offset to GPI separations, and an additive offset to GPI position angles (corresponding to calibration errors in planet scale and true north, respectively). The fit values for these offsets are $\rho_S/\rho_G = 1.001 \pm 0.003$, and $\theta_S - \theta_G = -0.47 \pm 0.14$, suggesting no offset in plate scale, but a true north offset of about half a degree between the two instruments. Figure 13 compares the relative astrometry from GPI and SPHERE, indeed showing SPHERE position angles systematically ~0.5 smaller than GPI data at the same epoch.

Maire et al. (2019) present new astrometry of the planet 51 Eri b from SPHERE and an independent reduction of GPI data. From their analysis, they find a systematic PA offset of $1.0^\circ \pm 0.2^\circ$, with SPHERE having larger values than GPI. With the revised astrometric calibration, we found an offset of $\Delta\theta = -0.16 \pm 0.26$ from a joint fit to our reduction of our GPI data and the SPHERE astrometry published in Maire et al. (2019; De Rosa et al. 2019), consistent with the offset found for β Pic b in this work. Using the old astrometric calibration and data reduced with the same version of the DRP used by Maire et al. (2019), we calculated an offset of $\Delta\theta = 0.28 \pm 0.26$,
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Planet mass, eccentricity, and period posteriors for different data sets. Data sets with GPI but not SPHERE data tend to favor smaller planet masses, and lager eccentricity and period (Cases 3 and 4). Combinations of GPI and SPHERE data (5 and 6) have more probability at larger masses, and a fit that excluded GPI data (7) moves to the largest planet masses.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Comparison of the GPI and SPHERE astrometry between 2013 and 2020, with the lowest $\chi^2$ orbit subtracted off (Case 3). We find no evidence for a systematic offset in plate scale (1.001 ± 0.003), but significant evidence of a position angle offset of -0.47 ± 0.14.</figcaption></figure>

which is closer to the value in Maire et al. (2019) but still significantly different. This lends further evidence to the conclusion that the culprit is not a single constant offset between the two instruments, but perhaps an algorithmic difference in how astrometry is extracted. Indeed, Maire et al. (2019) note that when they refit GPI data on 51 Eri, they find ~0.35 larger values of PA than those presented by De Rosa et al. (2015) for the same data sets. Further analysis is ongoing to determine the precise cause of these offsets, as well as their impact on derived orbital parameters.

### 3.3.6. Comparison to Previous Orbit Fits

In Figure 14, we compare a modified version of our Case 4 to the results from Snellen & Brown (2018), who examined a similar relative astrometric data set. For consistency in this comparison, here we use the published astrometry from Wang et al. (2016), rather than the updated astrometry presented here. We also do not use the CRIRES RV or the 2009 NaCo M-band point for this fit, to match the Wang et al. (2016) orbit fitting. A key difference between the methods is that Snellen & Brown (2018) did not simultaneously fit the relative astrometry and

Hipparcos and Gaia data as we did, but rather they took the Wang et al. (2016) orbital element posteriors as the constraints from the relative astrometric fit. Additionally, while we use the Hipparcos IAD as constraints on the orbit in the plane of the sky, Snellen & Brown (2018) converted the IAD into one-dimensional measurements along the orbital plane given by Wang et al. (2016). When reporting the mass posterior, Snellen & Brown (2018) restricted the fit to the $1\sigma$ period range of Wang et al. (2016), the most circular orbits. However, as seen in Table 3, the addition of the Hipparcos and Gaia data push the visual orbit toward longer periods and higher eccentricity than the Wang et al. (2016) fit to the relative astrometry alone.

As a result, while Snellen & Brown (2018) find a well-constrained mass for the planet of $11 \pm 2 M_{Jup}$, we find a broader range of $12.7_{-3.1}^{+6.4} M_{Jup}$ when analyzing the same data set. A key factor in the reported smaller uncertainty is that Snellen & Brown (2018) fixed a number of parameters, including the mass of the star ($M_{tot}$) and position angle of nodes ($\Omega$), as well as the orbital period of the planet. In the covariance panel between mass and period within the triangle plot of Figure 14, we show the Snellen & Brown (2018) 1, 2, and $3\sigma$ contours, extracted from their Figure 3, against ours. By restricting the period range to the $1\sigma$ range of Wang et al. (2016) of <28 yr, the planet mass appears more constrained than it actually is given the full data set. Including the GPI 2018 astrometry, as well as updating the astrometry following fixes to the pipeline, (Case 3) produces a somewhat more constrained planet mass compared to our modified Case 4, $12.8_{-3.2}^{+5.5} M_{Jup}$, but with error bars still a factor of two larger than reported by Snellen & Brown (2018). This offset illustrates the importance of a simultaneous fit of relative and absolute astrometry, given the complicated covariant structure of such orbits.

Recently, Dupuy et al. (2019) presented a fit to the $\beta$ Pic b orbit based on relative astrometry from the literature, including the Lagrange et al. (2019a) SPHERE measurement from 2018 as well as the *Hipparcos–Gaia* Catalog of Accelerations (Brandt et al. 2019). Their analysis differs from ours in a number of ways; most significantly, they utilize the *Hipparcos* catalog values rather than the *Hipparcos* IAD and their fit includes the *Gaia* proper motion for $\beta$ Pic, though with inflated errors. Additionally, our analysis benefits from the more precise relative astrometry from GPI in 2018. Dupuy et al. (2019) also fit the RV of the star (Lagrange et al. 2012) and of the planet (Snellen et al. 2014), though given the large jitter in the stellar RVs and the moderate error bars on the planet RV, we do not expect the inclusion of RVs to have a significant difference in the two fits. We also find a more constrained parallax for the system (51.44 ± 0.13 mas from our Case 3 fit, largely based on the *Hipparcos* IAD), compared to their inflated *Hipparcos* parallax error, a linear combination of the original ESA (1997) catalog and the rereduced van Leeuwen (2007a) catalog.

Similar to the comparison to Snellen & Brown (2018), we produce a modified Case 3 fit, before the correction to the GPI astrometry, to compare the two methods. Dupuy et al. (2019) find orbital parameters generally similar to our modified Case 3 orbit fit. They find a planet mass of $13.1_{-3.2}^{+2.8} M_{Jup}$, period of $29.9_{-3.2}^{+2.9}$ yr, and eccentricity of 0.24 ± 0.06, compared to our values of $11.1_{-2.3}^{+2.7} M_{Jup}$, 27.1 ± 2.0 yr, and 0.19 ± 0.05. Thus, we find a somewhat lower planet mass, shorter period, and smaller eccentricity, with slightly smaller error bars. We find a
<figure><img src="image_18.png" /><figcaption>Figure 14. Posteriors from the combined imaging and astrometric fit for period, planet mass, and eccentricity, but without the 2018 GPI data (Case 4). Overplotted in the mass/period covariance plot are 1, 2, and 3σ contours extracted from Figure 3 of Snellen & Brown (2018) for the same data set. While we find generally good agreement with covariance contours for periods less than 28 yr, there is significant probability at larger periods and masses, creating a mass measurement ($12.7^{+0.4}_{-3.1} M_{Jup}$) more uncertain than that reported by Snellen & Brown (2018) ($11 \pm 2 M_{Jup}$).</figcaption></figure>

stellar mass of $1.81 \pm 0.03 M_\odot$, similar to the Dupuy et al. (2019) value of $1.84 \pm 0.05 M_\odot$; these two estimates are the first time planet mass and stellar mass have been measured simultaneously from the same fit for a directly imaged planet. Comparing our triangle plot in Figure 7 to their Figure 3, Dupuy et al. (2019) do not reproduce our U-shaped covariance between semimajor axis and planet mass, nor between eccentricity and planet mass; rather, they see a roughly linear relationship for both covariances. The intersection of the two sets of covariances includes short-period lower-mass planets and long-period higher-mass planets, while our results include another family of short-period higher-mass planets not seen by Dupuy et al. (2019). The source of this discrepancy is not clear: given the large error bars on both the RVs and the recomputed Gaia proper motion of Dupuy et al. (2019), neither should significantly move the fit. It is not likely that the GPI 2018 data is to blame, since our Case 7 fit—which is an imaging data set similar to the one used by Dupuy et al. (2019)—also has these U-shaped covariances. For our final Case 3 fit, using updated GPI astrometry, we find a mass of $12.8^{+5.5}_{-3.2} M_{Jup}$, eccentricity of $0.12^{+0.04}_{-0.03}$, and period of $24.3^{+1.5}_{-1.0}$ yr: a larger uncertainty on planet mass, and smaller values of period and eccentricity.

### 3.4. Effects of a Second Giant Planet in the β Pic System

After this paper was submitted, Lagrange et al. (2019b) presented RV measurements of β Pic and an orbit fit for an inner giant planet, β Pic c, orbiting at 2.7 au. In particular, by fitting for the δ Scuti pulsations of the star, they were able to detect the ~4 yr signal of the inner planet. The orbit fit performed took 219 sets of orbital elements from a chain of a separate MCMC orbit fit to the astrometry, and used these

elements as the basis for fitting the RVs of the star, with the added assumption that the two planets are coplanar.

Here, we perform a joint fit to four types of data simultaneously: imaging data from NaCo, NICI, MagAO, and GPI, absolute astrometry from Hipparcos and Gaia, RV of the planet from CRIRES (the three of which constitute Case 3 above), as well as the δ Scuti-corrected RVs and errors of the star taken from supplementary Table 1 of Lagrange et al. (2019b). We expand our thirteen-element orbit fit with an additional eight parameters: semimajor axis, eccentricity, inclination angle, argument of periastron, position angle of nodes, epoch of periastron passage, and planet mass of β Pic c, along with the RV offset of the star. As with β Pic b, we place priors on β Pic c orbital parameters that are uniform in log(a), in eccentricity (e), in cosine of the inclination angle (cos i), in argument of periastron (ω), in position angle of nodes (Ω), in epoch of periastron passage (T₀), in planet mass ($M_c$), and in RV offset (γ). In addition to these priors, we perform a second “coplanar” fit, where the mutual inclination ($i_m$) between the two planets given by:

$$ \cos i_m = \cos i_b \cos i_c + \sin i_b \sin i_c \cos(\Omega_b - \Omega_c), \quad (11) $$

which is constrained to be a Gaussian centered on zero, with standard deviation of 1°. In both cases, the two planets are assumed to not interact with each other, so that the position and RV of the star are just the linear combination of the reflex motion from each planet's orbit. We also note that the RV offset γ does not represent the system velocity, since Lagrange et al. (2019b) have subtracted off the δ Scuti pulsations, and thus any RV offset. Here, γ represents an additional RV correction beyond this.

We give the posteriors to the unconstrained mutual inclination fit in Figure 15 and in Table 4. Despite having no constraint on mutual inclination, the inclination angle and position angle of nodes for c ($i_c$ and $\Omega_c$) differ from the priors, and follow the orbit of β Pic b, but with larger uncertainties: $i_b = 88^{\circ}8 \pm 1^{\circ}0$ and $\Omega_b = 32^{\circ}02 \pm 0^{\circ}08$ for the outer planet, compared to $i_c = 98^{\circ+12}_{-14}$ and $\Omega_c = 36^{\circ} \pm 15^{\circ}$ for the inner planet. Since radial velocities do not constrain either of these parameters, the absolute astrometry of Hipparcos and Gaia must be supplying these constraints.

Other than inclination angle and position angle of nodes for the inner planet, there are no significant differences in the derived posteriors for the parameters of β Pic c between the two fits. The mass of β Pic c changes slightly: in the fit that does not constrain mutual inclination, it is $M_c = 9.4^{+1.1}_{-0.9} M_{Jup}$, compared to $M_c = 9.2^{+1.0}_{-0.9} M_{Jup}$ in the coplanar fit (Figure 16). This is true also for the parameters of the outer planet, β Pic b, as shown in Figure 17, where the two fits incorporating two planets have similar posteriors on β Pic b.

Similarly to Lagrange et al. (2019b), we find the presence of the c planet results in a lower mass for the b planet. In our one-planet Case 3 fit, we found a mass for β Pic b of $M_b = 12.8^{+5.5}_{-3.2} M_{Jup}$, which drops to $M_b = 8.0 \pm 2.6 M_{Jup}$ in the coplanar fit (Figure 17). In this coplanar fit, the semimajor axis, period, and eccentricity posteriors also shift to lower values compared to the one-planet Case 3 fit. A possible explanation for this could be that the evidence for a nonzero eccentricity of β Pic b came from the absolute astrometry of the star, and that this astrometric motion can be equally well explained with a more circular outer planet and a second inner planet.
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Posteriori from the fit to Case 3 and the stellar RVs from Lagrange et al. (2019b), including a second planet (β Pic c), with no additional constraints on mutual inclination between the two planets.</figcaption></figure>

A comparison of our fit values to those of Lagrange et al. (2019b) reveals that their results are generally similar to ours, but with noticeable differences, likely resulting from performing a joint fit for all data and both planets, rather than using MCMC chains from a fit to β Pic b to fit the RVs. Lagrange et al. (2019b) found values of [$a_c, e_c, P_c, M_c$] of [2.69 ± 0.003 au, 0.24 ± 0.02, 3.335 ± 0.005 yr, 8.93 ± 0.14 $M_{Jup}$], compared to values from our coplanar fit of [2.72 ± 0.019 au, 0.24$^{+0.1}_{-0.09}$, 3.39 ± 0.02 yr, 9.18$^{+1.0}_{-0.9} M_{Jup}$]. That these measurements are in such good agreement, but with errors several times larger from the joint fit, suggests that extracting a limited number of orbits from the posterior, as

done in Lagrange et al. (2019b), underestimates the errors on the derived parameters.

In the one-planet fit, we found the *Gaia* proper motion to be significantly offset, i.e., ~2σ from the predicted astrometric motion of the star in R.A. (Figure 9). Considering the proper motion of the star in the coplanar two-planet fit does not resolve this, as Figure 18 shows this offset remains the same. The two planets are similar in mass, but the inner planet accounts for more of the proper motion signature, since stellar orbital velocity scales as $a^{-0.5}$ and β Pic c is ~3.6 times closer to the parent star than β Pic b. The shorter period suggests that future *Gaia* data releases could detect the astrometric motion of
<table><caption>**Table 4**
Two-planet Fit Posteriors</caption><thead><tr><th rowspan="3"></th><th colspan="4">Unconstrained Mutual Inclination</th><th colspan="4">Coplanar</th></tr><tr><th rowspan="2">Median</th><th colspan="2">68% CL</th><th colspan="2">95% CL</th><th rowspan="2">Median</th><th colspan="2">68% CL</th><th colspan="2">95% CL</th></tr><tr><th>Min.</th><th>Max.</th><th>Min.</th><th>Max.</th><th>Min.</th><th>Max.</th><th>Min.</th><th>Max.</th></tr></thead><tbody><tr><td>a<sub>b</sub> (au)</td><td>9.68</td><td>9.43</td><td>9.98</td><td>9.16</td><td>10.40</td><td>9.65</td><td>9.39</td><td>9.95</td><td>9.10</td><td>10.34</td></tr><tr><td>e<sub>b</sub></td><td>0.076</td><td>0.049</td><td>0.105</td><td>0.020</td><td>0.143</td><td>0.072</td><td>0.043</td><td>0.101</td><td>0.012</td><td>0.138</td></tr><tr><td>i<sub>b</sub> (deg)</td><td>88.824</td><td>88.726</td><td>88.922</td><td>88.627</td><td>89.019</td><td>88.826</td><td>88.729</td><td>88.923</td><td>88.630</td><td>89.019</td></tr><tr><td>ω<sub>b</sub> (deg)</td><td>-159.79</td><td>-166.32</td><td>-152.88</td><td>-174.29</td><td>-139.61</td><td>-160.13</td><td>-167.09</td><td>-152.40</td><td>-175.46</td><td>166.13</td></tr><tr><td>Ω<sub>b</sub> (deg)</td><td>32.011</td><td>31.934</td><td>32.087</td><td>31.858</td><td>32.162</td><td>32.008</td><td>31.932</td><td>32.084</td><td>31.855</td><td>32.159</td></tr><tr><td>T<sub>0b</sub></td><td>2013.81</td><td>2013.41</td><td>2014.20</td><td>2012.89</td><td>2014.81</td><td>2013.77</td><td>2013.32</td><td>2014.17</td><td>2012.47</td><td>2015.06</td></tr><tr><td>P<sub>b</sub> (yr)</td><td>22.7</td><td>21.8</td><td>23.7</td><td>20.9</td><td>25.2</td><td>22.5</td><td>21.6</td><td>23.6</td><td>20.7</td><td>25.0</td></tr><tr><td>M<sub>b</sub> (M<sub>Jup</sub>)</td><td>8.35</td><td>5.76</td><td>10.91</td><td>3.19</td><td>13.42</td><td>8.03</td><td>5.41</td><td>10.64</td><td>2.80</td><td>13.20</td></tr><tr><td>a<sub>c</sub> (au)</td><td>2.72</td><td>2.70</td><td>2.74</td><td>2.68</td><td>2.76</td><td>2.72</td><td>2.70</td><td>2.74</td><td>2.68</td><td>2.76</td></tr><tr><td>e<sub>c</sub></td><td>0.248</td><td>0.156</td><td>0.359</td><td>0.063</td><td>0.489</td><td>0.241</td><td>0.151</td><td>0.348</td><td>0.062</td><td>0.476</td></tr><tr><td>i<sub>c</sub> (deg)</td><td>97.860</td><td>84.044</td><td>110.352</td><td>69.863</td><td>121.177</td><td>88.852</td><td>88.139</td><td>89.568</td><td>87.424</td><td>90.280</td></tr><tr><td>ω<sub>c</sub> (deg)</td><td>94.59</td><td>77.54</td><td>118.47</td><td>57.28</td><td>170.61</td><td>95.61</td><td>77.88</td><td>120.71</td><td>56.41</td><td>172.77</td></tr><tr><td>Ω<sub>c</sub> (deg)</td><td>36.471</td><td>21.411</td><td>51.757</td><td>4.743</td><td>69.401</td><td>32.016</td><td>31.307</td><td>32.726</td><td>30.601</td><td>33.435</td></tr><tr><td>T<sub>0c</sub></td><td>2013.17</td><td>2013.00</td><td>2013.36</td><td>2012.73</td><td>2013.79</td><td>2013.18</td><td>2013.00</td><td>2013.38</td><td>2012.71</td><td>2013.81</td></tr><tr><td>P<sub>c</sub> (yr)</td><td>3.39</td><td>3.36</td><td>3.41</td><td>3.34</td><td>3.44</td><td>3.39</td><td>3.36</td><td>3.41</td><td>3.34</td><td>3.44</td></tr><tr><td>M<sub>c</sub> (M<sub>Jup</sub>)</td><td>9.37</td><td>8.44</td><td>10.43</td><td>7.56</td><td>11.79</td><td>9.18</td><td>8.31</td><td>10.14</td><td>7.47</td><td>11.31</td></tr><tr><td>M<sub>*</sub></td><td>1.76</td><td>1.73</td><td>1.79</td><td>1.70</td><td>1.82</td><td>1.76</td><td>1.74</td><td>1.79</td><td>1.71</td><td>1.82</td></tr><tr><td>α<sub>H0</sub><sup>*</sup> (mas)</td><td>-0.095</td><td>-0.253</td><td>0.065</td><td>-0.410</td><td>0.228</td><td>-0.025</td><td>-0.151</td><td>0.101</td><td>-0.275</td><td>0.228</td></tr><tr><td>δ<sub>H0</sub> (mas)</td><td>0.021</td><td>-0.143</td><td>0.183</td><td>-0.306</td><td>0.346</td><td>-0.021</td><td>-0.172</td><td>0.128</td><td>-0.322</td><td>0.277</td></tr><tr><td>π (mas)</td><td>51.413</td><td>51.281</td><td>51.544</td><td>51.149</td><td>51.676</td><td>51.396</td><td>51.265</td><td>51.526</td><td>51.134</td><td>51.657</td></tr><tr><td>μ<sub>α*</sub></td><td>4.965</td><td>4.947</td><td>4.981</td><td>4.928</td><td>4.997</td><td>4.963</td><td>4.948</td><td>4.977</td><td>4.932</td><td>4.992</td></tr><tr><td>μ<sub>δ</sub></td><td>83.967</td><td>83.948</td><td>83.986</td><td>83.925</td><td>84.004</td><td>83.969</td><td>83.950</td><td>83.987</td><td>83.930</td><td>84.004</td></tr><tr><td>γ (m s<sup>-1</sup>)</td><td>-26.1</td><td>-51.9</td><td>-3.7</td><td>-88.7</td><td>19.1</td><td>-25.3</td><td>-51.2</td><td>-2.8</td><td>-86.9</td><td>20.0</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>**Figure 16.** Comparison of posteriors on the orbital parameters on the inner planet, β Pic c, with unconstrained mutual inclination angle (blue) and a coplanar fit (red). While the two fits differ greatly in the derived inclination angle and position angle of nodes, the other parameters are very similar. The coplanar fit favors slightly smaller planet masses: $M_c = 9.4^{+1.1}_{-0.9} M_{Jup}$ for the unconstrained mutual inclination fit, and $M_c = 9.2^{+1.0}_{-0.9} M_{Jup}$ for the coplanar fit.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 17.** The parameters of the outer planet β Pic b, for the unconstrained mutual inclination fit (blue), the coplanar fit (red), and the regular Case 3 fit assuming only one planet in the system (black). The addition of the radial velocities and a second planet push the mass of β Pic b to lower values, in addition to slightly decreasing eccentricity, period, inclination angle, and position angle of nodes.</figcaption></figure>

β Pic due to the inner planet, as this orbit fit predicts significant acceleration of the star.

### 3.5. Comparison to Evolutionary Models

Chilcote et al. (2017) analyzed the SED of β Pic b, and found a model-dependent mass of $12.9 \pm 0.2 M_{Jup}$ using the
<figure><img src="image_22.png" /><figcaption>**Figure 18.** Predicted proper motion from our two-planet coplanar fit. As in the one-planet fit (Figure 9), the *Gaia* proper motion (which is not included in our fit) is a ~2σ outlier in R.A. The much shorter period of the oscillation, coupled with the larger amplitude, indicates that future *Gaia* data releases could detect the astrometric signature of β Pic c.</figcaption></figure>

bolometric luminosity of the planet, though this error bar does not include model uncertainty. Figure 19 compares the luminosity determined by Chilcote et al. (2017), $log L/L_⊙ = -3.76 \pm 0.02$, to predictions from the COND (Baraffe et al. 2003) and Sonora (M. Marley et al. 2020 in preparation) model grids, as well as the predicted luminosity given our dynamical mass measurement. As expected, the Chilcote et al. (2017) luminosity is significantly more precise than the uncertainty on our model-based prediction, given the ~30% errors on the dynamical mass; nevertheless, the estimate is consistent with the measurement. We compare our dynamical mass measurement to model predictions from this luminosity for the COND, Sonora, and SB12 (Spiegel & Burrows 2012) model grids in Figure 20, showing that the model-dependent luminosity estimates are well within the range implied from our one-planet fit mass measurement of $12.8^{+3.5}_{-3.2} M_{Jup}$. While the three model-dependent mass estimates all have exquisite precision relative to the dynamical mass, they are in significant disagreement with one another. The hot-start models (COND and Sonora) predict a significantly lower mass for the planet than the highest entropy models from the Spiegel & Burrows (2012) warm-start grid, and so the hot-start mass PDFs reach a maximum closer to the peak of the dynamical mass PDF than does the warm-start PDF, though all three model PDFs have peaks within the $1\sigma$ range of our dynamical mass measurement. The two-planet fit mass for β Pic b is significantly lower, with less than 5% of orbits corresponding to a mass larger than 12.5 $M_{Jup}$, more in tension with the model masses.

*Gaia* DR 3 proper motions and accelerations, along with continued monitoring of the relative orbit by direct imaging, will likely further constrain the orbit and the mass, and the DR 4 intermediate data will allow for a full fit including individual absolute astrometric measurements from Hipparcos and Gaia, as well as ground-based relative astrometry and radial

<figure><img src="image_23.png" /><figcaption>**Figure 19.** Comparison of the luminosity of β Pic b (Chilcote et al. 2017) and the age (Nielsen et al. 2016) to the Sonora (M. Marley et al. 2020 in preparation) and COND (Baraffe et al. 2003) models. Luminosity is given in solar units, and red and blue numbers mark the masses of the tracks in $M_{Jup}$ for the Sonora and COND models, respectively. Light red shaded region represents the $1\sigma$ region for the mass of the planet from our one-planet Case 3 fit and the Sonora models, consistent with the expected mass given the luminosity measurement.</figcaption></figure>

velocities. A precise determination of the mass of the planet using these data will allow β Pic b to be used as an empirical calibrator for evolutionary models at young ages where planets are still significantly radiating away their formative heat. We note that the luminosity-derived masses discussed previously assume prompt planet formation. A delay between star and planet formation may lead to a significantly younger age for β Pic b than its host star (e.g., Currie et al. 2009). Our current constraints on the dynamical mass of the planet do not allow us to distinguish between a prompt and delayed formation scenario assuming a given evolutionary model. The $8.0 \pm 2.6 M_{Jup}$ mass from the two-planet fit would require a significantly delayed epoch of planet formation to bring the luminosity in line with evolutionary models. A precise, model-dependent measurement of the entropy of formation will greatly constrain formation models for wide-separation giant planets as well.

In the meantime, more ground-based relative astrometry will also increase the mass precision. Figure 7 shows significant covariance between eccentricity, period, and planet mass. Thus, further constraints on the orbital parameters will reduce the mass errors. Figure 21 shows a significant divergence in the orbit tracks beyond ~2022, with higher masses generally corresponding to the shortest orbital periods.

To highlight this dependence, in Figure 22 we subtract off the lowest $\chi^2$ orbit from each of the tracks. The prediction for separation at 2020.0 has an uncertainty of 1.8 mas, which rises to 3.5 mas at 2021.0, and 8.2 mas at 2022.0. In comparison, Wang et al. (2016) demonstrated the ability to reach relative astrometric precision of less than 1 mas on β Pic b with GPI when the separation was above ~230 mas, a separation the planet should have reached again in 2019 June. Thus,
<figure><img src="image_24.png" /></figure>

**Figure 20.** Mass posteriors for β Pic b from our orbit fit for the one-planet Case 3 fit (green solid curve) and the two-planet coplanar fit (green dotted line), compared to model-derived masses from COND (Baraffe et al. 2003), Sonora (M. Marley et al. 2020 in preparation), and SB12 (Spiegel & Burrows 2012) (hybrid clouds, solar metallicity) based on the luminosity and age of the planet. Our one-planet fit dynamical mass is consistent with the model predictions, though the current precision does not allow us to differentiate between the different models. The two-planet fit mass is more discrepant with these model predictions, with less than a 5% probability that the mass is larger $12.5 M_{Jup}$.

<figure><img src="image_25.png" /></figure>

**Figure 21.** Future tracks of the Case 3, one-planet fit β Pic b orbit, again color-coded by mass, drawn from the posterior of the orbit including Hipparcos, Gaia, and all imaging data except SPHERE. There is a general trend where higher-mass planets result in a faster turnaround in ~2024.

continued monitoring with GPI and SPHERE will further refine the orbit determination. While there is not a direct correlation between mass and separation between 2020 and 2022, if the

<figure><img src="image_26.png" /></figure>

**Figure 22.** Same as Figure 21, but with the lowest $\chi^2$ orbit subtracted from the tracks to give more detail. Further monitoring of the system between 2020 and 2022 at the 1 mas level will greatly reduce the uncertainty in the orbital parameters, particularly in period.

shortest orbital periods are ruled out, this will also rule out the largest values of planet mass, leading to a more precise mass measurement.

### 3.6. Disk and Hill Sphere

Our new constraints for the orbital geometry of β Pic b are also relevant to the ongoing investigations of planet–disk dynamical interactions. Periastron occurs near maximum elongation, close to the sky plane in the SW at epoch 2013.72. With $a = 10.18$ au and $e = 0.122$, the projected periastron and apastron separations are 8.94 au and 11.42 au, respectively. If the planet is secularly forcing the eccentricities of the nearby material (Wyatt et al. 1999), then the inner cavity cleared by the planet should have a stellocentric offset similar to that of Fomalhaut’s dust belt (Kalas et al. 2005), i.e., roughly 4.18 au or 215 mas. The current scattered light data (Golimowski et al. 2006; Apai et al. 2015) and millimeter continuum maps (Matrà et al. 2019) do not have the required angular resolution to directly detect the hypothesized offset. However, the ~20% stronger mid-infrared thermal emission from the SW side of the disk compared to the NE is consistent with the offset (Lagage & Pantin 1994; Wahhaj et al. 2003).

Given the edge-on nature of the β Pic debris disk and the planet’s orbit, as well as evidence for a transit-like event in 1981 (Lecavelier Des Etangs & Vidal-Madjar 2009), determining whether the planet transits has become of great interest (e.g., Nielsen et al. 2014; Wang et al. 2016). Wang et al. (2016) ruled out the prospect of the planet itself transiting at $10\sigma$, and from our Case 3, we find the closest approach by the planet to be a projected separation of $22.7 \pm 1.9 R_*$, where we take the radius of the star to be $1.8 R_\odot$, from the interferometric measurement of Di Folco et al. (2004) of $1.8 \pm 0.2 R_\odot$ (Figure 23). From our MCMC chain, the smallest projected separation reached is $13.6 R_*$ (0.11 au, 0.09 $r_H$), similar to all of our other cases (minimum values of 13.5–14.3 $R_*$) except
<figure><img src="image_27.png" /><figcaption>**Figure 23.** Posteriors on the Hill sphere crossing and conjunction (closest approach) between 2017 and mid-2018, for orbit fits including data from 2018. For most models, the Hill sphere crosses in front of the star in 2017 mid-April, and the crossing lasts until 2018 early-February. In all fits including the 2018 GPI data, conjunction occurs in a ten-day window between 2017 September 11 and 2017 September 20 ($2\sigma$).</figcaption></figure>

Case 7, which utilized only imaging data sets from NACO and SPHERE. In this orbit fit, we find a minimum projected separation of $17.3 \pm 3.1 R_*$, with a single value out of $10^6$ having a separation $<5 R_*$, at $4.6 R_*$. Thus, we concur with Wang et al. (2016) that the astrometry strongly disfavors transit, at the $12.2\sigma$ level for Case 3. This conclusion is the same for the two-planet fit and the three-planet fit, with the smallest projected separation barely changing, from $22.7 \pm 1.9 R_*$ for Case 3 to $23.0 \pm 1.9 R_*$ for both the unconstrained mutual inclination and coplanar fits.

While transit of the planet itself is ruled out, the Hill sphere of the planet passes in front of the star, as noted by Wang et al. (2016). From 2017 to 2018, this offered a rare opportunity to probe the circumplanetary environment of a young Jovian exoplanet at large orbital separations where the influence of the star is minimal. There have been numerous observational efforts to monitor $\beta$ Pic both photometrically and spectroscopically during this Hill sphere crossing (e.g., Měkarnia et al. 2017; Stuik et al. 2017; Kalas et al. 2019; Mellon et al. 2019), so pinpointing the timeframe of the crossing is of prime interest to put these monitoring programs into context.

The Hill sphere radius ($r_H$) is given by

$$r_H \approx a(1 - e) \left(\frac{M_P}{3M_*}\right)^{\frac{1}{3}}, \qquad (12)$$

a function of the semimajor axis ($a$), eccentricity ($e$), planet mass ($M_P$), and star mass ($M_*$), from Hamilton & Burns (1992). For our Case 3 orbit fits, we find a value of the Hill sphere radius of $1.18_{-0.11}^{+0.15}$ au. We find that the Hill sphere first crosses in front of the star 2017 April 11 ($\pm 18$ days), and the crossing lasts until 2018 February 16 ($\pm 18$ days). Conjunction

is more tightly constrained, taking place on 2017 September 13 ($\pm 2.8$ days). The uncertainty in the timing of the Hill sphere crossings is dominated by the error on the planet mass; fixing the planet mass results in timing windows $\sim 6 \times$ smaller, more in line with the precision on time of conjunction. The predictions for Hill Sphere crossing do change in the two-planet fit, largely because the size of the Hill sphere is reduced thanks to a smaller planet mass for $\beta$ Pic b, starting in 2017 May 5 ($\pm 18$ days) and ending 2018 January 24 ($\pm 18$ days). Compared to the one-planet Case 3 fit, time of conjunction is similar for both two-planet fits. We find conjunction to take place on 2017 September 15 ($\pm 2.9$ days) for the unconstrained mutual inclination fit, and 2017 September 14 ($\pm 2.8$ days) for the coplanar fit.

The dates of these events vary from the different orbit fits, as shown in Figure 23. For Case 2, where the planet mass is not determined from the fit, we randomly sample from the Case 3 planet mass posterior. Generally, posteriors are similar for data sets that include the GPI 2018 astrometry (Cases 2, 3, 5, and 6). For example, time of conjunction has a median date of 2017 September 12 and 2017 September 13 for Cases 2 and 3, 2017 September 17 for Cases 5 and 6, and 2017 September 28 for Case 7. The large offset for the Case 7 fit, using only NACO and SPHERE data, suggests there is a bias between GPI and SPHERE relative astrometry, either due to instrumental calibration or data pipeline systematics. Indeed, for time of conjunction, the fits combining GPI and SPHERE data are in between fits excluding one of the two instruments. Our Case 6 combines data from the two instruments with offset terms in separation and position angle, but the posteriors on time of conjunction and Hill sphere crossings are identical whether these offsets are applied (case 6) or not (Case 5). Thus, a single offset between the two instruments does not appear to address the issue, suggesting that either a different parameterization is needed or the bias is time-variable. Further work is needed to understand how this offset arises between the two instruments. M. Nguyen et al. (2020 submitted) are currently analyzing multiple epochs of the same calibration field taken over GPI's lifetime as a validation of the astrometric calibration presented in De Rosa et al. (2019).

## 4. Conclusion

We combine relative astrometry of the planet $\beta$ Pic b with Gaia position and Hipparcos IAD to refine the orbit and measure the mass of the planet. We find a model-independent mass for the planet of $12.8_{-3.2}^{+5.5} M_{Jup}$, consistent with predictions from hot-start evolutionary models given the luminosity of the planet and age of the system. We find significant evidence for nonzero, but low, eccentricity for the planet, finding a value of $0.12_{-0.03}^{+0.04}$. Our comparison to previous work by Snellen & Brown (2018) and Dupuy et al. (2019) underscores the importance of performing a joint fit to the space-based absolute astrometry and ground-based relative astrometry. The reason for the offset between our fit and that of Dupuy et al. (2019) is less clear, and could be a combination of new relative astrometry, their fitting additional RV data of the star and planet, and their use of recalculated Hipparcos and Gaia catalog values.

When including the radial velocities of the star from Lagrange et al. (2019b) and adding an additional planet to the fit, $\beta$ Pic c, we find a significantly lower mass for $\beta$ Pic b,
$8 \pm 2.6 M_{\text{Jup}}$, and no significant difference in the orbital parameters whether the planets are assumed to be coplanar or not. We predict significant astrometric motion of the star from the orbit of $\beta$ Pic c, and future Gaia data releases may be able to detect the signature of this inner planet.

We have constrained the time of conjunction of the planet to an accuracy of 2.7 days, and the Hill sphere entrance and exit to 18 days. These values will guide analysis of the photometric monitoring of the star over the last two years to search for circumbinary material transiting in front of the star.

Future monitoring of $\beta$ Pic by both ground-based imaging and Gaia should further improve the precision on the measurement of the planet mass. As the planet moves further from the star, GPI and SPHERE will be able to determine the relative astrometry with increasing precision. Similarly, if Gaia astrometry on bright stars can be improved, the reflex motion of the star over the Gaia mission can be used to directly constrain the planet mass.

The combination of directly imaged short-period substellar companions and precision Gaia astrometry represents a new opportunity to directly measure the masses of these objects. As shown in the case of $\beta$ Pic, strong constraints on the orbital parameters allow us to connect the motion from the ~1991 Hipparcos IAD directly to the ~2015 Gaia astrometry. Further Gaia data releases and ground-based imaging will allow us to measure, or set upper limits on, other directly imaged substellar companions with shorter orbital periods, including 51 Eridani b (Macintosh et al. 2015), HR 2562 (Konopacky et al. 2016), and HD 984 (Meshkat et al. 2015). When coupled with James Webb Space Telescope (JWST) mid-IR observations of these objects that will sample the part of the SED with the bulk of the flux, we can directly compare the luminosity predictions of evolutionary models to the measured masses from absolute and relative astrometry.

This work is based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), and Ministério da Ciência, Tecnologia e Inovação (Brazil). This research has made use of the SIMBAD and VizieR databases, operated at CDS, Strasbourg, France. This work has made use of data from the European Space Agency (ESA) 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. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562.

J.R. and R.D. acknowledge support from the Fonds de Recherche du Québec. J.R.M.'s work was performed in part under contract with the California Institute of Technology

(Caltech)/Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. Support for M.M.B.'s work was provided by NASA through Hubble Fellowship grant #51378.01-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. Supported by NSF grants AST-1411868 (E.L. N., K.B.F., B.M., J.P., and J.H.), AST-141378 (G.D.), and AST-1518332 (R.D.R., J.J.W., T.M.E., J.R.G., P.G.K.). Supported by NASA grants NNX14AJ80G (E.L.N., B.M., F. M., and M.P.), NNX15AC89G and NNX15AD95G (B.M., T. M.E., R.J.D.R., G.D., J.R.G., P.G.K.), NN15AB52I (D.S.), NSSC17K0535 (E.L.N., R.D.R., B.M., J.B.R.). J.J.W. is supported by the Heising-Simons Foundation 51 Pegasi b postdoctoral fellowship. Portions of this work were performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. This work benefited from NASA's Nexus for Exoplanet System Science (NExSS) research coordination network, sponsored by NASA's Science Mission Directorate.

*Facility:* Gemini:South (GPI).

*Software:* Astropy (Astropy Collaboration et al. 2013), emcee (Foreman-Mackey et al. 2013), GPI DRP (Perrin et al. 2016), IDL Astronomy Library (Landsman 1993), pyKLIP (Wang et al. 2015).

## ORCID iDs

Eric L. Nielsen https://orcid.org/0000-0001-6975-9056

Robert J. De Rosa https://orcid.org/0000-0002-4918-0247

Jason J. Wang https://orcid.org/0000-0003-0774-6502

Johannes Sahlmann https://orcid.org/0000-0001-9525-3673

Paul Kalas https://orcid.org/0000-0002-6221-5360

Gaspard Duchêne https://orcid.org/0000-0002-5092-6464

Julien Rameau https://orcid.org/0000-0003-0029-0258

Mark S. Marley https://orcid.org/0000-0002-5251-2943

Didier Saumon https://orcid.org/0000-0001-6800-3505

Bruce Macintosh https://orcid.org/0000-0003-1212-7538

Maxwell A. Millar-Blanchaer https://orcid.org/0000-0001-6205-9233

Meiji M. Nguyen https://orcid.org/0000-0002-9350-4763

S. Mark Ammons https://orcid.org/0000-0001-5172-7902

Vanessa P. Bailey https://orcid.org/0000-0002-5407-2806

Travis Barman https://orcid.org/0000-0002-7129-3002

Joanna Bulger https://orcid.org/0000-0003-4641-2003

Jeffrey Chilcote https://orcid.org/0000-0001-6305-7272

Tara Cotten https://orcid.org/0000-0003-0156-3019

Thomas M. Esposito https://orcid.org/0000-0002-0792-3719

Michael P. Fitzgerald https://orcid.org/0000-0002-0176-8973

Katherine B. Follette https://orcid.org/0000-0002-7821-0695

Benjamin L. Gerard https://orcid.org/0000-0003-3978-9195

Stephen J. Goodsell https://orcid.org/0000-0002-4144-5116

Alexandra Z. Greenbaum https://orcid.org/0000-0002-7162-8036

Pascale Hibon https://orcid.org/0000-0003-3726-5494

Li-Wei Hung https://orcid.org/0000-0003-1498-6088

Patrick Ingraham https://orcid.org/0000-0003-3715-8138

Quinn Konopacky https://orcid.org/0000-0002-9936-6285

James E. Larkin https://orcid.org/0000-0001-7687-3965
Franck Marchis https://orcid.org/0000-0001-7016-7277

Christian Marois https://orcid.org/0000-0002-4164-4182

Stanimir Metchev https://orcid.org/0000-0003-3050-8203

Rebecca Oppenheimer https://orcid.org/0000-0001-7130-7681

David Palmer https://orcid.org/0000-0001-7128-0802

Marshall Perrin https://orcid.org/0000-0002-3191-8151

Abhijith Rajan https://orcid.org/0000-0002-9246-5467

Fredrik T. Rantakyrö https://orcid.org/0000-0002-9667-2244

Jean-Baptiste Ruffio https://orcid.org/0000-0003-2233-4821

Dmitry Savransky https://orcid.org/0000-0002-8711-7206

Adam C. Schneider https://orcid.org/0000-0002-6294-5937

Anand Sivaramakrishnan https://orcid.org/0000-0003-1251-4124

Inseok Song https://orcid.org/0000-0002-5815-7372

Remi Soummer https://orcid.org/0000-0003-2753-2819

Sandrine Thomas https://orcid.org/0000-0002-9121-3436

J. Kent Wallace https://orcid.org/0000-0001-5299-6899

Kimberly Ward-Duong https://orcid.org/0000-0002-4479-8291

Sloane Wiktorowicz https://orcid.org/0000-0003-4483-5037

Schuyler Wolff https://orcid.org/0000-0002-9977-8255

## References

Apai, D., Schneider, G., Grady, C. A., et al. 2015, ApJ, 800, 136

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701

Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42

Barrado y Navascués, D., Stauffer, J. R., Song, I., & Caillault, J.-P. 1999, ApJL, 520, L123

Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, MNRAS, 454, 593

Binks, A. S., & Jeffries, R. D. 2014, MNRAS, 438, L11

Brandt, T. D., Dupuy, T., & Bowler, B. P. 2019, AJ, 158, 140

Butkevich, A. G., & Lindegren, L. 2014, A&A, 570, A62

Chauvin, G., Lagrange, A.-M., Beust, H., et al. 2012, A&A, 542, A41

Chilcote, J., Pueyo, L., De Rosa, R. J., et al. 2017, AJ, 153, 182

Currie, T., Lada, C. J., Plavchan, P., et al. 2009, ApJ, 698, 1

Currie, T., Thalmann, C., Matsumura, S., et al. 2011, ApJL, 736, L33

De Rosa, R. J., Nguyen, M. M., Chilcote, J., et al. 2019, arXiv:1910.08659

De Rosa, R. J., Nielsen, E. L., Blunt, S. C., et al. 2015, ApJL, 814, L3

Di Folco, E., Thévenin, F., Kervella, P., et al. 2004, A&A, 426, 601

Dupuy, T. J., Brandt, T. D., Kratter, K. M., & Bowler, B. P. 2019, ApJL, 871, L4

ESA 1997, The Hipparcos and Tycho Catalogues, ESA SP-1200

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

Golimowski, D. A., Ardila, D. R., Krist, J. E., et al. 2006, AJ, 131, 3109

Hamilton, D. P., & Burns, J. A. 1992, Icar, 96, 43

Kalas, P., Graham, J. R., & Clampin, M. 2005, Natur, 435, 1067

Kalas, P., & Jewitt, D. 1995, AJ, 110, 794

Kalas, P., Zwintz, K., Kenworthy, M., et al. 2019, AAS Meeting Abstracts, 233, 218.03

Kervella, P., Arenou, F., Mignard, F., & Thévenin, F. 2019, A&A, 623, A72

Konopacky, Q. M., Rameau, J., Duchêne, G., et al. 2016, ApJL, 829, L4

Lagage, P. O., & Pantin, E. 1994, Natur, 369, 628

Lagrange, A. M., Boccaletti, A., Langlois, M., et al. 2019a, A&A, 621, L8

Lagrange, A.-M., Boccaletti, A., Milli, J., et al. 2012, A&A, 542, A40

Lagrange, A.-M., Bonnefoy, M., Chauvin, G., et al. 2010, Sci, 329, 57

Lagrange, A.-M., Gratadour, D., Chauvin, G., et al. 2009, A&A, 493, L21

Lagrange, A. M., Meunier, N., Rubini, P., et al. 2019b, NatAs, 3, 1135

Landsman, W. B. 1993, in ASP Conf. Ser. 52, Astronomical Data Analysis Software and Systems II, ed. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes (San Francisco, CA: ASP), 246

Lecavelier Des Etangs, A., & Vidal-Madjar, A. 2009, A&A, 497, 557

Lindegren, L., Hernandez, J., Bombrun, A., et al. 2018, A&A, 616, A2

Macintosh, B., Graham, J. R., Barman, T., et al. 2015, Sci, 350, 64

Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, PNAS, 111, 12661

Maire, A.-L., Rodet, L., Cantalloube, F., et al. 2019, A&A, 624, A118

Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006a, ApJ, 641, 556

Marois, C., Lafrenière, D., Macintosh, B., & Doyon, R. 2006b, ApJ, 647, 612

Martín-Fleitas, J., Sahlmann, J., Mora, A., et al. 2014, Proc. SPIE, 9143, 91430Y

Matrà, L., Wyatt, M. C., Wilner, D. J., et al. 2019, AJ, 157, 135

Mékarnia, D., Chapellier, E., Guillot, T., et al. 2017, A&A, 608, L6

Mellon, S. N., Stuik, R., Kenworthy, M., et al. 2019, AAS Meeting Abstracts, 233, 140.22

Meshkat, T., Bonnefoy, M., Mamajek, E. E., et al. 2015, MNRAS, 453, 2379

Michalik, D., Lindegren, L., Hobbs, D., & Lammers, U. 2014, A&A, 571, A85

Millar-Blanchaer, M. A., Graham, J. R., Pueyo, L., et al. 2015, ApJ, 811, 18

Morzinski, K. M., Males, J. R., Skemer, A. J., et al. 2015, ApJ, 815, 108

Nielsen, E. L., De Rosa, R. J., Rameau, J., et al. 2017, AJ, 154, 218

Nielsen, E. L., Liu, M. C., Wahhaj, Z., et al. 2014, ApJ, 794, 158

Nielsen, E. L., Rosa, R. J. D., Wang, J., et al. 2016, AJ, 152, 175

Perrin, M. D., Ingraham, P., Follette, K. B., et al. 2016, Proc. SPIE, 9908, 990837

Perrin, M. D., Maire, J., Ingraham, P., et al. 2014, Proc. SPIE, 9147, 91473J

Pueyo, L., Soummer, R., Hoffmann, J., et al. 2015, ApJ, 803, 31

Sahlmann, J., Lazorenko, P. F., Ségransan, D., et al. 2013, A&A, 556, A133

Sahlmann, J., Lazorenko, P. F., Ségransan, D., et al. 2016a, A&A, 595, A77

Sahlmann, J., Martín-Fleitas, J., Mora, A., et al. 2016b, Proc. SPIE, 9904, 99042E

Sahlmann, J., Ségransan, D., Queloz, D., et al. 2011, A&A, 525, A95

Sivaramakrishnan, A., & Oppenheimer, B. R. 2006, ApJ, 647, 620

Smith, B. A., & Terrile, R. J. 1984, Sci, 226, 1421

Snellen, I. A. G., Brandl, B. R., de Kok, R. J., et al. 2014, Natur, 509, 63

Snellen, I. A. G., & Brown, A. G. A. 2018, NatAs, 2, 883

Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJL, 755, L28

Sparks, W. B., & Ford, H. C. 2002, ApJ, 578, 543

Spiegel, D. S., & Burrows, A. 2012, ApJ, 745, 174

Stuik, R., Bailey, J. I., Dorval, P., et al. 2017, A&A, 607, A45

van Leeuwen, F. 2007a, A&A, 474, 653

van Leeuwen, F. 2007b, ASSL, 350, 1

Wahhaj, Z., Koerner, D. W., Ressler, M. E., et al. 2003, ApJL, 584, L27

Wang, J. J., Graham, J. R., Pueyo, L., et al. 2016, AJ, 152, 97

Wang, J. J., Perrin, M. D., Savransky, D., et al. 2018, JATIS, 4, 018002

Wang, J. J., Rajan, A., Graham, J. R., et al. 2014, Proc. SPIE, 9147, 914755

Wang, J. J., Ruffio, J.-B., De Rosa, R. J., et al. 2015, pyKLIP: PSF Subtraction for Exoplanets and Disks, Astrophysics Source Code Library, ascl:1506.001

Weinberger, A. J., Becklin, E. E., & Zuckerman, B. 2003, ApJL, 584, L33

Wolff, S. G., Perrin, M. D., Maire, J., et al. 2014, Proc. SPIE, 9147, 91477H

Wyatt, M. C., Dermott, S. F., Telesco, C. M., et al. 1999, ApJ, 527, 918

Zuckerman, B., Song, I., Bessell, M. S., & Webb, R. A. 2001, ApJL, 562, L87