SCExAO/CHARIS Near-infrared Integral Field
Spectroscopy of the HD 15115 Debris Disk

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

► To cite this version:

Kellen Lawson, Thayne Currie, John Wisniewski, Motohide Tamura, Glenn Schneider, et al. SCExAO/CHARIS Near-infrared Integral Field Spectroscopy of the HD 15115 Debris Disk. The Astronomical Journal, 2020, 160 (4), pp.163. 10.3847/1538-3881/ababa6 hal-03255536

HAL Id: hal-03255536

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

Submitted on 3 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>

# SCExAO/CHARIS Near-infrared Integral Field Spectroscopy of the HD 15115 Debris Disk

Kellen Lawson¹, Thayne Currie²,³,⁴, John P. Wisniewski¹, Motohide Tamura⁵,⁶,⁷, Glenn Schneider⁸, Jean-Charles Augereau⁹, Timothy D. Brandt¹⁰, Olivier Guyon²,⁵,⁸,¹¹, N. Jeremy Kasdin¹²,¹³, Tyler D. Groff¹⁴, Julien Lozi², Jeffrey Chilcote¹⁵, Klaus Hodapp¹⁶, Nemanja Jovanovic¹⁷, Frantz Martinache¹⁸, Nour Skaf²,¹⁹,²⁰, Eiji Akiyama²¹, Thomas Henning²², Gillian R. Knapp²³, Jungmi Kwon⁶, Satoshi Mayama²⁴, Michael W. McElwain¹⁴, Michael L. Sitko²⁵, Ruben Asensio-Torres²⁶, Taichi Uyama⁷,²⁷,²⁸, and Kevin Wagner⁸

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

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

³NASA-Ames Research Center, Moffett Boulevard, 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

⁸Steward Observatory, The University of Arizona, Tucson, AZ 85721, USA

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

¹⁰Department of Physics, University of California, Santa Barbara, Santa Barbara, CA, USA

¹¹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

¹⁶Institute for Astronomy, University of Hawaii, 640 North A'ohokū Place, Hilo, HI 96720, 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

²¹Department of Engineering, Niigata Institute of Technology, 1719 Fujihashi, Kashiwazaki, 945-1195, Japan

²²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

²⁴The Graduate University for Advanced Studies, SOKENDAI, Shonan Village, Hayama, Miura, Kanagawa 240-0193, Japan

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

²⁶Department of Astronomy, Stockholm University, AlbaNova University Center, SE-106 91, Stockholm, Sweden

²⁷Infrared Processing and Analysis Center, California Institute of Technology, 1200 East California Boulevard, Pasadena, CA 91125, USA

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

Received 2020 May 28; revised 2020 July 23; accepted 2020 July 31; published 2020 September 11

## Abstract

We present new, near-infrared (1.1–2.4 µm) high-contrast imaging of the debris disk around HD 15115 with the Subaru Coronagraphic Extreme Adaptive Optics (SCExAO) system coupled with the Coronagraphic High Angular Resolution Imaging Spectrograph (CHARIS). The SCExAO/CHARIS resolves the disk down to $\rho \sim 0''2$ ($r_{\text{proj}} \sim 10$ au), a factor of ~3–5 smaller than previous recent studies. We derive a disk position angle of PA ~279°–4–280° and an inclination of i ~85°–3–86.2°. While recent SPHERE/IRDIS imagery of the system could suggest a significantly misaligned two-ring disk geometry, CHARIS imagery does not reveal conclusive evidence for this hypothesis. Moreover, optimizing models of both one- and two-ring geometries using differential evolution, we find that a single ring having a Hong-like scattering phase function matches the data equally well within the CHARIS field of view ($\rho \lesssim 1''$). The disk’s asymmetry, well evidenced at larger separations, is also recovered; the west side of the disk appears, on average, around 0.4 mag brighter across the CHARIS bandpass between 0''25 and 1''. Comparing Space Telescope Imaging Spectrograph (STIS) 50CCD optical photometry (2000–10500 Å) with CHARIS near-infrared photometry, we find a red (STIS/50CCD–CHARIS broadband) color for both sides of the disk throughout the 0''4–1'' region of overlap, in contrast to the blue color reported at similar wavelengths for regions exterior to ~2''. Further, this color may suggest a smaller minimum grain size than previously estimated at larger separations. Finally, we provide constraints on planetary companions and discuss possible mechanisms for the observed inner disk flux asymmetry and color.

Unified Astronomy Thesaurus concepts: Circumstellar disks (235); Debris disks (363)

## 1. Introduction

Gas-poor, dusty debris disks around stars are key laboratories for studying planetary system structure and the late stages of their formation (Wyatt 2008; Hughes et al. 2018). Scattered-light imagery of debris disks around young stars clarifies the disks’ structures and can identify the signatures of sculpting planets and

in situ formation and erosion of icy Kuiper Belt objects (Kalas et al. 2005; Kenyon & Bromley 2008). Further, these studies enable analysis of the composition and scattering properties of the material within the disks, potentially providing reference points for the evolution of the Kuiper Belt (Currie et al. 2015b). High-contrast imaging produced using so-called “extreme
adaptive optics" (exAO) facilities (e.g., SPHERE, Beuzit et al. 2019; GPI, Macintosh et al. 2015; and the Subaru Coronagraphic Extreme Adaptive Optics (SCExAO) system, Jovanovic et al. 2015; Lozi et al. 2018; Currie et al. 2019a) provides the opportunity to study these systems to smaller inner working angles than was possible with conventional AO. This enables the assessment of the disks' scattering phase functions at previously inaccessible angles, as well as placing more significant constraints on the presence of embedded planets.

The debris disk around HD 15115, an F2V star at a distance of $49.0 \pm 0.1$ pc (Gaia Collaboration et al. 2018) and an estimated age of $\sim 10–100$ Myr (e.g., Moór et al. 2006; Rhee et al. 2007; Gagné et al. 2018), could be a particularly good target for studying planetary system structure and the results of the initial formation stages. Discovery optical scattered-light imagery from the Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS) reveals a highly inclined disk with an “extreme” east-west length asymmetry in the HST/ACS F606W bandpass ($\lambda_{\text{pivot}} = 5886$ Å, FWHM = 2325 Å) resolving the eastern extent out to $\sim 7''$, while the western extent is recovered to the edge of the field of view at $12''/38$ (Kalas et al. 2007). Their measurements of the disk's surface brightness on either side show an approximately symmetric brightness at $2''/0$, with the west becoming brighter than the east at larger separations ($\Delta m \sim 1$ at $6''$). Follow-up H-band Keck/NIRC2 AO imagery resolved the disk at $\rho \sim 1''–3''/3$ and revealed a blue (F606W-H) color on both sides and a brightness asymmetry beyond $2''$. The combination of the reported color and its highly inclined orientation led to HD 15115's disk being informally referred to as “the Blue Needle.” Follow-up observations from ground-based AO and space expanded the wavelength range over which HD 15115's disk is resolved and further clarified its properties at separations beyond $1''$. The HST/NICMOS 1.1 µm data revealed evidence of more complicated color gradients and a wavelength dependence for the disk's asymmetry and an angular separation dependence for its colors (Debes et al. 2008).

Subsequent studies found evidence of a bow-like shape in the disk at $\sim 1''–2''$ (Rodigas et al. 2012; Mazoyer et al. 2014; Sai et al. 2015), consistent with a ringlike disk at $\sim 90$ au. Using archival near-infrared (NIR) imaging from the Gemini Observatory, Mazoyer et al. (2014) concluded that, while the system's ring is asymmetrical in brightness, the geometry of the ring itself is symmetric about the parent star. Schneider et al. (2014) reported HST Space Telescope Imaging Spectrograph (STIS) data that significantly improved upon the visible-light telemetry and morphology of the disk. These data showed that the bowing and asymmetry in visible wavelengths continue down to $0''/4$. Additionally, this HST/STIS imaging revealed a previously unseen morphological bifurcation on the east side of the outer disk (Schneider et al. 2014).

Recent results suggest the possible existence of multiple debris ring components. Engler et al. (2019) reported the first exAO observations of HD 15115, consisting of Very Large Telescope/SPHERE total intensity data in the *J* and *H* bands and polarized intensity data in the *J* band. They recovered the disk over stellocentric separations of $\rho \sim 1''/0–5''/5$ and suggested, from peaks in their polarized intensity profiles, the possibility of a distinct noncoplanar inner disk having a fiducial radius of $\sim 1''/3$. Attempting to investigate this using their total intensity imagery, they were unable to reveal conclusive evidence regarding the disk's geometry. MacGregor et al. (2019)

reported 1.3 mm Atacama Large Millimeter/submillimeter Array (ALMA) observations of the system with a synthesized beam size of $0''/58 \times 0''/55$. From these observations, evidence exists for either a distinct inner ring with a radius of $\sim 0''/95$ or a significant gap in the canonical disk at a separation of $1''/2$. Notably, however, they reported a lack of evidence in their data to support the misalignment of the inner disk hypothesized by Engler et al. (2019). Additionally, MacGregor et al. (2019) found an absence of the east-west brightness asymmetry typically reported in previous NIR and optical imagery (e.g., Kalas et al. 2007; Mazoyer et al. 2014). They suggested that the large-grain dust population probed by ALMA was unaffected by the mechanism responsible for the asymmetry reported by other studies over similar separations. As perturbations from planetary-mass companions are often used to explain disk asymmetries and more complex, multiringed disk geometries (e.g., MacGregor et al. 2019), the details of these occurrences in the HD 15115 disk are significant. To better clarify the presence or absence of additional ring components and brightness asymmetries over a wide wavelength range, high-contrast imaging data interior to $1''$, matching the coverage of STIS data from Schneider et al. (2014), are needed.

In this work, we report new NIR scattered-light imagery of the HD 15115 system using the SCExAO system and the Coronagraphic High Angular Resolution Imaging Spectrograph (CHARIS) integral field spectrograph (IFS) in broadband (spanning the NIR *J*, *H*, and *K* bands, 1.13–2.39 µm) mode (Groff et al. 2016). This imagery provides a view of the disk to separations a factor of $\sim 3–5$ smaller than previous recent studies ($\rho \sim 0''/2$). We conduct an analysis of the disk's color in NIR and optical wavelengths by combining CHARIS/IFS data with prior HST/STIS imagery. Through both spine tracing and forward modeling, we investigate the details of the system's geometry and offer constraints for the presence of additional rings or planet companions within CHARIS's $2'' \times 2''$ field of view.

## 2. Data

### 2.1. *Observations*

