# The Orbit and Transit Prospects for β Pictoris b
Constrained with One Milliarcsecond Astrometry

Jason J. Wang, James R. Graham, Laurent Pueyo, Paul Kalas, Maxwell A.

Millar-Blanchaer, Jean-Baptiste Ruffio, Robert J. de Rosa, S. Mark Ammons,
Pauline Arriaga, Vanessa P. Bailey, et al.

► To cite this version:

Jason J. Wang, James R. Graham, Laurent Pueyo, Paul Kalas, Maxwell A. Millar-Blanchaer, et al.,

The Orbit and Transit Prospects for β Pictoris b Constrained with One Milliarcsecond Astrometry.

The Astronomical Journal, 2016, 152, 10.3847/0004-6256/152/4/97 . insu-03691493

HAL Id: insu-03691493

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

Submitted on 7 Sep 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_1.png" /></figure>

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

THE ORBIT AND TRANSIT PROSPECTS FOR β PICTORIS b CONSTRAINED WITH ONE MILLIARCSECOND ASTROMETRY

JASON J. WANG¹, JAMES R. GRAHAM¹, LAURENT PUEYO², PAUL KALAS¹, MAXWELL A. MILLAR-BLANCHAER³,
JEAN-BAPTISTE RUFFIO⁴, ROBERT J. DE ROSA¹, S. MARK AMMONS⁵, PAULINE ARRIAGA⁶, VANESSA P. BAILEY⁴, TRAVIS S. BARMAN⁷,
JOANNA BULGER⁸, ADAM S. BURROWS⁹, ANDREW CARDWELL¹⁰,¹¹, CHRISTINE H. CHEN², JEFFREY K. CHILCOTE¹², TARA COTTEN¹³,
MICHAEL P. FITZGERALD⁶, KATHERINE B. FOLLETTE⁴, RENÉ DOYON¹⁴, GASPARD DUCHÊNE¹,¹⁵, ALEXANDRA Z. GREENBAUM²,¹⁶,
PASCALE HIBON¹⁷, LI-WEI HUNG⁶, PATRICK INGRAHAM¹⁸, QUINN M. KONOPACKY¹⁹, JAMES E. LARKIN⁶, BRUCE MACINTOSH⁴,
JÉRÔME MAIRE¹², FRANCK MARCHIS²⁰, MARK S. MARLEY²¹, CHRISTIAN MAROIS²²,²³, STANIMIR METCHEV²⁴,²⁵,
ERIC L. NIELSEN⁴,²⁰, REBECCA OPPENHEIMER²⁶, DAVID W. PALMER⁵, RAHUL PATEL²⁷, JENNY PATIENCE²⁸, MARSHALL D. PERRIN²,
LISA A. POYNEER⁵, ABHIJITH RAJAN²⁸, JULIEN RAMEAU¹⁴, FREDRIK T. RANTAKYRÖ⁹, DMITRY SAVRANSKY²⁹,
ANAND SIVARAMAKRISHNAN², INSEOK SONG¹³, REMI SOUMMER², SANDRINE THOMAS¹⁸, GAUTAM VASISHT³⁰, DAVID VEGA²⁰,
J. KENT WALLACE³⁰, KIMBERLY WARD-DUONG²⁸, SLOANE J. WIKTOROWICZ³¹, AND SCHUYLER G. WOLFF²,¹⁶

¹ Astronomy Department, University of California, Berkeley, Berkeley, CA 94720, USA; jwang@astro.berkeley.edu

² Space Telescope Science Institute, 3700 San Martin Drive, Baltimore MD 21218, USA

³ Department of Astronomy & Astrophysics, University of Toronto, Toronto ON M5S 3H4, Canada

⁴ Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, Stanford, CA 94305, USA

⁵ Lawrence Livermore National Lab, 7000 East Avenue, Livermore, CA 94551, USA

⁶ Department of Physics and Astronomy, UCLA, Los Angeles, CA 90095, USA

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

⁸ Subaru Telescope, NAOJ, 650 North Aohoku Place, Hilo, HI 96720, USA

⁹ Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

¹⁰ Gemini Observatory, Casilla 603, La Serena, Chile

¹¹ Large Binocular Telescope Observatory, 933 North Cherry Avenue, Tucson AZ 85721, USA

¹² Dunlap Institute for Astronomy & Astrophysics, University of Toronto, 50 St. George Street, Toronto ON M5S 3H4, Canada

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

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

¹⁵ Université Grenoble Alpes/CNRS, Institut de Planétologie et d'Astrophysique de Grenoble, F-38000 Grenoble, France

¹⁶ Physics & Astronomy Department, Johns Hopkins University, 3600 North Charles Street, Baltimore, MD 21218, USA

¹⁷ European Southern Observatory, Alonso de Cordova 3107, Vitacura, Santiago, Chile

¹⁸ Large Synoptic Survey Telescope, 950 North Cherry Avenue, Tucson AZ 85719, USA

¹⁹ Center for Astrophysics and Space Sciences, University of California, San Diego, La Jolla, CA 92093, USA

²⁰ SETI Institute, Carl Sagan Center, 189 Bernardo Avenue, Mountain View, CA 94043, USA

²¹ NASA Ames Research Center, MS 245-3, Moffett Field, CA, 94035, USA

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

²³ University of Victoria, 3800 Finnerty Road, Victoria, BC, V8P 5C2, Canada

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

²⁵ Stony Brook University, 100 Nicolls Road, Stony Brook, NY 11794, USA

²⁶ American Museum of Natural History, New York, NY 10024, USA

²⁷ Infrared Processing and Analysis Center, California Institute of Technology, 700 South Wilson Avenue, Pasadena, CA 91106, USA

²⁸ School of Earth and Space Exploration, Arizona State University, PO Box 871404, Tempe, AZ 85287, USA

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

³⁰ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena CA 91109, USA

³¹ The Aerospace Corporation, 2310 East El Segundo Boulevard, El Segundo, CA 90245, USA

Received 2016 May 27; revised 2016 July 1; accepted 2016 July 16; published 2016 October 3

ABSTRACT

A principal scientific goal of the Gemini Planet Imager (GPI) is obtaining milliarcsecond astrometry to constrain exoplanet orbits. However, astrometry of directly imaged exoplanets is subject to biases, systematic errors, and speckle noise. Here, we describe an analytical procedure to forward model the signal of an exoplanet that accounts for both the observing strategy (angular and spectral differential imaging) and the data reduction method (Karhunen–Loève Image Projection algorithm). We use this forward model to measure the position of an exoplanet in a Bayesian framework employing Gaussian processes and Markov-chain Monte Carlo to account for correlated noise. In the case of GPI data on β Pic b, this technique, which we call Bayesian KLIP-FM Astrometry (BKA), outperforms previous techniques and yields 1σ errors at or below the one milliarcsecond level. We validate BKA by fitting a Keplerian orbit to 12 GPI observations along with previous astrometry from other instruments. The statistical properties of the residuals confirm that BKA is accurate and correctly estimates astrometric errors. Our constraints on the orbit of β Pic b firmly rule out the possibility of a transit of the planet at 10-σ significance. However, we confirm that the Hill sphere of β Pic b will transit, giving us a rare chance to probe the circumplanetary environment of a young, evolving exoplanet. We provide an ephemeris for photometric monitoring of the Hill sphere transit event, which will begin at the start of April in 2017 and finish at the end of January in 2018.

Key words: astrometry – planets and satellites: individual (β Pictoris b) – techniques: image processing
# 1. INTRODUCTION

Astrometry is an essential tool for characterizing directly imaged exoplanets and their physical relationship to other elements of the planetary system in which they reside. The Gemini Planet Imager (GPI; Macintosh et al. 2014) was designed with a goal of achieving $\le 1.8$ mas astrometric accuracy (Graham 2009), which is necessary for characterizing the eccentricity distribution of exoplanet orbits from the GPI Exoplanet Survey (Konopacky et al. 2014). To do so, the astrometric calibration of GPI has continually been benchmarked to well calibrated astrometric fields (Konopacky et al. 2014). This had led to some of the most precise astrometry on directly imaged exoplanet systems to date (De Rosa et al. 2015; Millar-Blanchaer et al. 2015; Rameau et al. 2016), allowing us to constrain or fit the first ever orbit of some of these directly imaged exoplanets. However, limited by either the signal-to-noise ratio (S/N) of these exoplanets or by biases in the various data analysis algorithms, so far no astrometric study with GPI has achieved the design goal of 1.8 mas precision.

The importance of understanding planetary orbits is highlighted by the β Pictoris system, a young (~23 Myr; Mamajek & Bell 2014; Binks & Jeffries 2016) and nearby (19.3 pc; van Leeuwen 2007) system that has been extensively studied. β Pic harbors a near edge-on debris disk that was first imaged by Smith & Terrile (1984) and which was subsequently observed to have a warp thought to be induced by a planet whose orbit is inclined relative to the debris disk (Burrows et al. 1995; Mouillet et al. 1997; Heap et al. 2000). Additional indirect signatures of a planet were derived from variable spectral features modeled as infalling comets (Beust & Morbidelli 2000) and a peculiar light curve anomaly detected in 1981 (Lecavelier Des Etangs et al. 1997).

