The PHANGS–JWST Treasury Survey: Star Formation, Feedback, and Dust Physics at High Angular Resolution in Nearby GalaxieS

Our NIRCam primary observations are divided into four separate sequences to maximize observing efficiency (Table 5). This allows for MIRI sky background observations to be taken in parallel with each NIRCam sequence, with one sequence for each of the four MIRI filters in our program. With NIRCam, we use Module B (FULL subarray) and an INTRAMODULEBOX-4 primary dither pattern with no additional subpixel dither. ⁷¹ We obtain four F200W exposures, one at each INTRAMODULEBOX-4 dither position (for a total of 16 exposures), using the BRIGHT1 readout pattern with one integration per exposure. Of these independent sequences, two have five groups/integration, contributing 386.5 s each to the total F200W exposure time, and the other two have three groups/integration, yielding 214.7 s each. The total F200W exposure time is thus 1202.5 s pointing^–1.

Simultaneously with F200W primary imaging in the short-wavelength channel (and coordinated parallel MIRI imaging of sky background), we observe the target with F300M, F335M, and F360M in the long-wavelength channel. The F300M and F335M observations are paired with the two longer F200W exposure sequences. The exposure time for F300M and F335M is consequently 386.5 s (each). During both of the two shorter F200W exposure sequences, F360M imaging is obtained, for a total exposure time of 429.5 s.

For our MIRI primary observations, we use a four-point dither pattern optimized for extended sources. We observe with the F770W, F1000W, F1130W, and F2100W filters and the FASTR1 readout pattern, yielding total exposure times per pointing of 88.8, 122.1, 310.8, and 321.9 s, respectively (Table 6). All of the galaxies are larger than the MIRI field of view, so we obtain a background measurement using this same MIRI imaging sequence for a blank region of sky near the galaxy in parallel with NIRCam, as just discussed.

We include the following special requirements in our observation plan: (1) a timing requirement (sequence observations, noninterruptible) that ensures the MIRI background imaging is obtained directly following the MIRI observation of the target and (2) position angle requirements to maximize the overlap of the primary MIRI observations with existing HST, MUSE, and ALMA coverage, while placing the parallel MIRI background observation on blank sky. In some of the galaxies with the largest angular extents in our sample, the parallel sky background measurement may not be completely beyond the outskirts of the galaxy disk. Finally, whenever possible, we attempted to widen the orient constraints to provide a minimum observation window of 2 weeks to provide greater opportunities for scheduling. It should be noted that the range in allowable orient angles can result in imaging that is not optimally aligned with previous data or astrophysically interesting features in the galaxy. For example, in NGC 1365, the observed orientation led to the omission of portions of the eastern spiral arm.

We start with the uncalibrated raw ("uncal") files obtained from the Mikulski Archive for Space Telescopes (MAST) and use Detector1Pipeline, which converts from "ramps" to "slopes" and applies basic detector-level corrections to the data, including the flagging of saturated pixels and corrections for chip persistence. ⁷⁴ We ran the Detector1Pipeline with the default parameters, with one exception. To attempt to recover some of the saturated pixels, we removed the restriction that during ramp fitting, at least two groups must not reach saturation. These saturated pixels are primarily an issue in the centers of NGC 1365 and NGC 7496, which both host bright active galactic nuclei (AGN). Removing the restriction did not recover the very central pixels coincident with the AGN themselves, but it did allow us to recover a number of the brighter pixels surrounding the central point source. However, some of these sources, especially the brightest compact, massive star-forming regions in the starburst ring of NGC 1365, still show artifacts that have structure suggesting that the source brightness remains underestimated by the imaging.

Stage two of the pipeline produces calibrated individual exposures. ⁷⁵ During this stage, the pixel coordinates are translated into WCS coordinates, which includes application of distortion corrections. This stage also involves flat-fielding and background subtraction (the background subtraction is used only for the MIRI data). Because of our parallel observing visit sequence (Section 3), the pipeline does not automatically associate the appropriate MIRI off (sky background) observations with the on-galaxy images. Instead, we manually implement this association during this processing step. For IC 5332, only the on-galaxy MIRI observations were obtained during our initial visit. Therefore, we used the off-galaxy observation from NGC 7496 for this purpose in IC 5332. This appears to work reasonably well but does result in a few visible imperfections in the background subtraction for IC 5332. Also for MIRI, we masked the part of the image associated with the Lyot coronagraph. In this early reduction, the coronagraphportions of the image consistently showed imperfect background subtraction that caused issues with our mosaicking.