On 2017 August 30 and 2017 September 7, HD 15115 was observed using the Subaru Telescope's SCExAO paired with the CHARIS/IFS operating in low-resolution ($R \sim 20$) broadband (1.13–2.39 µm) mode and utilizing SCExAO's Lyot coronagraph with a 217 mas diameter occulting spot. CHARIS has a nominal pixel scale of $0''/0164$ pixel$^{-1}$, which has been revised to $0''/0162$ pixel$^{-1}$ (Currie et al. 2018). Both sets of data were collected in angular differential imaging (ADI) mode (Marois et al. 2006), achieving total parallactic angle rotations of $76^\circ$ and $56^\circ$ with total integration times of $t_{\text{int}} = 81$ and 55 minutes, respectively. Each set is made up of 80 individual exposures, with the August 30 images having exposure times of 60.48 s and the September 7 images having exposure times of 41.3 s. Sky frames were obtained for both data sets.

For the September 7 data, the conditions were good, with the “slow” (long coherence time) seeing having an FWHM in the *V* band of $\theta_V \sim 0''/5$. The real-time telemetry data of SCExAO estimated *H*-band Strehl ratios of $\sim 80\%$. The conditions for the August 30 data were comparable. No telemetry data were recorded for the August observation, but the point-spread function (PSF) quality appeared slightly superior by eye.
<table><caption>**Table 1**
PSF Subtraction Algorithm Settings</caption><thead><tr><th rowspan="2">Data</th><th rowspan="2">Parameter Tuning</th><th colspan="4">A-LOCI</th><th colspan="4">KLIP</th></tr><tr><th>g</th><th>N<sub>A</sub></th><th>δ<sub>FWHM</sub></th><th>Δr<sub>sub</sub></th><th>N<sub>PCA</sub></th><th>N<sub>zones</sub></th><th>δ<sub>FWHM</sub></th><th>Δr<sub>sub</sub></th></tr></thead><tbody><tr><td>Aug 30</td><td>Disk</td><td>0.1</td><td>1000</td><td>1.50</td><td>4</td><td>5</td><td>1</td><td>1.25</td><td>4</td></tr><tr><td>Sep 7</td><td>Disk</td><td>0.1</td><td>1000</td><td>1.50</td><td>4</td><td>5</td><td>1</td><td>1.25</td><td>4</td></tr><tr><td>Aug 30<sup>a</sup></td><td>Companion</td><td>1.0</td><td>50</td><td>0.5, 1</td><td>10</td><td></td><td></td><td></td><td></td></tr><tr><td>Sep 7<sup>a</sup></td><td>Companion</td><td>1.0</td><td>50</td><td>0.75, 1</td><td>10</td><td></td><td></td><td></td><td></td></tr></tbody></table>

**Note.** Algorithm settings for A-LOCI and KLIP utilized for PSF subtraction of each of the three data sets. Here g refers to the aspect ratio of the optimization regions, with g < 1 producing azimuthally elongated sections and g > 1 producing radially elongated sections; N<sub>A</sub> refers to the area of optimization regions in units of PSF cores; δ<sub>FWHM</sub> indicates the minimum rotation gap in units of PSF FWHM (for both A-LOCI and KLIP); Δr<sub>sub</sub> gives the radial size of subtraction regions in units of pixels (for both A-LOCI and KLIP); N<sub>PCA</sub> indicates the number of principal components utilized in construction of the model PSF; and N<sub>zones</sub> is the number of subsections into which each KLIP optimization annulus was divided (with a value of 1 corresponding to full annuli). In all A-LOCI reductions, we also truncated the covariance matrix to zero out (normalized) singular values smaller than 1.25 × 10<sup>-6</sup> and constructed a reference PSF from only the 50 most correlated images.

<sup>a</sup> The two entries for the rotation gap refer to the ADI rotation gap and the SDI radial movement gap.

Additionally, we make use of HST/STIS analysis-quality imagery (STIS/50CCD, 2000–10500 Å, λ<sub>pivot</sub> = 5752 Å), originally reported and analyzed in Schneider et al. (2014), to better explore the colors of HD 15115’s disk (see Section 5).

## 2.2. CHARIS Data Reduction

The CHARIS data were extracted from raw CHARIS reads using the CHARIS Data Reduction Pipeline (Brandt et al. 2017). Extracted data take the form of image cubes with dimensions (N<sub>λ</sub>, N<sub>x</sub>, N<sub>y</sub>) = (22, 201, 201) (i.e., 201 × 201 pixel images for each of 22 wavelength channels). Subsequent basic image processing—e.g., sky subtraction, image registration, and spectrophotometric calibration—was carried out as in Currie et al. (2011, 2018).

The PSF subtraction was performed by application of both the Karhunen–Loève Image Projection (KLIP; Soummer et al. 2012) and the Adaptive, Locally Optimized Combination of Images (A-LOCI; Currie et al. 2012, 2015a) algorithms independently.

We performed PSF subtraction with settings geared toward the detection of (a) the HD 15115 debris disk and (b) companions plausibly responsible for sculpting the disk. Table 1 lists our parameter choices for each reduction, with the motivations for these choices summarized hereafter.

Disk detection—The HD 15115 debris disk is oriented nearly edge-on in the plane of the sky (e.g., Kalas et al. 2007; Mazoyer et al. 2014). To detect the disk, we performed PSF subtraction exploiting ADI only, not Spectral Differential Imaging (SDI). For A-LOCI, tuning the geometry of the optimization and subtraction regions—i.e., their relative azimuthal and radial widths—is essential for the recovery of disk flux. Combining the minimum rotation gap with azimuthally elongated regions allowed LOCI coefficients to be computed with less perturbation by the radially extended disk flux while still producing a strong reconstruction of the speckle noise. For KLIP, we performed PSF subtraction in full annuli. To limit self-subtraction of the disk, we imposed a minimum rotation gap of δ ~ 1.25–1.50 λ/D when selecting suitable reference frames for both A-LOCI and KLIP reductions.

Companion detection—To achieve the deeper contrast limits needed to detect faint planets, we used a combination of ADI and then SDI (on the ADI residuals) using A-LOCI following Currie et al. (2018). For the ADI component, optimization regions were a factor of 20 smaller (50 PSF footprints), and the rotation gap was reduced to δ = 0.5–0.75. For a given section of the science image, up to the 50 most correlated sections from

the reference image library were used to build a reference PSF (with the number of available images depending on the portion of the 80 exposures satisfying the minimum rotation gap requirement for the section). For the SDI component, the optimization zone covers an annular region with the same depth Δr<sub>sub</sub> as the subtraction zone, but the smaller annular wedge-shaped subtraction zone is masked.

## 2.3. Results

Both PSF subtraction techniques yield strong detections of the disk to ρ ~ 0.′15–0.′25 in CHARIS data (Figures 1 and 2 and), improving upon the 0.′4 angular separation achieved with optical HST/STIS data (Schneider et al. 2014). The CHARIS data mark a substantial improvement over previous ground-based, NIR scattered-light imaging of the disk, with conventional AO data limited to ρ ≳ 1" and exAO imaging from Engler et al. (2019) detecting the disk exterior to ρ ~ 0.′75–1". This improvement is owed in part to the significant field rotation achieved, with field rotations of 76° and 56° for our two sets of observations versus, e.g., 23° for the data from Engler et al. (2019).

The quality of the detection varies from the *J* band, where the disk detection is contaminated by residual speckles, to the *H* and *K* bands, where the images are free of strong residuals exterior to 0.′25 (Figure 3). The detection in the broadband images is strongest in the August 30 data (especially for the A-LOCI reduction), with a signal-to-noise per resolution element (SNRE) along the trace of the disk of ~5–7 for most regions exterior to 0.′25<sup>29</sup> (Figure 4). In the *H* and *K* bands, the disk detection is strong, generally achieving an SNRE of ~3–5 along the disk and peaking around 5.6 to the west; the *J*-band detection of the disk is considerably weaker (though still definitive), with an SNRE of ~2–3 over the same regions and peaking around 4 in the west (see Figure 5).

CHARIS imagery reveals a strongly asymmetrically scattering disk whose maximum intensity is unambiguously offset (with a projected semiminor axis of ~0.′12) from the system’s major axis throughout the ~2" × 2" field of view. This indicates a view of the system entirely inside the bow-like feature originally described by Rodigas et al. (2012). For the assumption of preferential forward scattering (Hughes et al. 2018), the

<footnote n="29"><sup>29</sup> In the S/N calculation, we use a software mask to reduce the amount of disk signal included in the noise estimation. This approach increases the finite-element correction penalty (Mawet et al. 2014), so the gain in S/N is small.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The HST/STIS imagery of HD 15115 (originally reported in Schneider et al. 2014) with lower and upper inset images corresponding to CHARIS imagery and the same HST/STIS imagery scaled to the CHARIS field of view. The CHARIS imagery presented is the average of the results for the A-LOCI reductions of the August 30 and September 7 data using settings for disk detection (see Table 1). The CHARIS image's central mask has a radius of $0''_{15}$. The STIS, STIS inset, and CHARIS images have the same orientation and are depicted with linear display stretches spanning 0–0.02, 0–4.0, and 0–5.21 mJy arcsec$^{-2}$, respectively.</figcaption></figure>

brighter (~northern) side of the disk observed clearly in our data would be presumed as the near side. Signal to the west and slightly south of the center in Figures 2 and 4 (annotated in the latter) may constitute marginal detections of the disk's dimmer (presumably far) side, which has been recovered in previous ground-based imagery (e.g., Mazoyer et al. 2014; Engler et al. 2019).

### 3. Disk Morphology

To make estimates of the disk's geometric parameters, we seek to identify the position of peak brightness along the disk's bright (~northern) edge for each of the wavelength-collapsed (CHARIS broadband) final images resulting from both A-LOCI and KLIP reductions of the August 30 and September 7 data (four images in total). We begin by rotating the image based on the PA reported in Engler et al. (2019; PA = 278°.9 ± 0°.1, which is then adjusted to account for the CHARIS PA correction discussed below) to orient the disk's major axis along the x-axis.<fnref n="30" /> For each x-direction integer pixel position, we identify the approximate y-position of the brightness peak by taking the median of the locations of the brightest 5% of pixels in that column. To determine the precise peak position, we then fit a Lorentzian profile to the array of flux values within 15 pixels (~0''25) of the approximate peak, seeking the y-position of the profile's peak (and taking the standard error from this fit as the uncertainty). For these fits, each flux value is weighted by the inverse of the corresponding noise levels from noise maps (described in Section 4.1). The y-position of the spine as a function of x-position is computed this way for each of the four images being analyzed. The spine profiles for the

four images are then combined by taking the weighted average of the four values at each x-position as the nominal average spine y-position, with uncertainty corresponding to the standard error for a weighted average.<fnref n="31" /> The resulting average spine profile is then fit with an ellipse (which is centered on the star), described by projected semimajor axis (a), nominal position angle (PA₀, which accounts for the initial image rotation applied previously), and inclination (i). Our best fit (over the region $0''_{25} \le \rho \le 1''_{0}$) is achieved for an ellipse with parameters a = 80.1 ± 3.3 au, PA₀ = 277°.82 ± 0°.05, and i = 85°.76 ± 0°.22. The nominal position angle measurement is then corrected for the CHARIS PA offset of −2°.20 ± 0°.27³² (Currie et al. 2018). This results in a final measurement of PA = 280°.02 ± 0°.27. The average spine positions and the best-fit ellipse are visualized in Figure 6.

