JWST/NIRCam Detection of the Fomalhaut C Debris Disk in Scattered Light

Observations of debris disks offer important insights into the formation and evolution of planetary systems. Though M dwarfs make up approximately 80% of nearby stars, very few M dwarf debris disks have been studied in detail—making it unclear how or if the information gleaned from studying debris disks around more massive stars extends to the more abundant M dwarf systems. We report the first scattered-light detection of the debris disk around the M4 star Fomalhaut C using JWST's Near Infrared Camera (NIRCam; 3.6 and 4.4 μm). This result adds to the prior sample of only four M dwarf debris disks with detections in scattered light and marks the latest spectral type and oldest star among them. The size and orientation of the disk in these data are generally consistent with the prior Atacama Large Millimeter/submillimeter Array submillimeter detection. Though no companions are identified, these data provide strong constraints on their presence—with sensitivity sufficient to recover sub-Saturn mass objects in the vicinity of the disk. This result illustrates the unique capability of JWST to uncover elusive M dwarf debris disks in scattered light and lays the groundwork for deeper studies of such objects in the 2–5 μm regime.


INTRODUCTION
Studies of debris disks, the dust-dominated remnants of the planet formation process, offer a number of compelling insights into how planetary systems form and evolve.Analysis of debris disk morphology provides a window into the dynamical history of these systems -revealing past stellar interactions (e.g., Ardila et al. 2005) or uncovering the "sign-posts" of yet-unseen companions (e.g., Kalas et al. 2005;Kenyon & Bromley 2008).Observational signatures in disks -such as sharp edges and gaps, warps, kinks, and filaments, brightness asymmetries, dust clumps, and other deviations from a smooth radial profile -have all been used to suggest the existence and properties of potential planets (e.g., Hashimoto et al. 2012;Gáspár et al. 2023).Simulations of these disk signatures have shown that they might be explained by direct planetary interactions (e.g., Dong et al. 2012;Lee & Chiang 2016), indirect planetary interactions (e.g., Yelverton & Kennedy 2018), or other phenomena relevant to understanding planet formation processes (e.g., Dong et al. 2018).Indeed, such morphological signatures have successfully predicted the presence of previously-unseen planets (e.g., β Pictoris b; Lagrange et al. 2009).Analysis of debris disks' brightness as a function of wavelength via multiwavelength studies enables assessment of the composition and size distribution of the constituent dust (Rodigas et al. 2015;Ballering et al. 2016;Arnold et al. 2019).This provides the opportunity to study the composition of planetesimals left over from the formation of planets, while also helping us to assess the mechanisms driving the evolution of post-accretion circumstellar material (e.g., Meyer et al. 2007).
Though M dwarfs make up approximately 80% of nearby stars (Reylé et al. 2021), few M-dwarf debris disks have been identified, and fewer studied in detail (e.g., Avenhaus et al. 2012;Luppe et al. 2020;Ren et al. 2021;Cronin-Coltsmann et al. 2023).Of the 10 known M-dwarf debris disk systems within 100 pc, only five have had their disks definitively spatially resolved (Cronin-Coltsmann et al. 2023), and only four have been detected in scattered light.This is in contrast to the ∼ 46 resolved debris disks around A-type stars 1 .The small sample of M-dwarf debris disks available for study precludes assessment of how the debris disks of the lowmass majority of stars, and thus their planetary systems, differ from those of more massive stars.While numerous studies have proposed mechanisms that might lead to Mdwarf debris disks being intrinsically less common than those of higher-mass stars (e.g., Lestrade et al. 2011), recent results suggest that they are equally abundant, but are often simply too faint to detect (Luppe et al. 2020;Cronin-Coltsmann et al. 2023).Moreover, simply detecting an M-dwarf debris disk is not sufficient for accessing the wealth of information that debris disks can provide.Deep morphological studies require strong, well-resolved detections to meaningfully disentangle dynamical scenarios.For example, precise measurements of disk eccentricity can test the possibility of past stellar interactions (e.g., Cronin-Coltsmann et al. 2021), while precise measurement of spiral arm motion can distinguish between their likely drivers -either gravitational instability or a planetary perturber (e.g., Dong et al. 2018).Meanwhile, studies of debris disk composition and dust properties require detections at multiple wavelengths -spanning both scattered light and thermal emission -to break degeneracies (e.g., Rodigas et al.   1 Based on the catalog of resolved disks at circumstellardisks.org 2015; Ballering et al. 2016).In both cases, the difficulty of resolving M-dwarf debris disks introduces significant barriers for understanding their planetary systems.
JWST's Near Infrared Camera (NIRCam) is a stable space-based high-contrast imager operating at nearinfrared (NIR) wavelengths (0.6-5 µm, Rieke et al. 2005).Since light scattered by dust around intrinsically redder M dwarfs can result in higher surface brightness in the near-IR than at optical wavelengths (even for disks with intrinsically blue scattering; e.g., Ren et al. 2021), observations with NIRCam can ease the requisite sensitivity for detection of M-dwarf debris disks in scattered light.Additionally, NIRCam provides the capability to characterize debris disks in the 2−5 µm range (e.g., Lawson et al. 2023), where ground-based observations are limited by telluric absorption and background levels.NIRCam's ability to use broader filters than the common ground-based filters at these wavelengths may likewise provide an advantage for efforts to study these disks.With many important scattered-light spectral features occurring at these wavelengths (e.g., the 3 µm water-ice feature; Kim et al. 2019), detections here are particularly powerful for diagnosing debris disk composition.Combined with the utility of ALMA for providing thermal detections of these objects (Luppe et al. 2020;Cronin-Coltsmann et al. 2023), JWST detections in scattered-light may provide a compelling new avenue for understanding the planetary systems of M dwarfs.Lawson et al. (2023) demonstrated the efficacy of NIRCam coronagraphy for studying M dwarf debris disks at 3-5 µm in application to the debris disk of AU Mic.Herein, we build upon these results to present the first scattered-light detections of the Fomalhaut C debris disk using coronagraphic observations from JWST/NIRCam in F356W (3.6 µm) and F444W (4.4 µm).Based on these data, we provide analysis of the disk's scattered-light morphology, brightness, and color, and also conduct a deep search for companions in the data -reaching sensitivities sufficient to detect masses below that of Saturn.

