# Multiband GPI Imaging of the HR 4796A Debris Disk

Christine Chen, Johan Mazoyer, Charles Poteet, Bin Ren, Gaspard Duchêne, Justin Hom, Pauline Arriaga, Maxwell Millar-Blanchaer, Jessica Arnold, Vanessa Bailey, et al.

► To cite this version:

Christine Chen, Johan Mazoyer, Charles Poteet, Bin Ren, Gaspard Duchêne, et al. Multiband GPI Imaging of the HR 4796A Debris Disk. The Astrophysical Journal, 2020, 898 (1), pp.55. 10.3847/1538-4357/ab9aba. hal-03038969

HAL Id: hal-03038969

https://hal.science/hal-03038969v1

Submitted on 27 Sep 2024

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>

Multiband GPI Imaging of the HR 4796A Debris Disk

Christine Chen<sup>1,2</sup>①, Johan Mazoyer<sup>3,4,26</sup>①, Charles A. Poteet<sup>1</sup>①, Bin Ren<sup>2,5,6</sup>①, Gaspard Duchêne<sup>7,8</sup>①, Justin Hom<sup>9</sup>①,
Pauline Arriaga<sup>10</sup>①, Maxwell A. Millar-Blanchaer<sup>5,11</sup>①, Jessica Arnold<sup>12</sup>①, Vanessa P. Bailey<sup>11</sup>①, Juan Sebastián Bruzzone<sup>13</sup>①,
Jeffrey Chilcote<sup>14</sup>①, Élodie Choquet<sup>15</sup>①, Robert J. De Rosa<sup>16</sup>①, Zachary H. Draper<sup>17</sup>①, Thomas M. Esposito<sup>7</sup>①,
Michael P. Fitzgerald<sup>10</sup>①, Katherine B. Follette<sup>18</sup>①, Pascale Hibon<sup>16</sup>①, Dean C. Hines<sup>1</sup>①, Paul Kalas<sup>7,19,20</sup>①,
Franck Marchis<sup>19</sup>①, Brenda Matthews<sup>19</sup>, Julien Milli<sup>21</sup>①, Jennifer Patience<sup>9</sup>, Marshall D. Perrin<sup>1</sup>①, Laurent Pueyo<sup>1</sup>,
Abhijith Rajan<sup>1</sup>①, Fredrik T. Rantakyrö<sup>22</sup>①, Timothy J. Rodigas<sup>12</sup>①, Gael M. Roudier<sup>11</sup>①, Glenn Schneider<sup>23</sup>①,
Rémi Soumer<sup>1</sup>①, Christopher Stark<sup>1</sup>, Jason J. Wang<sup>5,27</sup>①, Kimberly Ward-Duong<sup>9,18</sup>①, Alycia J. Weinberger<sup>12</sup>①,
David J. Wilner<sup>24</sup>①, and Schuyler Wolff<sup>25</sup>①

<sup>1</sup> Space Telescope Science Institute (STScI), 3700 San Martin Drive, Baltimore, MD 21218, USA; cchen@stsci.edu

<sup>2</sup> Department of Physics and Astronomy, The Johns Hopkins University, 3701 San Martin Drive, Baltimore, MD 21218, USA

<sup>3</sup> NASA Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA

<sup>4</sup> LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Jules Janssen, F-92195 Meudon, France

<sup>5</sup> Department of Astronomy, California Institute of Technology, 1200 East California Boulevard, Pasadena, CA 91125, USA

<sup>6</sup> Department of Applied Mathematics and Statistics, The Johns Hopkins University, 3400 North Charles Street, Baltimore, MD 21218, USA

<sup>7</sup> Astronomy Department, University of California, Berkeley, CA 94720, USA

<sup>8</sup> Université Grenoble-Alpes, CNRS Institut de Planétologie et d'Astrophysique (IPAG), F-38000 Grenoble, France

<sup>9</sup> School of Earth and Space Exploration, Arizona State University, P.O. Box 871404, Tempe, AZ 85287, USA

<sup>10</sup> Department of Physics & Astronomy, 430 Portola Plaza, University of California, Los Angeles, CA 90095, USA

<sup>11</sup> NASA Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

<sup>12</sup> Department of Terrestrial Magnetism, Carnegie Institution for Science, 5421 Broad Branch Road, Washington, DC 20015, USA

<sup>13</sup> Department of Physics and Astronomy, The University of Western Ontario, London, ON N6A 3K7, Canada

<sup>14</sup> Department of Physics, University of Notre Dame, 225 Nieuwland Science Hall, Notre Dame, IN 46556, USA

<sup>15</sup> Aix Marseille Univ., CNRS CNES, LAM, Pôle de l'Etoile Site de Château-Gombert 38, rue Frédéric Joliot-Curie F-13388 Marseille Cedex 13 France

<sup>16</sup> European Southern Observatory, Alonso de Córdova 3107, Vitacura, Santiago, Chile

<sup>17</sup> University of Victoria, 3800 Finnerty Rd., Victoria, BC V8P 5C2, Canada

<sup>18</sup> Physics and Astronomy Department, Amherst College, 21 Merrill Science Drive, Amherst, MA 01002, USA

<sup>19</sup> SETI Institute, Carl Sagan Center, 189 Bernardo Avenue, Mountain View, CA 94043, USA

<sup>20</sup> Institute of Astrophysics, FORTH, GR-71110 Heraklion, Greece

<sup>21</sup> European Southern Observatory (ESO), Alonso de Córdova 3107, Vitacura, Casilla, 19001, Santiago, Chile

<sup>22</sup> Gemini Observatory, Casilla 603, La Serena, Chile

<sup>23</sup> Steward Observatory, The University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA

<sup>24</sup> Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

<sup>25</sup> Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands

Received 2020 March 6; revised 2020 June 5; accepted 2020 June 8; published 2020 July 23

Abstract

We have obtained Gemini Planet Imager (GPI) J-, H-, K1-, and K2-Spec observations of the iconic debris ring around the young, main-sequence star HR 4796A. We applied several point-spread function (PSF) subtraction techniques to the observations (Mask-and-Interpolate, RDI-NMF, RDI-KLIP, and ADI-KLIP) to measure the geometric parameters and the scattering phase function for the disk. To understand the systematic errors associated with PSF subtraction, we also forward-modeled the observations using a Markov Chain Monte Carlo framework and a simple model for the disk. We found that measurements of the disk geometric parameters were robust, with all of our analyses yielding consistent results; however, measurements of the scattering phase function were challenging to reconstruct from PSF-subtracted images, despite extensive testing. As a result, we estimated the scattering phase function using disk modeling. We searched for a dependence of the scattering phase function with respect to the GPI filters but found none. We compared the H-band scattering phase function with that measured by Hubble Space Telescope STIS at visual wavelengths and discovered a blue color at small scattering angles and a red color at large scattering angles, consistent with predictions and laboratory measurements of large grains. Finally, we successfully modeled the SPHERE H2 HR 4796A scattered phase function using a distribution of hollow spheres composed of silicates, carbon, and metallic iron.

Unified Astronomy Thesaurus concepts: Coronagraphic imaging (313); Debris disks (363); Circumstellar matter (241); Planetary system formation (1257)

1. Introduction

During the past 20 years, visual, infrared, and millimeter
observations have provided incontrovertible evidence that most
stars are surrounded at birth by circumstellar accretion disks, and

at least some disks have built the ~3000 extrasolar planets that
have been discovered thus far (Williams & Cieza 2011). In
general, high angular resolution visual to near-infrared images
show starlight that is scattered off of predominantly submicron-
sized dust grains while those at mid-infrared to millimeter
wavelengths show thermal emission from increasingly larger and
larger dust grains (Hughes et al. 2018). At visual to near-infrared
wavelengths, observations suffer from the complication that host

<footnote n="26"><sup>26</sup> NFHP Sagan Fellow.</footnote>

<footnote n="27"><sup>27</sup> 51 Pegasi b Fellow.</footnote>
stars are bright compared to circumstellar disks; thus, these
observations benefit from the application of high-contrast
imaging techniques to reveal faint scattered light (Schneider
et al. 1999). High-contrast imaging techniques typically use
point-spread function (PSF) subtraction to remove light in
conjunction with coronagraphs that suppress light from the
bright central star (Lafrenière et al. 2009; Soummer et al. 2012).
Indeed, high-contrast imaging studies have discovered planetary
mass companions in a handful of protoplanetary and debris disks
thus far (Lagrange et al. 2010; Marois et al. 2010; Meshkat et al.
2015).

In disk imaging, disk geometry and surface brightness
measurements can shed light on the dynamics and composition
of the circumstellar dust and therefore the planets in the
underlying planetary systems. High-resolution images of debris
disks have revealed that the dust in many systems appears
sculpted into rings. Some rings are elliptical. In these systems,
undetected planets on elliptical orbits may force the
eccentricities of the parent bodies and dust into elliptical orbits
(Wyatt et al. 1999). High-resolution images have also revealed
local brightness enhancements, or clumps, and warps in some
disks, the properties of which have been used to effectively
predict the presence of an exoplanet (Mouillet et al. 1997).
However, these clumps may also be generated by giant
collisions. In our solar system, giant collisions between
forming planets and Mars-sized impactors are believed to have
ejected Mercury's mantle, formed the Moon, and created the
large crater on Mars' northern hemisphere. Such events are
expected to produce copious quantities of dust initially in very
localized areas. Finally, high-resolution imaging enables
measurements of the scattering phase function (SPF). Measure-
ments of the SPF have been used to constrain grain size (Milli
et al. 2017, hereafter M17). Unfortunately, PSF subtraction
techniques impact the field of view differently depending on
the angular distance from the central star and the position
angle. For companions, the impact on point-source photometry
and astrometry is relatively minor. However, for disks that are
spatially extended (and can fill the field of view), the impact
can be significant. For example, Angular Differential Imaging
(ADI; Marois et al. 2006) produces regions of self-subtraction
around a disk that strongly impact the disk shape and local
photometry (Milli et al. 2012).

In 2014, Gemini Observatory commissioned the Gemini
Planet Imager (GPI), a second-generation high-contrast ima-
ging instrument that provides near-infrared integral field
spectroscopy and polarimetry (Macintosh et al. 2014; Perrin
et al. 2015) on the Gemini South Telescope. Although GPI was
primarily designed to search for and characterize Jovian-mass
planets, it has spatially resolved a handful of debris disks for
the first time, particularly in the nearby Sco-Cen OB
association (Kalas et al. 2015; Draper et al. 2016). To date,
GPI studies have focused on characterizing the gross morph-
ology of disks because understanding the impact of PSF
subtraction on high-contrast imaging data is challenging. We
are carrying out a multifilter, integral field spectroscopy and
polarimetry study of approximately one dozen bright debris
disks that have been spatially resolved using Hubble Space
Telescope (HST) as part of the 2015B Gemini Large and Long
Program "Characterizing Dusty Debris in Exoplanetary
Systems" (PI Chen). Our goal is to not only extract

measurements of the disk geometry and surface brightness
but also to combine multifilter observations to provide the best
constraints on the dust grain properties (e.g., size, porosity,
shape, and composition). In 2014, the European Southern
Observatory commissioned the Spectro-Polarimetric High-
contrast imager for Exoplanets REsearch (SPHERE; Beuzit
et al. 2019) instrument, a second-generation high-contrast
imaging instrument for the Very Large Telescope. SPHERE
addresses many of the same science goals as GPI with slightly
different instrument capabilities.

The A0V star HR 4796A at a distance *d* ~ 72.8 pc possesses a spectacular narrow, inclined ring (*i* = 76°) imaged in scattered light (Schneider et al. 2009, 2018; Perrin et al. 2015; Rodigas et al. 2015) and thermal emission (Telesco et al. 2000; Kennedy et al. 2018) with a semimajor axis *a* ~ 72 au. An age of ~8 ± 2 Myr has been estimated for the central star based on the lithium abundance and isochrone fitting of its M-type companion (Stauffer et al. 1995) that is located 7"/7 from the primary star. Since the collisional and Poynting-Robertson drag lifetimes of the circumstellar dust are substantially shorter than the age of the star, the dust in this system is believed to be replenished by collisions among parent bodies (Jura et al. 1993). The lack of circumstellar dust close to the star has led to speculation that there are planetary-mass companions sculpting the disk (Jura et al. 1995). However, no companions have yet been detected to a completeness limit of ~4 *M*<sub>Jup</sub> (Milli et al. 2017); for reference, a Neptune-mass planet is sufficient to maintain the sharp inner edge of this disk. Recent, deep HST STIS imaging has revealed that the HR 4796A ring is located within a substantially more spatially extended population of circumstellar dust that may be interacting with the local interstellar medium (Schneider et al. 2018). Recent SPHERE IRDIS and ZIMPOL observations (M17; Milli et al. 2019) have provided exquisite measurements of the disk geometry, infrared SPF, and visual polarized intensity phase function. Most interestingly, Milli et al. (2019) show that the north-south asymmetry in the HR 4796 disk cannot be explained exclusively by a geometrical effect such as pericenter glow (introduced by Wyatt et al. 1999), but it is probably due to a dust density enhancement at pericenter (Olofsson et al. 2019). Finally, no circumstellar gas has yet been detected (Chen & Kamp 2004; Kennedy et al. 2018) around HR 4796A.

We report the results from a GPI integral field spectroscopy
study of the iconic HR 4796A debris ring. We PSF-subtract our
observations using multiple algorithms leveraging both ADI
and Reference Differential Imaging (RDI) techniques to
explore the benefits and weaknesses of each algorithm. For
comparison, we also fit our total intensity images using a
simple disk model and a Markov Chain Monte Carlo (MCMC)
approach. We recover measurements of the disk geometry and
SPF that are consistent with results from other groups using
other high-contrast imaging instruments and PSF subtraction
techniques. We find that the HR 4796A disk may be slightly
less forward scattering than previously reported, as well as
evidence that the near-infrared to visual color of the HR 4796A
SPF varies as a function of phase angle with a slightly blue
color at relatively small scattering angles (<60°) and a slightly
red color at relatively large scattering angles (>120°),
consistent with measurements of large particles in the
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Slices of calibrated GPI exposures of the HR 4796A disk obtained using *J*-Spec (left), *H*-Spec (middle), and *K1*-Spec (right). The contrast between the disk and the bright speckles is poorer in *J* and *H* bands, where the speckles are brighter. The faint ring visible in the calibrated *K1*-Spec exposure is the bright HR 4796A ring. The *J*-band diffraction pattern is complicated with the first-order satellite, waffle, and second-order satellite spots visible as three spots near the corners of the detector. At longer wavelengths, the angular size of the diffraction pattern increases. Only the first-order and waffle spots are visible in *H*. Only the first-order spots are visible in *K1* and *K2*.</figcaption></figure>

laboratory. We demonstrate that a population of large, irregular particles can reproduce the total intensity SPF and the measured near-infrared to visual color of the SPF as a function of phase angle.

## 2. Observations and Data Reduction

### 2.1. *Observations*

