PHANGS–JWST First Results: Mid-infrared Emission Traces Both Gas Column Density and Heating at 100 pc Scales

We consider mid-IR emission from the four MIRI bands imaged by PHANGS–JWST, F770W (7.7μm), F1000W (10μm), F1130W (11.3μm), and F2100W (21μm). To first order, the F770W and F1130W bands are expected to be dominated by emission from strong PAH bands while the F2100W filter covers mainly continuum emission from dust (e.g., see spectra in Draine & Li 2007; Smith et al. 2007; Galliano et al. 2018; Hensley & Draine 2022; Lai & Armus 2022). As discussed in Section 4.4, the situation with F1000W appears more ambiguous. Though expected to reflect primarily continuum, its behavior mirrors the PAH-tracing bands, suggesting that either (1) weaker PAH features or the wings of the nearby strong PAH features contribute to the band or (2) that the small grains contributing the emission mimic PAHs in many respects. Silicate absorption near 10 μm may also be important at high column densities, but should make only a minimal contribution for the column densities that cover most of the area in our targets.

Our expectations and calculation of the background levels (Appendices A and B) also make frequent reference to Spitzer's 8 μm and 24 μm bands and WISE's 12 μm band (Werner et al. 2004; Wright et al. 2010). The 8 μm band is dominated by the same PAH band as F770W, while the 24 μm band should reflect mainly continuum emission. The very broad WISE 12 μm filter integrates over both PAH features and the continuum. In SINGS H ii regions with mid-IR spectroscopy, Whitcomb et al. (2022) find that ≲50% of the overall emission in the WISE 12 μm band arises from PAH features, but this value might be somewhat higher in more diffuse regions.

Table 1 reports a series of typical band ratios for the first four PHANGS–JWST target galaxies, which can be used to translate between bands. Appendix A describes how we estimate these ratios from comparing our images to one another and to previous mid-IR imaging of our targets at 17 and 15'' resolution. We quote the ratios in three ways:

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal R }}_{7.7}^{X} & \equiv & \left(\displaystyle \frac{{I}_{\nu }^{X}}{{I}_{\nu }^{{\rm{F}}770{\rm{W}}}}\right)\,,\\ {{ \mathcal R }}_{24}^{X} & \equiv & \left(\displaystyle \frac{{I}_{\nu }^{X}}{{I}_{\nu }^{\mathrm{MIPS}24}}\right)\,,\\ {{ \mathcal R }}_{X}^{Y} & = & \displaystyle \frac{{{ \mathcal R }}_{7.7}^{Y}}{{{ \mathcal R }}_{7.7}^{X}},\end{array}\end{eqnarray} \tag{ 1 }$$

where X and Y denote some mid-IR bands, and we work with I_ν in units of MJy sr⁻¹. The two sets of ratios in Table 1, ${{ \mathcal R }}_{7.7}^{X}$ and ${{ \mathcal R }}_{24}^{X}$, are self-consistent results of a single calculation, not independent. We quote both because F770W is our highest-resolution image and because a large amount of pre-JWST work has centered on Spitzer's MIPS 24 μm band. The final line of Equation (1) simply notes that Table 1 can be used to express any pair of bands that we consider.

Table 1. Typical Observed Band Ratios in the First Four PHANGS–JWST Targets

Band	${{ \mathcal R }}_{7.7}^{X}$	${{ \mathcal R }}_{24}^{X}$	σa
			(dex)
F770W	1.0	1.31	⋯ b
F1000W	0.38	0.49	0.06
F1130W	1.35	1.76	0.05
F2100W	0.61	0.80	0.11
WISE3 12 μm	1.57	2.06	0.06
WISE4 22 μm	0.80	1.06	0.20
IRAC4 8 μm	1.0	1.31	0.04
MIPS24 24 μm	0.77	1.0	0.14

Notes. Observed ratios estimated from band–band comparisons in our first four targets as described in Appendix A and Figure 15. See definitions of ${ \mathcal R }$ in Equation (1). Note that ${{ \mathcal R }}_{24}^{X}$ are scaled versions of ${R}_{7.7}^{X}$ normalized to Spitzer at 24 μm.

^a Robustly estimated scatter in the log of the measured ratio between the listed band and F770W, measured at 15'' resolution (i.e., scatter in the blue points in Figure 15). ^b In Appendix A we use F770W as the reference band so no scatter is defined.

Download table as:  ASCIITypeset image

In framing our expectations, we will assume that all bands can be simply linearly converted to one another, which appears reasonable to first order based on Appendix A. To second order, the ratios among these bands do vary, especially in bright, high-intensity regions like massive H ii regions or the starburst ring at the center of our target NGC 1365. These ratio variations offer clues to how the physical properties and abundance of the dust grains producing mid-IR emission change as a function of environment. Measuring and interpreting these variations represent a key topic of other papers in this series (Chastenet 2023a, 2023b; Dale et al. 2023; Egorov 2023; Sandstrom 2023a). For our work, the most relevant effect will be that PAHs appear to be selectively destroyed in regions of intense star formation. This should suppress the PAH-dominated F770W and F1130W relative to the other bands in regions of intense extinction-corrected Hα emission. This effect is indeed present in our data and seen as a secondary correlation in Sections 4.2, 4.4, and 5.2. Perhaps surprisingly, F1000W shows almost identical behavior to F770W and F1130W in this regard, so that to some extent, those three bands contrast with F2100W in overall behavior (Section 4.4).

We work with CO (2–1) intensities in units of K km s⁻¹, I_{CO 2−1}. For a typical R₂₁ ≡ I_{CO 2−1}/I_{CO 1−0} = 0.65 (den Brok et al. 2021; Yajima et al. 2021; Leroy et al. 2022) and a Galactic CO (1–0)-to-H₂ conversion factor ${\alpha }_{\mathrm{CO}\ 1-0}^{\mathrm{MW}}=4.35$ M_⊙ pc⁻² (K km s⁻¹)⁻¹ (e.g., Bolatto et al. 2013),

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{mol}}}{{M}_{\odot }\,{\mathrm{pc}}^{-2}} & \approx & 6.7\left(\displaystyle \frac{{\alpha }_{\mathrm{CO}\ 1-0}}{{\alpha }_{\mathrm{CO}\ 1-0}^{\mathrm{MW}}}\right)\left(\displaystyle \frac{0.65}{{R}_{21}}\right)\left(\displaystyle \frac{{I}_{\mathrm{CO}\,2-1}}{{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)\\ \displaystyle \frac{N({\rm{H}})}{{10}^{20}\,{\mathrm{cm}}^{-2}} & \approx & 6.2\left(\displaystyle \frac{{\alpha }_{\mathrm{CO}\ 1-0}}{{\alpha }_{\mathrm{CO}\ 1-0}^{\mathrm{MW}}}\right)\left(\displaystyle \frac{0.65}{{R}_{21}}\right)\left(\displaystyle \frac{{I}_{\mathrm{CO}\,2-1}}{{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right),\end{array}\end{eqnarray} \tag{ 2 }$$

where Σ_mol is the molecular gas mass surface density including a factor of 1.4 to account for helium and heavy elements and N(H) is the total hydrogen column density, with no helium included. We phrase Equation (2) in terms of N(H) not N(H₂) = 0.5N(H) because dust emission is frequently normalized to N(H).

We also work with extinction-corrected Hα intensity in units ⁴⁹ of erg s⁻¹ cm⁻² sr⁻¹. We provide reference conversions adopting the Murphy et al. (2011) conversion from Hα luminosity to SFR, which is almost identical to the value advocated by Calzetti et al. (2007). In intensity units,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{SFR}}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-2}} & \approx & 642\,\left(\displaystyle \frac{{I}_{{\rm{H}}\alpha }}{\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{sr}}^{-1}}\right)\\ \displaystyle \frac{{\rm{\Sigma }}\left({Q}_{0}\right)}{{{\rm{s}}}^{-1}\,{\mathrm{pc}}^{-2}} & \approx & 8.8\times {10}^{49}\left(\displaystyle \frac{{I}_{{\rm{H}}\alpha }}{\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{sr}}^{-1}}\right)\\ \displaystyle \frac{{\rm{\Sigma }}\left({M}_{\mathrm{ZAMS}}\right)}{{M}_{\odot }\,{\mathrm{pc}}^{-2}} & \approx & 2200\left(\displaystyle \frac{{I}_{{\rm{H}}\alpha }}{\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{sr}}^{-1}}\right),\end{array}\end{eqnarray} \tag{ 3 }$$

where Σ_SFR is the SFR per unit area for a Kroupa IMF truncated at 100 M_⊙. In the second version, Q₀ refers to the production rate of hydrogen-ionizing photons with h ν > 13.6 eV and I_Hα .

The ionizing photon production rate may be more concretely related to the recombination line emission than the somewhat abstract SFR, but we also note that at our resolution and since we are studying whole galaxies, Hα can more accurately be thought of as being driven by the number of local ionizations or recombinations than the production rate. The final expression provides a reference conversion to the mass surface density of zero-age main-sequence stars needed to produce that density of ionizing photons for a fiducial Starburst99 run (Leitherer et al. 1999). See Belfiore et al. (2022a) for more discussion on the diffuse ionized medium and Belfiore et al. (2022b) and Murphy et al. (2011) for translations among different SFR tracers.

