# The JWST Early Release Science Program for Direct Observations of Exoplanetary Systems. VI. Evidence for Radially Evolving Icy Grains in the HD 141569A Disk via NIRCam Coronagraphic Imaging

Maxwell A Millar-Blanchaer, Élodie Choquet, Kellen Lawson, Sebastián Marino, Jens Kammerer, Aarynn L Carter, Isabel Rebollido, Jarron M Leisenring, Minjae Kim, Paul Kalas, et al.

► To cite this version:

Maxwell A Millar-Blanchaer, Élodie Choquet, Kellen Lawson, Sebastián Marino, Jens Kammerer, et al. The JWST Early Release Science Program for Direct Observations of Exoplanetary Systems. VI. Evidence for Radially Evolving Icy Grains in the HD 141569A Disk via NIRCam Coronagraphic Imaging. The Astrophysical Journal, 2025, 994 (2), pp.199. 10.3847/1538-4357/ae0615. hal-05486944

HAL Id: hal-05486944

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

Submitted on 31 Jan 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# The JWST Early Release Science Program for Direct Observations of Exoplanetary Systems. VI. Evidence for Radially Evolving Icy Grains in the HD 141569A Disk via NIRCam Coronagraphic Imaging

Maxwell A. Millar-Blanchaer¹, Élodie Choquet², Kellen Lawson³, Sebastián Marino⁴, Jens Kammerer⁵,⁶, Aarynn L. Carter⁷, Isabel Rebollido⁸, Jarron M. Leisenring⁹, Minjae Kim¹⁰,¹¹, Paul Kalas¹²,¹³,¹⁴, Karl R Stapelfeldt¹⁵, Sasha Hinkley¹⁶, Mark Booth¹⁷,¹⁸, Carol A. Grady¹⁹, Elisabeth C. Matthews²⁰, Beth A. Biller²¹,²², Andrew Skemer²³, Julien H. Girard⁷, Schuyler G. Wolff²⁴, Kimberly Ward-Duong²⁵, Michael R. Meyer²⁶, Anthony Boccaletti²⁷, Eric Pantin²⁷, Brenda C. Matthews²⁸, Stanimir Metchev²⁹, Marshall D. Perrin⁷, Christine H. Chen⁷,³⁰, Katie Crotts⁷, Olivier Absil³¹, William O. Balmer⁶,³⁰, Per Calissendorff³², Gabriele Cugno³³, Thayne Currie³⁴,³⁵, Camilla Daniëlski³⁶, Kielan K. W. Hoch⁷, Markus Janson³⁷, Elena Manjavacas³⁸, Pierre-Olivier Lagage³⁹, Ben J. Sutlieff²¹,²², Shrishmoy Ray⁴⁰, Bin B. Ren⁴¹, Emily Rickman⁴², Genaro Suárez⁴³, Christopher A. Theissen⁴⁴, Taichi Uyama⁴⁵, Andreas Quirrenbach⁴⁶, Jason J. Wang⁴⁷,⁴⁸, Niall Whiteford⁴⁹, Mark C. Wyatt⁵⁰, and Alice Zurlo⁵¹,⁵²

(The JWST ERS Collaboration)

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

²Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

³NASA-Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA

⁴Department of Physics and Astronomy, University of Exeter, Stocker Road, Exeter EX4 4QL, UK

⁵European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748, Garching, Germany

⁶Space Telescope Science Institute, Baltimore, MD 21218, USA

⁷Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

⁸European Space Agency (ESA), European Space Astronomy Centre (ESAC), Camino Bajo del Castillo s/n, 28692 Villanueva de la Cañada, Madrid, Spain

⁹Steward Observatory, University of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA

¹⁰Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK

¹¹Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV 4 7AL, UK

¹²Astronomy Department, University of California, Berkeley, CA 94720, USA

¹³SETI Institute, Carl Sagan Center, 189 Bernardo Avenue, Mountain View, CA 94043, USA

¹⁴Institute of Astrophysics, FORTH, GR-71110 Heraklion, Greece

¹⁵Jet Propulsion Laboratory, California Institute of Technology, Mail Stop 321-100, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

¹⁶University of Exeter, Astrophysics Group, Physics Building, Stocker Road, Exeter, EX4 4QL, UK

¹⁷UK Astronomy Technology Centre, Royal Observatory Edinburgh, Blackford Hill, Edinburgh, EH9 3HJ, UK

¹⁸Astrophysikalisches Institut und Universitätsternwarte, Friedrich-Schiller-Universität Jena, Schillergäßchen 2-3, D-07745 Jena, Germany

¹⁹Eureka Scientific, 2452 Delmer Street, Suite 1, Oakland, CA 96402, USA

²⁰Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

²¹Scottish Universities Physics Alliance, Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, UK

²²Centre for Exoplanet Astronomy, University of Edinburgh, EH9 3HJ, UK

²³University of California Santa Cruz, Santa Cruz, CA, USA

²⁴Steward Observatory and the Department of Astronomy, The University of Arizona, 933 N Cherry Avenue, Tucson, AZ 85719, USA

²⁵Department of Astronomy, Smith College, Northampton, MA 01063, USA

²⁶Department of Astronomy, University of Michigan, 1085 S. University, Ann Arbor, MI 48109, USA

²⁷LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Jules Janssen, 92195 Meudon, France

²⁸Herzberg Astronomy & Astrophysics Research Centre, National Research Council of Canada, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada

²⁹Western University, Department of Physics & Astronomy and Institute for Earth and Space Exploration, 1151 Richmond Street, London, Ontario N6A 3K7, Canada

³⁰Department of Physics and Astronomy, The Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA

³¹STAR Institute, Université de Liège, Allée du Six Août 19c, 4000 Liège, Belgium

³²Department of Astronomy, University of Michigan, 1085 S. University, Ann Arbor, MI 48103, USA

³³Department of Astrophysics, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland

³⁴Department of Physics and Astronomy, University of Texas-San Antonio, 1 UTSA Circle, San Antonio, TX, USA

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

³⁶INAF-Oservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy

³⁷Department of Astronomy, Stockholm University, AlbaNova University Center, SE-10691 Stockholm, Sweden

³⁸AURA for the European Space Agency (ESA), ESA Office, Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

³⁹Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France

⁴⁰School of Mathematics and Physics, University of Queensland, St Lucia, QLD 4072, Australia

⁴¹Université Côte d'Azur, Observatoire de la Côte d'Azur (OCA), CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice cedex 4, France

⁴²European Space Agency (ESA), ESA Office, Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

⁴³Department of Astrophysics, American Museum of Natural History, Central Park West at 79th Street, NY 10024, USA

⁴⁴Department of Astronomy & Astrophysics, University of California, San Diego, La Jolla, CA 92093, USA

⁴⁵Department of Physics and Astronomy, California State University Northridge, 18111 Nordhoff Street, Northridge, CA 91330, USA

⁴⁶Landesfernwarte, Zentrum für Astronomie der Universität, Heidelberg, Königstuhl 12, D-69117 Heidelberg, Germany

⁴⁷Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA

⁴⁸Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA

⁴⁹American Museum of Natural History Department of Astrophysics, Central Park West at 79th Street, New York, NY 10034, USA

⁵⁰Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK

⁵¹Instituto de Estudios Astrofísicos, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército Libertador 441, Santiago, Chile
⁵² Millennium Nucleus on Young Exoplanets and their Moons (YEMS), Chile

Received 2025 June 5; revised 2025 August 11; accepted 2025 August 12; published 2025 November 26

## Abstract

We present JWST NIRCam coronagraphic observations of the HD 141569A circumstellar disk, obtained as part of the JWST Early Release Science program. The observations recover the multi-ringed structure seen in previous shorter-wavelength observations, but at filters centered on the ~3 µm water ice absorption feature and a complementary continuum region (F300M and F360M, respectively). The observations reveal apparent absorption between the F300M and F360M filters that decreases with radius, with a notable change around 200 au, between the innermost and outermost two rings. These results are consistent whether the data is reduced via deconvolution or through a forward-modeling approach. We demonstrate that these changes suggest a radial decrease in the water ice mass fraction by a factor of ~3–10 and possibly a change in minimum grain size. We do not detect any point sources within the system and can exclude planetary companions 2 Jupiter masses and greater beyond 1" radius (~111 au). These observations and the subsequent analysis illustrate a robust pathway for using JWST/NIRCam to characterize the distribution of water ice in other circumstellar disks. We highlight some of the early lessons learned from this work that we hope will be useful for future circumstellar disk observation planning and analysis.

Unified Astronomy Thesaurus concepts: James Webb Space Telescope (2291); Coronagraphic imaging (313); Circumstellar disks (235); Direct imaging (387); High contrast techniques (2369); Infrared observatories (791); Space observatories (1543); Exoplanet astronomy (486)

## 1. Introduction

Protoplanetary and debris disks, collectively referred to as circumstellar disks, are important constituents of planetary systems, as they can directly inform us about ongoing or past planet formation processes (e.g., M. C. Wyatt et al. 2015; A. M. Hughes et al. 2018; M. Benisty et al. 2023). Scattered light imaging of circumstellar disks has largely been conducted to date at wavelengths within and shortward of the K-band atmospheric transmission window, where the low sky and telescope thermal backgrounds are most conducive to the detection of faint circumstellar nebulosity and where coronagraphic instruments have historically operated (e.g., E. A. Rich et al. 2022; A. Garufi et al. 2024; C. Ginski et al. 2024; P. G. Valegård et al. 2024). Longward of the K band, only a small handful of disks have been imaged in scattered light, and prior to JWST, all from the ground (e.g., C. McCabe et al. 2003; G. Duchêne et al. 2004, 2010; C. McCabe et al. 2011; J. Milli et al. 2014; T. J. Rodigas et al. 2015; D. Mawet et al. 2017). Due in part to high sky backgrounds, and the complications associated with observing through the Earth's atmosphere, these observations have only been possible so far for edge-on disks, or the brightest inner regions of non-edge-on disks.

With the launch of the James Webb Space Telescope (JWST; J. P. Gardner et al. 2023), it is now possible to observe at these longer wavelengths outside the Earth's atmosphere and with a much lower telescope thermal background, providing a major opportunity to expand the scope of disk scattered light and thermal imaging. Long-wave (≥2 µm) near-infrared (NIR) and mid-infrared (MIR) observations open up many exciting new prospects: e.g., (a) constraining dust properties through color measurements over a wide wavelength baseline from the MIR to the optical (when combining data from the Hubble Space Telescope, HST, and ground-based observatories; I. Rebolledo et al. 2024), (b) probing for the signature of water ice on dust grain surfaces through the

strong 3 µm absorption feature (e.g., M. Kim et al. 2019, 2024); and (c) conducting more sensitive searches for self-luminous planets, which may be sculpting the disk substructure (e.g., A. L. Carter et al. 2021).

The JWST Early Release Science (ERS) Program for Direct Observations of Exoplanetary Systems (Program ID 1386; S. Hinkley et al. 2022) has the primary goal of demonstrating these new scientific opportunities enabled by JWST's observing capabilities relevant to the exoplanet and circumstellar disk direct imaging community. Previous results from this program include: (a) the observation and characterization of planet HIP 65426b (A. L. Carter et al. 2023), the first extrasolar planet to be imaged with JWST; (b) the spectroscopic characterization of planetary-mass object VHS 1256b using the NIRSpec IFU and Mid-Infrared Instrument (MIRI) MRS modes, resulting in the detection of the first silicate cloud features in a planetary-mass companion, and delivering the most detailed spectrum yet for an object in this mass range (B. E. Miles et al. 2023); and (c) aperture masking interferometry of the star HIP 65426 (S. Sallum et al. 2024; S. Ray et al. 2025). These results have already led to follow-up analyses. S. Petrus et al. (2024) followed-up the VHS 1256b data from B. E. Miles et al. (2023) in an exploration of how well various self-consistent atmospheric models fit the extremely high-quality data from JWST, and N. Whiteford et al. (submitted) have carried out a detailed retrieval analysis from the same data set. Finally, early lessons learned and recommendations developed from the ERS program have been published in S. Hinkley et al. (2023).

In order to demonstrate JWST's ability to coronographically image faint extended emission, the ERS 1386 program includes both NIRCam and MIRI coronagraphic imaging of the HD 141569A disk. This paper focuses on the NIRCam observations while the MIRI observations will be reported in E. Choquet et al. (2025, in preparation).

### 1.1. HD 141569A's Disk