The position angle measured by this methodology falls above recent measurements made by MacGregor et al. (2019; PA = 278° ± 1°) and Engler et al. (2019; PA = 278°.9 ± 0°.1). The measured inclination is consistent with the values of both works (MacGregor et al. 2019 measured i = 86°.3 ± 0°.4, and Engler et al. 2019 measured i = 85°.8 ± 0°.7). Evaluating $\chi^2_ν$ for the spine parameters from Engler et al. (2019) with our data suggests that the difference we measure is significant, with these parameters giving $\chi^2_ν = 36.6$ versus our best fit of $\chi^2_ν = 1.1$. Given that our measurements of the disk are made in the region of $\rho \sim 0''_{2-1}''_{0}$, it is possible that we are measuring the overlapping signal of the canonical outer ring and the inner ring proposed in both works. If the disk profile observed is the result of

$$ \sigma'_{i} = \sqrt{\frac{\sum_{j=1}^{4} (y_{ij} - \mu_{i})^{2} \cdot w_{ij}}{\sum_{j=1}^{4} w_{ij}}}. $$

<footnote n="30">³⁰ The stated PA is assumed here only for the purpose of this initial rotation, which is carried out to simplify the fitting of the spine across the narrowest part of the disk, where the peak will appear sharpest. This additionally allows meaningful measurements of the disk's projected FWHM for later use (see Section 5). Rotating the images based on the 1σ upper and lower limits from Engler et al. (2019) instead (279°.0 and 278°.8, respectively) changes the eventual PA measurement from this procedure by ~0°.01, likely largely as a result of the rotation interpolation itself.</footnote>

<footnote n="31">³¹ For the jth image's ith x-axis spine position, call the corresponding y-axis spine position $y_{ij}$ with an associated fit y-position uncertainty $\sigma_{ij}$. The standard error for the corresponding weighted average, $\mu_i$, with weight $w_{ij} = 1/\sigma_{ij}^2$, is then</footnote>

<footnote n="32">³² Hereafter, any values of PA presented (e.g., in the case of disk modeling in Section 4.5) are already corrected for this PA offset.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Wavelength-collapsed results for the August 30 and September 7 observations following PSF subtraction using either A-LOCI or KLIP techniques with settings for disk detection (see Section 2.2). The central mask in each subplot has a radius of 0"/15, and the image has a linear display stretch spanning -0.61 to 6.10 mJy arcsec⁻². In the August 30 products, a plainly visible flux enhancement appears just beyond the inner software mask to the southwest, but it is not evident in the September 7 data. This feature is likely residual speckle noise (likewise for the similar feature to the southeast in the September 7 imagery).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The A-LOCI PSF-subtracted imagery for August 30 using settings for disk detection, with wavelength channels combined to produce images comparable to the *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. The central mask in each subplot has a radius of 0"/15. Images are displayed with linear stretches spanning 0–24.13 (*J* band) or 0–4.39 (*H* and *K* bands) mJy arcsec⁻².</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Maps of SNRE for August 30 and September 7 data following PSF subtraction using A-LOCI or KLIP techniques with settings for disk detection (see Section 2.2). The central mask in each subplot has a radius of 0.′′15. Possible signal from the disk's fainter side is indicated in the subplot for the KLIP reduction of August 30 data, peaking at an SNRE of ~4. This feature is also visible in the A-LOCI reduction of the same data, albeit at a lower SNRE (~2.5).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Maps of SNREs for the *J*, *H*, and *K* bands from A-LOCI PSF subtraction of August 30 data (see Section 2.2). The central mask in each subplot has a radius of 0.′′15.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Stellocentric separation of disk spine fits along the major and minor axes. The spine positions are used for ellipse fitting and evaluating disk geometry. Points with errors correspond to the weighted average position from all four utilized CHARIS broadband images (see Section 3). The blue arc depicts the best-fit ellipse solution for the data shown and provides the disk major and minor axes against which the spine positions are plotted here.</figcaption></figure>

an architecture featuring a distinct noncoplanar inner ring (a possibility suggested by Engler et al. 2019), we should expect the fit values to be skewed somewhere between those of the inner and outer component. The difference in measured PA could also be explained by an imprecise calibration for either instrument (or

both). However, the CHARIS PA calibration utilized was performed using data collected only a day after our September observations (Currie et al. 2018); it should provide a reasonable assessment of the PA calibration for our data. A reevaluation of the CHARIS north PA and pixel scale using additional data
obtained at additional epochs reaffirms these results (T. Currie 2020, in preparation).

Attenuation of disk flux during PSF subtraction can also have an effect on the measured position of the spine and thus on the derived parameters as well. However, we note that carrying out the aforementioned measurements on attenuated models with known PAs from our forward-modeling procedure (see Section 4.1) indicates that this effect is small. For example, for our best-fit one-ring model with true PA = 279°.8, we measure PA = 279°.7 ± 0°.3.

# 4. Modeling the Debris Disk of HD 15115

## 4.1. Disk Forward Modeling

We implement a strategy of forward-modeling synthetic disks, as described in Currie et al. (2018, 2019b), to investigate the details of HD 15115’s debris disk. In this approach, coefficients (for A-LOCI) or Karhunen–Loève (KL) modes (for KLIP) retained from the science data reduction are applied to image cubes containing only the signal of a model disk that has been rotated to reproduce the array of observed position angles and convolved with the instrumental PSF.

For this procedure, we consider all three sources of flux annealing described by Pueyo (2016): oversubtraction (speckle noise being subtracted from the disk signal), direct self-subtraction (resulting from inclusion of the disk signal in basis vectors), and indirect self-subtraction (resulting from the perturbation of basis vectors by the disk signal). For our disk reductions, indirect self-subtraction is expected to be the smallest of the three terms. For our KLIP reductions, we retained a small number of KL modes compared to the total number of reference images available. For A-LOCI, we adopted a large optimization area ($N_A$ = 1000 PSF cores) that is azimuthally elongated, in contrast to the nearly edge-on debris disk. For both reductions, we adopted a large rotation gap of $\delta$ = 1.25–1.5 PSF cores. Thus, oversubtraction and self-subtraction likely dominate over indirect self-subtraction.

Once processed this way, the result for a given model can be compared to the result for science data to assess the relative strength of the model. The procedure for calculating $\chi_\nu^2$ that we implement is as previously described in Goebel et al. (2018) but briefly summarized here. First, each value in the model image, $f_{\text{model}}$, and the science image, $f_{\text{obs}}$, is replaced with the sum of values within an FWHM-sized aperture. Following this, a finite element-corrected noise map is computed from the science image as described in Currie et al. (2011), additionally utilizing a software mask as described above in Section 2.3. Then, the model image is rescaled to minimize the inverse variance-weighted residuals with the science image in a region of interest. The difference of the scaled model and the science image, weighted by the noise map, is squared to create a $\chi^2$ map. This map is then binned to the size of the instrumental PSF. Finally, the reduced $\chi^2$ metric is computed from this as $\chi_\nu^2 = \chi^2/\nu = \nu^{-1}\sum_i^N(f_{i,\text{obs}} - f_{i,\text{model}})^2/\sigma_i^2$, where the degree of freedom, $\nu$, is given by the difference between the number of bins in the binned optimization region, $N$, and the number of free parameters in the model, $M$ (Thalmann et al. 2013).

The region of interest considered is a rectangular box of unbinned dimensions 120 pixels × 30 pixels (~2.0′′0 × 0.5′′) centered on the star and oriented to fall along the disk’s

approximate major axis. The region interior to 6 pixels (~0.1′′1) is excluded. For $\rho \lesssim 6$ pixels, not-a-number (NaN) values begin to appear when no reference frames can meet the minimum rotation gap requirement during PSF subtraction. The region of interest described is overlaid as a white rectangle in Figures 9 and 10.

We delineate models that are acceptably consistent with our observations as in Thalmann et al. (2013), i.e., those having $\chi_\nu^2 \le \chi_{\nu,\text{min}}^2 + \sqrt{2/\nu}$.

## 4.2. Model Debris Disks

The forward-modeling procedure outlined above was applied to model debris disks generated using a version of the GRaTeR software (Augereau et al. 1999). The formalism and assumptions of the models are detailed in Augereau et al. (1999) but summarized briefly hereafter. The models assume an optically thin disk with a radial dust grain distribution described by a smooth combination of two power laws and a vertical distribution described by an exponential function. For simplicity, and in line with the analysis of Engler et al. (2019), we set the vertical exponential distribution to be Gaussian in shape ($\gamma$ = 2) and restrict the flaring of the disk to be linear ($\beta$ = 1). To describe the angular distribution of scattered light, it is common to adopt the Henyey–Greenstein (HG) phase function (Henyey & Greenstein 1941). The HG phase function is parameterized by a single variable, the asymmetry parameter $g$, defined as the average of the cosine of the scattering angle, weighted by the (normalized) phase function, over all directions. However, as noted in Hughes et al. (2018), this formalism is not physically motivated and may introduce misleading results. Moreover, a simple HG phase function fails at reproducing the surface brightness profile we observe for HD 15115’s disk (see Appendix A). Hong (1985) implemented a linear combination of three HG phase functions to describe the observed angular distribution of scattered light for zodiacal dust. To limit model freedom and avoid nonphysical solutions, we adopted the phase function of Hong (1985) with the asymmetry parameters and corresponding weights identified therein: $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. Though allowing the asymmetry parameters and weights to vary during the exploration of disk models may improve the eventual result, it would also massively increase the size and complexity of the parameter space. As early testing showed that the empirically derived parameters of Hong (1985) reproduced our observed disk surface brightness quite closely, we chose to adopt them as-is to allow a more thorough exploration of the remaining disk parameters. See Appendix A for comparisons of scattering phase functions with our data, including simple HG phase functions using commonly reported asymmetry parameters for HD 15115’s disk.

The model geometries investigated fall under two archetypes. The first is a single-ring model, defined by six parameters:

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

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

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

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

5. PA, the position angle of the disk in degrees; and

6. *i*, the inclination of the disk in degrees.
The second is a two-ring model, taken to be the linear superposition of two single-ring models<fnref n="33" /> plus an additional parameter: $F_{max,2}/F_{max,1}$, the ratio of the peak flux of ring 2 to that of ring 1. This results in a model described by 13 parameters (allowing inclination and position angle to differ between the inner and outer disks). To better explore the parameter space of the inner ring, we reduce these to seven parameters by setting the well-studied outer ring's parameters to approximately match the ones identified by prior studies of the disk (e.g., Engler et al. 2019): $R_{0,1} = 96$ au, $\alpha_{in,1} = 2$, $\alpha_{out,1} = -3$, $H_{0,1}/R_{0,1} = 0.03$, PA$_1$ = 278°,9, and $i_1$ = 86°0.

