# Multiband Imaging of the HD 36546 Debris Disk: A Refined View from SCExAO/CHARIS*

Frantz Martinache, Kellen Lawson, Thayne Currie, John Wisniewski, Motohide Tamura, Jean-Charles Augereau, Timothy Brandt, Olivier Guyon, N. Jeremy Kasdin, Tyler Groff, et al.

► To cite this version:

Frantz Martinache, Kellen Lawson, Thayne Currie, John Wisniewski, Motohide Tamura, et al. Multi-band Imaging of the HD 36546 Debris Disk: A Refined View from SCExAO/CHARIS*. The Astronomical Journal, 2021, 162 (6), pp.293. 10.3847/1538-3881/ac2823. hal-03518281

HAL Id: hal-03518281

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

Submitted on 2 Aug 2025

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

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

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

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

# Multiband Imaging of the HD 36546 Debris Disk: A Refined View from SCExAO/CHARIS*

Kellen Lawson¹, Thayne Currie²,³,⁴, John P. Wisniewski¹, Motohide Tamura⁵,⁶,⁷, Jean-Charles Augereau⁸, Timothy D. Brandt⁹, Olivier Guyon²,⁵,¹⁰,¹¹, N. Jeremy Kasdin¹²,¹³, Tyler D. Groff¹⁴, Julien Lozi², Vincent Deo², Sebastien Vievard², Jeffrey Chilcote¹⁵, Nemanja Jovanovic¹⁶, Frantz Martinache¹⁷, Nour Skaf²,¹⁸,¹⁹, Thomas Henning²⁰, Gillian Knapp²¹, Jungmi Kwon⁶, Michael W. McElwain¹⁴, Tae-Soo Pyo², Michael L. Sitko²², Taichi Uyama⁷,²³,²⁴, and Kevin Wagner¹⁰,<fnref n="25" />

¹Department of Physics and Astronomy, University of Oklahoma, Norman, OK, USA; kellenlawson@gmail.com

²Subaru Telescope, National Astronomical Observatory of Japan, 650 North A'ohōkū Place, Hilo, HI 96720, USA

³NASA-Ames Research Center, Moffett Blvd., Moffett Field, CA, USA

⁴Eureka Scientific, 2452 Delmer Street Suite 100, Oakland, CA, USA

⁵Astrobiology Center of NINS, 2-21-1, Osawa, Mitaka, Tokyo, 181-8588, Japan

⁶Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan

⁷National Astronomical Observatory of Japan, 2-21-2, Osawa, Mitaka, Tokyo 181-8588, Japan

⁸Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

⁹Department of Physics, University of California, Santa Barbara, Santa Barbara, CA, USA

¹⁰Steward Observatory, The University of Arizona, Tucson, AZ 85721, USA

¹¹The Wyant College of Optical Sciences, University of Arizona, Tucson, AZ 85721, USA

¹²College of Arts and Sciences, University of San Francisco, San Francisco, CA, USA

¹³Department of Mechanical Engineering, Princeton University, Princeton, NJ, USA

¹⁴NASA-Goddard Space Flight Center, Greenbelt, MD, USA

¹⁵Department of Physics, University of Notre Dame, South Bend, IN, USA

¹⁶Department of Astronomy, California Institute of Technology, 1200 East California Boulevard, Pasadena, CA 91125, USA

¹⁷Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, France

¹⁸LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Jules Janssen, F-92195 Meudon, France

¹⁹Department of Physics and Astronomy, University College London, London, UK

²⁰Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany

²¹Department of Astrophysical Science, Princeton University, Peyton Hall, Ivy Lane, Princeton, NJ 08544, USA

²²Space Science Institute, 475 Walnut Street, Suite 205, Boulder, CO 80301, USA

²³Infrared Processing and Analysis Center, California Institute of Technology, 1200 E. California Boulevard, Pasadena, CA 91125, USA

²⁴NASA Exoplanet Science Institute, Pasadena, CA 91125, USA

*Received 2021 May 25; revised 2021 September 2; accepted 2021 September 17; published 2021 December 7*

## Abstract

We present the first multiwavelength (near-infrared; 1.1–2.4 µm) imaging of HD 36546’s debris disk, using the Subaru Coronagraphic Extreme Adaptive Optics (SCExAO) system coupled with the Coronagraphic High Angular Resolution Imaging Spectrograph (CHARIS). As a 3–10 Myr old star, HD 36546 presents a rare opportunity to study a debris disk at very early stages. SCExAO/CHARIS imagery resolves the disk over angular separations of $ρ \sim 0''25–1''0$ (projected separations of $r_{proj} \sim 25–101$ au) and enables the first spectrophotometric analysis of the disk. The disk’s brightness appears symmetric between its eastern and western extents, and it exhibits slightly blue near-infrared colors on average (e.g., $J–K = -0.4 \pm 0.1$)—suggesting copious submicron-sized or highly porous grains. Through detailed modeling adopting a Hong scattering phase function (SPF), instead of the more common Henyey–Greenstein function, and using the differential evolution optimization algorithm, we provide an updated schematic of HD 36546’s disk. The disk has a shallow radial dust density profile ($α_{in} \approx 1.0$ and $α_{out} \approx -1.5$), a fiducial radius of $r_0 \approx 82.7$ au, an inclination of $i \approx 79°1$, and a position angle of $PA \approx 80°1$. Through spine tracing, we find a spine that is consistent with our modeling, but also with a “swept-back wing” geometry. Finally, we provide constraints on companions, including limiting a companion responsible for a marginal Hipparcos–Gaia acceleration to a projected separation of $\lesssim 0''2$ and to a minimum mass of $\lesssim 11 M_{Jup}$.

*Unified Astronomy Thesaurus concepts: Debris disks (363); Circumstellar disks (235); T Tauri stars (1681); Infrared astronomy (786)*

## 1. Introduction

Gas-depleted, dust-dominated *debris disks* are critical targets for understanding the composition, structure, and (typically late-stage) formation of planetary systems (Kenyon & Bromley 2008; Wyatt 2008; Hughes et al. 2018). Scattered-light imaging of debris disks can reveal and characterize signatures of planets responsible

for shaping the disks’ morphologies and help constrain the properties of the planets themselves (e.g., Kalas et al. 2005; Lagrange et al. 2010). “Extreme adaptive-optics” (exAO) facilities, such as SPHERE (Beuzit et al. 2019), GPI (Macintosh et al. 2015), and SCExAO (Jovanovic et al. 2015; Lozi et al. 2018; Currie et al. 2020a), now deliver high-fidelity images of numerous debris disks at subarcsecond separations (e.g., Milli et al. 2017; Duchêne et al. 2020; Esposito et al. 2020; Lawson et al. 2020).

Modeling of multiwavelength imaging can enable analysis of the properties and composition of the constituent dust (e.g.,

<footnote>* Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.</footnote>

<footnote n="25">²⁵ NASA Hubble Fellowship Program—Sagan Fellow.</footnote>
Boccaletti et al. 2003; Fitzgerald et al. 2007; Goebel et al. 2018; Lawson et al. 2020), providing a context for the chemistry and evolution of our own Kuiper Belt (Currie et al. 2015b; Milli et al. 2017; Chen et al. 2020). However, the assumed scattering phase function (SPF) critically affects resulting model disk images—and thus our ability to infer disk properties. The traditionally adopted Henyey–Greenstein (H-G) SPF (Henyey & Greenstein 1941) can reproduce scattered-light images at large scattering angles. However, it is not physically motivated, and it is often inconsistent with data at smaller scattering angles (e.g., Hughes et al. 2018). Combined with evidence suggesting that the best-fit H-G SPF from disk modeling is correlated with the scattering angles probed by the data (rather than any more meaningful physical properties of the system), Hughes et al. (2018) suggest that some other SPF—capable of explaining debris disk imagery irrespective of the scattering angles probed—is needed. Lawson et al. (2020) showed that the Hong (1985) SPF, which is derived from observations of solar system zodiacal dust, reasonably reproduced near-infrared (NIR) imagery of the debris disk of HD 15115. This finding is consistent with the suggestion in Hughes et al. (2018) that a near-universal SPF may exist to explain observed scattered light of a particular wavelength originating from nearly any circumstellar dust. Moreover, Lawson et al. (2020) found that the adopted SPF can fundamentally change the interpretation of scattered-light images, e.g., the number of separate dust populations required to reproduce the data. If scattered-light imagery of additional systems can be modeled well using this SPF, it could provide further evidence of such a universal SPF.

In this study, we present the first multiwavelength scattered-light imaging of the debris disk around HD 36546 (Currie et al. 2017), using the Subaru Coronagraphic Extreme Adaptive Optics (SCExAO) system and the Coronagraphic High Angular Resolution Imaging Spectrograph (CHARIS) integral field spectrograph in broadband (spanning NIR *J*, *H*, and *K* bands, 1.13–2.39 µm) mode (Groff et al. 2016). HD 36546’s disk is highly inclined, bright in scattered light, and has an extremely high fractional luminosity ($L_{IR}/L_* \sim 4 \times 10^{-3}$) comparable to long-studied debris disks around $\beta$ Pictoris and HR 4796A as well as new benchmarks such as HD 115600 (Smith & Terrile 1984; Schneider et al. 2014; Currie et al. 2015b; Millar-Blanchaer et al. 2015; Gibbs et al. 2019).