Lagrange et al. (2009, 2010) then discovered β Pic b, a planet at an appropriate mass (~10 $M_{\oplus}$) and semimajor axis (8–13 au) to be responsible for the previously observed indirect signatures of planets. A key focus of subsequent observations was determining the alignment of β Pic b relative to the main outer disk and the warp to determine if β Pic b is causing the warp. By observing the disk and planet simultaneously, Lagrange et al. (2012) concluded that the planet is misaligned from the main outer disk and consistent with being responsible for the warp. Additionally, Dawson et al. (2011) ruled out the possibility of having another giant planet in the system massive enough to cause the warp instead. Thus, β Pic b is responsible for the warp in the debris disk. This was confirmed in astrometric monitoring campaigns (Chauvin et al. 2012; NielSEN et al. 2014; Millar-Blanchaer et al. 2015) which used homogeneous data sets to limit systematics and constrain the orbit of β Pic b.

Refining the orbital elements of β Pic b is not only crucial for investigating the dynamical link between the planet and the disk warp, but also because β Pic b may transit its host star once every ~20 years. To date, there are no other known systems where the physical properties of an exoplanet can be characterized by using both the direct imaging and transit techniques. Currently, the tightest constraint on the probability of transit is ~0.06%, obtained with a dedicated astrometric monitoring campaign with GPI (Millar-Blanchaer et al. 2015). However, Lecavelier des Etangs & Vidal-Madjar (2016) point out that GPI measurements from Millar-Blanchaer et al. (2015)

have a higher position angle (PA) than astrometry from previous measurements, which could arise from a possible systematic calibration offset between GPI and previous instruments instead of actually due to β Pic b's orbit being slightly inclined away from edge on. We note there is currently no evidence of a PA offset with the GPI astrometric calibration, and the GPI astrometry of HD 95086 b are consistent with previous astrometry from other instruments (Rameau et al. 2016). However, it is important to more accurately compute the orbit of β Pic b because in late 2017 to early 2018 (Millar-Blanchaer et al. 2015) it will be at its closest projected separation from the star. The transit of the planet and/or any circumplanetary material orbiting around it could therefore be detectable.

One of the obstacles in characterizing directly imaged exoplanets is that even with the newest instrumentation, the glare of the host star covers the signal of the planet. In order to subtract the point-spread function (PSF) of the star and maximize the S/N of the planet, observing techniques such as Angular Differential Imaging (ADI; Marois et al. 2006a) and Spectral Differential Imaging (SDI; Marois et al. 2000) and data reduction algorithms like Karhunen–Loève (KL) Image Projection (KLIP; Soummer et al. 2012; Pueyo et al. 2015) are used in combination to disentangle the point-spread function (PSF) of the star from potential astrophysical sources. However, these techniques distort the planet signal and create data reduction artifacts, which are usually nuisance parameters that need to be calibrated out to obtain unbiased astrometry.

Forward modeling effects of observing techniques and data reduction algorithms on the PSF of the planet was first done in the context of ADI and LOCI, showing significant improvements in astrometry and photometry for simulated planets (Marois et al. 2010; Galicher & Marois 2011). In similar contexts with ADI and LOCI, Brandt et al. (2013) and Esposito et al. (2014) used forward modeling to correct for the flux loss of exoplanets and the flux and morphology of disks, respectively. For classical ADI (cADI), Cantalloube et al. (2015) showed that forward modeling techniques can improve the sensitivity of cADI and mitigate biases. The use of Markov-chain Monte Carlo (MCMC) in conjunction with forward modeling was presented in Bottom et al. (2015) for reference differential imaging.

Recently, Pueyo (2016) introduced a method called KLIP-FM to analytically forward model the degradation of a faint astrophysical signal that occurs when using least squares-based PSF subtraction algorithms such as KLIP that is also generally applicable to any observing strategy. Additionally, the computation of the forward model with KLIP-FM is much quicker than negative fake planet injection methods (Lagrange et al. 2010; Marois et al. 2010), as the stellar PSF subtraction algorithm, KLIP, needs to be run only once. In this paper, we demonstrate the advantages of KLIP-FM for precise astrometry and constraining planetary orbits by applying it to GPI observations of β Pic b reduced using KLIP and ADI+SDI. In Section 2, we describe our new astrometry technique, Bayesian KLIP-FM Astrometry (BKA), in which we forward model the PSF of the planet with KLIP-FM and then use the forward model in a Bayesian framework to measure the position of the planet while also modeling the correlated nature of the noise. In Section 3, we validate our technique by fitting an orbit to the data and showing that our astrometry and
uncertainties are consistent with Keplerian motion with no obvious systematics. Finally, in Section 4, we apply our new astrometry to constrain the orbit of β Pic b and place the tightest constraints yet on the transit of the planet and its Hill sphere.

## 2. OBSERVATIONS AND DATA REDUCTION

### 2.1. Observations

To obtain a large temporal baseline of GPI astrometric points, we compiled GPI data of β Pic from commissioning (Gemini program GS-ENG-GPI-COM), an astrometric monitoring campaign of β Pic b (Gemini programs GS-2015A-Q-21 and GS-2015B-Q-9; PI: Graham), the GPI Exoplanet Survey (GS-2014B-Q-500; PI: Macintosh), and a Gemini Large and Long Program to characterize debris disks (GS-2015B-LP-6; PI: Chen). All the data used in the following analysis are listed in Table 1. Most of these data were published and analyzed in a previous study characterizing the orbit of β Pic b by Millar-Blanchaer et al. (2015). We have reprocessed those data with the new astrometry methods presented in this paper and combined them with five additional new epochs.

There were three data sets in Millar-Blanchaer et al. (2015) that we did not use. We did not use the polarimetry data set as there was no instrumental PSF obtained with the data to forward model. The data set on 2014 March 23 was taken during tests of the adaptive optics system, causing the instrumental PSF to be highly varying and making it unsuitable for forward modeling. The 2015 January 24 data set contained artifacts in the construction of the forward model that we could not remove. We chose not to include the measured astrometry from Millar-Blanchaer et al. (2015) for the data sets we omitted in order to maintain homogeneity in astrometric data sets and reduce potential systematic errors. We also note that in the last data set taken on 2016 January 21, there was saturation on the edge of the occulting mask due to bad seeing. This affected the forward modeling of β Pic b, which was also near the occulting mask so a large portion of the frames taken could not be used in this analysis.

### 2.2. Reducing Raw GPI Data

The raw integral field spectrograph (IFS) data were reduced to form three-dimensional (3D) spectral data cubes using the GPI Data Reduction Pipeline (DRP) version 1.2.1 or 1.3 (Perrin et al. 2014). There were no significant changes between the two versions of the GPI DRP that impacted the astrometry of β Pic b. We used the same data reduction process as in Millar-Blanchaer et al. (2015) and will summarize them here. First, dark subtraction and bad pixel correction were applied to each two-dimensional (2D) image. For earlier data sets in which cryocooler vibration induced correlated noise on the detector, the frames were “destriped” to remove this noise (Ingraham et al. 2014). Then, we corrected for instrument flexure using an argon arc lamp taken before each sequence to align each individual spectrum for extraction (Wolff et al. 2014). Then each 2D frame was turned into a spectral data cube, corrected for any remaining bad pixels, and corrected for distortion (Konopacky et al. 2014). For the K-band data, thermal background frames were taken along the sequence. We constructed thermal background cubes in the same fashion and subtracted them from the K-band data.

To spatially register and photometrically calibrate our data, we used the GPI DRP to measure the flux and location of the “satellite” spots, which are centered on the occulted star and imprinted with its attenuated spectrum (Marois et al. 2006b; Sivaramakrishnan & Oppenheimer 2006; Wang et al. 2014). The satellite spot fluxes were used to derive a flux calibration in each spectral channel, which we used when constructing the forward models in Section 2.4.

The location of the occulted star at each wavelength in each data cube was found using a least squares fit to all of the satellite spots' positions and the magnitude of the atmospheric differential refraction. The occulted star center is used to align all the images together before PSF subtraction, and is crucial for determining the astrometry of β Pic b relative to its host star. The precision on the star center is 0.05 pixels (0.7 mas) for satellite spots with $S/N > 20$ (Wang et al. 2014), which is certainly the case for all of our data on bright stars like β Pic.

### 2.3. PSF Subtraction

To subtract the stellar PSF from each individual data cube, we used *pyKLIP* (Wang et al. 2015), a Python implementation of the KLIP algorithm. We used both ADI and SDI to decorrelate the stellar PSF from the PSF of β Pic b. As we were only concerned with β Pic b, we applied our PSF subtraction only on the annulus that included the planet, rather than the full image.

