PHANGS–JWST First Results: The 21 μm Compact Source Population

Dust grains and polycyclic aromatic hydrocarbons (PAHs) in the interstellar medium (ISM) play a central role in shaping the radiation field in galaxies, converting short-wavelength light from stars and other emitters into long-wavelength emission in the infrared (IR; Galliano et al. 2018). The mid-IR emission observed by the James Webb Space Telescope (JWST) is generated by small dust grains, which are heated stochastically to temperatures of ≳100–150 K to emit photons at mid-IR (5 ≲ λ/μm < 60) wavelengths (Draine & Li 2001). The PAHs reemit absorbed radiation in discrete spectral band features in the range 3 < λ/μm ≲ 21, which corresponds to the stretching and bending modes of bonds in the large molecules (Allamandola et al. 1989; Li 2020). Because the IR radiation is significantly less affected by absorption, dust and PAH emission provide a vital low-extinction view into the densest regions of galaxies.

Thanks to previous generations of mid-IR observatories, the properties of warm dust and PAHs have been broadly surveyed. Such studies have illustrated the crucial role of IR observations in understanding the process of star formation in galaxies, since star formation occurs in high-extinction regions (e.g., Perez-Gonzalez et al. 2005; Jarrett et al. 2013). In particular, embedded high-mass stars create strong radiation fields that enhance the mid-IR emission so that the mid-IR can be used as a tracer for star formation (Calzetti et al. 2007; Kennicutt & Evans 2012). Mid-IR data are usually paired with a short-wavelength tracer of unobscured star formation to approximate the full star formation rate (Kennicutt et al. 2009; Hao et al. 2011; Leroy et al. 2019). However, mid-IR emission is not a linear tracer of the star formation rate at high resolution. In addition to the decorrelation of different stages of the star formation process at small (<1 kpc) scales (e.g., Onodera et al. 2010; Schruba et al. 2010; Kruijssen & Longmore 2014), the mid-IR emission also arises from the diffuse, average interstellar radiation field (ISRF) heating the neutral ISM (Boquien et al. 2015; Leroy et al. 2023; Sandstrom et al. 2023a).

Understanding the nature of the correlation between star formation and the mid-IR has been challenging because the resolution of the previous generation of observatories had comparatively poor resolution: 65 for the Spitzer/MIPS 24 μm band (Werner et al. 2004) and 119 for the WISE 22 μm band (Wright et al. 2010). The study of individual star-forming regions at <100 pc scales was thus limited to the Milky Way, the Local Group, and the nearest (d ≲ 3 Mpc) galaxies (e.g., Peeters et al. 2002; Meixner et al. 2006; Verley et al. 2007; Carey et al. 2009; Chastenet et al. 2019). Whereas PAH emission is distributed throughout the neutral ISM, several authors have noted the spatial correspondence between Hα emission and the mid-IR continuum (Rice et al. 1990; Helou et al. 2004; Verley et al. 2007; Relano & Kennicutt 2009). These moderately resolved extragalactic studies suggested that the PAH emission tends to be found at the edges of H ii regions; this scenario is supported by the high physical resolution observations of individual Milky Way objects (e.g., Carey et al. 2009). Given the close link between star formation and IR emission, extragalactic observations have used the mid-IR to identify compact regions of star formation. Such studies characterize the population of stellar clusters and compact associations that are in the process of forming and estimate the properties of young stellar structures (Sharma et al. 2011) or the characteristic spacing between regions (Elmegreen & Elmegreen 2019).

With the launch of JWST, we now have the opportunity to make sensitive, higher-resolution observations at mid-IR wavelengths. In particular, the JWST/MIRI observations at 5–28 μm have an order-of-magnitude improvement in resolution over Spitzer/MIPS and superior sensitivity. This capability offers a new opportunity to push studies of resolved mid-IR emission to galaxies well beyond the Local Group.

