Star Formation and Molecular Gas Diagnostics with Mid- and Far-infrared Emission

In the following we test a variety of MIR features and their joint correlations with SF and gas; here we briefly review some key points about them and their origins.

1.1.1. Tracing SFR

1.1.2. Tracing H₂

The SINGS program observed nuclear and extranuclear H ii regions in 75 nearby galaxies with the Short-Low (SL2 5.2–7.7 and SL1 7.4–14.7 μm), Short-High (SH 9.9–19.4 μm), Long-Low (LL 14–38 μm), and Long-High (LH 19.1–37.1 μm) modules of the InfraRed Spectrograph (IRS) on the Spitzer Space Telescope (Kennicutt & Armus 2003; Houck et al. 2004; Smith et al. 2007; Dale et al. 2009). We use data from SINGS Data Release 5. ⁴ The LL spatial coverage of extranuclear regions in SINGS is sparse compared to SH and LH, so we make use of LH for the long wavelength coverage. Our sample is therefore composed of the overlapping areas between the SL, SH, and LH data for each nuclear and extranuclear target. We note that this data set is focused on optically bright Hii regions in local spiral galaxies. This biases our sample toward more massive and evolved H ii regions.

A similar data set was defined in Roussel et al. (2007) in order to study the H₂ rotational lines. Our selection is nearly the same as that in Roussel et al. (2007) except we use the Data Release 5 of the SINGS data and the apertures are defined independently. The resulting spectral range includes all the major PAH features (except 3.3 μm, which contributes ∼2% of the total PAH power; Lai et al. 2020) and the 12.8 μm [Ne ii] and 15.6 μm [Ne iii] lines that are needed for our SFR tracer (see Section 2.6). We extract spectra in all IRS orders from the same apertures. The full data set of extraction apertures consists of 154 polygonal apertures with numbers of vertices ranging from 4–7. A table with the vertices of each region and all data used in the following analysis is available in machine-readable form. A description of the contents of this table is shown in Appendix A.1.

We use Spitzer, Wide-field Infrared Survey Explorer (WISE), and Herschel IR photometric observations from 3.4–250 μm from SINGS, the "z = 0 Multiwavelength Galaxy Synthesis" (z0MGS; Leroy et al. 2019), and "Key Insights on Nearby Galaxies: a Far-Infrared Survey with Herschel" (KINGFISH; Kennicutt et al. 2011). We use the SINGS photometric data from the four InfraRed Array Camera (IRAC) bands at 3.6, 4.5, 5.8, and 8 μm, and the 24 μm Multiband Imaging Photometer for Spitzer (MIPS) band. We did not apply extended source correction factors to this data since they do not affect our correlation analysis. The z0MGS data includes the four WISE bands at 3.4, 4.6, 12, and 22 μm. We also include FIR KINGFISH data from the three Photodetector Array Camera and Spectrometer (PACS) photometric bands at 70, 100, and 160 μm, as well as the Spectral and Photometric Imaging Receiver (SPIRE) 250 μm photometric band (Kennicutt et al. 2011).

We convolve all photometric data to 15'' to match the resolution of our CO data, except SPIRE 250 μm which has 18'' resolution (see Section 2.4). We perform a simple background subtraction on each photometric image using the median value of the image outside the R₂₅ isophotal radius of the galaxy with a 3σ mask to filter out background galaxies and foreground stars. We multiply each photometric band by its characteristic frequency in order to match the units of integrated intensities for spectral emission features (10⁻⁷ W m⁻² sr⁻¹) to facilitate direct comparisons between photometry and spectroscopy. We also use the photometric bands longer than 20 μm to calculate TIR using the calibration adapted from Galametz et al. (2013; see Section 2.6). We assume a constant 10% of the measured surface brightness in each aperture for all photometric bands. We find our results do not vary significantly if this assumption is doubled.

The molecular gas data used in this work consists of CO J = (2 − 1) measurements from the "Physics at High-Angular resolution in Nearby GalaxieS"—Atacama Large Millimeter Array (PHANGS-ALMA) and the "Heterodyne Receiver Array CO-Line Extragalactic Survey" (HERACLES) from the compilation in Leroy et al. (2022). These CO maps have been convolved to a Gaussian point-spread function (PSF) with FWHM of 15''. We leave the maps at this resolution in order to match the resolution of FIR maps from KINGFISH. Of the 75 galaxies in the SINGS sample, we have a CO J = (2 − 1) map for 41. Each map has a corresponding error map from which we extract the same IRS overlap apertures to determine the uncertainty on the corresponding CO emission. Of the 154 regions in our IRS spectroscopy sample, 104 have a CO measurement, and 89 of these have detections of CO with a signal-to-noise ratio (S/N) >3.

We use CO emission as our reference molecular gas tracer and we do not convert to H₂ gas masses with a CO-to-H₂ conversion factor. Since the CO-to-H₂ conversion is approximately the same for most of our regions, this choice will not affect the correlation coefficients. We explore the expected effect of translating CO emission to M_H2 in Section 4.2.

We extract a single average spectrum from all Spitzer-IRS orders in our apertures based on the overlap of SL1, SL2, SH, and LH spatial coverage. The size of these apertures is comparable to the angular resolution of the CO data (e.g., ∼15'' in diameter). To ensure matched extractions for the photometry, we convolve each IRAC, WISE, MIPS, and PACS image to a Gaussian PSF with FWHM of 15'' matched to the CO maps using convolution kernels from Aniano et al. (2011). We note SPIRE 250 μm has a slightly lower resolution (∼18'') that cannot be matched to 15'' but is close enough that we do not expect major issues, so we leave it at the native PSF.

For the spectroscopic data, the coverage of the spectral maps is small, so we cannot convolve to matched 15'' resolution without introducing artifacts. Our apertures are larger than the PSF of the IRS data at all wavelengths (Pereira-Santaella et al. 2010), which should minimize any PSF effects, so we proceed to extract and average the spectra in these regions without PSF matching. To check that PSF effects are minimal and no aperture corrections are needed for the IRS spectral extraction, we compare the native resolution Spitzer photometric maps (which have similar PSFs to the spectroscopy) to 15'' resolution convolved images at the IRAC 5.8 and 8 μm bands and the MIPS 24 μm band. We extract photometry in our apertures from the native and convolved images and determine the native-to-convolved ratio for each region. The native-to-convolved ratio is found to vary little with wavelength: 1.1 ± 0.1 for each photometric band. We then interpolate the ratio as a function of wavelength from 5.8–24 μm to estimate the effect of using the native resolution spectra in comparison with the 15'' photometric maps. We find our final results are not altered if we use this aperture correction on the IRS data, so we proceed without applying any correction.