Throughout this paper, we work with Hα emission corrected for the effects of extinction using the Balmer decrement. This approach combines observations of Hα and Hβ with an assumed wavelength-dependent extinction (e.g., see Osterbrock & Miller 1989). It has a long history (e.g., Miller & Mathews 1972), and Balmer decrement-corrected Hα measurements, including from SDSS, SINGS, CALIFA, and MaNGA, underpin much of our knowledge of extragalactic SFRs, including the calibration of IR-based SFR tracers (e.g., Brinchmann et al. 2004; Moustakas et al. 2006, 2010; Kennicutt & Hao 2009; Hao et al. 2011; Catalán-Torrecilla et al. 2015; Belfiore et al. 2022b). Although the Balmer decrement has limitations, particularly in dense, high-extinction, mixed media (e.g., Melnick 1979; Lequeux et al. 1981; Wong & Blitz 2002), at the moderate Σ_SFR and moderate extinction that we study here, we have every expectation that the Balmer decrement works well. This expectation appears to be borne out by numerical simulations (Tacchella et al. 2022), and studying SINGS Prescott et al. (2007) found no evidence for a substantial highly obscured population pervading normal star-forming disk galaxies. Indeed, in the brightest regions of PHANGS–MUSE (Belfiore et al. 2022b; F. Belfiore et al., in preparation) find good agreement between the estimates of Hα extinction based on the Balmer decrement and those obtained from "hybrid" Hα+IR prescriptions calibrated by Calzetti et al. (2007) to match extinction-corrected Paschen α. Our results in Section 4.2 comparing the extinction-corrected Hα with JWST mid-IR emission also find no evidence that the Balmer decrement significantly underestimates the extinction (see also Hassani et al. 2023, in this Issue).

Even in the presence of a relatively weak radiation field, dust in the ISM should produce mid-IR emission, primarily through stochastic, single-photon heating. The intensity of emission depends on the column density of dust and the intensity of the illuminating radiation field. For relatively weak radiation fields, a reasonable fiducial expectation based on the Draine & Li (2007) dust models and also following Compiegne et al. (2010) is

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{I}_{\nu }^{X}}{\mathrm{MJy}\,{\mathrm{sr}}^{-1}} & \approx & 0.022\,{{ \mathcal R }}_{24}^{X}\,\left(\displaystyle \frac{N({\rm{H}})}{{10}^{20}\,{\mathrm{cm}}^{-2}}\right)\\ & & \times \left(\displaystyle \frac{D/G}{0.01}\right)\left(\displaystyle \frac{U}{{U}_{0}}\right)\\ & \approx & 0.020\,{{ \mathcal R }}_{24}^{X}\,\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{gas}}}{1\,{M}_{\odot }\,{\mathrm{pc}}^{-2}}\right)\\ & & \times \left(\displaystyle \frac{D/G}{0.01}\right)\left(\displaystyle \frac{U}{{U}_{0}}\right),\end{array}\end{eqnarray} \tag{ 4 }$$

where ${I}_{\nu }^{X}$ is the mid-IR intensity at band X, N(H) is the line-of-sight column density of hydrogen, and Σ_gas is the gas mass surface density. The U/U₀ term expresses the mean local ISRF illuminating the dust in units of the solar neighborhood Mathis et al. (1983) field adopted by Draine & Li (2007). The formulae apply at 24 μm but the factor ${{ \mathcal R }}_{24}^{X}$ can be drawn from Table 1 to express a prediction for any of the mid-IR bands of interest, X.

In Equation (4), the factor D/G is the dust-to-gas mass ratio, where a value of 0.01 represents a typical result applying the Draine & Li (2007) model to star-forming massive galaxies (e.g., Draine et al. 2007; Sandstrom et al. 2013; Aniano et al. 2020). When considering a PAH-dominated band, one should expand this part of the equation with an additional term proportional to q_PAH, the fraction of the dust mass in PAHs. This allows one to account for variations in the dust composition that may make PAHs either more or less abundant. Draine et al. (2007) and Draine & Li (2007) note q_PAH ≈ 0.046 as appropriate for Milky Way–like galaxies, and adding a factor

$$\begin{eqnarray*}&&\times \,\left(\displaystyle \frac{{q}_{\mathrm{PAH}}}{0.046}\right)\end{eqnarray*}$$

to Equation (4) offers a reasonable first-order approach at, e.g., F770W, F1130W, or 8 μm.

For U ∼ U₀, dust emissivity in the mid-IR depends approximately linearly on U in the Draine & Li (2007) model (see their Figure 13), reflecting the fact that U sets the rate of stochastic, single-photon heating. ⁵⁰ Because the mass surface density of the dust is given by the product Σ_dust = D/G × Σ_gas, Equation (4) simply predicts that for low-intensity, mid-IR emission tracks the dust column times the illuminating radiation field.

The Σ_gas in the second line of Equation (4) includes helium, but N(H) does not. In general, Σ_gas includes both atomic and molecular phases of the ISM, but for much of this paper, we will focus on molecular-gas-rich regions and will approximate Σ_gas ≈ Σ_mol. The ability to generalize a relationship calibrated in the molecular-gas-dominated ISM to the atomic-gas-dominated regime will depend on how D/G, q_PAH, and U vary as a function of ISM phase or density. This is discussed in Sandstrom (2023b) and in Section 7.1 of this paper.

The mid-IR is also often used to trace star formation. In massive star-forming galaxies, dust absorbs and then re-emits much of the radiation from massive young stars. Some of this emission emerges in the mid-IR, which indicates the UV heating of dust grains. Based on observed excellent correlations between recombination line and mid-IR emission, the mid-IR has become a widely used SFR indicator. A variety of calibrations exist in the literature. Almost all of them have a fundamentally empirical calibration (e.g., see Murphy et al. 2011; Leroy et al. 2012; Kennicutt & Evans 2012; Calzetti 2013). We note here a common formulation for this relation for emission at λ = 24 μm:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{SFR}}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}{\mathrm{kpc}}^{-2}} & \approx & 2.97\times {10}^{-3}{\left({{ \mathcal R }}_{24}^{X}\right)}^{-1}\\ & \times & \left(\displaystyle \frac{{C}_{24}}{{10}^{-42.7}}\right)\left(\displaystyle \frac{{I}_{\nu }^{X}}{\mathrm{MJy}\,{\mathrm{sr}}^{-1}}\right),\end{array}\end{eqnarray} \tag{ 5 }$$

where C₂₄ is an empirically anchored conversion factor that has units of M_⊙ yr⁻¹ (erg s⁻¹)⁻¹ (e.g., following Kennicutt & Evans 2012). The fiducial conversion in Equation (5) is the one suggested for λ = 24 μm by Kennicutt & Evans (2012) and Jarrett et al. (2013) and shown by Leroy et al. (2019) to match the integrated-galaxy population synthesis modeling results of Salim et al. (2016, 2018) well.

In a detailed analysis combining PHANGS–MUSE and WISE data, Belfiore et al. (2022b) show that C_WISE4 ∼ C₂₄ varies as a function of the local conditions in a galaxy disk, likely reflecting heating due to sources other than star formation contributing to the mid-IR. They consider mid-IR in combination with Hα and UV emission, and we note that their maximum C for UV emission is ${\mathrm{log}}_{10}{C}_{\mathrm{WISE}4}\,\sim \,$−42.7, similar to Equation (5), but that for regions with significant heating by old stellar populations or when combining mid-IR and Hα, they find even lower ${\mathrm{log}}_{10}{C}_{\mathrm{WISE}4}\,\sim \,$−42.9. We will return to compare to their findings in Section 6.

As above, the factor ${{ \mathcal R }}_{24}^{X}$ in Equation (5) translates from the fiducial 24 μm to other bands assuming the fixed band ratios in Table 1. However, one should bear in mind that these ratios have been calculated for extended regions of moderate intensity emission (Appendix A) and that the strongest band ratio variations observed in other papers in this series are found at high intensity in star-forming regions (see Section 2.1).

The conversions in above can be equated to express predictions for the relationships between CO and mid-IR emission and Hα and mid-IR emission. First, when gas is mostly molecular and the dust is illuminated by a relatively weak radiation field, Equations (2) and (4) together yield

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{I}_{\nu }^{X}}{\mathrm{MJy}\,{\mathrm{sr}}^{-1}} & \approx & \,0.134\,{{ \mathcal R }}_{24}^{X}\,\left(\displaystyle \frac{{I}_{\mathrm{CO}\,2-1}}{{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)\\ & & \times \left(\displaystyle \frac{{\alpha }_{\mathrm{CO}\ 1-0}}{{\alpha }_{\mathrm{CO}\ 1-0}^{\mathrm{MW}}}\right)\left(\displaystyle \frac{{R}_{21}}{0.65}\right)\left(\displaystyle \frac{D/G}{0.01}\right)\left(\displaystyle \frac{U}{{U}_{0}}\right),\end{array}\end{eqnarray} \tag{ 6 }$$

so that for this case, which might be thought of as "molecular IR cirrus," the mid-IR tracks CO emission with additional dependence on the CO-to-H₂ conversion factor, CO 2–1 to 1–0 line ratio, dust to gas mass ratio, and ISRF. For a PAH-dominated band, the equation should include an additional factor $\times \left({q}_{\mathrm{PAH}}/0.046\right)$.

Meanwhile, if we consider mid-IR driven only by heating due to star formation and require consistency between extinction-corrected Hα and mid-IR results, then Equations (3) and (5) together imply

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{I}_{\nu }^{X}}{\mathrm{MJy}\,{\mathrm{sr}}^{-1}} & \approx & 2.16\times {10}^{5}\,{{ \mathcal R }}_{24}^{X}\\ & & \times \left(\displaystyle \frac{{I}_{{\rm{H}}\alpha }}{\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{sr}}^{-1}}\right)\left(\displaystyle \frac{{10}^{-42.7}}{{C}_{24}}\right)\end{array}\end{eqnarray} \tag{ 7 }$$

for the case of "mid-IR produced directly by reprocessed light from young stars." We caution that while Equation (7) appears simple, depending only on C₂₄ and ${{ \mathcal R }}_{24}^{X}$, this largely reflects that Equation (5) combines all of the physics behind the empirically calibrated mid-IR to SFR conversion into a single number. Moreover, because C₂₄ for the mid-IR is empirical and often calibrated against Balmer decrement-corrected Hα or other recombination line emission, Equation (7) is somewhat circular. We still find this a useful way to frame the current field. For a C calibrated based on star-forming regions or overwhelmingly star-forming galaxies, we might expect Equation (7) to hold for bright star-forming regions within our target galaxies.

In our analysis we will be interested in "bright" mid-IR emission, by which we mean emission that we might plausibly expect to emerge from regions dominated by molecular gas or recent star formation. The equations in Section 2.5 allow us to roughly define two plausible thresholds.