We obtained GPI Integral Field Spectrograph (Spec) observations of the debris disk around HR 4796A using the *J*-Spec, *H*-Spec, *K1*-Spec, and *K2*-Spec observing modes. The *J*-Spec data were obtained as part of the “Debris Characterization in Exoplanetary Systems” Gemini Large and Long Program (PI C. Chen; GS-2016A-LP-6); the *H*-Spec data were obtained as part of the “GPI Exoplanet Survey” (PI B. Macintosh; GS-2016A-Q-500); and the *K1*- and *K2*-Spec data were obtained as part of a Principal Investigator Program “Does the HR 4796 Debris Disk Contain Icy Grains?” (PI C. Chen; GS-2015A-Q-27). GPI provides 14 mas sampling for all coronagraphic observing modes. The IFS provides observations at 1.12–1.35 µm, 1.50–1.80 µm, 1.9–2.19 µm, and 2.13–2.4 µm, respectively, with spectral resolutions $R(=λ/Δλ) \sim 35–39, 44–49, 62–70, \text{ and } 75–83$, respectively.

For each observation, the bright, primary star was centered behind the focal plane mask. On-axis light from the star was not only reduced but also diffracted by a grid imprinted on the pupil plane mask, generating a diffraction pattern including astrometric reference or satellite images of the primary star that are ~10,000 times fainter than the star itself. These spots are used not only to determine the exact location of the occulted star but also to flux-calibrate the astronomical scene (Sivaramakrishnan & Oppenheimer 2006; Wang et al. 2014). The distance of the satellite spots from the primary star scales with wavelength such that the spots appear closer to the star at shorter wavelengths. For *J*-band observations, the coronagraphic images are complex, revealing not only lower Strehl ratios but also additional PSF structure such as AO “waffle-mode” spots, second-order satellite spots, and diffractive spots from the Deformable Mirror (DM) actuator print-through (see Figure 1). To facilitate ADI analysis of the data, the instrument field of view was allowed to rotate during an observation, producing diversity in the orientation of the science target with respect to the instrument reference frame. In some of our

individual exposures, the satellite and waffle spots overlapped with the disk, making PSF subtraction more challenging.

