PHANGS-JWST First Results: The Dust Filament Network of NGC 628 and its Relation to Star Formation Activity

NGC 628 is one of the initial targets observed for the PHANGS-JWST Cycle 1 Treasury project (GO 2107, PI J. Lee). Our observations of NGC 628 target the central region ($R_{gal}$ $\lesssim$ 5 kpc) of the star-forming disk, to overlap areas where HST , ALMA, and MUSE data have been obtained.²²2Footprint maps of the HST , ALMA, and MUSE observations are available at https://archive.stsci.edu/hlsp/phangs-hst. Science-ready PHANGS-HST images are also available for download. The dataset includes imaging with NIRCam (F200W, F300M, F335M and F360M) and MIRI (F770W, F1000W, F1130W and F2100W). Photospheric emission from resolved stellar populations is major component in the four NIRCam bands (excepting F335M which primarily probes the 3.3µm PAH feature), whereas MIRI traces the ISM (both PAHs and hot, small dust grains). Resolution varies from 0.066 to 0.12″(NIRCam) and 0.25 to 0.67″(MIRI). For MIRI, this corresponds to 12–32 pc at the distance of NGC 628. A detailed description of the PHANGS-JWST observations and data reduction is presented by Lee et al. (in prep.). Here, we focus on MIRI imaging, deferring analysis of filamentary structure seen in the 3.3µm PAH feature to future work in anticipation of improved astrometric alignment among the PHANGS-JWST NIRCam and MIRI imaging. Of particular importance to the analysis in Sec. 4.2 is the sky background adopted for the MIRI data. We use background-corrected images which have been tied to the sky level measured in wide-field Spitzer and WISE archival imaging, as described in the Appendix of Leroy et al. (in prep.). The MIRI background levels are currently uncertain by $\pm$0.1 MJy sr^-1.

The HST NUV-U-B-V-I (F275W, F336W, F435W, F555W, F814W) observations of the central NGC 628 disk we use were obtained by LEGUS (GO 13364, Calzetti et al. 2015) using WFC3/UVIS for NUV and U, and by R. Chandar (GO 10402) using ACS/WFC for B, V and I. All data were reprocessed by PHANGS-HST. Full details are given in (Lee et al., 2022). The ACS/WFC B-band images we use to identify dust lanes in attenuation have resolution ($\sim$ 0.09$\arcsec$, 4.3 pc), approximately 2.5$\times$ finer than MIRI F770W.

We use the nebular catalog of H ii regions derived from the integral field unit (IFU) spectroscopy of the PHANGS-MUSE survey Emsellem et al. (2022). For NGC 628, the ”convolved, optimized” resolution in PHANGS-MUSE DR 2.2 is $0.92\arcsec$, corresponding to a spatial resolution of $44\text{\,}\mathrm{pc}$. Santoro et al. (2022) and Groves et al. (subm.) used PHANGS-MUSE data to create a catalog of H ii regions and provide fluxes corrected for Milky Way and internal extinction. Only star-forming regions classified using the BPT diagram (Baldwin et al., 1981) are retained in our analysis.

PHANGS-HST resolved stellar photometry has been used to identify and characterize stellar associations as summarized in (Lee et al., 2022) and described in detail by Larson et al. (subm.). Stellar clusters in NGC 628 have been studied by Thilker et al. (2022). PHANGS-HST catalogs are publicly available³³3https://archive.stsci.edu/hlsp/phangs-cat. For both associations and clusters, fluxes for the five available HST bands were measured (using upper-limits in non-detected photometric bands) and then age, mass and reddening were estimated⁴⁴4Degeneracy between age and reddening is apparent for a subset of objects, and could be more relevant in our dust filaments. via fitting of observed SEDs (Turner et al., 2021) to solar metallicity stellar population models using cigale (Boquien et al., 2019).

We use the PHANGS-HST associations catalog based on local over-densities of NUV point-like detections at a scale of 32 pc for our analysis. The stellar associations have ages ranging up to $\sim$10² Myr. We use a subset selected to have age less than five Myr in order to limit the population to the most recent star formation activity. Clusters for our analysis were also selected with the same upper limit on age.