This paper explores the relationship between compact sources seen in MIRI 21 μm observations and other tracers of star formation. We use JWST observations of the first four galaxies to be observed for the Physics at High Angular resolution in Nearby GalaxieS (PHANGS) Treasury program (02107; PI: J. Lee) in the near- and mid-IR (Lee et al. 2023). We combine these new observations with the rich set of supporting data about the star formation process gathered as part of the PHANGS surveys. We aim to understand the nature of the compact 21 μm sources seen in the JWST imaging data, extending the studies from the Local Group to a more diverse set of star-forming environments and determining if these sources include a set of truly embedded star-forming regions. But even with the 067 FWHM of the 21 μm filter with JWST, these maps are not resolving individual stars or even stellar clusters (e.g., Rodriguez et al. 2023). We thus face the challenge of how to identify and extract a uniform set of regions.

To achieve these goals, we first develop a method to find compact sources in images with large amounts of diffuse emission (Section 3). Then, we use the colors (flux ratios) of the sources in JWST bands to reject objects that are unlikely to be associated with the ISM in the galaxies (Section 4). Finally, we characterize the populations of the sources in these different systems in Section 5.

We analyze MIRI data from four galaxies in the PHANGS–JWST survey described by Lee et al. (2023): IC 5332, NGC 0628, NGC 1365, and NGC 7496, with properties summarized in Table 1. The data presented in this paper were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute. ²⁸ We focus our analysis on the sources detected in the MIRI 21 μm data (F2100W filter; FWHM: 067), though we consider all data from the PHANGS–JWST filter set, which includes NIRCam (F200W, F300M, F335M, and F360M) for NGC 0628, NGC 1365, and NGC 7496 and MIRI (F770W, F1000W, F1130W, and F2100W) for all four targets. ²⁹ Lee et al. (2023) described the data reduction process, including modifications of the default JWST pipeline and postprocessing for the final images. Of note, the modifications correct the processing for strong 1/f noise in the NIRCam data, correct for off-galaxy MIRI background images where available, and match background levels between individual MIRI fields while anchoring the overall background levels to low-resolution archival data. This survey is being carried out in the context of the broader PHANGS survey, and we include maps of molecular gas content from the Atacama Large Millimeter/submillimeter Array (ALMA) ¹²CO(2−1) integrated intensity maps presented in Leroy et al. (2021), the Hα and Hβ maps extracted from observations with ESO's Very Large Telescope (VLT) using the Multi Unit Spectroscopic Explorer (MUSE; Emsellem et al. 2022), and broadband optical imaging data from the Hubble Space Telescope (HST) as presented in Lee et al. (2022). Figure 1 shows the MIRI images of our four targets.

Zoom In Zoom Out Reset image size

Table 1. Summary of Galaxy Properties and Observational Parameters

Galaxy	Da	${\mathrm{log}}_{10}({M}_{\star }/{M}_{\odot })$ b	Re b	ℓF2100W c	ℓALMA c	ℓMUSE c
	(Mpc)		(kpc)	(pc)	(pc)	(pc)
IC 5332	9.01	9.7	3.6	29	32	38
NGC 0628	9.84	10.3	3.9	32	53	44
NGC 1365 d	19.57	10.0	2.8	64	130	110
NGC 7496 d	18.72	11.0	3.8	61	150	77

Notes.

^a Distances to targets based on observations aggregated in Shaya et al. (2017) and Anand et al. (2021). ^b Galaxy properties including stellar mass (M_⋆) and effective radius (R_e ) based on the PHANGS compilation presented in Leroy et al. (2021). ^c Projected linear resolutions at distances of the different galaxies based on observations with JWST (067) and the varying resolutions of data from PHANGS–ALMA (07–17; Leroy et al. 2021) and PHANGS–MUSE (08–12; Emsellem et al. 2022). ^d Objects NGC 1365 and NGC 7496 host IR-bright AGN that saturate several of the JWST images.

Download table as:  ASCIITypeset image

To extract compact sources from the MIRI 21 μm data, we use the astrodendro software package (Robitaille et al. 2019). We adopt this approach instead of a point-source (star) identification approach because of the extended diffuse structure present throughout the maps (Figure 1) and because the ISM-tracing nature of the 21 μm band can lead to irregular morphologies even for compact sources. The highly structured shape of the F2100W point-spread function (PSF) adds some additional complexity that confounds many point-source identification algorithms. Dendrograms provide a nonparametric description of the contour structure in the emission (Rosolowsky et al. 2008), which we then filter to identify our sources of interest.

