# Multi-epoch Direct Imaging and Time-variable Scattered Light Morphology of the HD 163296 Protoplanetary Disk

Evan Rich, John Wisniewski, Thayne Currie, Misato Fukagawa, Carol A Grady, Michael Sitko, Monika Pikhartova, Jun Hashimoto, Lyu Abe, Wolfgang Brandner, et al.

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

To cite this version:

Evan Rich, John Wisniewski, Thayne Currie, Misato Fukagawa, Carol A Grady, et al. Multi-epoch Direct Imaging and Time-variable Scattered Light Morphology of the HD 163296 Protoplanetary Disk. The Astrophysical Journal, 2019, 875 (1), pp.38. 10.3847/1538-4357/ab0f3b. hal-02303546

HAL Id: hal-02303546

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

Submitted on 12 Aug 2025

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

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

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

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

# Multi-epoch Direct Imaging and Time-variable Scattered Light Morphology of the HD 163296 Protoplanetary Disk

Evan A. Rich<fnref n="1" />, John P. Wisniewski<fnref n="1" />, Thayne Currie<fnref n="2" />,<fnref n="3" />,<fnref n="4" />, Misato Fukagawa<fnref n="5" />,<fnref n="6" />, Carol A. Grady<fnref n="2" />,<fnref n="7" />,<fnref n="8" />, Michael L. Sitko<fnref n="9" />,<fnref n="10" />, Monika Pikhartova<fnref n="9" />, Jun Hashimoto<fnref n="11" />, Lyu Abe<fnref n="12" />, Wolfgang Brandner<fnref n="13" />, Timothy D. Brandt<fnref n="14" />, Joseph C. Carson<fnref n="15" />, Jeffrey Chilcote<fnref n="16" />, Ruobing Dong<fnref n="17" />, Markus Feldt<fnref n="13" />, Miwa Goto<fnref n="18" />, Tyler Groff¹⁹, Olivier Guyon<fnref n="3" />,<fnref n="20" />,<fnref n="21" />, Yutaka Hayano<fnref n="3" />, Masahiko Hayashi<fnref n="5" />, Saeko S. Hayashi<fnref n="3" />, Thomas Henning<fnref n="13" />, Klaus W. Hodapp<fnref n="22" />, Miki Ishii<fnref n="5" />, Masanori Iye<fnref n="5" />, Markus Janson<fnref n="23" />, Nemanja Jovanovic<fnref n="24" />, Ryo Kandori<fnref n="5" />, Jeremy Kasdin<fnref n="25" />, Gillian R. Knapp<fnref n="26" />, Tomoyuki Kudo<fnref n="3" />, Nobuhiko Kusakabe<fnref n="5" />, Masayuki Kuzuhara<fnref n="3" />,<fnref n="11" />, Jungmi Kwon<fnref n="27" />, Julien Lozi<fnref n="3" />, Frantz Martinache<fnref n="28" />, Taro Matsuo<fnref n="29" />, Satoshi Mayama<fnref n="30" />, Michael W. McElwain<fnref n="7" />, Shoken Miyama<fnref n="31" />, Jun-Ichi Morino<fnref n="5" />, Amaya Moro-Martin<fnref n="32" />,<fnref n="33" />, Takao Nakagawa<fnref n="34" />, Tetsuo Nishimura<fnref n="3" />, Tae-Soo Pyo<fnref n="3" />, Eugene Serabyn<fnref n="35" />, Hiroshi Suto<fnref n="5" />, Ray W. Russel<fnref n="36" />, Ryuji Suzuki<fnref n="5" />, Michihiro Takami<fnref n="37" />, Naruhisa Takato<fnref n="3" />, Hiroshi Terada<fnref n="3" />, Christian Thalmann<fnref n="38" />, Edwin L. Turner<fnref n="14" />, Taichi Uyama<fnref n="27" />, Kevin R. Wagner<fnref n="21" />, Makoto Watanabe<fnref n="39" />, Toru Yamada<fnref n="40" />, Hideki Takami<fnref n="5" />, Tomonori Usuda<fnref n="27" />, and Motohide Tamura<fnref n="5" />,<fnref n="27" />

Received 2018 November 16; revised 2019 March 4; accepted 2019 March 7; published 2019 April 12

## Abstract

We present H-band polarized scattered light imagery and JHK high-contrast spectroscopy of the protoplanetary disk around HD 163296 observed with the High-Contrast Coronographic Imager for Adaptive Optics (HiCIAO) and Subaru Coronagraphic Extreme Adaptive Optics (SCExAO)/Coronagraphic High Angular Resolution Imaging Spectrograph (CHARIS) instruments at Subaru Observatory. The polarimetric imagery resolve a broken ring structure surrounding HD 163296 that peaks at a distance along the major axis of 0.65 (66 au) and extends out to 0.98 (100 au) along the major axis. Our 2011 H-band data exhibit clear axisymmetry, with the NW and SE side of the disk exhibiting similar intensities. Our data are clearly different from 2016 epoch H-band observations of the

<footnote n="1">¹Homer L. Dodge Department of Physics, University of Oklahoma, Norman, OK 73071, USA; erich66210@ou.edu, wisniewski@ou.edu</footnote>

<footnote n="2">²Eureka Scientific, 2452 Delmer, Suite 100, Oakland, CA 96002, USA</footnote>

<footnote n="3">³National Astronomical Observatory of Japan, Subaru Telescope, National Institutes of Natural Sciences, 650 North Aohoku Place, Hilo, HI 96720, USA</footnote>

<footnote n="4">⁴NASA-Ames Research Center, Moffett Field, CA, USA</footnote>

<footnote n="5">⁵National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo, 181-8588, Japan</footnote>

<footnote n="6">⁶Graduate School of Science, Osaka University, 1-1 Machikaneyama, Toyonaka, Osaka 560-0043, Japan</footnote>

<footnote n="7">⁷Exoplanets and Stellar Astrophysics Laboratory, Code 667, Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="8">⁸Goddard Center for Astrobiology, Goddard Center for Astrobiology, Greenbelt, MD 20771, USA</footnote>

<footnote n="9">⁹Department of Physics, University of Cincinnati, Cincinnati, OH 45221, USA</footnote>

<footnote n="10">¹⁰Space Science Institute, 475 Walnut Street, Suite 205, Boulder, CO 80301, USA</footnote>

<footnote n="11">¹¹Astrobiology Center of NINS, 2-26-1, Osawa, Mitaka, Tokyo, 231-8588, Japan</footnote>

<footnote n="12">¹²Laboratoire Lagrange (UMR 7293), Université de Nice-Sophie Antipolis, CNRS, Observatoire de la Cote d'Azur, 28 avenue Valrose, F-06108 Nice Cedex 2, France</footnote>

<footnote n="13">¹³Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="14">¹⁴Astrophysics Department, Institute for Advanced Study, Princeton, NJ 08540, USA</footnote>

<footnote n="15">¹⁵Department of Physics and Astronomy, College of Charleston, 58 Coming St., Charleston, SC 29424, USA</footnote>

<footnote n="16">¹⁶Department of Physics, University of Notre Dame, 225 Nieuwland Science Hall, Notre Dame, IN 46556, USA</footnote>

<footnote n="17">¹⁷Department of Physics & Astronomy, University of Victoria, Victoria, BC V8P 1A1, Canada</footnote>

<footnote n="18">¹⁸Universitäts-Sternwarte München, Ludwig-Maximilians-Universität, Scheinerstr. 1, D-81679 München, Germany</footnote>

<footnote n="20">²⁰Astrobiology Center, National Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo, Japan</footnote>

<footnote n="21">²¹Steward Observatory, University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="22">²²Institute for Astronomy, University of Hawaii, 640 N. A'ohoku Place, Hilo, HI 96720, USA</footnote>

<footnote n="23">²³Department of Astronomy, Stockholm University, AlbaNova University Center, SE-10691 Stockholm, Sweden</footnote>

<footnote n="24">²⁴Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="25">²⁵Department of Mechanical and Aerospace Engineering Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="26">²⁶Department of Astrophysical Science, Princeton University, Peyton Hall, Ivy Lane, Princeton, NJ 08544, USA</footnote>

<footnote n="27">²⁷Department of Astronomy, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan</footnote>

<footnote n="28">²⁸Observatoire de la Cote d'Azur, Boulevard de l'Observatoire, F-06300 Nice, France</footnote>

<footnote n="29">²⁹Department of Astronomy, Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo-ku, Kyoto 606-8502, Japan</footnote>

<footnote n="30">³⁰The Center for the Promotion of Integrated Sciences, The Graduate University for Advanced Studies (SOKENDAI), Shonan International Village, Hayama-cho, Miura-gun, Kanagawa 240-0193, Japan</footnote>

<footnote n="31">³¹Hiroshima University, 1-3-2, Kagamiyama, Higashi-Hiroshima 739-8511, Japan</footnote>

<footnote n="32">³²Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="33">³³Center for Astrophysical Sciences, Johns Hopkins University, Baltimore, MD 21218, USA</footnote>

<footnote n="34">³⁴Department of Space Astronomy and Astrophysics Institute of Space & Astronautical Science (ISAS) Japan Aerospace Exploration Agency (JAXA) 3-1-1 Yoshinodai, Chuo-ku, Sagamihara, Kanagawa 252-5210, Japan</footnote>

<footnote n="35">³⁵Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA</footnote>

<footnote n="36">³⁶The Aerospace Corporation, Los Angeles, CA, USA</footnote>

<footnote n="37">³⁷Institute of Astronomy and Astrophysics, Academia Sinica 11F of Astronomy-Mathematics Building, AS/NTU No.1, Sec. 4, Roosevelt Rd, Taipei 10617, Taiwan, R.O.C</footnote>

<footnote n="38">³⁸Institute for Astronomy, ETH Zurich, Wolfgang-Pauli-Strasse 27, 8093 Zurich, Switzerland</footnote>

<footnote n="39">³⁹Department of Cosmosciences, Hokkaido University, Kita-ku, Sapporo, Hokkaido 060-0810, Japan</footnote>