The disk around HD 141569A (B9.5V, 111.6 ± 0.4 pc; C. Jaschek & M. Jaschek 1992; Gaia Collaboration 2020) is commonly quoted to have an age of 5 ± 3 Myr

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
(A. J. Weinberger et al. 2000). However, we note that the age determined in A. J. Weinberger et al. (2000) used the Hipparcos distance of ~100 pc (M. A. C. Perryman et al. 1997) rather than the distance determined by Gaia (Gaia Collaboration et al. 2016), reported here. A farther distance would lower the absolute magnitude used for age estimation and push the system to a slightly younger age. The disk has an infrared (2–200 µm) excess luminosity, $L_{IR}/L_* = 8.7 \times 10^{-3}$ (N. Pascual et al. 2016), which lies in between the typical values associated with protoplanetary disks and debris disks. The excess level, in combination with the fact that significant gas has been detected in the system (B. Zuckerman et al. 1995; W. F. Thi et al. 2014; K. M. Flaherty et al. 2016; J. A. White et al. 2016; J. M. Miley et al. 2018; E. Di Folco et al. 2020) has resulted in it being classified as a “hybrid” disk (e.g., Á. Kóspál et al. 2013; J. Péricaud et al. 2017), in between the protoplanetary and debris phases.

The disk was first imaged in scattered light with HST by NICMOS (J. C. Augereau et al. 1999; A. J. Weinberger et al. 1999), and subsequently with the Space Telescope Imaging Spectrograph (STIS; D. Mouillet et al. 2001) and the Advanced Camera for Surveys (M. Clampin et al. 2003). These observations revealed two ring-like components at ~200 and ~325 au radius, and significant nonaxisymmetric structure possibly caused by dynamical interactions of an unseen planetary perturber (M. C. Wyatt 2005), or the pair of nearby M-stars (D. R. Ardila et al. 2005; A. C. Quillen et al. 2005). A third disk component, inside of ~100 au, was first imaged in thermal emission at 11 µm by M. M. Moerchen et al. (2010), and subsequently detected from ~46–116 au in scattered light using STIS by M. Konishi et al. (2016).

Ground-based follow-up has traced the scattered light emission from the H band to the L band down to angular separations of ~0.25, revealing the innermost component to be a highly asymmetric ring with surface brightness peaking at ~44 au and with tentative hints of faint spiral structures (T. Currie et al. 2016; C. Perrot et al. 2016; J. S. Bruzzone et al. 2020; G. Singh et al. 2021; C. Xie et al. 2022). A possible interpretation of the detected asymmetry is a recent collision between large planetesimals that causes the dust density to vary azimuthally (G. Singh et al. 2021). Alternatively, the asymmetric disks could be the result of dust–gas interaction (W. Lyra & M. Kuchner 2013; A. J. W. Richert et al. 2018).

This work focuses on JWST/NIRCam observations of the system, at 3 and 3.6 µm. The two primary science goals of the JWST/NIRCam observations are to: (a) leverage NIRCam's access to the 3.0 µm water ice feature to search for evidence of icy disk grains and measure their spatial distribution, and (b) search for planets that may be driving the dynamical features seen across the disk. A programmatic goal of the ERS program as a whole is to offer recommendations to the community on the most effective observational practices, as well as to provide "science enabling products." In this respect, we present a repository where all of the data products and analysis presented herein can be found in github<fnref n="53" /> and in Zenodo at doi: 10.5281/zenodo.16734248.

## 1.2. Water Ice in Circumstellar Disks

Given the critical role of water in the development and evolution of life on Earth (e.g., F. Westall & A. Brack 2018),