For the purpose of $\chi^2$ calculation, the overall scaling factor applied to the model (see Section 4.1) is considered to be an additional free parameter for the model. This results in $M = 7$ for the one-ring model and $M = 8$ for the two-ring model.

### 4.3. Model Optimization Using Differential Evolution

Though HD 15115's disk has been studied extensively in the region beyond $1''0$, the small separations observed with CHARIS provide a look at the disk in the $0''2-1''0$ regime, where the parameters of the posited inner ring could potentially be studied in much greater detail. However, given the nearly edge-on orientation of the system, a parameter space with significant model degeneracies and multiple local minima is possible. A broad but detailed search of the parameter spaces outlined in Section 4.2 is necessary to offer a meaningful assessment of any such degeneracies and ensure that a unique and globally optimum solution is identified. A grid search for the seven-parameter two-ring model quickly reaches an intractably large size; e.g., a coarse grid examining only five values of each parameter would require 78,125 models to be propagated through the time-consuming forward-modeling procedure (typically ~5–20 minutes per model). While Markov Chain Monte Carlo (MCMC) techniques are commonly used for similar purposes (e.g., MacGregor et al. 2019), they can become trapped in local minima, and their results can be dependent on the initialization. Moreover, MCMC exploration typically requires a number of model evaluations that is effectively unapproachable for ADI forward-modeling procedures comparable to ours (e.g., the MCMC procedure of MacGregor et al. 2019 evaluates ~10⁶ models).

Instead, we make use of the differential evolution (DE) algorithm (Storn & Price 1997) to explore possible solutions for the parameters of each model. It requires no initial assumptions about a solution beyond the boundaries within which to explore and is capable of efficiently probing large, correlated parameter spaces by quickly evolving a population of trial solutions away from regions that offer inferior solutions and allowing population members to move between local minima. Though DE has not seen widespread use in the optimization of scattered-light disk models, it has been used elsewhere in the study of astrophysics with noteworthy efficacy, e.g., to explore optimal capture trajectories for Jovian orbiters by the European Space Agency (Labroquère et al. 2014), identify edge-on galaxies with abundances of extra-planar dust (Shinn 2018), or search for flaring stars in sparsely sampled time-series data (Lawson et al. 2019). While DE does not enable the robust determination of parameter likelihood

distributions in the same manner as MCMC, it is extremely effective at quickly reaching a global solution with little to no tuning of algorithm control variables or specific experience with the algorithm itself and with relatively few function evaluations needed (e.g., see comparison benchmarks in Storn & Price 1997). Additionally, the relative simplicity of DE makes it trivial to add into existing frameworks (the C-style psuedocode for the algorithm presented in Storn & Price 1997 requires only 19 lines of code). For groups currently exploring model parameters using grid searches, in all but the rarest cases, DE will tend to identify a superior final model while evaluating many fewer models overall. For groups that might be interested in adopting this technique, we include a simple Python implementation of DE, whose procedure is outlined below, in Appendix B.

In the DE procedure, we initialize a population of $N_{pop}$ random model parameter sets restricted to fall between boundaries set for each parameter. For each type of model, $N_{pop}$ is set to be 10 times the number of free parameters<fnref n="34" /> (60 for the single-ring model and 70 for the two-ring model). The initial model population is run through the forward-modeling routine to evaluate the fitness of the models ($\chi^2$). Following this, a mutation ($v_i$) for each member of the population ($x_i$) is created by adding a scaled difference of the parameters of two distinct, random population members ($x_j, x_k$) to the parameters of the current best solution ($x_{best}$): $v_i = x_{best} + m(x_j - x_k)$ (this “strategy” is called “best/1/bin” in the notation of Storn & Price 1997). For our purposes, the value of the mutation constant, $m$, is randomly selected in the range [0.5, 1.0] for every generation.<fnref n="35" /> From a given mutation, a trial replacement ($u_i$) is formed by setting $u_i = x_i$ and allowing the probability $P = 0.7$<fnref n="36" /> for each parameter value in $u_i$ to be replaced by the corresponding value from $v_i$. The fitness of the set of $N_{pop}$ trial replacements is then evaluated using the same forward-modeling routine. Finally, each member of the population, $x_i$, is replaced by its corresponding trial vector, $u_i$, if the fitness of the trial vector is superior to that of the population member. The procedure of creating and evaluating trial replacements for the population is repeated until the population becomes stagnant or converges to a single solution.

Fitness is evaluated at each stage by computing the combined $\chi^2$ for both A-LOCI reductions as

$$ \chi^2 = \frac{\chi_1^2 + \chi_2^2}{N_1 + N_2 - M} = \frac{1}{\nu} \sum_{i=1}^{2} \chi_i^2, $$

where subscripts 1 and 2 correspond to the values for A-LOCI reductions of the August 30 and September 7 data. Once the DE procedure is completed, the subset of models meeting the threshold $\chi^2 \le \chi^2_{min} + \sqrt{2/\nu}$ is then propagated through the forward-modeling procedure for the KLIP reductions. The final

<footnote n="33">³³ For example, by coadding the synthetic images for the individual models, as in Boccaletti et al. (2019). This assumes that the rings are sufficiently optically thin that single scattering dominates over multiple scattering.</footnote>

<footnote n="34">³⁴ Storn & Price (1997) suggested 5–10 population members per free parameter, though larger values are often used in recent implementations as well, e.g., the default value of 15 per free parameter in the implementation from the Python package SciPy (Virtanen et al. 2020).</footnote>

<footnote n="35">³⁵ Storn & Price (1997) introduced a single mutation constant and suggest a value of 0.5. The “dithered” mutation constant implemented here is adopted from the default setting in the SciPy module for Python, where the authors suggested that a dithered mutation constant will typically speed convergence substantially (Virtanen et al. 2020).</footnote>

<footnote n="36">³⁶ A crossover probability, $P$, anywhere from 0.1 to 0.9 is recommended in Storn & Price (1997), and $P = 0.7$ is adopted from the default value in the SciPy module for Python.</footnote>
<table><caption>**Table 2**
One-ring Model Optimization Results</caption><thead><tr><th>Parameter</th><th>Bounds</th><th>Best</th><th>Acceptable<sup>a</sup></th></tr></thead><tbody><tr><td>R<sub>0</sub> (au)</td><td>86.0–106.0</td><td>93.2</td><td>87.7–99.6</td></tr><tr><td>α<sub>in</sub></td><td>2.0–10.0</td><td>2.0</td><td>2.0–3.6</td></tr><tr><td>α<sub>out</sub></td><td>−10.0 to −2.0</td><td>−9.6</td><td>−10.0 to −7.5</td></tr><tr><td>H<sub>0</sub>/R<sub>0</sub></td><td>0.01–0.05</td><td>0.05</td><td>0.04–0.05</td></tr><tr><td>PA (deg)</td><td>277.0–281.0</td><td>279.8</td><td>279.4–280.5</td></tr><tr><td>i (deg)</td><td>85.0–87.0</td><td>85.3</td><td>85.3–85.7</td></tr><tr><td>χ<sup>2</sup><sub>ν,tot</sub></td><td>...</td><td>1.166</td><td>≤1.196</td></tr></tbody></table>

**Notes.** Optimization bounds, best-fitting value, and range of acceptable values for each varied parameter of the single-ring model following propagation through A-LOCI and KLIP forward-modeling procedures for the August 30 and September 7 data. Here $\chi^2_{\nu,tot}$ indicates the combined measure for all four reductions (see Section 4.3).

<sup>a</sup> These ranges give the smallest and largest value of each parameter that resulted in an acceptable solution. Given the possibility of complicated correlations between parameters and the lack of perfect sampling, it cannot be stated conclusively that every value within these ranges can produce an acceptable solution.

<table><caption>**Table 3**
Two-ring Model Optimization Results</caption><thead><tr><th>Parameter</th><th>Bounds</th><th>Best</th><th>Acceptable</th></tr></thead><tbody><tr><td>R<sub>0.2</sub> (au)</td><td>20.0–70.0</td><td>40.9</td><td>36.3–62.1</td></tr><tr><td>α<sub>in,2</sub></td><td>2.0–10.0</td><td>3.0</td><td>2.0–8.2</td></tr><tr><td>α<sub>out,2</sub></td><td>−10.0 to −2.0</td><td>−5.3</td><td>−9.4 to −4.6</td></tr><tr><td>H<sub>0.2</sub>/R<sub>0.2</sub></td><td>0.01–0.10</td><td>0.03</td><td>0.02–0.08</td></tr><tr><td>PA<sub>2</sub> (deg)</td><td>277.2–285.0</td><td>281.6</td><td>280.2–283.8</td></tr><tr><td>i<sub>2</sub> (deg)</td><td>70.0–87.3</td><td>80.1</td><td>79.0–85.0</td></tr><tr><td>F<sub>max,2</sub>/F<sub>max,1</sub></td><td>0.5–10.0</td><td>2.0</td><td>1.2–2.7</td></tr><tr><td>χ<sup>2</sup><sub>ν,tot</sub></td><td>...</td><td>1.151</td><td>≤1.181</td></tr></tbody></table>

**Note.** Same as Table 2 but for two-ring model optimization.

best model is taken from among these as the one that minimizes the combined $\chi^2_{\nu}$ for all four reductions,

$$ \chi^2_{\nu,tot} = \frac{\sum_{i=1}^{4} \chi_i^2}{(\sum_{i=1}^{4} N_i) - M} = \frac{1}{\nu} \sum_{i=1}^{4} \chi_i^2, $$

with $i = 3$ and 4 indicating the two KLIP reductions. While it may be preferable to evaluate each set of trial models for all four reductions, this procedure cuts the total model optimization time nearly in half by assuming that the overall best model will be contained within the “acceptable” fitness bounds of the first two reductions. See the second column of Tables 2 and 3 for the bounds adopted for each parameter.

### 4.4. Model Results

**One-ring models**—The DE procedure for single-ring models yields an optimal combined fit for the two A-LOCI reductions of $\chi^2_{\nu} = 1.145$. A visualization of the full sample of models explored using DE is provided in Figure 7. Some parameters ($\alpha_{in}$ and $\alpha_{out}$) have converged to the boundaries. While this could indicate that the bounds are too restrictive, we note that

the adopted boundaries include all values of these parameters explored by other recent studies of HD 15115 utilizing GRaTeR (e.g., Mazoyer et al. 2014 used models with $-6 \le \alpha_{out} \le -4$ and $2 \le \alpha_{in} \le 10$, and Engler et al. 2019 used models with $-8 \le \alpha_{out} \le -2$ and $2 \le \alpha_{in} \le 10$). This result is discussed in more detail in Section 4.5. However, we note here that these particular parameters are unimportant for our overarching conclusions.