The HD 36546 system lies at a distance of $101.35 \pm 0.73$ pc (Gaia Collaboration et al. 2018), slightly foreground to the Taurus-Auriga star-forming region and is a member of the 3–10 Myr old 118 Tau association (99.1% probability; Gagné et al. 2018).<fnref n="26" /> Due to the system’s extreme youth, the HD 36546’s debris disk represents a rare opportunity to study debris disk structure and composition at very early stages where protoplanetary disk accretion ceases, colder gas dissipates, and collisions between icy planetesimals tracing Kuiper Belt–analog formation are predicted to yield detectable debris at peak luminosity (Currie et al. 2008, 2009; Cloutier et al. 2014; Ribas et al. 2015; Kenyon & Bromley 2008).

Previous work on HD 36546 found evidence for a disk with an extremely large scale height and strong forward-scattering: an H-G g value of ~0.7–0.85, among the highest of any known debris disk (Currie et al. 2017; Hughes et al. 2018). Revisiting HD 36546’s disk with our new, comparable-quality data obtained over a wider wavelength range will provide improved constraints on its structure and the first look at its composition. Much like for the case with HD 15115 (Lawson et al. 2020),

modeling these data with a different SPF may lead to revisions in our understanding of the HD 36546 disk’s properties.

Our new data resolve the disk down to separations as small as ~0.′′25 (~25 au) and enable the first spectrophotometric analysis of the disk. From these data, we provide surface brightness, surface contrast, and asymmetry profiles in CHARIS broadband and *J*, *H*, and *K* bands—as well as disk color profiles in *J*–*H* and *H*–*K*. Through detailed forward modeling, we provide a refined schematic for HD 36546’s disk. Finally, we provide constraints for the presence of companions within the CHARIS field of view (FOV; angular separations of 0.′′11 $\lesssim \rho \lesssim 1.′′0$).

## 2. Data

### 2.1. Observations

We observed HD 36546 (A0–A2V, $V=6.95$, $H=6.92$; Currie et al. 2017; Lisse et al. 2017) on 2019 January 12 using the Subaru Telescope’s SCExAO running at 2 kHz paired with the CHARIS integral field spectrograph (Groff et al. 2016) operating in low-resolution ($R \sim 20$), broadband (1.16–2.37 µm) mode, and utilizing SCExAO’s Lyot coronagraph with 217 mas diameter occulting spot. Conditions were photometric with good 0.′′4–0.′′5 seeing and 16 km hr$^{-1}$ winds.<fnref n="27" /> The data were collected in *angular differential imaging* mode (ADI; Marois et al. 2006) and achieved total parallactic angle rotation of 99° with total integration time of $t_{int} = 50$ minutes. The data set is made up of 66 individual exposures, with 65 having exposure times of 45.73 s and one having an exposure time of 30.98 s. Sky frames were also obtained following the science data to enable sky subtraction.

A second target, HR 2466 (A2V; $V=5.21$, $H=5.07$), was observed on the same night with the same instrument configuration, to enable reference star differential imaging (RDI). These data include 143 individual exposures of 20.65 s each, for a total integration time of 49 minutes. For all reduction and analysis, we adopt the revised CHARIS pixel scale of 0.′′0162 and the North-up PA offset of 2.′.2 as reported in Currie et al. (2018).

### 2.2. CHARIS Data Reduction

CHARIS data cubes were extracted from raw CHARIS exposures using the CHARIS Data Reduction Pipeline (Brandt et al. 2017), using modifications to read noise suppression listed in Currie et al. (2020a). Extracted data take the form of image cubes with dimensions ($N_\lambda$, $N_x$, $N_y$) = (22, 201, 201) (i.e., $201 \times 201$ pixel images for each of 22 wavelength channels). Subsequent basic image processing—e.g., sky subtraction, image registration, spectrophotometric calibration—was carried out as in Currie et al. (2011, 2018). Inspection of registered data cubes showed that the quality of the AO correction varied between exposures taken before (good), during (worst<fnref n="28" />), and after (best) transit. Particularly low-quality frames were

<footnote n="26">²⁶ This association is labeled as “Mamajek 17” in Currie et al. (2017).</footnote>

<footnote n="27">²⁷ However, AO performance was compromised compared to typical performance due to a significant (4–5x) loss in the IR secondary mirror’s reflectivity at blue optical wavelengths where AO188 (which provides a needed low-order correction for SCExAO) does wavefront sensing. This loss resulted from oxidization of the mirror’s silver coating from sulfur dioxide stemming from the Kilauea eruption in 2018 May. For wavefront sensing, HD 36546 was effectively a magnitude 10.5 star. The mirrors were realuminized with a fresh coating in late 2019.</footnote>

<footnote n="28">²⁸ This is a result of HD 36546 transiting very near zenith from Subaru, making alt/az tracking through transit challenging.</footnote>
<table><caption>**Table 1**
PSF Subtraction Algorithm Settings</caption><thead><tr><th>Method</th><th>Param. Tuning</th><th>Radial Prof. Sub.</th><th>r<sub>min</sub></th><th>r<sub>max</sub></th><th>Δr<sub>sub</sub></th><th>Algorithm Parameters</th></tr></thead><tbody><tr><td>RDI-KLIP</td><td>Disk</td><td>False</td><td>15</td><td>70</td><td>55</td><td>N<sub>PCA</sub> = 10, N<sub>zones</sub> = 1,</td></tr><tr><td>ADI-KLIP</td><td>Disk</td><td>True</td><td>15</td><td>65</td><td>50</td><td>N<sub>PCA</sub> = 15, N<sub>zones</sub> = 1, δ<sub>FWHM</sub> = 2.5</td></tr><tr><td>ADI-ALOCI</td><td>Disk</td><td>True</td><td>15</td><td>70</td><td>5</td><td>N<sub>A</sub> = 800, g = 0.1, δ<sub>FWHM</sub> = 2.5, SVD = 0.02</td></tr><tr><td>SDI+ADI-ALOCI</td><td>Planet</td><td>True</td><td>5</td><td>70</td><td>2.5</td><td>N<sub>A</sub> = 50, g = 0.25, δ<sub>FWHM</sub> = 0.55, SVD = 1.05 × 10<sup>-6</sup></td></tr></tbody></table>

**Note.** Algorithm settings for PSF subtraction with each of the three techniques utilized. “Radial Profile Sub.” indicates whether each reduction was carried out on radial profile subtracted data. Here, $r_{min}$ and $r_{max}$ refer to the inner and outer radii of the subtraction region for each reduction (in pixels); $\Delta r_{sub}$ indicates the radial width of individual subtraction annuli; “g” refers to the aspect ratio of the optimization regions (which are annular subsections, defined in polar coordinates by radial and azimuthal extents), with $g < 1$ producing azimuthally elongated regions and $g > 1$ producing radially elongated regions; “$N_A$” refers to the area of optimization regions in units of PSF cores; “$\delta_{FWHM}$” indicates the minimum rotation gap in units of PSF FWHM (for both A-LOCI and KLIP); “$\Delta r_{sub}$” gives the radial size of subtraction regions in units of pixels (for both A-LOCI and KLIP); “$N_{PCA}$” indicates the number of principal components utilized in construction of the model PSF; “$N_{zones}$” is the number of subsections into which each KLIP optimization annulus was divided (with a value of 1 corresponding to full annuli); and “SVD” refers to the threshold for normalized singular values below which A-LOCI coefficients are truncated. For the ADI-ALOCI reduction, the algorithm is not restricted in the number of images that may be used in constructing the PSF model. For the SDI+ADI-ALOCI parameter tuning, we also employed a moving pixel mask over the subtraction zone (Marois et al. 2010; Currie et al. 2012a) and restricted the reference PSF construction for each target image to the 50 best-correlated frames (Currie et al. 2012b).

identified by first measuring the peak-to-halo flux ratio of the satellite spots in each exposure.<fnref n="29" /> These values were then averaged over the four satellite spots and 22 wavelength channels, and the median and median absolute deviation of the values across the full data sequence were computed. Any exposures falling more than three median absolute deviations below the median peak-to-halo ratio were removed. We removed three frames near transit that met this criteria.

PSF subtraction to recover HD 36546’s disk was performed using RDI with the *Karhunen-Loève* Image Projection (KLIP; Soummer et al. 2012) algorithm (RDI-KLIP), as well as ADI with both KLIP (ADI-KLIP) and the *Adaptive, Locally Optimized Combination of Images* (A-LOCI; Currie et al. 2012a, 2015a) algorithm (ADI-ALOCI), which is based on the LOCI algorithm (Lafrenière et al. 2007). To search for planets, we exploit both ADI and *spectral differential imaging* (SDI; Sparks & Ford 2002; see, e.g., Beuzit et al. 2019), using a more aggressive implementation of A-LOCI (as in, e.g., Currie et al. 2018, 2020a, 2020b; Lawson et al. 2020). Unlike prior implementations, we performed an ADI-based reduction on the post-SDI residuals instead of SDI on the post-ADI residuals. We found that, for data sets with somewhat variable AO corrections leading to temporally decorrelated images, the former approach is slightly superior. The utilized algorithm parameters are provided in Table 1.