We determine metallicities for our apertures using the central values and gradients from the Kobulnicky & Kewley (2004; KK04) calibration found in Moustakas et al. (2010). We find no significant difference in our observed trends in Section 3 if the Pilyugin & Thuan (2005; PT05) calibrated metallicities are used instead of the KK04 values. For galaxies with no calculated gradient we use the characteristic value, and for those with no characteristic we use the luminosity–metallicity (L-Z) value. These metallicities are used in combination with the 12.8 μm [Ne ii] and 15.6 μm [Ne iii] lines to determine Σ_SFR according to the prescription from Whitcomb et al. (2020; see Section 2.6), which was also calibrated with KK04 metallicities from Moustakas et al. (2010).

We use nuclear classifications from Murphy et al. (2018) and Moustakas et al. (2010) to identify apertures that include an active galactic nucleus (AGN). These regions are excluded from our results. For simplicity, we consider galaxies indicated as "SF/AGN" (see Table 1) as AGN and exclude their centers as well. From the 154 regions in our sample, only 10 do not have significant detections of either of the two neon lines required for our SFR calibration (see Section 2.6). Of the remaining 144, 32 are centers of galaxies containing an AGN. We remove these AGN from our sample, leaving a total of 112 SF regions and of these there are 75 that are covered by CO mapping and 67 that have CO detections at S/N >3. We use the full set of 112 regions for SFR correlations and CO correlations only use the subset of 67 regions. Including the weaker detections as upper limits does not significantly constrain the slope and intercept of our power-law fits, so we exclude them from the analysis.

We also determine the deprojected physical area of our apertures in order to investigate how spatial scale affects correlations between MIR, CO, and SF. This requires the galaxy inclinations, position angles, and distances listed in Table 1. The median area of the 112 regions is about 0.4 kpc² and likewise for the 67 of these regions with CO data. We compare the behavior of correlations involving the half of regions greater than 0.4 kpc² to those with smaller areas in Section 3.1. Since the SINGS apertures were placed based on peaks of Hα emission (Kennicutt & Armus 2003), we expect that in at least half the regions with area <0.4 kpc² the Σ_SFR is resolved distinctly from the CO and the two are spatially uncorrelated (Schinnerer et al. 2019).

From the above-described data, we derive TIR for each of our apertures using the calibration presented in Galametz et al. (2013):

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{TIR} & = & 2.013\times (\nu {S}_{24})+0.508\times (\nu {S}_{70})\\ & & \,+0.393\times (\nu {S}_{100})+0.599\times (\nu {S}_{160})\\ & & \,+0.680\times (\nu {S}_{250}),\end{array}\end{eqnarray} \tag{ 1 }$$

where the constants are taken from their Table 3, S_λ is the surface brightness at λ microns in MJy sr⁻¹, and ν is the frequency corresponding to the characteristic wavelength of each band. We use this calibration for all galaxies for simplicity instead of applying galaxy-specific calibrations where they are available.

We also derive SFR for each of our apertures using the calibration presented in Whitcomb et al. (2020). This tracer was calibrated to match free–free radio continuum emission using a similar data set based on SINGS extranuclear spectra. The SFR surface density Σ_SFR is given by

$$\begin{eqnarray}\begin{array}{rcl}\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{SFR}}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-2}}\right) & = & \,6.2\times {10}^{-3.25}\\ & & \times \,{\left({\rm{\Sigma }}\mathrm{Ne}\right)}^{0.87}\times {\left([{\rm{O}}/{\rm{H}}]\right)}^{-0.48},\end{array}\end{eqnarray} \tag{ 2 }$$

where constants are taken from their Table 5, Σ Ne is the sum of the 15.6 μm [Ne iii] and the 12.8 μm [Ne ii] integrated line intensities in 10⁻⁷ W m⁻² sr⁻¹, and [O/H] is from KK04-calibrated 12 + log₁₀[O/H] metallicities. We use this as our reference SFR tracer and refer to it as the ΣNe and Z tracer. This tracer is based on fine structure lines from high ionization potential ions (I > 13.6 eV) and it is calibrated to free–free emission. Our reference tracer is then tracing ionizing photons generated by recent star-forming events (≲5 Myr). Since our reference tracer responds to SF events on short timescales, an emission feature could be less correlated with our SFR if it responds directly to SF on longer timescales. The ΣNe and Z tracer is also independent of dust emission and attenuation. It is crucial that our SF tracer is independent of any dust attenuation or obscuration correction, since our goal is to separate the correlation of various MIR observables with SF and gas separately.

We also include 21 regions that have detections of only one neon line, with the undetected line set to zero. Our results are not altered by including these regions and they typically have the lowest Σ_SFR in our sample but do not deviate from observed trends (see Figure 3). Equation (2) gives the SFR surface density Σ_SFR, and we refer to this quantity as simply SF or SFR interchangeably for the remainder of this paper.

The extracted SL spectra were fit with the IDL program PAHFIT (Smith et al. 2007). We input heliocentric systemic velocities from the Lyon–Meudon Extragalactic Database (LEDA; Makarov et al. 2014) to correct for redshift variations. Constraints on the attenuation are weak with only SL coverage and previous studies of SINGS regions have found that, excepting specific cases such as edge-on targets, the attenuation in the MIR is generally minimal (Roussel et al. 2007; Smith et al. 2007), so we do not fit the attenuation. From PAHFIT we obtained the integrated intensities of emission features, as well as PAH, line, and continuum component spectra. The subfeatures of major PAH complexes are combined such that we obtain integrated intensities for the following bands and band complexes: 6.2, 7.7, 8.3, 8.6, 11.3, 12.0, 12.6, 13.6, and 14.2 μm. Uncertainties for these bands and emission lines are also returned by PAHFIT.