After stage two, we applied additional processing to suppress 1/f noise in the NIRCam data. This correlated read noise primarily manifests as striping across the NIRCam images along detector rows (in the fast-read direction). This can severely degrade the quality of the images, particularly at shorter wavelengths. Schlawin et al. (2020) studied this effect using simulated data and provided a number of suggestions for how to deal with this 1/f noise. We experimented with simple row-by-row median subtraction but found that it was often unfeasible to calculate a robust background median in the small range covered by each amplifier due to the presence of widespread diffuse emission in most of our images. Moreover, we found that though this median subtraction produced visually better individual images, it yielded poor results when multiple images were combined to estimate 3.3 μm PAH maps (Sandstrom et al. 2023a). In those cases, the 1/f noise-driven stripes reappeared and were severe in the star-subtracted PAH images.

Instead, as shown in Figure 3, we had success modeling the 1/f noise using robust principal component analysis (PCA) based on the technique developed in Wild & Hewett (2005) and Budavári et al. (2009) for SDSS spectra. We treat each column of the data as an individual noise spectrum and fit the results for each amplifier separately because the noise properties differ between amplifiers. We mask out bright sources and then apply a Butterworth (1930) filter to remove large-scale structure but retain the small-scale noise properties. We also shuffle in the x-direction and apply a different offset to each column, effectively rolling the column along the y-axis by a random amount. These shifts prevent the PCA from simply learning where the mask is and responding to the location of bright emission. After this processing, we use PCA to fit 50 components. Then we reconstruct a model describing the 1/f noise-induced stripes for that column using the five components with the largest eigenvalues. Using only five components limits the processing time, and we found that increasing the number of components beyond five did not significantly improve the noise model. Finally, we subtract the model of the 1/f noise-induced stripes from the data.

Zoom In Zoom Out Reset image size

Figure 3 illustrates this process. The left image shows visible artifacts, mostly horizontal striping. The PCA model successfully captures these features without including the extended bright emission from the source. Occasionally, corrupted columns and other features arise, which the PCA also robustly handles. The right image shows visibly less striping and fewer artifacts than the original image. We find that the PCA reduces the rms noise in the images between around 10% and 30%, with the amount of reduction depending mainly on the amount of area filled by emission in the image. For images with more source-free sky background, the PCA produces a more significant decrease in the noise level.

The level three stage of the pipeline ⁷⁶ is intended to combine individual tiles to produce a final rectified mosaic. This is achieved through relative alignment between tiles, background matching, and the final drizzling to the output pixel grid.

In our early science reduction, this step involved the most changes from the default pipeline parameters. For MIRI, we ran the level three pipeline on individual tiles, i.e., the individual pointings, which we later mosaicked together by hand. For NIRCam, we used the JWST pipeline to stitch the tiles.

We adjusted the stage three pipeline as follows. First, in the tweakreg step, which provides relative alignment between tiles, we found it necessary to limit the "roundness" of sources used for alignment to −0.5 < roundness < 0.5 to reject diffraction spikes from stars and the numerous extended sources in our field. We also turned off the 2dhist option, as our relative corrections were generally ≲1'', and leaving this on degraded the quality of the alignment. We fit the per-tile offsets using only (x, y) shifts, as we found that allowing for rotation produced negligible rotation angles and degraded the quality of the fits. For data acquired in 2022 July and August, we found relative corrections as large as ∼1''. This primarily reflects the quality of the guide star catalog during that period. The size of the correction will likely decrease as JWST observations continue and the guide star catalog is updated and improved.

For the NIRCam data, we calculate and subtract the background level for each detector individually. By default, the pipeline will subtract a single background for the four short-wavelength NIRCam detectors; however, we found that this produced clear "steps" in our mosaics.

In the final redrizzling step during stage three, we produced mosaics with the conventional north-up orientation to allow for more straightforward comparison with data at other wavelengths.

The level three pipeline produced mosaics that have good relative alignment between NIRCam tiles, but the alignment of the JWST images relative to other data with well-established absolute astrometry tended to be poorer, as high as ∼1'' in some cases. After stage three processing, we therefore performed an additional correction to improve the absolute astrometry of the NIRCam and MIRI images. We adopt the general strategy of bootstrapping the astrometric solution from images that are close in wavelength instead of directly using astrometric catalogs, as was followed for PHANGS–HST (Lee et al. 2022).

Our JWST images lack enough Gaia (Gaia Collaboration et al. 2016) sources to anchor our astrometry directly to Gaia (Lindegren et al. 2018, 2018). However, our NIRCam images detect many asymptotic giant branch (AGB) stars that are also visible in the PHANGS–HST images, which have astrometric solutions accurate to better than 10 mas. Thus, we extracted a set of AGB stars from the PHANGS–HST DOLPHOT catalog (Lee et al. 2022; Thilker et al. 2022) and cross-matched them to sources detected in the NIRCam bands in the level three source_catalog step (the final level three step) using the XYXYMatch function in tweakwcs. Typically, we identified a few hundred AGB stars per target from the PHANGS–HST catalogs, and the pipeline found 50–100 good matches in the NIRCam imaging. We then solve for a linear offset and rotation to match the images (though the rotation was typically very small). After applying these solutions, comparisons between the NIRCam and HST data show that the rms scatter in the astrometric accuracy across the NIRCam images is about ±0.1 NIRCam pixel.