OBSERVATIONS
Fomalhaut C was observed as part of a JWST NIR-Cam Guaranteed Time Observing (GTO) program.This program, GTO 1184 (PI J. Schlieder)2 , targeted 9 nearby, young M dwarfs with NIRCam coronagraphy to search for low-mass companions.A full description of the GTO 1184 survey and its results will be presented in a separate paper (Bogat et al. in prep).
All targets were observed at two roll angles separated by ∼10 • using two filters from the NIRCam long wavelength (LW) channel (average pixel scale of 63 mas/pixel): F356W (λ pivot = 3.57 µm, ∆λ = 0.77 µm) and F444W (λ pivot = 4.36 µm, ∆λ = 0.94 µm) (Rieke et al. 2005;Leisenring 2021).Observing at two roll angles not only permits the use of angular differential imaging (ADI; Marois et al. 2006) for companion searches, but also enables any identified companions or smallscale disk features (which are fixed to the sky frame) to be distinguished from residual speckle noise or diffraction artifacts (which are fixed to the detector frame).The MASK335R coronagraph mask, which has an inner working angle (IWA)3 of 0. ′′ 63, was used for observations in both filters (Krist et al. 2009).Using the SUB320 subarray mode, each integration results in an image of 320×320 pixels (20×20 ′′ ).
Instead of observing dedicated reference stars for each target, the survey permits reference star differential imaging (RDI) using a self-referencing strategywherein RDI is performed for each target using the other survey targets as a reference library.The reference library that we utilize for Fomalhaut C includes four GTO 1184 survey targets found to be free of (detectable) circumstellar flux: AP Col, 2MJ0944, G-7-34, and HIP 17695.These observations are summarized in Table 1.
Of the reference targets, AP Col and 2MJ0944 were observed very close in time to the observations of Fomalhaut C, with no intervening tilt events or wavefront corrections.The remaining reference targets, G-7-34 and HIP 17695, were observed approximately two months prior to Fomalhaut C.During this interval, one large and one moderate tilt event occurred (Lajoie et al. 2023).As such, these targets' images are unlikely to improve point spread function (PSF) subtraction at separations where speckle noise dominates, but may nevertheless help to suppress noise in the background-limited regime.

DATA REDUCTION
We reduce the data predominantly following the procedure of Lawson et al. (2023), but also adopting a number of more recent additions to the spaceKLIP package (Kammerer et al. 2022) 4 and other changes, as summarized hereafter.
Using the rampfit step from spaceKLIP, we process Stage 0 products (*uncal.fits) to Stage 1 (*rateints.fits).For this purpose, we skip dark current subtraction5 , adopt a border of 4 pixels as pseudo reference pixels, and use a jump detection threshold of 4. This step now also includes a 1 / f (frequency) noise correction procedure to eliminate instrumental striping that can otherwise manifest in the data (see description in Rebollido et al. 2024).
Using the imgprocess step from spaceKLIP, we then process the Stage 1 products to Stage 2 (*calints.fits),including identification of outlier pixels with the outlier detection step.Next, we identify additional bad pixels and correct them (using fix bad pixels and replace nans from spaceKLIP.imagetools).Following this, manual inspection revealed six additional pixels with seemingly spurious values throughout the data (likely unflagged static hot pixels).Each of these pixels was replaced with the median value within a 5×5 pixel box.
After correction of outlier pixels, we fit and subtract a background model from each science and reference exposure.In this model, the distribution of background flux in the coronagraphic data is simulated by multiplying a uniform image by the coronagraph transmission map (including the neutral density squares; generated using WebbPSF ext, Leisenring 2021) and then convolving the result with a grid of spatially varying PSFs from WebbPSF ext -following the general convolution procedure described in Lawson et al. (2023, Appendix A).Justification of the need for this step and the details of our implementation are provided in Appendix A.
Following background subtraction, image alignment is carried out as described in Lawson et al. (2023) -using a synthetic stellar PSF to determine the positions of the occulted star in each exposure.The results of this process (see Table 1) show that the Fomalhaut C target Note-A summary of the science and reference targets and their observations used in this study.All targets except Fomalhaut C ("FOMALHAUT-C") were used as reference targets.All exposures used the MEDIUM8 readout pattern with 10 groups per integration and with 8 integrations per roll for F356W and 16 integrations per roll for F444W -for total exposure times of 1676 seconds and 3562 seconds respectively.a The identifier used for each target in the GTO 1184 program.b ALLWISE photometry; Wright et al. (2010) c The measured offsets between the target star and coronagraph center in detector coordinates (see Section 3).d 2MASS J09445422-1220544 acquisition (TA) procedure incurred particularly large errors, resulting in significant misalignment between the star and the coronagraph.As none of the available reference targets had similar misalignment, this should be expected to result in more significant speckle residuals at small separations following PSF-subtraction -inhibiting the sensitivity of the data in the speckle-dominated regime.