## 2.3. Results

In broadband (wavelength-collapsed) imagery, all three PSF-subtraction techniques result in clear detection of the bright side of HD 36546’s disk to separations as small as ~0.′′25 (Figures 1 and 2). Each result also shows a partial, weak detection of the fainter side of the disk in the west, which is strongest in the ADI-ALOCI result, with signal-to-noise ratio per resolution element (SNRE) ~2.5. As in the discovery imagery of the disk from SCExAO/HiCIAO (Currie et al. 2017), SCExAO/CHARIS imagery of HD 36546 reveals a highly inclined, significantly asymmetrically scattering disk. ADI-based reductions reveal no obvious evidence of ansae, indicative of gaps or a ring-like structure, within the field of view. The RDI-KLIP reduction,

however, appears to show a change in the disk profile at $\rho \sim 0.′′75$, consistent with the fiducial radius of the best-fitting disk model identified in Currie et al. (2017).

The SNRE of the disk is calculated following the standard procedure (e.g., Goebel et al. 2018); the signal per resolution element is determined by replacing each pixel with the sum of pixels within an aperture of one PSF full width at half maximum (FWHM), while the noise per resolution element is the radial noise profile of the aforementioned signal with a finite-element correction applied (Mawet et al. 2014). The quotient of these is taken to be the SNRE. The detection in the broadband images (Figure 3) is strongest for the ADI-ALOCI reduction, with an SNRE along the trace of the disk of ~4–9 for regions exterior to 0.′′25.<fnref n="30" /> The ADI-KLIP reduction yields slightly weaker SNRE (typically ~1 less than the same locations in the ADI-ALOCI result). SNRE maps for the RDI-KLIP reduction suggest a somewhat weaker detection having SNRE of ~3–4 at small separations, dropping to ~2 beyond 0.′′4. We note, however, that the noise calculation for the RDI-KLIP reduction is more heavily affected by the disk’s signal because RDI subtraction tends to be more conservative of disk flux. Therefore, the resulting SNRE depends on the size of the noise mask utilized. Using a slightly larger mask results in SNRE of ~5–6 out to 0.′′6, and ~2.5 beyond (Figure 3, upper right). The detection is clearest in *H* band, where the disk appears much like the broadband result. Detections in *J* band and *K* band are weaker, with *J* band suffering especially from prominent residual speckles at small separations, and *K* band from weaker scattering and higher background flux (Figure 4).

# 3. Modeling the Debris Disk of HD 36546

## 3.1. Disk Forward Modeling

To probe the morphology of HD 36546’s disk, we adopt a strategy of forward modeling synthetic disk images as described in Currie et al. (2019) and Lawson et al. (2020), focused on KLIP instead of A-LOCI due to the former’s much faster computational speed.<fnref n="31" /> We create synthetic scattered-light images using a version of the GRaTeR software (Augereau et al. 1999). This

<footnote n="29">²⁹ The peak flux is defined as the flux within an aperture with a radius of one instrumental full width at half maximum (FWHM), centered on the satellite spot. The halo flux is the flux within an annulus centered on the satellite spot, with inner and outer radii of 2-FWHM-dr/2 and 2-FWHM+dr/2 respectively, where dr is the radial extent that results in an annulus enclosing the same area as the peak aperture.</footnote>

<footnote n="30">³⁰ In the S/N calculation, we use a software mask to reduce the amount of disk signal included in the noise estimation.</footnote>

<footnote n="31">³¹ In this case, KLIP is faster because reductions are performed in full annuli, requiring far fewer matrix operations than A-LOCI, which is performed in much smaller annular segments.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Wavelength-collapsed final images for all three reduction methods noted in Section 2.2. The inner and outer masks in each panel have radii of 0.′′25 and 1.′′0, respectively.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** As Figure 1, but with a symmetric logarithmic color scale. This scaling is linear in the range [−0.1, 0.1], and logarithmic otherwise.</figcaption></figure>

model image is then rotated to an array of parallactic angles matching those of the data, and convolved with the instrumental PSF for each of the 22 CHARIS wavelength channels. Here, the PSF is estimated empirically by median combining normalized cutouts of the satellite spots for each wavelength over the full data sequence. Using forward modeling (Pueyo 2016), we simulate annealing of the disk model for each CHARIS data cube in our sequence. In this procedure, oversubtraction (where the presence of disk signal results in an overly bright PSF model), direct self-subtraction (inclusion of the disk signal in basis vectors), and indirect self-subtraction (perturbation of basis vectors by the disk signal) terms are computed from the Karhunen–Loève modes used for the science data reduction and their perturbations by the disk.

Once this process is completed, we compare the disk model to the wavelength-collapsed science image to assess the goodness of fit. For this purpose, we utilize $\chi^2$, which is calculated here within a region of interest as described in Lawson et al. (2020). The region of interest considered is a rectangular box centered on the star with unbinned dimensions of 120 pixels × 50 pixels (~2.′′0 × 0.′′8) and rotated 11° clockwise to approximately align the region with the disk’s major axis. Models that are acceptably consistent with the overall best model’s fitness, $\chi^2_{\nu,\min}$, are defined as those having $\chi^2 \le \chi^2_{\nu,\min} + \sqrt{2/\nu}$ (Thalmann et al. 2013). Here, $\nu$ is the degree of freedom, defined as the difference between the number of bins in the region of interest when binned to

resolution and the number of free model parameters. The disk models are restricted to those with simple one-ring geometries with linear flaring ($\beta=1$), and a Gaussian vertical density distribution ($\gamma=2$).

It is common to describe the angular distribution of scattered light in debris disks using the physically-unmotivated Henyey–Greenstein (H-G) scattering phase function (Henyey & Greenstein 1941). The H-G SPF, parameterized by a single variable —the asymmetry parameter $g^{32}$, can reproduce scattered-light images at large scattering angles but often fails at smaller angles (e.g., Hughes et al. 2018). Moreover, best-fit single H-G scattering asymmetry parameters are well-correlated with the scattering angles probed by observations, while lacking correlation with the properties of the systems (Hughes et al. 2018). This appears to suggest that some other SPF, capable of reproducing observations irrespective of the scattering angles probed, is necessary. Linear combinations of two or more H-G phase functions of differing asymmetry parameter could provide a suitable solution specific to the details of a particular disk (e.g., its composition). However, optimization of a disk model with an SPF formed from the linear combination of N H-G phase functions will introduce either 2N additional free parameters, with an asymmetry parameter ($g_i$) and weight ($w_i$) for each, or 2N–1 free parameters for a constrained optimization problem (with the constraint that the weights sum to

<footnote n="32">32 The asymmetry parameter is defined as the average of the cosine of the scattering angle, weighted by the (normalized) phase function, over all directions.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Signal-to-noise ratio per resolution element for the reductions shown in Figure 1. For all panels, a mask is applied over the disk when the noise is calculated. The RDI reduction is shown using two different mask sizes (overlaid dashed white boxes), as the chosen size was found to have a strong impact on the resulting SNRE. Both ADI reductions use the smaller mask.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** ADI-KLIP reduction of SCExAO/CHARIS HD 36546 data, with wavelength channels combined to produce images comparable to *J* (channels 1–5, 1.16–1.33 µm), *H* (channels 8–14, 1.47–1.80 µm), and *K* (channels 16–21, 1.93–2.29 µm) bands. Images are displayed with a symmetric logarithmic color scale that is linear in the range [−0.1, 0.1] and logarithmic otherwise.</figcaption></figure>

unity, allowing the weight for one asymmetry parameter to be assumed). Without providing additional constraints (e.g., that $w_i \ge w_{i+1}$), these weights and asymmetry parameters will also have significant correlations or degeneracies with one another. Since the utilized SPF is expected to alter the best-fit value of other disk parameters, the SPF cannot be optimized independently. Thus, exploration of such phase functions would introduce significant, possibly intractable, model complexity for attempts to optimize disk models.

Another possibility, evidenced by measurements of both solar system and extrasolar dust, is that there exists a nearly universal

SPF to describe circumstellar dust in a given wavelength regime (Hughes et al. 2018). In such a scenario, this SPF could be measured independently, such as by observation of solar system dust, and then simply adopted for use by groups modeling scattered light from disks. Lawson et al. (2020) showed that imagery of HD 15115 could be reasonably reproduced by adopting the fixed Hong (1985) SPF—modeled as the linear combination of three separate H-G phase functions based on observations of solar system zodiacal dust. Here, we again adopt the Hong (1985) scattering phase function, with asymmetry parameters and weights fixed to the values identified in Hong
(1985): $g_1 = 0.7$, $g_2 = -0.2$, and $g_3 = -0.81$, with weights $w_1 = 0.665$, $w_2 = 0.330$, and $w_3 = 0.005$.