At present PAHFIT only operates on low-resolution spectra. To separate the underlying hot dust continuum from PAH emission in high-resolution SH spectra, we constructed a program similar to PAHFIT following Smith et al. (2007). We modeled fine structure and molecular line emission with Gaussian profiles. PAH features were modeled with the Drude emission profiles. We subtracted two linear estimates of the continuum between three points with no nearby emission features: one between 10 and 15 μm, and another between 15 and 18.2 μm. We then fit the spectrum from 10–15 μm as a linear combination of the individual Gaussian and Drude profiles using initial estimates and constraints for the parameters shown in Table 2 adopted from Smith et al. (2007), except for those that compose the 11.3 μm PAH complex. We find this complex is better fit by four Drude profiles in the higher spectral resolution SH data than the analogous two Drude profiles prescribed in Smith et al. (2007), which are designed for lower spectral resolution SL data. Notably, the 11.0 μm PAH feature can be fit distinctly from the rest of the complex at 11.3 μm, which we find is best described by PAH components at 11.2, 11.25, and 11.4 μm. The same technique is also applied to the spectral features between 15 and 19 μm.

We obtain uncertainties on the integrated intensities with Monte Carlo trials. In each trial, we perturb each point in the SH spectrum within its uncertainty assuming a Gaussian distribution, then run our fitting routine to find the integrated intensity of each emission feature. We then take the mean and standard deviation of 500 trial fits as the value and uncertainty of the feature. For convenience, we refer to this method and the program used to implement it as the Fitter for Aromatic, Atomic, and Molecular features in the Infrared (FARAMIR). Figure 1 (center) shows the FARAMIR fit for an SH observation of an extranuclear region in NGC 5194. In Appendix A.2 we verify the accuracy of FARAMIR by comparison with PAHFIT results for analogous SL spectra. FARAMIR reproduces the relevant PAH feature and emission line intensities but overestimates by ∼20% or less on average relative to PAHFIT. This overestimation is likely a result of the simple linear continuum model used by FARAMIR, whereas PAHFIT uses a set of several blackbody functions to model the dust continuum.

Zoom In Zoom Out Reset image size

The LH spectra contain only three emission lines that can be fit in most regions: H₂ S(0) at 28.2 μm, [S iii] at 33.5 μm, and [Si ii] at 34.8 μm. Beyond 35 μm the S/N is too low to reliably fit emission lines such as [Ne iii] at 36.0 μm. The other notable lines in the LH spectral range, [O iv] at 25.91 μm and [Fe ii] at 25.99 μm, are too faint relative to the LH periodic continuum noise (scalloping ⁵ ) to be fit separately in most regions. We fit the LH spectra with three individual Gaussian profiles: for each line we crop the spectrum to 1 μm on either side of the peak wavelength, then fit a linear continuum using the average surface brightness between 1 and 0.5 μm on either side of the peak wavelength. This continuum is subtracted from the cropped spectrum and we follow the same Monte Carlo procedure outlined above to obtain integrated intensities and their uncertainties.

We quantify the detection threshold of our fits to emission features in the SH and LH spectra using the residual scatter after the fit is subtracted from the original data. We calculate the standard deviation of the residual difference and compare it with the fitted amplitudes determined by our method. We define our non-detection threshold as thrice the standard deviation of the residuals, such that any Drude or Gaussian with a fitted amplitude below this threshold is considered a non-detection.

We find very few detections for the 15.9 and 17.9 μm PAH features, likely because these are intrinsically weak and immediately near our assumed continuum points at 15 and 18.2 μm. Our linear continuum assumption for fitting SH spectral features is sufficient for most bright spectra from star-forming regions, but some have a significantly nonlinear continuum and a small but non-negligible portion of the PAH features can be removed by the linear fit. Likely as a result of this linear fit, we find less than fifteen 3σ detections of the 13.5 and 14.2 μm PAH features with our FARAMIR method from SH spectra, but this is not the case for these same PAH features as measured by PAHFIT from SL spectra. When necessary we use the PAHFIT results from SL spectra for these low S/N PAH features in our analysis. We also find less than fifteen detections for the 17.9 μm PAH feature. All major emission features in the SL spectral overlap with SH measured with PAHFIT and FARAMIR respectively are in good agreement (see the Appendix).

In order to study the contributions of PAH features, emission lines, and dust continuum to the IRAC4 8 μm and WISE3 12 μm photometric bands, we perform synthetic photometry on decomposed IRS spectra. The WISE3 12 μm photometric band is sensitive to emission from 7.2–16.4 μm. This range overlaps with the IRAC4 band at 8 μm, which collects from 6.2–9.7 μm. The SL1 and SL2 combined spectra span the full IRAC4 band, but the SL spectral range ends at 14.7 μm. For the remaining portion of the WISE3 wavelength coverage outside the IRS SL range, we append the 14.7–16.4 μm range of the corresponding SH spectrum. The continuum level differs between SL and SH, so we subtract a linear continuum from the SH as determined by the average brightness at 10.1 and 14.4 μm. We then match the continuum-subtracted SH spectrum to the SL spectrum based on the average SL brightness at these same wavelengths. The final stitched spectra are composed of the full SL spectrum from 5.2–14.7 μm and the continuum-adjusted SH spectrum from 14.7–16.4 μm. These combined spectra are only used for our synthetic WISE3 photometry analysis.

We use PAHFIT to decompose the SL1 and SL2 spectra into PAH, starlight continuum, dust continuum, and line emission contributions in the 120 non-AGN spectra. For the wavelengths longer than the SL coverage (greater than 14.7 μm) we consider the appended 14.7–16 μm portion of the SH spectrum to be continuum. The WISE3 bandpass is sensitive out to about 16 μm and other than continuum, the 15.6 μm [Ne iii] line is the only emission feature included between 15 and 16 μm. We find it contributes a negligible amount to the total integrated flux from 15–16 μm compared to continuum emission for the majority of regions.