The MIRI images have fewer point sources than those from NIRCam and HST, so we derive the astrometric solution for MIRI using a cross-correlation approach. To do this, we take advantage of the fact that the F335M NIRCam band, which has been anchored using AGB stars to the PHANGS–HST and Gaia frame, contains both stellar and ISM-like emission because of the strong PAH band at 3.3 μm (see Figure 2). The stellar emission allows us to pin the astrometry to the AGB stars and the PHANGS–HST frame, while the morphology of the ISM at 3.3 μm resembles the ISM emission in the longer-wavelength MIRI bands. This allows us to use a cross-correlation analysis to solve for the relative astrometry of the NIRCam and MIRI bands. Specifically, we used the image-registration package to run a cross-correlation and solve for the linear shift (by identifying the maximum of the cross-correlation) needed to match the shortest MIRI band, F770W, to the F335M astrometry. Then, with the F770W astrometry established, we solved for the shift needed to align F1000W to F770W, F1130W to F1000W, and F2100W to F1130W. We adopted this stepwise approach to reflect the increasing dominance of the ISM emission over stellar emission as the wavelength increases. The cross-correlation technique generally produces "good" results but can be sensitive to artifacts, e.g., diffraction spikes around the AGN. By eye, this seems to yield images that are aligned at the level of ±1 MIRI pixel, i.e., about 01. While this is sufficient for the first results presented in this Issue (except where noted in specific papers), the alignment is not as good as would be expected based on alignment based on a more limited set of carefully selected, bright, isolated point sources (e.g., Lee et al. 2022). We expect to significantly improve this in future work.

After level three processing and alignment, we combined the individual MIRI tiles into a single image for each galaxy at each band. To do this, we defined one of the tiles for each galaxy as the reference tile. Then, we reprojected the other tiles onto the astrometric grid for the reference tile and compared the intensity of the two tiles in the region where they overlap. We solved for and applied the additive offset needed to yield a one-to-one match in intensities between the other tile and the reference file. These offsets were typically ≲0.05 MJy sr⁻¹ for F770W and F1000W, ≲0.15 MJy sr⁻¹ for F1130W, and as high as ∼0.3–0.8 MJy sr⁻¹ for F2100W. This yielded a set of individual tiles with the same background level but did not anchor that background level to the correct absolute value.

Next, we defined a new astrometric grid with a larger size but the same pixel scale and orientation as that of the reference field. We then reprojected all tiles onto this new grid, averaging intensities from different tiles with equal weighting where tiles overlapped. The result at this stage was a single combined MIRI image registered via cross-correlation to the NIRCam and other MIRI images but with an absolute background level that was still uncertain.

Note that we expect to deprecate most of these manual steps as the pipeline processing improves in the future. However, for this first processing during these early months, we found that this simple, by-hand approach yielded smoother backgrounds and a better match to previous infrared imaging of our targets than the pipeline.

As a final step, we adjust the overall background level of the MIRI images to match previous wide-field mid-infrared imaging of our targets. The JWST field of view for most of our targets contains relatively little empty sky, and the match between the off-source and on-source image was not sufficient to yield a precise background level. Fortunately, all of our sources have previous observations at least by the Wide-field Infrared Survey Explorer (WISE) at 12 and 22 μm (Wright et al. 2010) and for NGC 628, NGC 1365, and IC 5332 by Spitzer at 8 and 24 μm (Kennicutt et al. 2003; Armus et al. 2009; Dale et al. 2009). Though these data have much worse resolution than the JWST images, they cover a much larger area and extend to empty sky and so have well-established background levels.

We perform this background homogenization in two stages following a procedure detailed in the Appendices of Leroy et al. (2023). First, we establish a common background system across all JWST bands. Then, we anchor the overall system to an external band with a well-established background. To do this, we make a version of the MIRI images that all share the PSF of the F2100W image and a version of all JWST, Spitzer, and WISE data that share a common 15'' FWHM Gaussian PSF. That is, we match the backgrounds among the JWST bands at common resolution. These procedures leverage the fact that mid-infrared emission at different wavelengths shows strong, often nearly linear correlations even though the exact band ratios vary with filter combination (see the Appendix in Leroy et al. 2023). This procedure appears to yield good results with an uncertainty in the overall background level of better than ±0.1 MJy sr⁻¹ across all MIRI bands.

The reduction procedures described above yield data that are useful for the initial work presented in this Issue. Over the longer term, improvements will be made to address known limitations that we summarize here (as well as possible additional issues found during future analysis). In addition to refining our own procedures, we expect to take full advantage of the monthly refinements to the pipeline being released by STScI. ⁷⁷ The pipeline documentation already identifies some of these issues as areas of future development.

Where these limitations affect our first results scientific papers, they are noted, but most of our analyses are conducted to mitigate their impact. For example, aperture photometry with local background subtraction avoids uncertainties in measurement of the background (e.g., Rodriguez et al. 2023; Whitmore et al. 2023b).