Intensity for molecular-gas-dominated lines of sight: For metal-rich star-forming galaxies, molecular gas typically makes up most of the ISM in regions with Σ_gas ≳ 15 M_⊙ pc⁻² (e.g., Bigiel et al. 2008; Leroy et al. 2008; Schruba et al. 2011). Assuming U ∼ (1 − 2)U₀ and D/G ∼ 0.01 then this implies

$$\begin{eqnarray}&&{I}_{\nu }^{X}\gtrsim 0.3-0.6\,{{ \mathcal R }}_{24}^{X}\,\mathrm{MJy}\,{\mathrm{sr}}^{-1}\end{eqnarray} \tag{ 8 }$$

for "molecular cirrus" with the range affecting a plausible range of U and D/G, each of which affects the estimate linearly (Equation (4)).

Intensity for H ii regions: Meanwhile, though there is no single cutoff between H ii region emission and DIG, the 16th percentile for H ii region surface brightness in Belfiore et al. (2022a) is ${\mathrm{log}}_{10}{I}_{{\rm{H}}\alpha }\sim 38.9$ for I_Hα in units of erg s⁻¹ kpc⁻², which translates to ${\mathrm{log}}_{10}{I}_{{\rm{H}}\alpha }\sim -5.2$ in our units of erg s⁻¹ cm⁻² sr⁻¹. Following Equation (7), this implies

$$\begin{eqnarray}&&{I}_{\nu }^{X}\gtrsim 1.4\,{{ \mathcal R }}_{24}^{X}\,\mathrm{MJy}\,{\mathrm{sr}}^{-1}\end{eqnarray} \tag{ 9 }$$

for H ii regions.

We will return to an empirical version of the two thresholds below.

The mid-IR data were obtained using the MIRI instrument with the F770W, F1000W, F1130W, and F2100W filters. Details of the observations and data reduction appear in Lee et al. (2023). We used the PHANGS–JWST internal release "version 0.5," which uses pipeline version 1.7.0 and CRDS context 0968 and follows the procedure described in Appendix B to set the background level in the maps to be self-consistent among the four MIRI bands and to match previous wide-field observations at 8 μm by Spitzer or 12 μm by WISE.

We compare the mid-IR images to ALMA CO (2–1) maps obtained as part of the PHANGS–ALMA survey (Leroy et al. 2021). We use the combined 12-m+7-m+total power data cubes from the public data release ("v4"). Taking into account typical CO/Hα and CO-to-MIR ratios, the CO (2–1) data are significantly less sensitive than the other data in this work (Section 3.3 and Table 3). Therefore, we construct a special set of "flat" CO-integrated intensity maps, designed to allow simple, robust statistical averaging. To do this, after convolving the CO cube to our working resolution of 17, we shift each spectrum of the cube along its velocity axis, recentering the spectrum for each line of sight so that v = 0 km s⁻¹ now corresponds to the expected mean local rotation velocity. For NGC 0628, 1365, and 7496 we use the low-resolution velocity field derived from the CO as a reference, filling in with a predicted local velocity from the rotation curve in the few regions without detection. For IC 5332, which has lower signal to noise than the other galaxies in CO, we use an estimated rotation curve for the reference at all locations. After adjusting the cube so that all emission is centered in roughly the same channel, we integrate over a fixed velocity window picked to encompass all readily detected emission in the disk (δ v = 25, 35, 80, and 55 km s⁻¹ for IC 5332, NGC 0628, 1365, and 7496). The integration also includes all bright (signal-to-noise ratio, S/N > 3 in 2 channels) emission in extended line wings, which effectively captures the broad wings in the centers of NGC 1365 and NGC 7496.

Table 3. Characteristic Noise at 17 Expressed in Different Units

Band	Expected σ/σ7.7	σ1.7	${\sigma }_{1.7}^{7.7}$	${\sigma }_{1.7}^{{\rm{H}}\alpha }$	${\sigma }_{1.7}^{\mathrm{CO}\ 2-1}$	${\sigma }_{1.7}^{\mathrm{gas}}$
		(MJy sr−1)	(MJy sr−1 at F770W)	(erg s−1 cm−2 sr−1)	(K km s−1)	(M⊙ pc−2)
F770W	1.0	0.025	0.025	0.63 × 10−7	0.027	0.18
F1000W	1.1 ± 0.1	0.028	0.072	1.7 × 10−7	0.076	0.51
F1130W	1.35 ± 0.1	0.033	0.025	0.59 × 10−7	0.027	0.18
F2100W	2.5 ± 0.2	0.063	0.10	2.9 × 10−7	0.031	0.21
CO (2–1) a	⋯	⋯	⋯	⋯	1.0	6.7
Hα b	⋯	⋯	⋯	0.035 × 10−7	⋯	⋯

Notes. Approximate characteristic noises for our data before any inclination correction for our sample at our working resolution of θ = 17. Columns: Expected σ/σ^7.7—approximate ratio of pipeline noise estimate at this band to that at F770W at their native resolutions; σ_1.7—noise on the scale of intensity in that band; ${\sigma }_{1.7}^{7.7}$—noise scaled to units of 7.7 μm via ${{ \mathcal R }}^{X}$ for comparison; ${\sigma }_{1.7}^{{\rm{H}}\alpha }$—noise in equivalent extinction-corrected Hα units assuming the data set wide median ratio; ${\sigma }_{1.7}^{{CO}2-1}$—noise in CO (2–1) units assuming the data set wide median ratio; and ${\sigma }_{1.7}^{\mathrm{gas}}$—noise in gas surface density units assuming data-set-wide median CO-to-mid-IR ratio and a typical CO-to-H₂ conversion factor (Equation (2)).

^a Median noise across all four "flat" integrated intensity (i.e., moment 0) maps (Section 3.1) used for this analysis. Varies moderately from target to target. ^b Typical noise in the MUSE NGC 1365 and NGC 7496 Hα maps at θ ∼ 150 pc. The extinction correction introduces additional uncertainty but does not affect detection.

Download table as:  ASCIITypeset image

The resulting "flat" moment maps appear in the bottom-left corners of Figures 1, 2, 3, and 4. These maps include almost all of the CO flux in each data cube, similar to the standard PHANGS–ALMA high-completeness "broad" moment maps. However, those broad maps use a local velocity window based on multiresolution signal detection algorithms and do include empty lines of sight where no signal is detected. The "flat" maps that we use here have a more even noise distribution because they use a single fixed velocity window and have no blank pixels within the region of interest. As a result, they are moderately noisier in any given pixel than the fiducial PHANGS–ALMA products but can be averaged spatially in a simple way to produce results equivalent to spectral stacking (e.g., Schruba et al. 2011; Ianjamasimanana et al. 2012). Aside from the higher but more even noise, which can be seen particularly clearly in the images of NGC 1365 (Figure 2) and IC 5332 (Figure 4), these maps capture essentially the same features as the standard PHANGS broad maps (compare to the atlas images for the same galaxies in Leroy et al. 2021).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We compare both mid-IR and CO emission to extinction-corrected Hα emission from the PHANGS–MUSE survey (Emsellem et al. 2022). In this work, we primarily use maps of Hα emission corrected for internal extinction using the Balmer decrement method contrasting Hβ and Hα. These maps were produced from the "convolved and optimized" maps PHANGS–MUSE internal data release 2.2 (see Emsellem et al. 2022) and have been described in Belfiore et al. (2022a, 2022b) and Pessa et al. (2021, 2022).

We use kernels produced following the method of Aniano et al. (2011) and the preflight estimates of the JWST PSFs to convolve all of our data to share a common Gaussian PSF with FWHM θ = 17. Table 2 gives the corresponding physical beam size for each target, which ranges from 70 to 160 pc. After convolution, we reproject all data onto a common astrometric grid centered at the galaxy center with 083-wide pixels, i.e., 4 pixels per PSF area.

Both NGC 1365 and NGC 7496 have a bright active galactic nucleus (AGN) in the inner galaxy, which leads to diffraction spikes that extend well into the surrounding galaxy disks. Following Hassani et al. (2023) we use an image of the PSF centered at the bright peak to identify the diffraction spikes. We convolve this mask to the working resolution of our data and drop pixels where diffraction spikes are expected to cover >1% of the area from our analysis.

After this processing, we build a database consisting of pixel-by-pixel measurements of the intensities of CO (2–1), extinction-corrected Hα, mid-IR emission at 7.7 μm (F770W), 10 μm (F1000W), 11.3 μm (F1130W), and 21 μm (F2100W). We correct all intensities by a factor of cos i to the values expected had we observed the galaxies face on. We draw inclinations from Lang et al. (2020) and Leroy et al. (2021).

We mostly utilize intensities in our expectations and analysis, it can also be of interest to compute the flux, luminosity, mass, or SFR of an individual resolution element. For our common angular resolution of θ = 17, to convert from intensity to flux per beam, the beam area is ≈7.4 × 10⁻¹¹ sr. To allow easy conversion from inclination-corrected surface density to mass or SFR within a resolution element, Table 2 also quotes the physical beam size ${A}_{\mathrm{beam}}=\pi {\left(\theta d\right)}^{2}\,{\left(4\mathrm{ln}2\right)}^{-1}\,{\left(\cos i\right)}^{-1}$ for each target.

Table 3 provides characteristic noise levels for our θ = 17 images. We quote a single value for each band estimated before any inclination correction is applied. We give single, approximate values per band because we work with early versions of the data and our images lack significant "empty sky." Given that we have convolved the data to significantly coarser-than-native resolution, we view an empirical estimate as more appropriate than an aggressive extrapolation of the native resolution pipeline noise estimate.

To derive these estimates, we use a version of the procedure described in Appendix B. We consider the low-intensity region of each galaxy. As shown in Sandstrom (2023b), even these faint regions still often have significant real emission from the galaxy, which could masquerade as noise. To account for this, we scale another band by the typical band ratio and subtract that scaled image from the original then measure the noise in the difference, e.g., considering ${I}_{\nu }^{{\rm{F}}770{\rm{W}}}-{{ \mathcal R }}_{11.3}^{7.7}{I}_{\nu }^{{\rm{F}}1130{\rm{w}}}$. After running a broad median filter across the difference image, we solve for the implied rms noise, taking into account both the measured band ratio and the noise ratio expected based on the pipeline noise estimates. That is, we measure the noise in the difference image and use this, with an appropriate scaling, to estimate the noise in the original images.