From the 660 models evaluated, 13 resulted in acceptable values of $\chi^2_{\nu}$ ($\chi^2_{\nu} \le 1.188$, for $\nu = 1077$ with the two A-LOCI reductions included). After forward modeling this subset for the remaining reductions, the final best model results in $\chi^2_{\nu,tot} = 1.166$ with a revised acceptable limit of $\chi^2_{\nu,tot} \le 1.196$ (for $\nu = 2161$ when all four reductions are included). The parameters for the best overall model and ranges of acceptable parameters are included in Table 2. It should not be assumed that an acceptable solution can be produced for any arbitrary combination of parameter values falling within the acceptable ranges. A given value included in the acceptable range indicates that the value, paired with specific values of the other parameters, produces a model meeting the given threshold for acceptability.

**Two-ring models**—Following the initial DE procedure, the optimal fit for two-ring models results in $\chi^2_{\nu} = 1.169$ for the two A-LOCI reductions. A visualization of the full sample of models explored using DE is provided in Figure 8. The best solution at this stage has an inner-ring radius that falls at the upper boundary of the allowed values. We note, however, that the selected bounds include the best-fit inner rings of both MacGregor et al. (2019; $R_{0,2} \sim 48$ au) and Engler et al. (2019; $R_{0,2} = 64$ au).

From the 910 models evaluated here, 20 resulted in acceptable values of $\chi^2_{\nu}$ ($\chi^2_{\nu} \le 1.212$ with $\nu = 1076$). After KLIP forward modeling for this subset, the final best model resulted in $\chi^2_{\nu,tot} = 1.151$ with a revised acceptable limit of $\chi^2_{\nu,tot} \le 1.181$ ($\nu = 2160$). Of the 20 models evaluated for all four reductions, 12 pass this revised threshold for acceptability. For clarity, we note that the best model from the initial DE/A-LOCI-only reductions (the model whose parameters are indicated in Figure 8) is ultimately excluded by this final acceptability threshold. This model simply ends up being a worse explanation for the KLIP results than other acceptable models, which pushes its combined score up sufficiently to be eliminated (with $\chi^2_{\nu,tot} = 1.183$). The parameters for the best overall model and ranges of acceptable parameters are included in Table 3.

The models before and after forward-modeling, as well as the residuals with the data and $\chi^2$ maps are shown for the best one- and two-ring models in Figure 9 and Figure 10 respectively.

### 4.5. Modeling Discussion

With the optimal one- and two-ring models identified producing comparable $\chi^2_{\nu}$ metrics of 1.166 and 1.151, respectively, both geometries appear statistically consistent with our data, making it difficult to rule out either scenario. However, our exploration of the model parameter spaces using DE allows us to place some disambiguating constraints.

A few noteworthy observations can be made regarding the one-ring model results. First, the corner plot of our DE procedure (Figure 7) reveals a preference for a very slow
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Variation of the corner plot for optimization of the single-ring model to the A-LOCI reductions of HD 15115's 2017 August 30 and 2017 September 7 data. Each off-diagonal plot visualizes solutions as a function of two of the parameters, with each bin colored according to the quality of the best fit achieved with the values of the two parameters in that range (and any values of the other parameters). Darker bins indicate smaller values of $\chi_v^2$, where $\chi_v^2$ refers to the combined metric for the A-LOCI reductions of both data sets (see Section 4.3). Diagonal elements provide a one-dimensional view of each of the parameters, indicating the lowest $\chi_v^2$ value (y-axis) achieved for the binned range of the given parameter (x-axis). For each parameter, the area within the bounds provided in Table 2 is divided into six equally sized bins. The best-fit solution values are indicated by black crossshairs, and their values are given above the corresponding one-dimensional subplot. The threshold fitness for acceptable solutions with A-LOCI reductions ($\chi_v^2 \le \chi_{\text{min}}^2 + \sqrt{2/\nu}$) is indicated by a horizontal gray dashed line in each diagonal plot; bars that end below this threshold line resulted in models meeting the acceptable fitness criteria for some values of the other parameters.</figcaption></figure>

density change interior to $R_0$ (small $\alpha_{\text{in}}$), with a very rapid change exterior to $R_0$ (large $\alpha_{\text{out}}$). For both of these parameters, the optimization converges at or very near the boundaries (2 and -10, respectively). While this could suggest true values for these parameters beyond the boundaries we have set, this seems unlikely given the results of previous studies of HD 15115's disk; for example, the observed extent of the disk with wider fields of view is inconsistent with a radial density profile having $\alpha_{\text{out}} \sim -10$. We note, however, that models more consistent with prior results produce acceptable fit metrics as well (see Table 2). Further, given our narrow field of view and the disk's high inclination, we should expect that our results are not as sensitive to changes in these particular parameters

anyway. Additionally, the PA identified is seemingly distinct from values typically found in previous studies of the disk, with our procedure finding $\text{PA} = 279.8^{+0.737}_{-0.4}$ compared to $\text{PA} = 278.9 \pm 0.1$ from Engler et al. (2019). However, the value identified here is consistent with the value of $280.02 \pm 0.27$ identified from ellipse fitting of the spine in Section 3.

From the schematic of our best two-ring model (Figure 10, panel a), we see that our data are best explained by an inner ring with a projected semiminor axis similar to that of the outer

<footnote n="37">³⁷ Uncertainties here are roughly approximated as simply the upper and lower limits for acceptable models as presented in Table 2.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Same as Figure 7 but for optimization of the two-ring model.</figcaption></figure>

ring ($b_{\text{proj}} \sim 6.7$ au for $R_0 = 96$ au and $i = 86^\circ$). Looking at the corner plot for the two-ring optimization (Figure 8), perhaps the strongest apparent correlation occurs between the inner ring's radius (equivalently, its projected semimajor axis, *a*) and its inclination. Noting that the projected semiminor axis, *b*, of the inner ring is related to its inclination and projected semimajor axis by $b_{\text{proj}} = a_{\text{proj}} \cdot \cos(i)$, the correlation in the $R_{0,2}$ versus $i_2$ subplot falls very nearly along the line corresponding to $b_{2,\text{proj}} = 6.7$ au (see Figure 11). In fact, for the full set of acceptable two-ring models following DE, all have projected semiminor axes between 4.8 and 8.0 au, while values from 0.9 to 23.9 au are permitted by our parameter bounds. The revised set of acceptable models following modeling of all four CHARIS reductions reduces this range even further, to 5.4–8.0 au. Combined with our analysis of the disk's spine and surface brightness profiles, our results appear to suggest a lack of any statistically significant distinct inner-ring spine.

Overall, the results of our modeling procedure can be interpreted in a number of ways.

1. From the strong preference for $b_{1,\text{proj}} \sim b_{2,\text{proj}}$, a distinct inner ring exists, but its brightest features happen to roughly line up along our line of sight with the canonical outer ring. This would result in a two-ring geometry that is statistically indistinguishable from one with a single ring. A distinct inner ring that is coplanar (or nearly coplanar) with the outer ring is consistent only as the inner ring's radius approaches that of the outer, where $b_{2,\text{proj}}$ can near $b_{1,\text{proj}}$ while maintaining a matching inclination.

2. In the CHARIS field of view (~0"/2–1"/0), either the outer or inner ring is substantially brighter, such that the other is not recovered in our data. The differing PA identified for the one-ring model compared to literature could suggest that we are seeing the latter. However, we remark that our optimization procedure for two-ring
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Lowest $\chi^2_{\nu,\text{tot}}$ single-ring model ($\chi^2_{\nu,\text{tot}} = 1.166$) identified by the procedure outlined in Section 4.1, depicted for our KLIP reduction of August 30 data. The model parameters utilized can be found in Table 2. In all panels, the white rectangle indicates the “region of interest” for the purpose of $\chi^2$ calculation, with the white circle indicating the inner exclusion radius. (a) Input disk model convolved with the instrumental PSF. A schematic of the disk model is overlaid as a dashed blue ellipse having the same radius, inclination, and PA as the model disk. (b) Disk model after application of the forward-modeling procedure for the 2017 August 30 KLIP reduction, displayed exactly the same as the images of Figure 2 (linear stretch over -0.61 to 6.10 mJy arcsec$^{-2}$). (c) Residuals for the 2017 August 30 KLIP data product (upper left image in Figure 2) and the processed disk model, also displayed as Figure 2 (linear stretch over -0.61 to 6.10 mJy arcsec$^{-2}$). (d) Unbinned $\chi^2$ map for this model displayed with linear stretch over 0–19.46 in $\chi^2$ (roughly equivalent to 0–0.036 in $\chi^2_\nu$).</figcaption></figure>

models allowed values of $F_{\text{max},2}/F_{\text{max},1}$ as large as 10; if the underlying system is well described within CHARIS’s field of view as an extremely faint outer ring with a misaligned bright inner ring, solutions with large ring flux ratios should have manifested.

3. Perhaps the system is truly better described as a single ring. The tendency of the single-ring DE procedure toward parameter values that are seemingly at odds with prior results could be caused by inaccuracies in assumptions made by our models. For example, if the true scattering phase function departs slightly from the Hong phase function assumed, the disk is nonnegligibly eccentric, or the disk features a nonlinear flaring profile, models matching the true underlying parameters may not coincide with the minimum $\chi^2_\nu$ in our analysis. In this case, even if no second ring exists, a second-ring model might serve to mitigate these inaccuracies sufficiently to result in competitive fitness metrics.

4. The disk has a geometry distinct from any probed here and that is difficult to diagnose as a result of the nearly edge-on orientation (e.g., a debris disk with spiral arms or significant warping).

We also point out the presence of the very small separation ($\rho \sim 0.25''$) residual signal that is not fit by our models (bright residuals appearing just beyond the software mask in panel (c) of Figures 9 and 10). While this feature could be evidence of the ~4 au “warm-dust” disk suggested by Moór et al. (2011) through spectral energy distribution fitting, perhaps the more likely explanation is that it is simply residual speckle noise. While “aggressive processing” is often cited as a cause of spurious features in debris disks (Duchene et al. 2020), poor

speckle suppression due to insufficiently aggressive PSF subtraction can also cause spurious features at small angular separations. Indeed, separate tests with more aggressive A-LOCI and KLIP settings for our data appear to confirm that this signal is due to residual speckles, but these approaches compromise our detections of the disk at 0.25–1" and thus are not used.

## 4.6. Comparison with Recent Studies