REFERENCE STAR DIFFERENTIAL IMAGING
To remove the pattern of diffracted starlight in each image, we adopt three distinct strategies: classical RDI, model constrained RDI (MCRDI, Lawson et al. 2022Lawson et al. , 2023)), and a variation of synthetic RDI (SRDI, Greenbaum et al. 2023).In each case, a final image is created by derotating and median combining the starlightsubtracted integrations (*calints.fitsfiles) across both rolls.

Classical RDI
First, the classical RDI procedure is carried out for each filter by median combining the integrations of AP Col ("V-AP-COL"), which is the better spectral match of the two GTO 1184 targets observed nearcontemporaneously with Fomalhaut C (see Table 1).We then scale the brightness of the PSF model to minimize squared residuals with the data within a 20 pixel radius aperture centered on the star (chosen to exclude the disk signal as identified in Cronin-Coltsmann et al. 2021).We find no visual improvement in the subtrac-tion by manually tuning the determined scaling factor by small amounts in either direction.

Model Constrained RDI
In Model Constrained RDI (MCRDI), a synthetic image of the circumstellar scene is used to prevent systematic overestimation of the stellar PSF's brightness during construction of the stellar PSF model (Lawson et al. 2022(Lawson et al. , 2023;;Rebollido et al. 2024), which otherwise results in so-called "oversubtraction" (e.g., Pueyo 2016).To apply MCRDI, an initial unconstrained RDI reduction is carried out.The RDI residuals are then fit with disk models using standard forward modeling techniques to identify the best-fitting disk model image.The resulting disk model is then used to carry out an MCRDI reduction of the data (Lawson et al. 2022).In the MCRDI procedure, the optimal stellar PSF model is constructed by comparing the reference images to science images from which the best-fitting disk model has been subtracted (where, in our implementation, the PSF model is formed from a linear combination of reference images).The resulting stellar PSF model is then subtracted from the original science images.This technique provides final images which are effectively free of oversubtraction and enables more accurate disk photometry than application of model-based photometric throughput corrections (Lawson et al. 2022).Since the F356W data have better spatial resolution than the F444W data, we first model the disk at F356W and then restrict the F444W model geometry to match the geometry of the F356W result, while still varying the scattering phase function and overall brightness of the disk.Descriptions of the disk models, model convolution, the PSF subtraction algorithm settings, and handling of background sources are provided in Appendix B.

Synthetic RDI
Variation in the on-sky NIRCam PSF is predominantly a function of a) the alignment between the target star and the coronagraph mask, b) the evolving wavefront error, and c) the spectral type of the star (e.g., Girard et al. 2022).The WebbPSF tool (Perrin et al. 2014) provides the ability to generate synthetic NIRCam PSFs where these factors are accounted for by (respectively) tuning mask alignment, using empirical optical path difference (OPD) maps, and providing a source spectrum.In some scenarios, these synthetic PSFs could be used in place of or to supplement on-sky reference PSFs for carrying out RDI (see, e.g., Greenbaum et al. 2023).This is a particularly appealing prospect in this application, given the previously noted differences in the coronagraph alignment between the science and reference data.
Though providing the aforementioned options for tuning of synthetic PSF models, the resulting models from WebbPSF are nevertheless imperfect.This is both because of imperfect knowledge of the true state of these parameters (mask alignment, wavefront error, and target spectrum), and because of inaccuracies or simplifications in the underlying WebbPSF optical model.In the case of the Fomalhaut C data, performing RDI with a nominal WebbPSF model (using the nearest available OPD file to the observations, the measured mask offsets, and an approximate AMES-Cond stellar spectrum, Allard et al. 2001) leaves significant, bright PSF residuals well in excess of the residuals from classical RDI or MCRDI.
To improve the nominal synthetic PSF models, we proceed with a strategy of empirical model corrections -leveraging the available reference images to mitigate the recurring inaccuracies in the synthetic PSFs.This procedure is detailed in Appendix C. By introducing these corrections, some noise is added to the otherwise noiseless synthetic PSF models.However, contrasts show a net improvement of up to a factor of six -reaching sensitivities comparable to those of MCRDI.

