Probing the Distinct Extinction Law of the Pillars of Creation in M16 with JWST

Investigating the extinction law in regions of high dust extinction, such as the Pillars of Creation within the M16 region, is crucial for understanding the densest parts of the interstellar medium (ISM). In this study, we utilize observations from the Near-Infrared Camera and the Mid-Infrared Instrument onboard the James Webb Space Telescope (JWST) to analyze the color-excess ratios E(F090W − λ)/E(F090W − F200W) across a wavelength range of 0.9–7.7 μm. Our method involves performing linear regression on color–color diagrams to derive these ratios. The enhanced detection capabilities of JWST data allow us to probe the distinct extinction law to the densest regions in M16 corresponding to an extinction depth up to A V ∼ 60 mag. Remarkably, the resultant color-excess ratio curve exhibits a flatter profile than predicted by typical dust extinction models with R V = 5.5 for dense ISM environments. Moreover, we observe that the mid-infrared extinction law diverges from the near-infrared power law, showing a tendency for the slope to flatten as the wavelength increases. These findings have significant implications for our understanding of the dust properties in dense interstellar environments.


INTRODUCTION
The extinction law characterizes the wavelength-dependent absorption and scattering of light by interstellar dust grains.This critical astrophysical parameter reveals the composition, size distribution, and structural properties of dust within the interstellar medium (ISM) (Draine 2003).Dense regions within molecular clouds, where significant reddening occurs, are pivotal areas of study as they are the cradles of star and planet formation.The characteristics of the dust in these regions significantly impact the initial conditions for star formation and the evolution of protoplanetary disks.Thus, understanding the extinction law in these environments, although challenging, is essential.
Extinction laws have been extensively studied in regions of the ISM with lower density, where extinction effects are less pronounced.The Galactic extinction curve at ultraviolet (UV) and optical wavelengths, which is often characterized in terms of the ratio R V = A V /E(B − V ), has revealed non-universal dust properties that differ across various sightlines (Cardelli et al. 1989;Fitzpatrick et al. 2019;Massa et al. 2020).In contrast, the near-infrared (NIR) extinction law, particularly from 1.2 to 2.2 µm, was previously considered uniform across the Galaxy, typically modeled by a power-law A λ ∝ λ −α with α values around 1.6-1.8(Rieke & Lebofsky 1985;Cardelli et al. 1989).However, more recent studies have reported α values ranging from 2 to 2.6 (Stead & Hoare 2009;Alonso-García et al. 2017;Maíz Apellániz et al. 2020), and suggested a broader applicability of this power-law form across a wider wavelength range (0.8-5.0 µm) (Decleir et al. 2022;Wang & Chen 2024).Furthermore, variations in the NIR power-law index both with wavelength and across different sightlines have been observed, challenging the previous assumption of a uniform NIR Li et al. extinction law (Alonso-García et al. 2017;Nogueras-Lara et al. 2019).In the MIR range (3-8 µm), earlier studies posited that the NIR power-law could extend into the MIR wavelengths (Rieke & Lebofsky 1985;Martin & Whittet 1990).However, later observations, especially those from the Spitzer Space Telescope, have depicted a mostly flat extinction curve across these wavelengths in various interstellar environments (Indebetouw et al. 2005;Gao et al. 2009;Chen et al. 2013;Xue et al. 2016).This suggests that the MIR extinction law might significantly differ from the NIR laws, although observational challenges have historically limited a comprehensive understanding in this regime.
Notably, the extinction depths investigated in most previous studies are not particularly high.However, in dense interstellar regions, variations in the extinction law are expected due to factors such as dust grain growth and ice-mantle formation, as suggested by both theoretical models (Ossenkopf & Henning 1994;Ormel et al. 2011) and observations (Juvela et al. 2015;Li et al. 2024).Observational studies indicate that in regions of significant star formation, where extinction is higher, the extinction law tends to flatten in comparison to more diffuse regions (Moore et al. 2005;Flaherty et al. 2007;McClure 2009).This flattening aligns well with the Weingartner & Draine (2001, hereafter WD01) extinction model, which posits an R V = 5.5 for dense ISM environments.Additionally, there is accumulating evidence that suggests a trend toward a flatter IR extinction law with increasing extinction depth (Chapman et al. 2009;Cambrésy et al. 2011;Li & Chen 2023).
Characterizing the IR extinction law in dense environments presents significant challenges due to the high column densities and deep extinction that often obscure background stars.However, the advent of the James Webb Space Telescope (JWST) with its advanced instruments, Near Infrared Camera (NIRCam; Rieke et al. 2005Rieke et al. , 2023) ) and Mid-Infrared Instrument (MIRI; Rieke et al. 2015), offers new opportunities.These instruments provide unprecedented sensitivity and spatial resolution, enabling a detailed examination of the extinction law across a broad range of IR wavelengths.In this paper, we focus on the IR extinction law in the Pillars of Creation within the Eagle Nebula (M16), a well-known star-forming region located at a distance of 1.74 ± 0.13 kpc (Kuhn et al. 2019).
The structure of this paper is organized as follows: Section 2 details the observations and the methodology of data reduction.Section 3 discusses the results, focusing on the measurement and implications of the IR extinction law in M16.Section 4 provides a summary of the findings and conclusions.