In order to compare the sensitivity of different bands to gas column or recent star formation, Table 3 also re-expresses the noise in the bands scaled in several different ways. We first express the relative sensitivity of all bands to "typical" dust emission by scaling by the band ratios in Table 1 onto a common F770W intensity scale, ${\sigma }_{1.7}^{7.7}$. Then we adopt the median CO-to-mid-IR and Hα-to-mid-IR ratios found in Section 4 to express the mid-IR noises in equivalent I_Hα , I_CO, and Σ_gas units. In Section 6, we will use these to comment on the sensitivity of the JWST images to recent star formation and gas column density.

In addition to the noise levels quoted in Table 3, there is an overall ≲0.1 MJy sr⁻¹ uncertainty in the level of the background in the MIRI images. Realistically, based on a visual inspection of the images at high stretch, this background level clearly varies across the images at a fraction of this level. To assess this quantitatively, we note that the medium-scale median filter level in the difference images mentioned above imply an rms variation ∼ ± 0.04−0.07 MJy sr⁻¹ in the background level. Comparing to Table 3 highlights that at our working resolution the uncertainty in the background is often as important as the statistical noise.

We focus our analysis on regions with an inclination-corrected F770W intensity above 0.5 MJy sr⁻¹. We target this level because based on the calculations in Section 2.6, we expect that where molecular gas makes up most of the ISM along a line of sight it will likely produce at least this level of emission, even if it is only weakly illuminated. The threshold for emission associated with H ii regions is even higher, so we expect this selection to capture most of the area in our targets where either dust mixed with molecular gas or dust heated by H ii regions contributes. This level is also high enough that uncertainty in the background level and the statistical noise, which both have a maximum magnitude ∼±0.1 MJy sr⁻¹, are only secondary concerns for the mid-IR.

Table 4 reports the fractions of flux in each band and area within the PHANGS–JWST footprint associated with this "bright" emission. In our three more massive targets, this selection captures ≳90% of mid-IR band, CO (2–1), and Hα flux, though only a much lower ∼50% fraction of the area. This partially reflects that PHANGS–JWST specifically targets the regions of active star formation in our target galaxies. In less molecular-gas-dominated, less actively star-forming regions, lower intensity emission from dust mixed with atomic or CO-dark gas will contribute larger fractions of the flux (Section 2 and see Sandstrom 2023a). In our sample, the lower surface density IC 5332 exemplifies this case, and our selection captures only about half of the total flux and about 20% of the area (see Figure 4).

Zoom In Zoom Out Reset image size

Table 4. Correlation, Ratio, and Fitting Results

Quantity	Variable(s)	All Data	IC 5332	NGC 0628	NGC 1365	NGC 7496
Fraction Associated with "Bright" Emission
(Entries Report Fraction of Flux at the Indicated Band or Area Selected for Our Analysis)
Flux	F770W	0.94	0.44	0.97	0.95 a	0.90
Flux	F1000W	0.93	0.46	0.97	0.94 a	0.90
Flux	F1130W	0.93	0.46	0.97	0.94 a	0.91
Flux	F2100W	0.96	0.49	0.98	0.98 a	0.95
Flux	CO (2–1)	0.97	0.68	0.97	0.98 a	0.96
Flux	Hα	0.96	0.55	0.98	0.99 a	0.96
Area	⋯	0.60	0.17	0.87	0.60 a	0.54
Rank Correlations
(Entries Report Rank Correlation Coefficient Relating Indicated Variable Pair)
Rank correlation	CO (2–1) versus Hα	0.47	0.12	0.47	0.41	0.57
Rank correlation	CO (2–1) versus F770W	0.63	0.20	0.70	0.52	0.72
Rank correlation	CO (2–1) versus F1000W	0.62	0.19	0.70	0.52	0.71
Rank correlation	CO (2–1) versus F1130W	0.63	0.21	0.71	0.53	0.72
Rank correlation	CO (2–1) versus F2100W	0.61	0.17	0.71	0.50	0.66
Rank correlation	Hα versus F770W	0.78	0.63	0.74	0.80	0.85
Rank correlation	Hα versus F1000W	0.74	0.65	0.67	0.75	0.83
Rank correlation	Hα versus F1130W	0.74	0.54	0.71	0.76	0.81
Rank correlation	Hα versus F2100W	0.75	0.63	0.72	0.73	0.78
Median Ratios for Bright Emission
(Entries Report Median ${\mathrm{log}}_{10}y/x$ and Scatter for Indicated Ratio)
${\mathrm{log}}_{10}y/x$	CO (2–1)/Hα	5.63 ± 0.59	5.42 ± 0.66	5.63 ± 0.55	5.76 ± 0.64	5.46 ± 0.60
${\mathrm{log}}_{10}y/x$	CO (2–1)/F770W	0.04 ± 0.33	0.10 ± 0.36	0.03 ± 0.29	0.14 ± 0.40	−0.11 ± 0.30
${\mathrm{log}}_{10}y/x$	CO (2–1)/F1000W	0.44 ± 0.32	0.51 ± 0.37	0.44 ± 0.28	0.53 ± 0.38	0.28 ± 0.29
${\mathrm{log}}_{10}y/x$	CO (2–1)/F1130W	−0.10 ± 0.32	−0.03 ± 0.36	−0.11 ± 0.28	−0.03 ± 0.39	−0.25 ± 0.29
${\mathrm{log}}_{10}y/x$	CO (2–1)/F2100W	0.30 ± 0.34	0.31 ± 0.41	0.32 ± 0.29	0.37 ± 0.43	0.09 ± 0.33
${\mathrm{log}}_{10}y/x$	Hα/F770W	−5.60 ± 0.43	−5.39 ± 0.48	−5.62 ± 0.42	−5.62 ± 0.41	−5.57 ± 0.47
${\mathrm{log}}_{10}y/x$	Hα/F1000W	−5.21 ± 0.47	−4.97 ± 0.52	−5.22 ± 0.47	−5.23 ± 0.45	−5.20 ± 0.51
${\mathrm{log}}_{10}y/x$	Hα/F1130W	−5.76 ± 0.45	−5.51 ± 0.50	−5.77 ± 0.44	−5.79 ± 0.44	−5.73 ± 0.51
${\mathrm{log}}_{10}y/x$	Hα/F2100W	−5.34 ± 0.44	−5.13 ± 0.46	−5.33 ± 0.42	−5.38 ± 0.44	−5.36 ± 0.50
Power-law Fits to Binned Data
(Entries Report Slope m and Intercept b in Equations (10) and (11))
m, b	CO (2–1) versus Hα b	0.47, 0.30	0.56, − 0.29	0.39, 0.40	0.50, 0.29	0.51, 0.12
m, b	CO (2–1) versus F770W	1.10, −0.07	0.80, −0.08	0.84, 0.05	1.25, −0.10	1.08, −0.17
m, b	CO (2–1) versus F1000W	1.22, 0.37	0.71, 0.19	0.92, 0.39	1.33, 0.39	1.16, 0.27
m, b	CO (2–1) versus F1130W	1.21, −0.28	0.97, −0.21	0.93, −0.11	1.33, −0.34	1.16, −0.36
m, b	CO (2–1) versus F2100W	0.74, 0.24	0.22, − 0.08	0.61, 0.32	0.95, 0.20	0.79, 0.05
m, b	Hα versus F770W	1.51, –5.72	1.98, −5.39	1.67, −5.80	1.43, −5.68	1.53, −5.68
m, b	Hα versus F1000W	1.51, −5.14	1.96, −4.64	1.72, −5.13	1.44, −5.19	1.66, −5.07
m, b	Hα versus F1130W	1.55, −5.97	2.20, −5.70	1.71, −6.03	1.47, −5.98	1.71, −6.00
m, b	Hα versus F2100W	1.29, −5.35	1.66, −5.04	1.35, −5.34	1.24, −5.42	1.29, −5.34

Notes. See Section 5. Hα refers to extinction-corrected Hα in units of erg s⁻¹ cm⁻² sr⁻¹. CO (2–1) intensity has units of K km s⁻¹. All mid-IR intensities have units of MJy sr⁻¹.

^a Fraction of flux and area above the lower threshold. In NGC 1365, we also exclude bright emission from the center of the galaxy from our analysis. ^b Note the definition for the intercept location for CO (2–1) versus Hα in Equation (11).

Download table as:  ASCIITypeset image

We also exclude the bright center of NGC 1365, defined by a cut of ${I}_{\nu }^{{\rm{F}}770{\rm{W}}}\gt 30$ MJy sr⁻¹ from most of our statistical analysis. This region hosts an AGN, forms stars at a rate of order ∼10 M_⊙ yr⁻¹, and represents a distinct regime from our other targets in terms of gas surface density, star formation intensity, and dust properties. We focus this paper on normal galaxy disks and leave the contrast between disk and starburst environments for future work, but note that several other papers in this Issue focus specifically on this rich region in NGC 1365 (Liu et al. 2023; Schinnerer et al. 2023; Whitmore et al. 2023). For the correlations, we do plot the data from this higher intensity region but indicate the region excluded from statistical analysis by shading.

When analyzing CO (2–1) as a function of extinction-corrected Hα, we focus on an intensity range defined by scaling our bright emission definitions by the typical Hα^corr to mid-IR ratio (Table 4) and then adopting the same limits used for the mid-IR.