Altogether, our prescription results in a model having six free parameters:

1. $R_0$, the radius of peak grain density in au;

2. $H_0/R_0$, the ratio of disk scale height at $R_0$ to $R_0$;

3. $\alpha_{in}$, the power-law index describing the change in radial density interior to $R_0$;

4. $\alpha_{out}$, the power-law index describing the change in radial density exterior to $R_0$;

5. *i*, the inclination of the disk in degrees;

6. PA, the position angle of the disk in degrees.

To explore these parameters, we utilize the differential evolution (DE) optimization algorithm as described in Lawson et al. (2020). In DE, a population of trial solutions is iteratively improved by adding the weighted difference of two random solutions to a third and keeping only those “mutated” population members that result in an improvement in the fitness metric. Algorithmic parameters are set as follows: $N_{pop} = 15$, resulting in a population size of 90 (6 free parameters ×15), a dithered mutation constant $m = [0.5, 1.0]$, and a crossover probability $P = 0.7$. We test for convergence in the current population at the end of each generation. We consider the population to be converged if either of the following criteria is met:

1. For each varying parameter, the standard deviation among the set of normalized parameter values for the current population falls below a threshold value, $\delta_{\sigma X}$.

2. The standard deviation of the fitness of the current population, divided by the median fitness of the current population, falls below a threshold value, $\delta_{\sigma FX}$, and more than 10 generations have completed.

Effectively, the first criterion is met by a population whose parameters have converged to a single minimum position, while the second accommodates a population with diverse parameter values but comparable goodness of fit, as might occur with weak or degenerate parameters. The additional requirement that 10 generations have finished for the second criteria prevents the unlikely case in which convergence might otherwise be achieved by a population in the first few generations with approximately uniform but weak fitness. Both thresholds, $\delta_{\sigma X}$ and $\delta_{\sigma FX}$, are set to 0.005.<fnref n="33" />

We carry out several distinct procedures for optimizing the disk model. Our primary results focus on a procedure in which the model is optimized for both the RDI-KLIP and ADI-KLIP reductions simultaneously (i.e., seeking the model that minimizes the combined $\chi^2$ for both reductions). However, since the adopted noise profile for the RDI reduction suggests substantially larger noise levels than that of the ADI-KLIP reduction, its contribution to a combined $\chi^2$ metric will be minor; effectively, if the global best-fit solution differs considerably between the RDI and ADI reductions, this would not likely be apparent for such a procedure. As such, we also carry out procedures optimizing the model for each individual reduction. A secondary motivation for this is to validate the consistency of the solutions provided by DE. Toward this latter

point, we also run two separate optimization procedures with identical settings for the RDI-KLIP data.<fnref n="34" />

## 3.2. Model Results

The key results of the four modeling procedures (two for RDI-KLIP, one for ADI-KLIP, and one combined) are provided in Table 2. A corner plot for the combined procedure is provided in Figure 5. For analysis hereafter (e.g., for attenuation correction in Section 4), we adopt the best-fitting model resulting for the combined procedure, having parameters: $R_0 = 82.65$ au, $H_0/R_0 = 0.005$, $\alpha_{in} = 1.00$, $\alpha_{out} = 1.51$, $i = 79.06$, and PA = 80.09 (visualized in Figures 6 and 7). The results of the other modeling procedures are visualized and compared in Appendix A.

## 3.3. Modeling Discussion

The results of our forward-modeling procedure echo many of the themes of the results in Currie et al. (2017), favoring a highly inclined disk model with shallow radial density profiles (small $\alpha_{in}$, $\alpha_{out}$), suggestive of an extended disk as opposed to the ring-like geometry often found for debris disks.

We emphasize here that, while the nominal parameter values are expected to very closely correspond to the global minimum of the parameter space, the acceptable parameter ranges from DE provide only a rough approximation of parameter uncertainties. Unlike ensemble samplers (e.g., Markov Chain Monte Carlo), which draw samples from a parameter space in proportion to the statistical likelihood of those samples, DE is designed to rapidly approach the global minimum with little regard for suboptimal combinations of parameters; in cases where the algorithm identifies a strong value for a particular parameter early in the run, the values of that parameter that yield “acceptable” results may be poorly explored. If a stronger approximation of parameter uncertainties is desired, repeated DE runs (or follow-up runs near the best DE solution with “local” optimization algorithms, such as damped least-squares) can provide thorough exploration of acceptable values while still requiring orders-of-magnitude fewer samples than would be required by an ensemble sampler.

While highly inclined disks are expected to be relatively insensitive to changes in the $\alpha$ parameters (see, e.g., Lawson et al. 2020), these parameters are perhaps the best constrained of any in our modeling (see Figure 5). This might be explained by the particularly small scale height favored by our modeling: $H_0/R_0 = 0.005$, an order of magnitude smaller than the median value of 0.06 provided by Hughes et al. (2018) and the lower boundary for values permitted by our procedure. However, we note that scale heights up to $H_0/R_0 = 0.052$ are also acceptable solutions. In fact, the diagonal subplot for $H_0/R_0$ in Figure 5 shows that the minimum $\chi^2$ changes very little as a function of $H_0/R_0$ for $H_0/R_0 \lesssim 0.05$. Given that the corner plot shows no obvious strong correlations between $H_0/R_0$ and other parameters, this appears to suggest that our data simply cannot discriminate between $H_0/R_0$ values in this range. This is not particularly surprising, given that $H_0/R_0 = 0.05$ for a fiducial radius of $R_0 \sim 83$ au corresponds to $H_0 \sim 4.2$ au, which projects to $\sim 0.004$ or 2.8 pixels, compared to our resolution of 2.67 pixels. Considering this apparent degeneracy for $H_0/R_0$, our results for this parameter should be interpreted as $H_0/R_0 \lesssim 0.05$.

<footnote n="33">³³ Typically, larger values ~0.01 are chosen. However, we opt for a stricter convergence requirement here, purely for the sake of producing full corner plot visualizations—the solutions show no statistically significant changes between a threshold of 0.01 and 0.005.</footnote>

<footnote n="34">³⁴ Since DE involves a randomly initialized set of trial solutions that are then randomly mutated over subsequent generations, it is conceivable that two DE procedures could converge to distinct solutions.</footnote>
<table><caption>**Table 2**
Disk Model Optimization Results</caption><thead><tr><th></th><th>R<sub>0</sub> (au)</th><th>H<sub>0</sub>/R<sub>0</sub></th><th>α<sub>in</sub></th><th>α<sub>out</sub></th><th>i (deg)</th><th>PA (deg)</th><th>χ<sup>2</sup><sub>ν</sub></th></tr></thead><tbody><tr><td>RDI-KLIP (1)<br>N = 1980</td><td>75.4<br>[60.0, 103.2]</td><td>0.007<br>[0.005, 0.143]</td><td>1.0<br>[1.0, 5.0]</td><td>-1.4<br>[-2.1, -1.0]</td><td>78.9<br>[75.0, 81.4]</td><td>79.0<br>[75.1, 82.4]</td><td>0.45<br>&lt;0.50</td></tr><tr><td>RDI-KLIP (2)<br>N = 2340</td><td>76.4<br>[60.0, 102.7]</td><td>0.012<br>[0.005, 0.138]</td><td>1.0<br>[1.0, 5.4]</td><td>-1.4<br>[-2.0, -1.0]</td><td>78.9<br>[75.0, 81.9]</td><td>78.9<br>[75.5, 82.4]</td><td>0.45<br>&lt;0.50</td></tr><tr><td>ADI-KLIP<br>N = 1800</td><td>87.4<br>[70.3, 114.6]</td><td>0.005<br>[0.005, 0.050]</td><td>1.0<br>[1.0, 1.9]</td><td>-1.6<br>[-2.0, -1.3]</td><td>79.3<br>[76.9, 81.3]</td><td>80.6<br>[78.7, 82.6]</td><td>0.76<br>&lt;0.81</td></tr><tr><td>ADI & RDI<br>N = 2160</td><td>82.7<br>[64.1, 109.0]</td><td>0.005<br>[0.005, 0.052]</td><td>1.0<br>[1.0, 1.9]</td><td>-1.5<br>[-2.0, -1.2]</td><td>79.1<br>[77.3, 80.7]</td><td>80.1<br>[78.6, 82.0]</td><td>0.61<br>&lt;0.65</td></tr><tr><td>Currie et al. (2017)<sup>a</sup></td><td>75.6<br>[66.7, 84.5]</td><td>0.118<br>[0.070, 0.160]</td><td>3<br>[3, 10]</td><td>-3<br>[-3, -3]</td><td>75<br>[70, 75]</td><td>75</td><td>1.03<br>&lt;1.07</td></tr><tr><td>Bounds</td><td>[60, 120]</td><td>[0.005, 0.2]</td><td>[1, 10]</td><td>[-10, -1]</td><td>[75, 85]</td><td>[75, 85]</td><td>...</td></tr></tbody></table>