Detection of the Fomalhaut C Debris Disk in Scattered Light
In both F356W and F444W, and using all three starlight-subtraction techniques, faint extended signal is detected in the vicinity of Fomalhaut C (Figures 1, 2).Though manifesting at marginal signal-to-noise per resolution element (SNRE ∼ 3 in F356W), the distribution of the flux is generally consistent with the orientation and size of the Fomalhaut C disk as seen in the ALMA detection reported in Cronin-Coltsmann et al. (2021) (indicated by the silver contours in Figure 1 and by the dashed ellipses in Figure 2).For F444W, additional flux is seen extending to wider separations; this is discussed further in Section 6.
No plausible companion candidates are identified in the data.Though a number of sources are visible (see Figure 1), they are all either a) likely background sources -based on proper motion versus prior imagery, F356W−F444W color, or morphology (e.g., the apparently extended source to the south east) -or b) consistent with residual speckle noise (i.e., inconsistent between the two rolls).Assessment of companion detection limits for these data is provided in Section 5.2.
We also provide deprojected radial surface brightness profiles based on the procedure introduced in Marshall et al. ( 2018) and adopted in Cronin-Coltsmann et al. ( 2021), assuming a disk position angle of −63 • and an inclination of 44 • .Measurements for both filters are made in annuli of radial width equal to the instrumental FWHM for F444W (0. ′′ 17).The deprojected noise level is determined by propagating uncertainties from a radial noise map, generated using the MCRDI reduction with the best-fit disk model subtracted.These profiles are presented in Figure 3. Here, the detection in F356W manifests with peak signal to noise ratio of ∼ 17.These profiles are discussed further in Section 6.
Adopting stellar flux estimates of 387 and 291 mJy for F356W and F444W respectively (see Appendix D), these profiles can be used to estimate the 3-5 µm disk color (where disk color is the surface brightness color corrected for the color of the incident starlight).For this purpose, we adopt the peaks of the raw (unconvolved) model curves from Figure 3 as the disk photometry, as these avoid biasing the measurement by the differences between the filters' PSFs.These values are 1.45 ± 0.05 µJy / arcsec 2 at 3.6 µm (F356W) and 1.61 ± 0.06 µJy / arcsec 2 at 4.4 µm (F444W)6 .Computing disk color as S Mie theory and the agglomerated debris particle (ADP) models from Arnold et al. (2022), each calculated for the F356W and F444W filters, shows that such a red color is difficult to reproduce for typical size distributions and compositions (e.g., astronomical silicates, water ice, tholins, and olivine).At these wavelengths, we find the reddest colors for either standard astronomical silicates or water ice.For Mie theory models, the observed color is reproduced only for an extremely steep grain size distribution (power law q ∼ −6), with a minimum grain size of a min ∼ 2.5 µm or a min ∼ 4 µm for compositions of silicates or water ice, respectively; for the ADP models, we identify no solutions reasonably reproducing the measured color.Since dust grains around such low-luminosity stars are not subject to radiation pressure blowout (Arnold et al. 2019(Arnold et al. , 2022)), very small minimum grain sizes are expected -with Cronin-Coltsmann et al. ( 2022) assuming a value 0.1 µm for Fomalhaut C. As such, a value of a min ∼ 2.5 µm appears implausible, suggesting instead some other cause for the measured color.This result is discussed further in Section 6.

Companion Detection Limits
Before generating contrast curves, we first subtract the best-fitting MCRDI disk model from both the MCRDI and SRDI images to mitigate the impact of the disk on the radial noise estimation.Contrast curves are calculated as a function of angular separation using the meas contrast routine from the PyKLIP package (Wang et al. 2015) for the disk subtracted images and using stellar flux estimates based on TA images (see Appendix D).We then correct these raw contrast curves by dividing them by the coronagraph transmission profile.
For the MCRDI reductions, we make no additional corrections for PSF subtraction algorithm throughput.While the brightness of any companions would be decreased slightly by the MCRDI procedure (by ≲ 10%; see Figure 5 of Lawson et al. 2023), any companion candidates identified with marginal significance could be directly incorporated into the circumstellar model for MCRDI to avoid this (e.g., Rebollido et al. 2024).As such, the companion itself can be assumed to have no effect on the 5σ detection limits for MCRDI.
In the case of SRDI, the only comparison made between the PSF model and the science data is when scaling the model brightness to match the data within a narrow annulus spanning 3 to 15 pixels from the star.Any hypothetical companions with flux falling within  In each subplot, the profiles for the data, convolved model, and raw model are indicated by the solid blue line, solid orange line, and dotted orange line, respectively.The dashed gray line shows the 5σ noise level for the data.While the data and convolved model are in good agreement for F356W, for F444W the data shows a deficit at and interior to the disk peak (likely a PSF-subtraction effect) and is significantly more radially extended (discussed in Section 6).
this annulus would lead to some oversubtraction and would thus have less than 100% algorithmic throughput.To quantify this, we carry out forward modeling of this effect on synthetic data sequences containing only a companion PSF (generated using WebbPSF).For both filters, we find that throughput is typically high, ∼ 95%, but can dip as low as ∼ 90% for companions that are coincident with brighter stellar speckles and inside the 3 < r < 15 pixel annulus.We adopt the median throughput over the sampled position angles at each separation to correct the SRDI contrast curves.
To map the resulting sensitivities to companion masses, we assume a system age of 440 Myr and adopt synthetic photometry from Linder et al. (2019, BEX-HELIOS grid) for companions of 2 M J or less and from Phillips et al. (2020, ATMO-CEQ grid) for higher masses.The final detection limits are presented for both filters in the left panel of Figure 4.
Using these detection limits, we also compute a companion detection probability map for the F444W MCRDI reduction (the most sensitive reduction to lowmass companions), following the procedure described in Lawson et al. (2023).For simplicity, we assume companions on circular orbits oriented in the plane of the disk, with position angles of −63 • and inclinations of 44 • .The resulting probability map is presented in the right panel of Figure 4.These results indicate that the data should have revealed any 0.3 M J companions beyond ∼ 10 au, and any 1 M J companions beyond ∼ 5 au.These limits are discussed further in Section 6.