4.8.1. Astrometric Uncertainties

4.8.2. Flux Calibration and Backgrounds

4.8.3. Saturation and Poor Ramp Fitting

4.8.4. Angular Resolution

While past infrared missions have provided imaging of the nearby galaxies in the PHANGS–JWST sample, JWST's vast improvement in sensitivity and resolution reveals an extraordinary new view of the complex organization of the ISM. In Figures 4 and 5, we present our JWST imaging of NGC 628 and NGC 7496 in selected filters alongside previous Spitzer observations to illustrate data quality and the new science enabled by the observations. The two galaxies are at 9.84 and 18.7 Mpc, respectively, and span most of the distance range (5–20 Mpc) covered by our sample. In Figure 2, we show the filter curves of the Spitzer bands that overlap JWST wavelength coverage together with the eight PHANGS–JWST filters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Perhaps the most remarkable feature in the PHANGS–JWST imaging is the network of filaments, shells, and bubbles and the star-forming populations nested in these structures, as traced by the dust in emission. While the highest surface brightness features are apparent in the Spitzer imaging, the JWST imaging shows the pervasiveness of the network deep into the interarm regions and morphological details that suggest that the multiscale impact of star formation feedback on the ISM is ubiquitous. With the PHANGS–JWST imaging, we will now be able to develop catalogs of these structures and characterize their ensemble properties (e.g., Hassani et al. 2023; Rodriguez et al. 2023; Thilker et al. 2023; Watkins et al. 2023) even with the lowest-resolution image at 21 μm.

For both NGC 628 and NGC 7496, we show the NIRCam F335M and Spitzer/IRAC 3.6 μm images. The NIRCam F335M filter better isolates the 3.3 μm PAH feature, and observations in the flanking F300M and F360M bands (not shown) enable continuum subtraction. It is clear that the increase in resolution (PSF FWHM 16 versus 011), which now samples 5 pc at NGC 628 and 10 pc at NGC 7496, enables the study of individual star clusters and associations and allows embedded stellar populations to be identified. A slightly extended point source southeast of the galactic center in NGC 628 is resolved into two background galaxies and illustrates how JWST significantly improves the characterization of extended stellar structures and the identification of background "interlopers."

For NGC 628, the MIRI F770W and Spitzer/IRAC 8 μm images are shown. The bandpasses, which are similar, capture the complex of PAH features in the ∼8 μm spectral region (Figure 2). Overall, the same structures are apparent in both images, with the MIRI F770W image being eight times sharper (PSF FWHM 20 versus 024), sampling 12 pc at NGC 628.

Finally, the MIRI F2100W and Spitzer/MIPS 24 μm images, which capture flux from the warm dust continuum, are shown for both galaxies. Here, although the relative increase in the angular resolution is similar to that for the previous two sets of images, the improvement in the fidelity of the JWST data is the most dramatic, as the MIPS 24 μm PSF was only 6'', and the majority of emission appeared diffuse. With the MIRI F2100W PSF of 067 (30 pc at NGC 628 and 60 pc at NGC 7496), the dust emission is resolved into compact sources, mostly along the spiral arms, and shows the same intricate ISM structure seen in the MIRI F770W image.

In terms of surface brightness, the Spitzer images achieve slightly better sensitivities than the JWST data, as might be expected. The NGC 628 images from the Spitzer SINGS program (Kennicutt, Jr. et al. 2003) are mosaics with 1σ surface brightness limits of ∼0.02 MJy sr⁻¹ at 3.6 μm and ∼0.2 MJy sr⁻¹ at 24 μm. The corresponding sensitivities are 0.04 MJy sr⁻¹ for NIRCam F335M and 0.25 MJy sr⁻¹ for MIRI F2100W (see Table 4 and discussion in the next section). Of course, the principal gain is in the ∼tenfold improvement in resolution and the corresponding decrease in the solid angle subtended by the PSF, so that the JWST point-source sensitivity is ∼50 times better than Spitzer.

Table 4 provides sensitivity limits based on our initial reduction of PHANGS–JWST data. We evaluated the depth of our data with respect to detection of point sources at the resolution achieved in each filter and measurement of diffuse emission on a fixed larger scale.

For the point-source limit, we first identified the faintest area mapped that was also fully contained within the footprint of maximum exposure time. Inside each area, aperture photometry was performed for a set of 50 independent circular apertures sized to match the 50% encircled energy radius of each filter. Background estimation for each aperture was also filter-dependent, using an annulus spanning a radial range from two to three times the aperture radius. The sigma-clipped standard deviation of the background-subtracted flux for these "empty" apertures was then scaled by a factor of 10 to estimate the 5σ point-source detection limit in each filter (5σ_f ), which accounts for a factor of 2 for the choice of 50% enclosed energy radius and a factor of 5 for the choice of significance level. No correction was made for the contribution of the point-source flux to the background annulus flux, which is expected to be negligible for an unresolved source. Averages of these galaxy-specific limits per filter are tabulated. It is important to note that this method does not include loss of point-source sensitivity due to source crowding in more complex, brighter regions of each JWST target. Such loss of completeness due to crowding and the associated photometric bias will be quantified later when PSF-fitting photometry is undertaken for our entire set of targets.