<footnote n="40">⁴⁰Astronomical Institute, Tohoku University, Aoba-ku, Sendai, Miyagi 980-8578, Japan</footnote>
Very Large Telescope (VLT)/Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE), which found a strong 2.7× asymmetry between the NW and SE side of the disk. Collectively, these results indicate the presence of time-variable, non-azimuthally symmetric illumination of the outer disk. While our SCExAO/CHARIS data are sensitive enough to recover the planet candidate identified from NIRC2 in the thermal infrared (IR), we fail to detect an object with JHK brightness nominally consistent with this object. This suggests that the candidate is either fainter in JHK bands than model predictions, possibly due to extinction from the disk or atmospheric dust/ clouds, or that it is an artifact of the data set/data processing, such as a residual speckle or partially subtracted disk feature. Assuming standard hot-start evolutionary models and a system age of 5 Myr, we set new, direct mass limits for the inner (outer) Atacama Large Millimeter/submillimeter Array (ALMA)-predicted protoplanet candidate along the major (minor) disk axis of 1.5 (2) $M_J$.

*Key words:* planets and satellites: detection – protoplanetary disks

*Supporting material:* data behind figure

# 1. Introduction

Protoplanetary disks are dust and gas disks around young stars that guide the accretion of material onto forming stars and serve as the birthplace of planets. Direct imaging of protoplanetary disks reveals likely sites of active planet formation, may identify planets in the final stages of assembly (protoplanets), and probes the interaction between protoplanets and the disk material from which they form. Herbig Ae/Be stars (Herbig 1960), the intermediate mass analogs to T Tauri stars, are known to both host protoplanetary disks and often exhibit evidence of ejecting material via collimated, bipolar jets (Herbig 1950; Grady et al. 2000; Ellerbroek et al. 2014; Bally 2016). The protoplanetary disks around Herbig Ae/Be stars exhibit a variety of structures—with some hosting spiral arms (Hashimoto et al. 2011), and others that are flat and settled, causing self-shadowing of the disk (Meeus et al. 2001)—and may host some of the first directly imaged Jovian protoplanets (Quanz et al. 2013; Currie et al. 2015).

HD 163296 is a young—5.1$^{+0.3}_{-0.8}$ Myr old (Montesinos et al. 2009) to 7.6$^{+1.1}_{-1.2}$ (Vioque et al. 2018)—Herbig Ae protoplanetary disk system located at a distance of 101.5 ± 1.2 pc (Gaia Collaboration et al. 2016, 2018). The disk has been spatially resolved by ground- and space-based observing platforms at a multitude of wavelengths, including optical (Hubble Space Telescope (HST)/Space Telescope Imaging Spectrograph (STIS): Grady et al. 2000, HST/Advanced Camera for Surveys (ACS): Wisniewski et al. 2008), near-infrared (IR) (Very Large Telescope (VLT)/NACO: Garufi et al. 2014, 2017, Gemini/Gemini Planet Imager (GPI): Monnier et al. 2017, VLT/Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE): Muro-Arena et al. 2018, Subaru/Coronagraphic Imager with Adaptive Optics (CIAO): Fukagawa et al. 2010, Keck/NIRC2: Guidi et al. 2018), and radio wavelengths (Very Large Array (VLA): Guidi et al. 2016, ALMA: Guidi et al. 2016; Isella et al. 2016).

Spatially resolved imaging observations have revealed a complex circumstellar environment and evidence for active planet formation at wide separations around HD 163296. Its disk extends to at least to 4.′′4 (447 au) in optical scattered light (Wisniewski et al. 2008). While near-IR observations reveal a 64 au-scale inner dust ring (Garufi et al. 2014, 2017; Monnier et al. 2017; Muro-Arena et al. 2018), 1.3 mm continuum ALMA imaging (Isella et al. 2016) revealed three azimuthal gaps in the disk located at 0.′′49, 0.′′82, and 1.′′31 (50, 83, and 133 au respectively given Gaia-DR2 distance of 101.5 pc). The surface distribution of small dust grains in the outer disk appears low, owing to settling or partial-to-complete depletion (Muro-Arena et al. 2018). Keck/NIRC2 thermal infrared

imaging led to the discovery of a candidate 7 $M_J$ protoplanet just exterior to the inner ring (Guidi et al. 2018), while modeling of ALMA gas emission data suggest Jovian planets at 83 and 137 au (Teague et al. 2018) and/or a single Jovian on an even wider orbit (260 au Pinte et al. 2018).

Multi-epoch observations have revealed a wealth of variability in the HD 163296 system, likely traceable to dynamical processes in the inner disk region. Both IR spectra and visibilities from optical interferometry show variability possibly connected to changes in the inner disk or the system's wind component (Sitko et al. 2008; Tannirkulam et al. 2008). Long-term optical photometric and IR spectroscopic monitoring revealed suggestive evidence of a 16 yr periodicity, with optical fluxes dimming when the IR fluxes reach a maximum level (Sitko et al. 2008; Ellerbroek et al. 2014), on similar timescales as the ejection of Herbig-Haro objects (Ellerbroek et al. 2014). The star's accretion rate increased more than 1 dex over ~15 yr (Mendigutía et al. 2013). However, no clear correlation between these variations and the 16 yr optical infrared periodicity has yet been found. CO rovibrational emission lines exhibit variability possibly connected to changes in the disk wind or episodic accretion (Hein Bertelsen et al. 2016).

Spatially resolved imaging may also reveal evidence for variability—time-dependent changes in the disk's surface brightness and morphology potentially linked to variable illumination (Wisniewski et al. 2008). However, despite this plethora of variability observed, the lack of contemporaneous observations of both the inner and outer regions of the HD 163296 disk limits efforts to connect these phenomena to one another.

In this paper, we present multi-epoch near-infrared scattered light imaging of HD 163296, obtained at H-band in polarized light as part of the Strategic Exploration of Exoplanets and Disks with Subaru survey (Tamura 2009) and in total intensity in JHK using Subaru Coronagraphic Extreme Adaptive Optics (SCExAO) (Jovanovic et al. 2015) coupled with the Coronagraphic High Angular Resolution Imaging Spectrograph (CHARIS) integral field spectrograph (Section 2). To help parse and complement these data probing the outer disk, we acquired near-contemporary IR spectra to characterize the inner disk region of the system. We modeled the H-band scattered light images and near-IR spectra using a well-established 3D Monte Carlo Radiative Transfer (MCRT) code to create a more coherent, full picture of the system at this epoch (Section 3). Using our model, we preformed forward modeling of the SCExAO/CHARIS data to help constrain the presence of potential planets in the system (Section 4). Finally we discuss the implications of our results in Section 5 including deeper constrains on protoplanets around HD 163296 with the new SCExAO/CHARIS data.
## 2. Observations and Data Reduction

### 2.1. HiCIAO Imagery

We obtained high-contrast H-band imaging of HD 163296 using the High-Contrast Coronographic Imager for Adaptive Optics (HiCIAO) instrument (Hodapp et al. 2008) along with the AO-188 system (Hayano et al. 2008, 2010) at the Subaru Observatory on 2011 August 3. We used a circular occulting mask having a diameter of 0.3, and observed the system in standard Polarized Differential Imaging mode at four wave-plate positions (0°, 22.5°, 45°, 67.5°). We obtained 72 frames using 30 s exposures, yielding a total of 18 complete wave-plate sets. We determined that 8 wave-plate sets had lower AO performance, and discarded them during the reduction of the data. We also obtained a short, direct H-band photometric observation of HD 163296, and determined that the source's brightness at this epoch was 5.62 ± 0.05 mag.

We reduced our observations using standard double differencing techniques, as described in Hashimoto et al. (2011). To briefly summarize, the two subimages of each frame contain an ordinary and an extra-ordinary image, which can be summed and subtracted from their 90° counterparts to create stokes parameter $-Q, +Q, -U, and +U$ images. The Q and U frames were then rotated into a common orientation, corrected for instrumental polarization, and summed to create final Q and U images. We corrected these data for the presence of a residual polarized halo having the properties of $p = 1.00 \pm 0.05\%$ and $\theta = 42.5 \pm 1.5^\circ$. Final polarized intensity (PI) imagery was created from the total Q and U data, using $PI = \sqrt{Q^2 + U^2}$, as shown in Figure 1.

To further simplify the analysis of our imagery, we adopt the now common practice of assuming single scattering, and rotated all of the light that is polarized perpendicular to the star by the angle $\phi$ into a $Q_\phi$ image and all of the light that is polarized parallel to the star into a $U_\phi$ image as defined below (Schmid et al. 2006).

$$Q_{\phi} = Q \times \cos 2\phi + U \times \sin 2\phi \quad (1)$$

$$U_{\phi} = Q \times \sin 2\phi + U \times \cos 2\phi. \quad (2)$$

The final $Q_\phi$ and $U_\phi$ imagery for HD 163296 are shown in Figure 1. Little coherent signal appears present in the $U_\phi$ image, which helps confirm that little residual instrumental contaminants remain in these data. Next, we computed a signal-to-noise (SN) image, following the procedure outlined by Ohta et al. (2016). In summary, we computed the noise by measuring the standard deviation of every pixel in each of the Q and U frames used to construct the final imagery, then divided by the square root of the number of frames. The resultant SN image is shown in Figure 1.

### 2.2. Near-infrared Spectra from SpeX, Broadband Array Spectrograph System (BASS), and TripleSpec

We also observed HD 163296 multiple times with several near-IR instruments on NASA's Infrared Telescope Facility (IRTF) and at Apache Point Observatory (APO). We observed HD 163296 using the SpeX spectrograph (Rayner et al. 2003) at IRTF in its short-wavelength mode (0.8–2.4 µm) and long-wavelength mode (2.3–5.5 µm) on 2011 July 31, 2016 May 4, and 2018 June 24. These observations are contemporaneous with the HiCIAO 2011 observation (Section 2.1), the Gemini/ GPI observation (Section 5.1), and the second SCExAO/ CHARIS observations (Section 2.3) respectively. We observed HD 163296 using the TripleSpec spectrograph (Wilson et al. 2004) at the APO 3.5 m telescope, covering a spectral range of 0.95–2.46 µm, on 2018 May 16. This observation is

contemporaneous with the first SCExAO/CHARIS observation (Section 2.3). We observed the nearby A0V star HD 163336 to perform telluric corrections for both the SpeX and TripleSpec observations. These data were reduced and calibrated using the standard reduction packages Spextool and Triplespectool (Vacca et al. 2003; Cushing et al. 2004). We also observed HD 163296 with The Aerospace Corporation's BASS, which covers two wavelength bands from 2.9–6 µm and 6–13.5 µm respectively, on 2011 August 1. HD 163336 was observed with BASS to flux calibrate these data. The instrument and data reduction method are fully described in Wagner et al. (2015). These SpeX, TripleSpec, and BASS spectra are plotted in Figure 2.

