# JWST/NIRCam Coronagraphy of the Young Planet-hosting Debris Disk AU Microscopii

Kellen Lawson, Joshua Schlieder, Jarron Leisenring, Ell Bogat, Charles Beichman, Geoffrey Bryden, András Gáspár, Tyler Groff, Michael Mcelwain, Michael Meyer, et al.

► To cite this version:

Kellen Lawson, Joshua Schlieder, Jarron Leisenring, Ell Bogat, Charles Beichman, et al. JWST/NIRCam Coronagraphy of the Young Planet-hosting Debris Disk AU Microscopii. The Astronomical Journal, 2023, 166 (4), pp.150. 10.3847/1538-3881/aced08. hal-04268671

HAL Id: hal-04268671

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

Submitted on 3 Nov 2023

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 CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# JWST/NIRCam Coronagraphy of the Young Planet-hosting Debris Disk AU Microscopii

Kellen Lawson<sup>1,15</sup>, Joshua E. Schlieder<sup>1</sup>, Jarron M. Leisenring<sup>2</sup>, Ell Bogal<sup>1,3</sup>, Charles A. Beichman<sup>4</sup>, Geoffrey Bryden<sup>5</sup>, András Gáspár<sup>2</sup>, Tyler D. Groff<sup>1</sup>, Michael W. McElwain<sup>1</sup>, Michael R. Meyer<sup>6</sup>, Thomas Barclay<sup>1,7</sup>, Per Calissendorff<sup>6</sup>, Matthew De Furio<sup>6</sup>, Marie Ygouf<sup>5</sup>, Anthony Boccaletti<sup>8</sup>, Thomas P. Greene<sup>9</sup>, John Krist<sup>10</sup>, Peter Plavchan<sup>11</sup>, Marcia J. Rieke<sup>2</sup>, Thomas L. Roellig<sup>9</sup>, John Stansberry<sup>12</sup>, John P. Wisniewski<sup>13</sup>, and Erick T. Young<sup>14</sup>

<sup>1</sup> NASA-Goddard Space Flight Center, Greenbelt, MD, USA; kellenlawson@gmail.com

<sup>2</sup> Steward Observatory, University of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA

<sup>3</sup> Department of Astronomy, University of Maryland, College Park, MD 20782, USA

<sup>4</sup> NASA Exoplanet Science Institute/IPAC, Jet Propulsion Laboratory, California Institute of Technology, 1200 E California Boulevard, Pasadena, CA 91125, USA

<sup>5</sup> Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA

<sup>6</sup> Department of Astronomy, University of Michigan, 1085 S. University, Ann Arbor, MI 48109, USA

<sup>7</sup> University of Maryland, Baltimore County, Baltimore, MD 21250, USA

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

<sup>9</sup> NASA Ames Research Center, MS 245-6, Moffett Field, CA 94035, USA

<sup>10</sup> Jet Propulsion Laboratory, 4800 Oak Grove Drive M/S 183-900, Pasadena, CA 91109, USA

<sup>11</sup> George Mason University Department of Physics & Astronomy, 4400 University Drive, MS 3F3, Fairfax, VA 22030, USA

<sup>12</sup> Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

<sup>13</sup> NASA Headquarters, 300 Hidden Figures Way SW, Washington, DC 20546, USA

<sup>14</sup> Universities Space Research Association, 425 3rd Street SW, Suite 950, Washington DC 20024, USA

Received 2023 May 11; revised 2023 July 14; accepted 2023 August 1; published 2023 September 7

## Abstract

High-contrast imaging of debris disk systems permits us to assess the composition and size distribution of circumstellar dust, to probe recent dynamical histories, and to directly detect and characterize embedded exoplanets. Observations of these systems in the infrared beyond 2–3 µm promise access to both extremely favorable planet contrasts and numerous scattered-light spectral features—but have typically been inhibited by the brightness of the sky at these wavelengths. We present coronagraphy of the AU Microscopii (AU Mic) system using JWST’s Near Infrared Camera (NIRCam) in two filters spanning 3–5 µm. These data provide the first images of the system’s famous debris disk at these wavelengths and permit additional constraints on its properties and morphology. Conducting a deep search for companions in these data, we do not identify any compelling candidates. However, with sensitivity sufficient to recover planets as small as ~0.1 Jupiter masses beyond ~2" (~20 au) with 5σ confidence, these data place significant constraints on any massive companions that might still remain at large separations and provide additional context for the compact, multiplanet system orbiting very close-in. The observations presented here highlight NIRCam’s unique capabilities for probing similar disks in this largely unexplored wavelength range, and they provide the deepest direct imaging constraints on wide-orbit giant planets in this very well-studied benchmark system.

*Unified Astronomy Thesaurus concepts: Exoplanet systems (484); Direct imaging (387); Debris disks (363); Coronagraphic imaging (313); M-dwarf stars (982)*

## 1. Introduction

Compact multiplanet systems are ubiquitous around the galaxy’s lowest-mass stars (Dressing & Charbonneau 2015; Hardegree-Ullman et al. 2019), so their formation and evolution, as well as the impact of giant, outer companions, are of great interest. This is particularly true in the context of galactic planet demographics (Barclay et al. 2017; Meyer et al. 2018; Gaudi et al. 2021), exoplanet habitability (Clement et al. 2022), and as analogs to Solar System formation (Childs et al. 2019). The presence or absence of giant planets in the outskirts of M-dwarf systems can have a substantial impact on the delivery of the volatiles necessary for habitability to terrestrial planets deep in their interiors (Alibert & Benz 2017; Bitsch

et al. 2019; Schlecker et al. 2021; Clement et al. 2022). To date, such planets remain generally elusive and particularly difficult to discover (Bowler et al. 2015; Lannier et al. 2016). Though microlensing surveys provide evidence for a significant population of sub-Jupiters at ~1–10 au separations (Shvartzvald et al. 2016; Suzuki et al. 2016), direct imaging studies have previously been limited in their ability to detect planets in this mass regime. Early work by Beichman et al. (2010) showed that the wavelength coverage and predicted contrast of JWST’s high-contrast imaging modes would be ideal to detect such planets around nearby M dwarfs. These capabilities would open an unprecedented region of the mass-separation parameter space for exoplanets and provide unique insights into the population of planets in the outskirts of M-dwarf systems. These predictions led to a deeper investigation into the optimal targets for coronagraphy from JWST’s Near-infrared Camera (NIRCam; Rieke et al. 2005, 2023) and new analyses of JWST’s high-contrast imaging capabilities in the context of M dwarfs using updated performance predictions and planet models (Schlieder et al. 2015; Brande et al. 2020; Carter et al.

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

<footnote n="15"><sup>15</sup> NASA Postdoctoral Program Fellow.</footnote>
2023). These studies revealed that particularly nearby and young M dwarfs provide access to the lowest-mass planets at the closest separations—with the potential to push below a Saturn mass at $\lesssim 10$ au.

AU Microscopii (AU Mic, HD 197481, GJ 803) is a very nearby ($d = 9.714 \pm 0.002$ pc; Gaia Collaboration et al. 2023), bright ($V = 8.6$, $K_s = 4.5$) M dwarf. Spectrophotometric observations from the late 1950s provided the first estimates of an early-M type (Vyssotsky 1956). Keenan & McNeil (1989) measured a spectral type of M1V using modern spectroscopic techniques. Indications of strong magnetic activity and a pre-main-sequence age were emerging as early as the 1970s, with the detection of chromospheric emission lines, flares, and spot-driven variability (Wilson & Woolley 1970; Kunkel 1973; Torres & Ferraz Mello 1973).

Barrado y Navascués et al. (1999) ultimately confirmed the young age of AU Mic by a combination of kinematic association with the young A-type star $\beta$ Pictoris ($\beta$ Pic), elevated HR diagram position, and strong activity. These diagnostics indicated an age of $20 \pm 10$ Myr. Further follow-up revealed that AU Mic and many other nearby, young M dwarfs, are members of the $\beta$ Pic moving group (βPMG) kinematic association (e.g., Zuckerman et al. 2001; Schlieder et al. 2010, 2012; Shkolnik et al. 2017; Gagné et al. 2018). The age of the βPMG has been determined to be $\approx 20-25$ Myr using numerous techniques that include color–magnitude diagram and isochrone analyses (e.g., Bell et al. 2015), lithium depletion boundary measurements (e.g., Binks & Jeffries 2014), binary dynamical mass constraints (e.g., Nielsen et al. 2016), and kinematic traceback estimates (e.g., Miret-Roig et al. 2020; Couture et al. 2023; see Table 6 of Miret-Roig et al. 2020 and references therein for a literature summary of βPMG age estimates). Bell et al. (2015) provide an isochronal age estimate of $24 \pm 3$ Myr for the βPMG, which we adopt as the age of AU Mic hereafter.

The Infrared Astronomical Satellite (IRAS) provided the first evidence of a debris disk around AU Mic via a $60 \ \mu\text{m}$ excess (Mullan et al. 1989; Mathioudakis & Doyle 1991, 1993). Subsequent reanalysis of IRAS data by Song et al. (2002) confirmed the infrared (IR) excess at $60 \ \mu\text{m}$ and refined the model of the spectral energy distribution (SED). At $850 \ \mu\text{m}$, the system was detected as a point source—indicating the presence of an unresolved population of cold ($\sim 40$ K) dust (Liu et al. 2004).

Around this time, visible-band coronagraphy provided the first resolved images of the nearly edge-on disk—extending from the inner edge of the coronagraph at $\sim 50$ au to separations of $\sim 210$ au in reflected light (Kalas et al. 2004). Follow-up observations using near-infrared (NIR) coronagraphy revealed the disk extending to smaller separations and identified substructure and asymmetries between the southeast and northwest sides (Liu 2004). Over the nearly two decades that followed, continuous high-contrast and high-resolution imaging spanning optical to millimeter wavelengths (e.g., Fitzgerald et al. 2007; MacGregor et al. 2013) have made the AU Mic disk one of the best studied. Notable features include a planetesimal belt at $\sim 40$ au (Augereau & Beust 2006), an extended halo (Matthews et al. 2015), and dynamical changes in the disk, including a color change (Lomax et al. 2018) and fast-moving substructures (Boccaletti et al. 2015; Sezestre et al. 2017; Boccaletti et al. 2018; Wisniewski et al. 2019; Grady et al. 2020).

More recently, a compact system of transiting exoplanets—well within the inner edge of its planetesimal belt—has been identified (Plavchan et al. 2020; Martioli et al. 2021; Gilbert et al. 2022; Wittrock et al. 2023). Two of the planets, AU Mic b and c, have radii and masses comparable to those of Neptune and Uranus (Cale et al. 2021; Klein et al. 2021; Zicher et al. 2022), but have orbital periods of only 8.46 and 18.86 days (semimajor axes of 0.0645 and 0.1101 au, respectively). These periods lead to a resonance that drives transit timing variations (TTVs) on the order of minutes (Wittrock et al. 2022), providing a deeper characterization of the system and evidence for additional perturbers. Wittrock et al. (2023) identify a third planet, AU Mic d, via TTVs, with a model-favored period of 12.74 days and a mass comparable to that of the Earth. The periods of the AU Mic planets are close to a 4:6:9 resonance and provide clues to planet formation and evolution in the deep interior of this young system (Wittrock et al. 2023).