the formation and long-term stability of water ice within exoplanetary systems is of high interest. The formation of water ice onto dust grains in protoplanetary disks likely plays a critical role in the formation of planetary cores by enhancing the abundance of solid material outside the ice line and increasing the stickiness of grains in collisions (e.g., M. Min et al. 2011). Water ice is ubiquitous in the rings and moons of the outer solar system (e.g., G. Filacchione et al. 2012). The delivery of water to forming or recently formed planets by material outside the snow line is thought to be a primary mechanism of accreting water onto terrestrial planets, including the Earth (e.g., D. P. O'Brien et al. 2018). In this respect, understanding the prevalence of water ice in circumstellar disks across their lifetimes is critical to understand the water life-cycle in planetary systems.

In scattered light, there are ice absorption features centered at roughly 2, 3, and 4.5 µm, where the 3 µm feature is strongest by more than an order of magnitude (A. K. Inoue et al. 2008). The 2.9–3.2 µm feature is most significant for smaller grain sizes (<1 µm), but it is also present for larger grain sizes. However, photosputtering by UV photons may reduce the prevalence of ice (A. Grigorieva et al. 2007; A. Oka et al. 2012). In the optically thick case, if the scale height of the disk increases as a function of radius, then this effect may be more prominent at larger radii.

Ground-based scattered-light water ice absorption measurements have been carried out on five protoplanetary disks (HK Tau B, HV Tau C, HD 142527, HD 100546, and AB Aur), the debris disk HR 4796A, and HD 141569A. The two edge-on disks HK Tau B and HV Tau C were first shown by H. Terada et al. (2007) to exhibit water ice absorption using Subaru IRCS spectroscopy. Interestingly HV Tau C showed significant time-variability, suggesting a localized dust cloud moving in and out of the line of sight. Using multiband imaging at K and L with Subaru/CIAO, in combination with a medium-resolution filter centered at 3.08 µm, M. Honda et al. (2009) reported the detection of water ice in one spiral arm of HD 142527. Similar observations were carried out by M. Honda et al. (2016), who also detected statistically significant absorption in the HD 100546 disk in a 3 µm filter using Gemini/NICI. Interestingly, they compared their measurements against the UV H₂O photodesorption model by A. Oka et al. (2012) and found that their radial profile slightly preferred the model with photo-desorption as compared to their model without it.

S. K. Betti et al. (2022) observed AB Aur with similar broadband filters and also found statistically significant absorption from H₂O ice. Their analysis revealed that grain mass fraction of water ice inferred upon comparison with MCFOST Mie-scattering models depends on whether one compared the ice band to the $K_s$ observations (0%–5% ice by mass), or the $L'$ data (~68%). At the distance where their measurements were made, they argue that the mass fraction is unlikely to be as high as 68%, and that the $L'$ scattering surface may be at a different height in the disk relative to the dominant water absorption.

The results are varied for the older disks surrounding HR 4796A and HD 141569A. T. J. Rodigas et al. (2015) presented a detailed dust grain compositional analysis of HR 4796A using 11 photometric points and concluded that water ice was either not present or present only at low abundance in the system. Recently, J. K. Kueny et al. (2024) published ground-based observations of HD 141569A with the Clio-2 camera

<footnote n="53">⁵³ https://github.com/maxwellmb/jwst_ers_NIRCAM_HD141569A</footnote>
<table><caption>**Table 1**
Instrument Configuration</caption><thead><tr><th>Target</th><th>Filter</th><th>Readout Mode</th><th>N<sub>ramps</sub></th><th>N<sub>ints</sub></th><th>Exp. Time (s) per Dither</th><th>Roll Angle</th><th>N<sub>Dthers</sub></th></tr></thead><tbody><tr><td>HD 141569A</td><td>F300M</td><td>DEEP8</td><td>10</td><td>20</td><td>4019.59</td><td>279.76</td><td>1</td></tr><tr><td>HD 141569A</td><td>F300M</td><td>DEEP8</td><td>10</td><td>20</td><td>4019.59</td><td>284.76</td><td>1</td></tr><tr><td>HD 141569A</td><td>F360M</td><td>DEEP8</td><td>12</td><td>17</td><td>4207.74</td><td>279.76</td><td>1</td></tr><tr><td>HD 141569A</td><td>F360M</td><td>DEEP8</td><td>12</td><td>17</td><td>4207.74</td><td>284.76</td><td>1</td></tr><tr><td>HD 140986</td><td>F300M</td><td>RAPID</td><td>40</td><td>9</td><td>384.85</td><td>282.46</td><td>9</td></tr><tr><td>HD 140986</td><td>F360M</td><td>BRIGHT1</td><td>24</td><td>9</td><td>436.17</td><td>282.46</td><td>9</td></tr></tbody></table>

behind the MagAO-X adaptive optics system. They detected the innermost ring at ~40 au (0.36) radius in three filters centered at 3.1 µm (narrowband), 3.3 µm (broadband), and 3.8 µm (broadband). Their measurements recover a significant deficit of emission in the 3.1 µm filter, consistent with substantial water ice absorption.

In this paper, we present new water ice measurements of the HD 141569A disk, leveraging JWST's stability and sensitivity to complement these previous observations by detecting all three previously imaged ring components of the system in and out of the ice absorption band. We present the observations in Section 2, followed by a description of our procedure to subtract a model of the nearby M-stars and the results of point-spread function (PSF) subtraction in Section 3.2. In Section 4 we present a deconvolution method to recover the uncon-volved scene in both filters. Our forward modeling method is presented in Section 5, where we derive a measure of the absorption in the F300M as compared to the F360M filter. We then present a toy model that can explain these measurements (consistent between the forward modeling and deconvolution processes) in Section 6.1. Finally, we present our sensitivity to perturbing planets in Section 7 before wrapping up our final conclusions in Section 8.

## 2. Observations

Observations of HD 141569A were obtained on 2023 February 17 using the JWST/NIRCam coronagraphic mode in the F300M (Pivot $\lambda = 3.00$ µm) and F360M (Pivot $\lambda = 3.62$ µm) filters in the long wavelength channel of NIRCam (J. E. Krist et al. 2007; J. H. Girard et al. 2022; M. J. Rieke et al. 2023). Exposures were obtained in each filter at two telescope roll angles (~280° and ~285°) to allow for angular differential imaging (ADI; Table 1). Although the observatory can roll up to 14° (at large solar elongations), the maximum available roll at any given time is limited due to observatory pointing constraints. The star HD 140986 ($V = 6.25$ mag; SpT = K0III; N. Houk & C. Swift 1999) was observed as a PSF reference star using a 9-grid dither pattern to sample possible pointing acquisition error offsets in placing HD 141569A behind the coronagraph mask. The exposure times of the reference star were chosen such that the signal-to-noise ratio (SNR) of the stellar speckles in each dither was comparable to a single science exposure. All observations were obtained using the MASK335R coronagraphic focal plane mask and the MASKRND pupil plane mask (J. E. Krist et al. 2010). We used the 320 × 320 pixel subarray, which provides a field of view $20''/16 \times 20''/16$ in the long wavelength channel given its pixel scale of 0.063.

## 3. Data Reduction

The uncalibrated JWST/NIRCam data of HD 141569 were downloaded from the Mikulski Archive for Space Telescopes<fnref n="54" /> and then reduced to flux-calibrated images using the spaceKLIP community pipeline (J. Kammerer et al. 2022; A. L. Carter et al. 2023). The specific observations analyzed can be accessed via doi 10.17909/ty1h-9x40. SpaceKLIP employs a customized version of the official *jwst* data reduction pipeline (H. Bushouse et al. 2023), here used in its 1.12.1 version together with calibration file *jwst_1146.pmap*. The *pmap* files define the mapping between a given science file and its reference files, such as darks, bad pixel maps, and PSF masks. By default, the pipeline always uses the most recent available *pmap* file. However, if the results from this paper shall be reproduced, past *pmap* files can be downloaded at https://jwst-crds.stsci.edu/ and manually applied for the data reduction. With *jwst_1146.pmap*, new photometric flux calibration files for NIRCam coronagraphy derived based on Cycle 1 commissioning data have been used. Briefly, the customization of the *jwst* pipeline consists of not flagging pixels located diagonally from saturated pixels, adding a 4 pixel wide border of “pseudo-reference pixels” at the edges of the coronagraphic subarrays, skipping the dark subtraction step, and applying a jump rejection threshold of 4. The motivation for these changes with respect to the default *jwst* pipeline can be found in A. L. Carter et al. (2023). Recently, an additional step to correct for 1/f noise has been added to spaceKLIP, which was applied to the HD 141569 data as well. A more detailed description of this step can be found in I. Rebollido et al. (2024). Here, we used a Savitzky-Golay filter to model and subtract the 1/f noise, including data points up to half of the saturation fraction in the fit.

Next, the flux-calibrated images were cleaned of bad pixels, also using spaceKLIP routines. First, a bad pixel map is constructed from the pixels flagged by the *jwst* pipeline and additional bad pixels identified by a sigma clipping routine. On average, ~0.2% additional bad pixels were identified. Then, the flagged pixels are cleaned using three different algorithms applied subsequently. The first one replaces bad pixels, which are only bad in some (but not all) of the individual integrations with the median value of the pixel in the good integrations. The second one replaces remaining bad pixels with the median of their left/right/bottom/top good neighbors, and the third one uses an image plane median filter of size 4 pixels to replace remaining bad pixels that were not corrected by the previous two algorithms. We note that this procedure is exactly the same as in A. L. Carter et al. (2023). Finally, all

<footnote n="54">⁵⁴ https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Left panel: a view of the entire coronagraphic F300M field of the HD 141569 system. The two unocculted M-star companions are located to the upper left. The color map has been chosen to highlight the coronagraphic PSF and as a result the M-star PSF cores visually appear saturated. In reality, the detector has not been saturated by the M-stars. Center panel: a zoom in on the central regions of HD 141569A, shown with a log image normalization. Right panel: a zoom in on the central regions of the PSF reference star.</figcaption></figure>

integrations are median-combined into a single image for each roll position (or dither position for the PSF reference images).

An example of one of the two rolls of the F300M calibrated data can be seen in Figure 1. Within the image are the two M-star companions HD 141569 B and C in the upper left of the image, and the coronagraphic PSF of the primary star centered in the frame. Comparing the center and right panels, it is apparent that extended disk emission can be seen in the raw data. The center panel also highlights the extended PSF wings of the two M-star companions, which, without subtraction, biases the disk emission. The second F300M roll image and the two roll images obtained at F360M all share similar qualitative features.

The next step in our analysis was to model the M-stars with synthetic PSFs, generated using *WebbPSF* (M. D. Perrin et al. 2014; recently renamed to *STPSF*). In order to optimize these models, we first carried out a multistep process to refine our instrument model, described below in Section 3.1. The resulting instrument model is not only used to generate the M-star models, but also a suite of off-axis PSFs that were used for forward modeling (Section 5) and deconvolution (Section 4).

### 3.1. Instrument Modeling

To improve the nominal synthetic PSF models used throughout this work (e.g., disk forward modeling in Section 5), we optimized the following tunable instrument parameters in *WebbPSF*: defocus, pupil shear, and pupil rotation. For this purpose, we used each filter's set of nine dithered reference star exposures. Here, we assumed that the position of the coronagraph mask in the TA filter (F335M for the long wavelength channel) is precisely known and that the filter-dependent displacements of that position for F300M and F360M are well characterized by prior measurements. Likewise, based on prior analysis (J. H. Girard et al. 2022), we assumed that the inaccuracy of the dither maneuver and any drifting of the telescope are both negligible—such that the position of the star throughout the dithered sequence is well described by a single position plus each exposure's prescribed dither. The FGS guiding data acquired during the observations confirms this assumption, indicating a guiding jitter stability of <1 mas over the course of all exposure, and dither offsets of <1 mas from nominal for each position. As the morphology of the instrumental response for NIRCam coronagraphy is significantly affected by the position of the source relative to

the coronagraph, the offset of the reference star from the coronagraph center must also be considered when evaluating the aforementioned instrument parameters. All together, this left six parameters to be optimized: defocus, x-axis pupil shear, y-axis pupil shear, pupil rotation, x-axis source offset, and y-axis source offset.

The parameters are fit for all nine reference dithers simultaneously (i.e., we seek the set of WebbPSF parameters that best explains all nine dithers with consideration for the dither offset of each image). We began by fitting the F300M data. To generate each *WebbPSF* model of HD 140986, we adopted an approximate AMES-Cond spectrum for the reference star HD 140986 ($T_{\text{eff}} = 4800$ K, $\log(g) = 4.0$; F. Allard et al. 2001), and used the nearest measured optical path difference (OPD) map to the time of the observations. As each trial solution for the considered parameters required generating a separate *WebbPSF* model for each of the nine dithers, this process was potentially very computationally intensive. However, this burden was eased significantly by the fact that many of these parameters are not significantly correlated with one another. This allowed the procedure to be separated into multiple smaller problems that each considered only a subset of the six parameters. First, we determined the source offset. To do this, we generated a PSF model cube having null values for all six parameters. Using cross correlation, we identified the shift necessary to align the perfectly centered model PSF with the reference star data within a circular region extending to $40\lambda/D$ from the assumed coronagraph center. The x- and y-axis source offsets were fixed to the determined offset for the remainder of the procedure.

Next, we optimized pupil rotation alone by considering the annular region spanning $20\lambda/D < r < 40\lambda/D$, where $r$ is the distance to the coronagraph center. Within this region, we sought to minimize the uncertainty-weighted squared residuals between the data and model cubes (where the uncertainty for each exposure is propagated from the integration-wise pixel uncertainties stored in the FITS "ERR" extension). As the pupil shear and defocus have almost no effect at these separations, their values could be safely fixed to zero while fitting pupil rotation. The pupil rotation was then fixed to the best-fitting value hereafter.

Finally, we fit the defocus, x-axis pupil shear, and y-axis pupil shear simultaneously, considering the region within $r < 40\lambda/D$ but otherwise adopting the same objective as for the pupil rotation fitting. Of the six parameters, only defocus is
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The median-collapsed calibrated data (left panel), the M-star model (center panel), and the model-subtracted data (right panel) for the first roll of the F300M data. All three panels are shown with the same logarithmic color stretch. The color stretch highlights the extent of the M-star PSF wings, which can be seen to extend to a ~10'' radius.</figcaption></figure>

expected to have wavelength dependence. As such, the procedure for the F360M data adopts the F300M parameter values with the exception of defocus, which was fit following the same procedure. The best-fit values for x-axis pupil shear, y-axis pupil shear, and pupil rotation were -1.34e-2, -0.86e-2, and -0.59, respectively (with pupil shear in units of fraction of primary). The best-fit defocus values were -0.062 and -0.069 for F300M and F360M, respectively (in units of waves of defocus).

## 3.2. M-star Fitting and Subtraction

Using the best-fit *WebbPSF* parameters, we proceeded with generating models of the M dwarf companions as follows. Each M dwarf companion PSF model was generated assuming an approximate AMES-Cond synthetic source spectrum ($T_{\text{eff}} = 3400$ K, $\log(g) = 4.5$; F. Allard et al. 2001), and using the nearest measured OPD map to the time of the observations. Throughout this analysis, we excluded any pixels that were replaced during the bad pixel cleaning procedure. For each source, we first determined an approximate position through visual inspection of the derotated images and then generate a PSF model cube corresponding to that position. The brightness of each initial PSF model was set to the brightness that minimizes the uncertainty-weighted squared residuals between the data and model cube when considering only pixels far from HD 141569A (at least 80 pixels) and where the PSF model has nonnegligible values (values above $10^{-4}$ times the maximum model value). Since PSF morphology is not expected to change meaningfully in the vicinity of the sources, we determine the final position of each source by cross-correlating the initial model with the data. For this final position, a new PSF model is generated and scaled in the same manner. No effort is made to mask the other source during this process; while the wings of the sources' PSFs overlap, the PSF cores are sufficiently separated that additional masking was found to have no effect on the result of cross correlation (with the sources being separated by approximately $11\lambda/D$).

While the determined positions were not affected by the proximity of the companions, the brightnesses were. To remedy this, we proceeded with a final optimization step to hone the brightness of both models simultaneously, before adding them to form a composite model cube. Here, we also corrected for

two additional effects. The first is an overall erroneous brightness offset for the data induced by the previous use of pseudo-reference pixels during ramp fitting and the subsequent background subtraction step performed by *spaceKLIP*. As the background offset differs for each roll, the brightness offset for each roll must be tuned separately, and is introduced by simply adding a background value to each roll's composite model image.<fnref n="55" /> Finally, we introduced a Gaussian blurring to the composite model to approximate detector effects that are not included in the *WebbPSF* models, in particular the effect of charge migration (sometimes also referred to as the brighter-fatter effect). This effect affects the HAWAII-2RG detectors used in all NIR JWST instruments, and we refer the reader to P. Goudfrooij et al. (2024) for a more detailed description. This resulted in an optimization procedure varying five parameters: a brightness scaling factor for each model, a brightness offset for each roll, and a Gaussian blurring kernel width. During the optimization, we again excluded the region within 80 pixels of HD 141569A, but now included any pixels where either model cube is above 1–6 times that model cube's maximum value. As before, the optimization procedure sought the values of these five parameters such that the error-weighted squared residuals between the resulting composite model and the data are minimized. After applying this technique for the F300M filter, the same source positions were used to generate the PSFs for the F360M model. Otherwise, following the same procedure as for the F300M models, the initial brightnesses of the models were then scaled, and the five-parameter model honing procedure was applied.

Figure 2 displays an example F300M M-star model and the residuals after subtraction. The M-star model represents the data well and the subtraction significantly reduces the contribution of the M-stars near the disk's extended emission. Systematic residuals remain in the image, but are typically a small fraction of the initial signal, in particular near the extended disk emission. Residuals are likely due to the incompleteness of the internal NIRCam wavefront error as modeled by *WebbPSF/STPSF*. Because the coronagraphic fields of view lie outside the traditional imaging field and contain additional optical

<footnote n="55">⁵⁵ While K. Lawson et al. (2024) considers the effects of the coronagraph and neutral density squares for background correction, these effects are negligible for the source brightnesses and background levels of the data considered here.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Roll-combined images after PSF subtraction. For each filter, the PSF subtraction was performed using MCRDI. The red dashed line delineates the radius of 50% transmission inside of the coronagraph.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison of disk structures between archival HST/STIS optical coronagraphic data (left panel) and the JWST/NIRCAM F300M observations (right panel). A key difference in the morphologies is that the JWST data do not show a well-defined deep gap between the A and B ring components. North is up, and east is left.</figcaption></figure>

components (field and pupil elements), these could not be measured during ground testing and rely on optical models rather than the as-built measurements.

### 3.3. Disk Structure in PSF-subtracted Images

After the M-star model was subtracted, we attempted to carry out reference star differential imaging (RDI) PSF subtraction with *spaceKLIP* using Karhunen–Loève image projection (KLIP via pyklip; R. Soummer et al. 2012; J. J. Wang et al. 2015) and nonnegative matrix factorization (B. Ren et al. 2018). However, a bright central disk component introduced bias into the PSF-subtraction procedure that resulted in significant oversubtraction of the PSFs. As a result, we carried out model-constrained RDI (MCRDI; K. Lawson et al. 2022) in order to simultaneously model the disk and optimize the PSF subtraction.

Figure 3 presents the PSF-subtracted image in both filters after the MCRDI process (described in detail in Section 5 below). Figure 4 shows archival HST/STIS coronagraphic observations (GO–13786) that we processed following the methods described in M. Konishi et al. (2016). These

are diffraction-limited optical images ($\lambda$ = 585.2 nm, $\Delta\lambda$ = 441.0 nm) with angular resolution 0.′10. The complex disk structure has three main components (A, B, and C), and the intermediate belt B is split along the eastern edge, and the outer belt A has tenuous spiral arms to the north and south.

The JWST images show disk structure from inside the coronagraph inner working angle (i.e., the radius with 50% throughput, 0.′64), out to a radius of ~4′′ (a similar extent as that seen with STIS). The NIRCam data also recover the split in the B ring seen by STIS. The emission inside the coronagraph inner working angle appears as a two-lobed structure, surrounded by a halo of nodes extending to a radius of 1.′5 along the disk’s semimajor axis. This clumpy structure corresponds to the bright inner disk C at ~44 au interacting with the radially changing coronagraphic throughput and the convolution with the complicated off-axis coronagraphic PSF. These structures extend to larger radii in the F300M than the F300M, as would be expected for a diffraction-related effect. A north–south asymmetry is marginally visible in the inner-most regions that may correspond to the asymmetries seen by various ground-based observations. However, the limited inner
working angle of the coronagraph prevents a direct comparison
with that feature.

The complex structures of the B ring, including the eastern
split feature, are clearly resolved in the JWST images. The A
ring does not have an inner edge that is as sharply defined in
the JWST data as in the HST data. Finally, the northern spiral
arm is marginally detected in the JWST data.

## 4. Deconvolution

We perform deconvolution on PSF-subtracted images using
a version of the iterative Richardson-Lucy (R-L;
W. H. Richardson 1972; L. B. Lucy 1974) algorithm that
has been modified to accommodate both the spatially varying
morphology of the NIRCam PSF and transmission of the
coronagraph. The key modification of this algorithm is
summarized hereafter, but the technique will be more formally
reported in an upcoming work (K. Lawson et al. in
preparation).

Let *I* be the observed 2D image array to be deconvolved and
let *B* be an operator that applies instrumental blurring (and any
throughput losses) to a 2D image, such that if *O* is the ground-
truth unblurred image, then *I* = *B*(*O*). Then, let the *i*-th
iteration of the deconvolution procedure be denoted *R*<sub>*i*</sub>, with
the subsequent iteration, *R*<sub>*i*+1</sub>, being computed by multiplying
the previous iteration by a correction term. Finally, let all
division and multiplication refer to element-wise (Hadamard)
operations.

Consider first the scenario of a noncoronagraphic, nonvary-
ing PSF. For conventional R-L deconvolution, the correction
term is computed as $\mathcal{B}[I/\mathcal{B}(R_i)]$—we blur the $i$th iteration,
divide the observed image by the result, and then blur that
quotient. In the case of a nonvarying PSF (and neglecting
noise), if the $i$th iteration is exactly the true unblurred image,
then the quotient is unity and blurring the quotient likewise
produces unity, such that the resulting correction term is unity.
In this case, the procedure has converged to the correct
solution, and no further corrections are made.

Now, consider a case comparable to NIRCam coronagra-
phy: where both the PSF morphology and throughput vary
with detector position. Here, when computing the R-L
correction term as above in the scenario where the *i*th iteration
is exactly the true unblurred image, we do not recover a
correction term of unity; while the quotient in the correction
term still evaluates to unity, the blurring of that quotient no
longer does.<fnref n="56" /> Denoting a 2D array of unity as $\mathcal{J}$, we account
for this simply by dividing the correction term at each iteration
by the result of blurring $\mathcal{J}$, such that the modified correction
term is $\mathcal{B}[I/\mathcal{B}(R_i)]/\mathcal{B}(\mathcal{J})$. With this modification, the conv-
ergence property of R-L deconvolution is recovered. This
procedure is run for 250 iterations on each roll, with the two
rolls then being averaged to produce the final deconvolved
image for each filter (Figure 5).

## 5. Disk Modeling

The purpose of our disk modeling was twofold: (i) to
provide a model to aid in PSF subtraction via MCRDI, and (ii)
to establish the ratio of scattered light between the F300M and
F360M filters as a function of stellocentric distance before
being convolved with the instrumental PSF. We assume that

the disk emission seen in our images is entirely due to
scattered light, rather than thermal emission from the dust. We
confirm that this assumption is reasonable at the end of
Section 5.1. Our modeling process involves two main steps.
First, a baseline disk model was established to carry out the
MCRDI PSF subtraction and to determine basic morphological
parameters (Section 5.1). In MCRDI, a disk model and a PSF
model are simultaneously fit to the data. Indeed, K. Lawson
et al. (2023) demonstrated that this technique provided
improved PSF subtraction in the AU Mic disk with NIRCam
observations as well. Here, the PSF model is composed of a
weighted average of the reference star images, where the
weights are free parameters in the fit. The disk model is
described below. Second, we freeze many of the disk
parameters, and carry out a Markov Chain Monte Carlo
(MCMC) analysis of the flux-related parameters in order to
evaluate the radial flux distributions in each filter
(Section 5.2).

### 5.1. Baseline Disk Model

We adopted the GRaTeR analytical ring prescription
(J. C. Augereau et al. 1999) as our disk model, and modified
it to adopt using the JAX framework (described below). The
GRaTeR model provides an analytical description of the radial
and vertical dust profile of a disk and can calculate scattered
light images based on a provided scattering phase function. We
refer the reader to J. C. Augereau et al. (1999) for a detailed
description of the analytical model.

To model our data, we generated images consisting of four
separate rings. Although this system is known to have
azimuthal asymmetries, here we focus on the radial dust
distribution and leave a study of azimuthal dust properties for
future work. We experimented with both three-, four-, and
five-ring models, but found that three rings provided
unsatisfactory images, leaving significant residuals, and that
using five rings provided nearly identical results to the four-
ring model.

As our primary goal is to estimate the pre-convolved surface
brightness, and not carry out detailed disk modeling, we opted
to fix several of the disk morphological parameters (using
symbols as defined in J. C. Augereau et al. 1999): the vertical
scale exponent (set to be a Gaussian vertical profile; $\gamma = 2$), the
reference scale height (set to be 3 au at the disk central radius;
$\zeta_{r0} = 3$), the flaring index (set to have linear radial flaring;
$\beta = 1$). As free parameters, each of the four rings had a unique
flux scaling parameter (related to $n_0$, but not strictly a dust
density parameter), central radius ($r_c$), and inner and outer
radial power laws ($\alpha_{in}$ and $\alpha_{out}$). For each ring, we adopted a
scattering phase function consisting of the weighted sum of
three Henyey-Greenstein (H-G) functions (L. G. Henyey &
J. L. Greenstein 1941) with three weights and three g
parameters. We attempted the fit with only two H-G functions,
but found the fits to be significantly worse. In this system, with
significant azimuthal structure in the disk, the scattering phase
functions are likely modeling both the true phase function and
also azimuthal changes in the dust grain density. Therefore,
each ring had 10 free parameters. We applied the same
inclination and on-sky position angle of the (projected)
semimajor axis of the disk to all four rings. The rings were
generated at the same resolution as the images, i.e., with pixel
sizes of 0.′′063.

<footnote n="56">⁵⁶ This is true even for noncoronagraphic data if the PSF morphology is shift variant.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Deconvolved images of the PSF-subtracted disk in the F300M (left panel) and F360M (center panel) filters, as compared to the pre-convolved model from Section 5.</figcaption></figure>

Once a given disk model was generated, we convolved it with a grid of spatially varying NIRCam off-axis PSFs that were generated with *WebbPSF* using the instrument parameters (defocus, x- and y-axis pupil shear and pupil rotation) that were fit for in Section 3.2. The off-axis PSF images were generated at six radial distances from the coronagraph center (plus one PSF at the coronagraph center), out to a maximum distance of 4", spaced logarithmically, such that the inner regions where the PSF is changing more rapidly had denser sampling than the outer regions. At each radial location, PSFs were generated at four azimuthal angles. To perform the convolution, each pixel in a generated disk image was individually multiplied by the off-axis PSF image generated at the nearest location, and the results were summed together. This process replicated that carried out by K. Lawson et al. (2023) with NIRCam images of AU Mic. For each of the two roll images, we included as free parameters the centering of the four-ring scene behind the coronagraph (four additional parameters; 2× [x, y] centers).

Finally, for each roll, we include as free parameters nine weight parameters that are used to combine the PSF reference star images (from nine dithers) as a weighted sum to create a model host star PSF that was then subtracted from the data. Each roll angle had its own weight parameters, adding a further 18 free parameters. Preliminary testing of this process indicated that re-centering of the PSF reference images was not necessary to subtract the PSF to negligible levels below the disk emission.

In total, our model included 64 parameters: 4 × 10 ring parameters, two common ring parameters (inclination and position angle), 2 × 2 star centers (one for each roll), and 18 for the PSF weights. In order to optimize over this large parameter space, we adapted the GRaTeR code included as part of the VIP software suite (C. A. Gomez Gonzalez et al. 2017; V. Christiaens et al. 2023), and the code used in K. Lawson et al. (2023) to be compatible with the JAX python library for high-performance numerical computing.<fnref n="57" /> Here, we leverage its GPU-compatibility, including its just-in-time compiler to quickly calculate a disk image and its ability to efficiently calculate the gradient of a function via auto-differentiation. We defined our objective function as the squared sum of the residuals of both roll images after subtracting the forward-modeled disk and the weighted PSF

images, divided by the errors output by the *spaceKLIP* pipeline. We minimized the objective function using the *jaxopt*<fnref n="58" /> (M. Blondel et al. 2022) wrapper for the SciPy (P. Virtanen et al. 2020) *minimization* function, using the 'L-BFGS-B' algorithm (D. C. Liu & J. Nocedal 1989; R. H. Byrd et al. 1995; C. Zhu et al. 1997), which takes advantage of the gradient enabled by JAX.

In the case of the F300M images, we included all 64 parameters as free parameters in the minimization. Under the assumption that the morphology and scattering parameters do not change significantly between our two filters,<fnref n="59" /> for the F360M images, we froze the ring morphological and scattering parameters, leaving the following parameters free (a total of 26 parameters): the flux scaling of each ring (four parameters), the centering behind the coronagraph (four parameters), and the 18 PSF weights. The F300M filter was chosen because the morphological features appear sharper in this image as compared to F360M. The minimization was run until convergence, and best-fit parameters are given in Table 2. This optimization process does not return error estimates on the best-fit parameters, but at this stage, we are only interested in a representative model, not necessarily a quantitative model.

Figure 6 displays the roll-combined data, the model, and the residuals for both F300M and F360M. Before being combined, the data were PSF subtracted using the weighted-sum PSF model with the best-fit weights taken from the minimization process. The data and models were derotated using the telescope roll angle obtained from the `GS_V3_PA'` header keyword and then mean-combined. The final deconvolved images are compared against the forward-modeled results in and can be seen in Figure 5.

Figure 7 displays the four rings of our disk model, along with the images of the data with the model rings sequentially subtracted to highlight the different components of the data that the rings are modeling. The innermost ring (Ring 1 in Table 2) is responsible for the two-lobed structure as well as the knobby structure surrounding it. The second-most inner ring (Ring 2) has a similar radius as the first ring, but models a much more diffuse component in the data. The third ring (Ring 3) models the middle ring component, labeled the “B ring” in M. Konishi et al. (2016) and Figure 4. The sharpness of the

<footnote n="58">⁵⁸ https://jaxopt.github.io</footnote>

<footnote n="59">⁵⁹ In reality, the presence of significant water ice would indeed change the scattering phase function (M. Kim et al. 2024), but for simplicity we keep these parameters fixed.</footnote>

<footnote n="57">⁵⁷ https://jax.readthedocs.io/</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Left column: roll-combined images of the F300M and F360M data after PSF subtraction. Center column: roll-combined images of the models convolved with the NIRCam off-axis PSF. Right column: the data minus the model.</figcaption></figure>

<table><caption>**Table 2**
Baseline Disk Parameter Optimization Results</caption><thead><tr><th>Parameter Name</th><th>Symbol (s)</th><th>Ring 1<br>Common Parameters</th><th>Ring 2</th><th>Ring 3</th><th>Ring 4</th></tr></thead><tbody><tr><td>Inclination</td><td><i>i</i></td><td></td><td></td><td>58.37</td><td></td></tr><tr><td>Position angle</td><td>PA</td><td></td><td></td><td>-356.4</td><td></td></tr><tr><td>Roll 1 center F300M</td><td>[x<sub>r1</sub>, y<sub>r1</sub>]</td><td></td><td></td><td>[146.94, 172.92]</td><td></td></tr><tr><td>Roll 2 center F300M</td><td>[x<sub>r2</sub>, y<sub>r2</sub>]</td><td></td><td></td><td>[147.25, 172.71]</td><td></td></tr><tr><td>Roll 1 center F360M</td><td>[x<sub>r1</sub>, y<sub>r1</sub>]</td><td></td><td></td><td>[147.25, 173.01]</td><td></td></tr><tr><td>Roll 2 center F360M</td><td>[x<sub>r2</sub>, y<sub>r2</sub>]</td><td></td><td></td><td>[147.53, 172.82]</td><td></td></tr><tr><td colspan="6" style="text-align: center;">Individual Disk Parameters</td></tr><tr><td>Flux scaling</td><td>...</td><td>170.73</td><td>48003.33</td><td>448.32</td><td>193.82</td></tr><tr><td>Disk central radius</td><td>r<sub>c</sub></td><td>31.71</td><td>31.04</td><td>223.01</td><td>297.71</td></tr><tr><td>Inner power-law Index</td><td>α<sub>in</sub></td><td>6.50</td><td>2.57</td><td>6.68</td><td>4.41</td></tr><tr><td>Outer power-law Index</td><td>α<sub>out</sub></td><td>-0.41</td><td>-1.29</td><td>-3.41</td><td>-0.77</td></tr><tr><td>Scattering parameters</td><td>[g<sub>1</sub>, g<sub>2</sub>, g<sub>3</sub>]</td><td>[0.47, -0.70, 0.44]</td><td>[0.84, 0.64, -0.90]</td><td>[0.61, -0.63, 0.06]</td><td>[0.02, -0.19, 0.80]</td></tr><tr><td>Scattering parameter weights</td><td>[w<sub>1</sub>, w<sub>2</sub>, w<sub>3</sub>]</td><td>[0.47, 0.38, 0.14]</td><td>[0.00, 0.00, 1.00]</td><td>[0.60, 0.40, 0.00]</td><td>[0.47, 0.11, 0.42]</td></tr></tbody></table>

disk component and the convolution with the NIRCam off-axis coronagraphic PSF—known to have significant flux in its satellite lobes—results in noticeable features in the disk model above and below the main ring. The outermost ring (Ring 4) models the more extended diffuse component seen after subtracting the first three ring models. This final component also attempts to model some of the emission from the “split” feature seen attached to the B-ring, which our simple model is not able to otherwise reproduce.

The model provides a good representation of the data, capturing many of the predominant features seen in the PSF-subtracted disk images, including the complicated inner features caused by the coronagraph throughput and PSF wings, the sharp middle disk, and the diffuse emission seen at larger radii. In general, our residuals are ≲10% of the original data values. The predominant residuals can be mostly

attributed to the simplicity of our axisymmetric ring model in the presence of the known asymmetries at all radial scales in the disk, and in the inner regions may be in part due to a mismatch between our instrument model and JWST’s as-built coronagraph or the lack of oversampling in our model. Despite these discrepancies, overall our model represents the data well.

Our results have several differences relative to prior NIR ground-based work. B. A. Biller et al. (2015) detected and modeled the “inner ring” (our Ring 3 model) at 245 au and the “outer ring” (our Ring 4 model) at 406 au. Our model places these rings at 223 and 298 au possibly because the gap between the inner and outer ring of B. A. Biller et al. (2015) and other work is filled with 3 µm emission in the JWST data. J. Mazoyer et al. (2016) and C. Perrot et al. (2016) derived even smaller radii of 208 au and 197 au, respectively, for Ring B. One partial explanation is that oversubtraction of the PSF in
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Top row: the F300M data from one roll angle after sequentially subtracting the different ring models. The images are shown at four different display levels to show the disk components from the bright inner dust component to the fainter outer ring. Bottom row: the individual ring components of our convolved model. All images are shown at the orientation of the raw data.</figcaption></figure>

this prior work shifts the apparent edge of the ring closer to the star, whereas the JWST PSF subtraction relies on reference differential imaging. On the other hand, the diffuse disk emission (our Ring 4 model) may shift the apparent peak of the narrow belt (Ring 3) to a greater radius in our analysis. Exploring the exact reason for this discrepancy is beyond the scope of this work, which focuses mostly on modeling the radial surface brightness profile rather the detailed disk morphology.

T. Currie et al. (2016) detected the inner disk at L' and reported a radius of 39 au, whereas we derived 32 au. They reported a “hook-like” morphology at the ansae in both Keck and Very Large Telescope (VLT) data. D. Mawet et al. (2017) and J. K. Kueny et al. (2024) also imaged HD 141569A at ~3 µm from the ground and also reported a 39 au radius, but they characterized the southern hook as a spiral arm extending eastward. We do not find evidence for this arm in our data, though this region lies inside our inner working angle of 71 au, where the throughput is <50%.

In both filters and both rolls, the disk centers are found to be approximately 2 pixels off from the expected center of the coronagraph ([x, y] = 149.1, 173.4]). This is potentially due to both miscentering behind the coronagraph, as well as the fact that the innermost disk is known to be asymmetric (e.g., J. S. Bruzzone et al. 2020; G. Singh et al. 2021), yet we model it as a simple ring. Exactly how the asymmetric inner region interacts with the coronagraph at these smaller angular separations will depend on the exact morphology of the inner ring. Such modeling could be explored in future work.