### 2.3. SCExAO/CHARIS High-contrast Near-infrared Spectroscopy

We observed HD 163296 on 2018 May 22 and July 1 at the Subaru Observatory with SCExAO coupled with the CHARIS integral field spectrograph operating in low-resolution ($R \sim 20$) and broadband (1.13–2.39 µm) mode, covering the JHK filters simultaneously (Groff et al. 2015). For the May observations, the conditions were stable with 0.4 seeing and 6–7 m s⁻¹ winds. Our observations consisted of co-added 60.4 s frames totaling ~30 minutes of integration time and covering a modest parallactic angle rotation (ΔPA = 14.8°). Due to highly variable conditions for the July observations, we obtained shorter exposures (30.9 s) and removed roughly 50% of the frames with poor AO correction, yielding ~40 minutes of data covering 30°9 of parallactic angle motion.<fnref n="41" />

We followed the standard setup used for SCExAO/CHARIS broadband observations (Currie et al. 2018a; Goebel et al. 2018a), using the Lyot coronagraph with the 217 mas occulting spot and bracketing our coronagraphic sequence with blank sky frames to remove sky emission and instrumental artifacts. We used satellite spots produced from a 25 nm modulation on SCExAO's deformable mirror for spectrophotometric calibration and image registration (Jovanovic et al. 2015). For data cube extraction, we utilized the least-squares algorithm from the CHARIS Data Reduction Pipeline (Brandt et al. 2017). Basic data processing, including sky subtraction, image registration, etc., follows methods used for recent SCExAO/CHARIS broadband studies (Currie et al. 2018a, 2018b; Goebel et al. 2018b).

Spectrophotometrically calibrating CHARIS data for pre- transitional disk sources like HD 163296 requires either observations of a separate spectral standard or contemporaneous near-IR spectra. We opt for the latter, using the IRTF/SpeX and APO/Triplespec data previously discussed in Section 2.2. The spectra show only minor differences between epochs.

We explored a range of point-spread function (PSF) subtraction approaches leveraging on angular differential imaging (ADI; Marois et al. 2006), spectral differential imaging (SDI; Sparks & Ford 2002), and combinations of the two (ASDI, e.g., Marois et al. 2014). We further considered a variety of PSF subtraction algorithms, including A-LOCI (Currie et al. 2012, 2018a), Karhunen–Loève Image Projection (KLIP; Soummer et al. 2012), and classical PSF subtraction (Marois et al. 2006). The approach

<footnote n="41">⁴¹ While a real-time estimate of the Strehl ratio (S.R.) was not recorded for these data sets, the raw contrast for the May data was just slightly poorer than that obtained for $\kappa$. And observations achieving S.R. ~0.90–0.92 in H band (Currie et al. 2018a). Raw contrasts for the July data considered in our study are roughly a factor of 2.5–3 worse at 0.4, more characteristic of performance at S.R. ~0.65–0.70.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** H-band scattered light from the HD 163296 disk is clearly seen in polarized intensity (PI) (panel A), the S/N map (panel B), and $Q_{\phi i}$ imagery (panel C). Little coherent signal is seen in the $U_{\phi i}$ image (panel D), indicating that these data are largely free from PSF residuals. The PI (panel A), $Q_{\phi i}$ (panel C), and $U_{\phi i}$ (panel D) images are displayed on a linear scale with units of mJy, and have not been filtered. We have applied a software mask having a radial size of 0.3 (gray circles) to match the effective inner working angle of these data. For all panels, north is up and east is to the left. The *Q* and *U* images shown in panels C and D of this figure are available as the data behind the Figure. The data used to create this figure are available.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Five epochs of flux calibrated IR spectra of HD 163296, taken with IRTF/SpeX, IRTF/BASS, or APO/TripleSec, are shown. A full description of these observations can be found in Section 2.2. The spectra are plotted in log-log space.</figcaption></figure>

implemented for $\kappa$ And in Currie et al. (2018a), using A-LOCI to subtract the PSF in ADI and then again to remove residuals in SDI mode, yielded the best speckle suppression while preserving the signal from the disk. Due to the limited parallactic angle motion of both data sets (especially in May) and the presence of the disk, we utilized large optimization zones for the ADI step, employed local masking in the SDI step, and imposed a rotation/magnification criterion of $\delta = 0.5–1.0$ PSF footprints in both steps to construct a reference PSF (see Lafrenière et al. 2007). For both steps, we used a singular value decomposition (SVD) cutoff of $10^{-6}$ to solve the set of linear equations that result in the weighted reference PSF for each region of each data cube slice (see Currie et al. 2015).

Figure 3 shows broadband (wavelength-collapsed) CHARIS images of HD 163296 from SCExAO/CHARIS for the May (left) and July (right) epochs after removing the stellar PSF through both ADI and SDI. Despite poor field rotation (May data) or variable conditions (July data), we clearly detect the outer ring of emission seen in polarimetry, which appears as a sharply defined crescent defining the forward-scattering edge of
<figure><img src="image_6.png" /><figcaption>**Figure 3.** SCExAO/CHARIS broadband (wavelength-collapsed) images from 2018 May (left) and 2018 July (right) after removing the stellar PSF through both ADI and SDI: the color scaling for both panels goes from −30 to 30 mJy arcsec$^{-2}$. In both data sets, self-subtraction footprints (dark regions) flank the disk signal, which is reduced due to processing. The throughput of the disk is slightly higher in the July data due to better field rotation; regions surrounding the disk show slightly less residual speckle noise in the May data due to better AO performance.</figcaption></figure>

the structure. Self-subtraction footprints due to both ADI and SDI flank the ring. In individual passbands, the disk is just marginally visible in *J*-band but is well separated from residual speckle noise in *H* and *K*.

We defined a conservative lower limit to the signal-to-noise ratio (S/N) of the trace of the disk in broad band, adopting the standard practice of replacing each pixel by the sum within its aperture, defining a radial-dependent noise profile, and applying a finite-element correction for the noise (Currie et al. 2011; Mawet et al. 2014). To be conservative, we include signal from the disk in our estimate of the noise profile. Except at the semiminor axis, where the disk signal is attenuated by self-subtraction, the disk trace is decisively detected, with a S/N per resolution element ranging from 3 to 8.5.

Our data do not reveal the candidate protoplanet identified in Keck/NIRC2 $L_p$ data from Guidi et al. (2018) nor the companions predicted from ALMA data (Teague et al. 2018). The inner disk seen by GPI polarimetry (Monnier et al. 2017) is also not visible, likely due to heavy self-subtraction due to poor field rotation. The position of the Guidi et al. candidate lies well separated from the ring and residual speckle noise; the S/N maps show no convolved pixel within one PSF footprint (~0.08) of this position with a significance greater than 1.3σ. More conservative reductions (e.g., larger rotation gap; higher SVD cutoff) may show slightly elevated residual emission consistent with additional extended structure at this separation (e.g., additional ring material). However, this signal is not statistically significant and is simpler to explain as residual speckle noise instead.

### 3. Analysis of the *H*-band Polarimetry Data

In this section, we characterize the distribution of scattered light in our *H*-band imagery, and construct a MCRT model to help interpret the contemporaneous *H*-band scattered light imagery and near-IR spectra.

#### 3.1. Geometry of the Disk

Figure 1 reveals the clear detection of scattered light surrounding the HD 163296 disk in our *H*-band imagery outside of the inner working angle of these data, 0.′′3 (30.5 au). The scattered light imagery reveals a broken ring structure that peaks at a distance along the major axis of 0.′′65 (66 au) and extends out to 0.′′98 (100 au) along the major axis (see Figure 4). Both the $Q_φ$ and SN imagery exhibit little coherent

signal between our inner working angle and the inner edge of the ring structure. We do not detect the inner disk component as previously detected by Monnier et al. (2017) due to our larger inner working angle. We therefore conclude that the small amount of scattered light interior to the ring in the PI image (panel A; Figure 1) could arise from a mixture of residual, uncorrected flux from the PSF and scattered flux from an inner disk component that is within our masked region (see, e.g., Takami et al. 2018). The NE side of the disk is known to be the near side (Rosenfeld et al. 2013), and the IR scattered light disk exhibits evidence of strong forward scattering (Guidi et al. 2018). The broken ring structure we observe is missing polarized intensity originating from the far side of the disk (SW region, along the minor axis; see Figure 1).

We fit an ellipse to the scattered light ring using a least-squares fitting code (Hammel & Sullivan-Molina 2019),<fnref n="42" /> assuming the ring is a perfect circle projected at inclination. Since a known bias of the code is to prefer a smaller ellipse by preferably fitting the inner points (Halir & Flusser 1998), we choose to fit the peaks of the ring to mitigate this effect. Due to the low signal along the SW minor axis and sporadic structure along the NE minor axis, we did not keep any vertical cuts between 70° < PA < 180° and between 270° < PA < 370°. We fit a Gaussian to each vertical crosscut, producing the peak x, y position of the ring, and input these positions into the ellipse code described above.

In order to estimate the error of our ellipse fit, we performed a Monte Carlo routine by randomly sampling the Gaussian xy-coordinate errors and adding them to the xy-coordinates found above. We additionally applied a random rotation of the image between 0° and 1° to constrain the error associated with the interpolation of the image due to rotation. We performed 500 iterations and used the average values of the 500 iterations as the best-fit ellipse. The errors were estimated by taking the standard deviation of the parameters found with the 500 iterations. The best-fit results are shown in Table 1. The best-fit ellipse is compared with the PI image in Figure 5, shown as the white oval, along with the center of the disk (small white circle) and the center of the star.

Our measured inclination of the disk (41°4 ± 0°3) and PA (132°2 ± 0°3) is in agreement with the values derived from ALMA data of 42° and 132° respectively (Isella et al. 2016). Additionally, the offset of the minor axis from the central star

<footnote n="42">⁴² https://github.com/bdhammel/least-squares-ellipse-fitting</footnote>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Crosscuts along the major axis of the 2011 H-band PI image (top row) and $Q_φ$ image (bottom row). The right column is the PI and $Q_{\text{phi}}$ images unscaled, and the left column is the PI and $Q_φ$ with a $r^2$ scaling applied. The gray shaded area represents $3\sigma$ error bars. The red point is the scaled flux from the 2016 VLT/SPHERE observation reported by Muro-Arena et al. (2018).</figcaption></figure>