To estimate the surface brightness sensitivity (σ_I ) for the imaging data, we consider the distribution of pixel values in an unsharp masked image. For each image, we compute a median filter with a circular top-hat filter with a diameter six times the PSF FWHM in the band. We subtract off the median smoothed version from the original map to remove large-scale structure and produce an unsharp masked image. We then use a median absolute deviation estimator to measure the local noise in the unsharp masked image. We iteratively reject all values larger than the local 3σ noise measurement in each filter. After the data rejection converges, the local noise estimates are smoothed with a top-hat filter that has a size of 15× the PSF FWHM and interpolating over missing values. We compare this local estimate of the noise to the noise estimates provided by the JWST imaging pipeline. The resulting noise maps are within 10% of the pipeline-generated noise maps. We also estimated the local noise by calculating the median autocorrelation of the data as a function of scale and examining the distributions of the data in signal-free parts of the images (though some bands show emission across the field). All methods arrive at similar noise estimates and reflect the spatial structure in the mosaicked fields.

Note that because of the extended JWST PSFs, there is less encircled light at a PSF FWHM than for a Gaussian source. Thus, the simple scaling of surface brightness noise to obtain a point-source sensitivity leads to an overestimate of the point-source sensitivity: σ_I Ω_PSF ≳ σ_f , where Ω_PSF is the solid angle subtended by the PSF.

Since the first observations were taken for the PHANGS–JWST Treasury survey in 2022 July, color composites that combine JWST imaging from multiple filters have been developed and released for all of our targets by astronomers, professional public outreach teams, the press, and the public. The color images have broadly captured the attention of both the science community and general public, not only because the structures portrayed are aesthetically stunning but also because the images vividly illustrate and make qualitatively accessible the physics of star formation, feedback, and the ISM, particularly when the UV-optical imaging from HST is included.

Figure 6 shows the first composite PHANGS–JWST image published for NGC 7496, which was observed shortly after the end of commissioning and among the first science data to be released from the mission. The image is a sum of two composites: one based on the HST UV-optical filters (red: F814W/F555W/F438W; green: F336W; blue: F275W) and a second based on the JWST MIRI F1000W, F1130W, and F2100W filters in red hue. The animation compares the emission from the HST and JWST data. The regions and structures that are dark in the optical due to dust obscuration are illuminated in the mid-infrared by the reradiated emission from small dust grains. The complex network of filaments, bubbles, shells, and compact sources described earlier can now be seen in the context of the visible young stellar populations that line the peripheries of the network and have provided the feedback energy that, together with galactic dynamics, shapes the ISM. Infrared compact sources without optical counterparts can be identified, revealing the sites of the earliest stages of star formation.

Figures 7 and 8 show RGB color composites, now based on only the MIRI imaging, for all four PHANGS–JWST first results galaxies. The diversity of the dust emission structure motivates the study of a representative sample of galaxies and our specific science goals (next section). The four galaxies bracket the full range of star formation properties in the sample (with NGC 1365 having the highest star formation activity and gas surface density and IC 5332 among the lowest) and feature varied dynamical structures (e.g., bars in NGC 1365 and NGC 7496, strong spiral arms in NGC 628, and weak spiral structure in IC 5332). The PHANGS–JWST team is collaborating with image processor Judy Schmidt to produce color composites as new data are received and with the public outreach offices at ESA and STScI to support the communication of new results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

PHANGS–JWST near- and mid-IR imaging will be used to detect dust-embedded clusters and associations, mitigate and measure the effects of extinction, and constrain other physical properties such as age and mass through SED fitting (e.g., Figure 9). We will test different selection strategies to systematically identify embedded cluster candidates and produce the first complete catalogs of these sites of the earliest stages of star formation. The catalogs will enable us to measure the fraction of stars that are born in compact clusters, characterize the distribution of cluster masses at birth, and investigate dependence on environment.

The ensemble properties of star clusters are highly debated in the literature (e.g., Larsen 2009; Chandar et al. 2017; Krumholz et al. 2019; Adamo et al. 2020a, 2020b). Due to dust extinction, it has been difficult to obtain a truly complete census of clusters across all galactic environments and ages. The PHANGS–JWST infrared imaging will now allow us to study the clusters in dusty galaxy centers, spiral arms, and bar-end regions. These are sites where super star clusters tend to be found, the cluster formation efficiency is expected to be highest (e.g., Kruijssen 2012; Krumholz 2014; Adamo et al. 2015; Ginsburg & Kruijssen 2018; Grudić et al. 2021), and molecular cloud populations show clear enhancements in both density and turbulence (e.g., Longmore et al. 2014; Beuther et al. 2017; Sun et al. 2020; Henshaw et al. 2022). Extrapolating from the visible cluster population, we expect to detect ≳100 new embedded clusters in each galaxy, or ∼2000 across the sample, with the greatest gain in the dustiest and most active parts of the sample. Analysis combining JWST with HST data will also allow us to break the age-reddening degeneracy for optically detected clusters and derive robust constraints on age, mass, and reddening for both embedded and visible cluster populations (e.g., Figure 2).