**Notes.** Results for each of the four distinct DE optimization runs outlined in Section 3.1, as well as the forward modeling of HD 36546 in Currie et al. (2017). Each parameter column provides the value for the best model for the run indicated in the leftmost column (with *N* indicating the total number of models evaluated for each procedure), with the following line providing the range of values resulting in an acceptably fitting model. The final column gives the corresponding run's minimum $\chi^2$ and the threshold for acceptable models. The permitted boundary for each parameter in modeling for CHARIS data (the modeling procedure of Currie et al. (2017) used different bounds) is given on the final line. We adopt the best model for the combined run ("ADI & RDI") for analysis hereafter. *a*.

a In Currie et al. (2017), models are discussed in terms of $H_0$ ($\xi_0$ or $ksi_0$ in the text) rather than $H_0/R_0$; for simplicity, we present these results in terms of $H_0/R_0$ instead. Additionally, we have updated relevant values for the difference in the assumed distance to the system—114 pc in Currie et al. (2017), versus 101.35 pc here.

Additionally, we note that the value of $\alpha_{in}$ converges to the lower bound for permitted values ($\alpha_{in} = 1$). This result, indicating an exceptionally slow decline in density interior to $r_0$, could suggest that the disk is better described by a geometry with a nearly flat density slope over some intermediate region, followed by typical power-law slopes interior and exterior to that region. An alternative explanation for the strong $\alpha_{out}$ constraint that manifests here is addressed briefly in Section 4.1.

The best-fit parameter values and acceptable ranges resulting from forward modeling in Currie et al. (2017) are included for comparison in Table 2. Notably, this modeling was conducted on a coarse grid; the acceptable ranges provided are much less likely to correspond to complete ranges of acceptably fitting parameters. Additionally, we note that, while the value for PA identified in Currie et al. (2017) is at the lower bound of allowed values for our presented DE runs, in earlier modeling runs not presented here, the range of allowed PA values extended as low as 60° but showed very poor fits for PA $\lesssim$ 75°. As such, the final procedure restricted PA to the provided range to speed convergence. As stated previously, our results are thematically consistent with the picture of the disk suggested by the analysis in Currie et al. (2017). Many of the specific parameter values identified, however, appear inconsistent; e.g., despite evaluating models with $\alpha_{out} = -1$, the best-fit in Currie et al. (2017) is achieved for $\alpha_{out} = -3$. Forward modeling the best-fit solution from Currie et al. (2017) for our reductions results in a combined $\chi^2$ value of 0.85 (well beyond the limit for acceptable models of $\chi^2 = 0.65$). This apparent difference in identified optimal parameters might be explained by some combination of a few key differences in the modeling procedures:

1. the data itself, e.g., with the CHARIS reduction making a partial detection of the fainter side of the disk where HiCIAO data did not;

2. the scattering phase functions, with our choice to adopt the Hong phase function, as opposed to fitting an asymmetry parameter as in Currie et al. (2017), which

was found, in prior work modeling the disk of HD 15115 (Lawson et al. 2020), to have a significant impact on the optimal parameters;

3. the choice in Currie et al. (2017) to not vary PA, instead adopting the PA resulting from spine tracing;

4. the parameter bounds, e.g., with values of $\alpha_{in}$ in Currie et al. (2017) not including the best-fit value that we identify;

5. the parameter optimization, with Currie et al. (2017) identifying best-fit parameters through a much less thorough grid search.

The best-fit model for the combined ("ADI & RDI") DE run results in a reduced $\chi^2$ score for each reduction that is very similar to the best scores achieved in the DE runs for individual reductions. The DE run for only the ADI-KLIP reduction results in a best-fit model having $\chi^2 = 0.758$ (with "acceptable" models having $\chi^2 \le 0.808$), while the best-fit model of the combined run results in $\chi^2 = 0.760$ when considering only the ADI-KLIP data well within the limit to be considered an acceptable solution. Similarly, the DE runs for only the RDI-KLIP reduction result in a best-fit model having $\chi^2 = 0.435$ (with $\chi^2 \le 0.504$ being acceptable), while the best-fit model of the combined run results in $\chi^2 = 0.459$ for the RDI-KLIP reduction. In other words, the adopted disk model from our combined modeling run is also a reasonable explanation when considering the imagery of each data reduction individually.

## 4. Disk Spine Trace and Surface Brightness

To measure the disk's spine and surface brightness (SB), we proceed largely as detailed in Lawson et al. (2020) and summarized hereafter. In place of the RDI-KLIP reduction utilized for modeling, we perform spine fitting and surface brightness measurements on an RDI-KLIP reduction in which radial profile subtraction is utilized (with all other parameters remaining unchanged). This results in a generally weaker detection of the
<figure><img src="image_7.png" /><figcaption>**Figure 5.** A variation of the corner plot for optimization of the disk model to both the RDI-KLIP and ADI-KLIP reductions simultaneously (corresponding to the row labeled "ADI & RDI" in Table 2). Each off-diagonal panel visualizes solutions as a function of two of the parameters. Each point in the subplots corresponds to a single sample, colored according to its $\chi^2$ score. The color mapping scales from $\chi^2_{\nu,\min}$ to $\chi^2_{\nu,\min} + 3 \cdot \sqrt{2/\nu}$, with darker colors indicating more optimal values of $\chi^2_{\nu}$. Diagonal elements provide a one-dimensional view of each of the parameters, while also serving as a color bar for the off-diagonals. The threshold fitness for acceptable solutions ($\chi^2_{\nu} \le \chi^2_{\nu,\min} + \sqrt{2/\nu}$) is indicated by a horizontal black dashed line in each diagonal plot. The black-outlined star (triangle) marker in each off-diagonal (diagonal) panel indicates the position of the overall best-fitting model. The text at the top of each diagonal subplot provides the best-fitting value for the corresponding parameter, with the range beneath indicating the smallest and largest value producing an acceptable $\chi^2$; these limits are also depicted using error bars for the "best-fit" marker in each diagonal panel.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** From left to right, for the RDI-KLIP reduction: the PSF-convolved input model, attenuated model, data, and residual for the overall best model resulting from the combined DE run. The gray dashed line overlaid in each panel is the ellipse corresponding to the radius, PA, and inclination of the model.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** As Figure 6, but for the ADI-KLIP reduction.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Stellocentric separation of disk spine fits along the major and minor axes (for a major axis position angle of 80.09°). The compass provided indicates both the orientation of the figure and the difference in scale between the axes (which results in $N \propto E$, as depicted). The spine positions are identified by fitting Lorentzian profiles to both wavelength-collapsed KLIP images (ADI and RDI) simultaneously, and may indicate a swept-back wing profile (see Section 4.1). The dashed blue ellipse corresponds to the adopted best-fitting model of Section 3.1 ($r_0 = 82.65$ au, $i = 79.06^\circ$, PA = 80.09°), and is shown to demonstrate that the ellipsoidal profile very poorly describes the spine profile measured.</figcaption></figure>

disk at very small separations and in wavelength-collapsed imagery, but provides a more consistent detection across the field of view and in individual wavelength channels. For the sake of attenuation correction, we forward model the adopted best-fitting model of Section 3.1 for this additional reduction. The wavelength-collapsed image for this reduction and the corresponding forward modeling results are provided in Appendix B.

## 4.1. Disk Spine

The positions of peak radial brightness along the disk are referred to as the disk's spine, and are typically used as locations for surface brightness measurements and to approximate the geometric parameters of the disk (e.g., fiducial radius, inclination, and position angle). We fit the spine position using wavelength-collapsed products and applying an initial image rotation of -11° to approximately horizontally align the disk's major axis. We then measure the spine at each pixel column position by fitting Lorentzian profiles to the image values weighted by the inverse of the radial noise profile. The fit Lorentzian peak and associated fit uncertainty are taken as the spine position and uncertainty. As a departure from the method of Lawson et al. (2020), the fitting is carried out on both the ADI-KLIP and RDI-KLIP images simultaneously, to identify the spine position most consistent with both images (rather than taking the weighted average of the results).

Notably, the fit spine (Figure 8) appears to grow further from the major axis at larger separations, suggesting a possible “swept-back wing” geometry (see, e.g., Hughes et al. 2018), in which a disk’s outer dust ring appears to bend (rather than remaining planar) at increasing separations. As a result of this, the measured spine is extremely poorly described by an ellipse and thus cannot easily provide a meaningful assessment of the disk’s inclination or radius. While it may be possible to estimate the position angle of the disk’s major axis from this spine profile, it is unclear how the value would be affected by the noted phenomenon. As such, we provide no estimate of the disk geometry based on spine fitting. For the locations of surface brightness measurements, we use the fit spine profile in place of the best-fit ellipse profile used for our analysis in Lawson et al. (2020)—correcting the initial rotation of -11° using a rotational transformation of the spine positions.