<table><caption>**Table 1**
Results of Ellipse Fitting to PI H-band Image</caption><thead><tr><td>Parameter</td><td>PI Image Value</td></tr></thead><tbody><tr><td>Major Axis of Disk (au)</td><td>58.01 ± 0.09</td></tr><tr><td>Minor Axis of Disk (au)</td><td>48.4 ± 0.3</td></tr><tr><td>Minor Axis offset (′′)</td><td>-0.0432 ± 0.0016</td></tr><tr><td>PA (deg)</td><td>132.2 ± 0.3</td></tr><tr><td>Inclination (°)</td><td>41.4 ± 0.3</td></tr></tbody></table>

that we find (−0.0432 ± 0.0016) is consistent with previous measurements, given their quoted errors when available (0.06, Garufi et al. 2014; 0.105 ± 0.045, Muro-Arena et al. 2018; 0.1, Monnier et al. 2017).

We applied an $r^2$ illumination correction to our data to better investigate the physical distribution of dust in the ring seen in Figure 1. We then azimuthally binned the average flux per area of the ring between two concentric ellipses. We adopted an inclination of 42°, and constructed each bin to be 8° wide and spanned a projected radial distance of 0.55–0.71 (55–72.5 au)

to encompass the majority of the disk flux. The binned disk flux is azimuthally symmetric along the major axis, with the NW and SE side of the disks exhibiting the same amount of polarized intensity (Figure 6). There is also a clear azimuthal asymmetry in the binned flux along the minor axis, with the near side of the disk (NE side) exhibiting substantially more flux than the far side (SW side). We observe a deficit in scattered light flux along the near side of the disk at a PA of 30° in both the binned imagery (Figure 6) and unbinned PI, $Q_{\text{rot}}$ images, and SN map images (Figure 1). This feature coincides with the position angle of the disk brightness enhancement and the position angle of the candidate point source noted by Guidi et al. (2018), and will be further discussed in Section 5.4.

## 3.2. Modeling of the HD 163296 Disk

To help interpret our imagery and contemporaneous IR spectroscopy of HD 163296, we utilized the 3D MCRT code, HOCHUNK3D (Whitney et al. 2013). HOCHUNK3D allows the user to characterize the radial dust distribution, dust composition, and disk illumination parameters, and outputs a

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

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 5.** Result of the best-fit ellipse to our *H*-band PI data, where the central white dot is the center of the ellipse, the white ellipse is the peak of the ellipse, the black x marks the location of the star, and the blue circle marks the inner working angle. The ellipse was fit to the peak points along the main elliptical ring by fitting Gaussians to the crosscuts along the ring. The best elliptical fit finds a minor axis offset of −0.055. This value is consistent with those reported by Garufi et al. (2014), Monnier et al. (2017), and Muro-Arena et al. (2018) given their quoted uncertainties.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 6.** Binned flux along the azimuthal ring located at 65 au. Each bin is 8° wide and extends from a projected distance of 55–71 au annulus along the ring seen in this figure.</figcaption></figure>

spectral energy distribution (SED) of the disk and imagery in a variety of user-defined bandpasses. The current version of HOCHUNK3D allows the user to decouple the disk into two dust distributions, allowing one to parameterize both a settled dust population toward the midplane and a different dust population in the upper surface layers of the disk. These two dust populations can either be co-spatial, or have different radial sizes. The dust distribution of each disk is characterized by several power-law parameters: the radial power law ($\alpha$), the vertical Gaussian distribution ($\beta$), and the height of the disk

from the midplane ($h$). Deviations from these power laws, such as a gap, spiral arms, warped disks, and walls, can all be included. The code also allows for the presence of a dusty envelope, which is parameterized by its minimum and maximum radius ($R_{minenv}$, $R_{maxenv}$), and a dust density power law (ENVEXP). The dusty envelope can also include gaps and a bipolar cavity. Following the techniques established by Sitko et al. (2008), Wagner et al. (2015), and Fernandes et al. (2018), we use the dusty envelope as a proxy to model material ejected from the disk, aka a disk wind.

We constrained our model-starting parameters by observations when possible, and adopted the parameters from M. Pikhartova et al. (2019, in preparation), who are using HOCHUNK3D to model the variations seen in two epochs of HD163296's SED, as a starting point for our parameter-space exploration. ALMA observations of HD 163296 revealed the presence of 3 gaps located at 0.49, 0.82, and 1.31 (50, 83, and 133 au respectively given Gaia-DR2 distance of 101.5 pc) (Isella et al. 2016). Because our HiCIAO imagery is only sensitive to the first dust ring and the near-IR SED is most sensitive to dust features closer to the star, we only include the inner gap in our model. We allowed the two components of the dust distribution to be vertically stratified, and chose the radial extent of these distributions to match those observed for grains populating the midplane (250 au from VLA and ALMA observations; Guidi et al. 2016) and surface layers (540 au, Isella et al. 2007; Wisniewski et al. 2008). We note that while ALMA observations of the system were best described by a radial power law multiplied by an exponential function (Isella et al. 2016), HOCHUNK3D only uses a power-law function. Nevertheless, we did adjust the large grain dust distribution to match, as closely as possible, the dust distribution as measured by ALMA in the inner portion of the disk (Isella et al. 2016). The dust parameters for the large grain disk that we used are adopted from Wood et al. (2002), and are composed of amorphous carbon and silicon dust particles ranging in size up to 1 millimeter. The small grain disk and envelope dust parameters are from Kim et al. (1994), which is the average galactic interstellar medium (ISM) dust grain model.

We adopted an interstellar extinction of $A_V = 0$ mag from Ellerbroek et al. (2014), who measured the level of extinction from the ejected Herbig-Haro (HH) knots. Note that Ellerbroek et al. (2014) concluded that the optical variability of the SED likely comes from on-source reddening. In our model, we utilize the dusty envelope, a proxy to model disk wind, to replicate the on-source reddening, which is further discussed in Section 5.2. We explored accretion rates ranging from $1.73 \times 10^{-7}$ to $4.35 \times 10^{-6} M_{\odot}$, calculated from contemporaneous Brγ emission lines in the SpeX 2011 data, first presented in Ellerbroek et al. (2014), but adjusted for the new distance of 101.5 pc.

We constrained these models using a SED (Figure 9) constructed from contemporaneous near-IR observations (Figure 2), along with non-contemporaneous photometry from the AllWISE catalog (Wright et al. 2010), 2MASS All Sky Survey (Cutri et al. 2003), IRAS point-source catalog (Helou & Walker 1988), and the historical variability of the *V*-band photometry as compiled in Ellerbroek et al. (2014). We also constrained these models using the surface brightness profiles along the major axis of our HiCIAO *H*-band scattered imagery (Figure 4).

We explored the parameter space of our models using a $\chi^2$ minimization scheme. Namely, we calculated the $\chi^2$ for the
<table><caption>**Table 2**
List of Key Best-fit Model Parameters and Estimates of the Upper and Lower Bounds of the Parameter</caption><thead><tr><th>Parameter (Units)</th><th>Best-fit Model</th><th>Lower Bound</th><th>Upper Bound</th></tr></thead><tbody><tr><td>Star Temperature (K)</td><td>9250</td><td>...</td><td>...</td></tr><tr><td>Star Radius (R<sub>☉</sub>)</td><td>1.4</td><td>1.2</td><td>1.6</td></tr><tr><td>Disk Mass (M<sub>☉</sub>)<sup>a</sup></td><td>0.05</td><td>...</td><td>...</td></tr><tr><td>Fraction of Mass in Large Grain Disk</td><td>0.9</td><td>0.8</td><td>0.95</td></tr><tr><td>Inner Gap Radius (au)</td><td>29</td><td>20</td><td>32</td></tr><tr><td>Outer Gap Radius (au)</td><td>59</td><td>55</td><td>62</td></tr><tr><td>Large Grain Disk Minimum Radius (R<sub>sub</sub>)<sup>b</sup></td><td>31.9</td><td>25</td><td>35</td></tr><tr><td>Large Grain Disk Maximum Radius (au)</td><td>250.1</td><td>...</td><td>...</td></tr><tr><td>Large Grain Disk Scale Height (R<sub>sub</sub>)<sup>b</sup></td><td>0.11</td><td>0.08</td><td>0.13</td></tr><tr><td>Large Grain Disk Radial Density Exponent</td><td>0.1</td><td>0.05</td><td>0.2</td></tr><tr><td>Large Grain Disk Scale Height Exponent</td><td>0.16</td><td>...</td><td>0.18</td></tr><tr><td>Small Grain Disk Minimum Radius (R<sub>sub</sub>)<sup>b</sup></td><td>1.22</td><td>1.0</td><td>1.5</td></tr><tr><td>Small Grain Disk Maximum Radius (au)</td><td>540.1</td><td>...</td><td>...</td></tr><tr><td>Small Grain Disk Scale Height (R<sub>sub</sub>)<sup>b</sup></td><td>0.11</td><td>0.08</td><td>0.13</td></tr><tr><td>Small Grain Disk Radial Density Exponent</td><td>0.05</td><td>...</td><td>...</td></tr><tr><td>Small Grain Disk Scale Height Exponent</td><td>1.25</td><td>...</td><td>...</td></tr><tr><td>Envelope Inner Radius (R<sub>sub</sub>)<sup>b</sup></td><td>0.41</td><td>...</td><td>...</td></tr><tr><td>Envelope Outer Radius (au)</td><td>2.38</td><td>...</td><td>...</td></tr><tr><td>Envelope Density (g cm<sup>-3</sup>)</td><td>4.0 × 10<sup>-17</sup></td><td>2.0 × 10<sup>-17</sup></td><td>6.0 × 10<sup>-17</sup></td></tr><tr><td>Accretion (M<sub>☉</sub>)</td><td>6.0 × 10<sup>-7</sup></td><td>...</td><td>...</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Disk mass value includes dust and gas. We assumed the gas to dust ratio is 100.

<sup>b</sup> R<sub>sub</sub> is the sublimation radius with 1 R<sub>sub</sub> = 0.36 au.