We adjusted three main parameters for the PSF subtraction, depending on the data set. The first was the number of modes used from the KL transform to model the stellar PSF. The second was an exclusion criteria for picking reference PSFs. The exclusion criteria is similar to the quantity $N_δ$ in Lafrenière et al. (2007) and is defined by the number of pixels that β Pic b would move azimuthally and radially in an observing sequence due to ADI and SDI. Third, we toggled an ~11 pixel wide spatial high-pass filter that was applied to some data sets before PSF subtraction. The high-pass filter was implemented using a Gaussian filter in Fourier space and the ~11 pixel size was chosen to remove low-frequency background without significantly distorting any point sources in the image. In general, increasing the number of KL modes, decreasing the exclusion criteria, and applying a high-pass filter improves the subtraction of the stellar PSF. However, taken to the extreme, all three options attenuate signal from the planet. Additionally, for a planet as bright as β Pic b, the forward modeling described in Section 2.4 may not be valid when the PSF subtraction becomes too aggressive (Pueyo 2016). Thus, to optimize the signal of the planet while maintaining the validity of the forward modeling, we varied these parameters for each data set. As β Pic b is bright, in most of the data sets we used an exclusion criteria of four pixels, which is slightly greater than 1 $λ/D$, where $λ$ is the wavelength and $D$ is the diameter of the telescope. When there was little field rotation for ADI, we decreased the exclusion criteria to not overly restrict our reference PSFs, but also decreased the number of KL modes used to avoid being too aggressive. When the speckle noise was too bright due to observing at a shorter wavelength or when the planet moved closer to the star, we applied a high-pass filter to remove some of the diffracted starlight. We list the chosen parameters for each data set in Table 1.

After PSF subtraction, all images were rotated so that north is up and east is left. Then, all the frames were mean combined
<table><caption>**Table 1**
Observations and Data Reduction Parameters for GPI Data on β Pic</caption><thead><tr><th>UT Date</th><th>Program</th><th>Filter</th><th>Exposure Time (s)</th><th>Field Rotation (°)</th><th>Average Seeing<sup>a</sup> (′′)</th><th>KL Modes</th><th>Exclusion Criteria (pixels)</th><th>High-pass Filtered</th><th>Fitting Box Size (pixels)</th></tr></thead><tbody><tr><td>2013 Nov 16</td><td>GS-ENG-GPI-COM</td><td>K1</td><td>1789</td><td>26</td><td>1.09</td><td>7</td><td>4</td><td>No</td><td>13</td></tr><tr><td>2013 Nov 16</td><td>GS-ENG-GPI-COM</td><td>K2</td><td>1253</td><td>18</td><td>0.93</td><td>7</td><td>4</td><td>No</td><td>13</td></tr><tr><td>2013 Nov 18</td><td>GS-ENG-GPI-COM</td><td>H</td><td>2446</td><td>32</td><td>0.68</td><td>7</td><td>4</td><td>No</td><td>13</td></tr><tr><td>2013 Dec 10</td><td>GS-ENG-GPI-COM</td><td>H</td><td>1312</td><td>39</td><td>0.77</td><td>7</td><td>4</td><td>No</td><td>13</td></tr><tr><td>2013 Dec 10</td><td>GS-ENG-GPI-COM</td><td>J</td><td>1597</td><td>19</td><td>0.70</td><td>7</td><td>4</td><td>No</td><td>13</td></tr><tr><td>2013 Dec 11</td><td>GS-ENG-GPI-COM</td><td>H</td><td>410</td><td>65</td><td>0.46</td><td>7</td><td>4</td><td>No</td><td>11</td></tr><tr><td>2014 Nov 8</td><td>GS-ENG-GPI-COM</td><td>H</td><td>2147</td><td>25</td><td>0.77</td><td>7</td><td>4</td><td>No</td><td>11</td></tr><tr><td>2015 Apr 2</td><td>GS-2015A-Q-21</td><td>H</td><td>1312</td><td>10</td><td>0.51</td><td>2</td><td>2</td><td>No</td><td>11</td></tr><tr><td>2015 Nov 6</td><td>GS-2014B-Q-500</td><td>H</td><td>2207</td><td>28</td><td>...</td><td>7</td><td>4</td><td>No</td><td>11</td></tr><tr><td>2015 Dec 5</td><td>GS-2015B-LP-6</td><td>J</td><td>4948</td><td>66</td><td>0.92</td><td>7</td><td>4</td><td>Yes</td><td>11</td></tr><tr><td>2015 Dec 22</td><td>GS-2015B-Q-9</td><td>H</td><td>2088</td><td>19</td><td>0.76</td><td>7</td><td>4</td><td>No</td><td>11</td></tr><tr><td>2016 Jan 21</td><td>GS-2015B-Q-9</td><td>H</td><td>954</td><td>17</td><td>1.18</td><td>10</td><td>2</td><td>Yes</td><td>7</td></tr></tbody></table>

**Notes.**

<footnote><sup>a</sup> Measured by the Gemini DIMM.</footnote>

<footnote><sup>b</sup> Seeing monitor data were not available for this observation.</footnote>
in both the spectral and temporal directions, resulting in a single PSF-subtracted frame for each data set.

## 2.4. Constructing the Forward Model

After stellar PSF subtraction, the PSF of the planet is distorted by over-subtraction, caused by the presence of the planet in the data we are subtracting, and self-subtraction, caused by the presence of the planet in the reference images. Over- and self-subtraction perturb the astrometric and photometric properties of the planet's PSF and prevents a straightforward measurement of its position. Typically, the biases and uncertainties in the astrometry caused by stellar PSF subtraction are estimated by injecting fake planets into the data at other PAs and comparing the retrieved position to the injected position. However, the over-subtraction and self-subtraction that distort the planet PSF are deterministically caused by the existence of a planet in the data and its apparent motion in the reference images induced by ADI and SDI. In turn, these features, if modeled, can inform us about the location of the planet and improve our astrometric precision and accuracy.

Recently, Pueyo (2016) introduced KLIP-FM, an analytic framework to compute the effect of a planet on the KL modes and use these perturbations to reconstruct the over-subtraction and self-subtraction features. Using this technique, we are able to generate the PSF of the planet after PSF subtraction (see Appendix A for a detailed description of the procedure). Briefly, we use the PSF of the planet before PSF subtraction, the apparent motion of the planet due to ADI and SDI, and a model of the planet spectrum to compute the distorted PSF of the planet after PSF subtraction. The result is a 2D broadband planet PSF, which we call $F$, centered at $(x_0, y_0)$, our initial estimate for the location of the planet.

For our GPI data, the forward models were generated using the implementation of KLIP-FM in `pyKLIP`. For each data set, we used the same parameters as our PSF subtraction to construct the forward model. To construct a model of the instrumental PSF at each wavelength, we used the average PSF of the satellite spots across all images for that wavelength. We used a basic center of light centroiding routine to measuring an approximate $x$ and $y$ positions of the planet $(x_0$ and $y_0)$ in the PSF-subtracted image. This initial estimate for the position is good to within a pixel for our GPI data. For the input planet spectrum of the forward model, we used the normalized best-fit model spectrum from Chilcote et al. (2015), which has an effective temperature of 1650 K and a $log(g) = 4.0$ (cgs units). We scaled the spectrum by eye to approximately match the contrast of $\beta$ Pic b in our data and used the satellite spot fluxes to convert from contrast to digital numbers. We found that varying the spectrum and photometry had negligible impacts on our measured astrometry. Even a spectrum differing by ~25% in shape changed the astrometry by <0.2 mas, significantly smaller than the uncertainties we find in the following analysis.

With this input information, we used `pyKLIP` to generate one forward-modeled PSF, $F$, for each data set. In Figure 1, we show an example forward-modeled PSF (not optimized to fit the data) and comparison to data for our 2014 November 18 data set on $\beta$ Pic b. Qualitatively, the forward model matches all the features seen in data, including the prominent negative self-subtraction lobes on either side of the planet. As `pyKLIP` parallelizes the computation, the generation of the forward

<figure><img src="image_3.png" /><figcaption>**Figure 1.** (Left) An image of $\beta$ Pic b from the 2013 November 18 GPI H-band data after stellar PSF subtraction. (Right) Unoptimized forward model of $\beta$ Pic b for the same data set. The forward model has not yet been optimized to fit the data using the MCMC procedure discussed in Section 2.5, but should already be accurate to within a pixel. To see the best-fit forward models and residuals for each GPI data set, see Appendix B. The star is to the upper left of the planet and far outside the region shown here.</figcaption></figure>

models is quick. On a 32 core machine with AMD Opteron 6378 processors clocked at 2.3 GHz, forward models for all 37 channels of a representative 37 cube GPI data set were generated in 15 minutes: 4 minutes of overheads for preprocessing and generating the instrumental PSFs and 11 minutes to execute KLIP-FM and create the forward model. We note that we chose to run the forward model on a large annulus to examine noise properties in the data and that the computation time for KLIP-FM decreases ~20%–30% if a small sector around the planet was used instead.

## 2.5. Locating the Planet with Bayesian Parameter Estimation

To use the forward-modeled PSF, $F$, to perform astrometry, we developed a Bayesian framework to fit $F$ to the data, account for correlated noise, and estimate our fitting uncertainties. First, ignoring the correlated noise, we can use three parameters to fit the forward model to the data: the location on the planet in $x$ ($x_p$), the location of the planet in $y$ ($y_p$), and a scale factor ($\alpha$) to match the flux of $F$ to the data. We can then write the posterior probability for $x_p$, $y_p$, and $\alpha$ given the data $D$ using Bayes' Theorem as