From the decomposed SL spectra we find the IRAC4 8 μm photometric band is dominated by emission from the PAH complex at 7.7 μm. But there is a roughly equal split between PAH and dust continuum portions for the full WISE3 12 μm band. Figure 2 shows the dust continuum contributes a larger fraction of emission to the WISE3 band than PAHs. The median fractional contribution to IRAC4 from PAHs is found to be 78%, that from continuum about 19%, and that from lines about 4%. The median fractional contribution from PAHs to WISE3 is 38%, that from continuum is 58%, and that from lines is about 4%. We note that the PAH contribution to WISE3 only exceeds 50% in a handful of regions.

Zoom In Zoom Out Reset image size

We quantify the correlation between two emission features with the Spearman rank correlation coefficient. We also fit a power law to the data by ordinary least-squares using the logarithmic form:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({\rm{Y}})=m\times {\mathrm{log}}_{10}({\rm{X}})+{\rm{b}}\end{eqnarray} \tag{ 3 }$$

where Y is CO [K km s⁻¹] or Σ_SFR [M_⊙ yr⁻¹ kpc⁻²], X is the comparison emission feature converted to surface brightness units 10⁻⁷ W m⁻² sr⁻¹ for all quantities (other than CO and Σ_SFR), m is the power-law slope, and b is the power-law intercept.

We use a Monte Carlo method to account for uncertainties in each axis. In each of the trials we randomly offset each point according to its uncertainty assuming a Gaussian distribution. We simultaneously perform a bootstrap that randomly resamples the whole data set with replacement. Then we fit the constants of the power law and calculate the standard deviation and Spearman correlation coefficient. We complete 500 such trials and take the mean and standard deviation as the value of each quantity and its uncertainty.

Figure 3 shows the strongest correlations with CO and SF as an example. We also include examples of weaker correlations by showing the CO correlation for the feature best correlated with SF, and the SF correlation for the feature best correlated with CO. We note that the SF correlations have a different number of points than the CO correlations. The number of points varies based on the number of detections of each emission feature, but the maximum possible is 112 for SF and 67 for CO correlations. We find that restricting the set to include the same points in both correlations does not significantly alter our results.

The color gradient in Figure 3 shows the ratio of 15.6 μm [Ne iii] to 12.8 μm [Ne ii]. This ratio is a tracer of ionization parameter, which is itself strongly dependent on local properties such as the hardness of the radiation field, and therefore metallicity. The relation between [Ne iii]/[Ne ii] and metallicity has been investigated in many other studies (e.g Madden et al. 2006; Ho & Keto 2007; Xie & Ho 2019). Whitcomb et al. (2020) showed there is a strong, negative correlation between [Ne iii]/[Ne ii] and metallicity for apertures extracted from the same SINGS SL spectra used here.

Zoom In Zoom Out Reset image size

In our effort to quantify how well IR emission features trace SF and molecular gas, the breakdown of the K-S relationship on small scales provides some unique insights into the correlations. On the smallest scales, H₂ and ionized gas from SF are not colocated. However, they are often found nearby in the dense star-forming complexes studied here, so on large scales the correlation between SFR and CO is strong and it becomes difficult to distinguish between tracers of SF and tracers of CO (Schruba et al. 2010; Leroy et al. 2013; Schinnerer et al. 2019). Our regions sit at intermediate size scales (∼100 s pc, but less than a kiloparsec), so there can be some distinction in the distributions. At the smallest sizes this distinction becomes even clearer.

Comparing the left and right plots in Figure 5, the dashed lines show the correlation of Σ_SFR with Σ_H2 is significantly stronger for the 34 regions with an area greater than 0.4 kpc² compared to the correlation for the other 33 smaller regions (ρ_s = 0.7 versus 0.2). The dashed lines for the small regions shown in Figure 5 (right) are less than the critical Spearman coefficient so the correlation between SF and CO in these regions is consistent with the null hypothesis. We note that the K-S correlation is likely weaker in general for our sample since our regions cover relatively small portions of many different galaxies.

Zoom In Zoom Out Reset image size

While dividing the sample by size reveals the behavior of various tracers more clearly, it does introduce some biases in the selection as well. In particular, since the SINGS apertures are all of similar angular size, the smallest regions will be from the closest sources, which are more likely to be dwarf and/or lower metallicity galaxies. The relative positions of the various tracers in Figure 4 are mostly unaffected by looking at smaller size scales in Figure 5 (right), indicating that larger apertures simply have more Σ_SFR–CO overlap. We therefore draw our conclusions from the more statistically robust sample of all SF regions. We define an emission feature as SF dominant if it has greater ρ_SFR than ρ_CO, and CO dominant for the opposite case.

All major PAH features and H₂ rotational lines are strongly correlated with CO (ρ_CO > 0.8). From the other plots of Figure 4 we find the PAH features and H₂ rotational lines also have about the same sublinear power-law slope and the same amount of residual scatter in the CO correlation (both being within their 1σ uncertainties). The exception is the PAH 17 μm complex, which has a slope consistent with linear and a greater amount of residual scatter in the CO correlation compared to the shorter wavelength PAH features.

All photometric bands from 3.4–250 μm are well correlated with CO. We find the bands at 24 and 70 μm have the lowest ρ_CO (0.7) while the MIR 5.8 and 8 μm and FIR 100, 160, and 250 μm bands have the highest ρ_CO (0.85). The TIR emission calculated from our photometric data has approximately the same ρ_CO and ρ_SFR as the 24 and 70 μm bands, implying that TIR for these regions reflects the warm dust in the vicinity of SF.

We find the weakest correlations with CO are fine structure lines from highly ionized species. For example, [Ne iii] and [S iv] have ρ_CO less than the critical Spearman coefficient 0.17, indicating the null hypothesis applies and they are not correlated with CO. These highly ionized lines and Humphreys-α are the only emission features with correlations weaker than the minimum K-S relation. This is possible because of our sample selection of regions focused on bright ionized gas associated with recent SF. This behavior is expected since emission from the most highly ionized gas is not correlated with molecular gas emission on small scales. We discuss this further in Section 4.3.