SED fit, the surface brightness along the major axis, and the minor axis offset, and added these values in quadrature to find the total $\chi^2$ value. Because some of the SED data were not contemporaneous, we also calculated a separate $\chi^2$ value that only incorporated comparisons of contemporaneously obtained data to the model. We began the iterative process with model runs of 5 million photons in order to find the best-fit SED to the SpeX and BASS spectra. Next, we increased the number of photons in each run to 50 million photons to obtain higher-quality model H-band images, and convolved the model image with the PSF of the H-band image. We explored parameter space to produce the best-fit $\chi^2$ value between model surface brightness along the major axis and the minor axis offset to the observed imagery. After finding the best chi-squared fit model image, we iteratively switched between the SED and the model until we found a model that optimized the combined chi-squared value, resulting in our best-fit model. We then reran this best-fit model using $10^9$ photons to produce the model SED and imagery used in all of our figures. We remind readers that MCRT models, like HOCHUNK3D, that employ a large family of parameters suffer from parameter degeneracy; thus, our best-fit model is not unique (Dong et al. 2012).

Table 2 lists the main parameters utilized in our best-fit model, and Figure 7 details the temperature and density profile of the disk in this model. Figures 8 and 9 show the SED and radial surface brightness profile along the disk major axis of our best-fit model as compared with our observations. We remark that our best-fit model parameters are generally similar to those previously reported in the literature. For example, our disk mass of 0.05 M<sub>☉</sub> (Table 2) is similar to that measured by Qi et al. (2011) (0.089 M<sub>☉</sub>) and Isella et al. (2007) (0.12 M<sub>☉</sub>).

Our best-fit model SED generally matches well with the contemporaneous spectroscopy and historical observations from optical to radio wavelengths (Figure 9). Since the optical

flux has been shown to be highly variable and we do not have contemporaneous optical photometry or spectroscopy, we do not know whether the modest model overestimation of the optical flux simply reflects that the star was at a high flux state in 2011. Additionally, our model reproduces the on-source extinction value of $A_V = 0.5$ mag from Ellerbroek et al. (2014) with a value of $A_V = 0.46$ mag. We note that the observed versus model imagery comparison matches well along the NW side of the disk (right-hand side of Figure 8), while the model imagery is marginally too narrow along the SE side of the disk (left-hand side of Figure 8). This could be due to slight geometrical variations in the wall of the disk, causing the illumination of the SE side of the disk to be broader. We provide a full comparison of the observed H-band PI imagery and model imagery in Figure 10. Our model imagery reveals little scattered light beyond the bright ring and little to no scattered light within the gap of the disk, which matches the observed PI and $Q_\phi$ images.

## 4. Analysis of SCExAO/CHARIS High-contrast Near-infrared Spectroscopy

### 4.1. Methodology: Disk and Planet Forward-modeling

Although none of the protoplanets/candidates reported from Keck/NIRC2 or ALMA are visible in our data, great care is needed to properly interpret these non-detections and their implications. For example, like HD 163296, HD 100546 has multiple imaged protoplanet candidates embedded in a bright, structured protoplanetary disk (Quanz et al. 2013; Currie et al. 2015). Follow-up claims of a spurious detection/non-detection of candidates around HD 100546 were faulty, as shown in Currie et al. (2017), in large part due to (1) incorrect spectrophotometric calibration and (2) a lack of forward-modeling of planet and disk signals.
<figure><img src="image_13.png" /><figcaption>**Figure 7.** Top row of panels present temperature profiles for three regions of our MCRT disk model. The bottom row of panels present the density profiles for these same three regions of the disk model.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 8.** Major axis crosscut of our 2011 H-band imagery data (PI image) compared with the best-fit model (red-dashed line). The vertical dashed lines represent the inner working angle of 0.′′28.</figcaption></figure>

Contemporaneous near-IR spectra of HD 163296 allowed us to spectrophotometrically calibrate CHARIS data cubes (see Section 2.1). To properly understand our non-detections and derive upper limits at the candidates' locations, we then performed forward-modeling of our images, investigating the reduction of the total source signal and the biasing of its spatial intensity distribution due to processing. This annealing results

<figure><img src="image_15.png" /><figcaption>**Figure 9.** The observed SED of HD 163296 is shown along with our best-fit model SED (black line). The SpeX 2011 (green line) and BASS 2011 (red line) data are from this work, as described in Section 2. The blue circles represent data from the AllWISE catalog (Wright et al. 2010), the green circles are from the 2MASS All Sky Survey (Cutri et al. 2003), and the purple circles are from IRAS point-source catalog (Helou & Walker 1988). The gray circles depict V-band photometry and represent the historical minimum, $1\sigma$ below median flux, median flux, and $1\sigma$ above the median flux as reported by Ellerbroek et al. (2014).</figcaption></figure>

from self-subtraction of the source by itself and over-subtraction of the disk in ADI and SDI. Our method follows that outlined in Currie et al. (2018a), where we save the A-LOCI coefficients $\alpha$ and model the disk and planet signals as
<figure><img src="image_16.png" /><figcaption>**Figure 10.** Our 2011 *H*-band polarized scattered light image (left panel), the best-fit model *PI H*-band scattered light image (middle panel), and the difference between the observed and model *PI* image (right panel) are shown. All three panels are displayed on the same linear scale and same spatial scale, and rotated such that north is up and east is left. The inner working angle is masked out with a white circle. Note that the *PI* image was binned to match the pixel scale of the model for the difference image.</figcaption></figure>

introducing a linear perturbation of value $\beta$, which provides an additional source of annealing (see also Brandt et al. 2013; Pueyo 2016). We focus on the 2018 May data due to its higher quality.

First, we explored the effect of the disk on the non-detections of planetary companions, using forward-modeling to determine its annealing due to processing and its effect on any point sources located exterior, like the proposed companions from Guidi et al. (2018) and Teague et al. (2018). We started with the best-fit scattered light disk model described in Section 3.2, which is drawn from our *H*-band scattered light imagery with Subaru/HiCIAO. We produced total intensity (not scattered light) images in *J*, *H*, and *K* passbands and interpolated the model images onto the CHARIS wavelength array and pixel scale. The model disk is slightly bluer than the combined light of the star+disk, with intrinsic colors of *J*−*H*, *H*−*K* of ~0.35 and ~0.35. Note that the model was constructed based on a single passband (*H*-band); thus, the model may not constrain the true color of the disk. The disk contrast with respect to the star on the forward-scattering side is typically $\Delta M = M_{\text{disk/arcsec}^2} - M_* \approx 3.5–4$. The visible trace of a disk may differ in total intensity versus scattered light. Therefore, we slightly adjusted the model parameters to provide a better match to the forward-modeled disk image, specifically increasing the semimajor axis by 5%.

Second, we verified that an object consistent with the 6–7 $M_J$ candidate from Guidi et al. would be detected in our data. We used standard hot-start evolutionary models from Baraffe et al. (2003), adopting a planet age equal to the system age (5 Myr). This approach is intermediate between possible extremes that would yield higher and lower luminosities for a given planet mass. While we assume a planet age of 5 Myr when estimating mass limits, the age of a superjovian planet is likely much younger than that of the host star (Currie et al. 2013). This is especially true for protoplanets, which are nearing the end of their formation and thus much closer to $t \lesssim 1$ Myr for any evolutionary model, where the planet luminosity is maximum. The inferred limits adopting would be then substantially lower than those we report. Conversely, we could adopt planet mass limits using the “cold-start” evolutionary models (e.g., Marley et al. 2007). However, recent literature casts serious doubt on the validity of the cold-start model formalism, which relies on specific assumptions about the entropy of accreted material. As shown by Berardo et al. (2017), classic cold-start conditions are

extremely difficult to reach because the protoplanet will be substantially heated by the accretion shock, which will increase its entropy, resulting in hot-start-like initial conditions. Furthermore, imaged planets for which we have derived dynamical masses—$\beta$ Pic b, HR 8799 bcde (Lagrange et al. 2010; Marois et al. 2010; Currie et al. 2011; Snellen & Brown 2018; Wang et al. 2018; Dupuy et al. 2019)—are inconsistent with a cold-start evolutionary model. At the candidate’s location in each data cube, we injected a planet whose temperature matches that expected for a 4 $M_J$, 5 Myr planet according to these models. Although such a planet is predicted to be near the L/T dwarf transition ($T_{\text{eff}} \sim 1300$ K), we assume a (cloudier) L dwarf spectrum drawn from the Bonnefoy et al. (2014) library, because annealing due to SDI will be stronger for such a spectrum. Integrated over the CHARIS wavelength array, the broadband contrast of this planet with respect to HD 163296 is $\sim 8 \times 10^{-6}$, about 2.5–3.5 times as high as the predicted contrast for the Guidi et al. companion using a cloudy planet atmosphere from Currie et al. (2011).

Finally, our forward-modeling calculation allowed us to compute radially averaged, throughput-corrected broadband contrast curves. As with our fake planet injection test, we used the Baraffe et al. (2003) models to map between planet mass and temperature. To map between temperature and spectrum, we further used atmosphere models drawn from Burrows et al. (1997), adopting cloud prescriptions that provide reasonable fits to near-IR photometry for HR 8799 bcde and ROXs 42Bb, whose temperature, gravity, and masses cover most of our range (Currie et al. 2011, 2014; Madhusudhan et al. 2011; T. Currie et al. 2018, in preparation).

## 4.2. Results: Limits on Planets

Figure 11 shows the wavelength-collapsed image of the input disk (left panel) and output image after forward-modeling the disk through ADI and SDI (right panel). While the disk in total intensity is more forward-scattering than the model based on polarimetry would predict, and its brightness is ~30% higher, the model otherwise reproduces the CHARIS data and is sufficient for investigating the impact of self-subtraction on the forward-scattering side. The proposed candidate from Guidi et al. (2018) lies exterior to the main trace of the disk (cyan cross). After processing, the candidate’s location is free of
<figure><img src="image_17.png" /><figcaption>**Figure 11.** (Left) Broadband image of the best-fit synthetic disk model derived from polarimetry interpolated onto the CHARIS pixel scale and wavelength array and (right) forward-model of the disk after propagating its signal through ADI and SDI. The location of the proposed protoplanet candidate from Guidi et al. (2018) lies well exterior to the azimuthal and radial self-subtraction footprints in the forward-modeled disk. The images have been smoothed with a top-hat filter to more clearly reveal the trace of the disk: localized emission exterior to the disk is an artifact of this smoothing.</figcaption></figure>

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

<figure><img src="image_19.png" /><figcaption>**Figure 12.** (Left) 2018 May broadband image with a 4 $M_J$, 5 Myr old planet injected into our observing sequence at the location of the candidate from Guidi et al. (2018) ($\Delta F \sim 8 \times 10^{-6}$) and propagating its signal through ADI and SDI. Even with signal from the disk contributing to an estimate of the noise, the injected companion is detected at S/N ~5. (Right) Broadband contrast curve for the 2018 May and 2018 June data compared with broadband contrasts for 2–10 $M_J$ planets assuming the Burrows atmosphere models. The 5$\sigma$ contrast at $0''49$ is in agreement with expectations based on our injected 4 $M_J$ planet in the left panel. The contrast for a 1 $M_J$ companion lies off the graph at $\Delta F \sim 3.7 \times 10^{-7}$.</figcaption></figure>