A selection of M-dwarf targets is currently being studied through a JWST NIRCam Guaranteed Time Observing (GTO) program. This program, GTO 1184 (PI J. Schlieder)<fnref n="16" /> targets nine of the closest, youngest M dwarfs using NIRCam coronagraphy in two filters spanning 3–5 $\ \mu\text{m}$. AU Mic’s fourfold combination of proximity, young age, dynamically active debris disk, and exoplanet system place it among the richest laboratories to study the formation and evolution of planetary systems around low-mass stars. Moreover, analysis in Sezestre et al. (2017) explored the origin of AU Mic’s fast-moving disk features and found plausible scenarios involving an orbiting companion with semimajor axis between roughly 5 and 25 au. For these reasons, AU Mic was considered a prime target for GTO 1184.

In this study, we present the results of JWST/NIRCam observations of AU Mic from GTO 1184. We report conclusive first detections of the famous debris disk at both 3.6 and 4.4 $\ \mu\text{m}$, and spanning separations from $\sim 0''3$ to $\sim 5''$ (2.9 to 49 au). With these data, we conduct a deep search for companions—effectively probing companion masses well below that of Saturn—and provide preliminary analyses of the disk itself. A follow-up study will perform further analysis focusing on the disk.

## 2. Observations

Rather than conducting dedicated reference observations to accompany each science target, GTO 1184 improves survey efficiency by adopting a self-referencing strategy for reference star differential imaging (RDI). In this strategy, RDI is performed for each target using the suite of other targets’ images as a reference library (with consideration for any off-axis sources that may be present in the library). All observations used the *MASK335R* coronagraph mask (Krist et al. 2009), which has an inner working angle (IWA)<fnref n="17" /> of $0''63$.

Observations were conducted using two filters from the NIRCam long wavelength (LW) channel (average pixel scale of 63 mas pixel$^{-1}$): F356W ($\lambda_{\text{pivot}} = 3.563 \ \mu\text{m}$, $\Delta\lambda = 0.787 \ \mu\text{m}$) and F444W ($\lambda_{\text{pivot}} = 4.421 \ \mu\text{m}$, $\Delta\lambda = 1.024 \ \mu\text{m}$) (Rieke et al. 2005). Using the *SUB320* subarray mode, each integration in the LW channel produces an

<footnote n="16">¹⁶ GTO 1184—A NIRCam Coronagraphic Imaging Survey of Nearby Young M Dwarfs</footnote>

<footnote n="17">¹⁷ Where the IWA is defined as the angular separation at which coronagraph transmission reaches 50%</footnote>
<table><caption>**Table 1**
JWST/NIRCam Observations</caption><thead><tr><th rowspan="2">Prop. ID<sup>a</sup></th><th rowspan="2">Spec. Type</th><th rowspan="2">W1 (mag)<sup>b</sup></th><th rowspan="2">W2 (mag)<sup>b</sup></th><th rowspan="2">Obs. Date</th><th colspan="2">F356W</th><th colspan="2">F444W</th></tr><tr><th>N<sub>int</sub><sup>c</sup></th><th>Exp. Time<sup>d</sup></th><th>N<sub>int</sub><sup>c</sup></th><th>Exp. Time<sup>d</sup></th></tr></thead><tbody><tr><td>V-AU-MIC</td><td>M1</td><td>4.45</td><td>4.01</td><td>2022 Oct 3</td><td>34</td><td>1708</td><td>70</td><td>3517</td></tr><tr><td>HIP 17695</td><td>M3</td><td>6.81</td><td>6.66</td><td>2022 Oct 3</td><td>16</td><td>1676</td><td>34</td><td>3562</td></tr><tr><td>G-7-34</td><td>M4</td><td>8.01</td><td>7.81</td><td>2022 Oct 3</td><td>16</td><td>1676</td><td>34</td><td>3562</td></tr><tr><td>TYC 5899<sup>e</sup></td><td>M3</td><td>6.78</td><td>6.58</td><td>2022 Oct 3</td><td>8</td><td>838</td><td>17</td><td>1781</td></tr><tr><td>2MJ0443<sup>f</sup></td><td>M9</td><td>10.83</td><td>10.48</td><td>2022 Sep 6</td><td>16</td><td>1676</td><td>34</td><td>3562</td></tr><tr><td>LP-944-20</td><td>M9</td><td>9.13</td><td>8.80</td><td>2022 Sep 6</td><td>16</td><td>1676</td><td>34</td><td>3562</td></tr></tbody></table>

**Notes.** A summary of the science and reference targets and their observations used in this study. All targets except AU Mic ("V-AU-MIC") were used as reference targets.

<sup>a</sup> The identifier used for each target in the GTO 1184 proposal.

<sup>b</sup> ALLWISE photometry; Wright et al. (2010).

<sup>c</sup> The total number of integrations for the corresponding filter; every integration had 10 groups.

<sup>d</sup> The total effective exposure time over all integrations in units of seconds.

<sup>e</sup> TYC 5899-193-1.

<sup>f</sup> 2MASS J0443376+000205.

image of 320 × 320 pixels (20" × 20"). The F356W and F444W filters were selected primarily to maximize sensitivity to planets and to enable rejection of background sources via color analyses. For NIRCam's non-coronagraphic imaging mode, wavelengths shorter than ~4 µm in the LW channel produce point-spread functions (PSFs) with full width at half maximum (FWHM) smaller than two pixels and are thus considered undersampled. However, for coronagraphy, the Lyot stop that is paired with the round coronagraph masks reduces the effective telescope diameter from 6.5 to 5.2 m. This results in PSF FWHMs of approximately 2.29 pixels (0"/14) and 2.84 pixels (0"/18) for F356W and F444W, respectively. As such, neither filter's PSF is undersampled when used for coronagraphy.

All targets were observed at two roll angles separated by ~10°, with the exception of TYC 5899 (for which only a single roll was executed on the indicated date, due to a target acquisition failure). To accommodate its relative brightness, observations of AU Mic used the *SHALLOW2* readout pattern, while all other observations used *MEDIUM8*. The signal-to-noise ratio in the AU Mic integrations is somewhat higher than for the reference integrations (by a factor of ~2 compared to HIP 17695 and TYC 5899—the closest targets in brightness). However, when the reference integrations are combined during RDI (see Section 3.3), the S/N for the resulting models of the stellar diffraction pattern in each science integration is much more comparable. Nevertheless, this should be expected to slightly diminish the sensitivity in the background-limited regime ($r \gtrsim 2''$) compared to a conventional RDI sequence in which the reference observations are tuned to achieve comparable S/N.

A summary of the GTO 1184 JWST/NIRCam observations used in this study, including observation dates and instrument settings, is provided in Table 1. Targets utilized as references were found to be free of extended circumstellar emission, with any nearby sources in the reference images being accounted for as described in Section 3.3. The spectral type mismatch between AU Mic and the reference targets is not expected to significantly impact sensitivity (e.g., Girard et al. 2022). Other available GTO 1184 observations excluded from Table 1 were not ultimately used, as they were not found to improve RDI subtraction for AU Mic. A full description of the GTO 1184

program and the survey results will be presented in a forthcoming publication.

### 3. Data Reduction and Post-processing

To reduce the data, we make use of the *spaceKLIP* package (Kammerer et al. 2022), which in turn uses the *JWST Pipeline* (Bushouse et al. 2022) and largely follow the procedure of Carter et al. (2023), summarized hereafter, with some exceptions as indicated.

Beginning from Stage 0 products (*uncal.fits*), we process the data to Stage 1 (*rateints.fits*) using the *rampfit* step of *spaceKLIP*. Following Carter et al. (2023), we (a) use the updated reference pixel definition described therein, (b) skip the dark current subtraction step in order to avoid the negative effects of the low-quality calibration data that are currently available, and (c) adopt a jump detection threshold of 5.

We then process the resulting Stage 1 products to Stage 2 (*calints.fits*) using the *imgprocess* step from *spaceKLIP*. However, we found that the included pixel cleaning procedures from both *spaceKLIP* and the *JWST Pipeline* left numerous problematic pixels uncorrected. The apparent differences in performance between application to the data of Carter et al. (2023) and our data may result from differences in the observing strategies—e.g., with Carter et al. (2023) observing a dedicated reference target with a small grid dither strategy versus our self-referenced survey lacking dedicated reference targets and dithers. Instead, we executed the *imgprocess* step without pixel cleaning, and performed subsequent pixel cleaning using the procedures outlined hereafter.

#### 3.1. *Outlier Rejection*

We rejected outlier pixels based on three criteria, initially replacing offending pixels with NaN values in each case:

1. pixels falling arbitrarily close to zero (absolute values less than 1e-8);

2. pixels having data quality (DQ) flag values above 100 (e.g., flagged cosmic ray events); and

3. pixels falling more than 7.5 median absolute deviations (MADs) from the median among the values of a given
pixel across all integrations of a single roll angle (e.g., unflagged cosmic ray events).

For each criterion, the numerical threshold was set ad hoc based on visual inspection of the data. After this procedure, we manually masked an additional 28 seemingly discrepant pixels that were not identified by the criteria above. Generally, these were pixels that were persistently much brighter or fainter than neighboring pixels (i.e., spatial outliers rather than temporal outliers). We remark that it may be feasible to automatically reject these pixels with careful application of something akin to a tophat filter. However, in limited testing, we were unable to identify a process that would not occasionally mask bright diffraction speckles as well. Once this process was finished, we replaced each masked pixel with the median of values within a 5 × 5 pixel window.

## 3.2. Image Registration

To identify the position of the occulted star in the data (and thus the misalignment between the star and the coronagraph), we largely follow the procedure of Carter et al. (2023). Because any pointing changes between integrations are expected to be negligible (~1 mas rms; Rigby et al. 2023), we median combine the integrations of each exposure for determining offsets before then applying corresponding shifts to the individual integrations.

To begin, we generate a synthetic coronagraphic image of a perfectly centered star using *WebbPSF* with the closest optical path difference (OPD) map to the observations (O2022100401-NRCA3_FP1-1.fits; measured approximately one day after the collection of the AU Mic data). We then use cross-correlation with this synthetic image to determine the position of the star in the first exposure of the observations of HIP 17695 (the target observed nearest in time to AU Mic)—hereafter referred to as the centering reference—using the *JWST Pipeline image registration procedure* (*imageregistration align_array*). We chose to use HIP 17695 for this purpose to avoid any impact from AU Mic's bright disk, which is faintly visible in the raw data. We then cross-correlated each additional exposure with the centering reference to determine the relative offsets between exposures.

Throughout this procedure, we consider only the pixels within a coronagraph-centered annulus spanning 12 < *r* < 32 pixels. While the methodology of Carter et al. (2023) uses an 11 × 11 pixel coronagraph-centered box for relative centering, our data span a much larger temporal baseline—making small-separation changes to the diffraction pattern more likely. Within the aforementioned annular region considered for relative centering, centering for the AU Mic images explicitly excludes a rectangular region approximately aligned with the disk's major axis (adjusted for the parallactic angle of the exposure)—assuming a disk position angle of 128°48 (Vizgan et al. 2022) and a width of 12 pixels. Once these offsets are computed, we shift all of the integrations such that the stars are aligned with the reference pixel for the NIRCam long wavelength (LW) target acquisition (TA) filter, F335M. This accounts for a recently identified filter-dependent offset for NIRCam, which shifts the entire scene relative to the TA filter (J. Leisenring, private communication). In detector pixels, these offsets are [*x*, *y*] = [0.751, −0.121] for F356W, and [*x*, *y*] = [0.157, −0.224] for F444W.

This procedure is carried out separately for each filter. For the two rolls of the F356W observations of AU Mic, we find [*x*, *y*] offsets from the coronagraph center to the star of [30 mas, −6 mas] and [18 mas, −6 mas]. For F444W, we find offsets of [24 mas, −6 mas] and [12 mas, −5 mas].

## 3.3. Reference Star Differential Imaging