$$P(x_p, y_p, \alpha|D) = P(D|x_p, y_p, \alpha)P(x_p, y_p, \alpha), \quad (1)$$

ignoring normalization constants. The first term on the right-hand side is the likelihood and the second is the prior.

To construct the likelihood, we must first use our input parameters and $F$ to generate a model to compare to the data. We scale $F$ by $\alpha$ and recenter it from its guessed location ($x_0$, $y_0$) to ($x_p$, $y_p$). The residual, $R$, between the model and the data, $D$, in fitting region $\mathcal{F}$ is calculated by

$$R \equiv R(x_p, y_p, \alpha) = (D - \alpha F(x_p, y_p))_{\mathcal{F}}, \quad (2)$$

where fitting region $\mathcal{F}$ is a fixed rectangular box centered at the approximate location of the planet in the data. We pick the size of the fitting region to be a few $\lambda/D$ to encompass the PSF and the self-subtraction wings. We varied the fitting box size for each of our data sets in order to keep the fit focused on the area where signal from the planet can be seen. We list the size of the fitting box for each data set in Table 1. Note that we do not fit any background term to the image, as one of the first steps of KLIP is to subtract off the mean of the image, removing any spatially constant background.
For data with uncorrelated errors, the log of the likelihood function of the data for a particular model is

$$ P(D|x_p, y_p, \alpha) = \ln \mathcal{L}(x_p, y_p, \alpha) \\ = -\frac{1}{2} \sum_{i}^{N_{\text{pix}}} \left[ \frac{R_i^2(x_p, y_p, \alpha)}{\sigma_i^2} + \ln(2\pi\sigma_i) \right], \quad (3) $$

where $\sigma_i$ is the uncertainty in pixel $i$ and $N_{\text{pix}}$ is the number of pixels in the fitting region. However, the assumption of uncorrelated errors does not hold for images limited by speckle noise. Except in the cases of very aggressive PSF subtraction or at much greater separations from the star, the residual noise after PSF subtraction in GPI data of a bright star like $\beta$ Pic is still dominated by correlated speckle noise, which has a correlation scale that depends on the aggressiveness of the PSF subtraction. Due to the bright nature of $\beta$ Pic b, we could not use a very aggressive reduction to ensure that the planet remained a perturbation on our KL modes, thereby preserving the validity of the analytical forward modeling technique. The conservative PSF subtraction combined with the close separation of $\beta$ Pic b from its bright host star required us to capture the correlated nature of the noise. We thus write our likelihood function instead as

$$ \ln \mathcal{L} = -\frac{1}{2}(R^T C^{-1} R + \ln(\det C) + N_{\text{pix}} \ln(2\pi)), \quad (4) $$

where $C$ is the covariance matrix of size $N_{\text{pix}} \times N_{\text{pix}}$ and $R$ is a $N_{\text{pix}} \times 1$ matrix.

We applied a Gaussian process framework to characterize the covariance in the noise (see Czekala et al. 2015 for an in-depth explanation of the application of Gaussian processes to astronomical data). We only aim to model the correlations within each individual speckle, which spans $\lambda/D$ in spatial extent. While in a single unprocessed frame speckle noise also has additional correlations on much larger scales, there are no significant correlations between speckles at larger spatial scales in our small fitting region due to stellar PSF subtraction with KLIP and averaging uncorrelated speckles together when collapsing the frames in our ADI+SDI sequence. Thus, within our fitting region, the dominant correlation in our noise is from pixels within a single speckle, and thus is the one correlation we modeled.

Following the procedure in Czekala et al. (2015) for fitting one-dimensional (1D) correlations in stellar spectra, we used the Matérn covariance function parametrized with $\nu = 3/2$ to model the correlated speckle noise. We chose the Matérn function with $\nu = 3/2$ as it better fits the correlations at larger separations compared to a simple squared exponential relation. However, we note that switching between covariance functions that have similar shapes does not significantly alter the error bars. We also chose to assume symmetric noise as we did not find any difference in the correlation scale of our noise in the radial and azimuthal directions. This is likely due to the fact that we used both ADI and SDI to model the stellar PSF and thus are better at subtracting speckle noise, which, before PSF subtraction, is more correlated radially than azimuthally due to the finite spectral bandwidth of the data. In instruments that are not able to utilize SDI, it might be necessary to separate the noise into radial and azimuthal components, each with its own characteristic correlation length. Thus, for our purposes, we

chose the symmetric Matérn covariance function with $\nu = 3/2$, which defines the covariance between two pixels $i$ and $j$ separated by a distance of $r_{ij} = \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2}$ as

$$ C_{ij} = \sigma_i \sigma_j \left( 1 + \frac{\sqrt{3} r_{ij}}{l} \right) \exp \left( -\frac{\sqrt{3} r_{ij}}{l} \right), \quad (5) $$

where $\sigma_i$ and $\sigma_j$ are the uncertainties for each pixel and $l$ is a characteristic correlation length scale that increases when noise is correlated at larger spatial scales. We calculated the uncertainty for each pixel by computing the standard deviation of pixels in an annulus with a width of $\sim 2\lambda/D$ (6 pixels), centered on the star, and with a mean radius equal to the distance between that pixel and the star. Any pixels containing signal from the planet were masked and not used in estimating the noise in the annulus.

The correlation length $l$ is not fixed, but rather kept as a hyperparameter parameter in our Bayesian framework that we will marginalize over in the end. Thus, the final posterior we are trying to calculate is really

$$ P(x_p, y_p, \alpha, l|D) = P(D|x_p, y_p, \alpha, l)P(x_p, y_p, \alpha, l), \quad (6) $$

where $\alpha$ and $l$ are both hyperparameters as we are only interested in the astrometry.

Compared to our likelihood function, our prior is relatively straightforward. We allow for uniform priors in $x$ and $y$ within 2 pixels from our initial guess location. Similarly, we allow for an uniform prior between 0 and 5 for $\alpha$ to determine how much to scale $F$, which was already scaled to an approximate contrast of $\beta$ Pic b. The correlation length $l$ has a uniform prior between 0 and 10 pixels, which provides an ample range to explore the correlations within individual speckles of size $\lambda/D \approx 3$ pixels.

We used the Goodman & Weare (2010) Affine Invariant MCMC sampler implemented in the `emcee` Python package (Foreman-Mackey et al. 2013) to sample the posterior distribution and custom `cython` code to quickly generate the covariance matrix as we vary $l$. The MCMC sampler was run for 800 steps using 100 walkers, with a “burn-in” of 200 steps beforehand. In Figure 2, we show the posterior distributions from the MCMC fit for the 2014 November 18 H-band data set as a representative posterior distribution. The value of $l$ is close to $\lambda/D \approx 3$ pixels, indicating that we are accurately fitting the correlated speckle noise. We use the resulting posterior distributions to calculate the most likely values and uncertainties for the location of $\beta$ Pic b in our image at each epoch.

To convert our results to more useful physical values, we need to convert our measured location of $\beta$ Pic b to its positional offset from its host star in right ascension (R.A.) and declination (decl.). As our images are already rotated so that $-x$ is positive R.A. and $+y$ is positive decl., it is straightforward to convert from image fates to sky coordinates. We use the satellite spots to measure the location of the star behind the coronagraph, which has a precision of 0.7 mas (Wang et al. 2014). This allows us to derive the separation of the planet from the star in pixels. To convert from pixel separations to physical separations in R.A. and decl., we use the most recent astrometric calibration numbers from De Rosa et al. (2015): a plate scale of 14.166 ± 0.007 mas lenslet⁻¹ and a residual north angle offset of 0°10 ± 0°13 from the north angle value used in the GPI DRP (versions 1.2.1 to 1.3, the current version). These astrometric calibration numbers show no
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Posterior distribution of the four parameters in the MCMC fit for the astrometry for the 2014 November 18 epoch. The vertical dashed lines in the marginalized posterior distribution plots indicate the 16th, 50th, and 84th percentile values.</figcaption></figure>

significant change over time, so we apply them to all our epochs of data. Then, we assume all of these error terms are uncorrelated and add them in quadrature with our measurement errors from our MCMC fit to determine our full astrometric precision.

The combination of the forward-modeled PSF and the Bayesian framework makes up the Bayesian KLIP-FM Astrometry technique we introduce in this work. We apply BKA for all 12 GPI data sets and report the measured astrometry and error budgets in Table 2. The best-fit forward models and residuals to the fit are shown in Appendix B. On most of our data sets, we are not limited just by the uncertainty in the location of the planet, which was as low as 0.3 mas. The uncertainties in the location of the star and north angle also make significant contributions to the error budget. Typically, we achieved ~1 mas precision on the relative astrometry between β Pic b and its host star. This is a factor of ~2–4 improvement over previous techniques (Millar-Blanchaer et al. 2015) using the same data, indicating that this technique can be useful in reanalysis of archival data to obtain better astrometry in cases where the limiting factor is the uncertainty on the planet position. In two of the later data sets where the planet is observed closer in, we were limited by the S/N of the planet and unable to achieve 1 mas precision. In the 2015 December 5 data set, the noise was higher due to the planet being fainter relative to the star in J-band. In the 2016 January 21 data set, a combination of poor seeing and a small amount of usable data limited our astrometric precision.