negative self-subtraction footprints. Inspection of the individual data cubes containing the disk model processed through ADI and SDI likewise show a flat background. At wider separations overlapping with the proposed candidate from Teague et al. (2018), the disk likewise leaves negligible residual effects.

As shown in Figure 12 (left panel), a 6–7 $M_J$ candidate similar to the one proposed in Guidi et al. should have been detected in our data. The fainter, even lower-mass (4 $M_J$) candidate injected into our data is clearly visible. While its S/N is formally ~4.8, our inclusion of disk signal contributions leads our estimate of the noise to be conservative. A planet corresponding to the Guidi et al. candidate ($\Delta F \sim 2.5 \times 10^{-5}$) would be even more decisively detected (S/N ~15).

Broadband contrast limits in the right-hand panel of Figure 12 provide stringent limits on protoplanets covering the range probed with Keck/NIRC2 and ALMA. At $\rho \sim 0''49$, the azimuthally averaged 5$\sigma$ contrast limit is ~8.5 $\times 10^{-6}$, in agreement with our expectations from the fake planet injection. If the Guidi et al. companion is real, it would then have to be redder than $H-L_p \sim 3.5$ to escape detection: redder than all directly imaged planets except for the extreme L/T transition object HD 95086 b (De Rosa et al. 2016). Over the separations

just interior or close to the visible trace of the disk and comparable to the separation of the Guidi et al. companion— $\rho \sim 0''4$ (0''7) along the minor (major) axis—we can exclude planets with masses of 2–5 $M_J$, assuming standard hot-start evolutionary models. The CHARIS field encloses the possible location of the innermost companion proposed by Teague et al. (2018), which would lie at a projected separation of $r_{\text{proj}} \sim 83$ au ($\rho \sim 0''82$) along the major axis, or $r_{\text{proj}} \sim 40$ au ($\rho \sim 0''4$) along the minor axis. At these locations, our data rule out planets more massive than 5 $M_J$ and ~1.5 $M_J$, respectively. If located along the minor axis, the outermost proposed companion from Teague et al. (2018) would be at $\rho \sim 0''65$ with a mass less than ~2 $M_J$, according to our data.

# 5. Discussion

## 5.1. Previous Optical–IR Disk Imaging

HD 163296 has been observed numerous times across optical–IR bandpasses. Here we briefly summarize some of the major results of those investigations, to compare and contrast with our new imagery.
Space-based optical imagery has been obtained in both white light (HST/STIS; Grady et al. 2000) and broadband filters (HST/ACS; Wisniewski et al. 2008), tracing the disk out to 4"/4 (447 au) and detecting HH knots. Comparison of these data revealed evidence for significant variation (~1 mag) in the disk surface brightness, changes in the number of disk ansae visible over time, and changes in the relative brightness of features located in the NW and SE disk regions (Wisniewski et al. 2008). Unfortunately, none of these optical observations fully overlapped in wavelength coverage.

Ground-based AO imagery of the system can be generally summarized into 3 categories. First, a subset of observations clearly reveal the detection of the disk in scattered light, but the presence of residual AO speckle noise in the disk vicinity prevents a robust characterization of the surface brightness or detailed morphological structure of the disks (e.g., 2012 H-band imagery Garufi et al. 2014; 2014 Ks-band imagery Garufi et al. 2017). Second, a subset of observations (e.g., 2012 Ks-band imagery; Garufi et al. 2014) reveal the detection of an inclined ring structure extending out to 1"/03 (103 au), where the intensity of scattered light is strongest along the major axis of the disk and is symmetrical about both sides of the disk major axis (NW side and SE side). Third, a subset of observations (e.g., 2014 J-band, Monnier et al. 2017; 2016 H-band imagery, Muro-Arena et al. 2018) reveal clear evidence of this same inclined ring structure whose flux is both azimuthally asymmetric and not the strongest along the major axis. In particular, the NW side of the major axis is brighter than the SE side of the disk in J-band GPI observations (see Figure 2, Monnier et al. 2017), and the maximum flux from the disk is north of the major axis peaking on the NW side of the disk in these data. The 2016 H-band VLT/SPHERE observations (Muro-Arena et al. 2018) also exhibit strong azimuthal asymmetry, with the NW side of the disk along the major axis exhibiting 2.7× more scattered light than the SE side of the disk along the major axis. Muro-Arena et al. (2018) used 3D radiative transfer modeling to suggest that this strong azimuthal asymmetry could be reproduced by including an inner disk component that was misaligned by 1° compared with the outer disk.

## 5.2. Evidence for Time-dependent Azimuthal Asymmetry

Our 2011 epoch H-band imagery is consistent with the second category of disk appearance we discussed in Section 5.1. Namely, we observe a broken ring structure in H-band scattered light whose flux peaks along the major axis and exhibits clear symmetry between the NW and SE side of the disk. Our 2011 epoch H-band data are thus clearly different from the 2016 epoch VLT/SPHERE H-band data, which show a 2.7× asymmetry between the NW and SE side of the disk (Muro-Arena et al. 2018).

To illustrate these differences, we scaled the peak flux along the major axis of the 2016 VLT/SPHERE data and present these data as dashed horizontal lines in our Figure 4. The 2.7× asymmetry about the major axis observed in the 2016 VLT/SPHERE data is clearly outside of the 3σ errors of our 2011 data. This obvious difference is also seen by comparing Figure 1 of Muro-Arena et al. (2018) with Figure 1 in this paper. We note that neither data set exhibits evidence of large-scale gradients in its $U_φ$ component, indicating that systematic artifacts are not the cause of this phenomenon. We suggest that this is clear evidence that the system exhibits large changes in

the appearance of its scattered light disk as seen in multi-epoch observations obtained with the same filter, and supports previous suggestions of this phenomenon as deduced from multi-epoch observations from similar, albeit not the same, filters (Wisniewski et al. 2008).

There are several mechanisms that could cause an azimuthal asymmetry of scattered light, including an asymmetrical distribution of dust (Muro-Arena et al. 2018), an inclined inner disk shadowing the outer disk (Muro-Arena et al. 2018), a warped inner disk structure shadowing the outer disk (Sitko et al. 2008), or dust ejected above the midplane of the disk that shadows the outer disk (Ellerbroek et al. 2014).

Muro-Arena et al. (2018) suggested that an asymmetric distribution of dust in the system was unlikely, as no asymmetry was observed with ALMA (Isella et al. 2016). Muro-Arena et al. (2018) were able to replicate the azimuthal asymmetry they observed in their scattered light imagery by inclining the inner disk by 1° compared with the outer disk, which is consistent with previous near-IR interferometric observations (Tannirkulam et al. 2008; Lazareff et al. 2017; Setterholm et al. 2018). However, our 2011 epoch data reveal the presence of no azimuthal asymmetry along the major axis in the same filter bandpass as the 2016 SPHERE observations. An inclined inner disk is unlikely to precess significantly over a 5 yr time frame; hence, an inclined inner disk alone is unlikely to produce the observed significant azimuthal variations in the scattered light disk. Moreover, we have shown that we can reproduce the basic properties of both our contemporaneously obtained near-IR SED and H-band imagery with a model that does not include an inner inclined disk. Thus, while the system could plausibly host an inclined disk, we suggest that this feature is unlikely to be responsible for producing the time-dependent azimuthal variations in the outer scattered light disk of the system.

We consider several other mechanisms that could explain the change in disk surface brightness seen in the system. First, a warped inner disk structure, such as a puffed-up inner disk wall (Turner et al. 2014), could be shadowing the outer disk (Sitko et al. 2008). If this disk warp were to dissipate or rotate azimuthally within a 2–3 yr timescale, this could cause a change in illumination of the outer disk similar to that observed between the 2011 and 2016 epoch H-band data sets. Dynamical simulations are needed to determine whether a substantial change in the appearance of a warped disk could occur on this short of a timescale and lead to the amplitude of variable disk illumination observed.