In order to validate our assumptions that the majority of the light that we are seeing at these wavelengths has come from scattered light and not thermal emission, we constructed a toy model of the disk based on our best-fit results using MCFOST (C. Pinte et al. 2006, 2009). The model included the four rings, using the same radial dust profiles as in Table 2 and orientation (inclination and PA). The dust model consisted of compact grains, where the two inner rings were a mix of silicates, with 10% water by volume, and the outer two rings had 1% water

ice. Propagation was carried out using Mie theory. We set the masses of the components based on the three-ring model from W. F. Thi et al. (2014), where we split the mass of their innermost disk between ring 1 and ring 2. We used the masses for their outer two rings for rings 3 and 4. After carrying out the MCMC radiative transfer simulations, we split the contributions from scattered light and thermal emission and found that the thermal emission at 3 µm was 10 orders of magnitude fainter than the scattered light. Thus, we conclude that even if our model assumptions are slightly off, the emission seem by JWST in these filters is heavily dominated by scattered light.

## 5.2. Flux Ratio Modeling

At this stage in the modeling process, we were interested in quantitatively understanding the pre-convolved surface brightness profiles of the system in our two wavelengths. With that in mind, we froze the disk morphological parameters, the scattering phase functions, the PSF weights and the disk centering parameters, and left as free parameters only the flux scaling parameters of each of the four rings. We constructed a log-likelihood for the data given the model parameters as:

$$ \ln p(\text{data}|\Theta) = -0.5 \sum_{\text{rolls}} \sum_{n} \left[ \frac{(I_{\text{data}} - I_{\text{model}})^2}{s_n^2} + \ln(2\pi s_n^2) \right] \quad (1) $$

where $I_{\text{data}}$ and $I_{\text{model}}$ are the intensities of the data and model, respectively, $\Theta$ are the model parameters (the four flux scaling factors—one for each ring), and the sums are over the two roll images and the n central 140 × 140 (8" × 8") pixels in the image. The $s_n^2$ term is a variance term equal to the squared noise output by the JWST Stage 2 pipeline (as implemented by spaceKLIP), $\sigma^2$, added to an additional fractional model
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Left panel: the normalized surface brightness as a function of deprojected stellocentric distance for the data and the models. The vertical gray dashed lines show the locations of the A, B, and C rings identified by M. Konishi et al. (2016), and the vertical red dashed line shows the radial separation at which the coronagraph throughput reaches 50%. Right panel: the scattering efficiency ratio between the F300M and F360M images as a function of deprojected stellocentric distance. The vertical dashed lines represent the same values as in the left plot. The MCMC samples in both plots, which are fit to the noisy data, closely follow the data (orange) and demonstrate that the increasing trend seen in the data as a function of radial distance is detected with high significance.</figcaption></figure>

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

variance term, defined by $f$:

$$s_n^2 = \sigma^2 + f^2 I_{\text{model}}^2, \qquad (2)$$

where $f$ was included as a fifth free parameter. We include the additional fractional model variance to account for the fact that our simple model does not include the complicated azimuthal/ spiral structure seen in the real data, and we expect this additional error to scale with the disk flux.

The posterior probability distributions were sampled using the `emcee` MCMC package (D. Foreman-Mackey et al. 2013). Here we opt for using an MCMC sampler, rather than the simple optimization in Section 5.1 in order to recover parameter uncertainties and covariances. The sampler was set up with 32 walkers and run for an initial burn-in phase of 300 steps, after which the chains no longer moved significantly and visually appeared to have converged. The full run was carried out using 3000 steps, about 60 times the mean autocorrelation time. This procedure was carried out for both the F300M and F360M images. An examination of the joint posterior distributions reveals a correlation between the flux scaling of the two first rings (with similar semimajor axes), but all of the marginalized posterior distributions appear to be Gaussian in shape and well constrained (Figures 12 and 13 in Appendix A).

Figure 8 displays normalized surface brightness profiles for the data, the deconvolved data (see Section 4), the convolved model, and the unconvolved model. In each case, the profiles were built by measuring the median surface brightness in concentric elliptical annuli of width 2 pixels (in the direction of the semimajor axis) centered on the star’s location and divided by our estimate of the stellar flux in each filter. The convolved model radial profile is shown as a set of 100 randomly sampled models from the MCMC fitting described in this section. However, the scatter in the MCMC samples is so small that it appears as a single black line that matches the orange line of the data’s radial profile very well.

The right plot of Figure 8 displays the scattering efficiency ratio between the F300M and F360M images as a function of stellocentric distance for the raw data, the deconvolved data, the pre-convolved model, and the convolved model. These radial profiles were calculated by dividing the F300M images

by the F360M images, scaling by the ratio of stellar fluxes at F300M and F360M, and measuring the median value in concentric elliptical annuli with widths of 2 pixels. We estimate the stellar flux in the F300M and F360M filters using the `pysynphot` package (STScI Development Team 2013). We simulate the stellar spectrum using a Phoenix model at a temperature $T = 10,000$ K, $M/H = -0.5$, $\log(g) = 4.28$ (stellar parameters from B. Merín et al. 2004), which we scaled to the HD 141569A Two Micron All Sky Survey J-band magnitude $J = 6.872$ (R. M. Cutri et al. 2003), to avoid contamination by the warm disk infrared excess. We use the NIRCam spectral throughput from the `WebbPSF` (M. D. Perrin et al. 2014) package and obtain stellar fluxes of 712 mJy and 533 mJy in the F300M and F360M filters, respectively.

A 2D map of the scattering efficiency ratio is shown in Figure 9, which compares the ratio of the images from the two filters to a similar image generated from the models of each filter. Overall, the forward-modeled radial profile and 2D map matches the data well (where there is good disk signal). The match between the deconvolved data and the unconvolved model is not as exact, yet they qualitatively agree; the ratio between the two filters is ~0.4 inside of the B ring, and increases to at least 0.8 in the outer reaches of the system.

## 6. Grain Models

To interpret the surface brightness change as a function of radius described in Section 5.2, here we assess how a change in the minimum grain size or a change in the water ice fraction in grains affects the surface brightness ratio between 3 and 3.6 µm. We build a radiative transfer model using `DISC2RADMC`<fnref n="60" /> (S. Marino et al. 2022), a package that allows for the mixing of optical constants, the calculation of opacities using Mie theory,<fnref n="61" /> and the simulation of disk images using `RADMC3D`.<fnref n="62" />

Our disk model consists of a set of axisymmetric Gaussian rings that resemble the inferred dust distribution. Note that the

<footnote n="60">⁶⁰ https://github.com/SebaMarino/disc2radmc</footnote>

<footnote n="61">⁶¹ Models of nonspherical grains could lead to different results, and thus, future work could explore the effect of using nonspherical grain models.</footnote>

<footnote n="62">⁶² https://github.com/dullemond/radmc3d-2.0</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The ratio of the F300M and F360M data (left panel) and convolved models (right panel), after dividing each image by the integrated stellar flux across the bandpass. The three white ellipses correspond to the A, B, and C rings from M. Konishi et al. (2016).</figcaption></figure>

exact dust spatial distribution is not relevant to the grain properties as it does not affect the surface brightness ratio. We fix the size distribution to a power law with the commonly assumed exponent value of −3.5 for a collisional cascade (J. S. Dohnanyi 1969). The exponent affects how many small grains there are relative to the large grains, potentially changing the strength of the water ice feature. However, we choose to keep this fixed and instead vary the minimum grain size, as these can be degenerate. We also fix the maximum grain size to 1 mm. Its value has little effect as long as it is much larger than the wavelength and thus does not contribute significantly to the scattered light. We compute a grid of models varying the minimum grain size from 0.1–10 µm and the water ice fraction (by mass) from 0.1%–100%, in a 30 × 30 logarithmically spaced grid. Grains are assumed to be composed of amorphous water ice and silicates (B. T. Draine 2003; R. M. Mastrapa et al. 2009), and their optical constants are combined using the Bruggeman rule.<fnref n="63" /> We compute the dust opacities using Mie theory, assuming compact spherical grains. We produce images at several wavelengths to sample the width of each filter, and we then compute F300M and F360M model images by integrating over each filter transmission. These scattered light images use the Henyey–Greenstein approximation for anisotropic scattering. Finally, we extract surface brightness profiles normalized by the stellar flux, and we calculate the average surface brightness ratio between the two filters.

The model surface brightness ratio as a function of the minimum grain size and water ice fraction is presented in Figure 10 (top panel). We find that any ratio can be reproduced by a wide range of parameters. Namely, the ratio decreases with water ice fraction and increases with minimum grain size. This is because the water ice feature near 3µm becomes more prominent as the amount of water ice in grains increases and as an inverse function of grain size. We also find a valley or local minimum in the surface brightness ratio near a minimum grain