In our analysis, all of the mid-IR bands correlate well with CO (2–1) at 70–160 pc resolution. Given the conventional view of the mid-IR as a star formation tracer, we emphasize the impressive agreement between the CO and mid-IR maps at high physical resolution as one of our key results. Despite the comparatively high noise level in the CO, CO and mid-IR typically exhibit correlation coefficients of ρ ≈ 0.62 and show a clear correlation across a wide range of mid-IR intensities. CO (2–1) as a function of mid-IR emission shows an approximately linear slope, i.e., near m ∼ 1, with slope m_CO−MIR in the range m_CO−MIR ∼ 0.8−1.3 for most cases, and all four galaxies are in fairly good agreement with one another (Figures 5, 6, and 10). Treating all galaxies together, we find a mildly superlinear slope, m_CO−MIR ∼ 1.1–1.2 for F770W, F1000W, and F1130W and a mildly sublinear slope for F2100W. In the simplest terms, this slope near m_CO−MIR ∼ 1 implies that these bands not only correlate with CO, but that the CO and mid-IR actually track one another with a nearly fixed constant of proportionality.

Visually, this agrees with the impression of excellent agreement when one "blinks" the CO and mid-IR maps (e.g., Figures 1 through 4). The structure of the low-intensity, extended emission in the mid-IR map resembles that seen in the CO map. Observationally, this reflects that the JWST mid-IR observations have reached the resolution and sensitivity where they capture the glow of dust in neutral ISM heated by the ISRF. Physically, the nearly linear slope could reflect that variations in the dust-to-gas ratio, radiation field, and PAH abundance are weak compared to overall column density variations, so that over a large part of the region mapped by PHANGS–JWST, the bright mid-IR simply traces the molecular gas. The slightly superlinear slope in the relation for F770W, F1000W, and F1130W might indicate that either the dust-to-gas ratio, the PAH abundance, or the radiation field have a modest correlation with column density (all of these are reasonable to expect; e.g., Jenkins 2009; Sandstrom et al. 2010; Roman-Duval et al. 2017; Chastenet et al. 2019). The slightly sublinear slope for F2100W may reflect that bright star-forming regions influence the observed CO versus mid-IR correlation more heavily for this band. We attempt to disentangle these regions from the cold ISM in Section 5.

Mid-IR emission should depend on the ISRF. We do note that CO (2–1) brightness may also track the ISRF, U, because the ISRF should affect molecular gas heating and the CO (2–1) brightness will reflect the underlying excitation and temperature of the molecular gas (e.g., Bolatto et al. 2013; Gong et al. 2020; Liu et al. 2021). This could enhance the correlation that we observe, and these temperature and excitation effects can render the interpretation of both CO and mid-IR as a simple gas tracer more challenging.

In our analysis, extinction-corrected Hα and mid-IR show an even stronger correlation than CO and mid-IR emission, with ρ ≈ 0.75. Though the difference in ρ can probably be ascribed to the much better S/N in the Hα compared to the CO (Section 3.3), the fact remains that Hα and mid-IR correlate very well in our data. Given that both mid-IR emission and Hα are widely viewed as tracers of massive, recent star formation, this correlation is expected, and the power-law fits in Table 4 can be applied to predict the extinction-corrected Hα from the mid-IR intensity with reasonable accuracy.

Figure 10 and Table 4 show that Hα as a function of mid-IR emission exhibits a consistently steeper-than-unity slope. For F770W, F1000W, and F1130 we find m_Hα−MIR ∼ 1.4 − 1.7 with m_Hα−MIR ∼ 1.5 on average. F2100W also shows a superlinear slope, but a shallower one, m_Hα−MIR ∼ 1.3 on average (see Section 4.4 for a discussion of the differences).

Zoom In Zoom Out Reset image size

These slopes with m_Hα−MIR > 1 have the sense that extinction-corrected Hα becomes brighter relative to mid-IR at high intensities and such a trend appears consistent between all four systems (Figure 8). Visually, this agrees with the impression from Figures 1–4 that there is more extended and low-intensity mid-IR emission compared to Hα emission. That is, the peaks of the Hα and mid-IR align well, but the mid-IR shows a much more significant "diffuse" or extended component. The existence of such a component is in good agreement with previous infrared observations of very nearby galaxies (e.g., Walterbos & Schwering 1987; Calzetti et al. 2005; Leroy et al. 2012; Li et al. 2013; Groves et al. 2012; Crocker et al. 2013; Calapa et al. 2014; Boquien et al. 2015, 2016; Kim et al. 2021; Kim 2023; Belfiore et al. 2022b; and clearly visible in Figures 1 through 4). Because the bright Hα emission tends to be concentrated in the H ii regions, this extended IR emission is associated with low Hα-to-mid-IR ratios while the most luminous regions show the highest Hα-to-mid-IR ratios.

Zoom In Zoom Out Reset image size

A number of recent studies using PHANGS–MUSE have shown that Balmer decrement extinctions show a positive correlation between Hα intensity and Hα extinction, A_Hα (see Belfiore et al. 2022a; Emsellem et al. 2022; Santoro et al. 2022). The same trend of high extinction associated with high Hα emission or SFR appears to hold for integrated star-forming galaxies (e.g., Garn & Best 2010; Ly et al. 2012; Dominguez et al. 2013). As a result, in this data set, the regions with high I_Hα that show the lowest mid-IR emission relative to Hα also show the highest Hα extinction, on average.

This result has the same sense as the trend found by Belfiore et al. (2022b) at coarser resolution for the whole PHANGS–MUSE sample, that mid-IR emission becomes brighter relative to other tracers of star formation at low intensities. They expressed this as a dependence of the mid-IR-to-SFR coefficient C (Section 2.4, Equation (5)) on specific SFR or star formation surface density and found that more intensely star-forming parts of galaxies required a smaller C (i.e., a more negative exponent) than more quiescent regions. Here we see that this result continues to hold line of sight by line of sight on scales of 70–160 pc.

Note that this behavior differs somewhat from results focusing directly on star-forming complexes or actively star-forming galaxy centers by Calzetti et al. (2007). In their work, the ratio of mid-IR-to-extinction-corrected Pα increases or remains approximately constant with increasing star formation activity (they find the equivalent of m_Pα−MIR ≈ 1.06 for 8 μm while m_Pα−MIR ≈ 0.81 for 24 μm). Any of the "hybrid" tracers that linearly use the mid-IR with an obscured term also predict a ratio of mid-IR to extinction-corrected Hα that increases with extinction. The difference seems naturally explained by our inclusion of all mid-IR emission in our analysis, including emission not directly associated with star-forming regions. In Section 5, we return to this question, attempting to describe the full mid-IR emission as a combination of "cirrus" emission associated with gas and emission tracing star-forming regions.

Thus the explanation for our steep Hα versus mid-IR slope appears to mix at least two effects: (1) the impact of "cirrus" emission by dust not immediately associated with the star-forming regions, which preferentially contributes mid-IR to low-intensity regions and (2) the destruction of PAHs and small dust grains (for the F770W, F1000W, and F1130W bands), which we discuss more in Section 4.4 and likely explains the difference between the F2100W and the other bands.

The overall nonlinear slope of the relationship between extinction-corrected Hα and the mid-IR intensity also shows that at these resolutions, the mid-IR cannot be translated into a local estimate of star formation activity at high precision by a single factor. Even if one ignores that much of the Hα emission arises from nonlocal ionizations, i.e., the DIG, the slope of ∼1.3–1.5 relating Hα-to-mid-IR emission still leads to a scatter of 0.4–0.5 dex, implying an FWHM of a full order of magnitude in the ratio of Hα-to-mid-IR across our data set. Our power-law fits or the two-component model in Section 5 offer a better approach, and we expect to produce even sharper prescriptions using the full PHANGS–JWST survey.

To help interpret the CO versus mid-IR and Hα versus mid-IR correlation, we apply an identical analysis of the correlation of CO (2–1) and extinction-corrected Hα and show the results in Figure 9 and Table 4.

Zoom In Zoom Out Reset image size

The relationship between CO and Hα is weaker, more scattered, and less linear than that between mid-IR and either CO or Hα. Over the whole data set, CO and extinction-corrected Hα show a rank correlation coefficient ρ ≈ 0.47, weaker than either the correlation between CO and mid-IR emission or Hα and mid-IR emission. The scatter in the CO-to-Hα ratio is larger than that in either the CO-to-mid-IR ratio or Hα-to-mid-IR ratio, and the scatter of CO intensity within individual Hα intensity bins is larger than that of CO in bins of mid-IR emission (compare Figures 5 and 9). Finally, the best-fitting slope of the CO versus Hα relation is very shallow, m_CO−Hα ≈ 0.5.

The shallow slope, high scatter, and weaker correlation coefficient all indicate that the global scaling between molecular gas and star formation activity traced by Hα is beginning to break down at the 70–160 pc scales of our data. For reference, at 1.5 kpc resolution and comparing extinction-corrected Hα and molecular gas estimated from CO in the full PHANGS–MUSE sample, Sun et al. (2022a) and Sun et al. (2022b) find a slope of m ∼ 1.1, scatter about the best-fit power law of σ ≈ 0.3 dex, and a rank correlation coefficient of ≈0.88 (see also Pessa et al. 2021, 2022). This agrees with previous work showing a separation of CO and Hα emission in this sample at similar scales (Kreckel et al. 2018; Schinnerer et al. 2019; Chevance et al. 2020; Kim et al. 2022; Pan et al. 2022). It also agrees with the wider literature demonstrating a spatial separation between tracers of recent massive star formation and molecular gas at high resolution (including Kawamura et al. 2009; Schruba et al. 2010; Leroy et al. 2013; Corbelli et al. 2017; Grasha et al. 2018; Kruijssen et al. 2019; Turner et al. 2022). ⁵¹

The results that the CO versus mid-IR and Hα versus mid-IR relationships each appear significantly stronger than the CO versus Hα relationship and that the CO versus Hα relationship has begun to break down in our data both support the idea that our observations can distinguish the impacts of column density and heating on mid-IR emission. Lending further support to this view, the residuals about the CO versus mid-IR and Hα versus mid-IR show clear secondary correlations. In Figures 5 and 7, we color the data points by the mean intensity of the unplotted variable (i.e., Hα in the CO–IR plot, CO in the Hα–IR plot). Both figures reveal a secondary dependence on this unplotted variable. That is, for a given CO intensity, the mid-IR still appears to correlate with Hα and vice versa.