Second, this phenomenon could be caused by dust ejected above the midplane of the disk, which partially shadows the outer disk, as proposed by Ellerbroek et al. (2014). These dust “clouds” could differentially obstruct the illumination of the outer disk while they are between the star and the outer disk, as shown in Figure 13. We do have IR spectra that were obtained at a similar epoch to both our 2011 HICAO data and the 2016 SPHERE data. The contemporaneous IR spectra cannot constrain the possible asymmetric nature of the dust clouds, but can constrain the total amount of dust in the disk wind when compared with our MCRT models. As shown in Figure 2, while both have the same flux around 0.9 µm, the 2011 epoch IR spectrum is brighter (~0.5 mag at K') around 2 µm than the 2016 epoch IR spectrum. We remark that we can best reproduce the 2016 SED in our model by adopting a ~2× lower envelope density, e.g., 9.0 × 10⁻¹⁸ g cm⁻³, which
<figure><img src="image_20.png" /><figcaption>**Figure 13.** Diagram of the disk wind model. (A) shows the disk wind, which is asymmetric and shadows the SE portion of the disk. (B) shows a symmetric disk wind where the both sides of the disk are equally illuminated. The left-hand side of the diagram shows the outer portion of the disk, where the right-hand side of the diagram shows a zoomed-in version of the disk. The outer disk as been rotated and inclined to match the observed orientation of HD 163296 shown in Figure 1.</figcaption></figure>

corresponds to a lower circumstellar extinction in 2016 of $A_V = 0.1$ mag. We predict that the 2016 epoch should be 0.4 mag brighter in the V-band compared with the 2011 epoch data, similar to the optical variability found by Ellerbroek et al. (2014). Because the interstellar extinction is consistent with $AV = 0$ mag, the observed reddening most likely originates from the system. Thus an asymmetric disk wind launching dust clouds can explain both the variable illumination of the outer disk (Figure 13) and the reddening optical variability observed by our disk wind models and Ellerbroek et al. (2014).

If the system does have an inclined inner disk as suggested by Muro-Arena et al. (2018) that during some epochs produces non-axisymmetric illumination of the outer disk (e.g., NW side brighter than SE side; 1998 HST/STIS, Grady et al. 2000, 2014 J-band, Monnier et al. 2017; 2016 H-band, Muro-Arena et al. 2018), the spatial distribution of any dust clouds elevated by a disk wind must also be non-axisymmetric to produce the observed epochs of axi-symmetric illumination of the outer disk (e.g., as seen in 2012 Ks-band imagery, Garufi et al. 2014; 2011 H-band, this study) and the sole-epoch of observed non-axisymmetric illumination with the SE side of the disk brighter than the NW side (2004 HST/ACS, Wisniewski et al. 2008). Future observations that simultaneously observe quiescent and wind events with contemporaneous optical and IR photometry and coronagraphic imagery could help to test whether shadowing by dust clouds could explain the observed behavior

of the inner and outer disk of the system, and better parameterize the azimuthal distribution of such dust clouds.

### 5.3. Model

We were able to reproduce the basic properties of our contemporaneous near-IR spectra and scattered light H-band imaging with a 3D MCRT disk model, which approximated the features of a disk wind via an envelope. As seen in Figure 9, our model SED is consistent with the highest observed V-band flux that was reported by Ellerbroek et al. (2014), but we caution that the robustness of this agreement is uncertain because we do not have contemporaneous optical photometry.

Muro-Arena et al. (2018) also performed MCRT modeling of HD 163296, and compared their models to the ALMA dust continuum image from Isella et al. (2016), their own VLT/SPHERE image, and historical photometry and spectroscopy. They modeled all three gaps that were observed in the ALMA continuum image and introduced an inclined disk to explain the asymmetric scattered light flux observed with the VLT/SPHERE image as noted above in Section 5.2. Their model images and SED are well matched to their observed images and historical photometry and spectroscopy. While they do not employ a disk wind model as we did (Section 3.2, see Figure 13), their model does not have a clear mechanism to explain the time-dependent azimuthal asymmetries seen in near-IR scattered light images (Section 5.2) or the optical-IR
photometric and spectroscopic variability that has been
observed (Sitko et al. 2008; Ellerbroek et al. 2014). We
caution that the inability of an inclined disk by itself to explain
the observed time-dependent azimuthal asymmetries in scat-
tered light does not exclude the possibility that the system does
in fact have an inclined disk. Due to limitations with
HOCHUNK3D, we leave applying our disk wind model to
the archival images and SEDs to future work.

5.4. *Scattered Light Features along the Minor Axis*

We note that a deficit of scattered light is seen in the near
side of our disk imagery at a PA of 30°, in both our binned
imagery (Figure 6) and our unbinned PI and Q<sub>rot</sub> images. We
caution that while this feature could be real, it is not uncommon
to observe depolarization along the minor axis due to the
residual presence of an uncorrected polarized halo. Interest-
ingly, this feature coincides with the disk brightness increase
observed with the Keck/NIRC2 L'-band vortex coronagraph by
Guidi et al. (2018) and is located at a similar position angle,
albeit closer to the host star, as the purported candidate
planetary mass object reported by Guidi et al. (2018). As noted
by Guidi et al. (2018), this disk feature is located where
forward scattering should be significant. If the feature we
observe at the similar disk position is astrophysical, the
decreased amplitude of the feature in polarized intensity
suggests that it could be polarized less than its neighboring
disk material.

5.5. *Limits on Protoplanets Orbiting HD 163296*

Our data improve the detection limits for protoplanets in
thermal emission around HD 163296 compared with Keck/
NIRC2 data from Guidi et al. (2018): from 5–7 *M*<sub>J</sub> to now
2–5 *M*<sub>J</sub>, assuming standard hot-start evolutionary models, near
the projected trace of the disk. At wider separations covering
the possible locations of the inner proposed candidate from
Teague et al. (2018) (*r*<sub>proj</sub> ~ 83 au/ρ ~ 0.82), the limits have
now improved from 4.5 *M*<sub>J</sub> to 1.5 *M*<sub>J</sub>, the latter of which is just
slightly higher than the predicted mass of the companion
(1 *M*<sub>J</sub>). Limits for the outer Teague et al. candidate along the
minor axis are likewise just slightly higher than the predicted
mass (a limit of 2 *M*<sub>J</sub> versus a predicted 1.3 *M*<sub>J</sub>). Thus, at least
for now, the ALMA-predicted protoplanet candidates are
consistent with direct imaging constraints.

Our data appear to rule out the proposed, marginally
significant candidate identified from thermal IR data in Guidi
et al. (2018). Using standard assumptions for planet atmo-
spheres, our forward-modeling demonstrates we could have
detected an even fainter planet at the location of the proposed
candidate. For an assumed age of 5 Myr and hot-start
evolutionary models, the candidate is predicted to be 6–7 *M*<sub>J</sub>,
while our radially averaged contrast limits are significantly
lower (~4–5 *M*<sub>J</sub>).<sup>43</sup>

The simplest explanation for our conflicting results is that the
NIRC2 candidate is instead residual, partially subtracted
speckle noise or partially subtracted disk emission left over

from processing. Figure 1 of Guidi et al. (2018) shows multiple
emission peaks with a similar or slightly smaller spatial scale as
the candidate (e.g., at the 2, 6, 7, and 8 o'clock positions just
exterior to the masked region). An even brighter, seemingly
point-source-like peak at nearly the same position angle in
these data appears to be an artificially enhanced region of the
disk, which could have been mistaken for a point in shallower
and/or higher background data. Convolving the image with a
Gaussian kernel may further accentuate the point-source-like
appearance of these features.<fnref n="44" /> The position of the candidate
also coincides with the minor axis of a second ring of emission
detected with ALMA. Forward-modeling as performed in
Currie et al. (2015) could better clarify whether the candidate's
morphology is consistent with an annealed point source or
residual disk emission.

Alternatively, the candidate could be extremely red/under-
luminous in the near-IR and thus difficult to detect. If
embedded in the disk, it would be preferentially extincted in
the near-IR compared with the thermal infrared, as has been
proposed for HD 100546 b (Currie et al. 2015; Quanz et al.
2015). It could also retain an extremely dusty/cloudy
atmosphere characteristic of some young exoplanets near the
L/T transition (Currie et al. 2011; De Rosa et al. 2016), making
it appear “underluminous” in the near-infrared. Follow-up
thermal infrared imaging at L<sub>p</sub> or M<sub>p</sub> could provide a more
decisive probe of these possibilities.

6. Conclusions

We report *H*-band polarimetric imagery of the HD 163296
system along with contemporaneous infrared spectra observa-
tions and near-IR extreme AO imaging in total intensity.
We find:

1. Our 2011 *H*-band polarimetric imagery resolve a broken ring structure surrounding HD 163296 that peaks at a distance along the major axis of 0.65 (66 au) and extends out to 0.98 (100 au) along the major axis. Our non-detection of the inner disk component is driven by our inner working angle (0.3, 30.5 au), and does not conflict with the detection of this component by Monnier et al. (2017).

2. Our 2011 epoch *H*-band imagery exhibits clear axisym-
metry, with the NW and SE side of the disk exhibiting
similar intensities. Our 2011 epoch *H*-band data are thus
clearly different from the 2016 epoch *H*-band data from
VLT/SPHERE reported by Muro-Arena et al. (2018),
which exhibit a strong 2.7× asymmetry between the NW
and SE side of the disk. These results indicate the
presence of time-variable, non-azimuthally symmetric
illumination of the outer disk.

3. We were able to reproduce the basic properties of our contemporaneous near-IR spectra and spatially resolved *H*-band polarimetric imagery of the HD 163296 disk with a 3D MCRT disk model that approximated the features of

<footnote n="43">43 Note that any new age estimates for HD 163296 drawn from its Gaia-revised distance do not change our results. Comparisons to some isochrones may imply an older age (e.g., 7.6 ± 1.1 Myr; Vioque et al. 2018). However, other (e.g., the MIST and PARSEC) isochrones imply ages comparable to or just slightly greater than 5 Myr. These differences do not change the fact that the proposed HD 163296 companion should have been detected in our data under standard assumptions for planet atmospheres.</footnote>

<footnote n="44">⁴⁴ The large spatial scale of the residuals may also be traced to the PSF subtraction method used, which leverages on the Karhunen–Loéve Image Projection (KLIP) algorithm with few KL modes retained (Soummer et al. 2012). Compared with standard implementations of A-LOCI, KLIP with few KL modes retained may yield larger spatial-scale residuals. This is especially true for KLIP implementations performing PSF subtraction in full annuli as in Guidi et al. instead of smaller wedge-shaped annular regions, because the subtraction is less local, in addition to constructing a low-rank approximation of the data set’s covariance matrix.</footnote>
a disk wind via an envelope and did not specifically require an inclined inner disk component. We suggest that, while the system could plausibly host an inclined disk as suggested by Muro-Arena et al. (2018), such a component is unlikely to be responsible for producing the observed time-dependent azimuthal variations in the outer scattered light disk of the system. We speculate that a variable, non-axisymmetric distribution of dust clouds elevated by a disk wind could produce the diversity of morphological appearances of the outer disk now reported in the literature for this system.

4. While our 2018 epoch SCExAO/CHARIS observations easily recover the disk, they fail to recover the candidate $6-7 M_J$ protoplanet identified from Keck/NIRC2 data (Guidi et al. 2018). The Keck/NIRC2 detection is likely a residual speckle or a partially subtracted piece of the disk; alternatively, this object could be a heavily embedded or particularly red/cloudy object only identifiable in the thermal infrared.

5. Assuming hot-start evolutionary models and a system age of 5 Myr, our SCExAO/CHARIS detection limits for protoplanets in thermal emission around HD 163296 near the projected trace of the disk are 2–5 $M_J$. At wider separations, covering the possible locations of the inner proposed candidate from Teague et al. (2018) ($r_{\text{proj}} \sim 83 \text{ au}/\rho \sim 0.82$), our data lower the mass limit for detections from 4.5 $M_J$ to 1.5 $M_J$, which is still slightly higher than the predicted mass of the companion ($1 M_J$). Limits for the outer Teague et al. candidate along the minor axis are likewise just slightly higher than the predicted mass (a limit of 2 $M_J$ versus a predicted 1.3 $M_J$). The ALMA-predicted protoplanet candidates are currently still consistent with direct imaging constraints.

We acknowledge support from the NASA XRP program via NNX-17AF88G. The authors recognize and acknowledge the significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

## ORCID iDs

Evan A. Rich https://orcid.org/0000-0002-1779-8181

John P. Wisniewski https://orcid.org/0000-0001-9209-1808

Thayne Currie https://orcid.org/0000-0002-7405-3119

Misato Fukagawa https://orcid.org/0000-0003-3500-2455

Michael L. Sitko https://orcid.org/0000-0003-1799-1755

Timothy D. Brandt https://orcid.org/0000-0003-2630-8073

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

Ruobing Dong https://orcid.org/0000-0001-9290-7846

Olivier Guyon https://orcid.org/0000-0002-1097-9908

Yutaka Hayano https://orcid.org/0000-0003-4937-4233

Klaus W. Hodapp https://orcid.org/0000-0003-0786-2140

Markus Janson https://orcid.org/0000-0001-8345-593X

Ryo Kandori https://orcid.org/0000-0003-2610-6367

Gillian R. Knapp https://orcid.org/0000-0002-9259-1164

Tomoyuki Kudo https://orcid.org/0000-0002-9294-1793

Jungmi Kwon https://orcid.org/0000-0003-2815-7774

Frantz Martinache https://orcid.org/0000-0003-1180-4138

Taro Matsuo https://orcid.org/0000-0001-7484-1222

Satoshi Mayama https://orcid.org/0000-0002-3424-6266

Shoken Miyama https://orcid.org/0000-0001-5017-180X

Amaya Moro-Martin https://orcid.org/0000-0001-9504-8426

Tae-Soo Pyo https://orcid.org/0000-0002-3273-0804

Michihiro Takami https://orcid.org/0000-0001-9248-7546

Hiroshi Terada https://orcid.org/0000-0002-7914-6779

Christian Thalmann https://orcid.org/0000-0002-1664-2177

Taichi Uyama https://orcid.org/0000-0002-6879-3030

Kevin R. Wagner https://orcid.org/0000-0002-4309-6343

Makoto Watanabe https://orcid.org/0000-0002-3656-4081

Tomonori Usuda https://orcid.org/0000-0001-9855-0163

Motohide Tamura https://orcid.org/0000-0002-6510-0681

## References

Bally, J. 2016, ARA&A, 54, 491

Baraffe, I., Chabrier, G., Barman, T. S., et al. 2003, A&A, 402, 701

Berardo, D., Cumminng, A., & Marleau, G.-D. 2017, ApJ, 834, 149

Bonnefoy, M., Chauvin, G., Lagrange, A.-M., et al. 2014, A&A, 562, 127

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

Brandt, T. D., Rizzo, M., & Groff, T. 2017, JATIS, 3, 8002

Burrows, A., Marley, M., Hubbard, W. B., et al. 1997, ApJ, 491, 856

Currie, T., Brandt, T. D., Uyama, T., et al. 2018a, AJ, 156, 291

Currie, T., Brittain, S., Grady, C., et al. 2017, RNAAS, 1, 40

Currie, T., Burrows, A., Itoh, Y., et al. 2011, ApJ, 729, 128

Currie, T., Burrows, A., Madhusudhan, N., et al. 2013, ApJ, 776, 15

Currie, T., Cloutier, R., Brittain, S., et al. 2015, ApJL, 814, L7

Currie, T., Daemgen, S., Debes, J., et al. 2014, ApJL, 780, L30

Currie, T., Debes, J., Rodigas, T., et al. 2012, ApJL, 760, L32

Currie, T., Kasdin, N. J., Groff, T., et al. 2018b, PASP, 130, 044505

Cushing, M. C., Vacca, W. D., & Rayner, J. T. 2004, PASP, 116, 362

Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, yCat, 2246

De Rosa, R., Rameau, J., Patience, J., et al. 2016, ApJ, 824, 121

Dong, R., Rafikov, R., Zhu, Z., et al. 2012, ApJ, 750, 161

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

Ellerbroek, L. E., Podio, L., Dougados, C., et al. 2014, A&A, 563, A87

Fernandes, R. B., Long, Z. C., Pikhartova, M., et al. 2018, ApJ, 856, 103

Fukagawa, M., Tamura, M., Itoh, Y., et al. 2010, PASJ, 62, 347

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

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

Garufi, A., Meeus, G., Benisty, M., et al. 2017, A&A, 603, A21

Garufi, A., Quanz, S. P., Schmid, H. M., et al. 2014, A&A, 568, A40

Goebel, S., Currie, T., Guyon, O., et al. 2018a, AJ, submitted

Goebel, S., Currie, T., Guyon, O., et al. 2018b, AJ, 156, 279

Grady, C. A., Devine, D., Woodgate, B., et al. 2000, ApJ, 544, 895

Groff, T. D., Kasdin, N. J., Limbach, M., et al. 2015, Proc. SPIE, 9605, 1

Guidi, G., Ruane, G., Williams, J. P., et al. 2018, MNRAS, 479, 1505

Guidi, G., Tazzari, M., Testi, L., et al. 2016, A&A, 588, A112

Halir, R., & Flusser, J. 1998, in Proc. Sixth Int. Conf. Computer Graphics and Visualization (Plzen: Univ. West Bohemia), 125

Hammel, B., & Sullivan-Molina, N. 2019, Bdhammel/least-squares-ellipse-fitting: Initial release, v1.0, Zenodo, doi:10.5281/zenodo.2578663

Hashimoto, J., Tamura, M., Muto, T., et al. 2011, ApJL, 729, L17

Hayano, Y., Takami, H., Guyon, O., et al. 2008, Proc. SPIE, 7015, 701510

Hayano, Y., Takami, H., Oya, S., et al. 2010, Proc. SPIE, 7736, 77360N

Hein Bertelsen, R. P., Kamp, I., van der Plas, G., et al. 2016, MNRAS, 458, 1466

Helou, G., & Walker, D. W. 1988, Infrared Astronomical Satellite (IRAS) Catalogs and Atlases, Vol. 7 (Washington, DC: US Gov. Printing Office)

Herbig, G. H. 1950, ApJ, 111, 11

Herbig, G. H. 1960, ApJS, 4, 337

Hodapp, K. W., Suzuki, R., Tamura, M., et al. 2008, Proc. SPIE, 7014, 701419

Isella, A., Guidi, G., Testi, L., et al. 2016, PhRvL, 117, 251101

Isella, A., Testi, L., Natta, A., et al. 2007, A&A, 469, 213

Jovanovic, N., Guyon, O., Martinache, F., et al. 2015, ApJL, 813, L24

Jovanovic, N., Martinache, F., Guyon, O., et al. 2015, PASP, 127, 890

Kim, S.-H., Martin, P. G., & Hendry, P. D. 1994, ApJ, 422, 164

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

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

Lazareff, B., Berger, J.-P., Kluska, J., et al. 2017, A&A, 599, A85

Madhusudhan, N., Burrows, A., & Currie, T. 2011, ApJ, 737, 34

Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007, ApJ, 655, 541
Marois, C., Correia, C., Veran, J.-P., & Currie, T. 2014, in Proc. IAU Symp. 299, Exploring the Formation and Evolution of Planetary Systems, ed. M. Booth, B. C. Matthews, & J. R. Graham (Cambridge: Cambridge Univ. Press), 48

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

Marois, C., Zuckerman, B., Konopacky, Q., Macintosh, B., & Barman, T. 2010, Natur, 468, 1080

Mawet, D., Milli, J., Wahhaj, Z., et al. 2014, ApJ, 792, 97

Meeus, G., Waters, L. B. F. M., Bouwman, J., et al. 2001, A&A, 365, 476

Mendigutía, I., Brittain, S., Eiroa, C., et al. 2013, ApJ, 776, 44

Monnier, J. D., Harries, T. J., Aarnio, A., et al. 2017, ApJ, 838, 20

Montesinos, B., Eiroa, C., Mora, A., & Merín, B. 2009, A&A, 495, 901

Muro-Arena, G. A., Dominik, C., Waters, L. B. F. M., et al. 2018, A&A, 614, 24

Ohta, Y., Fukagawa, M., Sitko, M. L., et al. 2016, PASJ, 68, 53

Pinte, C., Price, D. P., Menard, F., et al. 2018, ApJL, 860, L13

Pueyo, L. 2016, AJ, 824, 114

Qi, C., D'Alessio, P., Öberg, K. I., et al. 2011, ApJ, 740, 84

Quanz, S., Amara, A., Meyer, M., et al. 2015, ApJ, 807, 64

Quanz, S., Amara, A., Meyer, M. R., et al. 2013, ApJL, 766, L1

Rayner, J. T., Toomey, D. W., Onaka, P. M., et al. 2003, PASP, 115, 362

Rosenfeld, K. A., Andrews, S. A., Hughes, M. A., Wilner, D. J., & Qi, C. 2013, ApJ, 774, 16

Schmid, H. M., Joos, F., & Tschan, D. 2006, A&A, 452, 657

Setterholm, B. R., Monnier, J. D., Davies, C. L., et al. 2018, ApJ, 869, 164

Sitko, M. L., Carpenter, W. J., Kimes, R. L., et al. 2008, ApJ, 678, 1070

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

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

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

Takami, M., Fu, G., Liu, H. B., et al. 2018, ApJ, 864, 20

Tamura, M. 2009, in AIP Conf. Proc. 1158, Exoplanets and Disks: Their Formation and Diversity, ed. T. Usuda, M. Tamura, & M. Ishii (Melville, NY: AIP), 11

Tannirkulam, A., Monnier, J. D., Millan-Gabet, R., et al. 2008, ApJ, 677, 51

Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018, ApJL, 860, 12

Turner, N. J., Benisty, M., Dullemond, C. P., & Hirose, S. 2014, ApJ, 780, 42

Vacca, W. D., Cushing, M. C., & Rayner, J. T. 2003, PASP, 115, 389

Vioque, M., Oudmaijer, R. D., Baines, D., Mendigutía, I., & Perez-Martinez, R. 2018, A&A, 620, 128

Wagner, K. R., Sitko, M. L., Grady, C. A., et al. 2015, ApJ, 798, 94

Wang, J., Graham, J., Dawson, R., et al. 2018, AJ, 156, 192

Whitney, B. A., Robitaille, T. P., Bjorkman, J. E., et al. 2013, ApJS, 207, 30

Wilson, J. C., Henderson, C. P., Herter, T. L., et al. 2004, Proc. SPIE, 5492, 1295

Wisniewski, J. P., Clampin, M., Grady, C. A., et al. 2008, ApJ, 682, 548

Wood, K., Wolff, M. J., Bjorkman, J. E., & Whitney, B. 2002, ApJ, 564, 887

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868