size of 2 µm. As a reference, the white and black contours represent a brightness ratio of 0.4 and 0.8 as found near 100 and 400 au, respectively, representative values of the minimum and maximum values. The dashed orange line represents the expected minimum grain size estimated as those grains with a β value of 0.5, β being the ratio between the radiation and gravitational force, taking into account the scattering opacity and the degree of forward scattering (J. A. Burns et al. 1979). In the absence of gas, we would expect the minimum grain size to be close to this orange line, as smaller grains would be efficiently removed by radiation pressure.

To understand how these results could vary if we considered porous grains, the bottom panel of Figure 10 presents the surface brightness ratio assuming a porosity of 50%. We find that the porosity tends to decrease the surface density ratio, shift the blow-out size toward larger sizes and the valley from 2 to ~7µm.

Overall, we find that we can explain the increase in the brightness ratio from 0.4 at 100 au to 0.8 at 400 au by either decreasing the amount of water ice by a factor of 3–10 (e.g., from 10% to 1% in mass) or by increasing the minimum grain size from 2–10 µm while maintaining the water ice fraction at 10%, assuming nonporous grains. As the porosity increases, it quickly becomes more difficult to reproduce these results by a strict change in minimum grain size.

At first, it seems unlikely that the amount of water ice could decrease with distance, as any desorption process should act the other way, increasing the amount of water ice with distance. However, if the relative amount of other volatiles increased with distance as they freeze out, that could lower the abundance of water ice. In a second scenario, the minimum grain size may increase with radius if the inner regions are richer in gas (as may be expected in this system; J. M. Miley et al. 2018; J. A. White & A. C. Boley 2018; E. Di Folco et al. 2020) and grains smaller than the blow-out size could be dragged sufficiently. However, our model in Figure 10 shows that a change in the minimum grain size from below the blow-out size to the blow-out size (orange line) alone cannot explain the change in the surface brightness ratio. Both the water ice

<footnote n="63">⁶³ Almost the same result was obtained mixing these using the Maxwell Garnett equation with water ice as an inclusion.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Expected surface brightness ratio between F300M and F360M as a function of the minimum grain size and the water ice fraction in grains (by mass). The top panel presents a case of compact grains (0% porosity), while the bottom panel presents a porous case (50% porosity). The white and black contours represent a brightness ratio of 0.4 and 0.8 as found at 100 and 400 au, respectively. The orange dashed line represents the expected minimum grain size as a function of the ice fraction.</figcaption></figure>

fraction and minimum grain size may change with radius. The change in the surface brightness could also be reproduced by a change in the porosity of grains from 50% porous to 0% with radius while maintaining the minimum grain size at 10 µm and water ice at 10%. This could be the case if the origin of dust in the inner and outer regions was different (e.g., collisionally generated in a collisional cascade as in debris disks versus protoplanetary disk dust).

This is overall a degenerate problem if we only consider these two filters, a result that is also found through modeling by M. Kim et al. (2024). Considering the surface brightness at a third wavelength or the spectral energy distribution (SED) could lift these degeneracies, but doing so is beyond the scope of this paper. We recommend that future work trying to constrain the presence of water ice using scattered light images considers at least three wavelengths in their analysis, with two wavelengths bracketing the water ice feature. At the time of our observations, simultaneous observing with NIRCam short-wave and long-wave channels was not yet enabled. However, now it is possible to observe in this mode, facilitating a three-wavelength comparison more easily. Furthermore, future work incorporating information from the absolute scattered light flux, phase function, polarization fraction, MIR surface brightness using MIRI and VLT/VISIR (E. Choquet et al.

2025, in preparation), and the SED could lift these degeneracies between the composition and size distribution of grains.

We also consider how the presence of other contaminants would affect the brightness ratio in these two filters. Using the Mie theory calculations as our baseline model, we systematically add various materials to the silicate model. We find that even with exceptionally high fractions (up to 100%) of substances such as polycyclic aromatic hydrocarbons, ethanol, ethane, and tholin (I. Cherchneff et al. 1991; H. Imanaka et al. 2004; T. L. Myers et al. 2018), the surface brightness ratio between F300M and F360M consistently remained within a narrow range of 0.8–0.9. Additional exploration of compounds like benzene, toluene, acetone, and IS mix-like materials (D. M. Hudgins et al. 1993; K. Moutzouris et al. 2014) indicates their contribution to be minimal, as we found these components are present at such negligible levels that they fall below meaningful detection thresholds for the purposes of this study. Despite changes in the specific molecules under examination, this consistency strongly suggests that the base composition of silicates and water ice, rather than the varied materials introduced, is primarily responsible for the observed absorption. Thus, water ice potentially remains the dominant factor in the disk’s spectral behavior.

## 6.1. Discussion