Overall, though, this GPI β Pic b data is an excellent demonstration for Bayesian KLIP-FM Astrometry as the planet is bright enough that the extended PSF features, such as the negative self-subtraction lobes, are clearly seen and provide significant information to constrain the position of the planet. For fainter planets, the extended features are harder to distinguish from the noise. As one of the main advantages of BKA over techniques that do not forward model the PSF is being able to forward model the extended self-subtraction lobes, the astrometric improvement would not be as large for lower S/N planets. There still should be some improvement, though, due to accurately modeling the over-subtraction on the core of the PSF and small contributions from the extended features even if they are hard to distinguish from noise. Regardless, in addition to the improved precision, BKA should also more accurately estimate the uncertainties as it fits for the correlation scale of the noise at the location of the planet.

### 3. VALIDATION THROUGH ORBIT FITTING

To explore the validity of our new astrometry technique, we fit a Keplerian orbit to our 12 epochs of astrometry. Since each epoch of astrometry is fit independently, and since the Keplerian orbit is agnostic toward the exact data-reduction methods, having all 12 epochs of astrometry fit the Keplerian orbit would only be possible if all of the astrometry is accurate and precise. If there are errors in estimating the magnitude of the uncertainties or any remaining biases in our measurements, we expect them to become evident in the residuals of the Keplerian fit either as systematic trends or fit outliers.
<table><caption>**Table 2**

Astrometric Error Budget and Measured Astrometry of β Pic b</caption><thead><tr><th>Data Set</th><th>Planet x/y Uncertainty (mas)</th><th>Star x/y Uncertainty (mas)</th><th>Plate Scale Uncertainty (mas)</th><th>PA Uncertainty (°)</th><th>ΔR.A. (mas)</th><th>ΔDecl. (mas)</th><th>Radial Separation (mas)</th><th>Position Angle (°)</th></tr></thead><tbody><tr><td>2013 Nov<br>16 K1</td><td>0.6/0.7</td><td>0.7/0.7</td><td>0.3</td><td>0.13</td><td>-228.5 ± 1.3</td><td>-366.2 ± 1.1</td><td>431.6 ± 1.0</td><td>212.0 ± 0.2</td></tr><tr><td>2013 Nov<br>16 K2</td><td>0.5/0.4</td><td>0.7/0.7</td><td>0.3</td><td>0.13</td><td>-229.2 ± 1.2</td><td>-364.5 ± 1.0</td><td>430.6 ± 0.9</td><td>212.2 ± 0.2</td></tr><tr><td>2013 Nov 18 H</td><td>0.3/0.3</td><td>0.7/0.7</td><td>0.3</td><td>0.13</td><td>-229.1 ± 1.1</td><td>-364.7 ± 1.0</td><td>430.6 ± 0.8</td><td>212.1 ± 0.2</td></tr><tr><td>2013 Dec 10 H</td><td>0.4/0.4</td><td>0.7/0.7</td><td>0.3</td><td>0.13</td><td>-227.9 ± 1.2</td><td>-359.9 ± 1.0</td><td>426.0 ± 0.9</td><td>212.3 ± 0.2</td></tr><tr><td>2013 Dec 10 J</td><td>0.6/0.7</td><td>0.7/0.7</td><td>0.3</td><td>0.13</td><td>-227.2 ± 1.3</td><td>-361.1 ± 1.2</td><td>426.6 ± 1.1</td><td>212.2 ± 0.2</td></tr><tr><td>2013 Dec 11 H</td><td>0.5/0.4</td><td>0.7/0.7</td><td>0.3</td><td>0.13</td><td>-227.8 ± 1.2</td><td>-359.2 ± 1.0</td><td>425.4 ± 0.9</td><td>212.4 ± 0.2</td></tr><tr><td>2014 Nov 8 H</td><td>0.5/0.5</td><td>0.7/0.7</td><td>0.2</td><td>0.13</td><td>-194.0 ± 1.1</td><td>-299.1 ± 1.0</td><td>356.5 ± 0.9</td><td>213.0 ± 0.2</td></tr><tr><td>2015 Apr 2 H</td><td>0.4/0.5</td><td>0.7/0.7</td><td>0.2</td><td>0.13</td><td>-172.1 ± 1.0</td><td>-266.5 ± 1.0</td><td>317.2 ± 0.9</td><td>212.9 ± 0.2</td></tr><tr><td>2015 Nov 6 H</td><td>0.7/0.7</td><td>0.7/0.7</td><td>0.2</td><td>0.13</td><td>-137.8 ± 1.1</td><td>-207.2 ± 1.0</td><td>248.8 ± 1.0</td><td>213.6 ± 0.3</td></tr><tr><td>2015 Dec 5 J</td><td>1.2/1.3</td><td>0.7/0.7</td><td>0.2</td><td>0.13</td><td>-133.9 ± 1.5</td><td>-200.5 ± 1.5</td><td>241.1 ± 1.4</td><td>213.7 ± 0.4</td></tr><tr><td>2015 Dec 22 H</td><td>0.5/0.5</td><td>0.7/0.7</td><td>0.2</td><td>0.13</td><td>-130.0 ± 1.0</td><td>-194.7 ± 0.9</td><td>234.1 ± 0.9</td><td>213.7 ± 0.2</td></tr><tr><td>2016 Jan 21 H</td><td>1.8/1.6</td><td>0.7/0.7</td><td>0.2</td><td>0.13</td><td>-126.8 ± 2.0</td><td>-185.8 ± 1.8</td><td>225.0 ± 1.8</td><td>214.3 ± 0.5</td></tr></tbody></table>
<table><caption>**Table 3**
Orbital Parameters of β Pic b</caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">Unit</th><th rowspan="2">Prior Range</th><th rowspan="2">Prior Distribution</th><th colspan="3">Posterior Percentiles</th></tr><tr><th>16</th><th>50</th><th>84</th></tr></thead><tbody><tr><td>Semimajor axis (a)</td><td>au</td><td>4–40</td><td>Uniform in log a</td><td>9.02</td><td>9.66</td><td>10.78</td></tr><tr><td>Epoch of Periastron (τ)</td><td>...</td><td>-1.0–1.0</td><td>Uniform in τ</td><td>0.32</td><td>0.73</td><td>0.87</td></tr><tr><td>Argument of Periastron (ω)</td><td>°</td><td>-360–360</td><td>Uniform in ω</td><td>192.8</td><td>205.8</td><td>258.4</td></tr><tr><td>Position Angle of the Ascending Node (Ω)</td><td>°</td><td>25–85</td><td>Uniform in Ω</td><td>31.67</td><td>31.76</td><td>31.84</td></tr><tr><td>Inclination (i)</td><td>°</td><td>81–99</td><td>Uniform in cos i</td><td>88.70</td><td>88.81</td><td>88.93</td></tr><tr><td>Eccentricity (e)</td><td>...</td><td>0.00001–0.99</td><td>Uniform in e</td><td>0.027</td><td>0.080</td><td>0.171</td></tr><tr><td>Total Mass (M<sub>T</sub>)</td><td>M<sub>☉</sub></td><td>0–3</td><td>Uniform in M<sub>☉</sub></td><td>1.76</td><td>1.80</td><td>1.83</td></tr><tr><td colspan="7" style="text-align: left; font-weight: bold;">Derived Parameters</td></tr><tr><td>Period (P)</td><td>years</td><td>...</td><td>...</td><td>20.21</td><td>22.47</td><td>26.24</td></tr><tr><td>Hill Sphere Ingress</td><td>MJD</td><td>...</td><td>...</td><td>57,840</td><td>57,846</td><td>57,854</td></tr><tr><td>1/2 Hill Sphere Ingress</td><td>MJD</td><td>...</td><td>...</td><td>57,916</td><td>57,924</td><td>57,934</td></tr><tr><td>Closest Approach Date</td><td>MJD</td><td>...</td><td>...</td><td>57,986</td><td>57,996</td><td>58,008</td></tr><tr><td>1/2 Hill Sphere Egress</td><td>MJD</td><td>...</td><td>...</td><td>58,056</td><td>58,069</td><td>58,082</td></tr><tr><td>Hill Sphere Egress</td><td>MJD</td><td>...</td><td>...</td><td>58,132</td><td>58,147</td><td>58,162</td></tr></tbody></table>

We followed the same analysis as in Millar-Blanchaer et al. (2015) to obtain the orbital elements of β Pic b using a MCMC fit with the parallel-tempered sampler implemented in emcee. We combined the 12 GPI astrometric points presented here along with the data sets presented in Chauvin et al. (2012) and Nielsen et al. (2014). Unlike Millar-Blanchaer et al. (2015), we did not explicitly include the radial velocity measurement of β Pic b from Snellen et al. (2014) in order to limit potential systematics between instruments, but we do use it to constrain the prior on the PA of the ascending node and thus the direction of the orbit (i.e., we know that that β Pic b has been moving toward us since ~2007). Our model of the orbit fits the same seven parameters as Millar-Blanchaer et al. (2015). For convenience, the parameters are listed in Table 3.