All features and photometric bands have a stronger correlation with Σ_SFR than the CO:Σ_SFR correlation (ρ_SFR > 0.4, see Section 4). We find the [S iii] lines and their sum, total [S iii], are best correlated with Σ_SFR with Spearman coefficient ρ_SFR ∼0.95. Close behind at ρ_SFR = 0.8 are the 15.6 μm [Ne iii] and 12.8 μm [Ne ii] lines that compose our reference Σ_SFR, as well as the MIPS 24 μm band, which is a widely used SF indicator (Calzetti et al. 2007). The [Si ii] line and PACS 70 μm emission also have ρ_SFR ≈0.8. The weakest Σ_SFR correlations are the H₂ S(3) line, and the short-wavelength photometric bands WISE1 3.4 μm and WISE2 4.6 μm with ρ_SFR ∼ 0.5 for each. These features are less than 1σ above the CO:SF correlation, and less than 1σ below the PAH:SF correlation.

From Figure 4 (bottom left) we find the PAH:SF and H₂ S(3):SF correlations have the highest residual scatter. More than 1σ below these are the IR continuum bands: WISE 3.4–12 μm, PACS 100, 160 μm, and SPIRE 250 μm. About 1σ below these are the remaining H₂ rotational lines, fine structure lines, and the 24 and 70 μm IR bands. The [S iii] lines or their sum have the lowest scatter of any SF correlation.

From Figure 4 (bottom right) we find most IR features have a slope with SFR that is consistent with linear to within 1σ. The exceptions are all sublinear (m < 1), the lowest being m ∼ 0.5 for the H₂ rotational lines and Humphreys-α, but we note the latter is a weakly detected line with few detections as a function of Σ_SFR.

To quantify the difference between the SF- and CO-tracing features, we plot the difference between the correlation coefficients ρ_SFR−ρ_CO from Figure 4 as a function of the observed wavelength in Figure 6.

Zoom In Zoom Out Reset image size

Figure 6 (left) shows the difference in correlations ρ_SFR−ρ_CO of each photometric band as a function of wavelength from WISE1 3.4 μm to SPIRE 250 μm. We find the emission from ∼20−70 μm (i.e., WISE 22 μm, MIPS 24 μm, and PACS 70 μm) is SF dominant. The WISE3 12 μm and PACS 100 μm bands are found to be slightly CO dominant, but not by greater than the 1σ statistical uncertainty. All other bands, from WISE1 3.4 μm to IRAC4 8 μm, PACS 160 μm, and SPIRE 250 μm are CO dominant by greater than 1σ. The overall trend in this figure suggests a potential peak at about 40 μm, but we only have evidence that the continuum between 20 and 70 μm is likely a comparably strong SF tracer as the 24 and 70 μm emission.

We also use our synthetic photometry decomposition (see Section 2.8) to show the difference between SF and CO correlations for the continuum and PAH portions of the IRAC4 8 μm and WISE3 12 μm bands in Figure 6 (left). Notably, the continuum portion of the WISE band is SF dominant while the PAH portion is CO dominant, in good agreement with our expectations based on the behavior of the nearby 24 μm continuum and that of each individual PAH band in Figure 6. The distinction is less pronounced in the decomposed IRAC4 band. The continuum portion of both the IRAC and WISE bands better fits the general trend over wavelength seen in Figure 6 (left) than the total band with PAH emission included. This shows a clear difference in behavior between the PAH vibrational emission and the small, hot dust grain continuum in the MIR, which we discuss further in Section 4.

From Figure 6 (right) we find the emission features from MIR spectra that are most CO dominant are the molecular gas rotational lines H₂ S(0) through H₂ S(3). The correlation coefficients of PAH features are remarkably similar to the H₂ rotational lines and to each other; they are less CO dominant by only about 1σ due to their stronger correlation with SF compared to the H₂ rotational lines. However, the residual scatter and power-law slopes vary significantly between PAH and H₂ rotational line correlations (see Figure 4, bottom; discussed further in Section 4.1.1).

We find the fine structure lines from ions with ionization potential greater than hydrogen—[Ne ii], [Ne iii], [S iii], and [S iv]—are all SF dominant by at least 1σ, with exception of the blend of the H₂ S(5) 6.91 μm and [Ar ii] 6.99 μm lines, which traces SF and CO equally. Since we found the other H₂ lines are highly CO dominant, this blended line likely differs from the others since it includes a mix of CO-dominant H₂ emission in addition to the SF-dominant fine structure line emission from Ar⁺. We also see the Humphreys-α line from ionized hydrogen is highly SF dominant, but this line is weakly detected and has large uncertainties. The remaining fine structure line from Si⁺ has a lower ionization potential than hydrogen, and therefore may originate both in the ionized and neutral gas, but we find it is SF dominant by ∼1σ.

We find both of the lines from S²⁺ are highly SF dominant and behave indistinguishably from the reference SF tracer itself in Figures 4 and 5. This suggests either of the [S iii] lines can be used as a single-parameter MIR tracer of star formation (see Section 4.3.1). The correlation is equally strong for each [S iii] line individually as well as all combinations of summing them with the [S iv] line, as shown in Figure 7. For all individual and summed lines among [S iv], [S iii] 18.7 μm, and [S iii] 33.5 μm, the slope of their correlation with SF is about 1σ greater than linear (m ∼ 1.1, see Table 4). There is a small trend in the scatter of these correlations with the ratio of 15.6 μm [Ne iii] to 12.8 μm [Ne ii]. However, including this ratio, or its physical correlate, metallicity (Madden et al. 2006; Whitcomb et al. 2020), as an additional feature in the regression results in a negligible improvement to the already tight correlation with Σ_SFR, slightly reducing the scatter σ.

Zoom In Zoom Out Reset image size

This behavior is distinct from the similar pair of neon MIR fine structure lines, 12.8 μm [Ne ii] and 15.6 μm [Ne iii]. The neon lines individually as a function of Σ_SFR have a strong metallicity dependence in the scatter, which is partially corrected by summing the two lines (see Zhuang et al. 2019; Whitcomb et al. 2020).

Using the fit constants from Table 4 in Equation (3), the single-line tracer of Σ_SFR based on the 18.7 μm [S iii] line is given by