We identify dust filaments using FilFinder (Koch & Rosolowsky, 2015). This code applies an adaptive thresholding algorithm and graph-based medial skeleton analysis to isolate and then characterize filaments. Thresholding is conducted over local neighborhoods, allowing for the extraction of structure over a large dynamic range both in intensity and spatial scale (the latter when the code is run multiple times with different parameter choices). The potential effect of bright sources interspersed in the web of filaments is mitigated by an arctan intensity transform before thresholding. We use a slightly modified version of FilFinder in which the arctan transformed image is convolved to the filament extraction scale before each run. In the current analysis we only utilize the filament masks produced by the code, leaving FilFinder’s skeleton analysis capabilities for future work.

We apply FilFinder independently to the background-corrected image in each MIRI band and to a pre-processed version of the B-band HST image (in which attenuation features are most evident compared to NUV, U, V, and I).

For pre-processing the HST image, multi-scale median filtering is used to suppress peak-like features over a range in scale, while retaining small scale dips (concave areas). Our specific filtering method follows from Hoversten et al. (2011). At each pixel the output value is assigned to be that location’s minimum in a stack of circular median filtered images. Filter kernel diameters are taken from a ladder of physical scales (starting at the resolution limit and proceeding up to 32 pc). The result is that confusion by bright stars and stellar clusters is greatly reduced, emphasizing the dust lane structures. Fig. 1 (left) and Fig. 6 (left) show the pre-processed image. FilFinder nominally operates by finding positive filamentary features above the background. before passing the pre-processed B-band image to FilFinder, we invert the sense of the intensity (subtracting the image from a constant value equal to the maximum in the field of view)

For both attenuation and emission, we use FilFinder to identify potential filaments with narrow dimension (width) starting at 25 pc then stepping by factors of $\sqrt{2}$ (0.15 dex) up to 400 pc (25, 35, 50, 70, 100, 140, 200, 280, 400 pc). The minimum scale of 25 pc corresponds to approximately twice the PSF FWHM of MIRI F770W. For F2100W analysis, we begin at the 35 pc scale due to the larger F2100W PSF. FilFinder parameters are set as follows: size_thresh = 6$\pi$$w^{2}$, adapt_thresh = 2$w$, glob_thresh = 2$\sigma$ above sky level, smooth_size = 0.5$w$, fill_hole_size = 0.5$w^{2}$, where $w$ is the extraction scale (in pixels) and $\sigma$ is the standard deviation noise level expected in the arctan transformed, convolved images.

In addition to the filament masks generated for extraction at specific scales, we also produce masks representing the union of filaments detected cumulatively at different scales. For these cumulative masks, we sum the individual masks (for scales less than or equal to the current scale) and then flatten the result, such that it has a value of one anywhere a constituent scale contributes filament coverage. The cumulative summed mask before flattening is also retained, as it highlights the multi-scale nature of the filament network. Appendix A presents individual scale and cumulative multi-scale masks for F770W, F2100W, and HST B-band as a Figure Set.

The wide range of scales initially allowed for filament extraction is exploratory. In the second half of Sec. 4.1), we argue that all emission and attenuation filament masks for scales $>200$ pc not be used, although we do include them in plots allowing for the reader to make their own .

For Sec. 4.1, we measure filament mask overlap categorizing each pixel as belonging to one of four classes: (1) attenuation filament only, (2) emission filament only, (3) both attenuation and emission filaments, or (4) no filamentary features detected. Note that the sight line fractions we report in Sec. 4.1 for classes (1)-(3) are normalized to the total count of pixels with any detected filamentary feature, rather than the total count of pixels in the image. It is beyond the scope of the current study to quantify the dependence of pixel classes on the depth of the imaging – however, given the pervasive character of the detected filament network, we suggest that at least the F770W, F1000W, and F1130W observations are sensitive enough to make this a moot point. Our F2100W data is about 2$\times$ less sensitive in absolute terms of limiting surface brightness, $\sigma_{\mathrm{I}}$ [MJy sr^-1], and also suffers from a similar loss in resolution compared to F770W. These factors likely contribute to loss of some smaller scale filamentary structure in the F2100W imaging.