In a pure “classical” RDI procedure, the pattern of diffracted starlight in the science target’s sequence would typically be modeled by taking the median of the sequence of reference images and then scaling the result to match the brightness of the target based on prior knowledge of the two stars’ fluxes. With consideration for the variability of the targets in our data and the lack of precise knowledge regarding the stellar fluxes in these bandpasses, such a procedure is problematic. Instead, we estimate the ratio of the stellar fluxes empirically by median combining the science target and reference sequences and then determining the scaling factor that minimizes the squared residuals between the two median images within a region that explicitly excludes the vicinity of AU Mic’s disk (using the same region defined in Section 3.2). We then multiply the median of the reference images by this scaling factor and subtract it from the target sequence. For this purpose, the reference sequence includes only the images of HIP 17695—the star in our sample whose observations were closest in time and whose spectral type is the closest match for that of AU Mic. In addition to this classical RDI reduction, we also carry out an RDI/Karhunen–Loève Image Projection (KLIP) reduction as implemented in *SpaceKLIP* (Kammerer et al. 2022), utilizing the full frame and again using the images of HIP 17695 as the reference data.<fnref n="18" /> Because the NIRCam diffraction pattern changes based on the alignment between the target star and the coronagraph, techniques like KLIP are expected to more effectively eliminate starlight by combining multiple reference images such that residual flux in the science image is minimized.

To better mitigate stellar residuals at small separations, we additionally apply a Locally Optimized Combination of Images (LOCI; Lafrenière et al. 2007) RDI procedure utilizing just two optimization/subtraction regions and two separate sets of reference images. The inner optimization region covers *r* < 20 pixels (1"26), with subtraction being performed over the same region, and incorporates reference images from the sequences for HIP 17695, 2MJ0443, G-7-34, LP-944-20, and TYC 5899. Meanwhile, the outer optimization region spans 10 < *r* < 35 (0"63–2"21; from roughly the IWA to where a background source begins to infringe for AU Mic), with the corresponding subtraction region covering *r* ≥ 20 pixels (1"26), and includes only reference images from HIP 17695 and G-7-34, which were found to be free of off-axis sources within the optimization region. This strategy allows us to benefit from a more diverse selection of reference images at small separations—where differences in the wave front or coronagraph alignment change the coronagraphic image more significantly—without being affected by the bright off-axis sources present at larger separations for many of the survey targets.

<footnote n="18">¹⁸ The RDI/KLIP reduction with *SpaceKLIP* performs the nominal centering strategy described in Carter et al. (2023) as part of the starlight subtraction procedure, but otherwise uses the same outlier-corrected (“cleaned”) data as the other reductions.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Starlight-subtracted and roll-combined final JWST/NIRCam images of AU Mic in both F356W (top) and F444W (bottom). The images have been oriented such that the assumed disk major axis is parallel to the x-axis. To highlight the relative brightness of the disk at F356W, the images are displayed with the same linear color stretch. Both images have been smoothed with a σ = 0.5 pixel Gaussian for presentation and are displayed within a 10" × 3"/5 FOV. The approximate coronagraphic inner working angle is indicated by the dashed black circle, while the approximate position of the occulted star is marked with a white star symbol.</figcaption></figure>

With KLIP or LOCI, a model of the starlight in an image is constructed from some combination of reference images to minimize the residuals with a science image. When circumstellar signal (CSS) is present in the science image, these techniques effectively identify the combination of reference images that best nulls the entire scene, including both stellar and circumstellar signal. Unlike for classical RDI, this results in a model of the starlight that is systematically brighter than it should be and that ultimately suppresses the throughput of all circumstellar flux in the science data—an effect referred to as “oversubtraction.” The significance of oversubtraction depends on the relation between the spatial distributions of the circumstellar and stellar flux such that it is often negligible for point sources. However, the effect is often severe for extended sources—suppressing a significant fraction of flux while also introducing color offsets and inducing or obfuscating morphological features.

To attain image products free of significant oversubtraction, we also carry out model-constrained RDI (MCRDI), as described in Lawson et al. (2022), in which a model of the CSS is optimized alongside the model of the starlight. By subtracting the CSS estimate from the data during optimization of the stellar model, the contribution to the stellar model resulting from the projection of the underlying disk signal onto the reference images (or reference eigenimages, in the case of KLIP), the cause of RDI oversubtraction, can be eliminated. Effectively, this technique seeks the models of the stellar and circumstellar light that best explain the data—then subtracts the resulting stellar model to isolate the circumstellar flux. A description of the utilized synthetic disk model, its parameters, and the optimization procedure are provided in Appendix A. With the exception of the model-based “constraint,” this reduction is identical to the LOCI RDI reduction described previously.

The results of the MCRDI procedure are displayed in Figure 1, while the MCRDI images, the corresponding model constraints, and the residuals are shown in Figure 2. The results for each of the four reductions in both filters are compared in Figure 3.

# 4. Detection of the AU Mic Debris Disk

For both filters and all four starlight subtraction methods, the AU Mic debris disk is unambiguously recovered from the NIRCam observations (Figure 3)—marking the first reported images of the disk at these wavelengths (~3–5 µm). The appearance of the disk in these data is generally consistent with the appearance at shorter wavelengths. Extended signals above and below the disk midplane (e.g., Figure 3) are aliases of the disk produced by the off-axis “lobes” of the NIRCam diffraction pattern. Modeling of AU Mic’s spectral energy distribution in Matthews et al. (2015) found that the disk component is consistent with a 53K blackbody having a fractional luminosity of 3.5e-4, for which emission would peak at ~55 µm. Based on this, the observed light from the disk at the 3–5 µm can be assumed to be overwhelmingly composed of scattered starlight, with thermal emission from the disk being negligible.