$$\begin{eqnarray}\begin{array}{rcl}\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{SFR}}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-2}}\right)\, & = & \,{10}^{0.74\,\pm \,0.11}\\ & & \times \,{\left(\displaystyle \frac{[{\rm{S}}\,{\rm\small{III}}]18.7\,\mu {\rm{m}}}{{10}^{-7}\,{\rm{W}}\,{{\rm{m}}}^{-2}\,{\mathrm{sr}}^{-1}}\right)}^{1.13\,\pm \,0.08}.\end{array}\end{eqnarray} \tag{ 4 }$$

The most CO-dominant tracers of Figure 6 (right) are found to be the rotational lines of H₂ followed closely by the PAH features. We find a strong mutual correlation between these three types of emission (Spearman coefficient >0.8); each H₂ line intensity is well correlated with both PAH and CO emission intensities, and the brightest PAH features are just as well correlated with both H₂ and CO emission. The correlation between PAH and H₂ emission has been well established in the literature for regions with negligible nonradiative heating (i.e., heating not driven by UV photons) such as the star-forming regions in this study (Roussel et al. 2007; Naslim et al. 2015; Smercina et al. 2018). In regions with significant nonradiative heating (e.g., shocks, turbulent dissipation), studies have found enhanced H₂ rotational relative to PAH emission (Ingalls et al. 2011). Spatial and intensity correlations between PAH and CO emission have also been observed in many previous and recent studies (Regan et al. 2006; Cortzen et al. 2019; Chown et al. 2021; Leroy et al. 2023).

In Section 4.1.1 we discuss the PAH and H₂ correlations with Σ_SFR and CO in detail. The PAH and H₂ correlations with CO have about the same amount of residual scatter (except PAH 17 μm versus CO which has higher scatter). The PAH features, however, have more scatter than H₂ rotational lines do versus SF (except H₂ S(3)). The scatter in the PAH:SFR relation is dependent on metallicity, but metallicity does not introduce scatter in the PAH:CO relation (see center panels of Figure 3). All the H₂ lines have a sublinear slope as a function of both CO and SF, while the PAHs are approximately linear as a function of SF and the slope as a function of CO varies from sublinear for shorter wavelength features to superlinear for longer wavelength features.

The strong correlations between PAHs, CO, and H₂ rotational lines agree with expectations for emission generated by dense PDRs. In typical star-forming regions PAH vibrational bands are excited by UV photons, CO rotational lines are collisionally excited, and the H₂ rotational lines are excited by both collisions and UV pumping (Tielens 2008; Pereira-Santaella et al. 2014). The collisional excitation rates are set by the density and the gas temperature, which is the result of heating from the photoelectric effect. Photoelectric (PE) heating from PAHs and dust is the dominant heating mechanism in the warm portion of dense PDRs and it continues to dominate at optical depths where CO emits most strongly (Wolfire et al. 1993, 2022). The PE efficiency of PAHs is found to depend on the incident UV field (Croxall et al. 2012) and through the resulting photoelectrons the gas is heated indirectly by recent SF. Thus, the PAH, CO, and H₂ rotational line emission are all tracing the conversion of UV photons into excitation and collisions, which may explain the strong correlations in dense PDRs both spatially and in intensity.

We found PAHs and H₂ lines have a significantly weaker correlation with SF compared to the fine structure lines from ions. Comparing the correlations with SF for CO, H₂, and PAHs, we see H₂ is slightly less correlated with SF than PAHs. This agrees with a picture where the H₂ excitation is mainly collisional whereas PAHs respond directly to the local UV field. The warm portion of a dense PDR (10^2-3K) is ideal for excitation of the H₂ rotational lines. Closer to the H ii regions the UV flux is sufficient to dissociate H₂(Roussel et al. 2007; Togi & Smith 2016). Togi & Smith (2016) conclude that the warmer portion of gas traced by the MIR H₂ rotational lines is typically about 15% of the total molecular gas mass. The strong correlations found in our results at small and large scales between each of the H₂ rotational integrated line intensities and CO J = (2 − 1) surface brightness are evidence of both a spatial connection and a connection in how the warm and cold molecular gas emission lines are excited. The envelope of warm gas surrounding cold molecular clouds is bright in H₂ emission, but some of it is CO dark (Wolfire et al. 2010; Leroy et al. 2011; Smith et al. 2014; Lee et al. 2015; Madden et al. 2020). In addition, we found the H₂ and CO emission have a similar correlation strength with SF (with H₂ stronger than CO by about 1σ). This may be evidence that the H₂ rotational emission is tracing a fairly constant fraction of the total H₂, consistent with the findings of Roussel et al. (2007) and Togi & Smith (2016) for the SINGS sample.

The PAH emission features require UV excitation, and we find they have a better correlation with Σ_SFR compared to H₂ rotational lines. PAHs are relatively large molecules and/or small dust grains, and it is theorized that high-energy photons may destroy them (Micelotta et al. 2010b; Chastenet et al. 2019). Therefore, we do not expect PAH emission to be spatially correlated with ionized gas on the smallest scales of individual H ii regions (≲10 pc). The PAH features are CO dominant, but the correlation PAH:SF is slightly stronger than H₂:SF (though by less than the 1σuncertainty) while PAH:CO is the same as H₂:CO. Unlike the low-J H₂ rotational lines, PAH vibrational excitation is more dependent on radiation hardness (Draine et al. 2021), which may account for the small difference in SF correlation for PAH and H₂ emission.

Since PAH and H₂ rotational emission are both spatially correlated with molecular gas but only PAH emission is directly dependent on the UV intensity, we might expect the PAH:SF correlation to be significantly stronger than H₂:SF, but this is not what we observe. The large scatter in the PAH:SF correlation has a strong dependence on the ratio [Ne iii]/[Ne ii], but the scatter in H₂:SF has a negligible dependence on this ratio. This suggests PAH emission is more sensitive to changes in metallicity than H₂ emission, which has been established in many previous studies (Engelbracht et al. 2005; Madden et al. 2006; Gordon et al. 2008; Muñoz-Mateos et al. 2009; Sandstrom et al. 2012; Lai et al. 2020; Zang et al. 2022). The ratio [Ne iii]/[Ne ii] tracks changes in radiation hardness and metallicity, so the scatter in the PAH:SF correlation is likely driven by variations in these properties. PAHs respond directly to UV photons, but the resulting emission bands have a dependence on the abundance of PAHs and their size and charge distributions, which in turn are dependent on local environmental properties such as metallicity or radiation hardness (Draine et al. 2021). Our sample is an ensemble of star-forming regions from many different galaxies, so these dependencies on local conditions are likely contributing to the scatter that weakens the correlation between PAH emission and Σ_SFR.