Interestingly, applying this spine fitting procedure to our best-fitting models also results in a divergent spine trace (both before and after attenuation is induced by forward modeling). After comparison to the complete sample of models, this behavior manifests only if $|\alpha_{out}| < 2$ (Figure 9). This could suggest that the observed spine shape is simply the result of a disk with a dust distribution meeting this criteria and viewed nearly edge-on. Alternatively, models having such values of $\alpha_{out}$ may have been the only accessible means by which to reproduce this feature of the observations—with the true cause
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Disk models forward-modeled for our ADI-KLIP reduction with the parameters of the overall best model identified in Section 3.1 but with differing values of $\alpha_{\text{out}}$. To emphasize the apparent disk spine, a $1/r^2$ correction and arbitrary color scale have been applied to each image. The dotted black ellipse in each panel corresponds to the radius, PA, and inclination of the model, while the blue line shows the spine trace found for the model following the procedure used for the real data (Section 4.1). The spine trace for the real data is overlaid in white. For $\alpha_{\text{out}} \gtrsim -2$, the fit spine deviates considerably from the profile of the ellipse. The model having $\alpha_{\text{out}} = -1.51$ corresponds to the overall best model identified in Section 3.1.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Measurements of peak surface brightness for HD 36546's disk as a function of projected stellocentric separation in arcseconds and au along the major axis in CHARIS data. Here, “broadband” refers to CHARIS’s broadband. Since measurements are made using an aperture of radius $0''06$, comparable SB measurements can only be provided to within $0''06$ of the smallest separations at which the disk is recovered ($0''25$)—hence, SB is provided only to $\rho \sim 0''31$. We note that the measurements here and in other Section 4 figures are not statistically independent from neighboring measurements as a result of the $0''12$ diameter aperture. The apparent dip in western K-band surface brightness at small separations is not likely physical—given that no comparable behavior manifests in other bands, we conclude that this is probably the result of the somewhat compromised K-band data quality that was noted in Section 2.3.</figcaption></figure>

being more similar to the explanations summarized by Hughes et al. (2018; e.g., ISM interaction).

## 4.2. Disk Surface Brightness

Peak surface brightness is measured for both the RDI-KLIP and ADI-KLIP reductions at the aforementioned spine locations using an aperture of diameter $0''12$ (7.4 pixels). As in Lawson et al. (2020), (a) we apply wavelength-dependent attenuation corrections based on the adopted disk model (determined by comparison of the PSF-convolved input model cube and the attenuated model cube after forward modeling), and (b) we approximate the uncertainty for each SB measurement as the standard deviation of an array of like measurements taken at the same stellocentric radius but with azimuthal angles placed every $10^\circ$. These measurements are made on the “broadband” wavelength-collapsed images, as well as images corresponding to the *J*, *H*, and *K* bands. The final surface brightness profile for a given filter is taken to be the weighted average of the results for the two individual reductions. The results of this procedure are visualized in Figure 10.

Disk surface contrast (the surface brightness of the disk relative to the stellar flux) is then computed by combining disk SB with measurements of the stellar flux for each filter (Figure 11).

Any east-west asymmetry is measured by comparing the surface brightness of opposing sides in a given filter (Figure 12).

Finally, the disk color in isolation from the stellar color is measured by converting surface contrasts to magnitudes and taking the difference of measurements between filters (see Figure 13), e.g., $(J_{\text{disk}} - J_{\text{star}}) - (K_{\text{disk}} - K_{\text{star}}) = (J - K)_{\text{disk}} - (J - K)_{\text{star}} = \Delta(J - K)$, where $\Delta(J - K)$ is the *J* − *K* color of the disk without contribution from the star.

## 4.3. Surface Brightness Power Law

For the broadband (wavelength-collapsed) images, we make an additional set of measurements to facilitate fitting a power law to the disk’s radial surface brightness profile. For this purpose, we follow the SB measurement procedure of Section 4.2, except:
<figure><img src="image_13.png" /><figcaption>Figure 11. Measurements of surface contrast for HD 36546's disk in CHARIS data.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12. Relative surface brightness between the eastern and western extents of the disk for CHARIS *J*, *H*, *K*, and broadband. In each subplot, the horizontal colored dashed line indicates the weighted average of the full set of constituent flux measurements, while the colored region corresponds to the 3σ confidence interval for the average. The averages for all four filters are within 1σ of zero—i.e., consistent with an overall symmetric disk flux.</figcaption></figure>

1. measurements are made at positions with an array of radial separations but corresponding to a single scattering angle (25°<fnref n="35" />, assuming PA = 80°.09 and i = 79°.06 based on the results of Section 3.2);

2. images are binned to resolution (FWHM = 2.67 pixels = 0.′′043);

3. the aperture diameter is set to one resolution element.

The first change enables characterization of the radial dependence of the disk's SB for comparison with results in the literature, while the latter two changes result in approximately statistically independent SB measurements for fitting. The locations of these measurements are overlaid on the binned ADI-KLIP image in the top panel of Figure 14. We then fit power laws to these SB measurements as a function of deprojected stlocentric separation<fnref n="36" /> for the eastern and western extents of the disk separately (Figure 14). We find best-fitting power-law slopes of $k = -2.45 \pm 0.14$ and $k = -2.50 \pm 0.21$ for the west and east, respectively.

### 4.4. Surface Brightness Discussion

Disk Color—Our individual measurements of disk color along the apparent spine of the disk (Figure 13) appear generally consistent with a disk scattering profile that is gray in NIR. However, taking the inverse-variance weighted average of the individual measurements for each color suggests a slightly blue color for the disk in both cases: $J - H = -0.19 \pm 0.05$ and $J - K = -0.37 \pm 0.06$.<fnref n="37" /> Notably, the locations of the apparent disk spine at which this blue color was measured fall beyond the projected location of peak density for the disk. To assess the properties of the bulk population of dust in the disk, rather than just those of the dust pushed to larger separations by pressure from the parent star, we make additional measurements of the disk's color in the vicinity of the disk's peak radial density. For this purpose, we measure surface brightness in apertures of the same size (0.′′06) along the brighter side of the

<footnote n="35">³⁵ The particular scattering angle is chosen arbitrarily in an attempt to probe a scattering angle providing radial coverage with both generally high disk signal and well-defined attenuation estimates in the east and west. Values within ~5° provide comparable results.</footnote>

<footnote n="36">³⁶ Derived deprojected stlocentric separations consider only the determined PA and inclination of the system, and assume that our line of sight probes light scattered from the disk midplane. A more detailed assessment would take scale height, flaring, and similar effects into account, but is beyond the scope of this work.</footnote>

<footnote n="37">³⁷ This general result appears to be robust to the specific disk model adopted for the sake of attenuation correction. Calculating the average disk colors using each of the 691 acceptable disk models for attenuation correction instead resulted in color ranges of $J - H = [-0.23, -0.16]$ and $J - K = [-0.43, -0.32]$ with uncertainties comparable to the values found using the best-fit disk for attenuation correction.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Disk color as a function of stellocentric separation along the disk major axis (see Section 4), with $\Delta(J - K) = (J - K)_{disk} - (J - K)_{star}$, etc. Major axis separation <0 corresponds to the eastern side of the disk. Horizontal lines and shaded regions of the corresponding color represent the overall average and $1\sigma$ confidence, respectively. The averages suggest an overall slightly blue color with $\ge 3\sigma$ confidence. $\Delta(H - K)$ manifests similarly to $\Delta(J - H)$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Broadband radial surface brightness (SB) profiles created following the procedure outlined in Section 4.3. The top panel shows the uncorrected ADI-KLIP image binned to resolution, with markers indicating the locations of SB measurements. The lower panel shows the radial surface brightness measurements, with east/west markers offset very slightly from one another in the x-axis direction to avoid overlap. Power-law fits to the data are shown as dashed lines of the corresponding color (with shaded regions indicating the $1-\sigma$ confidence intervals). Separate sides of the disk are generally in good agreement, with confidence intervals overlapping significantly.</figcaption></figure>

disk at the projected location of the adopted disk model's density peak (i.e., along the southern edge of the ellipse overlaid in Figure 8—corresponding to a single deprojected stellocentric separation). From these measurements, we proceed to average color measurements as before, finding $J - H = -0.13 \pm 0.08$, $H - K = -0.30 \pm 0.09$, and $J - K = -0.40 \pm 0.11$. These results again suggest a slightly blue NIR color for the disk.

Comparison with the numerical simulations in Boccaletti et al. (2003), which model NIR disk colors as a function of the dust size distribution's minimum grain size ($a_{min}$) and porosity (P), can enable some constraints regarding the size distribution of the dust in HD 36546's disk. With $P=0$, for $1.1~\mu$m versus $2.2~\mu$m ($\sim J - K$), a blue disk color is achieved only for

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Broadband ($1.13-2.39~\mu$m) contrast curve (solid black line) for the planet detection reduction of CHARIS HD 36546 data, as outlined in Section 2.2. The $5\sigma$ contrast is given as a function of stellocentric angular separation (arcsec, lower x-axis) and projected separation (au, upper x-axis). The right edge of the figure gives contrast mapped to companion masses for 3 Myr (red) and 10 Myr (blue) planets, based on hot start, solar metallicity, hybrid cloud planet evolution models of Spiegel & Burrows (2012). Overplotted lines of the corresponding color indicate the minimum mass and projected separations corresponding to the companion implied by HGCA measurements, mapped to contrasts based on the same models. The displayed ages represent the lower and upper limits of age estimates for HD 36546.</figcaption></figure>