DISCUSSION
Prior to these results, only four M-dwarf debris disks had been imaged in scattered light: AU Mic, TWA 7, TWA 25, and GSC 07396-00759 (Kalas et al. 2004;Choquet et al. 2016;Sissa et al. 2018).Fomalhaut C is exceptional among these in having the lowest stellar luminosity by a factor of ∼ 24 (L ∼ 0.005 L ⊙ ; Stassun et al. 2019) and the oldest age by a factor of ∼ 18 (440 Myr; Mamajek et al. 2013;Kalas et al. 2004;Sissa et al. 2018).As such, Fomalhaut C represents an entirely unprecedented class of scattered-light debris disk host.
The value of these data toward making a more detailed assessment of the disk's scattered-light morphology and composition is limited by the modest significance of the detection and by the presence of numerous background objects and significant residual speckle noise.We can, however, comment tentatively on a number of key parameters in the context of the prior Cronin-Coltsmann et al. ( 2021) analysis.Modeling of the ALMA detection in Cronin-Coltsmann et al. ( 2021) identified a particularly narrow best-fitting disk scale width of ∼0.′′ 11 with a 3σ upper limit of 0. ′′ 6.Meanwhile, the spatial resolution of NIRCam F356W data should be capable of distinguishing scale widths as small as approximately 0. ′′ 06 (half of the PSF FWHM in the narrowest direction).The best fitting disk model we identify during the MCRDI procedure has a fiducial radius of r 0 = 24.5 au and a radial power law index of α = 8.6.This corresponds to a radial scale width of 0. ′′ 54 (4.2 au)just within the reported 3σ upper limit from Cronin-Coltsmann et al. (2021).Estimating the uncertainty on the fit parameters via Hessian matrix inversion7 , yields a 3σ lower limit on scale width of 0. ′′ 46.Combined with the smaller disk radius that we identify -24.5 ± 0.1 au versus 26.4±0.6 au from Cronin-Coltsmann et al. ( 2021) -this could indicate that the scattered light detection corresponds to a population of dust that is not colocated with the material seen by ALMA.We remark, however, that these uncertainties neglect a number of less tractable sources of uncertainty (e.g., the possible presence of unmasked background sources) and thus likely overestimate the significance of these differences.A more confident assessment of this possibility will require higher significance detections.
Unlike the F356W data, the extended emission observed in F444W is not well-explained by a narrow ringlike disk alone.After subtraction of the best-fitting disk model, significant residuals are seen spanning much of the field of view (see Figure 2, lower right panel, and Figure 3, right panel).This could be explained in a few ways.First, we may be seeing residual background emission.While our background model assumes a uniform background convolved with the instrumental response, astrophysical variations in the background could produce the observed residuals.In this case, the F444W brightness of the disk derived from forward modeling may also be inflated, which would help to reconcile the exceptionally red color noted in Section 5.1.Alternatively, this could be evidence of an extended outer halo.Some weak evidence of such a halo is seen in analysis of Herschel 160 µm PACS data -in the form of an elevated background level (Cronin-Coltsmann et al. 2021).The apparent absence of such a feature in the F356W data could be reconciled if the halo appears more extended at 3.6 µm -perhaps due to an increasing presence of smaller grains at wider separations.In this scenario, the halo could appear with a sufficiently shallow slope over the NIRCam field of view so as to be effec- tively nulled during the background subtraction process.This background oversubtraction would also induce an apparent reddening of the main ring-like component, which would help to explain the seemingly-implausible red disk color that we measure.Finally, we note that the departure of the radial profile for the F444W data from that of the model interior to the peak (≲ 3. ′′ 5) is consistent with PSF subtraction residuals.Applying PSF subtraction to a simulated disk-free dataset (having the same coronagraph alignments as the real science and reference data) and then repeating these measurements shows a similar negative trend for the F444W at these separations.
The companion detection limits of Section 5.2 place strong constraints on the presence of planetary-mass companions around Fomalhaut C. For a mid-M target, the constraints afforded by JWST's sensitivity are unprecedented.While ground-based observations of targets with similar spectral types reach contrasts sufficient to detect masses ≳ 6 M J at 2 ′′ (e.g., Uyama et al. 2023), these JWST NIRCam data reach contrasts sufficient to uncover ∼ 0.3 M J planets at the same separation.These limits can also inform predictions specific to the Fomalhaut C system.Cronin-Coltsmann et al. (2021) note that if the observed radial extent of the disk is the result of material confined between the 3:2 and 2:1 orbital resonances of an unseen planet, then the planet should have a semimajor axis of ∼ 17 − 20 au.If such a planet is present in the system, our results effectively constrain its mass to ≲ 0.3 M J (approximately the mass of Saturn).