Taken together, the results of our analysis suggest that we are observing distinct relationships between CO versus mid-IR and Hα versus mid-IR. These are probably amplified by the still-present, if weak, correlation between CO and Hα in our data. Beyond these first results, analysis of the partial correlation coefficients or principal component analysis represent logical ways forward. In this paper we adopt a simple two-component modeling approach as the next step (Section 5).

To first order, we find consistent results among all four mid-IR bands and Appendix A shows that a single scaling relating the bands represents a reasonable description at modest intensity and modest resolution. In more detail, Table 4 and Figures 5 through 10 suggest moderately different behavior for the F2100W 21 μm band compared to the other three PHANGS–JWST mid-IR bands. Hα versus F2100W shows a shallower, nearly linear slope of m_Hα−MIR ∼ 1.3 compared to the other bands, which show m_Hα−MIR ∼ 1.5. Meanwhile CO (2–1) versus F2100W also shows a shallower slope, m_CO−MIR ∼ 0.7, compared to m_CO−MIR ∼ 1.2 for the other bands. These differences hold within each galaxy, as well as across the whole data set.

Zoom In Zoom Out Reset image size

The differences between the Hα versus F2100W and Hα versus F770W that we observe agree well with the results of Calzetti et al. (2007) studying 8 μm and 24 μm in star-forming regions in SINGS galaxies. They used extinction-corrected Paschen α and found that m_Pα−MIR ≈ 1.06 for 8 μm while m_Pα−MIR ≈ 0.81 for 24 μm. ⁵² The difference in slopes between the 8 μm and 24 μm is almost identical to what we observe here between the F2100W and F770W. However, the values of the slopes themselves differ due to differences in experiment design and sample, with our analysis yielding significantly steeper slopes. As discussed above, the difference should be expected, since Calzetti et al. (2007) focuses on star-forming regions and galaxy centers, while we plot all lines of sight for our target galaxies.

The simplest explanation of the differences in slopes between Hα and the mid-IR for different bands is that F2100W appears to trace the extinction-corrected Hα, itself a tracer of recent star formation and heating, more directly than the other mid-IR bands trace Hα. This is consistent with the observation that the emission of the PAH-tracing bands, F770W and F1130W, relative to the dust continuum, diminishes in bright H ii regions. This result is the main topic of Chastenet (2023a), Dale et al. (2023), and Egorov (2023) in this Issue. Previously, many Spitzer and WISE studies of the Milky Way and the nearest galaxies observed similar effects, either as a decline in the 8 μm/24 μm (or WISE 12μm/22μm) ratio in H ii regions or as a visible ringlike morphology for the PAH-tracing band, suggesting that the brightest emission surrounds the H ii region while the 22 μm or 24 μm emission peaks in the star-forming region (e.g., Helou et al. 2004; Povich et al. 2007; Bendo et al. 2008; Gordon et al. 2008; Relaño & Kennicutt 2009; Anderson et al. 2014; Calapa et al. 2014; Chastenet et al. 2019).

The differing slopes of the CO versus mid-IR relation for different bands could arise from the same effect. At high mid-IR intensities, there will be more contribution from bright star-forming regions to F2100W than to the PAH-tracing bands. This might mean relatively lower CO-to-F2100W ratios compared to, e.g., CO-to-F770W ratios in these bright regions, and help explain why the PAH-tracing bands tend to show m_CO−MIR ∼ 1.2 while the F2100W shows m_CO−MIR ∼ 0.7.

This relative suppression of PAH emission compared to ∼24 μm emission in star-forming regions is frequently ascribed to PAH destruction in ionized gas. Given this explanation, it may be surprising that F1000W appears better matched to F770W and F1130W than to F2100W. In theory, the F1000W band should be primarily continuum rather than PAH-dominated, with the main feature of note being the silicate absorption band, which we do not expect to be strong at these column densities (e.g., Draine & Li 2007; Smith et al. 2007). The closer coupling of F1000W to the two PAH-dominated bands rather than F2100W is a somewhat surprising result of the first PHANGS–JWST science. Apparently, either the smaller, hotter, more easily destroyed dust grains responsible for 10 μm compared to 21 μm heating mimic PAHs in their behavior, or alternatively weaker PAH features or the line wings of the adjacent strong bands contribute significantly to the band.

This apparent suppression of PAH emission in bright regions, combined with significant metallicity trends observed in the PAH abundance (e.g., Engelbracht et al. 2005; Calzetti et al. 2007; Gordon et al. 2008; Chastenet et al. 2019; Li 2020), led to a preference for the continuum-dominated 24 μm emission over 8 μm as a star formation tracer during the Spitzer era (e.g., Calzetti et al. 2007; Murphy et al. 2011; Kennicutt & Evans 2012; Catalán-Torrecilla et al. 2015). Despite this, PAH-dominated bands, including the 8 μm emission and the 12 μm emission from WISE still often yielded good quantitative correspondence with independent SFR estimates, especially when used as part of "hybrid" tracers with Hα or UV (e.g., Calzetti et al. 2007; Kennicutt & Hao 2009; Jarrett et al. 2013, though see also Calzetti et al. 2005).

Finally, we note that we are finding F2100W to be more linearly associated with extinction-corrected Hα, not only associated with Hα. The slopes relating Hα to F2100W and CO (2–1) to F2100W lie about equidistant from unity, and F2100W still correlates quite well with CO emission, especially at low intensities. The overall local correlation of F2100W with CO still appears excellent and in Section 5 we will find a substantial CO-associated component to contribute to the emission. Similarly, the PAH-tracing bands still correlate exceptionally well with Hα emission; they simply do so nonlinearly. Indeed a main conclusion of this paper is that all four of the bands act simultaneously as tracers of recent star formation and the neutral ISM.

Table 5 reports the best-fit c and h and the fractional contribution of each component to the best-fitting model for each target, f_c and f_h . Figure 11 shows the mean absolute residual in ${\mathrm{log}}_{10}h$ versus ${\mathrm{log}}_{10}c$ space (i.e., the grid of plausible coefficients) for the full data set. Then, Figures 12 and 13 show the resulting model mid-IR intensities as a function of the observed mid-IR intensities using the best-fitting c and h for each galaxy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5. Template-fitting Results

Predicted Band	Quantity	All Data	IC 5332	NGC 0628	NGC 1365	NGC 7496
F770W	Mean $| {\mathrm{log}}_{10}\tfrac{\mathrm{model}}{\mathrm{data}}|$	0.204	0.224	0.180	0.221	0.196
F770W	c	0.468	0.417	0.525	0.288	0.776
F770W	h/105	1.318	0.741	1.202	1.862	0.871
F770W	fc	0.50	0.46	0.51	0.40	0.62
F770W	fh	0.50	0.54	0.49	0.60	0.38
F1000W	Mean $| {\mathrm{log}}_{10}\tfrac{\mathrm{model}}{\mathrm{data}}|$	0.223	0.245	0.199	0.239	0.213
F1000W	c	0.214	0.162	0.234	0.141	0.339
F1000W	h/105	0.447	0.316	0.389	0.692	0.316
F1000W	fc	0.58	0.45	0.59	0.47	0.67
F1000W	fh	0.42	0.55	0.41	0.53	0.33
F1130W	Mean $| {\mathrm{log}}_{10}\tfrac{\mathrm{model}}{\mathrm{data}}|$	0.218	0.230	0.190	0.240	0.213
F1130W	c	0.724	0.589	0.794	0.490	1.202
F1130W	h/105	1.660	0.912	1.445	2.570	0.977
F1130W	fc	0.56	0.50	0.56	0.46	0.70
F1130W	fh	0.44	0.50	0.44	0.54	0.30
F2100W	Mean $| {\mathrm{log}}_{10}\tfrac{\mathrm{model}}{\mathrm{data}}|$	0.203	0.230	0.166	0.239	0.211
F2100W	c	0.269	0.263	0.282	0.178	0.479
F2100W	h/105	0.676	0.479	0.603	0.977	0.550
F2100W	fc	0.55	0.43	0.55	0.49	0.61
F2100W	fh	0.45	0.57	0.45	0.51	0.39

Notes. Best-fit model coefficients, c and h, to match the observed mid-IR intensity using a linear combination of scaled CO and Hα (Section 5, Equations (12)). c is the coefficient to scale CO and has units of MJy sr⁻¹ (K km s⁻¹)⁻¹. h is the coefficient to scale extinction-corrected Hα and has units of MJy sr⁻¹ (erg s⁻¹ cm⁻² sr⁻¹)⁻¹. Following Equation (13), f_c and f_h report the fraction of the total model flux associated with the CO-tracing component and Hα-tracing component.

Download table as:  ASCIITypeset image

Table 5 and Figures 11 through 13 show that this simple procedure does a fairly good job of fitting the data. We find coefficients that clearly minimize the residuals at ≈ ±0.2 dex. Figures 12 and 13 show that the resulting models do a good job of matching the observations over a wide range of intensities. The color in Figure 12 indicates the relative contribution of the Hα-tracing component and CO-tracing component in that bin. As expected, the Hα-tracing component contributes more emission at the higher mid-IR intensities, while the CO-tracing component makes up most of the model at lower mid-IR intensities.

Figure 11 and Table 5 show that the best-fitting coefficients c and h vary from galaxy to galaxy. Following Section 2, we expect such variations if the ISRF, U; D/G; and mean extinction vary across galaxies. All of these quantities do vary from galaxy to galaxy (e.g., Draine et al. 2007; Aniano et al. 2020), so we suggest to interpret these variations as indicative of real changes in the emissivity and amount of dust in these galaxies. Supporting this view, in Figure 11 and Table 5 all galaxies show similar relative c and h in all four bands. This would be expected if the drivers of the variations are local conditions in the galaxy. For instance, variations in U or D/G should affect c similarly for all four bands.

The main defect in the model, visible in Figures 12 and 13, is that it overpredicts F770W, F1000W, and F1130W emission at high intensity. Figure 13 shows that the overprediction appears present for all galaxies, and Figure 12 shows that it coincides with regions where the Hα-tracing component makes up most of the model. This feature reflects the same relative PAH emission drop in bright H ii regions seen by Chastenet (2023a) and Egorov (2023) and discussed in Section 4.2 and Section 4.4. The F2100W model, on the other hand, appears quite linear across the full intensity range.