4.1.1. Connections between Individual PAH and H₂ Correlations with CO and Σ_SFR

In Figure 8 (top left) we show our power-law fit with CO for all four studied H₂ emission lines. We find the power-law slope of each correlation is approximately the same (about 0.7) but the intercept varies between the lines, in order from dimmest: H₂ S(0), (2), (3), and S(1) is the brightest of the four. The trend between these warm molecular gas tracers and CO is found to be nonlinear and the brightest of the warm gas lines S(1) has the strongest correlation with CO, and with all the major PAH features as well. Figure 8 (top right) shows these same correlations with Σ_SFR for comparison.

Zoom In Zoom Out Reset image size

The lack of variation in the slope of each correlation in Figure 8 (top left) suggests that the heating contribution by UV pumping does not change between the lines. If UV pumping was significant, we would expect the slope as a function of Σ_SFR for higher-energy rotational lines, such as H₂ S(3), to be greater than the slope of the ground state S(0) line. Instead, we find the slopes are approximately equal for all four rotational lines (∼0.55 ± 0.1). Since the correlation with CO for each line is also much stronger, we conclude direct excitation of H₂ via UV pumping is unlikely to be dominant compared to collisional excitation for the H₂ rotational lines in our sample. This agrees with the analysis of the SINGS sample by Roussel et al. (2007).

We also include the correlations for the four brightest PAH features in Figure 8 (bottom). In contrast with the H₂:CO correlations, the power-law slope increases for CO correlations with longer wavelength PAH features, ranging from 0.8 at 7.7 μm, 0.9 at 11.3 μm, and 1.2 at 17 μm (see Table 4). In general, the PAH features at longer wavelengths originate more from larger PAHs, while shorter wavelength PAH features are more effectively emitted by small PAHs (Draine & Li 2007; Maragkoudakis et al. 2020; Draine et al. 2021). The longer wavelength features are brighter relative to the 7.7 μm feature at higher CO brightness, which may suggest PAH growth and/or more neutral PAHs in molecular gas regions. PAH growth in molecular clouds has been proposed by previous studies in the literature (Sandstrom et al. 2010; Chastenet et al. 2019). The increase in the 17/7.7 μm ratio could also result from a decrease in the average hardness of the radiation field as a function of CO. PAHs exposed to harder radiation fields, i.e., with more UV relative to lower-energy photons, emit significantly more in shorter relative to longer wavelength features (Draine et al. 2021). If dust content is higher in regions with more CO, more of the photons will be attenuated which could decrease the average hardness of the radiation field that is exciting the PAHs, potentially accounting for the enhanced 17/7.7 μm PAH ratio as a function of CO emission.

In Figure 8 (bottom right) we find the slope of the correlation with Σ_SFR is the same within 1σ for each major PAH feature. This implies that band strength variations are not driven by changes in the incident UV field intensity. The fact that the slope as a function of Σ_SFR does not vary between PAH features could be explained if the radiation field hardness is not correlated with Σ_SFR, or if the PAH population is changing in such a way that negates the expected band variations. If it is the PAH population changing, then the size and charge distribution must be shifting to larger and more neutral PAHs as a function of Σ_SFR. This is possible if small PAHs are destroyed or coagulated into larger PAHs, as the PAH:CO slopes could be interpreted to suggest.

Throughout this work, we have studied correlations with CO J = (2 − 1) emission. The issues with converting to molecular gas mass with a CO-to-H₂ conversion factor are explored in-depth in Bolatto et al. (2013). The CO-to-H₂ conversion factor is generally used with CO J = (1 − 0) emission, but our results are based on correlations with CO J = (2 − 1) emission. CO J = (2 − 1) can be converted to CO J = (1 − 0) with a ratio R₂₁ that is ∼0.65 on average and can vary between ∼0.4 and 1.0 (Leroy et al. 2022). This introduces some uncertainty in translating our correlations with CO J = (2 − 1) emission into correlations with M_H2. However, R₂₁ is correlated with Σ_SFR, so we expect the values for our sample to be relatively constant and near the high end of this range (R₂₁ ∼ 1.0). The CO-to-H₂ conversion factor also has a strong dependence on metallicity (Wolfire et al. 2010; Leroy et al. 2011; Glover & Clark 2012) and on CO excitation and optical depth (Israel 2020; Teng et al.2022). Our sample of H ii regions in nearby galaxies has a relatively small range of metallicities, with the exception of a few dwarf galaxies.

In the centers of a few of our targets there are also large deviations in CO-to-H₂ conversion factor due to excitation and optical depth effects (Sandstrom et al. 2013; Israel 2020). In Figure 3 we circled in red three such regions for which the CO-to-H₂ conversion factor is significantly different from the Milky Way value; these regions are in NGC 3184, NGC 3351, and NGC 6946, which have α_CO lower than the Milky Way value by 0.39, 0.79, and 1.04 dex, respectively (determined using dust as a tracer; Sandstrom et al. 2013). Nevertheless, the three points are not outliers in the H₂:CO correlation despite their different conversion factor compared to the other points. If the measured CO-to-H₂ conversion factors were applied, we would expect the red points to shift to lower H₂ mass than the others of similar CO brightness, likely increasing the scatter compared to the correlation with CO. The fact that these two points are not significant outliers in the trend between CO and H₂ rotational lines may suggest that the "H₂ rotational-to-H₂" conversion factor differs in those regions in the same way that CO-to-H₂ does. This could result from changes in the temperature distribution of the H₂ (either a change in the warm H₂ fraction or a change in the temperature distribution) coupled with changes in the temperature distribution and/or optical depth of the CO. If the α_CO changes are due to higher gas temperature, we might expect to observe warmer H₂ rotational temperatures in these regions, traced by H₂ S(3)/H₂ S(1). Togi & Smith (2016) also find that changes in the warm H₂ fraction and H₂ temperature may be correlated. With the three measurements currently available, we do not see clear differences in H₂ S(3)/H₂ S(1). Future studies of H₂ rotational ladders with JWST-MIRI in regions where α_CO varies may provide constraints on what processes can alter the CO-to-H₂ conversion factor.