Sec. 4.2 (flux fractions) requires local estimation of the diffuse emission to obtain background-subtracted filament flux. the photutils (Bradley et al., 2022) background2D code, supplying filament masks to indicate which pixels the procedure should ignore. A mesh of bins (each having size two-thirds the filament extraction scale) is defined, and the mode of unmasked pixels is determined in each bin. These mesh modes are median filtered with a 3$\times$3 boxcar (ignoring bins with too many pixels masked as within a filament) and the resulting values are interpolated across the image grid. from the input image and non-filament pixels are set to zero. Integrating this result gives the background subtracted flux of the filament structures, which is then divided by the total flux in the MIRI footprint to obtain flux fractions.

We start by highlighting the consistency between filamentary attenuation (dust lane) features and web-like MIR dust emission to illustrate the potential of using HST -detected features as a high resolution proxy for the dense dusty ISM morphology (and perhaps even molecular gas⁵⁵5The association between the DFN structures and ¹²CO(2-1)-traced molecular gas will be investigated in a future study.).

Fig. 2 shows the results of filament overlap analysis. For the PAH-dominated bands (F770W, F1130W) and 10$\micron$ thermal emission traced by F1000W, in the top panel of Fig. 2, we find that the percentage of sight lines in common between visible attenuation and MIR emission filaments is nearly 40%, at 25 pc and declines smoothly with increasing scale (filament width). This decline is due to more area becoming traced by attenuation only for larger individual filament extraction scales. The percentage of emission-only filament sight lines declines a small amount from 25 pc to 200 pc. Overlap statistics generated on the basis of cumulative multi-scale masks show a different picture, in which the percentage of attenuation plus emission sight lines grows with scale from just over 40% at 35 pc to 55% up to 200 pc. This is a consequence of different extraction scale masks picking up varied portions of the overall web-like DFN. We return to this point further in the current subsection (see Fig. 4 and Fig. 5).

In the bottom panel of Fig. 2, we show the overlap for F2100W emission filaments with the dust lanes from HST . At 21$\micron$ the percentage is 24% for 35 pc individual extraction scale, substantially less than for the three other bands, declining to $\sim$16% for 200 pc. For all separate scales, attenuation-only sight lines amount to more than 65% and emission-only $\lesssim$10%. Given that the attenuation filament mask remains constant F770W and F2100W, this could suggest that either we are sensitivity limited for F2100W or the filamentary emission detected in 21$\micron$ imaging is less consistently recovered into coherent structures. The latter interpretation is supported by mask inspection in the figures of Appendix A, and by the fact that the morphology of the F2100W image is more dominated by compact sources (e.g. IR-bright star-forming regions, see Hassani et al. in prep.) than F770W, F1000W, and F1130W. This serves as a reminder that we are tracing warm (140 K) dust at 21$\micron$, whereas the dust attenuation provides a more complete inventory with respect to dust over a wide range of [cooler] temperatures. Nevertheless, the overlap with the 200 pc cumulative multi-scale mask is 40% (only $\sim$1/3 less than for F770W).

In summary, the 40% level of pixel-by-pixel agreement of the attenuation and emission filament masks at 25 pc and the trend for the greatest agreement on the smallest scales that the correspondence may be even tighter if smaller physical scales are probed, such as in Local Group galaxies. We stress that, although 40% overlap may sound low, inspection of Fig. 1 shows that the detailed morphology (e.g. extent, shape) of filaments that are detected in attenuation and dust emission is frequently rather well-matched. We also note that regions of the filament network only found as dust emission are expected due to line-of-sight effects: some filaments will be positioned on the ‘back’ side of the face-on galaxy disk, and could account for the majority of ‘emission-only’ areas. Some attenuation-only regions also lie in false-negatives for MIR dust emission on the scale of interest, or lie in relatively MIR-faint, intermediate surface brightness regions of the visible galaxy disk. False negatives occur when a MIRI filament structure has a width that is somewhat different than the visible attenuation feature, or when scale-matched emission in the region is biased against detection at the scale of interest (e.g. a filament centered between brighter neighboring emission features spaced by about the width of the FilFinder adaptive thresholding box).