The MCMC sampler was run for 30,000 steps using 1024 walkers at each of the 20 temperatures after 30,000 steps of “burn-in” to allow the walkers to converge. We thinned the chains to remove any remaining correlations, keeping every 75th step to result in effectively 400 samples per walker for a total of 409,600 samples to construct our posterior.

The posteriors on the seven parameters in our model are shown in Figure 3 and the 16%, 50%, and 84% percentiles of the marginalized posterior distribution for each parameter is listed in Table 3. For comparison, we plot 500 randomly chosen possible orbits along with the measured astrometry in Figure 4. The fit residuals of our measured GPI points are shown in Figure 5. The residuals in R.A. and decl. offset (which are the parameters we use in our MCMC) are consistent with zero and do not show any systematics. The residuals in radial separation and PA, neither of which are used in our fit, are slightly further away from zero, but do not indicate any obvious errors in either our astrometry or our error estimates.

To quantitatively assess the validity of our measurements, we use the reduced chi-squared ($\chi^2_{red}$) statistic, which measures how consistent our astrometry is with a Keplerian orbit. Ideally, $\chi^2_{red}$ would be unity if all measurements and uncertainties were accurate. However, biases and improperly estimated errors would cause it to deviate from unity. Due to systematic astrometric calibration offsets between instruments that have not been characterized, we expect $\chi^2_{red}$ to be slightly above unity. For example, Millar-Blanchaer et al. (2015) reported a $\chi^2_{red}$ of 1.55$^{+0.09}_{-0.05}$ when combining GPI astrometry, measured

using standard techniques in the field, with previous astrometric measurements

Even with ~2–4 times smaller error bars on the GPI points, the $\chi^2_{red}$ of our accepted orbits is $1.53^{+0.08}_{-0.06}$, unchanged from the $1.55^{+0.09}_{-0.05}$ reported in Millar-Blanchaer et al. (2015). With BKA contributing 12 out of the 30 astrometry measurements, if the BKA technique contained biases larger than one milliarcsecond, they would have caused a significant increase in $\chi^2_{red}$. Likewise, if we had been too optimistic with our error estimates, $\chi^2_{red}$ should have also increased as the reported uncertainties would not have matched the scatter in the measurements. The lack of increase in $\chi^2_{red}$ indicates that the more precise astrometry from BKA is not biased and has accurate uncertainties.

In addition to the Keplerian orbit fit, we can examine the accuracy of our estimated one milliarcsecond uncertainties by looking at the astrometry measured on the same or consecutive days. As we do not expect the planet's position to change significantly in the span of a single day, the three measurements in 2013 November and the three measurements in 2013 December ought to be consistent with each other. Indeed, our measurements in Table 2 indicate that in both sets of astrometry, the measurements agree at the milliarcsecond level, confirming that our estimated uncertainties are accurate

Thus, the well behaved residuals of our GPI measurements and the lack of change in $\chi^2_{red}$ from Millar-Blanchaer et al. (2015) even with significantly smaller error bars lead us to conclude the measured astrometry using BKA are accurate and free from biases. The fact the residuals are consistent with zero, the lack of change in $\chi^2_{red}$, and the consistency of repeated measurements taken around the same time all indicate that the one milliarcsecond uncertainties estimated from BKA are also accurate. Together, these multiple assessments of the validity of BKA all indicate that this technique is both accurate and precise.

# 4. DISCUSSION

## 4.1. The Orbit of β Pic b

Having demonstrated the accuracy and precision of this new technique, we now analyze the new constraints on the orbit of β Pic b. Not surprisingly, the estimates for a and e have not
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Posterior distributions for the seven orbital elements in our Keplerian orbit model along with inferred distributions of possible dates for transit events in the top right corner. The vertical dashed lines in the marginalized posterior distribution plots indicate the 16th, 50th, and 84th percentile values. For the transit events, the dotted line corresponds to the ingress and egress of the full Hill sphere, the dashed line corresponds to the ingress and egress of the half Hill sphere, and the solid line corresponds to the date of closest approach.</figcaption></figure>

changed significantly from Millar-Blanchaer et al. (2015) since the GPI points reanalyzed in this paper do not have a sufficiently long time baseline to constrain these parameters. As seen in Figure 4, all but one astrometric measurement is on one-half of the orbit curve. The other half of the orbit is not as well constrained, leaving a degeneracy in *a* and *e*. This degeneracy can be broken with more measurements obtained when the planet appears on the other side of the star. Better constraints on *a* and *e* will provide better insight on how β Pic b interacts with the debris disk and potential unseen

planets in the system. The new total mass of the system, $M_T = 1.80_{-0.04}^{+0.03} M_\odot$, is significantly higher than the $M_T = 1.61 \pm 0.05 M_\odot$ from Millar-Blanchaer et al. (2015). This new total system mass, which effectively measures the mass of the star at this precision, is in better agreement with the stellar mass of $1.75 M_\odot$ derived from stellar photometry (Crifo et al. 1997). Our measurement of the PA of the ascending node, $\Omega$, slightly improves upon the value obtained by Millar-Blanchaer et al. (2015), but the overall value remains consistent. Thus, β Pic b is still consistent with being the
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (Top) Offset of β Pic b in R.A. (blue) and decl. (red) with respect to β Pic as a function of time. We have plotted the measured data and 500 randomly chosen accepted orbits from the MCMC sampler. (Middle) Radial separation of β Pic b from the star as a function of time. The same 500 randomly chosen orbits have are also plotted (Bottom) PA as a function of time for the data and the 500 randomly chosen orbits. To keep the data compact, we have wrapped PA by 180° to only consider PAs between 180° and 360°. This allows for easy comparison of the 2003 point, which is nominally at a PA of 34.24° but here displayed at a PA of 214.4°. The dashed black line indicates a constant PA of 212.4°, the weighted mean of all points. If the planet were to transit, we would not be able to see a significant deviation from constant PA in time. For all the plots, error bars are also plotted but many are too small to be seen on this scale.</figcaption></figure>

planet responsible for the known warp in the debris disk (e.g., Dawson et al. 2011). Compared to Millar-Blanchaer et al. (2015), the argument of the periastron, $\omega$ in our fit has increased from $156^{\circ+33}_{-76}$ to $206^{\circ+52}_{-13}$. This new value is consistent with $\omega = 200^{\circ} \pm 20^{\circ}$ that is required for the falling evaporating bodies (FEBs) scenario proposed by Thébault & Beust (2001) to explain redshifted absorption features in β Pic spectrum. Note that under previous definitions of the orbital parameters, this has been expressed as $\omega = -70^{\circ} \pm 20^{\circ}$ from the line of sight. This scenario also requires a slightly eccentric orbit, which is consistent with our derived orbital parameters. For a more in-depth discussion of β Pic b's relationship to the debris disk and the FEB scenario, we direct the reader to Millar-Blanchaer et al. (2015).

The biggest improvement in our understanding of the orbit of β Pic b is the improved constraint on the inclination of the orbit. We find the inclination to be $i = 88.81^{\circ}_{-0.11}^{+0.12}$, which allows us to place the tightest constraints on the probability that β Pic b will transit its host star. Assuming an angular diameter of the star of 0.736 mas (Defrère et al. 2012) and considering

the range of $a$ from our orbit fit, we find that in most cases we need $|i - 90^{\circ}| < 0.05$ in order for the planet to transit. With our current constraints on the inclination, we have ruled out the possibility that β Pic b will transit at $10\sigma$ significance. This tight constraint on the inclination and transit probability is due to the slightly longer time baseline and the improved precision in the measured PA of the GPI astrometry compared to Millar-Blanchaer et al. (2015). For an edge-on orbit that transits the star, we should see no significant change in PA over time. However, Figure 4 shows that the GPI points alone reveal a significant increase in PA over time. Thus, regardless of systematic astrometric calibration errors between instruments, we conclude that β Pic b will not transit its star.

## 4.2. Hill Sphere Transit

Unlike the planet, β Pic b's Hill sphere, the region around the planet that could contain gravitationally bound circumplanetary material, will transit the star. We define the radius of the Hill sphere as $r_H \approx a(1 - e)^{\frac{1}{2}} \sqrt{\frac{m}{3M}}$ using the approximate
<figure><img src="image_7.png" /></figure>

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

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

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

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

<figure><img src="image_12.png" /><figcaption>**Figure 5.** Residuals to the orbit fit for the average of 500 randomly chosen accepted orbits. The top row shows the residuals in ΔR.A. (blue) and Δdecl. (red) offset, which are the coordinates used in the MCMC analysis to fit the orbit. We also plot radial separation (middle row) and position angle (bottom row). The separation and position angle residuals were not optimized in the MCMC fit. The right column is a zoom-in of the left column showing only the residuals of the GPI astrometry.</figcaption></figure>