Based on SPHERE polarimetry, Engler et al. (2019) suggested the possibility of a misaligned/noncoplanar inner ring with a radius of ~1.3, ultimately finding a slightly better fit to their total intensity data for this geometry than for a one-ring or coplanar two-ring geometry. They further investigated the merit of the additional parameters of the two-ring model by comparing the Bayesian information criterion (BIC) for the two models, concluding that the BIC metric supports their best-fit noncoplanar two-ring model. While our modeling shows that the CHARIS data are consistent with a misaligned inner ring for some combinations of inclination and PA (see Sections 4.4 and 4.5), we find no clear evidence indicating the presence of an inner ring oriented as hypothesized from SPHERE imagery (PA = 276°, i = 80°; see bottom panel of Figure 12). Carrying out forward modeling on the best-fitting misaligned two-ring model identified in Engler et al. (2019) appears to reinforce this, with the model producing a $\chi^2_{\nu,\text{tot}}$ of 2.13 (for the same assumption of M = 10 free model parameters that they indicate) after forward modeling for our four reductions, compared with 1.151 for the overall best two-ring model we identify (see fourth row of Figure A1). Even allowing freedom for the other parameters, our two-ring optimization identified
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Lowest $\chi^2_{\nu,\text{tot}}$ two-ring model ($\chi^2_{\nu,\text{tot}} = 1.151$) identified by the procedure outlined in Section 4.1. The schematic overlaid in panel (a) shows the inner ring as a solid blue ellipse and the outer as a dashed blue ellipse. The $\chi^2$ map shown in panel (d) is displayed with identical scaling to that of Figure 9. Otherwise, same as Figure 9.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Reproduced from Figure 8: best $\chi^2$ for two-ring models as a function of inner disk radius and inclination. Values corresponding to the best model and acceptable models following the DE procedure of Section 4.1 (i.e., $\chi^2_\nu$ accounting for only the A-LOCI reductions) are shown in black and gray, respectively. The red line indicates values of radius and inclination that produce a projected semiminor axis equal to that of the adopted outer ring. Acceptable two-ring models are identified only for inner disk parameters resulting in nearly the same projected semiminor axis (see Section 4.5).</figcaption></figure>

no strong solutions with the inner ring oriented similarly (see Figure 8). The best coplanar two-ring model they identified fits our data somewhat better, resulting in a $\chi^2_{\nu,\text{tot}}$ of 1.7 (see bottom row of Figure A1). Models with a similar inner-ring radius and roughly coplanar orientation manifest in our final set of

acceptable models (Table 3) when using the Hong scattering phase function instead.

From ALMA observations, MacGregor et al. (2019) favored a disk model composed of two coplanar rings or a single ring with a Gaussian gap. Given the differences in parameterizations between our models, it is difficult to unambiguously translate their results for direct application to our data. However, the gap suggested by their models in either case is small enough (~14 au) that the profile should manifest consistently with the profile we observe (e.g., with the appearance of a single spine in our imagery; see Section 4.5 for relevant discussion regarding this constraint). We also note that our spine trace (Figure 6) shows an ~1σ shift around 0.8 on either side that is roughly coincident with the inner edge of the inner ring they proposed. By applying a simple 3 pixel rolling weighted average to the spine trace and disk-projected FWHM measurements, this feature becomes clearer (see Figure 13). While we find no significant evidence to support the presence of the ~14 au gap that they favored (given that it falls outside of our field of view, with $\rho \sim 1''2$), their interpretation appears generally consistent with CHARIS imagery.

The smoothed disk spine in Figure 13 also appears remarkably similar to the single profile spine trace of NZ Lup's disk (another highly inclined debris disk) reported in Boccaletti et al. (2019; their Figure 3). Boccaletti et al. (2019) ultimately favored a mutually inclined ($\Delta i \sim 5^\circ$) two-ring model for NZ Lup with an apparent gap that is roughly coincident with the dip in spine position seen in their spine trace. If a comparable explanation is assumed for the ~0.75–0.80 feature in our smoothed spine trace, the result would be a two-ring disk with a gap similar to the one suggested by MacGregor et al. (2019) but at a somewhat smaller separation than their best fit. However, given that the ALMA observations trace a significantly different dust
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Top: image of A-LOCI-processed August 30 CHARIS data with overlaid ellipses corresponding to the outer (dashed) and inner (solid) rings of our best-fit two-ring solution (see Section 4.2). While our solution's inner ring has an inclination offset of 5°9 and a PA offset of 2°7 with respect to the outer ring, the best-fit inner ring's radius results in an inner ring whose features predominantly coincide with those of the outer ring along our line of sight. Bottom: same as above but with the schematic depicting the inner and outer rings of the best-fit misaligned two-ring geometry from the results of Engler et al. (2019; see also the second row from the bottom of Figure A1). An inner ring oriented as posited by Engler et al. (2019) is not evident in the CHARIS data; to exist at such an orientation and still be consistent with our data, it would need to be substantially dimmer than the outer ring at similar projected separations.</figcaption></figure>

population than ours, these results may be fully consistent with one another. Notably, this interpretation manifests similarly to the overall best two-ring model identified in Section 4.4 (see Table 3), with a fiducial inner-ring radius of 40.9 au (0.83) and a mutual inclination of 5.9. The slight difference in the location of the peak FWHM measured for the disk between the east and west sides in Figure 13 (~0.75 and ~0.82, respectively) might be explained by the difference in PA between the rings suggested by our best two-ring model (ΔPA = 2°7).

# 5. Disk Surface Brightness and Color

## 5.1. Disk Photometry

To analyze the brightness and color of HD 15115's disk, we follow the general procedure of Goebel et al. (2018) to produce surface brightness profiles in the CHARIS broad band and *J*, *H*, and *K* bands. However, instead of fitting a fourth-order polynomial to the identified disk spine (see Section 3), we adopt the positions of the best-fit ellipse as the location of the disk spine for all imagery utilized.<fnref n="38" /> For measuring surface brightness, imagery is not rotated as it was in Section 3. Rather, the spine locations measured for the rotated images are

transformed to the native (north-up) image orientations, eliminating the possibility of image rotation interpolation affecting surface brightness measures. Flux attenuation cubes are created by dividing the PSF-subtracted best disk model by the pre-PSF-subtracted disk model (see Section 4.1). Bhowmik et al. (2019) showed that a comparable procedure results in erratic attenuation estimates for their KLIP reduction of SPHERE data for the highly inclined debris disk system HD 32297. However, this behavior does not manifest in our case (see Figure 14). These attenuation measurements are then used to correct the PSF-subtracted cubes produced with the A-LOCI and KLIP reductions of the August 30 and September 7 data. As in Goebel et al. (2018), we see fractional attenuation that tends to increase at smaller separations and further from the disk spine. Along the spine, disk flux in the CHARIS broad band is attenuated by ~45%–55% at 0.25 separation and ~30%–35% at 0.75. This attenuation varies by wavelength, with the *J* band typically being most attenuated (~65% and 45% at 0.25 and 0.75, respectively), followed by the *H* band (with values comparable to those in the broad band), and with the *K* band suffering the least attenuation (~45% and 30% at 0.25 and 0.75, respectively).

Following this, the channels of the attenuation-corrected reduction products are merged to create images corresponding to the *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. For each (x, y) pixel position along the spine, we take the nominal surface brightness to be the average flux value within a

<footnote n="38">³⁸ Though different observing wavelengths may trace distinct dust populations, resulting in different spine positions, testing showed that the utilized positions fall very near to locations we fit for HST/STIS imagery where meaningful fitting was feasible. For the purpose of surface brightness measurements, the spine identified from CHARIS broadband imagery appears to be a reasonable approximation of the spine for all bands we analyze.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Disk spine position and projected FWHM measurements made in Section 3 with a 3 pixel wide rolling weighted average applied. At ~0''75–0''80 on either side, a ~1σ–2σ dip in minor-axis separation is coincident with an increase in the measured disk FWHM. Further, a peak in FWHM can be seen on both sides but occurs slightly asymmetrically.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Estimates of fractional disk flux attenuation as a function of separation along the major axis. Solid and dashed lines show the values along the western and eastern extents of the disk, respectively. Values manifest similarly in all four reductions but are shown here for our A-LOCI reduction of August 30 data. These values correspond to $1 - \frac{1}{c}$, where c is the effective attenuation correction applied during measurement of surface brightness (see Section 5).</figcaption></figure>

circular aperture with diameter 0''12 (approximately the narrowest observed disk FWHM in CHARIS broadband imagery; this aperture size is used for all imagery). To accommodate the inclusion of partial pixels, we take the average in an aperture to be the weighted mean of the values with weights equal to the exact fraction of each pixel that is included in the aperture.

The uncertainty for each surface brightness measurement is determined as follows using non-attenuation-corrected images (attenuation maps become extremely noisy at the faint edges of the disk, where both attenuated and unattenuated models have values near zero). If the surface brightness, $F_s$, is measured at the spine position $(x_s, y_s)$ with corresponding stellocentric polar coordinates $(r_s, \theta_s)$, we make additional measurements the same way within the apertures at an array of positions $(r_i, \theta_i)$, with $\theta_i$ placed every 10°. Any of these measurements whose aperture overlaps with any spine aperture are removed. We then compute the standard deviation of this array of measurements as $\sigma_{0,s}$. Since $\sigma_{0,s}$ is representative of the uncertainty in the average surface brightness for the aperture at position $(x_s, y_s)$ before applying the attenuation correction, we additionally compute the nominal surface brightness at $(x_s, y_s)$ in the uncorrected image, $F_{0,s}$. The effective attenuation correction applied is then $C_s = F_s/F_{0,s}$. From this, our final uncertainty for $F_s$ is taken to be $\sigma_s = C_s \cdot \sigma_{0,s}$.

The procedure above is repeated to get arrays of surface brightness and corresponding uncertainties for each filter and reduction. For a given filter, the final nominal surface brightness at each position is taken to be the inverse variance–weighted average of the corresponding surface brightness

measurements for each of the four reductions utilized. The corresponding uncertainty is taken to be the standard error on the weighted average (see footnote, Section 3). Surface brightness is measured in the same manner for HST/STIS imagery of the system, except that no attenuation correction is necessary, and only one reduction is used. The results of this procedure are depicted in Figure 15. We point out here that, although the disk appears to be recovered only marginally in J-band imagery (Figure 3), making surface brightness measurements over a large (0''12) aperture and averaging measurements for multiple reductions results in J-band surface brightness measurements with a reasonable signal-to-noise ratio (S/N).

The surface brightness measurements for each bandpass are then combined with measurements of the stellar flux (for CHARIS data, from analysis of satellite speckles during spectro-photometric calibration of the cubes, and for HST/STIS data, as reported in Schneider et al. 2014) to compute the local surface brightness of the disk relative to the stellar flux (see Figure 16).

The disk color in isolation from the stellar color can then be analyzed by taking the difference of the relative magnitudes computed above (see Figure 18), while the east-west asymmetry can be assessed by comparing measurements of opposing sides in a particular bandpass (see Figure 17). These results are discussed in Section 5.2.