CONCLUSIONS
Based on coronagraphic imagery from JWST / NIR-Cam, we have presented the first scattered-light detections of the Fomalhaut C debris disk -marking the latest spectral type star with a debris disk detected in scattered light to date.We summarize our key findings hereafter.
1.In both F356W (3.6 µm) and F444W (4.4 µm), the orientation and size of the disk are largely consistent with the prior sub-mm detection with ALMA (Cronin-Coltsmann et al. 2021).In F444W, residual circumstellar flux not detected in F356W could suggest a faint extended halo or localized background variations.
2. By making empirical corrections to synthetic stellar PSFs, we improve contrasts for PSF subtraction with synthetic PSFs by factors of ∼ 6 in the speckle-limited regime and ∼ 2 in the backgroundlimited regime.Compared with other state-ofthe-art PSF subtraction methods, this approach achieves comparable contrasts in the backgroundlimited regime, and contrasts within a factor of ∼ 2 at smaller separations.
3. Modeling of the disk in scattered light favors a slightly smaller disk radius (24.5 au) and broader scale width (0. ′′ 54) than was identified from the sub-mm ALMA detection.This may indicate that the material probed in scattered light is not colocated with the material seen by ALMA.However, geometric constraints from these data are limited by the significance of the detection and the presence of coincident background sources -such that these modeling results should be considered tentative.
4. No companion candidates are identified; companion detection limits at 4.4 µm effectively rule out companions more massive than Saturn beyond approximately 10 au and more massive than Jupiter beyond approximately 5 au.
Deeper follow-up observations of Fomalhaut C with JWST/NIRCam would allow better assessment of the disk's scattered-light morphology while more confidently rejecting additional coincident background sources.Combined with improvements to the accuracy of NIR-Cam target acquisition centroiding since Cycle 1 (Girard et al. in prep), such observations should also enable a more meaningful probe of the environment interior to the ring as well.This would not only provide clues regarding any yet-unseen companions, but would also better inform the dynamical origins of the Fomalhaut triple system through precise measurement of Fomalhaut C's eccentricity (Shannon et al. 2014;Kaib et al. 2018;Feng & Jones 2018;Cronin-Coltsmann et al. 2021).
Overall, these results highlight a key utility for JWST: its ability to study debris disks around the lowest mass stars in scattered light.Future studies leveraging this capability -along with the capability of ALMA to study these elusive disks in thermal emission -stand to make transformational contributions to our understanding of planetary systems around the most common stars.

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 (among others) enabled this study through contributions to either the 2002 NIRCam instrument proposal, the development and commissioning of the NIRCam instrument, or the commissioning of the NIRCam coronagraphy mode: Martha Boyer, Alicia Canipe, Eiichi Egami, Daniel Eisenstein, Bryan Hilbert, Klaus Hodapp, Scott Horner, Doug Kelly, John Krist, Don McCarthy, Karl Misselt, George Rieke, John Stansberry, and Erick Young.We are grateful for support from NASA through the JWST NIRCam project, contract number NAS5-02105 (M.Rieke, University of Arizona, PI).
The authors thank G. Kennedy for sharing the reduced ALMA data for use in reproducing the Cronin-Coltsmann et al. ( 2021) contours in our Figure 1.
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 https://doi.org/10.17909/he30-tx49.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.
This research made use of POPPY, an open-source optical propagation Python package originally developed for the James Webb Space Telescope project (Perrin, 2012).
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: Astropy (Astropy Collaboration et al. 2013Collaboration et al. , 2018Collaboration et al. , 2022)), CuPy (Okuta et al. 2017), LMFIT (Newville et al. 2022), Matplotlib (Hunter 2007;Caswell et al. 2021), NumPy (Harris et al. 2020), POPPY (Perrin et al. 2012), pyKLIP (Wang et al. 2015), SciPy (Virtanen et al. 2020), SpaceKLIP (Kammerer et al. 2022), Vortex Image Processing (Gomez Gonzalez et al. 2017), WebbPSF (Perrin et al. 2014), WebbPSF ext (Leisenring 2021) APPENDIX A. BACKGROUND MODELS Treatment of the background in these data requires particular care for a number of reasons.First, the lack of a dedicated reference star is likely to lead to differences in background level between the science and reference images.Typically, a science and dedicated reference target will be near one another on-sky and will be observed in a non-interruptible sequence -leading to comparable backgrounds levels that are ultimately eliminated during PSF subtraction (so long as the stars are also of comparable brightness).Since the targets used as references here were much further on-sky than is typical and sometimes had sizeable offsets in time from the science observations, significantly distinct background levels could manifest.For the observations of Fomalhaut C, the jwst backgrounds tool (Rigby et al. 2023) predicts backgrounds of 4.4 µJy / arcsec 2 for F356W and 12.4 µJy / arcsec 2 for F444W.For the utilized reference target observations, jwst backgrounds predicts backgrounds of 2.3-3.4 µJy / arcsec 2 for F356W and 5.3-9.5 µJy / arcsec 2 for F444W.Second, the background will be affected by the coronagraph's transmission as well as by the neutral density squares in the subarray corners.As such, treating the background as a uniform value (e.g., subtracting a median value from the images) will effectively over-subtract the vicinity of the coronagraph and neutral density squares by differing amounts between images with differing background levels (while under-subtracting elsewhere).This will, in turn, degrade eventual PSF subtraction and likely introduce global photometric inaccuracy.Finally, while these effects may be negligible in scenarios where differences in background levels are substantially smaller than expected sensitivity, the background levels among these targets vary easily on the order of the expected sensitivity and thus cannot be safely ignored.
Fitting of the background model is carried out prior to alignment of the images, such that the brightness of the convolved background image described in Section 3 can simply be scaled in brightness to reproduce the observed background.In addition to this scaling factor, fitting of the background in these data considers two other com-ponents.First, the use of pseudo reference pixels during ramp fitting effectively subtracts a uniform background level from the data.Since the subtracted background level is not known, we simply optimize for a best-fit uniform offset alongside the scalar for the convolved background image (effectively correcting the incidental background subtraction that occurs during ramp fitting).To avoid any contribution from the wings of the stellar PSF, we also subtract a synthetic stellar PSF from the data during the fitting.Since the data have not yet been aligned, cross-correlation is used to align the synthetic PSF with each exposure being considered (as is typically done during image alignment).The overall brightness of the PSF model is then tuned with the other two parameters during the background optimization procedure.
The three parameter model is optimized for each individual exposure (whose integrations are mediancombined for this purpose), excluding the region within 60 pixels of the coronagraph center to avoid the portions of the stellar PSF that are poorly described by the cross-correlated PSF model (e.g., due to coronagraph misalignement).We additionally exclude a) a 5 pixel border around the image edge where rows of apparently spurious (significantly positive or negative) values sometimes manifest after ramp fitting, and b) the vicinity of two bright diffraction features that manifest at the subarray edges in the y-axis direction from the target star, but which do not appear in the WebbPSF models.To mitigate the effects of any faint background objects, uncorrected cosmic rays, or detector artifacts, the goodnessof-fit calculation excludes the 10% most outlying pixels above and below the median of the model-subtracted residuals for each set of trial parameters.
Once the optimization is completed for each exposure, we add the best-fit value for the uniform offset back to each individual integration, and then subtract the best-fit convolved background image to reach the final background-subtracted images.Figure 5 shows the results of this procedure for both a science and a reference exposure.The initial unconstrained RDI reduction for each filter considers two annular optimization regions spanning stellocentric separations of 4 ≤ r < 25 pixels and 25 ≤ r < 100 pixels, with corresponding subtraction regions of r < 25 pixels and r ≥ 25 pixels, respectively.The optimization regions additionally exclude the vicinity of four off-axis (likely background) sources in the science target's field of view.The radial boundaries for these regions were selected to approximately isolate the inner speckle-dominated regime.At larger separations, a total of 13 reference integrations featured significant uncorrected cosmic ray halos.By splitting the subtraction into two regions, we were able to exclude the affected frames for only the outer region, while allowing the inner region to benefit from a larger reference library.For each subtraction region, the PSF model is constructed as the linear combination of the reference images that minimizes squared residuals with the science image within the corresponding optimization region.
Following the initial RDI procedure, disk forward modeling is carried out.Raw disk models are generated using a simplified version of the GRaTer code (Augereau et al. 1999) implemented in Vortex Image Processing (VIP; Gomez Gonzalez et al. 2017).For the scattering phase function (SPF) of the disk, we adopt a simple Henyey-Greenstein (H-G) SPF (Henyey & Greenstein 1941).Though such an SPF is not expected to be physically informative (e.g., Hughes et al. 2018), the primary intent of the disk model here is to superficially reproduce the observed disk in order to suppress RDI oversubtraction.For this purpose, the adopted SPF is suitable.
The disk's inclination and position angle are fixed to the values from the overall best model reported in Cronin-Coltsmann et al. (2021): 44 • and −63 • respectively.The varied parameters for the F356W model include fiducial radius, scale height, density power-law index (α = α in = −α out ), H-G asymmetry parameter, eccentricity, and argument of periapsis.Additionally, an optimal overall brightness scaling factor is determined analytically for each model (i.e., this parameter is not varied by the optimizer).For the F444W model, we fix the parameters governing the disk's geometry to the best-fit values from the F356W procedure -tuning only the asymmetry parameter and the overall brightness scalar.
Convolution of disk models follows the procedure of Lawson et al. (2023), using a grid of synthetic PSF models and coronagraph transmission maps generated with consideration for the mask position in each roll after im-age alignment.The PSF grid samples the origin as well as 12 logarithmically-spaced radial positions at each of four linearly-spaced azimuthal positions (for a total of 49 spatial samples).Each model is then forward-modeled as typical for RDI, evaluating goodness of fit with a simple χ 2 metric within the outer optimization region described above (to mitigate the impact of residual speckle noise in the inner region).
After the best-fit model is identified, that model is used to carry out an MCRDI reduction of the data.In this process, the stellar PSF model is constructed by comparing the reference images with science images from which the optimal disk model has been subtracted.The resulting stellar PSF model is then subtracted from the original science frames.Besides the disk model "constraint", the PSF subtraction for the MCRDI reduction is identical to the initial unconstrained RDI procedure.