form proposed in Hamilton & Burns (1992) where *m* is the mass of the planet and *M* is the mass of the star. Assuming *β* Pic b follows a “hot-start” evolutionary track with a mass of 12.7 ± 0.3 *M*<sub>Jup</sub> (Morzinski et al. 2015) and using the range of semimajor axes, eccentricities, and stellar masses from our MCMC orbit fit, *r*<sub>H</sub> = 1.165<sup>+0.013</sup><sub>-0.016</sub> au (59.9<sup>+0.7</sup><sub>-0.9</sub> mas). Given that our prediction for the closest approach of *β* Pic b will be 9.9<sup>+0.9</sup><sub>-0.8</sub> mas (0.19 ± 0.02 projected au) from the star, the Hill sphere of *β* Pic b will transit the star, as shown in Figure 6.

From our MCMC orbit fit, we can compute when the transit of the Hill sphere will occur. We pick five notable events to focus on: two are the ingress and egress of the Hill sphere which are the extrema in time between which any circumplanetary material could transit; another is the date of closest

approach, which gives the opportunity to probe material closest in to the planet; and the last two are the transit of the sphere that is 1/2 r<sub>H</sub> in radial extent. Almost all stable prograde circumplanetary orbits reside within 1/2 r<sub>H</sub> (Shen & Tremaine 2008), so it is more likely to find material within half a Hill sphere. In Table 3, we list our constraints on the date of these events, and in Figure 3, we plot the posterior distribution for these events. The duration of the Hill sphere transit will be long: ingress is ~2017 April 3 and egress is ~2018 January 29. The 1/2 Hill sphere begins transit ~2017 June 20 and ends transit on ~2017 November 12, with closest approach on ~2017 August 31.

Given the R.A. of β Pic (05h47m17s), the star will not be visible from most ground-based observatories during almost
<figure><img src="image_13.png" /><figcaption>**Figure 6.** Region in the sky that the Hill sphere of β Pic b will sweep across during the closest approach of the planet in 2017. One-hundred randomly chosen accepted orbits (blue) are plotted along with the angular extents of their Hill spheres (light orange) and 1/2 Hill spheres (dark orange). The star, β Pic, is shown in its true angular size (small yellow dot). Both the Hill sphere and 1/2 Hill sphere will pass in front of the star. Note that for clarity we are not plotting the orbital path for β Pic b when it comes back in its orbit and passes behind the star.</figcaption></figure>

the entire time between ingress of the Hill sphere and closest approach. Ground-based telescopes in Antarctica, airborne observatories capable of travelling to Antarctica, and space-based observatories would provide the only opportunities to observe β Pic during this time period. The second half of the Hill sphere transit will be visible from most southern hemisphere ground-based observatories, making β Pic a well suited candidate for photometric monitoring in late 2017.

As it is a rare opportunity to probe circumplanetary material, it is not certain what will be seen when the Hill sphere of β Pic b transits the star. One possibility is that satellites could reside in the Hill sphere. To approximate the photometric transit depths, the 0.736 mas angular diameter of the star corresponds to $1.53 R_\odot$ or $1.065 \times 10^6$ km. Moons as large as Ganymede ($r = 2630$ km) or Io ($r = 1820$ km) would give transit depths of $2.473 \times 10^{-3}$ (2.7 mmag) and $1.711 \times 10^{-3}$ (1.9 mmag), respectively. Detecting these photometric signatures will require a high cadence, as any single satellite orbiting β Pic b will have a transit duration of ~2 days. Additionally, β Pic is a variable star with pulsation timescales of ~0.5 hours and variability amplitudes in B-band of <5 mmags (Koen et al. 2003), so careful modeling of stellar activity is necessary to be sensitive to these transit depths.

Another possibility is that, as β Pic b is still young and evolving, it may harbor a circumplanetary disk or ring system comprised of leftover material from planet formation. Such a hypothesis is not unprecedented as Kenworthy & Mamajek (2015) found evidence for a large circumplanetary disk around an unseen planet in the 1SWASP J140747.93-394542.6 (hereafter J1407) system, which has a similar age (~16 Myr; Mamajek et al. 2012) and thus likely at a similar stage in its evolution. Kenworthy & Mamajek (2015) interpreted the series of complex and deep eclipses in the J1407 light curve as from a circumplanetary disk 0.6 au in radial extent in the process of forming rings due to newly formed satellites. It is plausible that β Pic b can harbor a similar disk as a 0.6 au disk would be 1/2 $r_H$ in extent and consistent with where we would expect stable orbits to reside around β Pic b. Additionally, since β Pic b is young, it is plausible that there is a large amount of circumplanetary material which has yet to be cleared out dynamically. Such a large disk would transit the star and be suitable for detection through photometric monitoring of the host star if the disk is inclined by $\gtrsim 18^\circ$ with respect to the

orbital plane of β Pic b. To maintain such an inclination, the disk needs to be massive enough to prevent stellar tidal forces from aligning the disk to the orbital plane (Zanazzi & Lai 2016).

So far, there have not been many observational constraints of circumplanetary material around β Pic b. Lecavelier Des Etangs et al. (1997) reported a photometric event in 1981 and hypothesized it could be due to the transit of a planet that cleared out a hole in the debris disk around its Hill sphere. However, Millar-Blanchaer et al. (2015) show that the planet is not embedded in the debris disk and we definitely show the planet itself will not transit. Still, our orbit models give a 8% and 4% chance the photometric event coincided with the transit of the Hill sphere and 1/2 Hill sphere respectively, during which time material around the planet could have passed in front of the star. Additionally, a circumplanetary disk or ring may also be detectable through the planet's spectral energy distribution (SED). The planet's near-infrared spectrum would experience extinction if the dust resides between us and the planet, with the magnitude of extinction depending on the amount of dust. In the near-infrared, the extinction would be greater at shorter wavelengths due to the increased scattering and absorption by dust and would produce a spectral slope in the planet's near-infrared SED. Dust around the planet would also scatter starlight, causing the planets SED to appear brighter in the optical, as has been postulated for Fomalhaut b (Kalas et al. 2008). Furthermore, the dust will produce millimeter emission. However, detecting circumplanetary material in the planet's SED will require being able to distinguish it from the circumstellar disk with precise spectral data.

## 5. CONCLUSIONS

In the first part of this work, we have presented a new technique for more precise and accurate astrometry of directly imaged exoplanets using a new analytical forward modeling approach in a robust statistical framework.

1. Using the KLIP-FM framework presented in Pueyo (2016), we are able to analytically forward model the PSF of the planet through the data reduction process, giving us better information on the location of the planet. We apply KLIP-FM to GPI data on β Pic and forward model the PSF of β Pic b using the open source *pyKLIP* package.

2. For a close-in planet orbiting a bright star like in the case of β Pic b, we are limited by correlated speckle noise in our data. We developed a Bayesian framework utilizing Gaussian processes and MCMC to account for the correlated noise and to find the position of the planet simultaneously.

3. With this technique, we have achieved the most precise astrometry on β Pic b to date. On most of our GPI data sets, we achieve ~1 mas precision on the relative separation between β Pic b and its host star, a ~2–4 fold improvement over previous techniques using the same data (Millar-Blanchaer et al. 2015).

4. In data sets where the astrometry is limited by noise and not by astrometric calibration uncertainty, Bayesian KLIP-FM Astrometry approach should improve astrometric precision.

In the second part of this work, we apply our Bayesian KLIP-FM Astrometry technique to the orbit of β Pic b.
1. To validate this new astrometric technique, we used it to measure the position of β Pic b in 12 epochs of GPI data. We combined these 12 astrometric points with 2 previous astrometric monitoring campaigns and fit a Keplerian orbit using MCMC methods. We find the residuals to the fit are consistent with zero and show no apparent systematic trends, indicating that our fit is accurate and the uncertainties we estimate are reliable.

2. Due to the improved PA measurements from our technique, we have the tightest constraints on the inclination of the orbit and can exclude a possible transit of β Pic b at 10σ significance.

3. While the planet will not transit, we are confident the Hill sphere around β Pic b will transit. The Hill sphere will begin transit at the start of April in 2017 and finish transiting at the end of January in 2018 with closest approach in the end of August in 2017. The transit of β Pic b's Hill sphere should be our best chance in the near future to investigate young circumplanetary material.

In the future, this MCMC forward modeling technique can be applied to photometry and spectral extraction alongside of astrometry of directly imaged exoplanets, allowing for improved characterization of their atmospheres. For β Pic b, continued monitoring of its orbit will yield more insight into the dynamics of the star system, although the planet will soon be too close to its star to be seen with current direct imaging instrumentation. However, once the planet appears on the other side, continued astrometric monitoring should be able to constrain the semimajor axis and eccentricity of the orbit much better, which will improve our understanding of how β Pic b perturbs the disk and if there are other planets perturbing β Pic b.

The GPI project has been supported by Gemini Observatory, which is operated by AURA, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the NSF (USA), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia), MCTI (Brazil) and MINCYT (Argentina). We thank the National Science Foundation (NSF AST-1518332) and NASA (NNX15AC89G and NExSS program NNX15AD95G) for contributing to support of this research. 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 research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. The posterior distribution plots were made with corner (Foreman-Mackey et al. 2016).