We first estimate the 1σ noise level in the map by inferring the difference in brightness between data at the 2.5th percentile in brightness compared to the 16th percentile in brightness, which would correspond to the interval between the −2σ and −1σ intervals in a Gaussian distribution, leading to noise levels of ≈0.3 MJy sr⁻¹. This noise interval is likely an overestimate, since even these percentiles of brightness distribution may contain real flux. We resort to these empirical measures because the pipeline noise maps in the First Results data include the effects of the MIRI coronagraph pixels, which are not actually included in the real maps. We assume that a single value of noise characterizes the whole map, though there is some spatial inhomogeneity in the noise structure at the 30% level.

We mask the data and consider only regions above 5σ in brightness. Then, considering the region within the emission mask, we generate a dendrogram representation of the data, which finds local maxima in the map and identifies the contour levels below which each pair of local maxima is connected. We select local maxima that are >2σ above the contour level at which that maximum merges with another maximum. We further require that the maxima be spatially separated by at least one half-width of the PSF at 21 μm (033). This approach yields a set of local maxima that are well defined and significant with respect to the local background. This approach also identifies the diffraction structure from the PSF around bright sources as separate sources, so we further filter this set of local maxima by computing a roundness statistic around each peak inspired by daophot (Stetson 1987). Our statistic computes the absolute difference between a 2'' square subimage and a version of the subimage that has been rotated by 180° around the local maximum. We compare the sum of the absolute difference to the sum of the image, rejecting all local maxima for which this ratio is greater than 2, with this threshold value being chosen to eliminate diffraction features.

For each source, we measure the broadband optical-IR spectral energy distribution (SED) using the JWST and HST data. We convert the HST maps to surface brightness units (MJy sr⁻¹) and use circularized convolution kernels generated using the methods in Aniano et al. (2011) to match the HST data to the F2100W resolution (067). We sample each map at the location of the 21 μm local maximum, and we subtract an estimate of the median background calculated in an annulus with radii between 2× and 3× the width of the PSF. Note that these background levels are somewhat sensitive to the size of the annulus selected because of the strongly varying backgrounds and the possibility that the 21 μm peaks are extended. We estimate the local surface brightness uncertainty as the standard deviation of the data in this annulus, calculated using a median-absolute-deviation–based estimator. In the F2100W band, this leads to uncertainties from 4.1 to 16 μJy (16th to 84th percentiles). We will revisit this approach in generating a more robust catalog in future work spanning the full PHANGS–JWST sample.

Contrasting the peak intensity with this local background estimate, we also reject any source where the surface brightness at F2100W has a peak-to-background ratio of <5. The net effect of these filters is to find pointlike and slightly resolved bright sources in our images with typical diameters of ≲2'' × 067.

We convert surface brightness to flux density by measuring the solid angle subtended by the MIRI PSF in the F2100W band (Ω_F2100W = 1.61 × 10⁻¹¹ sr), calculated using the models from the WebbPSF package (Perrin et al. 2014). We note that we complete our analysis at a common angular resolution, but distances to our targets vary by a factor of 2. Hence, our analysis implicitly includes changing physical scales (Table 1) and different minimum recoverable luminosity between the different objects, and we note where this effect shapes our results below.

The active galactic nuclei (AGN) at the centers of NGC 1365 and NGC 7496 saturate the detector and scatter light over a wide area in the map. To identify structures that could be caused by this scattered light, we match the F2100W PSF models to the location of the AGN and the orientation of the diffraction pattern, rejecting sources where the PSF amplitude is >0.1% of the PSF maximum. This ends up covering a region of ∼15'' in diameter (∼1.3 kpc in these targets). Some other bright sources in the center of NGC 1365 are also saturated and not included in this analysis, though these are further explored in Schinnerer et al. (2023) and Liu et al. (2023).