C. EMPIRICAL CORRECTIONS FOR SYNTHETIC RDI
Empirical corrections to the nominal synthetic reference PSF model for the science data are determined as follows.For each reference exposure, we generate a nominal PSF model image using WebbPSF as we did for Fomalhaut C. The brightness of each PSF model is scaled to minimize squared residuals with the data within an annular region spanning 3 ≤ r < 15 pixels from the star.The model is then subtracted from the data and the residuals are divided by the total model flux to produce a normalized residual map for each reference exposure.
Inspection of the residual maps shows that the model inaccuracy changes little between exposures.This suggests that the bulk of the PSF model inaccuracy results from characteristics that do not change across the sample of reference exposures -such as underlying simplifications in the WebbPSF optical model.The most significant changes that do occur in the residuals from exposure to exposure appear to be temporally correlated; for example, the residual maps for AP Col and 2MJ0944 (November 2022) are more similar to one-another than to those of G-7-34 and HIP 17695 (October 2022).Temporally correlated errors could occur as a result of differences in OPD file accuracy.As the data are distributed roughly into two groups of observations -October 2022 and November 2022 -such a correlation could manifest if a tilt event occurred between the time of one set of observations and the utilized OPD file, resulting in a less accurate OPD measurement for those exposures.However, comparison of the OPD measurements flanking each observation show no evidence of such an event -with ∆WFE = 6.4 nm RMS between the preceding and succeeding OPD files for both the October and November data -below the Cycle 1 median of 9 nm RMS (Lajoie et al. 2023).
Since the GTO 1184 program was observed roughly in order of spectral type, the underlying driver of this phenomenon could also be related to the spectra of the targets.For example, synthetic PSF inaccuracies correlated with spectral type could manifest if a) the utilized synthetic spectra are inaccurate as a function of spectral type, or b) the inaccuracies in the WebbPSF models have wavelength dependence.Ultimately, these data are insufficient to confidently disentangle these possibilities.Conducting this analysis on a broader and more diverse set of observations is necessary to better understand the relevant factors.
To determine the model correction map for Fomalhaut C, we take the exposure-wise median of the normalized residual maps for exposures of both AP Col and 2MJ0944.These targets are selected because they were observed closest in time to the science target, and should thus yield the most accurate correction based on the analysis above.To produce the final corrected PSF model for each science exposure, we multiply the normalized correction map by the total science model flux and add this to the nominal PSF model.In Figure 6, we compare the resulting contrasts between the nominal and corrected synthetic RDI reductions and provide example images for one roll of the F356W data.The apparently worse performance for SRDI in F444W than in F356W (see Figure 2) could be evidence that the nominal PSF model deficiencies are indeed spectrally dependent.
The cost of this procedure is the introduction of noise to the otherwise noiseless synthetic PSF model.For this particular application, the noise introduced this way is significantly outweighed by the improved suppression of the diffraction pattern at all separations.In broader applications, it may be beneficial to adopt model corrections only at separations where this balance is favorable.