OBSERVATIONS AND DATA REDUCTION
We utilize publicly available observations of M16 captured using both NIRCam and MIRI on board the JWST, as part of the Director's Discretionary (DD) program (ID: 2739; PI: Klaus Pontoppidan).The observations employed the NIRCam filters F 090W , F 187N , F 200W , F 335M , and F 444W , along with the MIRI filters F 770W , F 1130W , and F 1500W .These observations spanned the main structure of the Pillars of Creation in M16.The fields-of-view (FoVs) are 7.4×4.4arcmin2 for NIRCam and 4.3×3.8arcmin 2 for MIRI, respectively.
We have retrieved the Level-2b data products from the Mikulski Archive for Space Telescopes (MAST)1 at the Space Telescope Science Institute for data reduction.The specific data analyzed can be accessed via the MAST DOI 10.17909/3w1e-qp71.The NIRCam images were processed using JWST pipeline version 1.10.1 with Calibration Reference Data System (CRDS) version 11.16.22 and context jwst 1100.pmap.MIRI images were processed using JWST pipeline version 1.13.3 with CRDS version 11.17.14 and context jwst 1216.pmap.Initially, the Level-2b *i2d.fitsfiles were corrected for 1/f noise using the Python routine image1overf.py 2 (Willott 2022).Subsequently, we use the F 200W Level-3 drizzled i2d.fits file as the astrometric reference to extract point sources.The astrometric alignment is applied to all Level-2b *i2d.fitsfiles using the JWST/Hubble Alignment Tool (JHAT3 ; Rest 2023).
For photometric analysis, we have performed PSF photometry on the NIRCam and MIRI data using the starbugII photometric tool (Nally 2023), which is optimized for JWST observations in crowded fields with complex backgrounds.Source detection on the aligned Level-2b images is executed using starbug2 --detect, with a detection threshold of 5σ for the NIRCam F 090W , F 187N , and F 200W sources, and 3σ for the NIRCam F 335M , F 444W , and MIRI sources.We have applied criteria related to geometric parameters such as sharpness and roundness to filter out contaminants including cosmic rays, dust structures, and galaxies.The parameters used for starbugII detection and photometry are summarized in Table 1.
Aperture photometry is also performed on all NIRCam and MIRI images to facilitate source detection and zero-point correction.For NIRCam, we have employed an aperture radius of 1.5 pixels with a sky annulus extending from 3.0 to 4.5 1.1 0.9 0.9 0.9 0.9 0.9 pixels.In contrast, for MIRI, aperture radii are set at 2.5 pixels for the F770W and F1130W bands and 3 pixels for the F1500W band, with background annuli ranging from 4 to 5.5 pixels and 4.5 to 6 pixels respectively.Aperture corrections are derived from jwst nircam apcorr 0004.fits for NIRCam and jwst miri apcorr 0005.fits for MIRI.Background emission in the images is modeled using starbug2 --background to produce clean images suitable for further analysis.Subsequently, PSF photometry is applied to these background-subtracted images using starbug2 --psf, in conjunction with WEBBPSF (Perrin et al. 2014) version 1.1.1.Instrumental zero points are estimated using starbug2 --calc-instr-zp, which involves comparing PSF-based magnitudes with those obtained from reliable aperture photometry.The PSF photometry results are initially calibrated to the AB magnitude system using these zero points.These are then converted to the Vega magnitude system using the offsets provided in jwst nircam abvegaoffset 0001.asdf for NIRCam and jwst miri abvegaoffset 0001.asdf for MIRI.Subsequently, individual catalogs are merged using the routine starbug2-match.