## 5.2. Surface Brightness Results

Disk color—Though the nearly unprecedented field of view probed by our CHARIS observations precludes numerical
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Measurements of surface brightness for HD 15115's disk as a function of stellocentric separation along the major axis in CHARIS and HST/STIS data. Here (and for other figures from Section 5) "broadband" refers to CHARIS's broad band. We note that the measurements here and in other Section 5 figures are not contiguously independent as a result of the choice to use a single aperture size for all photometric bands analyzed.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Measurements of surface reflectance for HD 15115's disk in CHARIS and HST/STIS data.</figcaption></figure>

comparisons of photometry with most prior studies, more quantitative comparisons can be made. Our red (STIS/50CCD —BB) color and neutral NIR color measured for the disk (Figure 18) appear generally consistent with prior literature that diagnosed the inner disk region in the optical and NIR. For example, combining the original discovery observations of Kalas et al. (2007) with new HST/NICMOS imagery, Debes et al. (2008) showed the disk's optical–NIR color becoming redder toward smaller separations.

Numerical simulations in the Boccaletti et al. (2003) model disk colors for infrared bandpasses as a function of the dust size distribution's minimum grain size ($a_{min}$) and porosity ($P$). With $P = 0$, for 1.6 µm ($\sim$H) versus 2.2 µm ($\sim$K), they show

1. a blue color for $a_{min} \lesssim 0.25$ µm,

2. a red color for $0.25 \mu m \lesssim a_{min} \lesssim 2 \mu m$, and

3. a neutral color for $a_{min} \gtrsim 2 \mu m$ (and briefly for $a_{min} \sim 0.25 \mu m$, as the color changes from blue to red).

Rodigas et al. (2012) found a predominantly gray ($K_s - L'$) color (2.1 and 3.8 µm, respectively) across the disk from 1.1"1 to 1.2"45. Comparing this result with grain-color models, they suggest a distribution comprised of grains from ~3 to 10 µm. The results of Boccaletti et al. (2003) show that a minimum grain size of ~3–10 µm should also produce a neutral color for ($H - K$), consistent with our measurements (Figure 18) and the suggestions of Rodigas et al. (2012).

However, by comparing measurements in the STIS/50CCD bandpass with our CHARIS broadband measurements, we find a definitively red color throughout the region of overlap
<figure><img src="image_19.png" /></figure>

**Figure 17.** Relative surface brightness between the eastern and western extents of the disk for STIS/50CCD and the CHARIS broad band. In each subplot, the horizontal colored line indicates the weighted average of the constituent flux measurements. The asymmetry in *J*, *H*, and *K* manifests similarly to the CHARIS broad band, albeit with larger uncertainties.

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

**Figure 18.** Disk color as a function of stellocentric separation along the disk major axis as identified from Section 3 (see Section 5), with $\Delta(H - K) = (H - K)_{disk} - (H - K)_{star}$, etc. Here “STIS” refers to measurements in STIS/50CCD, while “BB” refers to measurements in CHARIS’s broadband (wavelength-collapsed) imagery. The horizontal dotted line in each subplot indicates the weighted average of the constituent flux measurements for that side of the disk. Major axis separation <0 corresponds to the eastern extent of the disk. Other combinations of NIR filters (e.g., $(J - K)$) show a predominantly neutral color similar to our $(H - K)$.

($0''4-1''0$). While a wide range of minimum grain size values can produce a neutral color, a much smaller range results in a strong red color. Given that the redder filters analyzed in Boccaletti et al. (2003) predict no significantly red colors for minimum grain sizes larger than ~1 µm, the measurement of a very red (STIS/50CCD-CHARIS broadband) color suggests that a smaller minimum grain size, ~0.25–1.0 µm, is needed to simultaneously produce the red (STIS/50CCD-CHARIS broadband) and gray IR colors that we observe.

Rodigas et al. (2012) computed a blowout size, $a_{BO}$, of ~1–3 µm for HD 15115, indicating that grains of the minimum grain size that we estimate above would likely be expelled from the system. While a larger porosity would increase the estimated minimum grain size, with $a_{min} \propto (1-P)^{-1}$ (Boccaletti et al. 2003), it should also increase the blowout size by a comparable factor, with $a_{BO} \propto (1-P)^{-1}$ (Arnold et al. 2019). A minimum grain size below the theoretical blowout size can be explained in a number of ways.

Hughes et al. (2018) noted that this phenomenon is commonly observed and suggested that it is likely the result of a change in grain collision physics near the limit of small grains. Alternatively, this could manifest if some mechanism for continually replenishing smaller grains is present, such as planetesimal collisions (e.g., Hahn 2010).

Disk asymmetry—Numerous mechanisms were proposed previously to explain the observed flux asymmetry of HD 15115’s ringlike disk and extended halo. In the debris disk’s discovery paper, Kalas et al. (2007) suggested the possibility of a past encounter with nearby star HD 12545 perturbing planetesimal orbits to cause the asymmetry. However, MacGregor et al. (2019) pointed out that the spatial motion of the two objects makes this encounter unlikely. Debes et al. (2009) explored the possibility of disk sculpting through interaction with the interstellar medium (ISM) to explain the bluer western color observed at large separations as well as the observed east-west brightness asymmetry. Since HD 15115’s
<figure><img src="image_21.png" /><figcaption>**Figure 19.** Broadband (1.13–2.39 µm) contrast curves for planet detection reductions of CHARIS HD 15115 data of August 30 and September 7 outlined in Section 2.2. A 5σ contrast is given as a function of stellocentric angular separation (arcsec; lower x-axis) and projected separation (au; upper x-axis). Triangles along the right edge of the figure and corresponding colored lines give the planet detection limits for 10 (green), 25 (blue), and 100 (red) Myr planets based on hot-start, solar metallicity, hybrid cloud planet evolution models of Spiegel & Burrows (2012). The displayed ages are chosen to correspond to various determinations for the HD 15115 system (see Section 6).</figcaption></figure>

motion lies primarily in the direction of its apparently truncated eastern extent, pressure from clumps of ISM gas might redistribute dust from the eastern side to the western side (assuming the motion of the ISM gas itself is favorable for this scenario); with smaller grains being more susceptible to this mechanism, this should result in both a bluer and brighter western disk. The substantial asymmetries in the outer disk halo uniquely revealed by Schneider et al. (2014) STIS imaging provide further evidence that the outer disk and halo are being perturbed. However, using the equations and parameters provided in Debes et al. (2009) gives an approximate “deflection radius” (the stellocentric radius beyond which dust grains are likely to be significantly affected by the interaction) of 100–200 au. Given CHARIS’s ~5–50 au field of view, it seems unlikely that ISM interactions could be responsible for the asymmetry we observe. Moreover, this interaction should preferentially redistribute smaller grains to the western side, which is not supported by our disk color and surface brightness asymmetry measurements (see Figures 17 and 18); the fact that we measure a similar overall asymmetry in the STIS and CHARIS broadband data within the CHARIS field of view (Δ*m* ~ 0.6 and 0.4 mag, respectively) suggests that the phenomenon at work changes the overall dust density between the east and west, without significantly affecting the shape of the grain size distribution.

The results of Mazoyer et al. (2014) showed that while the eastern and western extents are significantly asymmetrical in brightness, the system features a symmetrical ring. This casts doubt on explanations of the brightness asymmetry that would necessitate an observable geometric asymmetry. More recently, Sai et al. (2015) reported an eccentricity for HD 15115’s disk of *e* = 0.06, which could contribute to the asymmetry we observe. By itself, this eccentricity does not appear capable of producing an asymmetry of the observed size, with limited

testing showing an induced east-west asymmetry of ≲10%. However, beyond asymmetry resulting directly from the eccentricity, Hahn (2010) noted that such a system may manifest with asymmetric dust distributions as a result of the difference in orbital velocity between apsides effectively enforcing differing ejection criteria. It is unclear if this mechanism would be capable of producing asymmetry of the observed magnitude.

A number of studies have proposed the possibility of asymmetry resulting from dynamical interactions with an embedded planet, both for HD 15115 (e.g., Sai et al. 2015) and for similar nearly edge-on systems (e.g., HD 111520; Draper et al. 2016). Sai et al. (2015) suggested the possibility of planetesimals being trapped in the Lagrange points of an embedded planetary companion. Such an embedded planet might also induce other disk structures; as mentioned in Section 4.5, complicated disk structures such as spiral arms could be present here but self-obscured by the system’s steep inclination. Such a geometry might result in asymmetries similar to those we observe.

Asymmetries might also be induced by major collisions within the disk (Hahn 2010). The possibility of a minimum grain size below the blowout size from our color analysis could be explained by this. In contrast to suggestions of planet “signposts” in similar systems, Thebault et al. (2012) simulated the interactions of debris disks and planets and concluded that for edge-on systems, only weak asymmetries will typically result from planet interactions.

## 6. Limits on Planets

Following the procedure described in Currie et al. (2018) for planet forward modeling, we computed 5σ contrast limits in the CHARIS broad band for the planet detection reductions outlined in Section 2.2. We then mapped these contrasts to
planet detection limits using the hot-start, solar metallicity, hybrid cloud, synthetic planet spectra provided by Spiegel & Burrows (2012).

Model planet spectra corresponding to an array of distinct determinations for the system's age are utilized. These age determinations include possible membership in the TW Hydrae association from Banyan Σ (98% likelihood; Gagné et al. 2018) with age $10 \pm 3$ Myr (Bell et al. 2015), possible membership in the β Pictoris moving group (Moór et al. 2006) with age $24 \pm 3$ Myr (Bell et al. 2015), and various other methods summarized in Rhee et al. (2007) that yield an age of ~100 Myr. Each planet spectrum was convolved with the filter profile for CHARIS's broadband filter and integrated to determine the photometric bandpass flux. The flux measured for HD 15115 was then converted to an absolute flux (to match the planet spectra) to determine the contrast at which each planet model would manifest. These values are indicated along the right axis of Figure 19. While planet contrast is intrinsically more favorable in the K band, the contrasts achieved are superior for CHARIS broadband imagery to the extent that the broad band offers the strongest constraints on the presence of planets.

The results of this procedure show that our August 30 data reduction is sensitive to $10 M_J$ companions at the lower and upper suggested ages to separations of ~7.5 and ~16 au, respectively. We note that, given the small mass of the possible companion proposed by MacGregor et al. (2019), $0.2 M_J$, we can place no constraints regarding its appearance anywhere within our field of view. On the other hand, for the scenario of a $12 M_J$ companion at a separation of 45 au discussed in Sai et al. (2015), we can rule out the planet over the majority of its orbit (e.g., ~93% of its orbit for an age of 25 Myr).

## 7. Conclusions and Future Work

CHARIS imagery of the HD 15115 system has revealed the inner regions of the disk in remarkable detail and probed substantially further than any previous scattered-light data (to $\rho \sim 0''2$). This imagery revealed no direct evidence of planetary-mass companions and has allowed for new constraints to be placed on the possibility of a yet-unseen substellar companion in the disk. Combined with the DE algorithm, CHARIS imagery has enabled us to conduct a thorough exploration of the recently proposed inner ring through forward modeling. In doing so, we find a poor fit for a significantly noncoplanar inner ring but reasonable fits for both a single ring and two rings aligned along our line of sight (either coplanar or manifesting with similar projected semiminor axes). These data, combined with HST/STIS imagery, have allowed for measurement of the disk's color and asymmetry at separations from 0''25 to 1''0 and spanning wavelengths from 0.6 to 2.3 µm. These measurements suggest a minimum grain size in

the CHARIS field of view of $\lesssim 1.0$ µm and thus smaller than previous estimates at larger separations.

The CHARIS observations presented here provide the first clear view of the system within $\rho \sim 0''4$. In general, follow-up observations probing this region of the system will better substantiate the results of our analysis. Follow-up observations with CHARIS would enable further constraints to be placed on the presence of inner disk features or companions, as well as gauging the significance of the $\rho \lesssim 0''25$ feature we note in Section 4.5. The use of CHARIS's new polarimetric integral field spectroscopy mode would allow for measurement of the disk's fractional polarization in CHARIS's field of view, a key diagnostic of the disk's dust properties (Perrin et al. 2015), while also allowing more rigorous assessment of any planet candidates that might be identified. High-S/N mid-IR spectra of HD 15115 (e.g., with the James Webb Space Telescope/MIR) could better constrain the dust composition within the disk by identifying the signatures of both silicates and nonsilicate species using spectral decomposition software (Hughes et al. 2018).

We thank our referee for providing helpful comments that improved this manuscript.

The authors wish to acknowledge the very 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. We wish to acknowledge the critical importance of the current and recent Subaru telescope operators, day crew, computer support, and office staff employees. Their expertise, ingenuity, and dedication are indispensable to the continued successful operation of Subaru.

We acknowledge funding support from the NASA XRP program via grants 80NSSC20K0252 and NNX17AF88G. T.C. was supported by a NASA Senior Postdoctoral Fellowship.

Based partially on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. These observations are associated with program No. 12228.

E.A. is supported by MEXT/JSPS KAKENHI grant No. 17K05399.

## Appendix A
### Scattering Phase Function Comparison

Figure A1 shows the results of forward modeling for models of various phase functions for our August 30 data. Overall, a simple HG phase function seems to very poorly describe the brightness profile that we observe.
<figure><img src="image_22.png" /></figure>

**Figure A1.** Each row's images depict (left to right) an initial two-ring disk model, the model following attenuation by the forward-modeling procedure for our August 30 A-LOCI reduction and rescaling to minimize $\chi^2$ (see Section 4.1), the residual after subtracting the model from the data, and the corresponding $\chi^2$ map. The value of $\chi^2$ given in the last panel is for the reduction shown only. In each case, images are shown at the same linear display stretch as the corresponding images in Figure 9. The first three rows show models in which only the phase function changes; the parameters are otherwise identical and correspond to the best overall two-ring disk model identified in Section 4.1. The model in the first row utilizes the same phase function implemented by Engler et al. (2019) for their two-ring models: a simple HG phase function with asymmetry parameters $g = 0.4$ for the inner ring and $g = 0.3$ for the outer. The second row's model changes the asymmetry parameter of both rings to $g = 0.7$, matching the highest-weighted term in the Hong phase function. The third row's model utilizes the phase function of Hong (1985), as adopted for our modeling procedure (see Section 4.2). The models in the fourth and fifth rows adopt the noncoplanar and coplanar (respectively) best-fitting two-ring models reported in Section 5.2 of Engler et al. (2019), which feature the same phase function as the model of the first row but with differing parameters elsewhere. Note that many of the models utilizing a simple HG phase function appear especially dim in the “processed model” panels as a result of the rescaling applied at the end of the forward-modeling procedure; this is simply the scaling of the model that best minimizes the weighted residuals.

## Appendix B
Python Implementation of DE

Here we provide a simple Python implementation of the DE algorithm (Storn & Price 1997) as described in Section 4.3. This code favors readability and simplicity over perfect computational efficiency and has no dependencies besides the NumPy<fnref n="39" /> module. In comparison with grid searches, this implementation of DE will result in both a superior fit and orders of magnitude fewer model evaluations for the overwhelming majority of cases.

```python
import numpy as np
def differential_evolution(objective_fn, converged, bounds, mutation = (0.5,1.), P = 0.7, popsize = 10):
    ```

A simple implementation of the differential evolution algorithm using the 'bestlbin'' strategy and allowing a 'dithered'' mutation constant.

### Parameters

objective_fn: callable

A function that takes the model parameters (1d array of length K) as its argument and returns the value to be minimized (typically chi-squared). This function should: generate the appropriate model from the list of parameters, propagate the model through your forward modeling routine, and then compare the model to your data to determine its fitness. You will probably want to have this function save the input and output models to disk as well.

<footnote n="39">³⁹ https://numpy.org/</footnote>
converged : callable

A function that takes the current (normalized)
population (2d array of shape (popsize\*K,
K) for K parameters) and their fit metrics
(1d array of length (popsize\*K)), returning
True if some convergence criteria has been
met and fitting should cease and False
otherwise.

bounds : numpy array of shape (K,2) where K is
the number of model parameters

Each entry of `bounds'`, `bounds[i.:]`, should
provide the lower and upper bound for a
parameter.

mutation : float or tuple (float, float), optional

The mutation constant to utilize. Storn &
Price (1997) suggested that values in the
range [0.4, 1.0] are typically more favor-
able. If given as a tuple, the mutation
constant is randomly selected each genera-
tion from the uniform distribution spanning
the two values given.

P : float, optional

The crossover probability to utilize. The
value of P should be in the range (0,1].

popsize : int, optional

The number of population members per free
parameter to utilize.

Returns

—: array, float

The set of best-fitting parameters and the
associated fitness metric.

"

N, K = bounds.shape[0]*popsize, bounds.shape
[0] # number pop. members and parameters

bmin, brange = bounds[:,0], np.diff(bounds.
T, axis = 0) # lower lims and range for each
param

x = np.randomrand(N, K) # Generate initial
(normed) population array

fx = np.array([objective_fn(xi) for xi in x
*brange+bmin]) # The initial pop's fitness

indices = np浓arange(N) # Define indices corresp-
onding to population members

while not converged(x,fx): # Loop until con-
verged(x,fx) returns True

if type(mutation) == tuple: m = np.random.
uniform(*mutation) # For dithered m

else : m = mutation

xtrial = np zeros_like(x)

j = np.argmin(fx) # For best1bin method, j is
the index of the best member

for i in indices:

k,l = np.random.choice indixes np.isin
(indixes, [i,j])], 2, replace = False)