Comparing the results from the different starlight subtraction methods, there are some notable differences. For both the RDI/ KLIP and RDI/LOCI results, oversubtraction is evident from the predominantly negative background values and from the negative alias of the coronagraphic stellar diffraction pattern at small separations (r ≲ 2"). Though the classical RDI procedure avoids these systematic oversubtraction effects, at separations ≲ 2" the results suffer from the nonoptimized model of the diffracted starlight, which cannot account for changes to the pattern between the science and reference data (largely due to differences in coronagraph alignment). MCRDI provides the same high fidelity as classical RDI at larger separations while permitting a much cleaner subtraction of the stellar diffraction pattern at small separations.<fnref n="19" />

We note that the differences between the KLIP and LOCI products are predominantly the result of the distinct optimization zones utilized (with KLIP using the full frame and LOCI using two narrow annular regions). A LOCI procedure

<footnote n="19">¹⁹ See Lawson et al. (2022) for discussion of additional considerations that may be relevant for broader applications of the MCRDI technique.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** For each filter (rows), the roll-averaged MCRDI final image (left), the best-fitting convolved model used as the constraint for MCRDI (middle), and the difference between the MCRDI result and the model image (right). The residuals use the same linear color stretch as the corresponding images to highlight that any differences are small compared to the total disk signal. A background source is visible in the lower left of the MCRDI images, at roughly (−1", −3").</figcaption></figure>

utilizing the same regions as KLIP results in nearly indistinguishable results.

The relative performance of classical RDI might be improved with the use of a dedicated reference sequence that adopts the small grid dither strategy (SGDS). Because the changes to the diffraction pattern between the science and reference sequences are generally dominated by differences in coronagraph alignment between the targets (e.g., Girard et al. 2022), observing the reference target at multiple dithers may increase the likelihood that the reference data includes a pointing that well-matches that of a given science pointing. However, as the dither offsets for NIRCam's SGDS are fixed (and the dithering is precise; Girard et al. 2022), the value of this approach will be dictated by the accuracy of NIRCam's pointing; if the size of the dithers is significantly larger than the typical pointing error, it is unlikely that any of the dither offsets will provide an improved match for the science image. As Girard et al. (2022) note the possibility of further improvements for NIRCam target acquisition, the utility of the SGDS for typical classical RDI reductions is currently uncertain. Alternatively, SGDS reference data might be better leveraged by some misalignment-aware classical RDI procedure—in which a model of the starlight is constructed by combining multiple reference dithers to best emulate the misalignment of the science data.

## 5. Companion Detection Limits

Inspecting the final images and disk-model-subtracted images, we identify no candidate companions near the plane of the disk. For reference, disk-model-subtracted images for MCRDI and classical RDI reductions of the F444W

observations are shown in Figure 4. A number of point sources appear within the field of view, but these are consistent with previously observed background sources based on AU Mic's proper motion (assuming Gaia DR3 values of $\mu_{\alpha} = 281$ mas yr$^{-1}$ and $\mu_{\delta} = -360$ mas yr$^{-1}$; Gaia Collaboration et al. 2023). For example, adjusting the locations of the two off-axis sources appearing in Figure 4 for proper motion, the results fall within ~50 mas of the locations of sources appearing in 2017 HST/STIS coronagraphy of AU Mic (originally published in Wisniewski et al. 2019; Grady et al. 2020). With no plausible companions detected, we proceed with analysis of the sensitivity of the data to assess both the companions we would likely have detected and what companions might still remain.

Computation of conventional contrast curves is complicated in this case by the fact that the flux of the occulted target star, AU Mic, cannot be directly measured from these data. While a synthetic stellar spectrum could be adopted for AU Mic's photosphere or otherwise used to convert existing photometry (e.g., for similar WISE filters) to NIRCam photometry, this nevertheless introduces an additional layer of model dependence. In light of this, we focus on the more direct calculation of companion detection limits in terms of companion brightness.

In this process, we determine limiting companion fluxes as a function of projected separation using forward modeling. Because companions are generally expected to orbit within or near the plane of the disk, we exclusively consider coplanar companions here (assuming a position angle of 128°.48 and an inclination of 90° for simplicity). For an array of projected separations spanning approximately 0".1–12" along both the southeastern and northwestern extents of the disk, we compute the corresponding position for each roll of the AU Mic data. Then, we use the *WebbPSF_ext* package to generate a PSF at
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Starlight-subtracted and roll-combined final images of AU Mic for the F356W (top row) and F444W (bottom row) filters using a variety of starlight subtraction techniques (as indicated by the column labels and described in 3.3). The black plus sign and dashed circle indicate the approximate coronagraph center and IWA, respectively. Interior to the IWA, coronagraph transmission decreases and is at least partially responsible for the intensity deficit. The images for a given filter are displayed with the same linear color stretch, as indicated by the colorbar at the end of each row. A mapping symmetric about zero is chosen to highlight the effects of RDI oversubtraction and speckle mismatch where present.</figcaption></figure>

each position using the closest OPD map and configured to match the observing configuration of the data. Each synthetic companion image is normalized to have a total flux of 1 mJy times the coronagraph transmission at the companion’s location (i.e., normalized to correspond to a source whose unocculted flux is 1 mJy). For each roll angle, our procedure considers the fit coronagraph offset from Section 3.2 when generating the model images. Subsequently, we simulate the effects of RDI on these sequences—as typical for RDI/LOCI forward modeling (e.g., Currie et al. 2019)—ultimately producing an oversubtracted and roll-combined model image. Additionally, we consider a scenario in which MCRDI is used to improve throughput (primarily at small separations, where RDI oversubtraction more significantly attenuates companion flux). For simplicity, we assume that MCRDI is able to effectively suppress all oversubtraction introduced by the companion, such that the coronagraph is the only significant limitation on companion throughput.

To determine the 1σ noise level in our data, we generate radial noise profiles for the MCRDI reductions of both filters. For this, we first replace the value of each pixel in the image with the sum of the values within an FWHM-diameter aperture (2.29 and 2.84 pixels for F356W and F444W, respectively). Then, for each radial separation from the star, $r_i$, the corresponding noise level is taken to be the finite-element-corrected standard deviation (Mawet et al. 2014) of all pixels with stellocentric separations, $r$, falling in ($r_i - 0.5 \cdot \text{FWHM}$) < $r$ < ($r_i + 0.5 \cdot \text{FWHM}$). To mitigate the impact of the disk, which would otherwise serve to artificially inflate the noise levels, we subtract the best-fitting MCRDI disk model (see Section A) before measuring the noise. Additionally, we exclude the vicinity of five background sources from this calculation. The signal-to-noise per resolution element (SNRE) map for the forward-modeled synthetic companion

image is then taken to be the aperture-summed companion image divided by the noise map computed for the data. Since the maximum SNRE achieved across a companion’s PSF may not be coincident with the injected candidate position—because the PSF spans a range of radial separations (and thus a range of noise levels)—we do not simply adopt the value at the injected position as the throughput SNRE. Instead, we adopt the largest SNRE value within 12 times the FWHM of the candidate position.<fnref n="20" /> This value is the throughput SNRE for a companion of 1 mJy at the prescribed location. Because the forward-modeled flux for a companion model scales linearly with the input flux, the 5σ limiting flux in mJy is taken to be 5 divided by the peak throughput SNRE of the 1 mJy model. When computed this way, the limits along the northwestern and southeastern extents of the disk are very similar (with small differences resulting from the asymmetry of the PSF and the misalignment of the coronagraph). Thus, for simplicity, we average the two values for a given separation when presenting sensitivity curves.

We note that the procedure outlined above does not account for residual disk flux resulting from inaccuracy of our disk model. Though these residuals affect our noise calculation (increasing our measured noise), they are not considered in the measurement of the “signal” component of the forward-modeled SNRE. In practice, where our disk model is brighter than the true underlying disk in the data (negative residuals), the flux of any coincident companions would appear to be diminished—and vice versa. Counterintuitively, this could result in a companion below the limits we report apparently manifesting with sufficient SNRE to be “recovered.” Given the faintness of the disk residuals (see Figures 2 and 4), we make

<footnote n="20">²⁰ This radius is chosen to include the six bright off-axis “nodes” of the NIRCam PSF.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The absolute value of disk-model-subtracted images for MCRDI and classical RDI reductions in F444W, oriented as in Figure 1. Each image has been multiplied by the projected radial separation in units of arcsec and is presented with a linear color stretch. The dotted black contour outlines the extent of the disk model, with two previously identified background sources indicated by the blue and red dotted circles. As a reference for the depth of the imagery, based on Linder et al. (2019) models, the brightness of the background source circled in red would correspond to a foreground companion with a mass comparable to Saturn's, while the brightness of the source circled in blue would correspond to a mass of ~0.12 $M_J$.</figcaption></figure>

no effort to account for this contribution quantitatively in this work.

The use of the MCRDI reduction for measuring the noise profile means that the detection limits achieved for “RDI” are representative of a reduction using MCRDI to mitigate oversubtraction from the disk, but making no effort to suppress oversubtraction for the considered off-axis companion. For a LOCI or KLIP reduction that does not account for the presence of the disk, simply forward modeling point sources would not be appropriate for assessing the faintest recoverable companions. If not suppressed, the disk will serve to significantly reduce the throughput circumstellar flux throughout the field of view (FOV)—including that of any companions (e.g., Lawson et al. 2022). To accurately account for the effects of starlight subtraction using any technique in which the stellar model is optimized by comparison with the data itself or in which the target data are used to build the PSF model (e.g., any common technique besides classical RDI), the entire circumstellar scene must be considered holistically in forward modeling. To address this, an “injection-recovery” approach (e.g., Carter et al. 2023) could be used in some scenarios. However, for an edge-on disk system like AU Mic, any position at which a companion would likely manifest is coincident with a non-negligible quantity of disk flux. So, while injection-recovery would account for the effects of the rest of the circumstellar scene on oversubtraction (and self-subtraction, in the case of angular differential imaging), we would still require some means of disentangling companion and disk signal. At some point, a model of the disk must be assumed in order to assess the sensitivity of these data to companions.

To map fluxes in the JWST/NIRCam filters to companion masses, we use the *species* Python package (Stolker et al. 2020). In *species*, companion properties (luminosity, surface gravity, effective temperature, etc.) are interpolated from grids of synthetic isochrones for a given age (assuming an age of 24 Myr for AU Mic) and companion mass. Synthetic spectra are then interpolated to produce a spectrum corresponding to these companion properties. Finally, *species* extracts

photometry from the resulting spectrum. For masses of 0.6 $M_J$ and above, we use the AMES-Cond isochrones from Allard et al. (2001). For lower-mass objects, none of the isochrone sets built into *species* have coverage at the age of AU Mic. Instead, we manually introduce the cloud-free petitCODE isochrones for low-mass objects provided by Linder et al. (2019) for masses of 0.5 $M_J$ and below. In either case, the isochrone properties are paired with AMES-Cond synthetic spectra (Allard et al. 2001) to estimate would-be companions' photometry. This ultimately provides valid synthetic photometry for masses as low as ~0.03 $M_J$. Comparing the photometry from *species* using Linder et al. (2019) isochrones and AMES-Cond spectra with the precomputed photometry for the Linder et al. (2019) isochrones shows F444W fluxes that are comparable but F356W fluxes that are much fainter for the precomputed photometry. As such, the companion detection limits for F444W are likely the more robust of the two filters. As the F444W filter is the more sensitive of the two to low-mass companions, this does not directly impact the faintest recoverable companions.

To convert flux limits to approximate contrasts, we scale a synthetic photosphere approximating AU Mic's mass (0.5 $M_\odot$) and age (24 Myr) to match 2MASS J- and H-band photometry of AU Mic (Skrutskie et al. 2006), and then extract photometry from the resulting spectrum. This yields approximate F356W and F444W fluxes of 4428 mJy and 3195 mJy respectively. Using photometry from the F335M target acquisition images (multiplied by a factor of 516 for use of the neutral density square; WebbPSF_ext; Leisenring 2021) to rescale the same synthetic model (in place of 2MASS photometry) yields very similar F356W and F444W fluxes of 4447 mJy and 3209 mJy, respectively. The 5$\sigma$ flux, mass, and contrast limits are presented in Figure 5.

We then compute maps of companion detection rates (“tongue plots”) for the MCRDI reductions using this information and assuming companions following circular orbits that are coplanar to the disk—which is assumed to be optically thin. For each synthetic companion mass, we consider a
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Top: 5σ companion detection limits in mJy for both filters—computed as described in Section 5. The right y-axis provides the corresponding mass limits for each filter in Jupiter masses—using AMES-Cond and cloud-free petitCODE isochrones as outlined in Section 5. Middle: 5σ contrast curves for the same reductions, adopting the model-based stellar fluxes for AU Mic indicated in the upper right corner (see Section 5). Bottom: the ratio of the companion throughput for the RDI-LOCI and MCRDI procedures. The improvement for MCRDI persists up to approximately the edge of the outer optimization region used for construction of the starlight model (r = 2.2/205 = 35 pixels). For all panels, the x-axis indicates the projected stellocentric separation of the measurement—in arcsec on the lower x-axis and in au on the upper x-axis. The relatively gentle slope of the coronagraph’s transmission profile permits some sensitivity to sources at smaller separations than the IWA of 0.63 (dotted gray line).</figcaption></figure>

logarithmic grid of semimajor axes with 300 values from 1 to 1000 au. For each semimajor axis and mass, we generate a sample every 0.01 over a full orbit, determine the projected separation for each sample, and then linearly interpolate the SNRE for that sample from the previously computed 5σ detection limits. This SNRE is then mapped to a detection probability assuming a normal distribution (e.g., that a 3σ detection corresponds to a 99.7% chance of detection). The overall detection rate for a companion at a given separation and mass is taken to be the average of the 36,000 samples over the full orbit. The companion detection maps for both filters are shown in Figure 6. We emphasize that these are not presented in terms of projected separations, but rather the true semimajor axes for companions on circular orbits.

5.1. Validation with Companion Injection

As a proof of concept for these detection limits, we injected two low-mass companions into the unsubtracted Stage 2 F444W products for AU Mic: a 0.1 *M*<sub>J</sub> planet at a projected separation of 2.71 along the southeastern extent of the disk, and a 0.2 *M*<sub>J</sub> planet at a projected separation of 1.72 along the northwestern extent of the disk. For this purpose, we use PSFs from *WebbPSF_ext* and again assume fluxes based on Linder et al. (2019) isochrones with AMES-Cond spectra (via *species*). The positions were drawn from circular coplanar orbits with semimajor axes of 30 and 20 au, respectively, in order to place the companions at projected separations where
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Maps of orbit-averaged companion detection probabilities as a function of companion semimajor axis (au) and companion mass (Jupiter masses on the left axis and Earth masses on the right). For the F444W observations, the map indicates detection probabilities for analogs of Jupiter, Saturn, and Neptune (i.e., matching semimajor axis and mass) of 0.73, 0.76, and 0.70, respectively.</figcaption></figure>

they would be slightly brighter than the detection limits of Figure 5.

We then reran the full MCRDI disk model optimization procedure exactly as for the real data, to verify that these companions would be recoverable in the resulting disk-model-subtracted results. During MCRDI optimization, we made no effort to mask the companions—meaning that the resulting disk model slightly overestimates the true disk’s brightness (by ~2%) in order to compensate for the companion flux. As both companions are visible in an initial “unconstrained” RDI/LOCI reduction, the final result could be improved by masking their locations throughout the procedure. Rerunning the procedure while masking the region within five times the FWHM of each companion had no perceptible impact on the final MCRDI image (<0.1% average change in brightness; the overbright disk in the unmasked scenario compensates for the companions’ flux in the least-squares optimization of the stellar

PSF model), but it does mitigate the overbrightening of the disk model. As such, there is a small difference in the brightness of the companions when the disk model is subtracted from the final MCRDI image in each scenario. For this demonstration, we proceed with the simpler unmasked version. Computing SNRE as before, the 0.1 *M*<sub>J</sub> and 0.2 *M*<sub>J</sub> companions have peak SNRE values of 6.6 and 9.5, respectively—manifesting with >5σ significance, as anticipated. The resulting MCRDI image and the disk-model-subtracted SNRE map are shown in Figure 7.

When applying MCRDI to data with visible companion candidates in nominal reductions, point sources could also be explicitly added to the MCRDI model instead. Though this increases the complexity of the optimization problem compared to simply masking them, it will likely improve the result when a candidate lies well within the speckle-limited regime (i.e., near the IWA), and it has the benefit of providing model-based measurements of companion fluxes as a byproduct. In application to data with companions that are much brighter relative to the disk, if not somehow addressed, the effect of the companions on the disk model will be more significant. In such a scenario, either masking or explicitly including the companions in the MCRDI CSS model is recommended.

## 5.2. Discussion of Companion Detection Limits

The limiting companion flux of ~0.004 mJy at larger separations in Figure 5 corresponds to 5σ flux contrasts of approximately 9.0 × 10<sup>-7</sup> and 1.3 × 10<sup>-6</sup> for F356W and F444W, respectively. These contrasts are consistent with the contrasts reported in, e.g., Carter et al. (2023) and Girard et al. (2022). In terms of companion masses, these contrasts permit exceptionally deep limits as a result of the relatively low luminosity of the parent star—with F444W showing sensitivity to <0.1 *M*<sub>J</sub> companions. The detection maps computed using forward modeling (Figure 6) show, for example, that for the F444W filter, our observations have a ~92% chance of detecting any Saturn mass (~0.3 *M*<sub>J</sub>) companion with a semimajor axis of 25 au (falling just inside the disk’s fiducial radius).

Based on modeling of the AU Mic disk in ALMA dust continuum emission, Pearce et al. (2022) and Vizgan et al. (2022) each provide estimates for a yet-unseen companion that might be responsible for shaping the inner edge of the disk. Pearce et al. (2022) find a companion mass and semimajor axis of ~0.44 *M*<sub>J</sub> and ~21.9 au (from an inner disk radius of 28.7 au), while Vizgan et al. (2022) find a companion mass and semimajor axis of ~0.34 *M*<sub>J</sub> and ~17 au (from an inner disk radius of 22.1 au). Permitting the aforementioned assumptions, our analysis indicates that we would have overall detection probabilities of ~91% and ~87% (respectively) for these planets in F444W. If such a planet remains, it would most likely be at very small angular separations at the epoch of these data—as each of the proposed objects would manifest with SNRE ≥ 5 until r ≲ 0.6 au.

Sezestre et al. (2017) conducted a dynamical analysis of a possible unseen “parent body” responsible for the fast-moving structures reported in Boccaletti et al. (2015)—exploring models for both an orbiting parent body (e.g., a planet) and a static parent body (e.g., a dust cloud resulting from a large planetesimal collision). The best-fit solution in the case of the orbiting parent body scenario would manifest at a projected separation of ~4.6 au (0.6"47) to the northwest at the epoch of
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Final F444W MCRDI products following the injection of two low-mass companions into observations of AU Mic (see Section 5.1). The images are oriented such that the disk major axis is aligned with the horizontal axis. Left: the final MCRDI image, which has been multiplied by the projected radial separation in units of arcsec to improve the visualization. Right: the SNRE map for the residuals between the MCRDI image and the MCRDI disk model. For both panels, the location of each injected companion is indicated by a white caret marker and annotation. Major and minor axis ticks delineate separation in increments of 1" and 0.5" respectively. A linear color stretch is used for both images.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Surface brightness measurements along the spine of AU Mic's disk in F356W (solid lines) and F444W (dashed lines) measured in circular apertures ($r = 4$ pixels $= 0.0252$) for images produced by three distinct starlight subtraction techniques: MCRDI (blue), RDI/LKLIP (orange), and classical RDI (green). The x-axis indicates the stellocentric separation projected onto the assumed disk major axis, with negative values corresponding to the southeastern extent of the disk. The roll-averaged coronagraph transmission profile is displayed as a light gray line (measured at the same locations and using the same aperture size), with values given along the right y-axis. To improve readability, error bars are plotted for alternating surface brightness measurements.</figcaption></figure>

the NIRCam data presented here (see right panel of Figure 8 in Sezestre et al. 2017). For F444W, a companion at this location would have been detected at $\ge 5\sigma$ for masses $\gtrsim 0.54 M_J$—placing a strong upper limit on the mass of the parent body in this scenario. Overall, future studies combining an analysis similar to that of Sezestre et al. (2017) with the detection limits provided herein should enable a much clearer understanding of the source of the fast-moving features observed in AU Mic's disk.

The nondetection of outer giant planets and these deep detection limits suggest the likely absence of outer giant planets in the AU Mic system. This echoes the results of Daley et al. (2019), who analyzed the vertical dust distribution of AU Mic's disk and ultimately inferred that its thin vertical size should preclude the presence of perturbers more massive than $\sim 2 M_\oplus$ in the outer disk. Given the lack of outer giant planets, which would otherwise prevent volatiles from streaming inward (e.g., Clement et al. 2022), transmission spectra for the inner planets should be relatively likely to reveal atmospheres that are rich in water or other volatiles.

## 6. Disk Surface Brightness and Color

Here, we provide measurements of surface brightness (SB) and color for AU Mic's disk from reductions using classical RDI, RDI/LKLIP, and MCRDI. A subsequent study focusing on analysis of the disk in these data will further investigate these measurements and their implications.

We measure SB within circular apertures with radii of 4 pixels (0.0252) placed along the spine of the disk. For simplicity, the spine positions are taken to be the analytic spine for a ring-like forward-scattering disk with a fiducial radius of 35 au, an inclination of 89°, and a position angle of 128°48 (i.e., falling along an ellipse for $r_{proj} < 35$ au, and along the major axis otherwise). These measurements are made using the derotated and roll-averaged final image products.<fnref n="21" />

Uncertainties on these measurements are estimated by making a set of like measurements at the same projected separation but differing position angles (in increments of 3° in the range [0°, 357°], for 119 total sets), and then calculating the median absolute deviation among the values at each position.

<footnote n="21">²¹ The results showed no statistically significant differences when instead averaging surface brightness measurements from individual rolls without derotation.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Surface brightness profiles for the MCRDI reductions (F356W in blue and F444W in red) presented in a log-log scale and including separations to 6". The roll-averaged coronagraph transmission profiles are displayed in light gray, with values given along the right y-axis. For both surface brightness and coronagraph transmission, the northwestern (NW) and southeastern (SE) profiles are indicated by solid and dashed lines (respectively). Error bars are excluded to improve readability.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** A comparison of the surface brightness profiles for the northwestern (NW) and southeastern (SE) extents of the disk for the MCRDI reductions (blue) in magnitudes. In the upper panel, profiles for the convolved disk models (black), whose unconvolved images are perfectly symmetric, are included to highlight effects artificially induced by coronagraph misalignment and by the asymmetry of the PSF itself. In the lower panel, the artificial asymmetry has been corrected by subtracting the model asymmetry measurements of the upper panel from those of the MCRDI measurements. For context, the dashed-dotted purple line shows the expected asymmetry based on the exponential fits to HST/STIS surface brightness profiles for AU Mic's disk from Schneider et al. (2014). In both panels, error bars are offset slightly in the x-axis direction to improve readability.</figcaption></figure>

To mitigate the effect of the disk itself on the uncertainty estimates, the disk model used with MCRDI is subtracted from each image before the measurements are made (in the case of the RDI/KLIP reductions, the disk model is first forward modeled to induce appropriate oversubtraction). Generally, the resulting uncertainties are much larger than the differences between the MCRDI and classical RDI measurements beyond the small separations where classical RDI suffers from stellar residuals. Meanwhile, the RDI/KLIP measurements are affected by oversubtraction—which is not considered in these

uncertainties—leading to measurements that are significantly discrepant with those of the other reductions.

SB profiles for both filters are shown in Figure 8. For the MCRDI reductions, we also provide (a) a log-log plot of surface brightness in Figure 9, and (b) a comparison of the brightness between the southeastern (SE) and northwestern (NW) extents of the disk in Figure 10.

To measure the disk's F356W versus F444W color, we adopt the model-based stellar fluxes used to estimate contrast in Section 5: $F_{356}^* = 4428$ mJy and $F_{444}^* = 3195$ mJy. Denoting a
<figure><img src="image_14.png" /><figcaption>Figure 11. Measurements of AU Mic's disk color as described in Section 6. Measurements for the MCRDI disk model (black) are included for both the convolved (solid) and unconvolved (dotted) model images, to provide a qualitative estimate of the effect of the differing PSFs on the measured color. To improve readability, (a) error bars are plotted for every other measurement and (b) points for the three data reductions are offset slightly from one another in the x-axis direction.</figcaption></figure>

measurement of the disk SB with superscript *d*, the disk color in magnitudes is computed as the difference between the SB color and the stellar color:

$$ \begin{aligned} \Delta(m_{356} - m_{444}) &= -2.5 \left[ \log_{10} \frac{F_{356}^d}{F_{444}^d} - \log_{10} \frac{F_{356}^*}{F_{444}^*} \right] \\ &= -2.5 \log_{10} \frac{F_{356}^d F_{444}^*}{F_{444}^d F_{356}^*}. \end{aligned} \quad (1) $$

The uncertainties for disk color measurements assume zero uncertainty for the stellar fluxes, such that all color uncertainty is contributed by the uncertainty for the disk SB. Disk color as a function of projected separation is presented in Figure 11.

We remark that these measurements have not been corrected for the effects of the PSF and coronagraph beyond the typical JWST flux calibration strategy. As such, a brightness measurement at a given position may be contaminated significantly by the wings of the PSFs for other parts of the disk (e.g., as in the Early Release Science observations of HD 141569 A with JWST MIRI; Hinkley et al. 2023).

## 6.1. Discussion of Disk Surface Brightness and Color

The surface brightness profiles presented in Figures 8 and 9 manifest with shapes similar to those presented in shorter-wavelength scattered-light studies—albeit smoothed by the wider 3–5 µm PSF. Most visibly in Figure 9, the profiles show a modest slope from $r \sim 1''–3''$, which steepens beyond $r \sim 3''–4''$ (the ring’s peak density radius)—as is noted in many prior studies of the system (e.g., Liu et al. 2004; Krist et al. 2005; Fitzgerald et al. 2007). While the brightness profiles show a marked flattening within $\sim 1''$, we note that this region is significantly affected by the coronagraph—such that much or all of this apparent change is likely instrumental rather than astrophysical (see Appendix B for additional context).

Taken at face value, the measurements presented in Figures 9 and 10 seem to suggest that the NW side of the disk is brighter for $r \lesssim 2''$—qualitatively consistent with measurements at shorter wavelengths (e.g., Fitzgerald et al. 2007; Schneider et al. 2014; Boccaletti et al. 2018). However, this apparent asymmetry is significantly smaller than the asymmetry in HST/STIS imagery of the disk reported by Schneider et al. (2014). For comparison, the asymmetry for the exponential fit to the SB measurements for the NW and SE extents of the disk, provided in Figure 39 of Schneider et al. (2014), is included in Figure 10. Additionally, comparison with the measurements of the truly symmetric MCRDI disk model,

also displayed in Figure 10, reveals that the two sets of measurements are generally consistent with one another at small separations.<fnref n="22" /> This suggests, instead, that much of the apparent asymmetry is induced by instrumental effects. Additional testing summarized in Appendix B shows that the misalignment of the coronagraph (see Section 3.2) can explain the majority of the asymmetry at small separations, with some additional asymmetry resulting from the asymmetry of the NIRCam PSF itself. The lower panel of Figure 10 shows the result when correcting the on-sky asymmetry measurements for the asymmetry induced in the model. Measuring the overall surface brightness of the inner disk using a large rectangular aperture spanning 0".63–2" in the major axis direction and $-0''$.5–$-0''$.5 in the minor axis direction—and correcting for instrumentally induced asymmetry in the same manner—yields overall inner disk asymmetry measurements of $-0.007 \pm 0.019$ mags for F356W and $-0.014 \pm 0.023$ mags for F444W.<fnref n="23" /> Overall, it appears likely that the enhanced brightness of the NW disk observed at shorter wavelengths is absent or significantly diminished at these wavelengths.

At larger separations, a seemingly localized increase in brightness on the SE side for both filters pushes the SE brightness above that of the NW at $r \sim 2''$.5, with a second lower-significance enhancement at $r \sim 3''$.25. Notably, the former roughly coincides with the projected location of the SE4 feature from Boccaletti et al. (2018) for the epoch of these data. While these data do not generally have sufficient spatial resolution to study these features in detail, this may be suggestive that NIRCam can identify them through such brightness enhancements. The recurrence of this feature in both filters lends further credence to this notion.

The colors measured from the classical RDI and MCRDI reductions (Figure 11) are generally consistent—indicating an approximately flat or slowly changing blue color of approximately $-0.3$ mags over the field of view, with a low-significance blue dip near 4". The anticipated small-separation stellar residuals for the classical RDI data manifest here as noisy/erratic measurements within $\sim 2''$. A key risk for using

<footnote n="22">²² A low-significance deviation at $r \sim 1''$.25 is apparent in the F356W profile of Figure 10. However, given the lack of recurrence in F444W, we conclude that it is not likely to be astrophysical—or, if astrophysical, it is much brighter in F356W and thus unlikely to be a companion (which should be much brighter in F444W).</footnote>

<footnote n="23">²³ The uncertainties for the corrected asymmetry measurements do not account for additional (but difficult to quantify) uncertainties introduced by this correction strategy—such as the accuracy of the determined coronagraph misalignment, the accuracy of the synthetic PSFs, or the accuracy of the disk model itself—and thus should be considered optimistic estimates.</footnote>
<table><caption>**Table 2**
MCRDI Disk Model Parameters</caption><thead><tr><th>Filter</th><th>Incl (deg)</th><th>PA (deg)<sup>a</sup></th><th>r<sub>0</sub> (au)</th><th>h<sub>0</sub>/r<sub>0</sub></th><th>α<sub>in</sub></th><th>α<sub>out</sub></th><th>g<sub>1</sub></th><th>g<sub>2</sub></th><th>w<sub>1</sub></th></tr></thead><tbody><tr><td>F356W</td><td>89.20</td><td>128.48</td><td>33.49</td><td>0.038</td><td>2.82</td><td>-4.92</td><td>0.95</td><td>0.46</td><td>0.44</td></tr><tr><td>F444W</td><td>89.20</td><td>128.48</td><td>30.26</td><td>0.039</td><td>3.39</td><td>-4.85</td><td>0.98</td><td>0.47</td><td>0.65</td></tr></tbody></table>

**Notes.** Disk model parameter values used in the Model-constrained RDI reductions presented herein. These are provided solely for the sake of reproducibility and should not be considered robust measurements of disk parameters. See Appendix A for details.

oversubtracted results to study disks is highlighted by the RDI/ KLIP disk color profile in Figure 11, which shows measurements deviating significantly from those of classical RDI and MCRDI alongside suggestive sharp radial trends. Notably, (a) the disk flux lost to oversubtraction is the projection of the underlying disk image onto the reference images and (b) the width of the stellar diffraction pattern scales linearly with the observing wavelength, while the true distribution of disk flux very likely does not. As a result, oversubtraction typically varies both spatially and spectrally, and so will tend to induce varying offsets in measurements of disk color.

Comparison of the color profiles for the raw and convolved models indicates a color offset introduced by the convolution as a result of the differences between the two filters' PSFs. In galactic astronomy, when measuring colors and color gradients of galaxies, this effect is often circumvented by convolving each filter's images with the PSF of the other filter before measuring surface brightness profiles ("cross-convolution"; e.g., den Brok et al. 2011). This produces images at differing wavelengths in which flux from sources is smeared comparably —but at the cost of degrading the effective spatial resolution of the data. Applying this approach to our MCRDI-reduced images, the result is much more similar to the measurements for the raw (unconvolved) model, but shows a more uniform and gentle slope: with disk color values spanning roughly -0.1 to -0.25 mags from 1" to 4". Deconvolution techniques, such as Richardson-Lucy deconvolution (Richardson 1972; Lucy 1974), may provide a solution that avoids degrading the spatial resolution of the data. However, in our case, these methods are significantly challenged by both the spatial variations of the NIRCam PSF and the use of a coronagraph. Ultimately, we leave a more detailed assessment of possible solutions for future work, and provide the current measurements as a preliminary assessment.

Though many prior scattered-light studies of the disk also found a blue disk color (e.g., Krist et al. 2005; Fitzgerald et al. 2007; Lomax et al. 2018), we emphasize that this is the first time the disk's 3-5 µm color has been measured. The NIRCam color probes a distinct region of the scattered-light spectrum and so is sensitive to different properties of the disk material. For example, optical color is much more effective for assessing the minimum grain size for AU Mic's disk (~0.2 µm; Arnold et al. 2022). Assuming a composition of standard astronomical silicates and a grain size distribution otherwise as the best-fit ADP solution in Arnold et al. (2022), the difference between a minimum grain size of 0.1 and 0.2 µm is a >1 mag change in the 0.4-0.9 µm color, while the 3.5-4.4 µm color is negligibly affected (changing on the order of millimages). On the other hand, the 3-5 µm regime is rich with scattered-light spectral features (e.g., tholins, CO₂ ice, and H₂O ice; Rodigas et al. 2015; Chen et al. 2019; Kim et al. 2019), while the optical regime is predominantly featureless. As such, it is possible—if

not likely—that the blue NIRCam color is the result of a distinct mechanism from the blue color at shorter wavelengths. As illustrated in Arnold et al. (2022), assessing dust composition and size distribution in debris disks via multi-wavelength observations requires particular care. The extensive degeneracies that exist between various compositional mixes and size distributions mean that identifying a single strong solution is insufficient. Rather, it is necessary to conduct a thorough investigation of all reasonable solutions, such that some subset of solutions that are consistent with the data can be identified. With consideration for this, we reserve further exploration of these details for a follow-up study that considers these NIRCam measurements alongside those of prior studies.

## 7. Vertical Disk Structure

The projected vertical size of the disk (i.e., its apparent width in the minor axis direction) is primarily a function of the disk's orientation, vertical dust distribution, and radial dust distribution (assuming $i \neq 90^\circ$). As such, detailed modeling is necessary before this observable can be distilled to more physically valuable information regarding the disk's vertical structure, and thus to a better understanding of the mechanisms shaping the disk (e.g., Daley et al. 2019). Notably, the disk's projected vertical size must be resolved for such an investigation to be meaningful. Whether this is the case for the NIRCam data is not immediately clear, as the spatial resolution is comparable to the projected vertical scale expected for the disk based on prior observations. For example, using HST/ACS imagery of the disk (FWHM ~0.063 at F606W), Krist et al. (2005) measure a projected disk FWHM as small as ~0.022 (occurring at a projected separation of ~1.5), while the FWHMs of the F356W and F444W filters are 0.14 and 0.18, respectively.

To assess whether the projected vertical extent of the disk is resolved in NIRCam imagery, we proceed as follows. Along several vertical slices (perpendicular to the major axis), we measure the disk's surface brightness in the F356W MCRDI reduction. These measurements are made using rectangular apertures with width 3 pixels (0.189) and height 1.5 pixels (0.095)—with each profile then being normalized to have the same peak value. For comparison, we make like measurements of a “thin model” of the disk having negligible projected vertical size ($h_0/r_0 = 10^{-5}$, $\alpha_{in} = 10$, $\alpha_{out} = -10$, $\gamma = 2$; otherwise as given in Table 2)<fnref n="24" /> and convolved with (a) coronagraphic images from WebbPSF_ext and (b) an empirical field PSF<fnref n="25" /> extracted from a bright background source in the near-contemporaneous GTO 1184 observations of TYC 5899. If the projected vertical size of AU Mic's disk is

<footnote><sup>a</sup> The value passed into VIP is PA-180° due to a difference of convention.</footnote>

<footnote n="24">²⁴ See Appendix A for descriptions of these parameters.</footnote>

<footnote n="25">²⁵ WebbPSF and WebbPSF_ext do not account for a number of detector-level effects that can result in underestimation of the PSF width.</footnote>
resolved in the NIRCam data, a PSF-convolved image for a disk much thinner than is supported by prior observations should yield substantially narrower vertical profiles than the on-sky data. Likewise, if the measurements of the aforementioned thin model yield profiles of comparable width to the on-sky data, it would mean that the projected vertical size is not resolved.

Alongside these measurements, we also provide profiles for the instrumental PSF for reference. At each separation, we used *WebbPSF_ext* to generate a PSF at the location of the disk major axis for each roll angle and then derotated and combined the images as we did for the data. This results in a sequence-averaged PSF having FWHM of ~0.18 in vertical profiles measured using the aforementioned aperture.

The measured vertical brightness profiles are presented in Figure 12. Uncertainties ($1\sigma$) for the measurements of the data are depicted as shaded gray regions and were computed following the same procedure outlined in Section 6.

## 7.1. Discussion of Vertical Structure

As is evident from comparison of the profiles for the PSF and the thin disk model using *WebbPSF_ext* for convolution, even a disk of unresolved vertical scale will manifest with much greater width than the PSF, as a result of the extended nature of the disk flux and the shape of the NIRCam PSF. As such, comparison with profiles of the PSF alone is insufficient for assessing whether the vertical extent of the disk is resolved. However, the profiles for the data are significantly wider than the profiles for the thin model—for both the *WebbPSF_ext* and empirical PSF versions. This suggests that we are seeing disk signal that is at least marginally resolved in the minor axis direction. As further evidence of this, we note that our best-fit MCRDI disk model appears consistent with the projected vertical size of the disk reported in Krist et al. (2005). At 2", Krist et al. (2005) report a projected disk FWHM of 0.25 for HST/ACS F606W imagery. Convolving our raw F356W MCRDI disk model with a Gaussian kernel with FWHM = 0.063 to approximate the F606W PSF, we measure a very similar projected disk FWHM value of 0.24 at this separation.

## 8. Conclusions and Future Work

We have presented high-contrast coronagraphic imaging of the AU Mic debris disk system from JWST/NIRCam. Our key findings are summarized hereafter.

1. The disk is unambiguously recovered in both F356W (3.563 µm) and F444W (4.421 µm) from separations as small as ~0.3 (well inside of the IWA) to separations as large as ~5". These detections mark the first images of the disk at 3–5 µm.

2. Using the model-constrained RDI (MCRDI) technique for the removal of starlight, final images are free of both the systematic oversubtraction that plagues RDI/KLIP or RDI/LOCI products and the significant small-separation stellar residuals apparent in classical RDI products.

3. No companions were identified, but analysis indicates that the data were capable of uncovering planets as small as ~0.1 $M_J$ with 5$\sigma$ confidence beyond ~2" (~20 au). These deep constraints on the presence of wide-orbit, massive planets are unique for AU Mic and provide relevant context for understanding planet formation and

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Normalized vertical brightness profiles for the F356W images described in Section 7. The profiles are normalized to have the same peak value and are offset along the y-axis according to the major axis position of the slice. Shaded gray regions indicate the $1\sigma$ uncertainties on the measurements for the MCRDI profiles. The top panel shows the locations and widths of the six slices in the context of the MCRDI image (which has been multiplied by the projected radial separation for this visualization only).</figcaption></figure>
evolution in a system with a dynamic debris disk and a
compact, multiplanet system at ≪1 au separations.

4. The significant brightness asymmetry favoring the north-western side of the inner disk at shorter wavelengths is not evident at 3–5 µm—with a 1σ upper limit of ~0.04 mags on any asymmetry (versus the best-fit model with asymmetry of ~0.2 mags from Schneider et al. 2014).

5. In both filters, we see evidence of one or more localized brightness enhancements to the southeast that may correspond to previously identified fast-moving disk features (e.g., Boccaletti et al. 2015).

6. A blue disk color of roughly −0.3 mags is measured between F356W and F444W—likely corresponding to a color closer to −0.2 mags when PSF effects are considered.

7. The projected vertical size of the disk is slightly resolved in these data, showing an apparent FWHM of 0.′′42 at 1′′ in F356W, compared to 0.′′25–0.′′28 for a disk model having infinitesimal vertical size before convolution with the instrumental PSF.

The high-significance detection of AU Mic's disk at these wavelengths demonstrates the unique suitability of JWST and NIRCam for studying this system—and motivates further studies facilitated by JWST's wavelength coverage. Follow-up observations in additional NIRCam filters would permit the study of the presence and location of ices in the system by further probing the 2–5 µm scattered-light spectral features of ices alongside the surrounding continuum (e.g., Honda et al. 2009; Debes et al. 2013; Tazaki et al. 2021; Betti et al. 2022)—thus testing the inferred presence of water ice based on polarimetric imagery from HST/ACS in Graham et al. (2007). As the theoretical ice line for AU Mic falls at ~2 au (Schüppler et al. 2015), ices should be observed throughout the bulk of the resolvable disk. With the planned addition of dual-channel coronagraphic imaging for NIRCam, the presence of ices could be tested with just one additional sequence using, e.g., F210M and F300M.

As demonstrated in Rodigas et al. (2015), the addition of thermal IR data can substantially alter the interpretation of a disk compared to scattered-light data alone. Following from Gáspár et al. (2023), who reported mid-infrared imaging of the Fomalhaut debris disk with JWST’s Mid-Infrared Instrument (MIRI), imaging of AU Mic with MIRI would likewise probe AU Mic’s disk in the thermal regime—further constraining its composition and helping to clarify the factors shaping it.

Reobservation of the system with NIRCam in the F444W filter would serve to place stronger constraints on the outer planets that might still remain hidden. In many cases, projected locations very near the coronagraph center are the only orbital positions that cannot be ruled out with high confidence from the existing data. A second epoch even just a year after the first would provide a sufficient baseline for the few degrees of orbital motion needed for many of these would-be planets to become detectable.

Overall, the results presented herein highlight the strength of JWST/NIRCam for studying circumstellar disk systems. While ground-based facilities observing at shorter wavelengths may achieve better spatial resolutions and reach higher contrasts, JWST is the only observatory capable of studying many compelling targets in the ~2–5 µm regime. As these wavelengths coincide with extremely favorable planet contrasts and

span numerous notable scattered-light spectral features, this
capability has considerable scientific value.

Acknowledgments

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

We acknowledge the decades of immense effort that enabled the successful launch and commissioning of the JWST; these results were possible only through the concerted determination of thousands of people involved in the JWST mission. In particular, we offer gratitude to a number of individuals who enabled this study through contributions to either the 2002 NIRCam instrument proposal or to the development and commissioning of the NIRCam instrument: Martha Boyer, Daniel Eisenstein, Klaus Hodapp, Scott Horner, Doug Kelly, Don McCarthy, Karl Misselt, and George Rieke. Additionally, we thank those that carried out commissioning of the NIRCam coronagraphy mode, including Julien Girard, Jens Kammerer, Mario Gennaroa, Armin Rest, Eiichi Egami, Ben Sunnquista, and many others. We are grateful for support from NASA through the JWST NIRCam project, contract number NAS5-02105 (M. Rieke, University of Arizona, PI).

The JWST data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific observations analyzed can be accessed via 10.17909/x9hs-pr32. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support to MAST for these data is provided by the NASA Office of Space Science via grant NAG5-7584 and by other grants and contracts.

This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.

K. Lawson's research was supported by an appointment to the NASA Postdoctoral Program at the NASA-Goddard Space Flight Center, administered by Oak Ridge Associated Universities under contract with NASA.

E. Bogat's work was supported by a grant from the Seller's Exoplanet Environments Collaboration (SEEC) at NASA GSFC, administered through NASA's Internal Scientist Funding Model (ISFM).

*Software:* Matplotlib (Hunter 2007; Caswell et al. 2021), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Astropy (Astropy Collaboration et al. 2013, 2018), CuPy (Okuta et al. 2017), LMFIT (Newville et al. 2022), WebbPSF (Perrin et al. 2014), WebbPSF_ext (Leisenring 2021), Space-KLIP (Kammerer et al. 2022), Vortex Image Processing (Gomez Gonzalez et al. 2017).

Appendix A
MCRDI Disk Models and Optimization

To generate models for use in the model-constrained RDI procedure, we utilize the *scattered-light_disk* module of the *Vortex Image Processing* (VIP; Gomez Gonzalez et al. 2017) Python package, which introduces a “lite” version of the GRaTer disk modeling code (Augereau et al. 1999) that assumes an optically thin disk dominated by single scattering. Because the goal of our application is merely to superficially
simulate the distribution of light from the disk for the purpose of preventing RDI oversubtraction, the physical validity of this assumption is inconsequential. Though the AU Mic disk is famous for structure and asymmetries observed at shorter wavelengths, inspection of the other NIRCam data reductions suggested that any such features are much less pronounced or absent in these data (due to wavelength dependence or spatial resolution). Moreover, to eliminate oversubtraction with MCRDI, it is not necessary to perfectly reproduce every nuance of the distribution of CSS in the data. Rather, the CSS can be overestimated in some areas and underestimated in others—so long as the two balance in the least-squares construction of the model of the stellar diffraction pattern (i.e., such that the residuals between the true CSS and the estimated CSS yield a negligible projection onto the reference images). As such, we adopt a simple axisymmetric ring-like disk geometry.

Using VIP, we model the disk as a function of nine parameters and with a scattering phase function (SPF) that is the linear combination of two Henyey–Greenstein (H–G) SPFs (Henyey & Greenstein 1941). The models assume linear flaring ($\beta = 1$) and a Gaussian vertical density distribution ($\gamma = 2$). The varied parameters include:

1. inclination (incl, in degrees);

2. position angle (PA, in degrees);

3. fiducial radius ($r_0$, in au);

4. the ratio of scale height to radius ($h_0/r_0$);

5. the density power-law indices interior and exterior to the fiducial radius ($\alpha_{in}$ and $\alpha_{out}$, respectively);

6. the H–G asymmetry parameters, $g_1$ and $g_2$, and the weight for the first SPF ($w_1$, with $w_2 = 1 - w_1$).

We create each raw model oversampled by a factor of two relative to the data. We then rotate the model to the appropriate roll angles for the data and convolve it with synthetic NIRCam coronagraphic images—sampled from an array of detector positions. These images are generated with the `WebbPSF_ext` package and using the OPD map measured closest in time to the observations.

In lieu of precise prior knowledge of the disk’s spectrum (which will affect the resulting diffraction pattern), we instead adopt a synthetic spectrum matching the spectral type of the parent star. Effectively, this assumes that scattered light is dominant at these wavelengths and that this scattering lacks wavelength dependence (i.e., “gray scattering”). Though not precisely accurate, this approximation is sufficient for the purpose of the desired superficial disk estimate.

At this point, the synthetic coronagraphic images to be used for convolution could be normalized as in Section 5—by normalizing to the transmission of the coronagraph at each sampled position. This would permit a convolved image that effectively samples the changes to the diffraction pattern’s morphology and to the coronagraph transmission at the resolution of the sampled grid. To permit finer sampling of the transmission, we instead normalize each of the synthetic coronagraph images to sum to one, but multiply each model image with a coronagraph transmission map just before convolution—resulting in effective transmission sampling at the resolution of the oversampled pixels. We remark that this significantly improves the consistency of the final disk models with the data in this case; we strongly recommended this

strategy for any convolution scenarios where circumstellar signal is present near the IWA.

We then create the convolved model image for each roll angle by convolving every pixel of the rotated input model with the nearest sample from the grid in polar coordinates (see Figure 13). For each roll angle, we consider the previously computed offset of the coronagraph center from the star for both the coronagraphic transmission and for matching pixels with samples from `WebbPSF_ext`. Once convolution is completed, we resample the resulting images to match the sampling of the data.

Once a convolved model sequence is created for a particular set of input disk parameters, we forward model it for our standard LOCI RDI reduction in the typical manner (e.g., Currie et al. 2019) and compare it with the observed LOCI RDI result within a region of interest. This approach is mathematically identical to the direct optimization of a CSS estimate using MCRDI,<fnref n="26" /> but allows for post facto analytic rescaling of the model’s brightness to minimize residuals with the data—ultimately eliminating the need to vary the brightness of the disk alongside the other parameters.<fnref n="27" /> The region of interest includes pixels with either $r < 35$ pixels or falling within a stellocentric rectangular region rotated to the approximate position angle of the disk (PA = 128°48; Vizgan et al. 2022) with a width of 12 pixels and a length of 94 pixels (where the detection in the LOCI RDI images begins to wane)—but excluding the region within three pixels of the star (where small differences in coronagraph alignment can produce significant stellar residuals that might otherwise impact disk model fitting). This region of interest is intended to include the majority of the disk flux while also including the inner background region where oversubtraction is most evident in the unconstrained LOCI result. We note that the region falling beyond the optimization region used for RDI subtraction has no effect on the amount of oversubtraction and thus does not need to be considered for the purpose of eliminating oversubtraction. However, as we needed a model that was superficially accurate at larger separations for use in Sections 5 and 6, we chose to extend the region of interest as described to cover both needs.

To explore the resulting disk model parameter-space, we use the differential evolution algorithm (Storn & Price 1997)—a genetic algorithm for global optimization that performs well in degenerate and/or multimodal parameter spaces and typically identifies a strong solution in only a few thousand samples (e.g., Lawson et al. 2020). For optimization, all of the aforementioned parameters are varied in wide ranges—with the exception of position angle, which is fixed to 128°48. The model for each of the two filters was optimized separately with this approach—except that the F444W model was restricted to the best-fit inclination of the F356W model.

After the optimal model is identified in this manner, the scaled and convolved input model is used to perform MCRDI on the data to achieve the final MCRDI result. The MCRDI

<footnote n="26">²⁶ That is, computing the RDI residuals after first subtracting the CSS estimate from the data; see Lawson et al. (2022) for more information regarding this equivalency.</footnote>

<footnote n="27">²⁷ Because oversubtraction is entirely the result of the presence of the disk itself, the brightness of the forward-modeled result scales linearly with that of the input. For example, if $I_M$ is the forward-modeled result for an input disk model $I_M$, then $sI_M$ will be the result for the input model $sI_M$ for any scalar constant s. The value of s producing the best fit can therefore be computed after forward modeling.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 13.** For both rolls of the F356W AU Mic observations: a map showing the vicinity of the star (gold star marker) with each oversampled pixel colored according to the coronagraphic image sample (colored round markers) with which it is paired for model convolution. Solid vertical and horizontal lines denote the boundaries of detector-sampled pixels, while dotted lines delineate the oversampled pixel boundaries. It should be noted that not all of the *WebbPSF_ext* samples in the grid are ultimately matched to any pixels.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** As Figure 10, but for measurements of differing disk models and with color now encoding the filter (F356W in blue and F444W in red). Solid lines correspond to measurements of the convolved MCRDI model (as the black lines in Figure 10), the dashed-dotted lines are for the same model but assuming the star and coronagraph are perfectly aligned for convolution, and the dotted lines are for the model before convolution.</figcaption></figure>

result, the best-fit model estimate, and the residuals for each filter are shown in Figure 2.

In Table 2, we provide the best-fitting model parameters. However, we emphasize that these are included purely for the sake of repeatability. These values were not optimized to be physically meaningful, but rather to superficially emulate the distribution of disk flux following convolution. These values should not be used other than to reproduce the MCRDI reductions for these or very similar data.

## Appendix B: Instrumental Asymmetry Effects

To diagnose the cause of the induced asymmetry for the otherwise symmetric input MCRDI model discussed in Section 6, we make additional asymmetry measurements using the same input model but performing convolution assuming perfect alignment between the star and coronagraph for both rolls. Figure 14 compares these results with those of the nominal convolved model and with the unconvolved input model. This reveals that the majority of the small-separation asymmetry results from the misaligned coronagraph, which effectively blocks more of the disk on the SE side. The reach of this effect to separations as large as ~2" is driven by the width
of the NIRCam PSF combined with the sharp decline in disk brightness as separation increases—resulting in a non-negligible contribution from the bright portion of the disk that is occulted to the southeast but not occulted to the northwest.

The remaining asymmetry for the “aligned” model appears to result from the asymmetry of the PSF itself. Data simulated with parallactic angles offset from those of the data (but with the same roll offset) show broadly varying induced asymmetries. These asymmetries were minimized when the major axis of the disk was aligned with the y-axis direction of the detector —where the PSF core is more symmetrical—and never became significantly larger than the asymmetries induced at the observed parallactic angles.

## ORCID iDs

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

Joshua E. Schlieder https://orcid.org/0000-0001-5347-7062

Jarron M. Leisenring https://orcid.org/0000-0002-0834-6140

Ell Bogat https://orcid.org/0000-0002-7325-5990

Charles A. Beichman https://orcid.org/0000-0002-5627-5471

Geoffrey Bryden https://orcid.org/0000-0001-5966-837X

Andrés Gáspár https://orcid.org/0000-0001-8612-3236

Tyler D. Groff https://orcid.org/0000-0001-5978-3247

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

Michael R. Meyer https://orcid.org/0000-0003-1227-3084

Thomas Barclay https://orcid.org/0000-0001-7139-2724

Per Calissendorff https://orcid.org/0000-0002-5335-0616

Matthew De Furio https://orcid.org/0000-0003-1863-4960

Marie Ygouf https://orcid.org/0000-0001-7591-2731

Anthony Boccaletti https://orcid.org/0000-0001-9353-2724

Thomas P. Greene https://orcid.org/0000-0002-8963-8056

Peter Plavchan https://orcid.org/0000-0002-8864-1667

Marcia J. Rieke https://orcid.org/0000-0002-7893-6170

Thomas L. Roellig https://orcid.org/0000-0002-6730-5410

John Stansberry https://orcid.org/0000-0003-2434-5225

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

Erick T. Young https://orcid.org/0000-0002-6395-4296

## References

* Alibert, Y., & Benz, W. 2017, A&A, 598, L5

* Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357

* Arnold, J. A., Weinberger, A. J., Videen, G., & Zubko, E. S. 2022, ApJ, 930, 123

* Augereau, J. C., & Beust, H. 2006, A&A, 455, 987

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

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

* Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

* Barclay, T., Quintana, E. V., Raymond, S. N., & Penny, M. T. 2017, ApJ, 841, 86

* Barrado y Navascués, D., Stauffer, J. R., Song, I., & Caillault, J. P. 1999, ApJL, 520, L123

* Beichman, C. A., Krist, J., Trauger, J. T., et al. 2010, PASP, 122, 162

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

* Betti, S. K., Follette, K., Jorquera, S., et al. 2022, AJ, 163, 145

* Binks, A. S., & Jeffries, R. D. 2014, MNRAS, 438, L11

* Bitsch, B., Raymond, S. N., & Izidoro, A. 2019, A&A, 624, A109

* Boccaletti, A., Sezestre, E., Lagrange, A. M., et al. 2018, A&A, 614, A52

* Boccaletti, A., Thalmann, C., Lagrange, A.-M., et al. 2015, Natur, 526, 230

* Bowler, B. P., Liu, M. C., Shkolnik, E. L., & Tamura, M. 2015, ApJS, 216, 7

* Brande, J., Barclay, T., Schlieder, J. E., Lopez, E. D., & Quintana, E. V. 2020, AJ, 159, 18

* Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2022, spacetelescope/jwst: JWST v1.6.2, Zenodo, doi:10.5281/zenodo.6984366

* Cale, B. L., Reefe, M., Plavchan, P., et al. 2021, AJ, 162, 295

* Carter, A. L., Hinkley, S., Kammerer, J., et al. 2023, ApJL, 951, L20

* Caswell, T. A., Droettboom, M., Lee, A., et al. 2021, matplotlib/matplotlib: REL: v3.5.1, Zenodo, doi:10.5281/zenodo.5773480

* Chen, C., Ballering, N., Duchene, G., et al. 2019, BAAS, 51, 342

* Childs, A. C., Quintana, E., Barclay, T., & Steffen, J. H. 2019, MNRAS, 485, 541

* Clement, M. S., Quintana, E. V., & Quarles, B. L. 2022, ApJ, 928, 91

* Couture, D., Gagné, J., & Doyon, R. 2023, ApJ, 946, 6

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

* Daley, C., Hughes, A. M., Carter, E. S., et al. 2019, ApJ, 875, 87

* Debes, J. H., Jang-Condell, H., Weinberger, A. J., Roberge, A., & Schneider, G. 2013, ApJ, 771, 45

* den Brok, M., Peletier, R. F., Valentijn, E. A., et al. 2011, MNRAS, 414, 3052

* Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45

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

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

* Gáspár, A., Wolff, S. G., Rieke, G. H., et al. 2023, NatAs, 7, 790

* Gaudi, B. S., Meyer, M., & Christiansen, J. 2021, in ExoFrontiers: Big Questions in Exoplanetary Science, ed. N. Madhusudhan (Bristol: IOP Publishing), 2

* Gilbert, E. A., Barclay, T., Quintana, E. V., et al. 2022, AJ, 163, 147

* Girard, J. H., Leisenring, J., Kammerer, J., et al. 2022, Proc. SPIE, 12180, 121803Q

* Gomez Gonzalez, C. A., Wertz, O., Absil, O., et al. 2017, AJ, 154, 7

* Grady, C. A., Wisniewski, J. P., Schneider, G., et al. 2020, ApJL, 889, L21

* Graham, J. R., Kalas, P. G., & Matthews, B. C. 2007, ApJ, 654, 595

* Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1

* Hardegree-Ullman, K. K., Cushing, M. C., Muirhead, P. S., & Christiansen, J. L. 2019, AJ, 158, 75

* Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

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

* Hinkley, S., Biller, B., Skemer, A., et al. 2023, arXiv:2301.07199

* Honda, M., Inoue, A. K., Fukagawa, M., et al. 2009, ApJL, 690, L110

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

* Kalas, P., Liu, M. C., & Matthews, B. C. 2004, Sci, 303, 1990

* Kammerer, J., Girard, J., Carter, A. L., et al. 2022, Proc. SPIE, 12180, 121803N

* Keenan, P. C., & McNeil, R. C. 1989, ApJS, 71, 245

* Kim, M., Wolf, S., Potapov, A., Mutschke, H., & Jäger, C. 2019, A&A, 629, A141

* Klein, B., Donati, J.-F., Moutou, C., et al. 2021, MNRAS, 502, 188

* Krist, J. E., Ardila, D. R., Golimowski, D. A., et al. 2005, AJ, 129, 1008

* Krist, J. E., Balasubramanian, K., Beichman, C. A., et al. 2009, Proc. SPIE, 7440, 74400W

* Kunkel, W. E. 1973, ApJS, 25, 1

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

* Lannier, J., Delorme, P., Lagrange, A. M., et al. 2016, A&A, 596, A83

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

* Lawson, K., Currie, T., Wisniewski, J. P., et al. 2022, ApJL, 935, L25

* Leisenring, J. 2021, WebbPSF Extensions, v1.0.4., https://github.com/JarronL/webbpsf_ext

* Linder, E. F., Mordasini, C., Mollière, P., et al. 2019, A&A, 623, A85

* Liu, M. C. 2004, Sci, 305, 1442

* Liu, M. C., Matthews, B. C., Williams, J. P., & Kalas, P. G. 2004, ApJ, 608, 526

* Lomax, J. R., Wisniewski, J. P., Roberge, A., et al. 2018, AJ, 155, 62

* Lucy, L. B. 1974, AJ, 79, 745

* MacGregor, M. A., Wilner, D. J., Rosenfeld, K. A., et al. 2013, ApJL, 762, L21

* Martioli, E., Hébrard, G., Correia, A. C. M., Laskar, J., & Lecavelier des Etangs, A. 2021, A&A, 649, A177

* Mathioudakis, M., & Doyle, J. G. 1991, A&A, 244, 433

* Mathioudakis, M., & Doyle, J. G. 1993, A&A, 280, 181

* Matthews, B. C., Kennedy, G., Sibthorpe, B., et al. 2015, ApJ, 811, 100

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

* Meyer, M. R., Amara, A., Reggiani, M., & Quanz, S. P. 2018, A&A, 612, L3

* Miret-Roig, N., Galli, P. A. B., Brandner, W., et al. 2020, A&A, 642, A179

* Mullan, D. J., Stencel, R. E., & Backman, D. E. 1989, ApJ, 343, 400
Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2022, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python, v1.1.0, Zenodo, doi:10.5281/zenodo.7370358

Nielsen, E. L., De Rosa, R. J., Wang, J., et al. 2016, AJ, 152, 175

Okuta, R., Unno, Y., Nishino, D., Hido, S., & Crissman 2017, in 31st Conf. on Neural Information Processing Systems (New York: ACM)

Pearce, T. D., Launhardt, R., Ostermann, R., et al. 2022, A&A, 659, A135

Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, Proc. SPIE, 9143, 1174

Plavchan, P., Barclay, T., Gagné, J., et al. 2020, Natur, 582, 497

Richardson, W. H. 1972, JOSA, 62, 55

Rieke, M. J., Kelly, D., & Horner, S. 2005, Proc. SPIE, 5904, 1

Rieke, M. J., Kelly, D. M., Misselt, K., et al. 2023, PASP, 135, 028001

Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001

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

Schlecker, M., Pham, D., Burn, R., et al. 2021, A&A, 656, A73

Schlieder, J. E., Beichman, C. A., Meyer, M. R., & Greene, T. 2015, arXiv:1512.00053

Schlieder, J. E., Lépine, S., & Simon, M. 2010, AJ, 140, 119

Schlieder, J. E., Lépine, S., & Simon, M. 2012, AJ, 144, 109

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

Schüppler, C., Löhne, T., Krivov, A. V., et al. 2015, A&A, 581, A97

Shkolnik, E. L., Allers, K. N., Kraus, A. L., Liu, M. C., & Flagg, L. 2017, AJ, 154, 69

Shvartzvald, Y., Maoz, D., Udalski, A., et al. 2016, MNRAS, 457, 4089

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Song, I., Weinberger, A. J., Becklin, E. E., Zuckerman, B., & Chen, C. 2002, AJ, 124, 514

Stolker, T., Quanz, S. P., Todorov, K. O., et al. 2020, A&A, 635, A182

Storn, R., & Price, K. 1997, J. Glob. Optim., 11, 341

Suzuki, D., Bennett, D. P., Sumi, T., et al. 2016, ApJ, 833, 145

Sezestre, É., Augereau, J. C., Boccaletti, A., & Thébault, P. 2017, A&A, 607, A65

Tazaki, R., Murakawa, K., Muto, T., Honda, M., & Inoue, A. K. 2021, ApJ, 921, 173

Torres, C. A. O., & Ferraz Mello, S. 1973, A&A, 27, 231

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

Vizgan, D., Hughes, A. M., Carter, E. S., et al. 2022, ApJ, 935, 131

Vyssotsky, A. N. 1956, AJ, 61, 201

Wilson, O., & Woolley, R. 1970, MNRAS, 148, 463

Wisniewski, J. P., Kowalski, A. F., Davenport, J. R. A., et al. 2019, ApJL, 883, L8

Wittrock, J. M., Dreizler, S., Reefe, M. A., et al. 2022, AJ, 164, 27

Wittrock, J. M., Plavchan, P., Cale, B. L., et al. 2023, arXiv:2302.04922

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

Zicher, N., Barragán, O., Klein, B., et al. 2022, MNRAS, 512, 3060

Zuckerman, B., Song, I., Bessell, M. S., & Webb, R. A. 2001, ApJL, 562, L87