Our final catalog consists of 1271 compact sources, with 188 in IC 5332, 502 in NGC 0628, 367 in NGC 1365, and 214 in NGC 7496. In addition to the JWST- and HST-based SEDs, we also sample the MUSE Hα and Hβ maps and the ALMA CO(2−1) integrated intensity and line width maps at the location of the 21 μm local maxima. For the ALMA and MUSE measurements, we do not perform any background subtraction. These maps are at coarser resolution (up to 17 depending on band and target; see Table 1) and thus yield slightly smoothed estimates of the Hα, Hβ, or CO(2−1) that would be measured at the JWST resolution. The subsequent CO-based analysis does not strongly depend on the differences in resolution, but measurements of Hα luminosity may be overestimated and could affect our results. In the case of a flat Hα map over the MUSE PSF area, the Hα brightness could be overestimated relative to the JWST beam by the ratio of PSF areas, which could be up to a factor of 3 in this extreme case (NGC 1365). However, the galaxies show similar resolutions between MUSE and JWST, so even for the extreme case, the overestimate will be less severe.

Of the 1271 compact sources, 1083 have NIRCam as well as MIRI data (recall that IC 5332 currently lacks NIRCam imaging). Figure 2 illustrates the SEDs for our sources. Most sources are bright at 21 μm by selection, with a characteristic "dip" at 10 μm defined by the gap between the bright PAH features at 7.7 and 11.3 μm. Sources with NIRCam imaging also show the PAH feature at 3.3 μm as a "bump" in the SED. There is modest variation in the median SEDs between sources in different galaxies, and the bright sources in the nuclear star-forming ring of NGC 1365 increase the median flux profile for that galaxy.

Zoom In Zoom Out Reset image size

To better understand the nature of the identified sources, two authors (H.H. and E.R.) visually examined the SEDs and an atlas of multiwavelength images showing the full set of native-resolution JWST, HST, ALMA, and MUSE data for each source. We find that most sources appear to be star-forming regions associated with slightly resolved clusters or associations of reddened stars, bright PAH features, and CO or Hα emission. We refer to these sources as "ISM" sources because their SED reflects the presence of strong PAH features, which emerge from dust mixed with ISM material. These ISM sources are typically embedded in diffuse filamentary structure visible in the F2100W, F1130W, and F770W bands (Thilker et al. 2023).

In Figure 2, we also show the SEDs of the visually identified non-ISM sources. Several of these appear to be dusty stars, likely extreme-AGB stars (Corbelli et al. 2011; D. Thilker et al. 2022, in preparation), which are recognizable by a smooth SED in the near- and mid-IR with weak PAH features and a pointlike morphology in F200W. We designate these as "STAR" sources. The remaining sources appear to be background galaxies. These have an extended shape in F200W images and no PAH features, likely because the PAH features are redshifted out of the corresponding JWST bands. We refer to these as high-z or HafeZ ³⁰ sources (hereafter HZ).

Leveraging these visual classifications, we attempt to identify color cuts that could effectively achieve the same assignments using ratios among the MIRI and NIRCam bands. Given the prominence of PAH features in the ISM sources and their absence or weakness in the HZ and STAR sources, we focus on the brightness of the bands with strong PAH features relative to the nearby PAH-free bands. The resulting color comparisons are illustrated in Figure 3. In the left panel, we plot r₇ ≡ f_F770W/f_F1000W versus r₁₁ ≡ f_F1130W/f_F1000W, which normalizes the bright PAH features in the F770W and F1130W bands by the nearby F1000W band. Here F1000W is expected to reflect mostly continuum or silicate absorption. The ISM-like features with strong PAH emission fall along a tight locus with r₁₁ ≈ 1.25r₇. These ISM features clearly separate from the STAR/HZ sources in this space, and we define a boundary

$$\begin{eqnarray}&&{r}_{11}\gt 2.5-{r}_{7},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{r}_{11}\lt {r}_{7}+1.5,\end{eqnarray} \tag{ 2 }$$

such that nearly all sources in this region are classified as lines of sight dominated by the ISM emission. Sources outside this region lack a dip in the 10 μm relative to the 7.7 and 11.3 μm PAH features that is characteristic of typical ISM emission from local star-forming galaxies.

Zoom In Zoom Out Reset image size

We construct a similar color–color diagram using the NIRCam bands (Figure 3, right). Here we focus on identifying lines of sight with bright PAH emission by contrasting the F335M band, which covers the 3.35 μm PAH feature, against the more continuum-dominated nearby F330M and F360M. For r_3.0 ≡ f_F335M/f_F300M and r_3.6 ≡ f_F335M/f_F360M, the regions with a significant PAH feature show