xmi = np .clip(x[j] + m* (x[k]-x[l]) , 0, 1) # ith
mutant vector, clipped to bounds

xtrial[i] = np .where(np .random .rand(K) < P,
xmi, x[i]) # Get trial pop. vector

fxtrial = np .array([objective_fn(xi) for xi
in xtrial*brange+bmin]) # Fitness of trial
pop.

improved = fxtrial < fx # Boolean array indi
cating which trial members were improvements
x [improved], fx [improved] = xtrial [improved],
fxtrial [improved] # Replace improved members
return x np .argmin(fx)]*brange+bmin, np .min
(fx) # Return the best params and fitness

The code as presented can be easily adapted for parallel
processing with minor changes to the two lines that evaluate
the fitness for a set of model parameters; e.g., using the Joblib
module,<fnref n="40" /> the fourth line of code in the function could be
replaced with (likewise for the 15th line):

from joblib import Parallel, delayed
fx = np.array(Parallel () (delayed(objecti
ve_fn) (xi) for xi in x*brange+bmin)).

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

Glenn Schneider https://orcid.org/0000-0002-4511-5966

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

N. Jeremy Kasdin https://orcid.org/0000-0002-6963-7486

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

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

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

Eiji Akiyama https://orcid.org/0000-0002-5082-8880

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

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

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

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

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

Ruben Asensio-Torres https://orcid.org/0000-0003-2990-0726

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

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

References

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

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

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

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

Bhowmik, T., Boccaletti, A., Thébault, P., et al. 2019, A&A, 630, A85

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

Boccaletti, A., Thébault, P., Pawellek, N., et al. 2019, A&A, 625, A21

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

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. 2012, ApJL, 760, L32

Currie, T., Guyon, O., Lozi, J., et al. 2019a, Proc. SPIE, 11117, 111170X

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

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

Debes, J. H., Weinberger, A. J., & Kuchner, M. J. 2009, ApJ, 702, 318

<footnote n="40">⁴⁰ https://joblib.readthedocs.io</footnote>
Debes, J. H., Weinberger, A. J., & Song, I. 2008, ApJ, 684, 1

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

Duchene, 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

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

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

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

Hahn, J. M. 2010, ApJ, 719, 1699

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

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

Hughes, A. M., Duchene, G., & Matthews, B. 2018, ARA&A, 56, 541

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

Kalas, P., Fitzgerald, M. P., & Graham, J. R. 2007, ApJL, 661, L85

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

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

Labroquère, J., Héritier, A., Riccardi, A., & Izzo, D. 2014, in Parallel Problem Solving from Nature-PPSN XIII, ed. T. Bartz-Beielstein et al. (Cham: Springer), 262

Lawson, K. D., Wisniewski, J. P., Bellm, E. C., Kowalski, A. F., & Shupe, D. L. 2019, AJ, 158, 119

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

MacGregor, M. A., Weinberger, A. J., Nesvold, E. R., et al. 2019, ApJL, 877, L32

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

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

Mazoyer, J., Boccaletti, A., Augereau, J. C., et al. 2014, A&A, 569, A29

Moór, A., Ábrahám, P., Derekas, A., et al. 2006, ApJ, 644, 525

Moór, A., Pascucci, I., Kóspál, Á, et al. 2011, ApJS, 193, 4

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

Pueyo, L. 2016, ApJ, 824, 117

Rhee, J. H., Song, I., Zuckerman, B., & McElwain, M. 2007, ApJ, 660, 1556

Rodigas, T. J., Hinz, P. M., Leisenring, J., et al. 2012, ApJ, 752, 57

Sai, S., Itoh, Y., Fukagawa, M., Shibai, H., & Sumi, T. 2015, PASJ, 67, 20

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

Shinn, J.-H. 2018, ApJS, 239, 21

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

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

Storn, R., & Price, K. 1997, Journal of Global Optimization, 11, 341

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

Thebault, P., Kral, Q., & Ertel, S. 2012, A&A, 547, A92

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

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