These investigations use the full set of PHANGS–JWST data for measuring the properties of embedded clusters. In particular, the near-IR bands provide key constraints on stellar properties. The F200W observations provide good visibility through the dust in the region of the SED still dominated by stellar photospheric emission. The excellent angular resolution (007) at 2 μm is required to resolve star clusters, which have effective radii between 0.5 and about 10 pc (Portegies Zwart et al. 2010; Ryon et al. 2017).

Rodriguez et al. (2023) and Whitmore et al. (2023b) took the first steps in identifying embedded clusters, estimating their physical properties, and examining the incompleteness in optical samples. Hoyer et al. (2023) studied the nuclear star cluster in NGC 628 and the nature of the mid-infrared structure surrounding it. Schinnerer et al. (2023) showed that massive, dusty clusters can provide a sensitive diagnostic of the very dusty star-forming environments in nuclear regions when comparing ALMA data and simulations. Some of the clusters detected in the PHANGS–JWST survey push current stellar and dust models used by CIGALE to their limits due to their very young ages, high dust content, and/or low mass. Including dust models computed in the optically thick case and stellar populations with a stochastically sampled initial mass function may prove necessary to model such clusters.

PHANGS–JWST data enable characterization of the poorly studied embedded phase of star formation before feedback disrupts the natal cloud. Observations that resolve galaxies into discrete star-forming regions, clusters, and molecular clouds have enabled statistical measurements of the timescales associated with the major phases of the star formation cycle (e.g., Kawamura et al. 2009; Kruijssen & Longmore 2014; Meidt et al. 2015; Corbelli et al. 2017; Grasha et al. 2018; Kruijssen et al. 2018, 2019a; Grasha et al. 2019; Hannon et al. 2019, 2022; Chevance et al. 2020a; Kim et al. 2022; Ward et al. 2022). Over the last decade, the (anti)correlation between clouds and H ii regions or optically identified clusters has constrained the total lifetimes of clouds to 10–30 Myr (where the range represents the physical variation rather than the uncertainty; see, e.g., Chevance et al. 2020b, 2022a). However, both the duration of the earliest, embedded phase and the length of time that a cloud remains truly devoid of star formation activity remain poorly known.

Measuring these timescales is essential for understanding the dominant physics of star formation. The time between cloud formation and the onset of star formation constrains the (im)balance between the driving and decay of turbulence (e.g., Gnedin et al. 2016; Padoan et al. 2017). The duration of the embedded phase reflects the time that it takes feedback to start dispersing a cloud (e.g., Olivier et al. 2021). Measuring the embedded timescale can also distinguish between models in which star and cluster formation "accelerate" over a cloud's lifetime (Murray 2011; Lee et al. 2016) and those that favor a more steady process (Krumholz & McKee 2005; Krumholz & Tan 2007; Krumholz & McKee 2020). Timescales for the embedded and optically visible stages can be used to infer which feedback mechanisms quench star formation (e.g., Lopez et al. 2014; Kruijssen et al. 2019a; Kim et al. 2022; Chevance et al. 2022b; Hannon et al. 2022) and whether clouds are magnetically sub- or supercritical (Kim et al. 2021b).

Accessing the embedded phase of star formation requires a tracer of obscured star formation that reaches the scale of individual clouds. The limited resolution of Spitzer (∼6'') limited the number of galaxies where such measurements were possible to five (Kim et al. 2021a). With the increased spatial resolution and sensitivity in the mid-infrared achieved with MIRI, the time evolution of the early stages of star formation can now be characterized in galaxies out to ∼30 Mpc.

By combining some of the first PHANGS–JWST MIRI observations with PHANGS–ALMA CO observations (Leroy et al. 2021), Kim et al. (2023) determined the relative duration of the early phases of star formation, from truly quiescent clouds, to embedded star formation, revealed star formation, and cold gas-free H ii regions in NGC 628. In combination with the stellar "clocks" offered by star clusters (see Section 6.1) and Hα emission (e.g., Leroy et al. 2013; Haydon et al. 2020), these relative measurements are converted into absolute timescales. PHANGS–JWST data will enable this measurement across the wide range of environments that are present in the sample. These observations will constrain how the physical mechanisms that regulate and drive the successive phases of star formation map to the environment.