$$\begin{eqnarray}&&{r}_{3.6}\gt 1.7\mbox{--}0.5{r}_{3.0},\end{eqnarray} \tag{ 3 }$$

which appears as a dashed line in the right panel of Figure 3. For reference, we also show the line defined by Lai et al. (2020) that indicates the locus of stellar continuum in the JWST bands and the locus of PAH emission defined by Sandstrom et al. (2023b; see also Rodriguez et al. 2023).

Our partition of this color space generally yields an accurate match to our manual source classification. The MIRI-based diagnostic (Figure 3, left) is 94% accurate on the 1271 sources with measured MIRI colors. Here accuracy is defined as the fraction of correct classifications (true negatives and true positives) compared to the by-eye SED and image classification. There still remain notable false positives and negatives, which stem partly from complex backgrounds affecting the extracted flux values. The NIRCam-based partition is only 88% accurate on the 1083 sources with NIRCam data (recall again that IC 5332 lacks NIRCam data). This lower accuracy results in substantially more false-negative classifications using NIRCam colors (128) than using MIRI colors (34). These false negatives are sources identified as not being ISM features because of the lack of a strong 3.35 μm PAH feature but that appear to be ISM sources in our manual assessment. If we classify a source as being an ISM feature based on having at least one diagnostic consistent with the color cuts, the accuracy of the classification rises marginally to 95% for the 1083 sources with both MIRI and NIRCam measurements. Misclassification could also occur when there is significant PAH (or ISM) emission in the bands we are using as "continuum" bands: F1000W, F300W, and F360W. For example, Lai et al. (2020) noted that the 3.47 μm PAH band may contribute to the F360W light, which is also seen in the PHANGS–JWST sample (Sandstrom et al. 2023b). Leroy et al. (2023) showed evidence that F1000W may have a significant contribution from the wings of the PAH features in adjacent bands. The 9.7 μm silicate feature is often found in absorption in evolved stars (e.g., Werner et al. 1980), which could make these stars appear more ISM-like.

To maximize the accuracy of our assessment on the largest number of sources while retaining the ability to apply these results to other targets, we proceed using the MIRI-based diagnostic as the primary source classifier. Future work will rely on such color classification, since the number of sources will be too large for human classification with continued JWST operations.

While our compact source catalog contains 187 non-ISM sources from background galaxies and stars, these sources do not contribute significantly to the 21 μm flux density of our targets. Using the MIRI classifier, we find that 97% of the total F2100W flux density in our catalogs comes from sources classified as "ISM" across all targets.

Young massive clusters and associations are thought to form in dark clumps of gas that the newly formed stars then disrupt with winds and radiation. There should thus be a phase of star formation where the young stars are embedded and visible primarily in the IR before they break out and become visible in the optical. This picture has been borne out in other systems (e.g., Kim et al. 2021). Our data broadly support this picture; from the visual inspection, the compact 21 μm bright sources detected by JWST in these targets appear to be primarily star-forming regions hosting high-mass star formation, though relatively few sources are associated with a truly optically dark phase (see also Kim et al. 2023; Whitmore et al. 2023). In this section, we quantify this impression by examining the associations between the 21 μm sources and other star formation tracers and then examining embedded sources. For this analysis, we restrict our analysis to 1085 compact sources that are classified as "ISM" by the MIRI-based flux ratio criterion.

5.1. SED Fitting and Derived Properties