$a_{min} \lesssim 0.25~\mu$m. Similarly, for $1.6~\mu$m versus $2.2~\mu$m ($\sim J - K$), a blue color is achieved for $a_{min} \lesssim 0.3~\mu$m. Finally, for $1.1~\mu$m versus $1.6~\mu$m ($\sim J - H$), a blue color is achieved for $a_{min} \lesssim 0.2~\mu$m.<fnref n="38" /> Based on these simulations, we can estimate a minimum grain size (assuming $P=0$) of $a_{min} \lesssim 0.2~\mu$m to produce the blue color that we observe. The minimum grain size producing a blue color tends to increase for higher

<footnote n="38">³⁸ The models of Boccaletti et al. (2003) are presented in terms of $\sim J - K$ and $H - K$. For our purposes, we combine the provided measurements to present these results in terms of $J - H$ as well.</footnote>
<table><caption>**Table 3**

HGCA Measurements of HD 36546</caption><thead><tr><th>Source</th><th>μ<sub>α</sub></th><th>σ<sub>α</sub></th><th>μ<sub>δ</sub></th><th>σ<sub>δ</sub></th><th>Corr</th><th>t<sub>α</sub></th><th>t<sub>δ</sub></th></tr></thead><tbody><tr><td>Hip</td><td>8.513</td><td>0.974</td><td>-41.876</td><td>0.472</td><td>0.257</td><td>1991.15</td><td>1990.76</td></tr><tr><td>HG</td><td>7.368</td><td>0.028</td><td>-41.226</td><td>0.014</td><td>0.087</td><td></td><td></td></tr><tr><td>Gaia</td><td>7.508</td><td>0.056</td><td>-41.305</td><td>0.040</td><td>-0.137</td><td>2016.00</td><td>2016.50</td></tr></tbody></table>

Note. HGCA calibrated proper motion measurements for HD 36546 (Brandt 2018) as updated with Gaia eDR3 measurements in Brandt et al. (2021, submitted).

porosity, e.g., with $a_{min} \lesssim 3$ µm producing a blue $J-K$ color for $P = 0.95$.

**Disk Asymmetry**—Both individual measurements and the overall averages for east-west flux asymmetry show no statistically significant evidence for a disk that is asymmetric in brightness in any of the four filters analyzed (Figure 12). While the processed imagery of the disk may initially suggest some asymmetry is present—with the western extent of the disk appearing brighter in all three presented reductions in Figure 2—no asymmetry is evident in SB measurements once attenuation corrections are applied. Thus, we conclude that the apparent asymmetry in CHARIS imagery is an artifact of PSF subtraction and is not truly astrophysical in nature.

**Disk Power Law**—The results of power-law fitting (Figure 14) are in good agreement with the results of our asymmetry measurements (Figure 12)—again suggesting that the disk’s flux is largely symmetrical at NIR wavelengths in the $\rho = 0.0''3-1.0''0$ region of CHARIS’s detection.

Comparing our best-fitting power-law slope of $k \sim -2.5$ for HD 36546 ($\rho \sim 25-100$ au) to values in the literature for other debris disks provides further evidence that HD 36546 has a particularly shallow density slope. For the debris disk of HD 32297 ($\rho \sim 30-280$ au), Currie et al. (2012b) find $k \sim -6$, and for the debris disk of HR 4796 A ($\rho \sim 80-120$ au), Thalmann et al. (2011) report $k \sim -10$.

## 5. Limits on Planets

For the planet detection reduction outlined in Section 2.2, we computed $5\sigma$ contrast limits in CHARIS broadband, following the procedure of Currie et al. (2018) for computing planet throughput corrections.

Contrasts are about a factor of 2–3 worse at $\rho < 0.0''75$ than good performance for bright diskless stars (Currie et al. 2020a) and slightly poorer than our previous limits for the debris-disk-bearing HD 15115 (Lawson et al. 2020). However, despite a loss of wavefront-sensing sensitivity (see Section 2), we reach contrasts at mid-spatial frequencies slightly better than those achieved in Currie et al. (2017) with SCExAO/HiCIAO (e.g., $4.4 \times 10^{-6}$ at $0.0''45$ with CHARIS, versus $\sim 8 \times 10^{-6}$ with HiCIAO). This gain is likely due to the use of SDI for additional speckle suppression.

To convert from contrast to planet mass limits, we integrated the hot start, solar metallicity, hybrid cloud, synthetic planet spectra provided by Spiegel & Burrows (2012) over the CHARIS bandpass and divided by the integrated stellar spectrum. We consider two mappings between contrast and mass, corresponding to the lower and upper age limits estimated for the system (3 Myr and 10 Myr).

The right axis of Figure 15 denotes these mass limits. For $r \gtrsim 35$ au ($\rho \sim 0.0''35$), our reduction is sensitive to $\sim 2 M_{Jup}$ planets at $\sim 3$ Myr, and $\sim 5 M_{Jup}$ planets at $\sim 10$ Myr. Despite achieving generally inferior contrasts compared to the best

reduction of CHARIS data for the debris-disk-bearing star HD 15115 (Lawson et al. 2020), we are able to probe slightly less massive planets at a given projected separation due to HD 36546’s relative youth.

While our data do not result in the detection of any candidate planets, HD 36546 may show tentative, indirect evidence for a massive and (in principle) imageable companion. The Hipparcos-Gaia Catalog of Accelerations (HGCA; Brandt 2018) provides calibrated proper motion data for $\sim 115,000$ nearby stars spanning the $\sim 24$ yr baseline between Hipparcos and Gaia astrometry. Accelerations based on this data have been used to measure dynamical masses of directly imaged exoplanets (e.g., Brandt et al 2019), as well as to provide evidence of gravitational interaction with as-yet unseen companions for targeted direct imaging searches (Currie et al. 2020b; Steiger et al. 2021). Brandt (2021) updates HGCA, replacing Gaia DR2 measurements with those from Gaia eDR3, which yield a factor of $\sim 3$ improvement in astrometric precision. Using these updated Gaia measurements, HGCA for HD 36546 (Table 3) indicates a marginally significant acceleration relative to a model of constant proper motion—$\chi^2 = 7.65$, or $2.3\sigma$ for two degrees of freedom—that could be indicative of a companion.

Following equations in Brandt et al. (2019) to estimate companion mass as a function of projected separation using HGCA accelerations, we can assess the completeness of our search for the HD 36546 companion that would induce this acceleration. Since radial velocity (RV) measurements are not available for HD 36546, we adopt an RV acceleration of zero to compute a minimum mass instead. Mapping of the resulting minimum masses to contrasts for 3 Myr and 10 Myr are overlaid on Figure 15. These results suggest that we would likely recover the HGCA companion for a projected separation beyond $\sim 0.0''15$ assuming an age of 3 Myr, or $\sim 0.0''18$ for an age of 10 Myr.

## 6. Conclusions and Future Work

Our study provides the first multiwavelength imaging of the HD 36546 disk and recovers disk signal to separations as small as $\sim 0.0''25$. Through detailed modeling implementing a Hong (1985) scattering phase function and using the differential evolution optimization algorithm, we provide an updated schematic of HD 36546’s disk. This schematic suggests a disk with a particularly shallow radial dust density profile, a fiducial radius of $\sim 83$ au, an inclination of $\sim 79^\circ$, and a position angle of $\sim 80^\circ$. Through spine tracing, we find a disk spine that is consistent with our modeling results, but also consistent with a “swept-back wing” geometry. Surface brightness measurements show no significant flux asymmetry between the eastern and western extents of the disk, and slightly blue NIR colors on average. While we report no evidence of directly imaged companions, we provide contrast limits to constrain the most plausible remaining companions.

Notably, prior work modeling the debris disks of HD 36546 and HD 15115 in NIR total intensity data identified very
distinct optimal HG SPF asymmetry parameters of $g = 0.70–0.85$ (Currie et al. 2017) and $g = 0.3–0.4$ (Engler et al. 2019) respectively. Alongside the results of Lawson et al. (2020), we have shown that imagery of both disks can be reasonably reproduced by models using the same Hong (1985) SPF instead. As is suggested in Hughes et al. (2018), this may support the existence of a nearly universal SPF for circumstellar dust. Additional data that can be used to study smaller scattering angles and additional disk systems will be necessary to clarify this possibility. If such an SPF is verified, it will enable both better exploration of model parameters (by reducing the number of varying parameters) and tighter constraints on parameter values (by mitigating parameter degeneracies).