Figure 3 shows the division of sight lines between those associated to a dust filament (here, either attenuation or emission) and those that do not lie in a filament for the particular extraction scale. The plots show very little change between bands. We find that approximately 50–55% of sight lines are attributed to dust filaments for scales 25–100 pc with no variation due to extraction scale, then a mild increase at larger scales (to $\sim$60–65% for 200 pc). Cumulative multi-scale mask measurements of the dust filament sight line fraction (dots and lines in Fig. 3) steadily rise across the 25–200 pc range, despite the constancy for individual extraction scales smaller than 140 pc. Our filament masks are detecting a morphologically diverse dusty ISM, Though beyond the scope of this paper, it will be important to investigate if the component is increasingly atomic-dominated compared to narrower filaments.

Figure 4 illustrates the tendency for dust filament masks to capture different morphological features when extracted using varied scales. In particular, we show the JWST F770W image in the top panel and HST B-band in the bottom panel. Using red and blue lines, overplotted on each of the images is a contour boundary of the 25 pc mask (dashed thin line) and 200 pc mask (solid thick line), with the variety of mask corresponding to the dust detection property of each image (emission on top, attenuation on bottom). Many of the single extraction scale 25 pc dust filaments are very long, with $l>$1 kpc and have rather high aspect ratio. However, there are also frequent cases of the 25 pc masks (dashed thin lines) only including localized substructure within larger coherent filaments left unjoined by small scale extraction. This selective property of the filament identification outcome can be seen in spiral arms (markedly less so in interarm regions) and in both emission and attenuation. Conversely, the 200 pc filament masks (solid thick lines) often entirely exclude areas with a significant population of narrow GMC-scale width filaments. We note that the 200 pc attenuation filament mask recovers continuous spiral structure more effectively than the F770W emission mask of the same scale. Additionally, Fig. 4 demonstrates that the recovered dust attenuation features can be quite modest in terms of apparent $A(B)$. Inspection in the peripheral areas of the HST panel nevertheless suggests the majority of such filaments are real. On the JWST image, we also plot a contour representing the cumulative (up to 400 pc) emission filament mask (dashed thin yellow). The 25 pc and 200 pc filament structures alone do not include all portions of the DFN, especially in dust emission. The yellow contour shows how using a cumulative mask addresses this issue, linking many smaller scale components. Using cumulative masks with large upper scales can excessively broaden the extent of filaments, and hence we urge caution in the choice of the upper limit on multi-scale integration.

Fig. 5 provides empirical justification for maximum scale on the basis of emission. The figure plots cumulative F770W mask boundaries for different maximum scales, starting at 100 pc (green), running through 140, 200, 280 pc (yellow, red, cyan), up to 400 pc (magenta). The top layer contour (100 pc, green) already structure of the DFN, but several regions of seemingly contiguous filamentary emission remain disconnected either internally or to the network. By looking at other colored contours emerging from under the green boundary, one can infer how adding progressively larger scales changes the network. The 200 pc cumulative mask (red) appears to provide reliable recovery of all filamentary emission structures without undue peripheral excess, although 140 pc and 280 pc are probably also acceptable. HST B-band filament extraction at 280 and 400 pc occasionally confuses interarm gaps as attenuation. Such large ($>200$ pc) scales also push the limit of what can be considered a filament in the sense of forming via turbulent Jeans scale gravitational instabilities (see Meidt et al. in prep.).