D. NIRCAM STELLAR FLUX ESTIMATION
To enable assessment of the sensitivity of these data in terms of contrast, we approximate the flux of Fomalhaut C in each filter as follows.We begin with the unocculted target acquisition images (one for each roll), which were taken in the F335M filter (3.35 µm, no neutral density square), and convert them from units of MJy/sr to mJy/pixel 2 .We measure initial fluxes for each image as the sum of pixels within a 15 pixel radius aperture centered on the PSF peak.To make aperture corrections to these fluxes, we use WebbPSF to compute the encircled energy at 15 pixels for a synthetic PSF normalized such that an infinite aperture would measure a

Figure 1 .
Figure 1.Left: NIRCam F356W detection of the Fomalhaut C debris disk in scattered light resulting from the MCRDI procedure (Section 4).The image has been smoothed by a gaussian kernel with FWHM of 6.7 pixels (3× the PSF FWHM).The inner software mask has a radius of 1. ′′ 5 and covers the region of significant residual speckle noise.Logarithmically-spaced black contours are drawn above the maximum for the color stretch to show the peak location of the bright south-western background source.Right: As the image on the left, but with silver contours showing the 1-5 σ levels for the ALMA detection of the disk reported in Cronin-Coltsmann et al. (2021).The gold arrows extend from the locations of the two assumed background sources modeled in Cronin-Coltsmann et al. (2021), showing their expected displacement due to the proper motion of Fomalhaut C (Gaia Collaboration et al. 2022).Both predicted positions are consistent with sources detected in the 2022 NIRCam data, supporting their identity as unassociated background objects.

Figure 2 .
Figure 2. PSF-subtracted final images of Fomalhaut C in F356W (top row) and F444W (bottom row), as well as the MCRDI disk model and disk model residuals (two rightmost columns).To reduce noise, each image has been median binned by a factor of 5 (such that each pixel in the binned image is the median of a 5 × 5 pixel square in the original image).The first three columns correspond to the three PSF subtraction strategies (Section 4): classical RDI, synthetic RDI, and model constrained RDI (MCRDI), respectively.The dashed ellipses mark the approximate shape and size of the disk identified by Cronin-Coltsmann et al. (2021).

Figure 3 .
Figure 3. Deprojected radial surface brightness profiles for MCRDI reductions of the F356W (left) and F444W (right) data.In each subplot, the profiles for the data, convolved model, and raw model are indicated by the solid blue line, solid orange line, and dotted orange line, respectively.The dashed gray line shows the 5σ noise level for the data.While the data and convolved model are in good agreement for F356W, for F444W the data shows a deficit at and interior to the disk peak (likely a PSF-subtraction effect) and is significantly more radially extended (discussed in Section 6).

Figure 4 .
Figure 4. Left: Radial 5σ companion detection limits for model constrained RDI (MCRDI) and synthetic RDI (SRDI)reductions as a function of angular separation in arcsec (lower x-axis) and projected separation in au (upper x-axis).Limits are provided in terms of contrast with the host star (upper left) and companion mass (lower left).Right: Companion detection probabilities for the F444W (4.4 µm) MCRDI reduction as a function of semimajor axis (for circular orbits coplanar with the disk) and companion mass (Jupiter-masses on the left axis, and Earth-masses on the right axis).White contours mark the 68%, 95%, and 99.7% probability thresholds.For context, the best-fit disk radius of 24.5 au is indicated by a vertical dashed line.These data effectively rule out companions more massive than ∼ 0.3 MJ in the vicinity of the disk.

Figure 5 .
Figure 5.A comparison of the initial and corrected images from the background model subtraction procedure (Appendix A) for one science exposure (top row) and one reference exposure (bottom row) in F356W.The images of each row are normalized to the initial stellar PSF brightness to highlight the difference between the backgrounds (visible as initially more negative regions within the neutral density squares at the corners).The right-most column shows the background model that is subtracted from the initial image..

Figure 6 .
Figure 6.Left: Contrast curves for the nominal (dotted lines) and corrected (dashed lines) synthetic RDI reductions of the Fomalhaut C data.The implemented corrections improve contrasts by up to a factor of six for both filters.Right: A comparison of synthetic PSF models for one roll of Fomalhaut C F356W data without (top row) and with (bottom row) the empirical corrections described in Appendix C. All images were divided by the radial average for the data to better compare the speckle patterns at all separations.The field of view for each panel is 2 × 2 ′′ ..