To characterize these compact sources, we use the cigale code (Boquien et al. 2019) to estimate the age and mass of stellar populations that could yield our measured JWST and HST luminosity densities, as well as some properties of the surrounding photodissociation region. For the fit, we consider a double exponential star formation history with an e-folding time ranging from 0.01 to 10 Myr and a late starburst with an e-folding time ranging from 0.001 to 0.1 Myr. Based on the characteristic timescales for star formation in our galaxies (e.g., Chevance et al. 2022), we consider a range of stellar population ages from 3 to 25 Myr. Based on the the Balmer decrement analysis of Hα and Hβ emission (Emsellem et al. 2022; Groves et al. 2022; Santoro et al. 2022; Belfiore et al. 2022), we restrict our models to 0.1 < E(B − V) < 3, and we use a Calzetti et al. (2000) attenuation law. We include a Draine et al. (2014) dust model with a mass fraction of PAHs ranging from 0.47% to 7.32%. Since individual H ii regions can have intense radiation fields, we consider $0.5\lt {U}_{\min }/{U}_{0}\lt 50$, where U_min is the parameter from Draine & Li (2007) that characterizes the ISRF, and U₀ is the ISRF in the solar neighborhood (Mathis et al. 1983). We consider a fraction of the dust mass heated by the power-law distribution of radiation fields with intensity up to ${U}_{\max }={10}^{7}\,{U}_{0}$ with a fraction of γ = 0.1. Furthermore, we allow the exponent of power-law distribution to span a range of 1–3 in intervals of 0.5. Finally, we choose a Chabrier initial mass function (IMF) for our model. We summarize our results in Figures 4 and 5. Our fits typically find that the 21 μm sources have SEDs consistent with stellar population masses of 10² < M_⋆/M_⊙ < 10^4.5. We find a range of mass-weighted ages spanning 2–25 Myr, though most objects are <8 Myr. Based on the cigale fits, we find a median ratio of $\nu {L}_{\nu }/{M}_{\star }=98\,{L}_{\odot }{M}_{\odot }^{-1}$ for our sources with a variation of ±0.3 dex.

Zoom In Zoom Out Reset image size

Figure 4 shows the luminosity distribution of sources in the four galaxies, ordered by distance to the target. All distributions show a turnover at the low-mass end that is due to completeness, since the value corresponds well to the typical 5σ luminosity uncertainties (blue dotted lines in Figure 4). The noise uncertainty is estimated from the typical errors in the flux measurements and thus includes both image noise and the complexity of the background. Above this limit, the distributions show good agreement with the power-law distribution function (red lines) shown in the figure with the form of ${dN}/{{dL}}_{\nu }\propto {L}_{\nu }^{\beta }$ running over 2 orders of magnitude in source luminosity. We use the median ν L_ν /M_⋆ ratio from the cigale fits to define the top axis of the graph. We fit each binned distribution with a simple linear regression ranging from the 5σ limit to 100× this value to provide an estimate of the index of the distribution. We assess uncertainties by rebinning the data to different widths and offsets and refitting the distribution. In NGC 1365, there is a tail of sources at high luminosity, which consists of the sources in the nuclear star-forming ring of this system that are not included in the fit. These high-mass young sources are explored further in Whitmore et al. (2023).

Our derived values of the power-law index (typically β ≈ −1.7 ± 0.1) are slightly shallower (top-heavy) compared to standard analyses of star cluster masses, which find β = −2.0 ± 0.3 (Mok et al. 2020), though analyses of young clusters sometimes find shallower indices (Whitmore et al. 1999; Adamo et al. 2017, 2020) for the cluster mass function. The slopes are essentially the same as the H ii region luminosity function in PHANGS targets (β_{H II} = −1.7; Santoro et al. 2022). Though careful source identification and detailed SED fitting remain the most precise way to assess the cluster mass function, the overall agreement between different results and reasonably stable mass-to-light ratios suggest that, similar to H ii region luminosities, the IR peaks may provide a useful indicator of the mass distribution of young associations and clusters. If we adopt a constant scaling of ν L_ν /M_⋆, the slope of the luminosity distribution would also trace the mass distribution of these objects. These values can be compared to the mass distribution of clusters from optical analyses, though we are not directly selecting for long-lived bound objects.

5.2. Association with Star Formation Tracers

5.3. A Dearth of Embedded Sources

There is an excellent correlation between the 21 μm luminosity of a source and the luminosity of the Hα line at the source location (see also Leroy et al. 2023). In Figure 5, we plot the correlation between the luminosity of the sources in the F2100W and the extinction-corrected Hα luminosity calculated using the Balmer decrement method applied to the MUSE Hα and Hβ maps (Belfiore et al. 2022). The Hα luminosity is calculated by sampling the Hα flux maps and scaling by the projected area of the JWST PSF: L_Hα = f_Hα Ω_F2100W d². The figure shows a nearly linear correlation between these two tracers throughout our sample with 〈L_Hα /ν L_ν,F2100W〉 = 0.05. A linear relation describes the population above L_ν,F2100W ≈ 10¹⁹ W Hz⁻¹, but below this threshold, there is a population of sources that have significant 21 μm emission but minimal or no associated Hα emission. We also investigated these results on a per-galaxy basis and found no significant differences in these behaviors despite the variations in luminosities sampled (Figure 4). These sources are our candidates for embedded sites of star formation based on the 21 μm data.