Most of the imaging bands used in PHANGS–JWST show emission features from dust continuum and PAHs. The PAHs reprocess as much as 20% of all absorbed stellar UV/optical light (Smith et al. 2007). As they do so, they eject electrons that heat the ISM and produce bright mid-IR emission features (e.g., Tielens 2008) that are among the most highly detectable star formation tracers across cosmic time (e.g., Riechers et al. 2014). However, PAH properties and abundance vary substantially within and between galaxies. Fewer PAHs are present in H ii regions or at low metallicities (e.g., Engelbracht et al. 2005; Chastenet et al. 2019), but they can make up >5% of the dust mass in the diffuse ISM of solar metallicity galaxies (Draine et al. 2007; Aniano et al. 2020). This variability hampers the use of PAH emission as a star formation tracer (Calzetti et al. 2007; Whitcomb et al. 2020) and leads to critical uncertainties when modeling their role in the ISM. The PHANGS–JWST Treasury data will yield the first highly resolved atlas of PAH properties—abundance, size distribution, and charge inferred from mid-IR band ratios—in the nearby galaxy population. These data can be used to build an empirical model to predict the abundance and characteristics of PAHs as a function of local metallicity, ISM phase, density, and radiation field (from HST, MUSE, and ALMA). The measurements underpinning this model must be highly resolved because radiation field, density, and ISM phase can vary by factors of ∼1000 on 10–100 pc scales, altering the PAHs (Chastenet et al. 2019; Draine et al. 2021). The PHANGS–JWST first results in this Issue from Chastenet et al. (2023a, 2023b), Dale et al. (2023), and Egorov et al. (2023) all show that PAH emission properties do indeed vary systematically with the local galactic environment.

The PAH emission offers a promising tracer for the structure of the neutral ISM. Traditionally, H i 21 cm emission and the lowest-energy rotational lines of CO (λ = 2.6 and 1.3 mm) have been used to map out the atomic and molecular phases of the ISM. These tracers suffer some limitations; the very long wavelength of 21 cm emission makes it intrinsically faint, and current interferometers can typically only achieve resolutions of 5''–10'' for nonprohibitive integration times. The (sub)millimeter CO emission lines are at shorter wavelengths, and thus a significantly better resolution of 1'' can be achieved with ALMA (Leroy et al. 2021), but they are an imperfect proxy for the H₂ that constitutes the molecular ISM. With JWST, the mid-IR PAH emission potentially provides a much higher-resolution view (<037) of the same interstellar gas that is traced by H i and CO. While the emissivity of PAHs likely varies with environment (Draine et al. 2007; Hensley & Draine 2022), the spatial distribution of the mid-IR PAH emission yields an exceptionally detailed, high dynamic range image of the structure of neutral gas in a galaxy. Based on studies in the Milky Way and the nearest galaxies, the older "spherical cloud" view of the ISM has been integrated with a newer filamentary description of ISM structure (Hacar et al. 2022). The enthusiasm for such a filamentary view has been driven by infrared imaging by the previous generations of instruments, but it was difficult to extend this view to a galactic scale in all but the nearest galaxies. The JWST imaging of PAH emission yields such high-resolution maps of neutral gas throughout the Local Volume.

Figure 7 shows the filamentary organization of the gas in NGC 628, which prompted initial investigations of ISM structure using PAHs as ISM tracers (e.g., Sandstrom et al. 2023b; Barnes et al. 2023; Watkins et al. 2023). Sandstrom et al. (2023b) demonstrated a new JWST band combination that produces high-quality maps of the PAH feature at 3.3 μm with excellent stellar continuum subtraction. These maps offer the highest-resolution (011) wide-area maps of the ISM in nearby galaxies. Sandstrom et al. (2023b) verified that PAH features are intimately connected to the neutral ISM and used the genus statistic for characterizing the emission topology. Thilker et al. (2023) showed how the filamentary network of 3.3 μm PAH emission relates to extinction lanes identified in the optical HST imaging and to the young cluster population of a galaxy. Meidt et al. (2023) showed that the spatial separation scales between filaments is consistent with large-scale gravitational fragmentation models. Williams et al. (2022) leveraged the high-resolution images to make a careful measurement of evolutionary timescales in a spiral arm spur.

The PAH imaging also shows us where the neutral ISM is not found. Rings of PAH emission with faint centers can be used to define the wall of "bubbles" or feedback-driven shells. While bubbles from feedback have been long studied (e.g., Oey & Clarke 1997; Weisz et al. 2009, and references therein), the exquisite resolution, sensitivity to structure on many scales, and wide dynamic range of the PHANGS–JWST images make them an excellent tracer of bubble walls. PHANGS–JWST data can be used to produce inventories of bubbles and shell features. These structures can be connected to the sources of feedback identified from PHANGS–MUSE and PHANGS–HST catalogs of star clusters and associations and to the kinematics of the cold ISM as traced by PHANGS–ALMA. As an early demonstration of such analyses, Watkins et al. (2023) created a bubble catalog for NGC 628 and examined the properties, including the size distribution for the ensemble population. Barnes et al. (2023) demonstrated the capabilities of the PHANGS multiwavelength data set for measuring the expansion rates and understanding the origins of large kiloparsec-scale cavities, focusing on the "Phantom Void" in the southeast quadrant of NGC 628.