RESULTS AND DISCUSSION
Fig. 1 presents color-magnitude diagrams (CMDs) for the JWST NIRCam and MIRI data, plotting various filters (λ) against the color index F 090W − F 200W .These filters include F 090W , F 187N , F 200W , F 335M , F 444W , and F 770W .The CMDs are based on the final merged photometric catalog.Photometric uncertainties are typically below 0.2 mag for NIRCam bands, while MIRI bands exhibit higher uncertainties due to their relatively lower signal-to-noise ratio (S/N).The diagrams distinctly show two prominent groups of stars, characterized by their varying color indices.The first group, comprising highly reddened young stars embedded within molecular clouds, exhibits a color index F 090W − F 200W exceeding 5 mag.The second group consists of foreground stars, with a color index approximately 2.5 mag.Notably, the diagrams indicate a larger population of sources detected in the NIRCam bands compared to the MIRI bands.

Analysis of Color-Color Diagrams and Reddening Slopes
In the NIR and MIR bands, color-color diagrams are pivotal in understanding the extinction law.In these diagrams, stars mainly share similar intrinsic colors and they are often aligned along a path defined by the reddening vector which quantifies the displacement due to dust extinction.In areas with heavy dust obscuration, the natural color variations among stars are minor relative to the color excess from dust.This significant difference means that the variability in intrinsic colors can generally be overlooked, streamlining the analysis.This method is especially useful in regions of high extinction, where the effects of reddening overshadow any other changes in stellar colors (Román-Zúñiga et al. 2007;Ascenso et al. 2013).When evaluating two color indices, m λ1 − m λ2 and m λ2 − m λ3 , the ratio of color excesses, , can be derived through a linear regression of the stellar distribution on the diagram.The color excess, denoted as E(m λ1 − m λ2 ), represents the deviation of the observed color from its intrinsic value, and is defined by:  The effectiveness of these diagrams in delineating reddening slopes is influenced by the observational constraints across different filters.The observations at longer wavelengths can detect the stars with larger dust column density thanks to the much smaller extinction at the longer wavelengths.This results in variations in the depth of extinction that can be probed in different colors, as well as the sensitivity of these colors to extinction effects.In the subsequent sections, we will explore multi-band color-color diagrams to assess the consistency and variability of color excess ratios across different wavelengths, thereby enhancing our understanding of the interstellar dust in various environments.
We first select F 090W − F 200W as the reference color to construct color-color diagrams for λ = F 187N, F 335M, F 444W, and F 770W (Fig. 2).Photometric uncertainties are constrained to less than 0.5 mag for each band.The MIRI bands F 1130W and F 1500W are excluded due to the limited number of stellar detections, which precluded a reliable determination of color excess ratios.
A linear regression is applied to the color-color diagrams using the LTS LINEFIT method, which accommodates errors in both color axes and intrinsic scatter within the diagram.This robust fitting technique efficiently identifies and excludes significant outliers (gray points in Fig. 2), focusing the fit on the retained black points.The computed color excess ratios, Table 2 also presents recent findings on extinction laws derived from observations by the JWST.Wang & Chen (2024) utilized NIRSpec spectroscopic data to explore the NIR extinction law in the young massive star cluster Westerlund 2. Concurrently, Fahrion & De Marchi (2023) analyzed NIRCam photometric data to probe the extinction law within the star-forming region 30 Doradus in the Large Magellanic Cloud (LMC).Our study reports color excess ratios β F 335M and β F 444W that are lower than those found in these studies, indicative of a comparatively flatter extinction law.
Fig. 3 illustrates the derived color excess ratios, β λ , plotted as a function of wavelength.Unlike previous studies that convert these ratios into relative extinction ratios, such as A λ /A K , which depend on uncertain NIR extinction ratios like A H /A K , we directly present β λ .This approach minimizes errors associated with uncertainties in A H /A K .For comparison, Fig. 3 includes reddening laws from from some important literatures.Our observed β λ curve is significantly flatter than the WD01 model for a diffuse environment at R V = 3.1, and even flatter than the WD01 model at R V = 5.5 for a dense environment, indicating a distinct extinction law in M16.Our results are more consistent with the extinction curve in diffuse ISM toward Cyg OB2-12 reported by Hensley & Draine (2020) and the recent "Astrodust" model by Hensley & Draine (2023).Specifically, the observed β λ is only slightly flatter than that reported by Hensley & Draine (2020), and slightly steeper than that of Hensley & Draine (2023).However, it is important to note that Cyg OB2-12 represents only a single case along a sightline through diffuse ISM, while M16 is a case for a dense region.More cases are required to confirm these results, and currently, we are conducting a systematic analysis of a sample of dense dark clouds to further substantiate our findings.
The flatter extinction curves observed in our study suggest variations in dust grain size and composition within the dense regions of M16.WD01 suggest that a flatter infrared reddening curve correlates with larger grain size distributions.Additionally, the formation of ice mantles on dust grains, as discussed by Ormel et al. (2011)   extinction effects in broadband filters, as noted in studies by Wang et al. (2013) and Xue et al. (2016).In our data, the F 335M filter marginally covers the H 2 O 3.0 µm absorption feature, and the F 444W band seems unaffected by the CO 2 4.27 µm absorption, indicating that other factors may contribute to the observed flat extinction curve.This flat curve implies either a predominance of larger dust grains or a higher proportion of such grains in the dense M16 environment.This evidence reinforces the hypothesis that dust grain characteristics significantly influence the observed NIR extinction laws in various interstellar environments.
The M16 cloud exhibits significant dust obscuration, rendering the F 090W band largely ineffective for penetrating the dense regions.By contrast, the color indices F 200W − F 335M and F 335M − F 444W demonstrate enhanced capability in detecting the more embedded background stars, thanks to their reduced sensitivity to extinction compared to F 090W − F 200W .On the other hand, this leads to broader dispersions in color-color diagrams when using the F 200W − F 335M and F 335M − F 444W colors.To mitigate this dispersion, we have adopted the method outlined by Cambrésy et al. (2011), setting a pixel size of 4 ′′ and applying a 3σ clipping algorithm to compute the average color within each pixel.Subsequently, we have generated color maps for are shown in Fig. 4.These diagrams are constructed using colors averaged from the maps, rather than relying on individual stellar observations.
As depicted in Fig. 4(a) and Fig. 2, the color index F 090W − F 200W spans a range of 0 to 9 mag.Considering a K0III star with an effective temperature of 4750 K, its intrinsic color index (F 090W − F 200W ) 0 is approximately 1.2 mag (Castelli & Kurucz 2003).Utilizing this intrinsic color and applying the WD01 R V = 5.5 model, the extinction depth A V for an embedded background star is calculated as ∼ 23 mag.This is substantial yet typical for such dense regions.A linear regression of F 200W − F 335M vs. F 090W − F 335M indicates a color excess ratio E(F 200W − Using the WD01 model with R V = 5.5 for a K0III giant star, the maximum depth of extinction reaches approximately A V ∼ 60 mag.The observed color excess ratio curve, as depicted in Fig. 3, is flatter than that predicted by this model, suggesting that the actual depth could be even higher.At this depth, the derived color excess ratio is E(F 335M − F 444W )/E(F 200W − F 335M ) ≈ 0.28.In contrast, using F 090W − F 200W as the reference, the color excess ratio is computed as (β F 444W − β F 335M ) / (β F 335M − 1) ≈ 0.347, which is higher than the 0.28 obtained from direct linear fitting of F 335M − F 444W versus F 200W − F 335M .This systematic difference suggests that discrepancies may arise due to changes in dust properties with increased extinction depth, potentially indicating larger dust grain sizes.Additional independent observational evidence for the growth of dust grains in dense interstellar environments comes from the observations of near-IR core shine (Jones et al. 2016), and of strongly forward-directed scattering phase functions, combined with a high long-wavelength albedo of dust in globules (Witt et al. 1990;Togi et al. 2017).

Variability of the power-law index α in the NIR extinction Law
Numerous studies have demonstrated that the NIR extinction law can be characterized as a power-law form A λ ∝ λ −α .Pioneering work by Naoi et al. (2006Naoi et al. ( , 2007) ) highlighted that the NIR extinction slope tends to flatten with increasing extinction depth.Further, Nogueras-Lara et al. (2019) documented variations in α between the NIR JH and HK bands, with α JH consistently exceeding α HK .These observations suggest a decrease in the extinction law slope at longer wavelengths.Recent spectroscopic studies extend these findings across the 0.8 to 5.0 µm range, confirming the power-law behavior with varying indices (Decleir et al. 2022;Wang & Chen 2024).
Despite the limitations of photometric data, our analysis leverages multi-wavelength color excess ratios to explore the α variability.Assuming a power-law model for NIR extinction, we derive the relationship between the color excess ratios and α as follows: Here, λ 1 , λ 2 , and λ 3 denote the effective wavelengths of the JWST NIRCam bands, ordered as λ 1 < λ 2 < λ 3 .lower than α F 090W ∼F 335M = 2.44.These findings affirm the non-uniformity of α across the spectrum from 0.9 µm to 4.4 µm, underscoring a trend toward flatter extinction curves at longer wavelengths.Moreover, in the MIR wavelengths, compared to the NIR, we are able to probe denser regions of molecular clouds that exhibit higher extinction values.The observed lower values of α in the MIR, indicative of flatter extinction curves relative to those in the NIR, further support our previous observations that denser regions within star-forming environments tend to exhibit flatter extinction curves.

SUMMARY
This study utilizes the sophisticated imaging capabilities of the JWST NIRCam and MIRI to conduct deep observations of the Pillars of Creation within the M16 star-forming region.Our analysis has yielded the color-excess ratio curve across a wavelength range of 0.9 to 8 µm.By conducting linear fits on the color-color diagrams, we are able to directly derive the extinction law in M16.The enhanced detection capabilities of JWST data allow us to probe extinction depths up to A V ≈ 60 mag.The derived color-excess ratio curve for the dense region in M16 exhibits a flatter profile than those typically obtained from previous works, and is even less steep than the WD01 R V = 5.5 extinction curve.
Furthermore, our findings reveal that the extinction law in the MIR range from 2 to 5 µm does not simply continue the NIR power-law relationship.Instead, we observed a decrease in the slope of the IR extinction law with increasing wavelength, suggesting a modification in dust characteristics or scattering processes in the MIR spectrum.
We would like to thank the anonymous referee for the very helpful comments and suggestions.This work is supported by the National Key R&D program of China 2022YFA1603102 and 2019YFA0405500, the National Natural Science Foundation of China 12133002, 12173034, and 12322304.B.Q.C. acknowledges the National Natural Science Foundation of Yunnan Province 202301AV070002 and the Xingdian talent support program of Yunnan Province.X.C. thanks to Guangdong Province Universities and Colleges Pearl River Scholar Funded Scheme (2019).

Figure 1 .
Figure 1.Color-magnitude diagrams for all identified sources within the Pillars of Creation in M16, depicted as black points.Each panel displays data obtained with a different filter λ plotted against the color index F 090W − F 200W , where λ includes F 090W , F 187N , F 200W , F 335M , F 444W , and F 770W .The total number of sources in each filter is indicated in the upper right corner of each panel.

Figure 2 .
Figure 2. Color-color diagrams of F 090W − λ versus F 090W − F 200W for all sources within the Pillars of Creation in M16.Panels display λ=F 187N (top left), F 335M (top right), F 444W (lower left), and F 770W (lower right).The red solid lines represent the best-fit linear models to the color indices for each band using LTS LINEFIT.The 3σ confidence intervals of these fits are indicated by red dashed lines.Outliers are marked as grey points.Best-fit slopes are annotated in the upper left corner of each panel.
and McClure et al. (2023), could further influence the extinction law characteristics.Recent observations by Ginsburg et al. (2023) demonstrated how ice absorption significantly alters star distributions on color-color diagrams in environments like the Galactic Center, as observed with JWST narrow-band photometry.These phenomena can lead to anomalous

F770WFigure 3 .
Figure 3.Comparison of the reddening law in M16 with literature results.Our derived color excess ratios (red points) are plotted alongside results from Weingartner & Draine (2001, WD01) for RV = 3.1 in blue dashed and RV = 5.5 in blue solid, and from Hensley & Draine (2020, HD20 Cyg OB2-12) in green dashed and Hensley & Draine (2023, HD23 Astrodust) in green solid, respectively.The inset is a zoom-in view of the MIR regime.

Figure 4 .
Figure 4. Color-color diagrams of F 200W −F 335M vs. F 090W −F 200W (panel a) and F 335M −F 444W vs. F 200W −F 335M (panel b), generated using color maps.The linear fits are shown with solid red lines, and the 3σ deviation ranges with dashed red lines.Outliers identified by LTS LINEFIT are marked with grey points.Reddening vectors for AV = 10 mag, calculated using the WD01 RV = 5.5 extinction law, are superimposed.

Table 1 .
StarbugII parameters used for photometry Parameter