Zoom In Zoom Out Reset image size

We construct a boundary at L_Hα /(ν L_ν,F2100W) = 1/80 (dashed black line in Figure 5) and examine the sources with L_ν,F2100W < 10¹⁹ W Hz⁻¹. We find that 105 of the 1085 (10%) ISM-classified sources are "embedded" under these criteria. These sources are typically found in extinction features in the HST data (Thilker et al. 2023). We find that the median aperture-subtracted flux densities in the HST bands are all typically negative at the 1σ–2σ level, whereas sources with Hα have positive optical fluxes, on average. Visual examination shows that several sources overlap with bright CO clouds in the ALMA data (∼1'' resolution); however, a comparable number of sources are found in small voids in the CO emission. This diversity of CO morphologies does not indicate a clear link to a specific phase of molecular cloud evolution. We also examined whether these embedded sources are distinct in several JWST color spaces but found no significant differences between embedded and optically visible sources. This result is similar to the results of Kim et al. (2021, 2023), who modeled the life cycle of a typical star-forming region in NGC 0628 and other nearby galaxies based on the location and scale dependence of the CO–to–21 μm and Hα–to–21 μm ratios. They found that the period of time when a region emits bright 21 μm emission only slightly precedes the time when it emits bright Hα (specifically, they inferred regions to be Hα-bright for 70% of the total time they emitted bright 21 μm), which agrees well with our observation. Note, however, that the current analysis does not apply a background subtraction to the extinction-corrected Hα maps, in contrast with the analysis of the H ii regions in Section 5.2 and the timescale analysis in Kim et al. (2023). Hence, the presence of diffuse emission may affect our results.

Based on our cigale fits, we see only a few high-mass clusters (>10⁴ M_⊙) in NGC 1365 that are deeply embedded (below the black dashed line). These objects are explored more in Whitmore et al. (2023), Schinnerer et al. (2023), and Liu et al. (2023). For low mid-IR luminosities, all of the young (<10 Myr) embedded sources have SED fits that also find low-mass measurements (M < 10³ M_⊙), reaching levels where stochastic sampling of the IMF can become important (da Silva et al. 2012) and a cluster may not form stars with significant ionizing flux. However, in smoothing the data to the relatively coarse resolution of the F2100W filter, we may obscure truly embedded sources by combining the light with brighter surrounding sources. Examining these sources in shorter-wavelength and thus higher-resolution data is another promising route for identifying embedded sources (Rodriguez et al. 2023).

In summary, we find only a few (<10%) sources visible at 21 μm that are candidate embedded star-forming regions traceable primarily in the mid-IR. These sources are fainter than the overall population of 21 μm emitters associated with young, low-mass stellar structures (M < 10³ M_⊙, age < 10 Myr) as a whole and contribute <1% of the 21 μm luminosity.

5.4. Spatial Associations with Gas versus H ii Regions

We have identified and investigated the compact 21 μm source population in four nearby galaxies observed with JWST as part of the PHANGS project. Using a dendrogram-based source identification algorithm, we identified 1271 sources in the survey area and measured their flux densities in the JWST bands (Section 3). By using the characteristic signature of PAH emission to make color cuts on the MIRI and NIRCam photometric data, we are able to identify and exclude a contaminating population of dusty stars and background high-z galaxies. This results in 1085 sources that show emission features consistent with being ISM in the target galaxies. For these ISM-like sources, we find the following.

• 1.  
  Multiband SED fitting of these compact objects shows that they are associated with stellar structures (clusters, associations) with stellar population mass 10² < M_⋆/M_⊙ < 10⁴. We find ν L_ν,F2100W/M_⋆ ≈ 98 L_⊙/M_⊙ with a range of ±0.3 dex.
• 2.  
  In each galaxy, the luminosity distribution is well described by a power-law relationship ranging over 2 orders of magnitude with the form ${dN}/{{dL}}_{\nu }\propto {L}^{\beta }$. The power law ranges from β = −1.8 to −1.6 (Figure 4), which agrees well with the luminosity distribution of H ii regions and the mass distribution of young clusters in previous work.