Table 5 also reports the flux associated with each component in the model fit. In almost all cases, f_c ≈ f_h , meaning that the model fit assigns about equal flux to the CO component and the Hα component. Taken at face value, this implies that about half of the mid-IR flux from our targets arises directly from star-forming regions traced by Hα, while about half of the flux can be attributed to emission from dust that is mixed with molecular gas and heated by the ISRF.

Extended "diffuse" or "cirrus" mid-IR emission not directly associated with local, recent star formation has been a longstanding topic of interest. Various studies have attempted to remove or model the effects of such emission using local background subtractions (e.g., Calzetti et al. 2005), physical dust models (e.g., Leroy et al. 2012), Fourier filtering (e.g., Kim et al. 2021; Kim 2023), multiwavelength comparisons (Crocker et al. 2013; Calapa et al. 2014; Whitcomb et al. 2022), or scale- (Li et al. 2013; Calzetti 2013; Boquien et al. 2015) or environment-dependent SFR calibrations (Boquien et al. 2016; Belfiore et al. 2022b).

Our experiment uses a physically motivated, template-based approach that leverages JWST's high resolution to estimate that in our target regions, ∼40%–60% of the mid-IR flux at 7.7–21 μm lies in a component that resembles the CO emission. This broadly agrees with the magnitude of previous estimates of the fraction of "diffuse" mid-IR emission ⁵⁴ (e.g., Leroy et al. 2012; Crocker et al. 2013; Kim et al. 2021; Kim 2023) and of the scale- or environment-dependent changes in the prefactor C (Equation (5)) to estimate the SFR (e.g., Boquien et al. 2015; Belfiore et al. 2022b). In addition to prototyping a new approach that yields an independent estimate, our data provide evidence that the "diffuse" mid-IR emission in the star-forming parts of galaxies has a distribution like that of the molecular gas. This could be expected theoretically, because dust mixed with molecular gas will dominate the overall dust budget in these molecule-rich inner parts of galaxies (see Leroy et al. 2012). Here we demonstrate using observations that a large fraction of the dust emission, including 40%–60% of the 21 μm emission often used to trace star formation, appears to resemble the CO.

Conversely, the strong, nearly linear correlation between mid-IR and CO emission at a coarser resolution has led to suggestions that mid-IR can trace molecular gas (e.g., Gao et al. 2019; Chown et al. 2021; Leroy et al. 2021; Leroy 2023; Zhang & Ho 2022). Our analysis supports such a potential use even at a quite high physical resolution. However, the modeling exercise here suggests that at high resolution, emission more directly associated with recent star formation will contribute ≈50% of the mid-IR flux. Because the distributions of recent star formation and neutral gas at these scales no longer perfectly match, this emission will represent an important contaminant that may need to be accounted for or masked out in quantitative gas estimates.

Our fits offer a good first-order description of the bright mid-IR, but several directions for refinement are also immediately obvious. First, to model the bright F770W, F1000W, and F1130W emission, the Hα component should include some nonlinearity to reflect decreased PAH or small-grain emission in the bright H ii regions (see references and discussion in Section 4.4).

At the faint-intensity end, a general form of the model must account for mid-IR emission associated with phases of the ISM other than molecular gas, including diffuse atomic gas or CO-dark molecular gas (see Sandstrom 2023b). Though only IC 5332 has significant flux associated with such regions in our current sample, atomic-gas-dominated regions cover most of the area in most galaxy disks. A fixed offset in the model may be a reasonable first approximation given the rough constancy of the atomic gas layer in inner galaxy disks. Alternatively, this approach could be used to isolate non-CO-emitting gas if all the other components were well understood.

We also note that our analysis treats the DIG and H ii regions together, positing that both types of emission may track dust reprocessing UV light from recent star formation. This may be reasonable to first order given that the DIG in our targets shows a morphology consistent with being created by light leaking from star-forming regions (Belfiore et al. 2022a). In the future, distinguishing the two components when comparing to the mid-IR should yield more insight into how UV light from young stars is reprocessed into dust emission.

Finally, we expect that working at even higher resolution will improve the physical insight gained from the modeling. For these first results, we used a fixed 17, and our physical resolution for NGC 1365 and NGC 7496 remains fairly coarse at ∼150 pc. At this resolution, there remains significant overlap between CO and Hα emission (e.g., see Schinnerer et al. 2019; Kim et al. 2022; Pan et al. 2022), while at higher resolution we expect to be able to distinguish the two components even more sharply. With the full combined PHANGS–JWST, MUSE, and ALMA sample, we will include a larger sample of nearby galaxies and include galaxies with higher -resolution ALMA maps. As a result, we expect to conduct a more extensive analysis at even higher physical resolutions, 50–100 pc, as the survey proceeds.

We measure typical CO-to-mid-IR ratios (blue points, Figure 14) of I_CO/I_MIR ≈ 0.8–2.7 K km s⁻¹ (MJy sr⁻¹)⁻¹. This agrees broadly with the "rule of thumb" observed for large parts of galaxies that a CO intensity of 1 K km s⁻¹ corresponds to about 1 MJy sr⁻¹ in the mid-IR (see also Leroy 2023). After accounting for the mean band ratios in Table 1, the typical values expressed in units of 7.7 μm intensity, ${I}_{\mathrm{CO}}/{I}_{\mathrm{MIR}}{{ \mathcal R }}_{7.7}^{X}\,\approx \,$ 0.8–1.3 K km s⁻¹ (MJy sr⁻¹)⁻¹.

The template-fit coefficients (green points) imply higher CO-to-mid-IR when we isolate the CO-tracing component. Here the ratios span c⁻¹ ≈ 1.4–4.7 K km s⁻¹ (MJy sr⁻¹)⁻¹, i.e., about twice the overall ratio, reflecting that only half the mid-IR flux in the model fit has been assigned to the CO-tracing component. When placed on a common 7.7 μm intensity scale yield, the template fits yield ${c}^{-1}{{ \mathcal R }}_{7.7}^{X}\,\approx \,$ 1.3–2.7 K km s⁻¹ (MJy sr⁻¹)⁻¹.

For comparison, Equation (6) suggests that for molecular gas heated by a solar neighborhood radiation field (U = 1) and a typical dust-to-gas ratio:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{I}_{\mathrm{CO}}}{{{\rm{I}}}_{\mathrm{MIR}}} & \sim & 5.7\,{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\left(\mathrm{MJy}\,{\mathrm{sr}}^{-1}\right)}^{-1}\\ & \times & {\left({{ \mathcal R }}_{7.7}^{X}\,\left(\displaystyle \frac{{\alpha }_{\mathrm{CO}\ 1-0}}{{\alpha }_{\mathrm{CO}\ 1-0}^{\mathrm{MW}}}\right)\left(\displaystyle \frac{{R}_{21}}{0.65}\right)\left(\displaystyle \frac{D/G}{0.01}\right)\left(\displaystyle \frac{{q}_{\mathrm{PAH}}}{0.046}\right)\left(\displaystyle \frac{U}{{U}_{0}}\right)\right)}^{-1},\end{array}\end{eqnarray} \tag{ 14 }$$

with a prefactor about five times higher than our measured ratio and ≈2.5 times higher than our typical template-fit coefficient.

For the template-fit coefficients, the ISRF factor, U, almost certainly explains most of the difference from the fiducial expectation. Studying the whole disk of a local sample of star-forming galaxies, Aniano et al. (2020) find a typical mean ISRF in the range U/U₀ ∼ 2 − 3. For the one target where we overlap their study, NGC 0628, they find both the diffuse and mean ISRF to be U/U₀ ∼ 2. Given that PHANGS–JWST targets the inner part of star-forming galaxies and the ISRF tracks star formation activity and stellar density, it seems very likely that in more diffuse regions U/U₀ ≳ 2 and plausibly higher. Given this, our reference comparison line in Figure 14 shows U/U₀ = 2.5 and this appears to be both physically reasonable and represent a good description of the best-fitting coefficients.

Logically, the difference between our higher c⁻¹ and the lower measured median ratios should stem from the inclusion of regions of recent star formation and intense heating in the median ratio. The template fit will assign much of the flux associated with these regions to the Hα-tracing component, but they are included in the overall ratio.

We measure extinction-corrected Hα-to-mid-IR ratios of ${\mathrm{log}}_{10}H\alpha /{I}_{\mathrm{MIR}}\sim -5.7$ to −5.2 (red points, Figure 14) in units of erg cm⁻² s⁻¹ sr⁻¹ (MJy sr⁻¹)⁻¹. If we adjust these for the band ratios in Table 1 all bands yield a fairly consistent median ${\mathrm{log}}_{10}{\rm{H}}\alpha /{I}_{\mathrm{MIR}}\,{{ \mathcal R }}_{7.7}^{X}\,\sim \,$ −5.6 to −5.4. As with the CO case, the template-fit coefficients yield higher values, ${\mathrm{log}}_{10}{h}^{-1}\,{{ \mathcal R }}_{7.7}^{X}\,\sim \,$−5.2 to −4.9, with ${\mathrm{log}}_{10}{h}^{-1}\,{{ \mathcal R }}_{7.7}^{X}\approx -5.0$ on average. This reflects that these template-fitting values do not include emission associated with the CO-tracing component.

Following Equations (5) and (7) the normalization of the Hα-to-mid-IR ratio can be expressed in terms of the coefficient C₂₄ to translate mid-IR to SFR, with ${\mathrm{log}}_{10}H\alpha /{I}_{\mathrm{MIR}}\,\sim \,$ −5.2 expected at 7.7 μm for a standard ${\mathrm{log}}_{10}{C}_{24}\,\sim \,$ −42.7 often used for whole galaxies or large parts of galaxies (e.g., see Kennicutt & Evans 2012; Belfiore et al. 2022b). The line that goes through the purple points in Figure 14 shows a value two times lower than this, ${\mathrm{log}}_{10}{C}_{24}\,\sim \,$ −43.0, implying that two times higher SFR-to-mid-IR ratio describes the template-fit values of h compared to whole galaxies or large regions. This result agrees qualitatively with the results of Belfiore et al. (2022b) for the brightest regions and for the overall literature, which finds lower C₂₄ when focusing on the most intense star-forming regions rather than whole galaxies (e.g., Kennicutt & Evans 2012; Calzetti 2013).

