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

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. 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 ≲10 au.

AU Microscopii (AU Mic, HD 197481, GJ 803) is a very nearby (d = 9.714 ± 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 β Pictoris (β Pic), elevated HR diagram position, and strong activity. These diagnostics indicated an age of 20 ± 10 Myr. Further follow-up revealed that AU Mic and many other nearby, young M dwarfs, are members of the β 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 ≈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 ± 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 μ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 μm and refined the model of the spectral energy distribution (SED). At 850 μm, the system was detected as a point source—indicating the presence of an unresolved population of cold (∼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 ∼50 au to separations of ∼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 ∼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) ¹⁶ targets nine of the closest, youngest M dwarfs using NIRCam coronagraphy in two filters spanning 3–5 μ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 μm, and spanning separations from ∼03 to ∼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.

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) ¹⁷ of 063.

Observations were conducted using two filters from the NIRCam long wavelength (LW) channel (average pixel scale of 63 mas pixel⁻¹): F356W (λ_pivot = 3.563 μm, Δλ = 0.787 μm) and F444W (λ_pivot = 4.421 μm, Δλ = 1.024 μm) (Rieke et al. 2005). Using the SUB320 subarray mode, each integration in the LW channel produces an 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 (014) and 2.84 pixels (018) 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 ≳ 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.

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

3.2. Image Registration

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. ¹⁸ 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 (126), 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 (063–221; from roughly the IWA to where a background source begins to infringe for AU Mic), with the corresponding subtraction region covering r ≥ 20 pixels (126), 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.

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

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

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 μ_α = 281 mas yr⁻¹ and μ_δ = − 360 mas yr⁻¹; 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.

Zoom In Zoom Out Reset image size

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 12848 and an inclination of 90° for simplicity). For an array of projected separations spanning approximately 01–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 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 · FWHM) < r < (r_i + 0.5 · 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. ²⁰ 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 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_⊙) 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σ flux, mass, and contrast limits are presented in Figure 5.

Zoom In Zoom Out Reset image size

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 logarithmic grid of semimajor axes with 300 values from 1 to 1000 au. For each semimajor axis and mass, we generate a sample every 001 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.

Zoom In Zoom Out Reset image size

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_J planet at a projected separation of 271 along the southeastern extent of the disk, and a 0.2 M_J planet at a projected separation of 172 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 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_J and 0.2 M_J 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.

Zoom In Zoom Out Reset image size

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

Here, we provide measurements of surface brightness (SB) and color for AU Mic's disk from reductions using classical RDI, RDI/KLIP, 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 (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 12848 (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. ²¹

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 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{eqnarray}\begin{array}{rcl}{\rm{\Delta }}({m}_{356}-{m}_{444}) & = & -2.5\left[{\mathrm{log}}_{10}\displaystyle \frac{{F}_{356}^{d}}{{F}_{444}^{d}}-{\mathrm{log}}_{10}\displaystyle \frac{{F}_{356}^{* }}{{F}_{444}^{* }}\right]\\ & = & -2.5{\mathrm{log}}_{10}\displaystyle \frac{{F}_{356}^{d}{F}_{444}^{* }}{{F}_{444}^{d}{F}_{356}^{* }}.\end{array}\end{eqnarray} \tag{ 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.

Zoom In Zoom Out Reset image size

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 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 ≠ 90°). 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 ∼ 0063 at F606W), Krist et al. (2005) measure a projected disk FWHM as small as ∼022 (occurring at a projected separation of ∼15), while the FWHMs of the F356W and F444W filters are 014 and 018, 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 (0189) and height 1.5 pixels (0095)—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₀/r₀ = 10⁻⁵, α_in = 10, α_out = −10, γ = 2; otherwise as given in Table 2) ²⁴ and convolved with (a) coronagraphic images from WebbPSF_ext and (b) an empirical field PSF ²⁵ 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 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 ∼018 in vertical profiles measured using the aforementioned aperture.

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

Zoom In Zoom Out Reset image size

7.1. Discussion of Vertical Structure

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 ∼03 (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σ 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 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 northwestern 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 042 at 1'' in F356W, compared to 025–028 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.

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), SpaceKLIP (Kammerer et al. 2022), Vortex Image Processing (Gomez Gonzalez et al. 2017).

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 (β = 1) and a Gaussian vertical density distribution (γ = 2). The varied parameters include:

• 1.  
  inclination (incl, in degrees);
• 2.  
  position angle (PA, in degrees);
• 3.  
  fiducial radius (r₀, in au);
• 4.  
  the ratio of scale height to radius (h₀/r₀);
• 5.  
  the density power-law indices interior and exterior to the fiducial radius (α_in and α_out, respectively);
• 6.  
  the H–G asymmetry parameters, g₁ and g₂, and the weight for the first SPF (w₁, with w₂ = 1 − w₁).

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.

Zoom In Zoom Out Reset image size

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, ²⁶ 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. ²⁷ 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 =12848; 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 12848. 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 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.

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.

Zoom In Zoom Out Reset image size

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.