The only SF-dominant tracers of Figure 6 (right) in MIR spectra are the fine structure lines from ions. Most of these have ionization potentials greater than that of hydrogen. The tight correlation of these ionized gas lines and our ionizing-photon-based reference SFR tracer is expected. The lack of correlation with CO emission also implies a spatial decorrelation between lines from the most highly ionized elements and cold gas. This is especially clear in Figure 4 for the case of 15.6 μm [Ne iii]. Despite being a component of our reference Σ_SFR tracer, we find the [Ne iii]:CO correlation is significantly weaker than the K-S relation. Doubly ionizing neon requires photon energies of at least 41 eV, so the lack of correlation with cold molecular gas and strong correlation with massive SF is expected for such species.

This correlation is also exaggerated by the combination of spatial and intensity dependencies: the [Ne iii] emission is not only spatially uncorrelated with CO in our smallest regions, but their intensities are driven by uncorrelated mechanisms. The [Ne iii] line strength depends directly on the spectrum of ionizing radiation >41 eV, the electron density, and the abundance of Ne²⁺ ions relative to hydrogen. The CO line strength depends on many factors, only one of which is an indirect dependence on the UV field intensity through PE heating of the gas by PAHs and dust.

This spatial separation between ionized and molecular gas is feasible since the median physical area of our apertures is about 0.4 kpc², corresponding to a scale length of about 600 pc. Schinnerer et al. (2019) found that in maps of nearby spiral galaxies with resolution smoothed to this scale, the molecular and ionized regions (as traced by CO and Hα emission, respectively) overlap in about 50% of pixels across the galaxy. Since the SINGS apertures were placed based on peaks of Hα emission (Kennicutt & Armus 2003), we expect that in at least half of each region with area <0.4 kpc² the Σ_SFR is resolved distinctly from the CO and the two are spatially uncorrelated.

In Figure 9 we show the difference between SF and CO correlation coefficients ρ_SFR–ρ_CO for the MIR fine structure lines as a function of their ionization potential. From this, we find a trend of increasing SF dominance with increasing ionization potential. All lines are SF dominant by at least 1σ, except the blended H₂ S(5) and [Ar ii] 6.9 μm line. Si⁺ has a lower ionization energy than hydrogen so it can exist and emit in ionized and neutral regions. This implies [Si ii] emission is SF dominant because the intensity is correlated and some of the emission may be spatially coincident with the H ii region.

Zoom In Zoom Out Reset image size

4.3.1. MIR [S iii] Lines as a Single-feature SF Tracer

4.3.2. Trends in IR Continuum Correlations with SF and CO

All PAH features are CO dominant but still correlated with SF (ρ_SFR ∼0.65), so PAH emission is associated with cold gas emission and associated with recent SF through the need for UV excitation. But there is also a significant amount of UV and optical radiation that can excite PAHs from the more numerous stars that are not recently formed O stars (Bendo et al. 2008; Crocker et al. 2013). In regions where the diffuse radiation field is stronger than those from SF, the excitation of the PAH vibrational emission will not be due to recent SF (Draine et al. 2007; Ingalls et al. 2011). In this case, the majority of PAH excitation will be due to stars less massive than O stars, implying the PAHs are not primarily tracing the most recent generation of SF but the SF averaged over long timescales.

Our investigation into the correlations between SF, CO, and the various portions of IR continuum, summarized in Figure 6 (left) has implications for the use of FIR wavelengths as SF tracers. We found that dust continuum beyond about 100 μm better traces CO emission than SF. Since dust content increases along with cold gas content, this result indicates the FIR emission is more sensitive to changes in the dust content than to the prevalence of recent SF events.

The H ii regions in this study are in galaxies within ∼20 Mpc. Galaxies at higher redshift are lower metallicity on average, so we hypothesize what the MIR correlations may look like for these high-redshift sources. Since we have few low-metallicity regions (only 12 between $\tfrac{1}{5}$ and $\tfrac{4}{5}$ Z_☉), we only make inferences based on the presence or absence of a metallicity trend in the scatter of each correlation.

The [S iii] 18.7 and 33.5 μm lines have no significant metallicity dependence in the scatter of their correlation with Σ_SFR. The SF tracer based on either of these lines is insensitive to changes in metallicity in our sample. Other MIR tracers of SF, such as PAH or 15.6 μm [Ne iii] and 12.8 μm [Ne ii] emission, require explicit metallicity corrections to reproduce Σ_SFR values as accurate as those derived from [S iii] emission alone.

The scatter in PAH emission correlations with CO is negligible, but the scatter in PAH:SF correlations has a significant metallicity dependence. Since the PAH:SF correlation fails at the low metallicities of our data set but the PAH:CO correlation does not, it is likely that PAH emission in high-redshift galaxies better traces CO than it traces SF.

The case is similar for H₂ rotational lines, but the H₂:SF scatter is less significant than in PAH:SF. This again implies that in high-redshift galaxies, the H₂ rotational lines will be better tracers of CO than SF, but the distinction will be less pronounced than for PAH emission at high redshift.

Table 5 shows the general form and content of the machine-readable table of data used in this work. This includes the coordinates of the vertices of the 154 apertures that define the overlap between SINGS SL, SH, and LH coverage, as well as all data derived for each: the aperture angular and physical area, the SFR surface density, integrated intensities of CO J = (2 − 1) and all IR photometry and spectroscopic features described in Section 2.

Figure 10 shows a comparison between PAHFIT and FARAMIR integrated intensities (see Section 2.7) for four emission features in the overlapping wavelengths of the SL and SH suborders of Spitzer-IRS (∼10 to 15 μm). PAHFIT is able to fit dimmer PAH features than FARAMIR, but the agreement is strong at higher intensities.

Zoom In Zoom Out Reset image size