Our median ratios correspond to somewhat lower SFR-to-mid-IR ratios, or higher C₂₄, than this typical value. This partially reflects our methodology. We plot median ratios, but the Hα-to-mid-IR relation (Section 4.2) is superlinear with much of the flux in a set of bright regions. Our approach down-weights these bright regions and emphasizes the lower intensity regions, which tend to be dominated by emission associated with molecular gas (see Figure 12). More generally the lower overall ratio reflects the fact that averaging over whole galaxies, emission from dust associated with molecular gas and other "diffuse" emission will be averaged in with emission directly from star-forming regions. This tends to be calibrated out by SFR conversion recipes at low resolution yielding the sort of higher C₂₄ and lower SFR-to-mid-IR ratios that we observe here.

IC 5332 has lower stellar mass and only 0.6 times the metallicity of our other three targets (Table 2). Consistent with the general trends for lower stellar mass galaxies to have lower star formation activity and a lower molecular to atomic gas fraction (e.g., Saintonge & Catinella 2022), the galaxy shows overall lower surface brightness in the mid-IR and in CO than the other three targets. We expect the dust-to-gas ratio, D/G, to correlate linearly with metallicity to first order (e.g., see Galliano et al. 2018). Therefore we expect D/G to also be ≈0.6 times lower in IC 5332 than the other three targets.

Likely reflecting this lower dust abundance, Table 4 and Figure 8 do show that IC 5332 exhibits a moderately higher Hα-to-mid-IR ratio than the other galaxies, as one would expect if the lower dust content leads to lower overall UV extinction and less IR emission. On the other hand, the CO-to-mid-IR ratio in IC 5332 closely resembles that in the other three galaxies, especially for F770W, F1000W, and F1130W. We note that this is effectively a stacking result because IC 5332 is only weakly detected in PHANGS–ALMA. Figure 6 shows that if we stack our "flat" CO cubes (Section 3.1) as a function of mid-IR intensity, we observe a linear relationship with a similar normalization to the other three galaxies. This result qualitatively resembles the observation of a good correlation between PAH and CO emission in the Magellanic Clouds (e.g., Sandstrom et al. 2010; Chastenet et al. 2019).

This might appear surprising given IC 5332's expected lower D/G compared to the other targets. In theory, a lower D/G should yield a lower IR emission per gas column (Equation (4)). The quantitative agreement between IC 5332 and the other galaxies may emerge because the CO-to-H₂ conversion factor α_CO (Equation (2)) also depends on metallicity (e.g., Glover & Low 2011; Leroy et al. 2011a; Bolatto et al. 2013; Hu et al. 2021). The sense of that relationship is that the CO emission per unit H₂ drops at low metallicity. Our results suggest that the drop in D/G may cancel with a rise in α_CO in IC 5332 to produce a CO versus mid-IR correlation similar to that found in higher-metallicity galaxies. With proper calibration this might offer a mid-IR-based approach to explore α_CO variations, something explored by Gratier et al. (2010) using Local Group Spitzer data. However, we note that such low surface densities, dust mixed with H i, or an extended CO-dark phase may also play a role (see Sandstrom 2023b for more discussion and quantitative exploration).

One main result from our analysis is that the mid-IR traces the molecular gas surprisingly well. This implies that one can use the mid-IR to trace gas column density, analogous to the way that longer-wavelength dust emission has long been used to trace the ISM in galaxies at low and high redshift (Israel 1997; Leroy et al. 2011a; Eales et al. 2012; Sandstrom et al. 2013; Scoville et al. 2014; Shi et al. 2014; Groves et al. 2015; Tacconi et al. 2020).

Using JWST mid-IR observations to trace the ISM represents an exciting prospect because these observations are sensitive, high resolution, and relatively phase-agnostic. As pointed out in Section 3.3 and Table 3, if we simply adopt our median CO-to-mid-IR ratio, then at matched resolution, our relatively short integration time JWST observations have significantly better mass sensitivity, σ_gas ∼ 0.2 M_⊙ pc⁻² at F770W, F1130W, or F2100W, than our (also relatively short integration time) ALMA observations, σ_gas ∼ 6.7 M_⊙ pc⁻². Second, JWST naturally achieves very high angular resolution. As an interferometer, ALMA often requires large time investments to achieve surface brightness sensitivity at high angular resolution. Mapping a single galaxy in the relatively bright CO lines at resolution matched to JWST at 11.3 μm, 7.7 μm, or even the 3.3 μm PAH feature (see Sandstrom 2023a) comes with prohibitive cost. Third, none of our expectations (Section 2) specifically require that mid-IR emission arise from dust mixed with molecular gas. We have simply focused our own comparison on CO because we have an independent tracer of that phase of the ISM. As pointed out and explored by Sandstrom (2023b) in this Issue, the dust traced by the mid-IR should have a broader phase sensitivity than just CO-emitting molecular gas, including both CO-dark gas and atomic gas.

Practically, what should one do to predict the gas distribution from the mid-IR? Before giving recommendations, we first emphasize again that these are early days for these kinds of observations, and we expect our approach to evolve. That said, to first order, one can adopt a purely empirical approach to estimate the gas column density. One can use either our measured median ratios or our power-law fits (both in Table 4) to predict I_CO from mid-IR emission (the median ratios may be slightly less accurate, but are also less scale dependent). Then one can convert to Σ_gas = Σ_mol from Equation (2), leveraging the fact that our data selection was designed to pick regions where most of the gas is likely to be molecular so that we can treat the result as a general relation between the total gas column density and the mid-IR.

We state the result in this way because we expect that the mid-IR traces the total neutral gas column. Dust is mixed with all gas, not only with molecular gas (e.g., see Bohlin et al. 1978; Sandstrom et al. 2013; Galliano et al. 2018), and we assume that to first order that the combination of the dust-to-gas ratio, PAH abundance, and ISRF strength remains about fixed for the gas outside star-forming regions. The ambient ISRF is indeed observed to vary smoothly based on resolved SED modeling of nearby galaxies (Aniano et al. 2012, 2020; Draine et al. 2014), though note that it does vary, showing, e.g., radial variations that track the overall stellar density. The PAH abundance outside star-forming regions also appears to vary relatively smoothly (e.g., Chastenet et al. 2019; Chastenet 2023a, the latter in this Issue), though again showing large scale variations. And the dust-to-gas ratio also shows coherent, slowly varying behavior across galaxies, tracking the metallicity well (e.g., Sandstrom et al. 2013; Draine et al. 2014; Chiang et al. 2018; Galliano et al. 2018). We note an important caveat here, that there is good evidence for a density or phase dependence of the dust-to-gas ratio (e.g., Jenkins 2009; Roman-Duval et al. 2014; Chiang et al. 2018). This has the sense that gas with very low densities could have lower dust-to-gas ratios, but the variations are relatively modest in magnitude and the dust-to-gas ratio in the dense, cold atomic phase is likely to more closely resemble that in molecular gas (see Jenkins 2009; Draine 2011). We leave a detailed treatment of this effect, along with a deeper investigation of the phase dependence of the ISRF or q_PAH, for future work, noting that even with an associated uncertainty of ∼50%, the mid-IR still represents a powerful new ISM tracer.

In this approach, bright star-forming regions will represent an important contaminant. One can mask them imprecisely by clipping on mid-IR brightness, e.g., using the approximate mid-IR cuts noted in Section 2.6 based on Belfiore et al. (2022a). Better, one can use a Hα map and the results of this paper to estimate where local star-forming regions make an important contribution to the mid-IR (i.e., scale the Hα, perhaps after decomposition into H ii regions, using our results and compare to the local mid-IR to gauge where local contamination rises above some user-specified tolerance). The downside of masking the star-forming regions is that one will then miss molecular gas associated directly with these regions. To address this, one could conduct a version of our template-fitting exercise and attempt to subtract the Hα-tracing component. With a broader set of JWST targets, it should also be possible to identify general prescriptions to subtract the Hα. Though subtracting bright emission or otherwise modeling the local radiation field to infer the local column density is the most appealing route forward, it remains to be seen if this approach can deliver accurate estimates of gas even in the presence of bright star-forming regions.

Physically, as outlined in Section 2.5 and Equations (4), the mapping between mid-IR emission and gas surface density depends on the dust-to-gas ratio and the ISRF, as well as the PAH abundance, and in the short term, we can hope to take a more physical and less empirical approach to this topic by directly estimating all of these quantities. The PHANGS data sets contain information on the resolved recent star formation, older stellar content, star formation history, and metallicity of our targets. The JWST data themselves constrain variations in the PAH abundance via the band ratios. In the near future, the combined PHANGS data sets should allow us to make predictions for how the ISRF, D/G, and q_PAH vary across each of our targets, and use these to calibrate local conversions from mid-IR to column density. In parallel, we expect our ability to identify and subtract emission from bright star-forming regions to improve, and both of these avenues should lead to higher-quality estimates of the gas column density.

The classical application of mid-IR emission is to trace star formation. Extinction-corrected Hα does correlate very well with mid-IR emission in our data. There remains a great deal of work to do on this topic with these data. Here we note that from this work, the power-law fits relating extinction-corrected Hα-to-mid-IR emission (Table 4) represent our most straightforwardly useful result. They capture the nonlinearity in the overall Hα versus mid-IR relations well, and the results for the three massive disk galaxies agree very well (though we clearly also expect metallicity effects; see Section 6.3). The fits have been derived for mid-IR intensities of ∼0.5–30 MJy sr⁻¹ and physical scales of 50–150 pc, and we expect them to perform best in this range. After estimating the extinction-corrected Hα, one can apply Equation (3) to estimate the SFR but should recall that our analysis used the entire Hα map and not explicitly separated the DIG and H ii regions (see Belfiore et al. 2022a; Groves et al. 2023).