# APPENDIX A
## ANALYTIC FORWARD MODELING OF THE PLANET PSF AFTER STELLAR PSF SUBTRACTION

Using the notation of Pueyo (2016), we can analytically forward model the PSF of the planet through the process of stellar PSF subtraction. Let us denote the target image as $T_{λ_p, t_p}(\mathbf{x})$, which is taken at wavelength $λ_p$ and time $t_p$ and contains the pixels $\mathbf{x}$ from which we want to subtract the stellar PSF.

To forward model the PSF of the planet, we need a model of the PSF of the planet prior to stellar PSF subtraction. For our GPI data, we use the satellite spots to generate a realization of

an unocculted and unprocessed point source as a function of wavelength. Using this model, we can generate $A_{λ_p, t_p}(\mathbf{x})$, a frame that consists solely of the unprocessed PSF of the planet in the target frame. The initial estimate for the planet's position in the frame is determined using an input separation and PA of the planet and the location of the star. Similarly, we can also generate an array of frames, each containing the unprocessed PSF of the planet for each image in the reference library, accounting for the fact that the planet position changes due to ADI field rotation and due to SDI rescaling to align speckles. In these frames, the PSF shape only depends on the wavelength of the frame and the position of the PSF depends on the apparent movement of the planet in the data due to ADI and SDI. This array of unprocessed PSFs of the planet will be used to calculate the perturbation of the KL modes due to the existence of the planet.

Assuming a spectrum, $f$, of the planet in the case of IFS data (for imaging data, $f$ is just a scalar), we can compute $\Delta Z_k(\mathbf{x})$, the perturbation on the $k$th KL mode, $Z_k(\mathbf{x})$, due to the existence of a planet, using equation E18 or E20 of Pueyo (2016). We will not write out explicitly the exact formula to compute $\Delta Z_k$ here since we have not made any modifications to it. We note that we use Equation E18 from Pueyo (2016) to compute $\Delta Z_k$ as it is computationally faster than Equation E20 for a fixed input spectrum. We calculate the perturbed $\Delta Z_k$ for the first $k_{Klip}$ KL modes, where $k_{Klip}$ is the number of KL modes we choose to use in our PSF subtraction.

From here, the forward-modeled PSF of the planet for this frame is computed as

$$F_{λ_p, t_p}(\mathbf{x}) = A_{λ_p, t_p}(\mathbf{x}) - \left[ \begin{aligned} & \sum_{k=1}^{K_{Klip}} \langle A_{λ_p, t_p}(\mathbf{x}), Z_k(\mathbf{x}) \rangle > Z_k(\mathbf{x}) \\ & - \sum_{k=1}^{K_{Klip}} \langle T_{λ_p, t_p}(\mathbf{x}), Z_k(\mathbf{x}) \rangle > \Delta Z_k(\mathbf{x}) \\ & + \langle T_{λ_p, t_p}(\mathbf{x}), \Delta Z_k(\mathbf{x}) \rangle > Z_k(\mathbf{x}) \end{aligned} \right], \quad (7)$$

where $\langle \cdot, \cdot \rangle$ is the inner product. Equation (7) is very similar to Equation F7 in Pueyo (2016), but is focused on generating the forward model with a fixed input spectrum and not concerned with keeping the planet's spectrum as a free parameter for spectral extraction. As mentioned in Pueyo (2016), the term in the first square bracket is the over-subtraction term that is due to the projection of the KL modes on data with a planet in it, and the terms in the second square bracket are the self-subtraction terms due to the presence of a planet in the reference images influencing the KL modes.

After generating forward models for each frame, we take the mean of all the forward models in time to give us a single forward-modeled PSF cube, $F_{λ_p}(\mathbf{x})$ that still has spectral information in the third dimension. In principle, this PSF cube can be used to retrieve spectral information from the data, but it is outside of the scope of this paper. For our astrometry purposes, we also take the mean in wavelength, generating a single 2D forward-modeled PSF, which we call $F$, for each data set.
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 7.** Data, best-fit forward model, and residual map after the model as been subtracted from the data for each of the 12 data sets. In each row, we plot two data sets. For each data set, we plot the data (left), best-fit forward model (center), and residual map (right) on the same color scale. While the scale of each data set is different, zero is mapped to the same color for all the data sets.</figcaption></figure>

## APPENDIX B
BKA FIT RESIDUALS

In Figure 7, we plot the residuals to the best-fit model for each epoch using BKA. As β Pic b moves closer to its host star, the magnitude of speckle noise increases relative to the signal of β Pic b. This is especially true in the last epoch when β Pic b was closest to its star. However, the fit was still accurate as the bright positive core of the PSF of the planet was successfully modeled and not seen in the residuals, which are consistent with noise.

## REFERENCES

* Beust, H., & Morbidelli, A. 2000, *Icar*, 143, 170

* Binks, A. S., & Jeffries, R. D. 2016, *MNRAS*, 455, 3345

* Bottom, M., Kuhn, J., Mennesson, B., et al. 2015, *ApJ*, 809, 11

* Brandt, T. D., McElwain, M. W., Turner, E. L., et al. 2013, *ApJ*, 764, 183

* Burrows, C. J., Krist, J. E., Stapelfeldt, K. R. & WFPC2 Investigation Definition Team 1995, *BAAS*, 27, 32.05

* Cantalloube, F., Mouillet, D., Mugnier, L. M., et al. 2015, *A&A*, 582, A89

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

* Chilcote, J., Barman, T., Fitzgerald, M. P., et al. 2015, *ApJL*, 798, L3

* Crifo, F., Vidal-Madjar, A., Lallement, R., Ferlet, R., & Gerbaldi, M. 1997, *A&A*, 320, L29

* Czekala, I., Andrews, S. M., Mandel, K. S., et al. 2015, *ApJ*, 812, 128

* Dawson, R. I., Murray-Clay, R. A., & Fabrycky, D. C. 2011, *ApJL*, 743, L17

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

* Defrère, D., Lebreton, J., Le Bouquin, J.-B., et al. 2012, *A&A*, 546, L9

* Esposito, T. M., Fitzgerald, M. P., Graham, J. R., & Kalas, P. 2014, *ApJ*, 780, 25

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

* Foreman-Mackey, D., Vousden, W., Price-Whelan, A., et al. 2016, corner.py, v1.0.2, Zenodo, doi: 10.5281/zenodo.45906

* Galicher, R., & Marois, C. 2011, in Second Int. Conf. on Adaptive Optics for Extremely Large Telescopes (Victoria, BC: Victoria Conference Centre), 25 http://ao4elt2.lesia.obspm.fr

* Goodman, J., & Weare, J., 2010, *Comm. App. Math. Comp. Sci.*, 5, 65

* Graham, J. R. 2009, Gemini Planet Imager Operational Concept Definition Document V4.1 (Hilo, HI: Gemini Observatory)

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

* Heap, S. R., Lindler, D. J., Lanz, T. M., et al. 2000, *ApJ*, 539, 435

* Ingraham, P., Perrin, M. D., Sadakuni, N., et al. 2014, *Proc. SPIE*, 9147, 914770

* Kalas, P., Graham, J. R., Chiang, E., et al. 2008, *Sci*, 322, 1345

* Kenworthy, M. A., & Mamajek, E. E. 2015, *ApJ*, 800, 126

* Koen, C., Balona, L. A., Khadaroo, K., et al. 2003, *MNRAS*, 344, 1250

* Konopacky, Q. M., Thomas, S. J., Macintosh, B. A., et al. 2014, *Proc. SPIE*, 9147, 914784

* Lafrenière, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, É. 2007, *ApJ*, 660, 770

* 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

* Lecavelier des Etangs, A., & Vidal-Madjar, A. 2016, *A&A*, 588, A60

* Lecavelier des Etangs, A., Vidal-Madjar, A., Burki, G., et al. 1997, *A&A*, 328, 311

* Macintosh, B. A., Graham, J. R., Ingraham, P., et al. 2014, *PNAS*, 111, 35
Mamajek, E. E., & Bell, C. P. M. 2014, MNRAS, 445, 2169

Mamajek, E. E., Quillen, A. C., Pecaut, M. J., et al. 2012, AJ, 143, 72

Marois, C., Doyon, R., Racine, R., & Nadeau, D. 2000, PASP, 112, 91

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

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

Marois, C., Macintosh, B., & Véran, J.-P. 2010, Proc. SPIE, 7736, 77361J

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

Mouillet, D., Larwood, J. D., Papaloizou, J. C. B., & Lagrange, A. M. 1997, MNRAS, 292, 896

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

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

Pueyo, L. 2016, ApJ, 824, 117

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

Rameau, J., Nielsen, E. L., De Rosa, R. J., et al. 2016, ApJL, 822, L29

Shen, Y., & Tremaine, S. 2008, AJ, 136, 2453

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

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

Thébault, P., & Beust, H. 2001, A&A, 376, 621

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

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, Astrophysics Source Code Library, pyKLIP, ascl:1506.001

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

Zanazzi, J. J., & Lai, D. 2016, arXiv:1605.02365