Fig. 6 presents a zoomed in view of a region east-southeast of the galaxy center, showing the filtered B-band data in comparison to an unprocessed version of the same image and to F1000W data from JWST . Contours of the 25 pc B-band and F1000W filaments are overplotted. This figure visually emphasizes the 40% sight line overlap and morphological agreement at this scale. Fig. 6 also demonstrates feasibility to probe even smaller scales in visible dust lanes across the entire PHANGS-HST sample. Significant substructuring (down to the $\sim 5$ pc WFC3/UVIS resolution limit in NGC 628) in attenuation of the features detected as dust emission filaments at 25 pc scales is apparent. We further expect HST to detect additional small scale filaments beyond the limit of MIRI – several super narrow attenuation filaments without corresponding larger scale intensity depression are apparent in other regions of NGC 628. Finally, Fig. 6 demonstrates that with HST we will be able to reveal candidate features we nickname ‘dust motes’ (examples marked with yellow circles), essentially compact ($\lesssim 10$ pc), dark clouds we cannot cleanly identify solely with JWST due to confusion with point-like dusty extreme AGB stars (Thilker et al., in prep.) lying outside the emission filament network (examples marked with green circles) in the short wavelength MIRI images. These dust motes could be individual molecular clouds in relative isolation. Their size is comparable to the Taurus Molecular Cloud, although Taurus is star-forming whereas the motes may often be quiescent. A better Milky Way might be the smallest scale clouds found in the 3D extinction maps of e.g. Leike et al. (2020). Dust mote clouds are challenging to confirm with MIRI, but appear in the HST B-band images before any pre-processing (Fig. 6 center). PHANGS-JWST F335M 3.3$\micron$ PAH imaging (Sandstrom et al., in prep.; Rodriguez et al., in prep.) and forthcoming PHANGS-HST H$\alpha$ imaging (GO 17126, PI R. Chandar) may prove useful to further vet dust candidates as a class.

The relative amount of structured and unstructured (diffuse) dust in a galaxy is fundamental metric of ISM morphology and offers a basis for comparison with simulations that aim to understand the impact of factors such as stellar feedback and dynamical influences (e.g. Smith et al., 2020). As noted in Sec. 4.3, this observable is also relevant to star formation rate (SFR) estimation. We measure the fraction of flux contained within the extracted filaments compared to the total flux in the image. As the filaments are brighter than their surroundings, this requires that we establish a local background estimate at all locations in the image, which we do using the procedure given in Sec. 3. Figure 7 presents our results concerning flux fraction as a function of filament extraction scale.

We find that at the smallest independent extraction scales considered (25 pc for F770W, F1000W, F1130W; 35 pc for F2100W), the background-subtracted filament flux fractions are nearly $30\%$ except for F2100W at $\approx 20\%$. The flux fraction is essentially constant versus extraction scale for F770W, F1000, and F1130W, but peaks at 35 pc and 200 pc for F2100W. The F2100W flux fractions are always significantly less than all three shorter wavelength MIRI bands. Given the present residual $\pm$ 0.1 MJy sr^-1 uncertainty in the background level of our MIRI images (see Sec. 2.1 and Leroy et al. in prep.), we assess the impact this has on filament flux fractions, by offsetting the image intensities by the uncertainty in a positive and negative sense. Dashes overplotted on the bars of Fig. 7 indicate the resulting perturbed fractions.

Cumulative multi-scale mask filament flux fraction results are also shown in Fig. 7, indicated by the solid lines and dots. As expected, the fraction increases as progressively more individual masks are combined and the maximum allowed extraction scale is increased. For our choice of maximum scale (200 pc) the background-subtracted filament flux is in the range $55$% to $60$% of the total F770W, F1000W, and F1130W flux in the field. At the same cumulative scale, the F2100W measurement is 43%. Scatter due to uncertainty in background level is also checked for cumulative values, and shown with dotted lines on Fig. 7. The impact is notably larger than for individual scales, but is generally less than systematic uncertainty due to our choice of maximum scale included in the DFN cumulative mask.

We note that if diffuse emission is neglected, rather than subtracted away as in our analysis, the integrated flux from the filament network becomes 40-100% higher boosting the cumulative filament flux fractions to $\sim$80% for F770W, F1000W, and F1130W and 60% for F2100W.

The results of Sec. 4.2 confirm that more than half of the MIR flux from NGC 628 originates in the intricately structured filament network, with the only exception being for F2100W (43%). This outcome is what we anticipated based on previous work for nearby galaxies (Liu et al., 2011; Leroy et al., 2012; Crocker et al., 2013; Calzetti, 2013; Boquien et al., 2016; Kumari et al., 2020; Belfiore et al., in prep.). It reinforces the necessity of correcting the integrated MIR luminosity of a galaxy for the presence of inherently diffuse emission heated by older stellar populations, when using MIR to measure current SFR (Lonsdale Persson & Helou, 1987; Boquien et al., 2016). One of the goals of PHANGS-JWST is to clarify systematics of this correction. We can positionally test the linkage of the DFN to age-dated markers of star formation events, checking whether all regions of the filament network can be attributed to heating by current star formation or if some filamentary structure is effectively quiescent and should be removed alongside diffuse emission when estimating SFR. Put differently, perhaps only filamentary structure up to a certain maximum scale, or down to a limiting dust surface density, is directly linked to the youngest stellar populations.