Deeper follow-up observations with SCExAO/CHARIS would provide (a) further constraints on companions within the disk and (b) higher SNR measurements of the disk color, enabling a more detailed assessment of the dust properties of the disk. The upcoming replacement of Subaru's AO188 system and additional commissioned improvements to wavefront sensing are expected to improve SCExAO/CHARIS contrasts at small separations by roughly a factor of 10 by late 2021 (Currie et al. 2020a). Assuming a system age of 3 Myr, these contrasts would enable detection of ~1 $M_{Jup}$-mass planets at ~25 au. A future nondetection would likely restrict the location of any super-Jovian-mass planet to $\rho \lesssim 0''1$.

Observations measuring the polarized intensity of the disk at comparable scales—or the total intensity of the disk in a wider field of view—could clarify the underlying cause of the noted diverging spine trace. CHARIS's integral field spectropolarimetry mode or space-based imaging of the system could provide these data. Further, polarized intensity imaging would enable even tighter constraints on disk model parameters by better probing regions of the disk where the polarized intensity phase function results in more favorable sensitivity.

We thank our referee, whose comments helped us to improve both the content and clarity of this manuscript.

This research is based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. We are honored and grateful for the opportunity of observing the Universe from Maunakea, which has cultural, historical, and natural significance in Hawaii.

We wish to acknowledge the critical importance of the current and recent Subaru telescope operators, daycrew, computer support, and office staff employees. Their expertise,

ingenuity, and dedication is indispensable to the continued successful operation of Subaru.

The development of SCExAO was supported by the Japan Society for the Promotion of Science (Grant-in-Aid for Research \#23340051, \#26220704, \#23103002, \#19H00703 & \#19H00695), the Astrobiology Center of the National Institutes of Natural Sciences, Japan, the Mt Cuba Foundation and the director's contingency fund at Subaru Telescope. We acknowledge funding support from the NASA XRP program via grants 80NSSC20K0252 and NNX17AF88G. T.C. was supported by a NASA Senior Postdoctoral Fellowship. M.T. is supported by JSPS KAKENHI grant Nos.18H05442, 15H02063, and 22000005. K.W. acknowledges support from NASA through the NASA Hubble Fellowship grant HST-HF2-51472.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.

## Appendix A
### Additional Model Results

The results of our additional differential evolution disk model optimization procedures are included here, with the information of Table 2 graphically summarized in Figure 16. Individual reductions are in close agreement with one another: while the optimal values for some parameters (e.g., $R_0$) differ noticeably compared to the permitted parameter bounds, each reduction's best-fit values fall within the acceptable ranges for the other reductions.

Comparison of the samples for $\alpha_{in}$ between the two RDI-KLIP DE runs (the first two rows of Figure 16) illustrates our prior warning regarding approximation of parameter uncertainties from DE (see Section 3.3). While the runs ultimately yield very similar acceptable ranges, the comparable upper limit achieved for “RDI-KLIP (2)” is the result of a single sample; except for the lone acceptable model with $\alpha_{in} \sim 5.5$, the acceptable range of $\alpha_{in}$ values for this run would manifest much more comparably to that of the ADI-KLIP run.

As a result of the larger noise levels in the RDI-KLIP reduction, the RDI-KLIP modeling procedure (Figure 17) shows much weaker constraints than the ADI-KLIP (Figure 18) modeling procedure. Despite some differences in the resulting parameter values, the best-fit solutions for these procedures manifest very similarly to those of the combined run (e.g., Figures 6 and 7).
<figure><img src="image_18.png" /><figcaption>**Figure 16.** The diagonal (one-dimensional) elements of corner plots for each of the four DE procedures (as indicated to the right of each row of subplots), as summarized in Table 2. Here, $\chi^2$ is presented normalized as $\chi^2_{\nu,\text{norm}} = (\chi^2_{\nu} - \chi^2_{\nu,\text{min}})/\sqrt{2/\nu}$, such that the best model for each run falls at $\chi^2_{\nu,\text{norm}} = 0$ and acceptable models have $\chi^2_{\nu,\text{norm}} \le 1$. Parameter bounds are truncated to better show the distributions in the vicinity of the minima.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 17.** As Figure 5, but for optimization of the disk model to only the RDI-KLIP reduction of HD 36546. This figure includes both runs for the RDI-KLIP data (corresponding to the entries for “RDI-KLIP (1)” and “RDI-KLIP (2)” in Table 2).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 18.** As Figure 5, but for optimization of the disk model to only the ADI-KLIP reduction of HD 36546 (corresponding to the entry for “ADI-KLIP” in Table 2). Compared to that of the RDI-KLIP data, the ADI-KLIP optimization results in stronger constraints for the model parameters by virtue of a higher signal-to-noise ratio.</figcaption></figure>
## Appendix B

### RDI-KLIP Reduction Using Radial Profile Subtraction

The results of forward modeling the adopted best-fitting model for the radial profile subtracted RDI-KLIP reduction (see Section 4) are provided in Figure 19.

<figure><img src="image_21.png" /><figcaption>**Figure 19.** As Figure 6, but for the RDI-KLIP reduction with radial profile subtraction that was utilized for spine tracing and surface brightness measurements (Section 4).</figcaption></figure>

## ORCID iDs

Kellen Lawson https://orcid.org/0000-0002-6964-8732

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

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

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

Jean-Charles Augereau https://orcid.org/0000-0002-2725-6415

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

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

Julien Lozi https://orcid.org/0000-0002-3047-1845

Vincent Deo https://orcid.org/0000-0003-4514-7906

Sebastien Vievard https://orcid.org/0000-0003-4018-2569

Nemanja Jovanovic https://orcid.org/0000-0001-5213-6207

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

Nour Skaf https://orcid.org/0000-0002-9372-5056

Thomas Henning https://orcid.org/0000-0002-1493-300X

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

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

Michael W. McElwain https://orcid.org/0000-0003-0241-8956

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

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

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

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

## References

Augereau, J. C., Lagrange, A. M., Mouillet, D., Papaloizou, J. C. B., & Grorod, P. A. 1999, A&A, 348, 557

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

Boccaletti, A., Augereau, J. C., Marchis, F., & Hahn, J. 2003, ApJ, 585, 494

Brandt, T. D. 2018, ApJS, 239, 31

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

Brandt, T. D., Rizzo, M., Groff, T., et al. 2017, JATIS, 3, 048002

Brandt, T. D. 2021, ApJS, 254, 42

Chen, C., Mazoyer, J., Poteet, C. A., et al. 2020, ApJ, 898, 55

Cloutier, R., Currie, T., Rieke, G. H., et al. 2014, ApJ, 796, 127

Currie, T., Brandt, T. D., Kuzuhara, M., et al. 2020b, ApJL, 904, L25

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

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

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

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

Currie, T., Guyon, O., Lozi, J., et al. 2020a, Proc. SPIE, 11448, 114487H

Currie, T., Guyon, O., Tamura, M., et al. 2017, ApJ, 836, L15

Currie, T., Kenyon, S. J., Balog, Z., et al. 2008, ApJ, 672, 558

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

Currie, T., Lisse, C. M., Kuchner, M., et al. 2015b, ApJL, 807, L7

Currie, T., Marois, C., Cieza, L., et al. 2019, ApJL, 877, L3

Currie, T., Rodigas, T. J., Debes, J., et al. 2012b, ApJ, 757, 28

Duchêne, G., Rice, M., Hom, J., et al. 2020, AJ, 159, 251

Engler, N., Boccaletti, A., Schmid, H. M., et al. 2019, A&A, 622, A192

Esposito, T. M., Kalas, P., Fitzgerald, M. P., et al. 2020, AJ, 160, 24

Fitzgerald, M. P., Kalas, P. G., Duchêne, G., Pinte, C., & Graham, J. R. 2007, ApJ, 670, 536

Gagné, J., Mamajek, E. E., Malo, L., et al. 2018, ApJ, 856, 23

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

Gibbs, A., Wagner, K., Apai, D., et al. 2019, AJ, 157, 39

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

Groff, T. D., Chilcote, J., Kasdin, N. J., et al. 2016, Proc. SPIE, 9908, 990800

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

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

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

Kenyon, S. J., & Bromley, B. C. 2008, ApJS, 179, 451

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

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

Lawson, K., Currie, T., Wisniewski, J. P., et al. 2020, AJ, 160, 163

Lisse, C. M., Sitko, M. L., Russell, R. W., et al. 2017, ApJL, 840, L20

Lozi, J., Guyon, O., Jovanovic, N., et al. 2018, Proc. SPIE, 10703, 1070359

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

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

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

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

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

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

Pueyo, L. 2016, ApJ, 824, 117

Ribas, Á., Bouy, H., & Merín, B. 2015, A&A, 576, A52

Schneider, G., Grady, C. A., Hines, D. C., et al. 2014, AJ, 148, 59

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

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

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

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

Steiger, S., Currie, T., Brandt, T. D., et al. 2021, AJ, 162, 44

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

Thalmann, C., Janson, M., Buenzli, E., et al. 2013, ApJL, 763, L29

Wyatt, M. C. 2008, ARA&A, 46, 339