Mid-IR emission from hot dust is widely used to trace star formation at both low and high z (e.g., Kennicutt & Evans 2012). Unfortunately, the calibration of mid-IR star formation tracers has a strong dependence on both spatial scale and environment (e.g., Boquien et al. 2015), which likely stems from a variable heating contribution by old stars and uncertain dust geometry (e.g., Groves et al. 2012; Nersesian et al. 2019). PHANGS–JWST imaging, combined with supporting PHANGS data sets, enables a data-driven description of the origin and variations in mid-IR emission. Mid-IR emission is due to the heating of dust and PAHs by both young and old stars, but the bright thermal dust emission in the mid-IR continuum bands (F2100W and F1000W) is expected to respond differently to the local dust-heating radiation field than the PAH emission (Draine et al. 2007; Hensley & Draine 2022), offering a potentially powerful diagnostic of local interstellar conditions. With the insights obtained by comparing mid-IR band ratios and physical dust models, we will build on previous work that calibrated the mid-IR emission as an SFR indicator on galaxy-wide and subkiloparsec scales (e.g., Calzetti et al. 2007).

Combining robust high-resolution measurements of the embedded star formation from PHANGS–JWST with PHANGS–ALMA and PHANGS–MUSE imaging of the molecular gas and optically visible tracers of star formation and feedback will enable the first cloud-scale inventory of all stages in the star formation cycle across the 19 galaxies in our sample. Previous work on ≲100 pc scales suggested that up to ∼60% of molecular gas in galaxies may be devoid of Hα emission (e.g., Onodera et al. 2010; Schruba et al. 2010; Kreckel et al. 2018; Kruijssen et al. 2019a; Schinnerer et al. 2019b; Chevance et al. 2020a; Kim et al. 2022; Pan et al. 2022). Yet highly resolved IR observations show that truly non-star-forming clouds are scarce in our galaxy (Mooney & Solomon 1988; Ginsburg et al. 2012), and Spitzer results suggest that molecular clouds in Local Group galaxies spend up to half their lifetime hosting only embedded massive star formation (Kim et al. 2021a). With PHANGS–JWST, we will confirm whether there is a large reservoir of molecular gas without star formation in nearby star-forming disk galaxies and, if so, whether it is due to a delayed onset (e.g., Gnedin et al. 2016; Padoan et al. 2017) or sustained suppression (e.g., Federrath et al. 2016; Semenov et al. 2017; Meidt et al. 2018, 2020; Kruijssen et al. 2019b) of star formation.

PHANGS–JWST high-resolution measurements will also enable important empirical tests of cloud-scale models of star formation. These models predict the integrated and per freefall time SFR and efficiency of individual clouds based on their density, turbulence, and dynamical state (e.g., Padoan et al. 2017; Kim et al. 2021b), quantities that have already been measured using the PHANGS–ALMA data (e.g., Sun et al. 2020, 2022). Cloud-scale star formation models are widely used to interpret observations (e.g., Barnes et al. 2017; Utomo et al. 2018) and implemented as subgrid models in galaxy simulations (e.g., Gensior et al. 2020; Kretschmer & Teyssier 2020). The models have never been rigorously tested beyond the Local Group, however, since the rapid action of feedback prohibits linking optical SFR estimates to the progenitor cloud's gas mass and hence the integrated SFR per unit gas mass. Lacking cloud-by-cloud estimates of the embedded star formation activity, studies in external galaxies have resorted to using population averages, obtaining results that suggest—but cannot currently prove—important inconsistencies with the models, such as inferred star formation efficiencies that decrease with the gravitational boundedness of the gas (e.g., Leroy et al. 2017; Schruba et al. 2019). The new PHANGS–JWST data are essential for unambiguous cloud-scale tests of these predictions.

The PHANGS–JWST papers in this Issue have begun to explore these questions. Hassani et al. (2023) studied compact sources in F2100W, finding that nearly all sources are compact star-forming regions. These regions show strong associations with both the molecular gas and Hα emission. Liu et al. (2023) evaluated how strong heating in the IR from dusty star-forming regions is connected to the density and temperature of the molecular ISM as seen by ALMA. Leroy et al. (2023) used a large set of archival data to show that the mid-IR has an excellent capacity to trace the properties of the star-forming ISM. This retrospective analysis provides context for the high-resolution view provided by PHANGS–JWST. Using the new data, Leroy et al. (2023) conducted a correlation analysis of the mid-IR data, showing strong correlations with both the Hα and the CO in galaxies, and these correlations are stronger than the correlation between Hα and CO at high resolution.