We take a first step toward the goal above by cross-correlating the filament masks with the PHANGS-MUSE H ii regions, plus PHANGS-HST young ($\leq 5$ Myr) clusters and associations, described in Sec. 2. With these populations, we measure unweighted inclusion fractions (number of objects contained within the filament mask divided by the total count in the image footprint) and also tracer-weighted inclusion fractions (summing weights corresponding to SFR rather than object counts). For weighting we use extinction corrected L(H$\alpha$), effectively unobscured SFR, for H ii region populations; and M_∗ / age $\approx$ effective SFR(5 Myr) for young clusters and associations.

Fig. 8 plots the measured inclusion fractions for each of the test populations. We begin by making general comments about Fig. 8 which are applicable to all panels. We find that the inclusion fractions (at individual scales or cumulatively) are consistently maximized for H ii regions but decline slightly when using young clusters and dramatically using young associations, regardless of scale or band. The inclusion fraction measured at individual filament extraction scales for H ii regions and young clusters ($\leq 5$ Myr), decline from a maximum at 25 pc until reaching 70–100 pc. For associations, there is far weaker, if any, dependence on scale over the same range, but an apparent enhancement at 140–200 pc scales. The random inclusion fractions for unweighted measurements (top panels) provide further insight. They indicate strong correlation between dust emission filaments and H ii regions / young clusters, moderate anti-correlation of associations with filaments at scales $<$100 pc, and increasingly strong correlation of associations with the DFN for extraction scales of 100, 140, and 200 pc. For measurements made with cumulative masks, the inclusion fractions (unweighted and weighted) increase with scale as would be expected. The slope of the cumulative curves in Fig. 8 (top) is consistent with the scale dependence of the randomized control, showing the influence of significant multi-scale mask covering fraction, but normalization is significantly higher (excess over random $20\pm 5\,\sigma$ for H ii regions, 3–6$\sigma$ for young clusters) except for associations. For associations, the small scale signal of anti-correlation dominates until our cumulative measurement finally reaches random equivalence by 200 pc.

The SFR-weighted inclusion fractions Fig. 8 (bottom) are perhaps more physically relevant to the questions raised at the start of this Section. Remarkably, our measurements show about 75% of the SF traced by H ii regions is occurring within the 25 pc scale of the F770W filament network (bottom left), whereas the peak SFR-weighted inclusion fraction for F2100W is $\sim$45% (at 35 pc, bottom right). We find the single scale SFR-weighted inclusion fractions are highest overall for F1130W and F1000W at nearly 80% for 25 pc filament extraction (not shown). Perhaps the most striking measurement linking the dust emission filaments to current star formation is the attainment of $\sim 95$% SFR-weighted inclusion by our uppermost cumulative scale, 200 pc, for each of F770W, F1000W and F1130W when using the multi-scale masks, and 72% for F2100W.

We anticipated seeing a reduction in inclusion fraction for stellar associations, since, although the associations are young, they are star formation events that have evolved from clusters via disruption/dissolution or were initially formed unbound, and either way seem more likely to have already cleared their environment of dust. However, the anti-correlation observed at small scales is an unexpected demonstration that such feedback actively helps to sculpt the dust into bubble/shell structures (e.g. Watkins et al., in prep.; Barnes et al., subm.).

Work is currently underway to obtain a catalog of PHANGS-ALMA GMCs associated with embedded ($t$ $\sim$ 0 Myr) star formation, indicated by compact F2100W sources coincident with a GMC (Lessing et al., in prep.). When ready, we will compare our dust emission filaments to that population. We emphasize that the F2100W source population does not drive the identification of filaments at $21\micron$ due to FilFinder’s transformation of the image intensity scale.