• 3.  
  The 21 μm sources are nearly always spatially coincident with an H ii region (>92% of sources). The spatial correlation with stellar associations identified in HST is also strong (74%), and both of these correlations are significantly stronger than random association.
• 4.  
  The luminosity of the 21 μm sources correlates linearly with the attenuation-corrected Hα emission over ∼3 orders of magnitude in luminosity (Figure 5). However, there is a tail of sources with low 21 μm luminosity that are underluminous in Hα emission. These may correspond to embedded sources. However, potential embedded sources are comparatively rare and not uniformly linked to bright CO emission, as would be expected if they were young massive clusters.
• 5.  
  While the 21 μm sources are strongly spatially associated with H ii regions, they are also associated with molecular clouds, as would be expected for these regions forming inside GMCs in a standard model of star formation.

As a "first-order" approximation, the 21 μm compact sources seen by JWST are H ii regions that are already optically visible. However, only 14% of the H ii regions have an associated 21 μm source. The 21 μm sources are preferentially associated with the bright H ii regions, and the low fraction of association seems to arise because the PHANGS–MUSE catalog of H ii regions is more sensitive than these initial 21 μm catalogs. This work demonstrates the capability to use broadband mid-IR imaging to efficiently identify the full set of compact star-forming regions in a galaxy.

This research was conducted as part of the PHANGS collaboration.

This work is based on observations made with the NASA/ESA/CSA JWST and Hubble Space Telescopes. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST and NASA contract NAS 5-26555 for HST. The JWST observations are associated with program 2107, and those from HST are associated with program 15454.

This paper makes use of the following ALMA data:

ADS/JAO.ALMA#2017.1.00886.L,

ADS/JAO.ALMA#2018.1.01651.S.

ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

Based on observations collected at the European Southern Observatory under ESO programs 094.C-0623 (PI: Kreckel), 095.C-0473, 098.C-0484 (PI: Blanc), and 1100.B-0651 (PHANGS–MUSE; PI: Schinnerer), as well as 094.B-0321 (MAGNUM; PI: Marconi), 099.B-0242, 0100.B-0116, 098.B-0551 (MAD; PI: Carollo), and 097.B-0640 (TIMER; PI: Gadotti).

H.H. and E.R. acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), funding reference No. RGPIN-2022-03499. M.C. gratefully acknowledges funding from the DFG through an Emmy Noether Research Group (grant No. CH2137/1-1). COOL Research DAO is a Decentralized Autonomous Organization supporting research in astrophysics aimed at uncovering our cosmic origins. J.M.D.K. gratefully acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program via the ERC Starting Grant MUSTANG (grant agreement No. 714907). T.G.W. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 694343). E.J.W., R.S.K., and S.C.O.G. acknowledge funding provided by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project-ID 138713538—SFB 881 ("The Milky Way System," subprojects A1, B1, B2, B8, and P1). M.B. acknowledges support from FONDECYT regular grant 1211000 and by the ANID BASAL project FB210003. J.K. gratefully acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the DFG Sachbeihilfe (grant No. KR4801/2-1). K.K. gratefully acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in the form of an Emmy Noether Research Group (grant No. KR4598/2-1; PI: Kreckel). F.B. would like to acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 726384/Empire). R.S.K. acknowledges support from the European Research Council via the ERC Synergy Grant "ECOGAL" (project ID 855130), the Heidelberg Cluster of Excellence (EXC 2181-390900948) "STRUCTURES," funded by the German Excellence Strategy, and the German Ministry for Economic Affairs and Climate Action for funding in project "MAINN" (funding ID 50OO2206). A.K.L. gratefully acknowledges support by grants 1653300 and 2205628 from the National Science Foundation, award JWST-GO-02107.009-A, and a Humboldt Research Award from the Alexander von Humboldt Foundation.

Facilities: JWST(NIRCam - , MIRI) - , HST(WFC3 - , ACS) - , ALMA - , VLT(MUSE). -

Software: astrodendro (Robitaille et al. 2019), astropy (Astropy Collaboration et al. 2013, 2018), CIGALE (Boquien et al. 2019), photutils (Bradley et al. 2020).