Our measurements provide an estimate of the relative scattering efficiency between F300M to F360M from roughly 50 au out to almost 500 au. These measurements reveal decreasing F300M absorption relative to F360M as a function of projected separation. The ratio that we measure for the innermost regions appears to be consistent with those from J. K. Kueny et al. (2024). In their work, they measure the ratio of the disk flux with respect to the star with a narrowband 3.1 µm filter and a broadband 3.8 µm (L') filter. The ratio between these two values is ~0.44, roughly consistent with our values of ~0.4 inside of 50 au (Figure 8), especially given the difference in filters used. Their observations provided a clearer picture of the inner ring in the system, undistorted by the NIRCam coronagraph IWA, whereas our new observations detect the B and C rings at these wavelengths for the first time.

In our modeling, we find that we can explain this absorption by varying the mass fraction of water ice and that in the inner regions of the disk (~100 au), it could be anywhere from 20%–100% depending on the minimum grain size (beyond the blow-out size). At the outer reaches of the system, closer to ~400 au, the mass fraction of water ice could range from ~1% to ~20%, again, depending on the minimum grain size.

Recent JWST spatially resolved, NIR scattered light spectroscopy has detected water ice scattering features in the outer edge of the debris ring and in the extended halo around the ~20 Myr old star HD 181327 (C. Xie et al. 2025). Indeed, the broad 3.0 µm water ice feature shows not only a wide bowl shape but also a narrow 3.1 µm Fresnel peak, similar to that observed toward trans-Neptunian objects (N. Pinilla-Alonso et al. 2025), believed to be created by large crystalline water ice grains. Furthermore, the IFU data reveal a gradient in the strength of the water ice feature as a function of stlocentric distance with estimated water ice fractions of 0.1%, 7.5%, and 21.0% at 80–90 au, 90–105 au, and 105–120 au, respectively, similar to that estimated for HD 141569. For HD 181327, the smaller water ice fraction in the inner disk may be the result of
efficient photodesorption of micron-sized water ice grain by stellar ultraviolet photons.

In AB Aur, using similar filters to ours, S. K. Betti et al. (2022) found mass fractions of up to 68% water ice, which they argued is likely unphysical. In the inner regions of HD 141569A, the measurements are consistent with values this high and potentially even higher. Assuming that these are also unphysical, this suggests that, at least in the inner regions of the disk, the minimum grain size is on the smaller end of our range. S. K. Betti et al. (2022) suggested that one mechanism to explain their measurements is that the surface probed by their longer wavelength measurements (L') is deeper into the disk due to a reduced opacity at longer wavelengths. In their case, this argument is supported by the difference in disk morphology they see between their L' observations and those from shorter wavelengths. In the case of HD 141569A, the disk morphology appears consistent between our observations and visible wavelengths (e.g., Figure 4). This may push the disk to the smaller dust grain regime in Figure 10, with minimum grains more in the 2–3 µm range. However, we note that our calculations were carried out assuming an optically thin disk.

In both HD 100546 (M. Honda et al. 2016) and AB Aur (S. K. Betti et al. 2022), mild trends of increasing water ice absorption with radial distance are seen. The stability of the JWST’s space-based platform and sensitivity at these wavelengths allows us to detect a robust trend with radius. However, here we measure a decreasing trend in the ice mass fraction with stellocentric distance. There are several potential explanations for this trend. One explanation could be that at larger separations, the disk is more flared and so the outer disk sees more UV flux from the host star (assuming the disk is at least marginally optically thick in the radial direction near the midplane), increasing the rate of photodissociation. Another explanation, as suggested above, is that other volatiles may freeze out onto the dust grains at larger separations, reducing the mass fraction of ice. Finally, J. A. White & A. C. Boley (2018) suggested that the outer regions of this disk (>95 au) may have a steeper/more evolved grain size distribution than the inner disk. This difference in the outer regions as compared to the inner regions could perhaps be due to the residual gas in the inner system. If this is the case, then it could be that (a) the process of collisional erosion has removed some of the water ice in these outer regions (A. Bonsor et al. 2023), (b) a difference in dust grain distributions is affecting the ratio between the JWST filters inside and outside of the ice feature (M. Kim et al. 2024), or (c) a combination of both effects.

The treatment of water ice as distinct inclusions mixed with silicates using the Bruggeman rule carried out here may not fully capture the complexity of real grain compositions. Recent laboratory work by A. Potapov et al. (2021) has demonstrated that physically mixed silicate/ice samples can show different optical properties than mathematically mixed optical constants. However, for our specific two-filter analysis at 3.0 and 3.6 µm, such differences appear relatively modest at the cold temperatures relevant to our disk environment. Nevertheless, comprehensive modeling—including systematic comparison of hydrated versus mixed silicates, crystalline versus amorphous ice phases, and their effects on the F300M/F360M brightness ratio—represents an interesting extension of this work that could be explored in dedicated follow-up work that may include multiwavelength observations.

## 7. Sensitivity to Planets

In order to search for directly imaged planets in our data set, we began by subtracting the M-star models (Section 3.2) and the disk models (Section 5) from the individual integrations of each roll for both the F300M and F360M data sets. The rest of the processing in this section was completed using *space-KLIP*. All of the integrations were cropped to a 10" × 10" box, in order to reduce the influence of the M-star residuals on subsequent processing. The integrations for a given filter were then all mutually aligned with each other via cross correlation (the ‘fourier’ method in *spaceKLIP*), before being median-combined for a given roll angle or dither position for the science and reference data, respectively. Although when subtracting the host PSF to image the disk we did not find it necessary to align the images, at this stage, we found it provided better results.

For each filter, KLIP PSF subtraction (R. Soummer et al. 2012) was performed with *spaceKLIP* (via *pyklip*; J. J. Wang et al. 2015) using the full data set. Before carrying out PSF subtraction, all of the collapsed images were highpass filtered with an FWHM of 3 and 4 pixels, for F300M and F360M, respectively (equal to ~twice the Airy FWHM at the filter center wavelengths), in order to minimize the impact of any residual extended disk emission. We attempted ADI, RDI, and a combination of both (ADI+RDI). For each method, we tried varying the number of annuli (1, 5, and 20) and the number of subsections (one, four, and six) over which the PCA algorithm is run independently, and the number of KL modes (1, 5, and 10). Ultimately, we found that using 20 annuli, four subsections, and 10 KL modes provided the best speckle suppression. Overall, RDI performed the worst of the three methods, and the ADI and ADI+RDI sensitivity curves were indistinguishable by eye. As a result, hereafter, we only considered the ADI+RDI and RDI measurements. This result is different from that of A. L. Carter et al. (2023), who found that RDI performed better than ADI (and that ADI+RDI was comparable to RDI) when imaging the HIP 65426 system. This may be due to the significant disk residuals in the RDI-only reduction, which become self-subtracted out when using ADI.

Raw contrast limits were calculated as a function of separation and then calibrated, following the methods described in A. L. Carter et al. (2023). The calibrated 5σ contrast curves were converted to planet masses using the ATMO2020 model grids (M. W. Phillips et al. 2020), with an age of 5 Myr and a distance of 111 pc, where the filter magnitudes for a given planet mass were calculated using the *species* Python package (T. Stolker et al. 2020).

No statistically significant point sources were detected. Figure 11 displays the final PSF-subtracted images for each filter, as well as the 5σ planet mass sensitivity curve. The F300M filter is more sensitive to lower planet masses, reaching a floor of ~2 *M*<sub>Jup</sub> outside of about 1". The F360M filter only reaches a sensitivity of 2 *M*<sub>Jup</sub> outside of 2", after which both filters' sensitivities remain relatively flat until they both turn up slightly around 4" (likely due to the influence of the M-star residuals).

The PSF-subtracted images show clear residuals from both the disk model subtraction and the M-star model subtraction, limiting the achievable point-source contrast. To evaluate the impact of these residuals, we also include in Figure 11 simulated contrast curves from *pyNRC*<sup>64</sup> (J. Leisenring 2025). In *pyNRC*, we simulated the exact instrument setup for both

<footnote n="64"><sup>64</sup> https://pynrc.readthedocs.io/</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Left and center panels: the final PSF-subtracted images for both filters. No point sources were detected. Subtracting off the M-star and disk models leaves significant residuals that affect the contrast floor. Right panel: the mass sensitivity of our observations, as measured from the data and simulated with *pyNRC*.</figcaption></figure>

filters, including customized PSF models for both the HD 141569A and HD 140986, the detector readout parameters, and the difference in roll angle between the two science images. We included realistic pointing offsets, and 2 and 5 nm wave front error drifts between subsequent rolls, and between the science and reference star, respectively (values that produce good agreement with on-sky results with HIP 65426 b; Leisenring, private communication). The simulated contrast curves demonstrate that in the absence of residuals and photon noise from the disk and M-stars, our observations should be sensitive to planets down to ~0.5*M<sub>Jup</sub>*, and that the F360M filter should be more sensitive, in contrast to our data. The superior performance of the F300M filter here may be due to the improved performance of our model at F300M relative to F360M (Figure 6), or possibly the lower disk brightness due to the ice absorption. We note that we did not attempt to correct for any potential extinction of a planet due to the disk’s dust, as this was beyond the scope of this work. With the current limits, we cannot, for example, rule out the 0.2–2*M<sub>Jup</sub>* planet predicted in the 235–250 au region by M. C. Wyatt (2005). These observations are also not sensitive to the kinds of planets predicted to be stirring and sculpting the cold debris disks found around slightly older stars, which are typically Neptune- to Jupiter-mass and orbiting in the 10–100 au region (T. D. Pearce et al. 2022). We note however, that the models from T. D. Pearce et al. (2022) do not consider the possibility of significant amounts of gas in the system (as is the case here), which means that there is not a straightforward correspondence between the distribution of small dust and that of planetesimals.

Despite the general expectation of improved mass sensitiv-
ity of JWST relative to ground-based observations at
separations >1", we find that in this scenario, ground-based
observations of HD 141569A are able to achieve deeper
sensitivities. In particular, planet search results from C. Perrot
et al. (2016) using VLT/SPHERE/IRDIS in the H band are
able to reach sensitivities down to ~1*M<sub>Jup</sub> outside of ~1".
Similar sensitivities were reached by M. Janson et al. (2013)
using Subaru/HiCIAO outside of the speckle limited regime.
Our poorer limits are likely due to the disk residuals in our
data, which may be affected by the larger PSF size and the
complicated PSF morphology of the NIRCam coronagraph.
Furthermore, in Figure 11, the M-star residuals can be seen to
cover nearly half the field in both filters, an effect not present
in the SPHERE data. JWST observations of the same system
might also be able to achieve lower mass sensitivities with a
filter better suited for a planet search, such as F444W, rather

than filters optimized for water ice detection. Indeed, the
results from A. L. Carter et al. (2023) demonstrate very good
agreement with simulated performance and suggest that the
disagreement here is related to the disk and the noise
contributions from the nearby M-stars.

8. Conclusion

Here we have presented a first analysis of new JWST
images of the HD 141569A disk, obtained in the F300M and
F360M filters. While hints of the disk can be seen in the raw
images, subtraction of the host star PSF and of the neighboring
M-stars is required in order to clearly isolate the disk features.
The M-stars were subtracted by fitting an instrument point-
source model (i.e., a *WebbPSF* model) to the data. We carried
out an MCRDI procedure in order to simultaneously constrain
the disk morphology and subtract the host star PSF. The PSF-
subtracted disk images are reminiscent of previous imaging
(e.g., with HST/STIS) but with subtle differences. A more
thorough comparison and broad multiwavelength analysis has
been left for future work.

The bright PSF wings of the M-star companions contam-
inate the disk emission in these images more than in previous
imaging, which may be due to a combination of the new
wavelength range and the unique JWST PSF structure. When
planning future observations with similar nearby stars, we
recommend observers be conscious of this effect and should
consider forward modeling their effects using software such as
*pyNRC*.

We constructed a disk model composed of four rings to fit
the data in both filters. The model matched the data reasonably
well and was able to recover the complicated inner structure
caused by a combination of the NIRCam off-axis PSF
morphology variations and the coronagraph mask throughput.
As with the extended M-star PSF wings, this behavior may not
be intuitive when compared with previous observations, and
we recommend that future observers planning to observe
bright disks near the coronagraph mask consider forward
modeling the expected behavior with *pyNRC* or similar
software.

By comparing our forward-modeled disk images in both
filters, we recovered a radial profile of the scattering efficiency
ratio between the F300M and F360M filters, which varied
from 0.4 up to ~0.8. A similar profile was extracted from a
deconvolved disk image based on the same PSF-subtracted
data. We demonstrated that this range of values could be
explained by changing water content by a factor of 3–10 and/
or changing minimum grain size by a factor of ~5. The exact values are also degenerate with porosity. However, unless the porosity is strictly 0%, it is hard to explain the data with only a radial change in minimum grain size. Future observations using either NIRSpec or NIRCam could be used to characterize the water ice around HD 141569. NIRSpec spectroscopy could be used to measure the shape of the 3.0 µm water ice feature and thereby break the degeneracy in the grain properties for the brightest regions of the disk at distances >1". NIRCam Coronagraphy could be used to examine the water ice content both interior to the NIRSpec Inner Working Angle and in the faint outer disk where high-SNR spectroscopy will be more difficult.

Our current two-filter analysis demonstrates that water ice features can plausibly explain the observed radial brightness variations. We recommend that future NIRCam imaging studies consider incorporating: (1) three or more wavelengths spanning the 3 m feature, (2) systematic comparison of different ice phases and silicate hydration states, and (3) integration with SED, polarization, and MIR data (e.g., 4.5, 6–7, and 13 µm broad water features) to break the degeneracies identified here.

We did not detect any point-source companions in the system. Our sensitivity analysis revealed that we were less sensitive to planets than previous ground-based observations, likely due to residual disk signal in our data and/or the residual PSF wings of the neighboring M-stars. Indeed, the observations did not achieve the sensitivity predicted by the pyNRC instrument simulator in the absence of disk and M-star companions. Improved planet sensitivity may be realized in the future with better M-star models and/or more sophisticated disk models that include azimuthal variations, and a more optimal filter for planet searches (e.g., F444W).

This work presents only a first look at this data set, which is of high SNR and ripe for more detailed analysis. Future work may include a broader multiwavelength grain characterization that combines the new data with archival data sets such as the STIS data included herein. Such an analysis could explore the influence of water ice content on the shape of the continuum (e.g., M. Kim et al. 2024), or include a different grain size distribution. Future work may also consider the azimuthal variation in the disk and search for high-order spatial variations in the absorption seen in this work. Additionally, an analysis that includes varying scattering phase functions between the two filters may also reveal a better fit to the data.

Overall, this work demonstrates JWST’s unique ability to search for and characterize the water ice content in circumstellar disks at unprecedented depths. In order to facilitate future analysis of these data, and realize the goals of the ERS program, we have released the analysis code described above here (github.com/maxwellmb/jwst_ers_NIRCAM_HD141569A).

## Acknowledgments

This project was supported by a grant from STScI (JWST-ERS-01386) under NASA contract NAS5-03127. This work benefited from the 2022 Exoplanet Summer Program in the Other Worlds Laboratory (OWL) at the University of California, Santa Cruz, a program funded by the Heising-Simons Foundation. J.M.L. acknowledges support from the JWST NIRCam project under NASA contract NAS5-02105 (M. Rieke, University of Arizona, PI). S.M. was supported by a Royal Society University Research Fellowship (URF-R1-221669). E.C. acknowledges funding from the European Union (ERC, ESCAPE, project No. 101044152). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. O.A. is a Senior Research Associate of the Fonds de la Recherche Scientifique—FNRS. A.Z. acknowledges support from ANID—Millennium Science Initiative Program—Center Code NCN2024_001.

**Facility:** JWST (NIRCam).

**Software:** astropy (Astropy Collaboration et al. 2013, 2018), JAX (J. Bradbury et al. 2018), NumPy (C. R. Harris et al. 2020), SciPy (P. Virtanen et al. 2020), spaceKLIP (J. Kammerer et al. 2022; A. L. Carter et al. 2023), pyNRC (J. Leisenring 2024), cmasher (E. van der Velden 2020), pyKLIP (J. J. Wang et al. 2015), Matplotlib (J. D. Hunter 2007), Vortex Image Processing (C. A. Gomez Gonzalez et al. 2017; V. Christiaens et al. 2023), MCFOST (C. Pinte et al. 2006, 2009).

## Appendix A

Figures 12 and 13 show the residuals (divided by the local noise) and the corner plots from the flux fitting procedures described in Section 5.2, demonstrating overall good behavior: the residuals in the disk images are typically <3σ, and the parameters in the corner plots are all well constrained.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Comparison of F300M and F360M residuals after the flux fitting procedure. The residuals have been divided by the local noise to show the signal-to-noise ratios. The fitting was carried out simultaneously on the individual roll data without being derotated and roll-combined. Only one roll angle is shown here, but the results are similar for both rolls.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Corner plots from the flux fitting procedure for the F300M (left panel) and F360M (right panel) filters.</figcaption></figure>

## Appendix B
Reconvolution of Deconvolved Images

Deconvolution is essentially a distinct method for solving the same problem as forward modeling of synthetic disk images: identifying some idealized image of the circumstellar scene that, when forward-modeled for instrumental effects, provides a good match for the observations. To validate the result of deconvolution, we can therefore simply treat the deconvolved image like a “pre-convolved” synthetic disk model image and apply the same forward modeling treatment (rotating to the appropriate position angles, applying mask

transmission, convolving with the PSF, and then derotating and combining). Applying this to our deconvolved images, we then compute the residuals between the real MCRDI images and the resulting “reconvolved” images. The results of this analysis are shown in Figure 14. Compared with the synthetic model residuals of Figure 6, the reconvolved residuals are generally significantly smaller. This exercise suggests that the implemented deconvolution algorithm is working as intended to produce a final image that is a good approximation of the ground-truth circumstellar scene (at least up to the accuracy of the synthetic PSF models used for convolution).
<figure><img src="image_17.png" /></figure>

**Figure 14.** The same as in Figure 6, but for reconvolution of the deconvolved images following the procedure outlined in Appendix B. Gray contours in the right column indicate ±3σ thresholds for the residuals. To approximate the noise level for this purpose, we adopt the propagated pixel-wise uncertainties from the ERR FITS extensions. As these estimates neglect residual speckle noise, they should be assumed to underestimate the noise where speckle residuals are nonnegligible (at very small separations and near the M dwarfs). The residual maps used for the contours were convolved with an FWHM-diameter tophat filter, which serves to mitigate the appearance of the single-pixel deviations expected from Gaussian noise but does not qualitatively alter the result.

## ORCID iDs

Maxwell A. Millar-Blanchaer https://orcid.org/0000-0001-6205-9233

Élodie Choquet https://orcid.org/0000-0002-9173-0740

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

Sebastián Marino https://orcid.org/0000-0002-5352-2924

Jens Kammerer https://orcid.org/0000-0003-2769-0438

Aarynn L. Carter https://orcid.org/0000-0001-5365-4815

Isabel Rebollido https://orcid.org/0000-0002-4388-6417

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

Minjae Kim https://orcid.org/0000-0001-6218-2004

Paul Kalas https://orcid.org/0000-0002-6221-5360

Karl R Stapelfeldt https://orcid.org/0000-0002-2805-7338

Sasha Hinkley https://orcid.org/0000-0001-8074-2562

Mark Booth https://orcid.org/0000-0001-8568-6336

Carol A. Grady https://orcid.org/0000-0001-5440-1879

Elisabeth C. Matthews https://orcid.org/0000-0003-0593-1560

Beth A. Biller https://orcid.org/0000-0003-4614-7035

Andrew Skemer https://orcid.org/0000-0001-6098-3924

Julien H. Girard https://orcid.org/0000-0001-8627-0404

Schuyler G. Wolff https://orcid.org/0000-0002-9977-8255

Kimberly Ward-Duong https://orcid.org/0000-0002-4479-8291

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

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

Eric Pantin https://orcid.org/0000-0001-6472-2844

Brenda C. Matthews https://orcid.org/0000-0003-3017-9577

Stanimir Metchev https://orcid.org/0000-0003-3050-8203

Marshall D. Perrin https://orcid.org/0000-0002-3191-8151

Christine H. Chen https://orcid.org/0000-0002-8382-0447

Katie Crotts https://orcid.org/0000-0003-4909-256X

Olivier Absil https://orcid.org/0000-0002-4006-6237

William O. Balmer https://orcid.org/0000-0001-6396-8439

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

Gabriele Cugno https://orcid.org/0000-0001-7255-3251

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

Camilla Daniëlski https://orcid.org/0000-0002-3729-2663

Kielan K. W. Hoch https://orcid.org/0000-0002-9803-8255

Markus Janson https://orcid.org/0000-0001-8345-593X

Elena Manjavacas https://orcid.org/0000-0003-0192-6887

Ben J. Sutlieff https://orcid.org/0000-0002-9962-132X

Shrishmoy Ray https://orcid.org/0000-0003-2259-3911

Bin B. Ren https://orcid.org/0000-0003-1698-9696

Emily Rickman https://orcid.org/0000-0003-4203-9715

Genaro Suárez https://orcid.org/0000-0002-2011-4924

Christopher A. Theissen https://orcid.org/0000-0002-9807-5435

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

Jason J. Wang https://orcid.org/0000-0003-0774-6502

Niall Whiteford https://orcid.org/0000-0001-8818-1544

Mark C. Wyatt https://orcid.org/0000-0001-9064-5598

Alice Zurlo https://orcid.org/0000-0002-5903-8316

## References

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

Ardila, D. R., Lubow, S. H., Golimowski, D. A., et al. 2005, ApJ, 627, 986

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

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

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

Benisty, M., Dominik, C., Follette, K., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 605

Betti, S. K., Follette, K., Jorquera, S., et al. 2022, AJ, 163, 145
Biller, B. A., Liu, M. C., Rice, K., et al. 2015, MNRAS, 450, 4446

Blondel, M., Berthet, Q., Cuturi, M., et al. 2022, in Advances in Neural Information Processing Systems, ed. S. Koyejo et al. (Red Hook, NY: Curran Associates Inc),

Bonsor, A., Wyatt, M. C., Marino, S., et al. 2023, MNRAS, 526, 3115

Bradbury, J., Frostig, R., Hawkins, P., et al. 2018 JAX: composable transformations of Python+NumPy programs v0.3.13, http://github.com/jax-ml/jax

Bruzzone, J. S., Metchev, S., Duchêne, G., et al. 2020, AJ, 159, 53

Burns, J. A., Lamy, P. L., & Soter, S. 1979, Icar, 40, 1

Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2023, JWST Calibration Pipeline v1.12.1, Zenodo, doi: 10.5281/zenodo.8380331

Byrd, R. H., Lu, P., Nocedal, J., & Zhu, C. 1995, SJSC, 16, 1190

Carter, A. L., Hinkley, S., Bonavita, M., et al. 2021, MNRAS, 501, 1999

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

Cherchneff, I., Barker, J. R., & Tielens, A. G. G. M. 1991, ApJ, 377, 541

Christiaens, V., Gonzalez, C., Farkas, R., et al. 2023, JOSS, 8, 4774

Clampin, M., Krist, J. E., Ardila, D. R., et al. 2003, AJ, 126, 385

Currie, T., Grady, C. A., Cloutier, R., et al. 2016, ApJL, 819, L26

Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, yCat, II/246

Di Folco, E., Pêricaud, J., Dutrey, A., et al. 2020, A&A, 635, A94

Dohnanyi, J. S. 1969, JGR, 74, 2531

Draine, B. T. 2003, ApJ, 598, 1017

Duchêne, G., McCabe, C., Ghez, A. M., & Macintosh, B. A. 2004, ApJ, 606, 969

Duchêne, G., McCabe, C., Pinte, C., et al. 2010, ApJ, 712, 112

Filacchione, G., Capaccioni, F., Ciarniello, M., et al. 2012, AGUFM, 2012, P53B

Flaherty, K. M., Hughes, A. M., Andrews, S. M., et al. 2016, ApJ, 818, 97

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Gaia Collaboration 2020, yCat, I/350

Gaia Collaboration, Prustyi, T., de Bruijne, J., H., J., et al. 2016, A&A, 595, A1

Gardner, J. P., Mather, J. C., Abbott, R., et al. 2023, PASP, 135, 068001

Garufi, A., Ginski, C., van Holstein, R. G., et al. 2024, A&A, 685, A53

Ginski, C., Garufi, A., Benisty, M., et al. 2024, A&A, 685, A52

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

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

Goudfrooij, P., Grumm, D., Volk, K., & Bushouse, H. 2024, PASP, 136, 014503

Grigorieva, A., Thébault, P., Artymowicz, P., & Brandeker, A. 2007, A&A, 475, 755

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

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

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

Hinkley, S., Carter, A. L., Ray, S., et al. 2022, PASP, 134, 095003

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

Honda, M., Kudo, T., Takatsuki, S., et al. 2016, ApJ, 821, 2

Houk, N., & Swift, C. 1999, MSS, 5, 0

Hudgins, D. M., Sandford, S. A., Allamandola, L. J., & Tielens, A. G. G. M. 1993, ApJS, 86, 713

Hughes, A. M., Duchêne, G., & Matthews, B. C. 2018, ARA&A, 56, 541

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

Imanaka, H., Khare, B. N., Elsila, J. E., et al. 2004, Icar, 168, 344

Inoue, A. K., Honda, M., Nakamoto, T., & Oka, A. 2008, PASJ, 60, 557

Janson, M., Brandt, T. D., Moro-Martín, A., et al. 2013, ApJ, 773, 73

Jaschek, C., & Jaschek, M. 1992, A&AS, 95, 535

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

Kim, M., Kennedy, G. M., & Roccatagliata, V. 2024, MNRAS,

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

Konishi, M., Grady, C. A., Schneider, G., et al. 2016, ApJL, 818, L23

Kóspál, Á., Moór, A., Juhász, A., et al. 2013, ApJ, 776, 77

Krist, J. E., Beichman, C. A., Trauger, J. T., et al. 2007, Proc. SPIE, 6693, 66930H

Krist, J. E., Balasubramanian, K., Muller, R. E., et al. 2010, Proc. SPIE, 7731, 77313J

Kueny, J. K., Weinberger, A. J., Males, J. R., et al. 2024, ApJ, 961, 77

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

Lawson, K., Schlieder, J. E., Leisenring, J. M., et al. 2023, AJ, 166, 150

Lawson, K., Schlieder, J. E., Leisenring, J. M., et al. 2024, ApJL, 967, L8

Leisenring, J. 2024, pyNRC: Python ETC and Simulator for JWST NIRCam v1.2.0, Zenodo, doi:10.5281/zenodo.5829552

Leisenring, J. 2025, pyNRC: Python ETC and Simulator for JWST NIRCam v2.0.1, Zenodo, doi:10.5281/zenodo.5829552

Liu, D. C., & Nocedal, J. 1989, MatPr, 45, 503

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

Lyra, W., & Kuchner, M. 2013, Natur, 499, 184

Marino, S., Cataldi, G., Jankovic, M. R., Matrà, L., & Wyatt, M. C. 2022, MNRAS, 515, 507

Mastrapa, R. M., Sandford, S. A., Roush, T. L., Cruikshank, D. P., & Dalle Ore, C. M. 2009, ApJ, 701, 1347

Mawet, D., Choquet, É., Absil, O., et al. 2017, AJ, 153, 44

Mazoyer, J., Boccaletti, A., Choquet, É., et al. 2016, ApJ, 818, 150

McCabe, C., Duchêne, G., & Ghez, A. M. 2003, ApJL, 588, L113

McCabe, C., Duchêne, G., Pinte, C., et al. 2011, ApJ, 727, 90

Merín, B., Montesinos, B., Eiroa, C., et al. 2004, A&A, 419, 301

Miles, B. E., Biller, B. A., Patapis, P., et al. 2023, ApJL, 946, L6

Miley, J. M., Panić, O., Wyatt, M., & Kennedy, G. M. 2018, A&A, 615, L10

Milli, J., Lagrange, A. M., Mawet, D., et al. 2014, A&A, 566, A91

Min, M., Dullemond, C. P., Kama, M., & Dominik, C. 2011, Icar, 212, 416

Moerchen, M. M., Telesco, C. M., & Packham, C. 2010, ApJ, 723, 1418

Mouillet, D., Lagrange, A. M., Augereau, J. C., & Ménard, F. 2001, A&A, 372, L61

Moutzouris, K., Papamichael, M., Betsis, S. C., et al. 2014, ApPhB, 116, 617

Myers, T. L., Tonkyn, R. G., Danby, T. O., et al. 2018, ApSpe, 72, 535

O'Brien, D. P., Izidoro, A., Jacobson, S. A., Raymond, S. N., & Rubie, D. C. 2018, SSRv, 214, 47

Oka, A., Inoue, A. K., Nakamoto, T., & Honda, M. 2012, ApJ, 747, 138

Pascual, N., Montesinos, B., Meeus, G., et al. 2016, A&A, 586, A6

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

Pêricaud, J., Di Folco, E., Dutrey, A., Guilloteau, S., & Piétu, V. 2017, A&A, 600, A62

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

Perrot, C., Boccaletti, A., Pantin, E., et al. 2016, A&A, 590, L7

Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 323, L49

Petrus, S., Whiteford, N., Patapis, P., et al. 2024, ApJL, 966, L11

Phillips, M. W., Tremblin, P., Baraffe, I., et al. 2020, A&A, 637, A38

Pinilla-Alonso, N., Brunetto, R., De Prá, M. N., et al. 2025, NatAs, 9, 230

Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

Potapov, A., Bouwman, J., Jäger, C., & Henning, T. 2021, NatAs, 5, 78

Quillen, A. C., Varnière, P., Minchev, I., & Frank, A. 2005, AJ, 129, 2481

Ray, S., Sallum, S., Hinkley, S., et al. 2025, ApJL, 983, L25

Rebollido, I., Stark, C. C., Kammerer, J., et al. 2024, AJ, 167, 69

Ren, B., Pueyo, L., Zhu, G. B., Debes, J., & Duchêne, G. 2018, ApJ, 852, 104

Rich, E. A., Monnier, J. D., Aarnio, A., et al. 2022, AJ, 164, 109

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

Richert, A. J. W., Lyra, W., & Kuchner, M. J. 2018, ApJ, 856, 41

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

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

Sallum, S., Ray, S., Kammerer, J., et al. 2024, ApJL, 963, L2

Singh, G., Bhowmik, T., Boccaletti, A., et al. 2021, A&A, 653, A79

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

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

STScI Development Team, 2013 pysynphot: Synthetic Photometry Software Package, Astrophysics Source Code Library, ascl:1303.023

Terada, H., Tokunaga, A. T., Kobayashi, N., et al. 2007, ApJ, 667, 303

Thi, W. F., Pinte, C., Pantin, E., et al. 2014, A&A, 561, A50

Valegård, P. G., Ginski, C., Derkink, A., et al. 2024, A&A, 685, A54

van der Velden, E. 2020, JOSS, 5, 204

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

Wang, J. J., Ruffio, J.-B., De Rosa, R. J., et al. 2015 pyKLIP: PSF Subtraction for Exoplanets and Disks, Astrophysics Source Code Library, ascl:1506.001

Weinberger, A. J., Becklin, E. E., Schneider, G., et al. 1999, ApJL, 525, L53

Weinberger, A. J., Rich, R. M., Becklin, E. E., Zuckerman, B., & Matthews, K. 2000, ApJ, 544, 937

Westall, F., & Brack, A. 2018, SSRv, 214, 50

White, J. A., & Boley, A. C. 2018, ApJ, 859, 103

White, J. A., Boley, A. C., Hughes, A. M., et al. 2016, ApJ, 829, 6

Wyatt, M. C. 2005, A&A, 440, 937

Wyatt, M. C., Panić, O., Kennedy, G. M., & Matrà, L. 2015, Ap&SS, 357, 103

Xie, C., Chen, C. H., Lisse, C. M., et al. 2025, Natur, 641, 608

Xie, C., Choquet, E., Vigan, A., et al. 2022, A&A, 666, A32

Zhu, C., Byrd, R., Lu, P., & Nocedal, J. 1997, ACM Transactions on Mathematical Software, 23, 550

Zuckerman, B., Forveille, T., & Kastner, J. H. 1995, Natur, 373, 494