Each of our observations was composed of several tens of individual exposures with integration times between ~30 and ~90 s. The frame times were selected to be long enough such that the detector readout noise did not dominate the disk signal but short enough to avoid saturation of residual speckles and angular smearing of the disk. *K1*-Spec and *K2*-Spec observations typically tolerate longer integration times because the star and therefore the speckles are fainter and the PSF is more oversampled. In general, most of our GPI observations were obtained in contiguous time blocks near transit to maximize the field rotation. However, our 2014 March 25 and 2015 April 2 *K1* data sets were not. In 2014, the target was observed for 10 minutes; the observations were paused 15 minutes before transit and resumed again 5 minutes after transit; after transit, the target was observed for an additional 10 minutes. In 2015, the target was observed for an hour; the observations were paused 2 hours before transit and resumed 2 hours after transit; after transit, the target was observed for an additional hour. Since the night sky is bright at *K* band, we obtained dedicated sky observations (20" offset from our target) for thermal/sky background subtraction in *K1* and *K2* bands. We typically took five sky exposures after every hour of *K1*- or *K2*-band observing.

A GPI observing mode usually corresponds to a choice of filter, apodizer, and focal plane, and Lyot plane masks. However, GPI lost the ability to reliably change apodizers from 2015 July until 2016 March; therefore, the majority of the observations obtained during this time used the *H*-band optimized APOD_H_G6205 independent of the wavelength of the observations made. For spectroscopic observations, using the APOD_H_G6205 apodizer for any mode other than *H*-Spec is considered “Nonstandard.” Nonstandard observations possess slightly different throughputs and inner working angles (IWAs). While our *H*, *K1*, and *K2* observations were obtained in Standard configurations, our *J* observations were not. Since we estimate the PSF from the satellite spots in our forward-modeling analysis, the nonstandard *J*-band observation does not affect our analysis. We list the integration times, the number of science exposures, the average air mass of the target, the seeing, and the amount of field rotation in Table 1.

We note that HR 4796 was observed using the *K1* IFS on three different nights spanning a period of more than 1 yr.
<table><caption>**Table 1**
Observations</caption><thead><tr><th>Target</th><th>UT Date</th><th>Program</th><th>Obs. Mode<sup>a</sup></th><th>Int.Time (s)</th><th>No. of Exps.</th><th>Air Mass</th><th>Seeing (arcsec)</th><th>Field Rot.<sup>b</sup> (deg)</th></tr></thead><tbody><tr><td>HR 4796A</td><td>2016 Mar 23</td><td>GS-2016A-LP-6</td><td>J-Spec<sup>c</sup></td><td>29.1</td><td>59</td><td>1.02–1.04</td><td>0.4–0.7</td><td>48.8</td></tr><tr><td>HR 4796A</td><td>2016 Mar 18</td><td>GS-2016A-Q-500</td><td>H-Spec</td><td>59.6</td><td>37</td><td>1.01–1.02</td><td>0.5</td><td>52.7</td></tr><tr><td>HR 4796A</td><td>2014 Mar 25</td><td>GS-ENG-GPI-COM</td><td>K1-Spec</td><td>59.6</td><td>20</td><td>1.01–1.02</td><td>0.3</td><td>15.1, 12.6</td></tr><tr><td>HR 4796A</td><td>2015 Apr 2</td><td>GS-2015A-Q-27</td><td>K1-Spec</td><td>88.7</td><td>47</td><td>1.13–1.24, 1.11–1.21</td><td>0.55</td><td>8.0, 8.4</td></tr><tr><td>HR 4796A</td><td>2015 Apr 3</td><td>GS-2015A-Q-27</td><td>K1-Spec</td><td>88.7</td><td>46</td><td>1.01–1.09</td><td>0.55</td><td>78.5</td></tr><tr><td>HR 4796A</td><td>2015 Apr 3</td><td>GS-2015A-Q-27</td><td>K2-Spec</td><td>88.7</td><td>23</td><td>1.02–1.04</td><td>0.55</td><td>33.8</td></tr></tbody></table>

**Notes.**

<sup>a</sup> GPI instrument mode corresponding to choice of filter plus corresponding focal plane and Lyot plane masks and apodizer.

<sup>b</sup> Change in parallactic angle over the course of the observation sequence.

<sup>c</sup> Observation made with nonstandard apodizer (APOD_H_G6205).

Thus, our observations provide some empirical insight into the repeatability of GPI *K1* IFS disk observations.

## 2.2. Data Reduction

Reduction and interpretation of ground-based coronagraphic observations are challenging because observing conditions change on rapid timescales. The disk around HR 4796A is so bright that it can be seen in the individual exposures before PSF subtraction (Figure 1). We used the GPI Data Reduction Pipeline<sup>28</sup> (DRP; Perrin et al. 2014) to reduce, wavelength-calibrate, and assemble spectral data cubes (which contain 37 2D images) at the individual exposure level. Specifically, for the *J* and *H* data, we used the DRP to subtract the dark background, update spot shifts for flexure, interpolate bad pixels in the 2D frame, assemble the IFS spectral data cube, interpolate bad pixels in the cube, correct distortion, and measure satellite spot locations. For the *K1* and *K2* data, we also subtracted dedicated thermal/sky background images and destriped the science images after dark subtraction but before interpolating bad pixels in the 2D frame. For *K1* and *K2* bands, with high thermal noise, we removed some spectral slices from the data cubes if the signal-to-noise ratio (S/N) of the satellite spots was lower than 2. We determined the location of the occulted star in each image cube using a least-squares fit to all of the satellite spot positions. We then used the star’s position to align all the images to a common frame before PSF subtraction. Our method allowed us to estimate the position of the star in each frame to a precision of 0.05 pixels or 0.7 mas (Wang et al. 2014). Next, we collapsed our IFS cubes in wavelength space to produce a single broadband image for each exposure and removed exposures in which the disk overlaps with one or two satellite spots. Next, we generated a single high-S/N empirical PSF for each observation, using all of the remaining satellite spots. Finally, we rotated the images so that north was up and east was to the left.

## 2.3. PSF Subtraction

Several methods have been developed to remove stellar speckles from coronagraphic images in post-processing. These techniques may severely impact circumstellar objects, especially extended objects (debris and protoplanetary disks), for which the photometry and the shape are impacted (Milli et al. 2012). In this work, we used four different methods to estimate the PSF for PSF subtraction. Our analysis allowed us to

quantify the relative strengths and weaknesses of each algorithm for measuring the geometry and phase function of a bright, narrow, inclined ring: (1) Mask-and-Interpolate (Perrin et al. 2015), (2) RDI-Non-negative Maxmatrix Factorization (NMF; Ren et al. 2018), (3) RDI-Karhunen-Loève Image Projection (KLIP; Soummer et al. 2012), and (4) ADI-KLIP (Marois et al. 2006).

The Mask-and-Interpolate method takes advantage of the high S/N and rather smooth nature of the residual PSF halo. For each exposure: (1) The disk is masked, assuming that it can be circumscribed by two nested ellipses with a maximum radial width of 15 pixels (210 mas) at the ansae. In addition, the first, second, and waffle-mode satellite spots are also masked. (2) The PSF is estimated under the mask using a fourth-order polynomial interpolation in both the *x*- and *y*-directions. The masked disk is replaced by the interpolated PSF. (3) The resulting PSF is smoothed using a median filter with a smoothing length of 7 pixels and subtracted from the original image. (4) The PSF-subtracted images are median-combined. (5) The remaining background in each PSF-subtracted image is fit in concentric annuli within a wedge ±85° from the major axis or excluding ±5° from the minor axis using a fourth-order azimuthal polynomial (see Figure 2). This background model is then subtracted from the smoothed PSF-subtracted image to produce the final image. We present the total intensity images for each of our observations and their corresponding S/N maps (formed by dividing our total intensity image by our uncertainty image) in Figure 2.

In RDI, independent observations of bright stars (ideally without companions or disks) are used to create a PSF library from which a representative PSF is estimated. The GPI Exoplanet Survey (GPIES) team is carrying out a search for planetary-mass companions around ~600 young stars using the GPI *H*-Spec mode. To date, the team has published observations of the first ~300 stars (Nielsen et al. 2019). The majority of the GPIES observations contain isolated stars, and all the observations for which no companion, disk, or background object has been detected can be used as a PSF reference library. We cross-correlate our HR 4796A Spec *H*-Spec observations with those in the library to select the 100 exposures that are most correlated with our target exposures. Then, we perform both an NMF and a KLIP principal component analysis on the empty reference observations to reconstruct the PSF. For RDI-NMF, the main parameter is the number of components used. We use 15 NMF components. For RDI-KLIP, the main parameter is the number of eigenmodes used in the Karhunen-Loève (KL) basis. We use five KL mode basis vectors. For both

<footnote n="28"><sup>28</sup> http://www.stsci.edu/~mperrin/software/gpidata/</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Mask-and-Interpolate total intensity images and S/N maps of the HR 4796A debris disk, oriented with north up and east to the left.</figcaption></figure>

RDI-NMF and RDI-KLIP, we perform the analysis on the whole image rather than concentric annuli to enable the best measurements of relative surface brightness within the disk. We present the *H*-band total intensity images for our RDI-NMF and RDI-KLIP reductions in Figure 3.

In ADI, the astrophysical scene is allowed to rotate with respect to the diffraction pattern of the occulted star. In this case, a representative PSF can be estimated directly from the observing sequence, including the astrophysical target. To minimize self-subtraction, we added an additional parameter compared with RDI-KLIP: the exclusion criterion (equivalent to $N_δ$ in Lafrenière et al. 2007), the minimum number of pixels a hypothetical point-like astrophysical source must move azimuthally to avoid overlap in two images. For all bands, we used three KL modes and an exclusion angle of $N_δ = 6°$,

which is a good compromise between S/N of detection and disk impact. To estimate the PSF using ADI-KLIP, we used the *pyKLIP*<sup>29</sup> package (Wang et al. 2015), a Python implementation of the KLIP algorithm. We applied *pyKLIP* globally (with the same parameters for the whole image) to avoid discontinuities in the reduction along the object.

We assumed that the uncertainty in our observations was dominated by speckle noise. Therefore, we estimated the uncertainty by masking out the disk and calculating the standard deviation of the remaining pixels in concentric annuli, assuming an annular width of 3 pixels, approximately the PSF FWHM for all observations and reduction methods. We present the total intensity images for each of our observations

<footnote n="29"><sup>29</sup> https://bitbucket.org/pyKLIP/pyklip</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** RDI-NMF (top left) and RDI-KLIP (top right) *H*-band total intensity images and S/N maps of the HR 4796A debris disk, oriented with north up and east to the left. The RDI-NMF reduction includes 15 NMF components. The RDI-KLIP reduction includes five KL mode basis vectors.</figcaption></figure>

(collapsed along the wavelength axis) and their corresponding S/N maps in Figure 4.

We generally detect the disk with S/N ~ 5–10. The S/N with which the disk is detected is higher in *H* and *K1* bands than in *J* and *K2* bands. The *J* data suffer from the largest speckle noise and the most overlap between the disk and PSF structures (e.g., satellite spots). The *K2* data suffer from a combination of lower throughput in the instrument (particularly at wavelengths 2.26–2.4 µm) and higher sky backgrounds. In addition, our *K2* observations have horizontal stripes, created by a misalignment of the focal plane mask. Consistent with previous ground- (Wahhaj et al. 2014; Rodigas et al. 2015, M17) and space-based observations (Schneider et al. 2009), we find that (1) the disk appears bright at the ansae, with the NE ansa ~10% brighter than the SW ansa and (2) the host star is offset from the center of the disk. Our *J* and *H* observations indicate that the disk surface brightness near the forward-scattering peak appears asymmetric, with the northern side of the disk brighter than the southern side.

The exact properties of the disk recovered are dependent on the PSF subtraction method used. For example, the disk recovered using “Mask-and-Interpolate” has the following features: (1) More spatially extended. ADI analyses produce disk self-subtraction for spatially extended disks with a finite spatial extent (Esposito et al. 2014). The detection of disk self-subtraction in our ADI-KLIP reduction is consistent with other

studies. (2) More affected by speckles. Interpolation is not expected to reproduce the complex GPI PSF structure with very high fidelity, particularly near the satellite and waffle spots. The lack of PSF fidelity is particularly problematic at *J* band, where the first- and second-order satellite spots, waffle spots, and DM spots are within the instrument field of view and overlap with the disk. To minimize the impact of satellite and waffle spots on the disk, we masked out areas of the disk where the satellite and waffle spots overlap when calculating our average PSF-subtracted image. However, we note that despite taking this precaution, some satellite and waffle spot artifacts may still remain. For example, the waffle spots are immediately adjacent to the disk ansae in approximately half of the *J*-band exposures, leading to the interpolation being based on more distant pixels. Thus, the bright spots near the ansae in the *J*-band Mask-and-Interpolate image, which do not appear in the KLIP image, are suspect.

### 3. Analysis Tools

Historically, scattered-light images have provided some of the highest angular resolution images of debris disks, enabling detailed measurements of disk geometry and SPF. The latter measures the change in total intensity of the scattered light as a function of scattering angle, $\theta$, the angle of deviation from forward scattering (M17). However, accurately extracting the
<figure><img src="image_6.png" /><figcaption>**Figure 4.** ADI-KLIP total intensity images and S/N maps of the HR 4796A debris disk, oriented with north up and east to the left. These reductions were made using one angular zone and three KL mode basis vectors.</figcaption></figure>

SPF for a disk is challenging because PSF subtraction techniques alter the disk surface brightness differently, depending on the apparent angular separation from the star and the position angle. For example, ADI produces regions of oversubtraction around inclined disks, resulting in flux losses (Milli et al. 2012). Therefore, we used two techniques to estimate the disk geometry and empirical SPF: (1) Reconstruction—we measure the disk geometry and the SPF from our PSF-subtracted images. For the SPF, we correct for oversubtraction. (2) Modeling—we model the disk using a simple prescription for the geometry and SPF and use an MCMC method to explore the disk parameters until the residuals between the disk model and observation are minimized.

## 3.1. Reconstruction

### 3.1.1. Geometric Parameters

We used the *Debris Ring Analyzer (DRA)* to measure the projected disk parameters using the techniques described in Stark et al. (2014). Briefly, the DRA divides the disk into pie-shaped wedges centered on the star and then iteratively measures the peak of the radial surface brightness distribution in each wedge using a polynomial fit. The peak coordinates are transformed to Cartesian coordinates and then fit with an ellipse using the `mpfitellipse` code (IDL). The peak coordinates and ellipse fit are shown in Figure 7 for the *H* band for all the PSF subtraction methods. The DRA requires a mask to select the
areas of the image to be included in the geometric fit. We
created one mask for each PSF-subtracted image of HR 4796A
approximately following the area of the ring peak surface
brightness, excluding wedges near forward and backward
scattering where the speckle noise is the highest. We overlaid
the measurements of peak surface brightness on our total
intensity images and excluded points that did not lie on the
disk. Our Mask-and-Interpolate images required the largest
exclusion angles around forward and backward scattering,
while our RDI-KLIP images required the smallest. For
example, for our H-Spec observations, Mask-and-Interpolate
required 60° exclusion regions around forward and backward
scattering; RDI-NMF required 0° and 15° exclusion regions
around forward and backward scattering, respectively; RDI-
KLIP required no exclusion regions; and ADI-KLIP required
10° and 30° exclusion regions (see Figure 7). We divided the
remaining disk area into 2.7 wedges in the sky plane, centered
on the star. For each wedge, we iteratively fit a third-order
polynomial to the disk surface brightness as a function of
distance to determine the radial location of the peak surface
brightness. The initial fit used the full radial extent of the
wedge and progressively selected smaller and smaller radial
sections of the wedge around the peak surface brightness until
its profile was well fit. We estimated the uncertainty in the peak
position by randomly varying the disk surface brightness
within a wedge according to the observed 1σ uncertainties
assuming a Gaussian distribution, refitting the radial peak, and
calculating the standard deviation of 50 unique instances.
Finally, we fitted an ellipse to the collection of peak surface
brightness coordinates and their uncertainties by minimizing
the perpendicular distance from the ellipse to the coordinates.

3.1.2. *Scattering Phase Function*

We hypothesize that we can reconstruct the intrinsic SPF by
correcting the SPF extracted from our PSF-subtracted images
for the effects of ADI subtraction at the one-dimensional level.
We note that the disk total intensity image appears bright near
the ansae as a result of limb brightening. This limb brightening
is generated by the large column density of scatterers near the
ansae in this optically thin disk, compared with that along
the minor axis, corresponding to the front and back sides of the
disk, as well as by the finite resolution of the telescope. When
extracting the SPF, we deproject the disk image and divide it
by the assumed model SPF to correct for limb brightening and
oversubtraction and therefore reveal the intrinsic SPF. This
approach has been used in previous studies of the HR 4796A
SPF (M17).

In our study, we use MCFOST (Pinte et al. 2006, 2009) to
simulate a disk with the same geometric parameters as HR
4796A and a known SPF. We match the disk surface brightness
to our observations of the HR 4796A ansa. We convolve our
idealized model disk with an empirical GPI PSF to mimic the
observed broadening of the disk ansae. Specifically, we use
unblocked Spec observations of the bright (K<sub>s</sub> = 7.6 mag),
nearby A-type star HD 118335 obtained using the J, H, K1, and
K2 filters and appropriate matching focal plane masks on 2014
March 26 to approximate the PSF. We inject our convolved
model disk into an empty GPI sequence (without a companion
or a disk) with similar noise properties, fixing the parallactic
angles in the empty data set to be the same as in our target
observations. We PSF-subtract the injected disk sequence using
the same PSF subtraction algorithm and parameters that we

<figure><img src="image_7.png" /><figcaption>**Figure 5.** SPFs extracted from a simulated observing sequence with an injected isotropic scattering disk ($g = 0$) and an HR 4796A-like geometry, corrected for limb brightening, using several different PSF subtraction techniques. The disk model is constructed so that the disk surface brightness $S/N$ is the same as the HR 4796A ansa. An isotropic disk is expected to have a constant SPF; thus, regions with values <1 are affected by oversubtraction.</figcaption></figure>

used for the HR 4796A sequence. We extract the SPF from the
PSF-subtracted, simulated disk image by performing aperture
photometry on the deprojected, $1/r^2$ corrected disk image.
Finally, we multiply the SPF extracted from the PSF-subtracted
HR 4796A image by the SPF assumed in our simulated disk
and divide by the SPF extracted from our PSF-subtracted,
simulated disk image. For simplicity, we simulate a disk with
an isotropic SPF or a Henyey-Greenstein (HG) scattering
parameter, $g = 0$. In this case, the SPF is constant as a function
of scattering angle.

We show how the various PSF subtraction algorithms impact the SPF extracted from a disk with an HR 4796A-like geometry by plotting the SPFs extracted from PSF-subtracted, isotropic disk sequences and corrected for limb-brightening effects in Figure 5. We find that Mask-and-Interpolate alters the SPF the least, with relatively little oversubtraction across a broad range of scattering angles (40°–145°). Mask-and-Interpolate produces the least oversubtraction because it uses an interpolation to estimate the PSF within the disk mask. RDI-NMF and RDI-KLIP produce more oversubtraction at a larger range of angles near forward and backward scattering. ADI-KLIP produces the most oversubtraction at a larger range of angles near forward and backward scattering and generates additional wavy structure beyond the oversubtraction at forward and backward scattering. ADI-KLIP relies on the observing sequence to estimate the PSF. In this case, large-scale, diffuse, disk scattered light could be misidentified as part of the PSF. We tested our method on a simulated data set in Section B.1 and demonstrated that SPF reconstruction should be feasible at scattering angles of 40°–140° for *H* and *K1* bands.

3.2. *Modeling*

An alternate method for estimating the SPF is fitting the total
intensity image using an MCMC approach combined with a
forward-modeling technique described in Pueyo (2016) and
assuming a simple prescription for the 3D dust density
distribution. Similar analyses have been used to estimate the
geometric properties of the several debris disks from GPI
observations (Millar-Blanchaer et al. 2016). We fit our ADI-KLIP total intensity images to estimate the SPFs from our GPI total intensity images.

### 3.2.1. Model Description

We assume that the dust density can be well described using a flared disk with a power-law surface density distribution

$$ \eta(r, z) \propto \left(\frac{r}{R_1}\right)^{-\beta} \exp\left[-\frac{1}{2}\left(\frac{z}{h_0 r}\right)^2\right] \quad (1) $$

(Millar-Blanchaer et al. 2015), where r is the radial distance of the dust from the star, z is the height of the dust above the disk midplane, and β is the power-law index of the dust density distribution. Since we study the relative brightness of the disk as a function of the disk and not the total disk mass, we do not solve for the porportionality constant in our analysis. Since scattered-light imaging of the HR 4796A debris disk indicates that it is a flat, narrow ring, we assume that the dust surface density $\eta(r, z) = 0$ interior to the ring ($r < R_1$) and exterior to the ring ($r > R_2$) and $h_0 = 0.01$. We note that $R_1$ (inner radius) and $R_2$ (outer radius) are distinct from the ellipse radius a measured in Section 4 that is extrapolated from points of maximum surface brightness calculated in pie-shaped wedges from ADI-KLIP, PSF-subtracted images. We initially left the scale height as a free parameter, and our best-fit models favored a very thin ring $h_0 < 0.001$ because the vertical extent of the disk was challenging to resolve. Unfortunately, the locations in the disk that provide the most leverage on the scale height measurements were along the disk minor axis where our data were most impacted by ADI self-subtraction. Since the disk scale height would be degenerate with the disk radial width parameters ($R_1$ and $R_2$) if the scale height were resolved, we decided to fix $h_0$ to a small value that was consistent with measurements on other disks. We note that small variations in $h_0$ would produce slightly different values for the disk radial width ($R_1$ and $R_2$).

Following M17, we combine the dust density distribution with a two-component HG SPF to generate a 2D projection of the disk from the observer's point of view. The two-component HG function

$$ p_2(g_1, \alpha, g_2, \theta) = \alpha HG(g_1, \theta) + (1 - \alpha) HG(g_2, \theta), \quad (2) $$

where $\theta$ is the scattering angle, HG(g, $\theta$) is the single-component HG function, $g_1$ and $g_2$ are the asymmetric scattering parameters for the two separate components, and $\alpha$ is the relative weight for the two components ($0 \le \alpha \le 1$). We note that the HG SPF function does not carry any physical meaning in our fits. Instead, we use HG SPFs to reproduce the essential behavior of complex SPFs using a small number of free parameters. For example, a single-component HG has one free parameter: the asymmetric scattering parameter g. A two-component HG function has three free parameters: the two asymmetric scattering parameters, $g_1$ and $g_2$, and the relative weight between them, $\alpha$. The HG SPF was first used to empirically model the SPF observed for interstellar dust (Henyey & Greenstein 1941). The majority of debris disk SPF studies in the literature use a single-component HG SPF

because most SPFs are not well measured over a large range of scattering angles. A two-component HG SPF has been successfully used to reproduce the SPFs of zodiacal dust (Hong 1985) and Saturn's rings (Hedman & Stark 2015) over a large range of scattering angles. We explore SPFs including additional complexity in Appendix A.

The intensity for a given pixel ($x', y'$) integrated along the line of sight is (Millar-Blanchaer et al. 2015)

$$ I(x', y') = I_0 + \cos i \int_{z'=-R_2}^{R_2} \frac{N_0}{r^2} \eta(r, z) p_2(g_1, \alpha, g_2, \theta) dz', \quad (3) $$

where $\eta(r, z)$ is the dust density distribution of the disk, inclined, i, from face-on and rotated so that the semimajor axis is PA degrees west of north. In this model, the disk offset is defined in the plane of the disk with $dx$ (in au) in the disk minor-axis direction toward the northwest and $dy$ (in au) in the disk major-axis direction toward the southwest. Using the disk-plane offset, we can measure the eccentricity ($e = \sqrt{dx^2 + dy^2}/R_1$) and argument of pericenter ($\omega = \arctan(dy/dx)$). Consistent with the definition, the SPF is implemented at the grain level, in the disk midplane before the image is inclined and offset. The observed azimuthal surface brightness distribution is dependent not only on the SPF but also on the inclination (creating limb brightening) and the stellar offset (more light shines on the pericenter than the apocenter, making the reflected light from the pericenter brighter than the apocenter). We assume that $N_0$ is the flux normalization and $I_0$ is a constant offset. If the ADI-KLIP algorithm is applied correctly, then the average sky background is expected to be zero (Soummer et al. 2012). As a result, we assume that $I_0 = 0$. In summary, our model has 11 free parameters, including 7 geometric parameters ($R_1, R_2, \beta, i, \text{PA}, dx$, and $dy$), 3 SPF parameters ($g_1, g_2, \alpha$), and the normalization $N_0$.

### 3.2.2. Bayesian Estimation of the Parameters

We used Bayesian parameter estimation to derive the best-fit values and the posterior distribution functions (PDFs) for our model of the HR 4796A total intensity images. Our model images (Figure 6, top left) are forward-modeled to include the instrument and PSF subtraction artifacts. This process, illustrated in Figure 6 for our GPI H-Spec total intensity image, is as follows:

1. First, we generated a model image of the disk (described in Section 3.2.1) with 11 free parameters.

2. Second, we estimated the PSF for the observing sequence being analyzed. Each GPI Spec observation included the four satellite spots. We approximated the PSF during each observation by averaging the images of the four satellite spots to increase the S/N. In some exposures, the disk rotated onto the satellite spots. In these cases, we omitted the affected satellite spot from the PSF estimate and masked the portion of the disk that overlapped with the satellite spot. The resulting PSF is shown in a small vignette in the bottom left panel of Figure 6. We convolved this PSF with the model.
GPI H band: Best Model and Residuals

<figure><img src="image_8.png" /><figcaption>Figure 6. Best-fit model resulting from the MCMC: GPI H-band (top left) best-fit model image; the red and green spots mark the position of the star and the disk pericenter, respectively. Bottom left: best-fit model image after convolution with an observed GPI PSF (shown in a small vignette at the bottom left). Top middle: best-fit model image after convolution and forward modeling to reproduce the ADI-KLIP effect. The white dashed ellipses indicate the inner and outer edges of the zone over which the likelihood for the MCMC is estimated. Bottom middle: Image showing the ADI-KLIP-reduced data set. Top left: residuals from the MCMC. Bottom right: S/N of the residuals of the MCMC.</figcaption></figure>

3. Third, we forward-modeled the observing sequence (Pueyo 2016) using the DiskFM pipeline<fnref n="30" /> of the pyKLIP package (Wang et al. 2015) to simulate the impact of the ADI-KLIP on the model. The forward model was constructed by simulating the convolved disk at each of the parallactic angles within our observing sequence and combining these images using the same pyKLIP parameters as the observations (Figure 6, top middle).

4. Fourth, we estimated the uncertainty map using a technique described in Gerard & Marois (2016). To estimate the noise map, we reduced the data using ADI-KLIP reduction with the same parameters but back-rotated to produce a final reduced image in which the disk is averaged out. We calculated the standard deviation of the pixel values within concentric annuli, centered on the star, with a constant width (3 pixels) in the combined image with the disk averaged out. We used this azimuthally symmetric image showing the standard deviation of the pixel values as in our noise map.

5. Finally, we compared the forward model to the reduced image (Figure 6, bottom middle) by measuring

$$ \chi^2 = \sum_s \frac{(\text{Data} - \text{ForwardModel})^2}{\text{Uncertainty}^2}, \quad (4) $$

where S is the zone, corresponding to the disk. S is defined as an ellipsoidal annulus with the same inclination and position angle as the HR 4796A disk and with inner and outer radii of 40 and 130 au, respectively (shown in a white dashed line in Figure 6, top middle). We applied an empirical scalar correction to the noise map to retrieve realistic error bars.

We performed these steps hundreds of thousands of times within an MCMC wrapper that maximizes $e^{-\chi^2/2}$ until the chains had converged, using the `emcee` package (Foreman-Mackey et al. 2013). We used 192 parallel walkers and removed some iterations during the burn-in phase. We noted the total number of iterations for each band at the top of each corner plot. After the MCMC converged, we plotted the residuals (`|Data - BestForwardModel|`; Figure 6, top right) and the S/N of the residuals (`|Data - BestForwardModel|`/Uncertainty; Figure 6, bottom right), as well as the posterior distribution functions. All corner plots are shown in Appendix C. We derived uncertainties based on the 16th ($-1\sigma$), 50th (median value), and 84th ($+1\sigma$) percentiles of the samples in the distributions (plotted as vertical lines in the corner plots). We demonstrated our method on a simulated data set in Section B.2. We found that the model priors could be recovered with Gaussian PDFs to within $1\sigma$ (with the exception of $R_2$) if the noise map was artificially multiplied by a scalar.

We found that the best value for $R_2$ often depends on the KLIP parameters used; higher KL modes, consistent with more

<footnote n="30">³⁰ https://pyklip.readthedocs.io/en/latest/diskfm_gpi.html</footnote>
<table><caption>**Table 2**
Debris Ring Analyzer Geometric Parameters (Sky Plane and Disk Plane)</caption><thead><tr><th></th><th>a' (mas)</th><th>e' (mas)</th><th>PA' (deg)</th><th>ΔR.A' (mas)</th><th>ΔDecl' (mas)</th><th>a (mas)</th><th>e (mas)</th><th>i (deg)</th><th>ω (deg)</th><th>Ω (deg)</th></tr></thead><tbody><tr><td colspan="11">Mask-and-Interpolate</td></tr><tr><td>J 20160323</td><td>1053 ± 6</td><td>0.970 ± 0.001</td><td>26.2 ± 0.5</td><td>6 ± 4</td><td>-15 ± 7</td><td>1053 ± 6</td><td>0.02 ± 0.01</td><td>75.9 ± 0.1</td><td>-22 ± 56</td><td>26.2 ± 0.5</td></tr><tr><td>H 20160318</td><td>1063 ± 16</td><td>0.970 ± 0.002</td><td>26.5 ± 0.4</td><td>0 ± 8</td><td>-6 ± 11</td><td>1063 ± 16</td><td>0.01 ± 0.02</td><td>76.5 ± 0.5</td><td>-61 ± 180</td><td>26.5 ± 0.4</td></tr><tr><td>K1 20140325</td><td>1079 ± 27</td><td>0.970 ± 0.002</td><td>26.5 ± 0.4</td><td>4 ± 10</td><td>-17 ± 19</td><td>1079 ± 28</td><td>0.02 ± 0.02</td><td>77.0 ± 0.6</td><td>-44 ± 62</td><td>26.5 ± 0.4</td></tr><tr><td>K1 20150402</td><td>1060 ± 12</td><td>0.972 ± 0.002</td><td>26.6 ± 0.4</td><td>-1 ± 6</td><td>-20 ± 10</td><td>1061 ± 12</td><td>0.04 ± 0.02</td><td>76.4 ± 0.5</td><td>-67 ± 19</td><td>26.6 ± 0.4</td></tr><tr><td>K1 20150403</td><td>1066 ± 22</td><td>0.973 ± 0.002</td><td>26.5 ± 0.4</td><td>8 ± 8</td><td>-19 ± 12</td><td>1066 ± 22</td><td>0.02 ± 0.01</td><td>76.5 ± 0.6</td><td>-14 ± 46</td><td>26.5 ± 0.4</td></tr><tr><td>K2 20150403</td><td>1075 ± 30</td><td>0.973 ± 0.003</td><td>26.5 ± 0.5</td><td>-6 ± 9</td><td>-13 ± 16</td><td>1076 ± 30</td><td>0.05 ± 0.02</td><td>76.7 ± 0.7</td><td>-79 ± 22</td><td>26.5 ± 0.5</td></tr><tr><td colspan="11">RDI-NMF</td></tr><tr><td>H 20160318</td><td>1069 ± 18</td><td>0.970 ± 0.002</td><td>26.5 ± 0.5</td><td>-1 ± 9</td><td>-14 ± 14</td><td>1070 ± 18</td><td>0.03 ± 0.02</td><td>76.5 ± 0.5</td><td>-68 ± 50</td><td>26.5 ± 0.5</td></tr><tr><td colspan="11">RDI-KLIP</td></tr><tr><td>H 20160318</td><td>1063 ± 15</td><td>0.972 ± 0.001</td><td>26.5 ± 0.4</td><td>0 ± 7</td><td>-15 ± 11</td><td>1063 ± 15</td><td>0.03 ± 0.01</td><td>76.3 ± 0.4</td><td>-65 ± 30</td><td>26.5 ± 0.4</td></tr><tr><td colspan="11">ADI-KLIP</td></tr><tr><td>J 20160323</td><td>1069 ± 14</td><td>0.972 ± 0.002</td><td>26.5 ± 0.5</td><td>2 ± 7</td><td>-18 ± 11</td><td>1069 ± 14</td><td>0.03 ± 0.02</td><td>76.5 ± 0.5</td><td>-58 ± 43</td><td>26.5 ± 0.5</td></tr><tr><td>H 20160318</td><td>1069 ± 15</td><td>0.972 ± 0.002</td><td>26.5 ± 0.4</td><td>0 ± 7</td><td>-15 ± 10</td><td>1069 ± 15</td><td>0.03 ± 0.02</td><td>76.5 ± 0.4</td><td>-65 ± 35</td><td>26.5 ± 0.4</td></tr><tr><td>K1 20140325</td><td>1070 ± 22</td><td>0.972 ± 0.002</td><td>26.4 ± 0.5</td><td>1 ± 9</td><td>-21 ± 16</td><td>1070 ± 22</td><td>0.04 ± 0.02</td><td>76.4 ± 0.5</td><td>-61 ± 23</td><td>26.4 ± 0.5</td></tr><tr><td>K1 20150402</td><td>1077 ± 16</td><td>0.973 ± 0.002</td><td>26.4 ± 0.4</td><td>-2 ± 8</td><td>-21 ± 13</td><td>1078 ± 16</td><td>0.05 ± 0.02</td><td>76.7 ± 0.5</td><td>-69 ± 23</td><td>26.4 ± 0.4</td></tr><tr><td>K1 20150403</td><td>1069 ± 13</td><td>0.972 ± 0.002</td><td>26.4 ± 0.5</td><td>1 ± 7</td><td>-24 ± 11</td><td>1069 ± 13</td><td>0.04 ± 0.02</td><td>76.4 ± 0.4</td><td>-62 ± 21</td><td>26.4 ± 0.5</td></tr><tr><td>K2 20150403</td><td>1074 ± 18</td><td>0.973 ± 0.002</td><td>26.5 ± 0.5</td><td>0 ± 8</td><td>-23 ± 13</td><td>1076 ± 18</td><td>0.05 ± 0.02</td><td>76.6 ± 0.5</td><td>-65 ± 24</td><td>26.5 ± 0.5</td></tr></tbody></table>

**Notes.** The uncertainties listed in this table are the statistical uncertainties from our fits and do not include the uncertainty either in the direction of true north or in the position of the star. The uncertainty in the measurement of true north is 0.03 for GPI (Konopacky et al. 2014). The uncertainty in the position of the star is 0.05 pixels (or ~0.7 mas) for GPI (Wang et al. 2014).

aggressive KLIP reduction, leads to smaller $R_2$. We hypothesized that ADI-KLIP has difficulty recovering $R_2$ because the diffuse extended structures were severely impacted by ADI self-subtraction. Indeed, the STIS coronagraphic image of this disk revealed diffuse, low surface brightness scattered light, extending hundreds of au from HR 4796A, extending well outside of both the GPI and SPHERE fields of view (Schneider et al. 2018). Therefore, we concluded that the values estimated for $R_2$ with this method are not physically relevant, a fact confirmed by our test on a simulated data set in Section B.2. As a result, we assumed a prior for $log R_2$ that was uniform from 82 to 100 au with an inverse exponential function $1/(1 + \exp(k(R - R_c))$ and $k = 40$ beyond 100 au to avoid wasting computational time on an irrelevant parameter. We assumed that all of the other priors were uniform between our minimum and maximum values. For example, we assumed uniform priors for $log R_1$ between 60 and 80 au, for $\beta$ between 1 and 30, for $g_1$ between 5% and 99.99%, for $g_2$ between -5% and -99.99%, for $\alpha$ between 1% and 99.99%, for $i$ between 70° and 80°, for $\theta$ between 20° and 30°, for $dx$ and $dy$ between -10 au and 10 au, and finally for $log N$ between 0.5 and 50,000.

## 4. Results: Disk Geometry

The dust in the HR 4796 disk is produced by collisions among parent bodies on Keplerian orbits. Thus, we idealize the narrow dusty ring using a Keplerian orbit, specified by five orbital elements: (1) the semimajor axis, $a$; (2) the eccentricity, $e$; (3) the inclination from face-on viewing, $i$; (4) the PA of the ascending node, $\Omega$; and (5) the argument of pericenter, $\omega$. We infer the orbital elements from our measurements of the projected semimajor axis, $a'$, the projected eccentricity, $e'$, the

position angle (measured east of north), $PA'$, and the projected offset of the host star from the center of the ellipse, $\Delta R.A'$ and $\Delta decl'$ (see Table 2). We adopted the same convention as M17: a positive $\Delta R.A'$ means that the ellipse center (blue cross in Figure 7) is located west of the star (red star in Figure 7), and a positive $\Delta decl'$ means that the ellipse center is located north of the star. The uncertainty in the measurement of true north is 0.03 for GPI (Konopacky et al. 2014).

### 4.1. Debris Ring Analyzer Analysis

We used the DRA code to estimate the projected parameters and then to apply the Kowalsky deprojection routine to transform the projected disk parameters to deprojected orbital parameters (Smart 1930; Stark et al. 2014). We measured the projected disk properties and inferred the orbital parameters for each of our observations using our Mask-and-Interpolate, RDI-NMF, RDI-KLIP, and ADI-KLIP PSF-subtracted images to better understand the uncertainties in our derived quantities.

We measured projected disk parameters (see Table 2) that are generally consistent from night to night and from filter to filter using all of our PSF subtraction techniques; however, we observed some minor systematic differences. In general, the standard deviation of our KLIP measurements is smaller than that of our Mask-and-Interpolate measurements. For example, the average and standard deviation of the KLIP and Mask-and-Interpolate are $\Delta R.A' = 0 \pm 2$ mas, $\Delta decl' = -20 \pm 3$ mas and $\Delta R.A' = 2 \pm 5$ mas, $\Delta decl' = -15 \pm 5$ mas, respectively. Thus, we conclude that the residual speckle noise within the Mask-and-Interpolate reductions makes it more difficult to measure the geometric parameters reliably. Despite this caveat, our measurements for the positional offset of the star $\Delta R.A'$ and $\Delta decl'$ are broadly consistent with those
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Demonstration of the projected ellipse fitting procedure for our H-band observation for each of the PSF subtraction techniques: Mask-and-Interpolate (top left), RDI-NMF (top right), ADI-KLIP (bottom left), and RDI-KLIP (bottom right). Diamonds show the positions of peak disk surface brightness calculated in 3° wedges centered on the star. The red star marks the position of the star, and the blue cross marks the center of the ellipse. We excluded wedges near forward and backward scattering because the residual speckle noise in this region was large. The size of the excluded wedges depends on the PSF subtraction technique, with Mask-and-Interpolate requiring the largest exclusion wedges and RDI-KLIP the smallest. The excluded angles WEST and EAST give the angular extent of the excluded region near forward and backward scattering, respectively.</figcaption></figure>

obtained from other observations (Schneider et al. 2009; Thalmann et al. 2011; Wahhaj et al. 2014; Rodigas et al. 2015; Perrin et al. 2015, M17).

Consistent with previous studies, we estimated a small eccentricity, suggesting that the ring is very circular. Specifically, we estimated an orbital eccentricity ($e \sim 0.03$) that is comparable to that reported from GPI commissioning observations ($e \sim 0.02$; Perrin et al. 2015), Magellan AO ($e \sim 0.06$; Rodigas et al. 2015),

and SPHERE ($e \sim 0.06$; M17). The almost circular nature of the ring makes it challenging to constrain the argument of pericenter, especially with this geometrical fit method. Comparing our measured value for $\omega$ to others published in the literature is further complicated by inconsistencies in the definition of this parameter. We adopted the same definition as M17 and Milli et al. (2019), who found $\omega = -74.3 \pm 6.2^\circ$ and $\omega = -74.2 \pm 11.9^\circ$, respectively. These results are in good agreement with
<table><caption>**Table 3**
Disk Fitting Geometry Parameters (Disk Plane)</caption><thead><tr><th></th><th><i>R</i><sub>1</sub><br>(mas)</th><th><i>e</i></th><th><i>i</i><br>(deg)</th><th>ω<br>(deg)</th><th>Ω<br>(deg)</th></tr></thead><tbody><tr><td>GPI J 20160323</td><td>1038 ± 5</td><td>0.023 ± 0.004</td><td>76.9 ± 0.1</td><td>-28<sup>+15</sup><sub>-13</sub></td><td>26.8 ± 0.1</td></tr><tr><td>SPHERE H2<sup>a</sup></td><td>1032 ± 1</td><td>0.047 ± 0.002</td><td>76.77 ± 0.03</td><td>-65 ± 1</td><td>26.94 ± 0.04</td></tr><tr><td>GPI H 20160318</td><td>1032 ± 2</td><td>0.028 ± 0.003</td><td>76.81 ± 0.07</td><td>-64 ± 4</td><td>26.7 ± 0.1</td></tr><tr><td>GPI K1 20150403</td><td>1031 ± 4</td><td>0.035<sup>+0.003</sup><sub>-0.002</sub></td><td>76.47<sup>+0.07</sup><sub>-0.08</sub></td><td>-46 ± 4</td><td>26.6 ± 0.6</td></tr><tr><td>GPI K2 20150403</td><td>1036<sup>+6</sup><sub>-7</sub></td><td>0.024 ± 0.004</td><td>76.6 ± 0.1</td><td>-42<sup>+12</sup><sub>-11</sub></td><td>26.5 ± 0.1</td></tr></tbody></table>

**Notes.** The uncertainties listed in this table are derived from the 16th, 50th, and 84th percentiles of the marginalized posterior distribution for each parameter and for each data set. They do not include the uncertainty either in the direction of true north or in the position of the star. The uncertainty in the measurement of true north is 0°03 for GPI (Konopacky et al. 2014) and 0°08 for SPHERE/IRDIS (Maire et al. 2016). The uncertainty in the position of the star is 0.05 pixels (or ~0.7 mas) for GPI (Wang et al. 2014) and 0.25 pixels (or 3 mas) for SPHERE/IRDIS (M17).

Olofsson et al. (2019), who adopted a different convention ($\omega = -254^{\circ}3 \pm 1^{\circ}8 = -180^{\circ}-74^{\circ}3 \pm 1^{\circ}8$), and with Rodigas et al. (2015), who adopted yet another convention ($\omega = 110.6^{\circ}6 \pm 13^{\circ}7 = 180^{\circ}-70.6^{\circ}6 \pm 13^{\circ}7$). Finally, Perrin et al. (2015) found $\omega = -16^{\circ}9 \pm 1^{\circ}9$. Since they did not discuss this value in detail, we could not confirm the convention they used; however, they found a relative position of the ring center (small offset to the southwest relative to the star) that would place the position of the pericenter closer to the north ansae, in slight disagreement with more recent studies. Using the DRA, we estimated $\omega$ between -14° and -80° with a typical uncertainty of a few tens of degrees. If we exclude the low-S/N mask and interpolate data set values, we are finding values between -59° and -69°, in very good agreement with the most recent analysis. The DRA is not precise enough to provide an accurate estimation of the ring center compared with the stellar position, especially for low-S/N observations, such as our J and K2 observations. For this reason, we used the modeling analysis to derive more precise geometric parameters that are described in the next section.

## 4.2. Disk Modeling Analysis

We fitted our GPI J-, H-, K1-, and K2-Spec total intensity images using the procedure described in Section 3.2. We plotted the posterior distribution functions for the parameters estimated in Figure C1 for H-Spec (see Figure B3 in Appendix B, and Appendix C for the other wavelengths). The best-fit model for the H-Spec observation and its forward model are shown in Figure 6, and the ones for the other observations are shown in Appendix C. We marked the position of the star with a red point and the position of pericenter with a green point. The S/N residuals for the individual pixels were typically less than 2; however, visual inspection of the residual S/N image revealed some disk structure indicating that the disk is not perfectly subtracted. Such residuals are not surprising because HR 4796A is observed with high S/N and the model for the dust geometry is simple with relatively few free parameters compared with the number of pixels over which the disk is resolved.

We extracted the best-fit disk geometric parameters (disk inner radius $R_1$, eccentricity $e$, inclination $i$, argument of the pericenter $\omega$, and principal angle $\Omega$) from the MCMC fit. We listed the best-fit parameters in Table 3. Since the time needed to optimize the fit for a single band was long, we decided to optimize only the fit for the highest-S/N K1-Spec data set (2015/04/03). The disk geometric parameter uncertainties were smaller but consistent with those estimated from DRA

fitting (in Section 4.1). The disk semimajor axis, $a$, estimated from DRA fitting should not be directly compared with the disk inner radius, $R_1$, estimated from model fitting because the two quantities are defined differently. The DRA semimajor axis $a$ is the position of the maximum of the disk in the convolved image, and model $R_1$ is the position of the inner radius before convolution.

We found that most of the GPI best-fit parameters are approximately consistent with one another. However, we noted that the stellar offset (eccentricity and argument of the pericenter) was slightly inconsistent because the uncertainties estimated in model fitting did not include the errors in the registration of the star and the direction of true north. Our geometric parameters are generally consistent with those previously measured by Olofsson et al. (2019) and Milli et al. (2017, 2019). Specifically, we also found that the disk pericenter was located on the front side of the disk (northwest of the star). However, we measured an eccentricity and pericenter that were slightly different than previously reported in Milli et al. (2019) ($e = 0.072 \pm 0.037$, $\omega = -74^{\circ}2 \pm 11.9$). We hypothesized that some of these differences may be the result of uncertainties in the measurement of stellar positions due to differences in the instrument and even the observing conditions. However, using the same SPHERE IRDIS H data, we also found differences between our measurements ($e = 0.047 \pm 0.002$, $\omega = -65^{\circ} \pm 1$) and those of Milli et al. (2017) ($e = 0.070 \pm 0.011$, $\omega = -72^{\circ}4 \pm 5.1$). As a result, we concluded that these small differences in the measured position of the disk center are the result of different measurement methods. Finally, in our disk modeling, we assumed that the disk was circular, while Olofsson et al. (2019) assumed that the disk was elliptical. Since the HR 4796A disk has such a small eccentricity, we do not expect that these differences in assumed eccentricity should impact the comparison of our analyses.

Multiwavelength observations indicate that the NE side of the disk is ~5%–20% brighter than the SW side of the disk (Telesco et al. 2000; Wahhaj et al. 2014; Olofsson et al. 2019). We generated a variety of disk models, varying the disk geometry and SPF, to explore the origin of this effect. Since the SPF is always symmetric with respect to the minor axis (and is studied in more detail in the next section), we did not anticipate that changes in the SPF could reproduce the observed brightness asymmetry. Instead, we expected that geometric effects would make one side of the disk appear brighter than the other if one side was closer to the star than the other. We

<footnote><sup>a</sup> Derived from a reanalysis of the data as presented in this manuscript and not from M17.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** SPFs for the HR 4796A disk obtained at *H* band (left) and *K1* band (right). The phase functions extracted using Mask-and-Interpolate (blue), RDI-NMF (pink), ADI-KLIP (yellow), and RDI-KLIP (green) are overplotted with SPHERE *H2* SPF (black; M17) for comparison. Unfortunately, the library of reference PSFs is substantially smaller in *J* and *K1* than in *H*; therefore, reference PSF subtraction for *J* and *K1* bands is currently not feasible. For *K1*, we plot the SPF extracted from the 2015 April 3 observation because it had a large number of individual frames and large field rotation under similar observing conditions and therefore the highest S/N. We use the standard deviation of the three measurements as a proxy for the uncertainty.</figcaption></figure>

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

created a hypothetical disk model with a stellar offset that produced a pronounced north-south asymmetry at the forward-scattering part of the disk in Figure B2 (bottom middle). Once this asymmetry was subtracted from the observations, we detected $2\sigma$ residuals at pericenter in the fits to the GPI *H* and GPI *K1* observations. These remaining residuals suggested that the disk pericenter may be bright not only because it is closer to the star but also because it has a higher dust density, consistent with the hypothesis of Olofsson et al. (2019).

Recent simulations of disks composed of collisionless particles show that the apocenter is expected to be $1+e$ brighter than the pericenter at far-infrared and submillimeter wavelengths because the larger number of particles at apocenter overcomes the brightness difference between the apocenter and pericenter at these wavelengths (Pan et al. 2016). However, existing ALMA observations have too low an S/N to detect the predicted apocenter glow for HR 4796 (Kennedy et al. 2018).

# 5. Results: Scattering Phase Function

## 5.1. Reconstruction Analysis

After verifying that reconstruction of the SPF from PSF-subtracted images should be possible, we reconstructed the SPF of the HR 4796A disk from our GPI total intensity images. First, we corrected the images for inverse-square illumination effects. We used our ADI-KLIP measured stellar offset to calculate the deprojected distance from the central star in pixels, *r*, and multiplied by the square of the deprojected distance from the central star (*r*²). Next, we computed the polar projection for each PSF-subtracted image by averaging the *r*²-scaled surface brightness of the disk within 3° wedge-shaped elliptical apertures at single-pixel intervals. Then, we extracted the SPF at each scattering angle from the polar projected images by summing the disk intensity over projected radii extending from ~0.′′18 to ~0.′′31. Finally, we corrected the SPF of the HR 4796 disk at the one-dimensional level using our isotropic disk model injected into an empty GPI observing sequence and PSF-subtracted with Mask-and-Interpolate, RDI-NMF, RDI-KLIP, and/or ADI-KLIP. More specifically, we

divided the HR 4796 disk SPF by that from a similarly extracted isotropic disk for each PSF subtraction technique.

We plot our corrected, empirical *H* and *K1* SPFs in Figure 8. Unfortunately, the library of reference PSFs is substantially smaller in the *J*, *K1*, and *K2* bands than in the *H* band; therefore, reference PSF subtraction (RDI) for *J*, *K1*, and *K2* is currently not feasible. We find that the *K1*-Spec Mask-and-Interpolate and ADI-KLIP SPFs are consistent with one another to within our uncertainties at scattering angles 30°–150°. In addition, we find that our GPI SPF is consistent with the SPHERE *H2* and *H3* phase functions previously reported by M17. Indeed the GPI *K1*-Spec phase function exhibits some of the same extremely forward-scattering behavior at small phase angles <45° and modest backward scattering at large phase angles >90°. However, the precision at which the phase function is measured at small (<40°) and large (>140°) phase angles is poorer in the GPI observations despite good observing conditions and a comparable range of parallactic angles. There are many factors that affect high-contrast imaging observations, including the weather conditions and the instrument performance.

Our GPI *H*-Spec Mask-and-Interpolate, RDI-NMF, RDI-KLIP, and ADI-KLIP SPFs are not consistent with one another or with the SPHERE SPF. Instead, all of our reconstructions appear to include an overcorrection for the SPF near forward and backward scattering. Since the shorter-wavelength data are more susceptible to speckle noise, one possibility is that our empty observing sequences are not sufficiently similar to our science observations to provide an accurate correction. Unfortunately, we cannot reconstruct *J* or *K2* SPFs using the PSF-subtracted images. The residual speckle noise in the *J*-band Mask-and-Interpolate image is too high. There is so little GPI data taken using the *K2* filter that we could not find a suitable empty observing sequence in which to inject our idealized disk model.

## 5.2. Modeling Analysis

For the GPI *H*-Spec image, we show the posterior distribution functions for our best-fit model, inferred from our MCMC analysis, assuming the two-component HG phase
<table><caption>**Table 4**
Disk Fitting Scattering Parameters</caption><thead><tr><th></th><th>g<sub>1</sub></th><th>g<sub>2</sub></th><th>α</th></tr></thead><tbody><tr><td>GPI J 20160323</td><td>0.80 ± 0.04</td><td>-0.21 ± 0.03</td><td>0.34 ± 0.03</td></tr><tr><td>SPHERE H2<sup>a</sup></td><td>0.90 ± 0.03</td><td>-0.17 ± 0.01</td><td>0.32<sup>+0.07</sup><sub>-0.04</sub></td></tr><tr><td>GPI H 20160318</td><td>0.79 ± 0.04</td><td>-0.21 ± 0.01</td><td>0.30<sup>+0.03</sup><sub>-0.02</sub></td></tr><tr><td>GPI K1 20150403</td><td>0.94<sup>+0.03</sup><sub>-0.04</sub></td><td>-0.15 ± 0.01</td><td>0.5<sup>+0.2</sup><sub>-0.1</sub></td></tr><tr><td>GPI K2 20150403</td><td>0.94<sup>+0.04</sup><sub>-0.06</sub></td><td>-0.15 ± 0.02</td><td>0.5 ± 0.2</td></tr></tbody></table>

Note. The uncertainties listed in this table are derived from the 16th, 50th, and 84th percentiles of the marginalized posterior distribution for each parameter and for each data set. They do not include the uncertainty either in the direction of true north or in the position of the star.

<sup>a</sup> Derived from a reanalysis of the data as presented in this manuscript and not from M17.

function model (Figure C1). We list the best-fit parameter values in Table 4. We estimate asymmetric scattering parameters, $g_1 = 0.78 \pm 0.04$ and $g_2 = -0.21 \pm 0.01$. We plot our best-fit model SPF in blue in Figure 9. This SPF is generally consistent with those inferred from previous ground- and space-based observations (Debes et al. 2008, M17). The HR 4796A dust SPF shows very forward-scattering behavior at scattering angles <40°, a minimum at scattering angles 30°–70°, and back scattering at scattering angles >70°. This shape is distinctive compared with SPFs for solar system dust and other debris disks that typically have a minimum around 90° (Hughes et al. 2018). Our MCMC-extracted SPF (blue solid line in Figure 9) shows small but significant differences in the shape of the SPF compared with that published by M17 (gray points), specifically at the position of the ansae (60° < θ < 120°). To understand this difference, we reanalyze the SPHERE H2 data with our disk fitting method.

### 5.3. Reanalysis of the SPHERE H2 Data

We applied our fitting code (already described in Section 3.2) to the already published SPHERE IRDIS H2 observations of HR 4796A (M17). These data were obtained in February 2015.<sup>31</sup> The raw frames were sky-subtracted, flat-fielded, and bad-pixel-corrected using the SPHERE data reduction and handling (DRH; Pavlov et al. 2008) pipeline and star-centered using the waffle spots. The uncertainty in the measurement of true north is 0°.08 for SPHERE/IRDIS (Maire et al. 2016). The uncertainty in the position of the star is 0.25 pixels (or 3 mas) for SPHERE/IRDIS (M17). We re-reduced the data using the *pyKLIP* ADI-KLIP algorithm and the same parameters as for the GPI data (three KL modes and an exclusion criterion of 6°), resulting in the image shown in Figure 10 (bottom middle).

One difference between the GPI Spec and SPHERE/IRDIS observations is the availability of a PSF reference. The SPHERE observations did not include satellite spots during the deep coronagraphic sequence; however, the host star was observed out of the coronagraphic mask right before and right after the ADI sequence to facilitate absolute flux calibration. We used these unocculted observations of the host star as a proxy for the PSF and convolved them with our model to estimate the observed astronomical scene. We presented the

<figure><img src="image_12.png" /><figcaption>**Figure 9.** SPHERE H2 SPF for the HR 4796A debris disk extracted by M17 from a corrected ADI-KLIP PSF-subtracted image (gray error bars) and by our MCMC disk fitting routine from SPHERE H2 (black line) and GPI H-Spec (blue line) observations. The shaded regions in the plots show the uncertainties and are estimated from 1000 randomly chosen SPFs generated by the MCMC sampler.</figcaption></figure>

original ADI-KLIP PSF-subtracted SPHERE H2 IRDIS image, the best-fit model, the best-fit model convolved with the PSF, the forward model of the observing sequence, the residuals for the observing sequence, and the S/N residuals for the observing sequence in Figure 10. We showed the posterior distribution functions for the model parameters estimated from our MCMC analysis, using the two-component HG phase function, in Figure C2.

We found that our SPHERE H2 and GPI H-Spec total intensity best-fit models were in very good agreement with one another. The geometric parameters (see Table 3) were consistent to within $1\sigma$ with the exception of the eccentricity (e). The smaller SPHERE uncertainties were the result of the higher S/N, especially at small angular separations. The uncertainties on the true north and star position, which were not folded into the reported uncertainties, were 2–3 times larger for SPHERE than for GPI. Even though the S/N was higher in the SPHERE image, the residuals in the SPHERE data were generally smaller than those in the GPI H observations.

The debris disk’s geometry and SPF were not expected to change in the intervening year between the observations, especially because the two sets of observations were made in approximately the same spectral band. From our measurements, we concluded that our 11-parameter model did not limit our ability to fit the disk because the fitting procedure produced similar results for the two different observations. We therefore concluded that the differing residuals amplitude between the IRDIS H2 and GPI H-Spec observations stemmed from differences in the fidelity of the PSFs used to simulate our observations. GPI and SPHERE observers used different methods to estimate the instrumental PSFs. GPI observers recreated the PSF from observations of the satellite spots that were spectrally collapsed. This method provides a simultaneous measurement of the PSF but produces low-S/N PSFs. SPHERE observers estimated the PSF from observations of the unocculted host star taken before and after the observational sequence. This method provides higher-S/N PSF measurements but not simultaneous. The smaller residuals in the SPHERE image indicate that, at least in the context of analysis

<footnote n="31"><sup>31</sup> Based on observations made at the Paranal Observatory under ESO program 095.C-0298(H).</footnote>
SPHERE H2 band: Best Model and Residuals

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Best-fit model resulting from the MCMC: same as Figure 6, but for SPHERE H2 band.</figcaption></figure>

of bright disks such as HR 4796, the second PSF generation method seems to produce better results.

Recent analysis suggests that the brightness asymmetry observed near pericenter cannot be entirely attributed to the fact that dust grains at pericenter receive, and therefore scatter, more light compared with those at apocenter. Indeed, Olofsson et al. (2019) required a density enhancement at pericenter to explain the brightness asymmetry observed in SPHERE observations. In our analysis, we assumed that the dust density was azimuthally symmetric. In this case, (1) any north-south brightness asymmetries in the model are the result of a stellar offset and (2) any additional brightness asymmetry beyond that predicted by our model must be the result of a local density enhancement. We did not find any large asymmetry in our residuals, suggesting that the bulk of the disk brightness asymmetry could be explained by star illumination differences. However, we did detect a marginal north-south ($2\sigma$) brightness asymmetry in the residuals (e.g., GPI H-Spec and K1-Spec), particularly at the pericenter of the disk (near forward scattering), that could be consistent with a dust density enhancement, similar to the one measured by Olofsson et al. (2019).

The fits to the SPHERE and GPI observations indicated that the SPFs were consistent to within $1\sigma$ (Figure 9), with the SPHERE SPF appearing slightly more forward scattering and slightly less backward scattering than the GPI SPF. This difference was seen in the best-fit values of $g_1$ and $g_2$ (see Table 4). For the SPHERE image, we measured asymmetric scattering parameters $g_1 = 0.90 \pm 0.03$ and $g_2 = -0.17 \pm 0.01$. Although the SPHERE and GPI measured values of $g_1$ and $g_2$ are inconsistent by >2$\sigma$, the plotted values of the SPFs are consistent to within the

uncertainties. As a result, we concluded that an SPF may be degenerate with values of $g_1$, $g_2$, and $\alpha$. Indeed, our corner plot (see Figure C2) showed a strong correlation between $g_1$ and $\alpha$. At present, no alternate prescription has been proposed that removes this degeneracy. As a result, we prefer to compare the SPFs directly and not the SPF fitting parameters $g_1$, $g_2$, and $\alpha$.

When we compared the SPF that we measured from the SPHERE IRSDIS H2 observations (black solid line in Figure 9) with that published by M17, we found significant differences in the shape of the two SPFs. This was surprising because the measurements were made from the exact same observations. When comparing an image of our best model for the SPHERE observations with that of an HR 4796A-like disk with the M17 SPF (left panel of Figure 11), we found only minor differences. Visually, the disk with the M17 SPF appeared only slightly more forward scattering. However, when we compared the residuals of these models after they were convolved with the PSF, forward-modeled, and PSF-subtracted (right panel of Figure 11), we found significant differences. Our SPF generally had much smaller residuals than the M17 SPF. For this reason, we preferred our modeling approach to measuring the SPF compared with extraction and correction from the ADI-reduced images, as described in M17 and in Section 5.1.

We attributed differences in the estimated SPFs to errors in the correction of ADI oversubtraction used by M17 and our team in Section 5.1. Previously, Milli et al. (2012) showed that ADI oversubtraction in disks was dependent on the disk geometry and surface brightness distribution. We concluded that correcting ADI oversubtraction on a complex SPF was challenging. Specifically, using an isotropic SPF model disk
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Comparison of the models (left) and residuals (right) for the SPHERE H2 data set obtained with our MCMC fit method or with the SPF extracted by Milli et al. (2017).</figcaption></figure>

that is not representative of the actual disk SPF biased the estimate of the SPF by overcorrecting the flux at the ansae. Indeed, we concluded that having a good estimate of the disk geometry and SPF is necessary in order to properly account for the complex effects of oversubtraction on disk shape and azimuthal brightness. Furthermore, we discovered that our MCMC sampler approach, using a random exploration of the parameter space with only a few assumptions for the priors, was well suited for this kind of problem.

In their analysis, M17 argued that their HR 4796A SPF was too complex to be accurately represented by a two-component HG phase function. However, in Appendix A we attempted to reproduce the M17 SPF using a more complex model (three-component HG phase function). We demonstrated that a more complex model was not necessary and that the SPF was well described using a two-component HG phase function.

## 5.4. GPI Color

We did not detect any statistically significant differences between the SPFs at *J*, *H*, *K1*, and *K2* (see Figure 12). The formal best-fit values for *g*₁, *g*₂, and *α* were statistically different; however, the curves were consistent with one another to within the uncertainties plotted. For example, the best-fit values of *g*₁ and *g*₂ indicate that the grains were more forward scattering at *K1* and *K2* (*g*₁ = 0.94 ± 0.05) than at *J* and *H* (*g*₁ = 0.80 ± 0.04) at the 3σ–4σ level. However, the SPFs overlapped to within the uncertainties over the entire range for which the SPFs were measured (~15°–155°; see Figure 12). Similar to our SPHERE H2 and GPI H-Spec analysis, we concluded that the fitting parameters did not accurately capture the consistencies and inconsistencies between the measurements and that multiple values of *g*₁, *g*₂, and *α* could be used to represent a single phase function.

The SPF uncertainties were the largest for *J* and *K2* bands and the smallest for SHERE H2 and GPI H and K1 bands. We attributed the large uncertainties in the *J*-band measurements to the large speckle noise at the shortest wavelengths, where the PSF contains the most substructure. We attributed the next largest uncertainties in *K2* band to the small throughput of the bandpass and large thermal background that limited the S/N at which the disk was detected. The best performance occurred in *H* band, where the PSF was still relatively fine but where the PSF was less complicated.

## 5.5. GPI/STIS Color

Visual scattered-light images of the HR 4796A disk and its interstellar environment have been obtained using the HST STIS coronagraph (Schneider et al. 2009, 2018).

<figure><img src="image_15.png" /><figcaption>**Figure 12.** HR 4796A SPFs estimated by fitting GPI *J* (brown), *H* (blue), *K1* (magenta), and *K2* (purple) images. Our analysis does not detect any differences between the shapes of the HR 4796A SPF as a function of near-infrared wavelength within the GPI data. The best values for all of the posteriors are given in Table 4.</figcaption></figure>

We compared the STIS data with the SPHERE near-infrared scattered-light images to estimate the near-infrared to visual color. For HR 4796A, the range of scattering angles over which the SPF could be measured was dependent on the instrument Inner Working Angle (IWA); the sampling of the SPF was dependent on the plate scale. HR 4796A was imaged using STIS twice: the first time using Wedge A at a location where the wedge is 0.63° wide (Schneider et al. 2009), and a second time using BAR5 that has a width 0.3° in conjunction with Wedge A at a location where the wedge is 1.0° wide to enable imaging near and far from the star at high dynamic range (Schneider et al. 2018). The STIS occultors were wider than the SPHERE Lyot N_ALC_YJH_S coronagraph mask (185 mas) used by M17. Unfortunately, STIS's occulter did not provide access to the disk's minor axis. In addition, the STIS plate scale (50.77 mas) was larger than the SPHERE IRDIS plate scale (12.25 mas). The SPHERE PSF was supersampled for an 8 m telescope operating at near-infrared wavelengths, while the STIS PSF was subsampled for a 2.5 m telescope operating at visual wavelengths.

We measured the HR 4796A visual SPF from the more recent STIS BAR5 and Wedge A-1.0 observations that detect the dust ring at very high S/N. Schneider et al. (2018) observed the disk at three roll angles during two epochs. During each epoch, they interleaved PSF observations between successive HR 4796A observations to provide the most representative
<figure><img src="image_16.png" /><figcaption>**Figure 13.** HR 4796A STIS, visual, total intensity SPF, extracted from high-S/N Wedge A-1.0 and BAR5 observations (Schneider et al. 2018). The 0.15 IWA does not permit measurement of the disk surface brightness along the minor axis.</figcaption></figure>

PSF. They performed reference PSF subtraction for each exposure, masking the diffraction spikes, and combining the images to provide a final PSF-subtracted image that is minimally impacted by diffraction artifacts. Their approach to PSF subtraction has been highly effective for STIS coronagraphic observations because HST and its instruments are extremely stable. We were able to extract the SPF directly from this combined, PSF-subtracted image without extensive modeling because Schneider et al. (2018) did not use ADI to construct their final STIS PSF-subtracted image of HR 4796A.

We extracted the STIS SPF for scattering angles 27°–153°, only a slightly smaller range than the 13.6–166.6° range accessible from our vantage point (Figure 13). We detected a hint of forward scattering at the smallest scattering angle and a predominantly flat SPF at larger scattering angles.

We calculated the ratio of the near-infrared to visual SPFs to search more carefully for trends in color as a function of scattering angle. Even though our measurement of the *H*-band SPF provided an improvement compared with previous measurements, we performed our analysis on the previously published SPHERE *H2* SPF derived by M17 because their analysis did not assume any functional form for the SPF. Further, to facilitate analysis of the two data sets, we interpolated the SPHERE *H2* measurements to the same scattering angles as the STIS measurements because the SPHERE measurements are well sampled and smooth. Laboratory measurements of large particles (>200 µm) suggested that the ratio of the reflectances measured at long wavelengths compared to short wavelengths should rise as a function of scattering angle, consistent with the prediction from geometric optics (Schröder et al. 2014). We plotted the ratio of the SPHERE *H2* and HST STIS SPFs in Figure 14. We found that the ratio of the SPHERE *H2* SPF with that of the HST STIS SPF has a slightly blue color at relatively small phase angles (<60°) and a slightly red color at relatively large phase angles (>135°), consistent with expectations based on measurements made in the laboratory (Schröder et al. 2014).

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Ratio of the VLT/SPHERE *H2* and the HST/STIS SPFs, revealing a slightly blue color for scattering at relatively small phase angles (<60°) and a slightly red color at relatively large phase angles (>120°). Overplotted in red and blue is the expected near-infrared to visual color predicted by our DHS model.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Observed VLT/SPHERE VBB polarized SPF in Milli et al. (2019) and the prediction by our DHS model. The DHS model is not able to reproduce the observed plateau at ~30° to ~90°.</figcaption></figure>

## 6. Dust Grain Properties

The SPF can be used to place detailed constraints on the properties of the dust. M17 modeled their SPHERE HR 4796A SPF using a grid of 7800 dust compositions and sizes assuming Mie theory and/or the Distribution of Hollow Spheres (DHS).
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Comparison of the SPF for our best-fitting DHS model (black line) and that extracted from the SPHERE H2 observations by M17.</figcaption></figure>

DHS is a theoretical construct that is used to approximate the behavior of aggregates (Min et al. 2005, 2016). M17 found that the local brightness enhancement on the back side of the disk (close to the ansa) could only be reproduced by DHS models with large silicate grains. Unfortunately, they were unable to reproduce the SPF over the full range of measured phase angles (13°–6–166°).6 using either Mie or DHS. Their smallest reduced $\chi^2 = 131$. They constructed a best-fit model to the forward-scattering part of the SPF (≤45°) using Mie theory assuming porous water ice grains (the fraction of void occupied by the ice, $p_{H_2O} = 90\%$, and the porosity without ice, $P = 0.10\%$) with a very steep power-law size distribution, $dn/ds \propto s^{-5.5}$, and a minimum grain size $s_{min} = 17.8 \mu m$. This model did not provide any back scattering. However, they noted that large grains modeled using Hapke theory reproduced not only the observed back-scattering behavior but also the forward-scattering behavior at scattering angles >30°, if they have a size of ~30 µm.

## 6.1. Modeling the SPHERE H2 Phase Function

We used MCFOST to model the SPHERE SPF, predicting the SPFs expected from various dust populations, assuming that the grains were well represented by the DHS model. We chose to model the SPHERE SPF so that the results from our modeling effort could be more directly compared with previous modeling efforts.

First, we experimented with fitting by eye. Building on the experience from M17, we assumed a large minimum grain size (~30 µm) to reproduce the observed back scattering from the more distant side of the disk. This requirement was consistent with the large minimum grain size expected from radiation pressure blowout (~10 µm). We found that we could reproduce the overall shape of the SPF by balancing the grain-size distribution power law with the minimum grain size. However, we could not reproduce the observed dip in the SPF at 30°–40° if the grains were composed of silicates or water ice only,

<table><caption>**Table 5**
MCFOST Dust Model</caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">Symbol</th><th colspan="2">Value</th></tr><tr><th>Prior<sup>a</sup></th><th>Posterior<sup>b</sup></th></tr></thead><tbody><tr><td>Scattering theory</td><td></td><td colspan="2">DHS</td></tr><tr><td>Maximum void fraction</td><td>V<sub>max</sub></td><td>(0%, 100%)</td><td>76<sup>+6</sup><sub>-4</sub>%</td></tr><tr><td>Porosity</td><td>P</td><td>(0%, 100%)</td><td>15<sup>+12</sup><sub>-11</sub>%</td></tr><tr><td>Minimum grain size</td><td>s<sub>min</sub></td><td>(0.5, 100) µm</td><td>25<sup>+4</sup><sub>-4</sub> µm</td></tr><tr><td>Maximum grain size</td><td>s<sub>max</sub></td><td>1000 µm</td><td></td></tr><tr><td>Power-law size distribution</td><td>ν</td><td>(-6, -3)</td><td>-3.74<sup>+0.12</sup><sub>-0.05</sub></td></tr><tr><td>Amorphous silicate volume<sup>c</sup></td><td>f<sub>Si</sub></td><td>(0%, 100%)</td><td>42<sup>+11</sup><sub>-13</sub>%</td></tr><tr><td>Amorphous carbon volume<sup>c</sup></td><td>f<sub>C</sub></td><td>(0%, 100%)</td><td>17<sup>+16</sup><sub>-12</sub>%</td></tr><tr><td>Metallic iron volume<sup>c</sup></td><td>f<sub>Fe</sub></td><td>(0%, 100%)</td><td>37<sup>+15</sup><sub>-9</sub>%</td></tr><tr><td>Reduced chi-squared value</td><td>χ<sub>ν</sub><sup>2</sup></td><td></td><td>2.3</td></tr></tbody></table>

**Notes.**

<sup>a</sup> The parameters are limited to three decimal digits, with uniform sampling in the prior range (except $s_{min}$, which is log-uniformly sampled).

<sup>b</sup> 16th, 50th, and 84th percentiles.

consistent with the findings of M17. We discovered that we could reproduce the dip in the SPF by adding amorphous carbon and metallic iron, common materials observed in the ISM and protoplanetary disks. As a result, we used amorphous silicate, amorphous carbon, and metallic iron in our grain models. This dust composition is consistent with those used in other studies that reproduced the color of the scattered-light and/or thermal emission from the disk (Debes et al. 2008; Rodigas et al. 2015).

After we found a reasonable fit by eye, we used a Bayesian analysis to locate the best fit and measure the posterior distribution functions. Specifically, we used the DebrisDiskFM (Ren et al. 2019) and emcee (Foreman-Mackey et al. 2013) packages to distribute the MCMC posterior calculation on different computer nodes. We allowed the maximum void fraction ($V_{max}$), porosity ($P$), amorphous silicate ($f_{Si}$), amorphous carbon ($f_{C}$), and metallic iron ($f_{Fe}$) volumes to vary between 0% and 100% with the constraint that $f_{Si} + f_{C} + f_{Fe} = 1$. Thus, only $f_{Si}$ and $f_{C}$ are explicitly sampled. We allowed the minimum grain size ($s_{min}$) to vary between 0.5 and 100 µm but fixed the maximum grain size ($s_{max} = 1000 \mu m$). We allowed the exponent for the power-law size distribution to vary between -6 and -3. The best-fitting parameters are listed in Table 5, and the resulting SPF is compared to the observed SPHERE SPF in Figure 16.

Our best-fit model is composed of 42% ± 4% amorphous silicate, 17+11/-13% amorphous carbon, and 37+15/-9% metallic iron by volume, suggesting that a substantial fraction of metallic iron was needed to reproduce the shape of the SPHERE H2 SPF. This composition has somewhat less silicates and more iron than estimated for protoplanetary disks (75% amorphous silicate, 10%–15% amorphous carbon, and 10%–15% iron sulfide; Min et al. 2016). However, we note that the amorphous carbon volume fraction posterior distribution function was not Gaussian, indicating that the porosity and the amount of amorphous carbon were not well constrained. In addition, our best-fit model requires a large minimum grain size, $s_{min} = 25 \pm 4 \mu m$, consistent with previous analyses of the SPHERE H2 SPF (M17) and somewhat larger than expected
from radiation pressure blowout if the grains are spheres. We find
that the minimum grain size ($s_{min}$), maximum void fraction ($v_{max}$),
and grain-size distribution power law ($\nu$) are correlated with
smaller minimum grains requiring grains with larger void fractions
and shallower grain-size distributions. The large minimum grain-
size estimates are consistent with new estimates for the minimum
grain size using the Discrete Dipole Approximation to more
accurately model irregular grains around luminous stars. Arnold
et al. (2019) found that the minimum grain sizes around A-type
stars are expected to be larger than previously predicted using Mie
theory. We expect that the difference between our measured GPI
SPF and our model SPF is small relative to the systematic
uncertainties in the dust models.

Once we located the best-fit for the SPHERE *H2* SPF, we
predicted the total and polarized intensity SPFs expected for the
same grain distribution in the GPI *J*, *H*, *K1*, and *K2* and the
STIS visual filter. We compared our predictions to the observed
(1) SPHERE/IRSDIS *H2* to STIS color (see Figure 14) and (2)
SPHERE/ZIMPOL visual polarized intensity phase function
(see Figure 15). We overlaid our predictions for the *H2* to STIS
color on the measurements in Figure 14. We found that our
model predictions are consistent with the observations to within
$1\sigma$. SPHERE/ZIMPOL polarized intensity images of HR
4796A in the very broad band indicate that the polarized
intensity phase function was relatively uniform over the
scattering angles observed (Milli et al. 2019). Our model
indicated that the HR 4796A dust was forward scattering at
scattering angles smaller than 13.6, the minimum scattering
angle accessible owing to the disk’s 76° inclination.

7. Discussion

In our analysis of the GPI HR 4796A observations, we (1) estimated the geometric parameters for the disk, (2) estimated the SPF for the grains, and (3) constrained the dust properties (e.g., size, porosity, and composition) assuming that the grains could be well described using a DHS. We tried two approaches for estimating the geometric parameters and the SPF: reconstruction from PSF-subtracted images (where we tried multiple PSF subtraction techniques), and recovery by forward-modeling the images, assuming a simple geometric model and a two-component Henyey–Greenstein SPF.

Our measurements of the geometric parameters estimated from
PSF-subtracted images and from forward-modeling are consistent
with one another, despite very different approaches. Our results
give us confidence that many measurement techniques are robust.
Indeed, our results are also broadly consistent with previously
published studies, albeit with significantly lower eccentricity. Our
results may also suggest that measuring the geometric parameters
for low-eccentricity disks may be disproportionately impacted by
uncertainties in the stellar position and/or PSF subtraction effects,
such as self-subtraction. We confirm the nonzero eccentricity and
therefore confirm the likelihood of an underlying planetary-mass
companion. We find tentative evidence of excess brightness at
periastron at the 2σ level, consistent with a density enhancement
there.

Consistent with M17, we measure strong forward scattering,
mild back scattering, and a minimum in the SPF at 50°–60°
scattering angle. However, the details of the SPF depend
sensitively on the technique used. For example, the SPF bump
at 90° reported by M17 is recovered in our Mask-and-Interpolate
SPF reconstruction (in K1-Spec) but not in our forward-modeling
analysis (of J-, H-, K1-Spec). The discrepancy between the two

techniques is unexpected because our model tests indicate that
both techniques can be used to robustly retrieve an input SPF.
However, our model test was carried out assuming extremely
stable observing conditions in which the PSF for the target and the
reference were the same. We did not detect any differences
between the SPFs measured at *J*, *H*, and *K1*; however, we did
detect differences between the near-infrared SPF and the visual
SPF measured from HST STIS observations. The near-infrared to
visual color of the SPF is consistent with measurements of very
large particles.

We modeled the HR 4796A SPF using MCFOST. Previous
studies struggled to fit the SPF using Mie theory and DHS. We
chose to focus our modeling efforts on DHS grains because they
better approximate aggregate grains that have been lofted off of
solar system minor bodies. Our best-fit MCFOST DHS model to
the SPHERE *H2* SPF broadly predicts (1) the visual SPF
measured from HST STIS observations (without using the HST
data in the model) and (2) the visual polarized intensity phase
function measured from SPHERE ZIMPOL observations,
although it does not predict the exact behavior in the ZIMPOL
data near 50° scattering angle correctly. Therefore, our analysis is
consistent with the presence of large aggregates composed of
smaller monomers. Unfortunately, modeling with true aggregate
scattering properties is beyond current computational capabilities.

In addition, our model predicts changes in the total intensity
and polarized intensity SPFs when moving from short to long
wavelengths within the near-infrared bands. For example, both
the total intensity and polarized intensity SPFs are expected to
have increasingly broader forward-scattering cones at longer
and longer wavelengths and smaller and smaller shoulders of
approximately constant polarized intensity at shorter and
shorter wavelengths (see Figure 17). The lack of detection of
these differences within the *J*, *H*, and *K1* SPFs is probably the
result of the low S/N with which we extract the SPFs. A new
analysis of the HR 4796A *K1*-Spec data using NMF with data
imputation, where the disk region is first ignored and then
empirically recovered using PSF-only signals to minimize
reduction bias, indicates that the observations did capture
changes in the SPF as a function of wavelength within the IFS
data cube (Ren et al. 2020).

Although not discussed here, GPI’s Polarimetric mode can
measure the polarized intensity at *J*, *H*, and *K1* bands to place
additional constraints on the grain properties. This mode uses a
Wollaston prism to simultaneously divide the light into the
Ordinary and Extraordinary beams and differencing the Ordinary
and Extraordinary images to obtain the PSF-subtracted polarized
light image. Thus, the polarized intensity image better preserves
the disk geometry and surface brightness distribution and achieves
a contrast performance ~10 times better than ADI (Perrin et al.
2015). This improved performance may allow the detection of the
trends in polarized intensity predicted by our model. GPI’s
Multiband Polarimetric images were recently analyzed in a
separate paper (Arriaga et al. 2020).

8. Conclusions

We have obtained Gemini Planet Imager *J*, *H*, *K1*-, and *K2*-
Spec observations of the iconic HR 4796A debris disk. We find
the following:

1. The *H*-band SPF is forward scattering at small scattering angles with a minimum near 60° and backward scattering at larger scattering angles, reaching a plateau at ~150°,
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Top left: DHS model total intensity SPFs shown for HST STIS and GPI J, H, K1, and K2. Top right: DHS model polarized intensity SPFs. Bottom left: ratio of the HST STIS and GPI J, H, and K1 SPFs. Bottom right: polarization fraction ($p = P/I$) SPFs.</figcaption></figure>

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

indicating that it is slightly less forward scattering than previously reported (M17).

2. The J, H, K1, and K2 SPFs are consistent with one another to within the uncertainties of the observations, suggesting there is no evidence for changes in the SPF as a function of wavelength at near-infrared wavelengths.

3. The SPF is well modeled assuming a smooth, offsetted, axisymmetric ring, indicating that there are no large asymmetries within the dust distribution. However, a marginal excess in the residuals at pericenter could be consistent with a pericenter dust enhancement recently reported by Olofsson et al. (2019).

4. The SPHERE H2 total intensity SPF is well modeled using large, irregular grains composed of amorphous silicate, amorphous carbon, and metallic iron. The large minimum grain size is consistent with expectations for irregular grains based on Discrete Dipole Approximation calculations.

5. This dust grain model also reproduces the ratio of the SPHERE H2 to STIS SPFs and the shape of the visual polarized intensity phase function, reinforcing the conclusion that the disk is composed of large, irregular grains.

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 Minis- tério da Ciência, Tecnologia e Inovação (Brazil). This material is based on work supported by the National Science Foundation under Astronomy and Astrophysics grant No. 1616097 (J.M., B.R.) and No. 1518332 (T.M.E., R.J.D.R., J.R.G., P.K., G.D.) and NASA grants NNX15AC89G and NNX15AD95G/NExSS (T.M.E., R.J.D.R., G.D., J.J.W., P.K.). J.M. acknowledges that support for part of this work was provided by NASA through the NASA Hubble Fellowship grant No. HST-HF2-51414 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. This work benefitted from NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate. This research has made use of the VizieR catalog access tool, CDS, Strasbourg, France (doi: 10.26093/cds/vizier). The original description of the VizieR service was published in A&AS 143, 23 (Ochsenbein et al. 2000). This research project (or part of this research project) was conducted using computational resources (and/or scientific computing services) at the Maryland Advanced Research Computing Center (MARCC). This paper is
dedicated to UCLA Professor Michael Jura, who not only discovered the very bright infrared excess associated with HR 4796A but also taught a generation of young astronomers how to think about debris disks.

**Software:** Gemini Planet Imager Data Pipeline (Perrin et al. 2014, 2016, http://ascl.net/1411.018), pyKLIP (Wang et al. 2015, http://ascl.net/1506.001), numpy, scipy, Astropy (Astropy Collaboration et al. 2018), matplotlib (Hunter 2007), emcee (Foreman-Mackey et al. 2013, http://ascl.net/1303.002), corner (Foreman-Mackey 2016, http://ascl.net/1702.002), DebrisDiskFM (Ren et al. 2019).

# Appendix A

## Two-component Henyey-Greenstein versus Three-component Henyey-Greenstein Scattering Phase Function

Our analysis clearly recovered a different SPF than the one published in M17 (see Figure 9). We were concerned that the main reason for this difference was that our model was too simple. The SPF extracted by M17 directly from the image may have been too complex to be reproduced by a two-component Henyey-Greenstein phase function. In this appendix we show that increasing the level of complexity of the SPF does not improve the fit and therefore that a two-component Henyey-Greenstein phase function is sufficient to describe the SPF in these data. Specifically, we introduced a third HG component to our SPF and compare to previous results. Equation (2) now becomes

$$p_3(g_1, \alpha_1, g_2, \alpha_2, g_3, \theta) = \alpha_1 \text{HG}(g_1, \theta) + \alpha_2 \text{HG}(g_2, \theta) + (1 - \alpha_1 - \alpha_2) \text{HG}(g_3, \theta). \quad (A1)$$

Using the SPF extracted by M17, we fitted a three-component Henyey-Greenstein function (green solid line in Figure A1). This shows that except for a few points, the general shape extracted by M17 can be reproduced by the five parameters ($g_1, g_2, g_3, \alpha_1, \alpha_2$) of a three-component Henyey-Greenstein function. We then ran the fitting extraction method on the SPHERE H2 data using the same model as previously but with a three-component Henyey-Greenstein function (13-parameter model instead of 11-parameter model). We chose as initial points the SPF parameters that best fitted the SPF extracted by M17 (green solid line) to be sure that the MCMC would consider this part of the parameter space.

The resulting best-fit three-component Henyey-Greenstein SPF is shown in Figure (A1; red dashed line). The black solid line is the two-component HG function extracted using this same method (see Section 5.3). First, this shows that the best fits from the two-component HG function MCMC and the three-component HG function MCMC almost totally overlap. The only difference is at very small scattering angle (<18°), where the three-HG-component SPF is very slightly more forward scattering than the two-HG-component SPF. However, considering the error bars in this region of the image, we considered that for this disk and for the range of scattering angle probed with this observation (13°5 to 166°5), a two-HG component (three parameters) can produce a function complex enough to reproduce the actual SPF. For this reason, we only used a two-HG component in this study for all the data. Second, the three-component HG function fitted to the SPF extracted by M17 (green solid line) is within the parameter range accessible (i.e., within the prior) of the MCMC sampler. However, once again, it is significantly different than the one

<figure><img src="image_22.png" /><figcaption>**Figure A1.** SPHERE H2 SPF: comparison of extraction methods. The gray points and error bars are the SPF previously extracted by M17 from a cADI PSF-subtracted image and then corrected by dividing by a forward model. We fit a three-component HG function to this SPF (green solid line). We overplot SPFs obtained by the MCMC fitting extraction method on the SPHERE H2 data with a three-component HG (red dashed line) and a two-component HG (black solid line). The shaded regions are plotted by measuring the SPF of 1000 randomly chosen accepted SPFs from the MCMC sampler.</figcaption></figure>

favored by the sampler a posteriori, which is why we favored ours for the theoretical analysis.

# Appendix B

## Scattering Phase Function: Extraction of Known Parameters

Accurately measuring the SPF is very challenging. In this paper, we used two different methods to estimate the SPF: a direct extraction from the image and an image-based forward-modeling fitting in an MCMC framework. To ensure that these techniques provide correct results on the SPF, we performed a series of tests by injecting a model disk in an empty GPI data set. This appendix shows the results of those tests.

### B.1. Validating the Reconstruction Analysis

We demonstrate that reconstruction of the HR 4796A SPF should be feasible by recovering the SPF from an idealized HR 4796A-like disk with strong forward scattering from a disk-injected, empty, observing sequence. We use MCFOST to simulate a disk with spherical grains composed of 70% silicates and 30% metallic iron by volume, a grain-size distribution, $dn/da \propto a^{-4}$, a minimum grain size $a_{min} = 5$ µm, and a maximum grain size, $a_{max} = 1$ mm. We inject the HR 4796-like disk into an empty GPI observing sequence using the same process described for the isotropic disk. As with the isotropic disk, we convolve the MCFOST disk image with an unblocked, Spec observation of HD 118335 and inject our convolved model into an empty GPI sequence. We PSF-subtract the injected, forward-scattering disk sequence and reconstruct the SPF from the PSF-subtracted image. Finally, we divide the SPF extracted from the injected, forward-scattering, disk sequence by that extracted from the injected, isotropic, disk sequence. We find that we can reconstruct the Mie disk SPFs in H and K1 bands (see Figure B1) at scattering angles $40° < \theta < 140°$. Our process fails to accurately recover the SPF at smaller and larger scattering angles, where the disk has a smaller angular offset from the coronagraphic spot and the image
<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>**Figure B1.** SPFs reconstructed from a forward-scattering, HR 4796A-like disk at *H* band (left) and *K*I band (right). The phase functions extracted using Mask-and-Interpolate (blue), RDI-NMF (pink), ADI-KLIP (yellow), and RDI-KLIP (green) are overplotted with the assumed SPF (black) for comparison. Unfortunately, the library of reference PSFs is substantially smaller in *J* and *K*I than in *H*; therefore, reference PSF subtraction for *J* and *K*I bands is currently not feasible.</figcaption></figure>

## GPI *H* band (injected model): Best Model and Residuals

<figure><img src="image_25.png" /><figcaption>**Figure B2.** Best-fit model resulting of the MCMC: same as Figure 6, but for an HD 48524 *H*-band data set with an injected HR 4796A-like disk with known parameters.</figcaption></figure>

suffers from larger speckle noise. Thus, we anticipate that reconstruction of the HR 4796A SPF at scattering angles $40^{\circ} < \theta < 140^{\circ}$ should be feasible.

### B.2. Validating the Modeling Analysis

We constructed an image assuming a disk with an inner radius $R_1 = 75$ au, an outer radius $R_2 = 100$ au, a power-law surface distribution $\beta = 12.4$, asymmetric scattering parameters $g_1 = 0.825$ and $g_2 = 0.201$, a relative weight $\alpha = 0.298$, an inclination $i = 76.8^\circ$, a position angle PA = 26.64$^\circ$, stellar positional offsets $dx = -2.0$ au and $dy = 0.94$ au, and a flux normalization $N = 80$

(hereafter “injected parameters”). Next, we injected the model into an empty GPI *H*-Spec observational sequence, with the same parallactic angles as in our HR 4796A data set. This sequence was obtained with GPI for the star HD 48525, on 2018 January 28. Finally, we reduced this sequence using the same ADI-KLIP method as for the GPI HR 4796A observations, producing the image shown in Figure B2 (bottom middle), and applied our MCMC method to extract parameters from this reduction (hereafter “recovered parameters”).

Our goal for this test was to show that disk modeling recovers the injected parameters to within $1\sigma$ of the values
<figure><img src="image_26.png" /><figcaption>**Figure B3.** MCMC posterior distributions recovered for an HR 4796A-like disk injected in an HD 48525 data set with a two-component HG phase function. The diagonal histograms show the posterior distributions of each parameter marginalized over all other parameters. In each plot, the dashed lines show the 16th, 50th, and 84th percentiles. The off-diagonal plots display the joint probability distributions with contour levels at the same percentiles. For each posterior, the parameter “true” value of the injected disk is overplotted in green.</figcaption></figure>

assumed. However, during initial testing, we discovered that the MCMC error bars were too small (the injected parameters were originally within $2\sigma$ or $3\sigma$ of the recovered parameters), suggesting that we underestimated the noise in our uncertainty maps (probably because we made an axisymmetric noise assumption). Therefore, we multiplied our noise map by a scalar factor of 4 to obtain PDFs with a Gaussian shape and parameters well within $1\sigma$ of the injected parameters. We estimated the scalar factor with which to multiply the noise by requiring that the amplitude of the noise in our residual maps has $S/N < 2$. For our data, we typically multiplied our uncertainty maps by factors of 3–5. Scalar multiplication of

the noise map does not change the maximum value of the likelihood but does produce larger error bars.

In this initial test, we recovered the values for all of the injected parameters (except $R_2$) to within $1\sigma$ (see Figure B3; injected parameter values are plotted as green lines), with a noise map multiplied by a scalar factor of 4. Only the value recovered for the outer radius ($R_2 = 96 \pm 1$ au) was $>4\sigma$ different from that assumed in the model, which was expected (see previous section). We perfectly recovered the injected SPF. In the left panel of Figure B4, we show the injected (red dashed line) and recovered SPFs (green solid line shows the SPF, and the shaded green area shows its associated
<figure><img src="image_27.png" /><figcaption>**Figure B4.** SPF injected into disk image (red dashed line) compared with the estimated phase function recovered from disk fitting (green solid line). The shaded regions in the plots show the uncertainties and are estimated from 1000 randomly chosen SPFs generated by the MCMC sampler.</figcaption></figure>

uncertainty). For this PSF plot and all SPF plots hereafter, the shaded regions represent the region covered by 1000 randomly chosen SPFs from the MCMC sampler, after convergence of the MCMC.

In this simple test, the model science images were generated using the same code that was used to simulate the disk in our retrievals, and the disk was completely removed in the residuals images. However, our simple 11-parameter disk model did not reproduce the complex structure observed at high S/N in our HR 4796A images. We attributed the challenges in our disk modeling to two main limitations. First, we assumed that the dust density was azimuthally symmetric within the disk. However, the residual maps showed significant local differences in the surface brightness along the major axis that could not be explained by the SPF because the SPF is symmetric with respect to the minor axis. Second, our satellite spot PSF had relatively low S/N. To build S/N, we averaged all of the satellite spot observations within an observational sequence. As a result, we did take into account changes in observing conditions on fast timescales.

## Appendix C: MCMC Detailed Results

This appendix shows the products of the MCMC for all data sets used in this paper. We show the corner plots for GPI J, H, K1, and K2 bands and SPHERE H2 band and the best model for GPI J, K1, and K2 bands (see Figures C1–C8).
GPI H band: 6,547 iterations (with 500 burn-in)

192 walkers: 1,257,024 models

<figure><img src="image_28.png" /><figcaption>Figure C1. Same as Figure B3, but for ADI-KLIP GPI H band of the HR 4796A debris disk.</figcaption></figure>
SPHERE H2 band: 3,325 iterations (with 2,000 burn-in)

192 walkers: 638,400 models

<figure><img src="image_29.png" /><figcaption>Figure C2. Same as Figure B3, but for ADI-KLIP SPHERE H2 band of the HR 4796A debris disk.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>GPI J band: Best Model and Residuals</figcaption></figure>

**Figure C3.** Best-fit model resulting from the MCMC: same as Figure 6, but for GPI J band.
GPI J band: 4,950 iterations (with 500 burn-in)

192 walkers: 950,400 models

<figure><img src="image_31.png" /><figcaption>**Figure C4.** Same as Figure B3, but for ADI-KLIP GPI J band of the HR 4796A debris disk.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>GPI K1 band: Best Model and Residuals</figcaption></figure>

**Figure C5.** Best-fit model resulting from the MCMC: same as Figure 6, but for GPI K1 band.
GPI K1 band: 5,775 iterations (with 1,500 burn-in)

192 walkers: 1,108,800 models

<figure><img src="image_33.png" /><figcaption>Figure C6. Same as Figure B3, but for ADI-KLIP GPI K1 band of the HR 4796A debris disk.</figcaption></figure>
GPI K2 band: Best Model and Residuals

<figure><img src="image_34.png" /><figcaption>**Figure C7.** Best-fit model resulting from the MCMC: same as Figure 6, but for GPI K2 band.</figcaption></figure>
GPI K2 band: 7,734 iterations (with 1,800 burn-in)

192 walkers: 1,484,928 models

<figure><img src="image_35.png" /><figcaption>Figure C8. Same as Figure B3, but for ADI-KLIP GPI K2 band of the HR 4796A debris disk.</figcaption></figure>

ORCID iDs

Christine Chen https://orcid.org/0000-0002-8382-0447

Johan Mazoyer https://orcid.org/0000-0002-9133-3091

Charles A. Poteet https://orcid.org/0000-0003-4845-7483

Bin Ren https://orcid.org/0000-0003-1698-9696

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

Justin Hom https://orcid.org/0000-0001-9994-2142

Pauline Arriaga https://orcid.org/0000-0001-6364-2834

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

Jessica Arnold https://orcid.org/0000-0001-7824-5372

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

Juan Sebastián Bruzzone https://orcid.org/0000-0002-2731-0397

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

Élodie Choquet https://orcid.org/0000-0002-9173-0740

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

Zachary H. Draper https://orcid.org/0000-0002-1834-3496

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

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

Dean C. Hines https://orcid.org/0000-0003-4653-6161

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

Franck Marchis https://orcid.org/0000-0001-7016-7277

Julien Milli https://orcid.org/0000-0001-9325-2511

Marshall D. 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

Timothy J. Rodigas https://orcid.org/0000-0002-7535-2997

Gael M. Roudier https://orcid.org/0000-0002-7402-7797

Glenn Schneider https://orcid.org/0000-0002-4511-5966
Rémi Soummer https://orcid.org/0000-0003-2753-2819

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

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

Alycia J. Weinberger https://orcid.org/0000-0001-6654-7859

David J. Wilner https://orcid.org/0000-0003-1526-7587

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

## References

Arnold, J. A., Weinberger, A. J., Videen, G., & Zubko, E. S. 2019, AJ, 157, 157

Arriaga, P., Fitzgerald, M. P., Duchêne, G., et al. 2020, arXiv:2006.06818

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155

Chen, C. H., & Kamp, I. 2004, ApJ, 602, 985

Debes, J. H., Weinberger, A. J., & Schneider, G. 2008, ApJL, 673, L191

Draper, Z. H., Duchêne, G., Millar-Blanchaer, M. A., et al. 2016, ApJ, 826, 147

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

Foreman-Mackey, D. 2016, JOSS, 1, 24

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

Gerard, B. L., & Marois, C. 2016, Proc. SPIE, 9909, 990958

Hedman, M. M., & Stark, C. C. 2015, ApJ, 811, 67

Henyey, L. G., & Greenstein, J. L. 1941, ApJ, 93, 70

Hong, S. S. 1985, A&A, 146, 67

Hughes, A. M., Duchêne, G., & Matthews, B. C. 2018, ARA&A, 56, 541

Hunter, J. D. 2007, CSE, 9, 90

Jura, M., Ghez, A. M., White, R. J., et al. 1995, ApJ, 445, 451

Jura, M., Zuckerman, B., Becklin, E. E., & Smith, R. C. 1993, ApJL, 418, L37

Kalas, P. G., Rajan, A., Wang, J. J., et al. 2015, ApJ, 814, 32

Kennedy, G. M., Marino, S., Matrà, L., et al. 2018, MNRAS, 475, 4924

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

Lafrenière, D., Marois, C., Doyon, R., et al. 2007, ApJ, 660, 770

Lafrenière, D., Marois, C., Doyon, R., & Barman, T. 2009, ApJL, 694, L148

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

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

Maire, A.-L., Langlois, M., Dohlen, K., et al. 2016, Proc. SPIE, 9908, 990834

Marois, C., Lafrenière, D., Doyon, R., et al. 2006, ApJ, 641, 556

Marois, C., Zuckerman, B., Konopacky, Q. M., et al. 2010, Natur, 468, 1080

Meshkat, T., Bailey, V. P., Su, K. Y. L., et al. 2015, ApJ, 800, 5

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

Millar-Blanchaer, M. A., Wang, J. J., Kalas, P., et al. 2016, AJ, 152, 128

Milli, J., Engler, N., Schmid, H. M., et al. 2019, A&A, 626, A54

Milli, J., Mouillet, D., Lagrange, A.-M., et al. 2012, A&A, 545, A111

Milli, J., Vigan, A., Mouillet, D., et al. 2017, A&A, 599, A108

Min, M., Hovenier, J. W., & de Koter, A. 2005, A&A, 432, 909

Min, M., Rab, C., Woitke, P., et al. 2016, A&A, 585, A13

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

Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, AJ, 158, 13

Ochsenbein, F., Bauer, P., & Marcout, J. 2000, A&AS, 143, 23

Olofsson, J., Milli, J., Thébault, P., et al. 2019, A&A, 630, A142

Pan, M., Nesvold, E. R., & Kuchner, M. J. 2016, ApJ, 832, 81

Pavlov, A., Feldt, M., & Henning, T. 2008, in ASP Conf. Ser. 394, Data Reduction and Handling for SPHERE, ed. R. W. Argyle, P. S. Bunclark, & J. R. Lewis (San Francisco, CA: ASP), 581

Perrin, M. D., Duchene, G., Millar-Blanchaer, M., et al. 2015, ApJ, 799, 182

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

Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

Pueyo, L. 2016, ApJ, 824, 117

Ren, B., Choquet, É., Perrin, M. D., et al. 2019, ApJ, 882, 64

Ren, B., Pueyo, L., Chen, C., et al. 2020, ApJ, 892, 74

Ren, B., Pueyo, L., Zhu, G. B., et al. 2018, ApJ, 852, 104

Rodigas, T. J., Stark, C. C., Weinberger, A., et al. 2015, ApJ, 798, 96

Schneider, G., Debes, J. H., Grady, C. A., et al. 2018, AJ, 155, 77

Schneider, G., Smith, B. A., Becklin, E. E., et al. 1999, ApJL, 513, L127

Schneider, G., Weinberger, A. J., Becklin, E. E., et al. 2009, AJ, 137, 53

Schröder, S. E., Grynkó, Y., Pommerol, A., et al. 2014, Icar, 239, 201

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

Smart, W. M. 1930, MNRAS, 90, 534

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

Stark, C. C., Schneider, G., Weinberger, A. J., et al. 2014, ApJ, 789, 58

Stauffer, J. R., Hartmann, L. W., & Barrado y Navascues, D. 1995, ApJ, 454, 910

Telesco, C. M., Fisher, R. S., Piña, R. K., et al. 2000, ApJ, 530, 329

Thalmann, C., Janson, M., Buenzli, E., et al. 2011, ApJL, 743, L6

Wahhaj, Z., Liu, M. C., Biller, B. A., et al. 2014, A&A, 567, A34

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

Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67

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