PHANGS–JWST First Results: A Global and Moderately Resolved View of Mid-infrared and CO Line Emission from Galaxies at the Start of the JWST Era

We compile a large set of galaxy-integrated estimates of CO luminosity, mid-IR luminosity, effective radius, stellar mass, and SFR. Given the sky-wide coverage of WISE, the limiting factor is the availability of CO data, so we construct the sample around CO surveys. We compile galaxy-integrated CO (1–0) luminosities from ALMaQUEST (Lin et al. 2020), AMIGA (Lisenfeld et al. 2011), CARMA-EDGE (Bolatto et al. 2017), MASCOT (Wylezalek et al. 2022), xCOLD GASS (Saintonge et al. 2017), and galaxy-integrated CO (2–1) luminosities from APEX-EDGE (Colombo et al. 2020), ALLSMOG (Bothwell et al. 2014; Cicone et al. 2017), and a survey by Jiang et al. (2015). We also use the galaxy-integrated values from the CO (2–1) and CO (1–0) mapping samples COMING (Sorai et al. 2019) and HERACLES (Leroy et al. 2009), with the additional IRAM 30 m maps described in Leroy et al. (2022), PHANGS–ALMA (Leroy et al. 2021b), and the Nobeyama Atlas by Kuno et al. (2007). The calculation of these integrated values for mapping surveys are described in Leroy et al. (2022) and when necessary includes an aperture correction to convert the mapped luminosity to the full galaxy CO luminosity (see, e.g., Leroy et al. 2021b).

We assign mid-IR luminosities at 12 and 22 μm, stellar mass (M_⋆), and SFR estimates by cross-matching each galaxy with a CO measurement through either the GSWLC database (Salim et al. 2016, 2018) or the z0MGS atlas (Leroy et al. 2019). In the handful of cases where we could not successfully cross-match a CO survey target to either of these mid-IR catalogs, we dropped the target from our analysis.

GSWLC compiled WISE catalog photometry and spectral-energy-distribution (SED)-modeling-based stellar masses and SFR estimates for more distant galaxies, while z0MGS constructed images, carried out photometry, and estimated SFR and stellar mass from UV, near-IR, and mid-IR for extended local systems. The z0MGS methods to estimate SFR and stellar mass were designed to yield values that match the GSWLC values, so all of the integrated values in this Letter should be self-consistent and anchored to the results of Salim et al. (2016, 2018), which use CIGALE (Boquien et al. 2019) modeling and the Bruzual & Charlot (2003) models.

We also assign an effective radius, R_e , to each galaxy. For galaxies at distance <50 Mpc, we prefer near-IR-based effective radius measurements by Muñoz-Mateos et al. (2015) or Leroy et al. (2021b). When these are not available, we estimate the effective radius from the optical R₂₅ in the LEDA database (Paturel et al. 2003; Makarov et al. 2014) using a scaling relation in which the ratio R_e /R₂₅ depends on stellar mass (derived from the data in Muñoz-Mateos et al. 2015). For more distant galaxies, we draw the effective radius from the survey itself. For the Sloan Digital Sky Survey (SDSS)-based surveys, these tend to be based on the r band, while for CALIFA these are based on stellar mass estimates from the modeling of integral field unit spectroscopy (Sánchez et al. 2016). Taking effective radius measurements from different literature sources adds a small systematic uncertainty that does not have a large impact on our results.

We use these data to estimate the mean I_CO, I_{12 μm}, and I_{22 μm} within R_e for each target. To do this, we assume that half of the luminosity in CO and mid-IR emission emerges from within the stellar effective radius, R_e , and calculate ${I}_{x}=0.5\,{L}_{x}/\pi {R}_{e}^{2}$ for each band x. Such an assumption is well justified, on average, for star-forming galaxies where the CO, stellar, and star formation distributions tend to track one another quite well (e.g., Young & Scoville 1991; Young et al. 1995; Regan et al. 2006; Leroy et al. 2008, 2009; Bolatto et al. 2017; Lin et al. 2020; Brown et al. 2021; though, see especially the latter for some important subtleties).

For the integrated galaxies, we treat observations with a signal-to-noise ratio (S/N) < 5 as nondetections and record their 5σ value as an upper limit. In practice, this affects CO much more than the mid-IR. For this part of the analysis, we also allow repeat measurements in the case that different surveys target the same galaxy.

We also assemble a large set of 17'' resolution measurements of CO emission and mid-IR emission that resolve nearby galaxies into individual regions. We target 17'' because this represents the common resolution achievable by WISE 22 μm data (θ ≈ 15''). ⁴¹ and much of the available single-dish CO mapping (θ ≈ 15''–17''). We smooth to a common Gaussian beam with a FWHM of 17'', which allows us to combine essentially all recent large CO (1–0) and CO (2–1) mapping surveys, Spitzer Infrared Array Camera (IRAC) 8 μm and Multiband Imaging Photometer (MIPS) 24 μm mapping, and WISE 12 and 22 μm mapping.

We use CO maps from BIMA SONG (CO (1–0); Helfer et al. 2003), COMING (CO (1–0); Sorai et al. 2019), HERACLES (CO (2–1); Leroy et al. 2009), and the follow up IRAM observations described in Leroy et al. (2022), the Nobeyama survey (CO (1–0); Kuno et al. 2007), PHANGS–ALMA (CO (2–1); Leroy et al. 2021b), and VERTICO (CO (2–1); Brown et al. 2021). Wherever possible we use the combined interferometric and total power data: all sources in BIMA SONG, PHANGS–ALMA, and many in VERTICO include total power (see table in Brown et al. 2021). These surveys sometimes target the same galaxy, and for this resolved part of the analysis we only use a single map per galaxy per transition, with our preferences set by the combination of calibration quality, field of view, and sensitivity. For CO (1–0) our priority is the COMING survey over the survey by Kuno et al. (2007) over BIMA SONG. For CO (2–1) we prefer PHANGS–ALMA over the IRAM/HERACLES survey over VERTICO.

We combine these CO maps with WISE 12 and 22 μm images from Leroy et al. (2019), based on Lang (2014). Whenever available, we also include IRAC 8 μm emission from either the Local Volume Legacy (LVL; Dale et al. 2009) or the SIRTF Nearby Galaxies Survey (SINGS; Kennicutt et al. 2003) and MIPS 24 μm emission from LVL, SINGS, or the data compilation by Bendo et al. (2012b).

We convolve all data to a common Gaussian resolution of 17'' at FWHM, using the kernels by Aniano et al. (2011) when necessary to convert from the point-spread function of the instrument to a Gaussian. At the distances to our targets, this 17'' resolution corresponds to a 16%–84% range of 0.5–1.9 kpc, with a median of 1.3 kpc, for the sample of azimuthal averages and a 16%–84% range of 0.3–1.8, with a median of 1.2 kpc, for the individual regions.

For the CO data, we mostly follow the procedures described in Leroy et al. (2021a, 2022) to convert from data cubes into maps with associated uncertainties. The only important difference is that for this work we create moment maps by integrating the line over a single velocity window for each galaxy tailored to encompass the emission. That is, we create simple "flat" moment maps that enable robust averaging. For the mid-IR maps, we estimate the uncertainties in the convolved data from empty regions of the sky near the galaxy. We project all data for each galaxy onto a common astrometric grid with 85 pixels, roughly Nyquist sampling the 17'' beam.

We exclude galaxies with inclination >70° from the resolved analysis (no inclination cut is applied in the integrated analysis). All resolved intensities are corrected for the effects of inclination by applying a factor of $1/\cos i$, with the inclinations drawn from Lang et al. (2020) when available or otherwise based on isophotal axial ratios used in Muñoz-Mateos et al. (2015) or from LEDA (Paturel et al. 2003).

We then use the convolved, inclination-corrected data to construct azimuthally averaged intensity profiles for each target. We average intensities in elliptical annuli of width equal to the FWHM beam size and spaced by one half-beamwidth, so that we oversample the beam radially by a factor of 2. The orientations for these annuli are drawn from the same sources as the inclination. The averages exclude data within ±30° of the minor axis to avoid the most severe projection oversampling effects.

We propagate errors from the original maps and mark rings and pixels with S/N > 5 as detections (as mentioned below, the S/N in CO is the limiting criterion for the overwhelming majority of the data). We retain information on the cases with S/N < 5 so that we can still use these to average and determine integrated trends, but consider these individual measurements unreliable.

We have compiled a large fraction of the local universe CO (1–0) and CO (2–1) observations of massive, star-forming galaxies. A detailed analysis of the effects of combining these different surveys is beyond the scope of this Letter, but Tables 1 and 2 and Figure 1 provide an overview of the resulting data set. The figures and Table 2 show the stellar mass and SFR for CO-detected galaxies and compare these to the "main sequence," relating SFR to M_⋆ for low-redshift star-forming galaxies (here we use the broken power-law form from Leja et al. 2022).

Our integrated-galaxy CO detections span from about ${\mathrm{log}}_{10}{M}_{\star }\sim \,$ 9.5–11 and mostly lie near the star-forming main sequence. Although there are significantly more galaxy-integrated CO (1–0) measurements than CO (2–1) measurements, both data sets do a good job of spanning this range. The CO (2–1) mapping covers the same range and also samples the star-forming main sequence well. ⁴² The CO (1–0) mapping data cover a narrower range of masses and show some bias toward high SFRs at fixed M_⋆. We discuss this in more detail and correct our results for this effect in Section 3.4.

We focus on massive, star-forming galaxies in our compilation. Although the data set includes a few early-type and low-SFR/M_⋆ galaxies, many of these are CO nondetections. We do not include CO mapping or galaxy-integrated surveys exclusively focused on early-type galaxies (e.g., Combes et al. 2007; Crocker et al. 2011; Young et al. 2011; Alatalo et al. 2013; Davis et al. 2014). Nor do we probe far into the dwarf galaxy regime. For example, we do not include measurements by Schruba et al. (2012), Hunt et al. (2015) or similar work. Finally, we also do not include dedicated surveys of major mergers or luminous/ultraluminous infrared galaxies (e.g., Herrero-Illana et al. 2019; Lisenfeld et al. 2019). Extending this work to such targets will be interesting in future work.

Our goal is to characterize the observed relationship between each CO line and each mid-IR band for integrated and moderately resolved galaxies. We aim to measure the strength of the correlation, typical ratio and scatter, and to provide basic functional relationships linking the two.

To do this, we construct a series of variable pairs, separately combining CO (1–0) and CO (2–1) intensity with each mid-IR band (see Tables 1 and 3). For integrated galaxies, we focus on WISE 12 and 22 μm, which are available for the whole sample. ⁴³ For these integrated galaxies, we also pair the CO-to-mid-IR ratio with the stellar mass, M_⋆, and with the specific star formation rate, SFR/M_⋆, estimated as described in Section 2.1. For resolved galaxies, we also pair each CO line with the available Spitzer 8 and 24 μm data.

Table 3. Summary of CO vs. Mid-IR Scaling Relations

Data Pair	Rank	Ratio	Scatter	Bin Range	m	b	σresid	Figures
(y versus x)	corr.	(${\mathrm{log}}_{10}y/x$)	(dex)	(in ${\mathrm{log}}_{10}x$)	(Equations (3) and (5))		(dex)
Integrated galaxies
CO (1–0) versus 12 μm	0.87	0.20	0.22	−0.5 to 1.9	0.86 ± 0.02	0.26 ± 0.02	0.20	2
CO (1–0) versus 22 μm	0.80	0.01	0.31	−0.3 to 2.1	0.67 ± 0.05	0.23 ± 0.05	0.25	2
CO (1–0)/12 μm versus M⋆ a	0.14	⋯	⋯	9.0 to 11.25	0.13 ± 0.05	0.14 ± 0.02	0.22	9
CO (1–0)/22 μm versus M⋆ a	0.23	⋯	⋯	9.25 to 11.25	0.19 ± 0.01	−0.06 ± 0.01	0.28	9
CO (1–0)/12 μm versus SFR/M⋆ a , b	−0.38	⋯	⋯	−11.5 to −9.25	−0.14 ± 0.05	0.14 ± 0.03	0.20	10
CO (1–0)/22 μm versus SFR/M⋆ a , b	−0.56	⋯	⋯	−11.5 to −9.25	−0.34 ± 0.11	−0.13 ± 0.07	0.23	10
CO (2–1) versus 12 μm	0.92	0.11	0.23	−0.5 to 1.5	1.02 ± 0.03	0.11 ± 0.02	0.21	2
CO (2–1) versus 22 μm	0.87	−0.09	0.28	−0.3 to 1.9	0.85 ± 0.05	0.01 ± 0.03	0.28	2
CO (2–1)/12 μm versus M⋆ a	0.45	⋯	⋯	9.0 to 11.0	0.16 ± 0.07	0.07 ± 0.03	0.17	9
CO (2–1)/22 μm versus M⋆ a	0.47	⋯	⋯	9.0 to 11.0	0.22 ± 0.04	−0.15 ± 0.02	0.24	9
CO (2–1)/12 μm versus SFR/M⋆ a , b	−0.34	⋯	⋯	−11.0 to −9.25	−0.14 ± 0.12	0.09 ± 0.04	0.19	10
CO (2–1)/22 μm versus SFR/M⋆ a , b	−0.57	⋯	⋯	−11.25 to −9.25	−0.41 ± 0.09	−0.16 ± 0.04	0.23	10
Rings in radial profiles
CO (1–0) versus 12 μm	0.90	0.38	0.29	−1.0 to 1.75	0.83 ± 0.02	0.40 ± 0.02	0.17	3
... with sample correction c	⋯	⋯	⋯	⋯	0.83	0.21	⋯
CO (1–0) versus 22 μm	0.86	0.17	0.39	−1.0 to 2.25	0.65 ± 0.05	0.30 ± 0.03	0.21	3
... with sample correction c	⋯	⋯	⋯	⋯	0.65	0.11	⋯
CO (1–0) versus 8 μm	0.91	0.23	0.19	−1.0 to 1.75	0.88 ± 0.05	0.29 ± 0.04	0.16	4
... with sample correction c	⋯	⋯	⋯	⋯	0.88	0.10	⋯
CO (1–0) versus 24 μm	0.90	0.28	0.39	−1.0 to 2.0	0.71 ± 0.06	0.33 ± 0.03	0.19	4
... with sample correction c	⋯	⋯	⋯	⋯	0.71	0.14	⋯
CO (2–1) versus 12 μm	0.94	0.13	0.22	−1.0 to 2.0	1.04 ± 0.05	0.09 ± 0.04	0.19	3
CO (2–1) versus 22 μm	0.91	−0.03	0.29	−1.0 to 2.25	0.91 ± 0.05	−0.09 ± 0.04	0.23	3
CO (2–1) versus 8 μm	0.94	−0.05	0.26	−1.0 to 1.75	1.14 ± 0.03	−0.11 ± 0.02	0.21	4
CO (2–1) versus 24 μm	0.91	0.01	0.36	−1.0 to 2.0	0.95 ± 0.07	−0.04 ± 0.02	0.23	4
Individual 17'' regions
CO (1–0) versus 12 μm	0.71	⋯	⋯	−1.0 to 2.25	0.79 ± 0.03	0.41 ± 0.03	0.19	5
... with sample correction c	⋯	⋯	⋯	⋯	0.79	0.22	⋯
CO (1–0) versus 22 μm	0.68	⋯	⋯	−0.5 to 2.75	0.59 ± 0.05	0.34 ± 0.07	0.20	5
... with sample correction c	⋯	⋯	⋯	⋯	0.59	0.15	⋯
CO (1–0) versus 8 μm	0.63	⋯	⋯	−0.75 to 1.5	0.89 ± 0.06	0.21 ± 0.03	0.22	6
... with sample correction c	⋯	⋯	⋯	⋯	0.89	0.02	⋯
CO (1–0) versus 24 μm	0.70	⋯	⋯	−0.75 to 1.5	0.73 ± 0.05	0.37 ± 0.03	0.21	6
... with sample correction c	⋯	⋯	⋯	⋯	0.73	0.18	⋯
CO (2–1) versus 12 μm	0.87	⋯	⋯	−1.0 to 2.25	0.96 ± 0.05	0.04 ± 0.04	0.18	5
CO (2–1) versus 22 μm	0.82	⋯	⋯	−0.5 to 2.25	0.80 ± 0.04	−0.01 ± 0.04	0.20	5
CO (2–1) versus 8 μm	0.81	⋯	⋯	−0.75 to 1.5	1.17 ± 0.08	−0.19 ± 0.05	0.21	6
CO (2–1) versus 24 μm	0.83	⋯	⋯	−0.75 to 1.75	0.90 ± 0.06	0.03 ± 0.04	0.19	6

Notes. See Section 2 for descriptions of data and methods.

^a Note that fits versus M_⋆ have their intercept at ${\mathrm{log}}_{10}{M}_{\star }=10$ and fits versus SFR/M_⋆ have their intercept at log₁₀ SFR/M_⋆ = −10. ^b Correlated axes. ^c The row reports the fit to the resolved data with an offset applied to the intercept to account for the +0.19 dex median offset of galaxies with CO (1–0) mapping from the integrated relations. The direct fits describe our resolved CO (1–0) data better, but these corrected fits should better describe the general population of massive, star-forming galaxies captured by the galaxy-integrated CO (1–0) data (Figure 1). See Section 3.4 for a more detailed description of the correction.

Download table as:  ASCIITypeset image

We treat the mid-IR emission as the independent (i.e., x-axis) variable in the analysis. Table 1 shows that our data often include high-S/N mid-IR data with low-S/N CO measurements but that the reverse is almost never true. This makes ordering the analysis around the mid-IR measurements convenient. When comparing the CO-to-mid-IR ratio to M_⋆ or SFR/M_⋆, we treat the latter as the independent variable, since we have SFR and M_⋆ estimates for all targets (Section 2.1).

Treating these variables as independent we bin in x, usually mid-IR emission, almost always using 0.25 dex wide bins, and then measure the median and 16%–84% range for the dependent (y-axis) variable, usually CO intensity. These calculations take into account upper limits for integrated galaxies. For resolved measurements we include the low-S/N data directly, and rely on the averaging performed in binning to yield a good S/N estimate of the median. In the figures we show estimates of the 16%–84% range when the lower bound is not a limit and there are at least 12 detected measurements contributing to the bin. We show the median estimate when it is not a limit and there are at least six measurements contributing. Overall the bins show an excellent match to the ridgeline and extent of the data (e.g., Figure 2). They offer a simple way to characterize general trends for a large data set, and a robust way to deal with the frequently modest S/Ns in the CO mapping data.

Zoom In Zoom Out Reset image size

The I_CO versus I_MIR relation appears to be reasonably well described as a power law, at least to first order (Section 3). We characterize this by fitting a line to the binned data in log–log space. In detail, we use an orthogonal distance minimization to fit

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{I}_{\mathrm{CO}}={m}_{\mathrm{CO \mbox{-} IR}}{\mathrm{log}}_{10}{I}_{\mathrm{IR}}+{b}_{\mathrm{CO \mbox{-} IR}},\end{eqnarray} \tag{ 3 }$$

where I_CO and I_MIR are the median CO line intensity in the bin and the mid-IR surface brightness at the bin center for the line and band in question in units of K km s⁻¹ and MJy sr⁻¹.

The subscript "CO-IR" denotes that we are fitting CO as a function of mid-IR intensity. This choice is driven by the practical considerations discussed above, but it also makes physical sense if one were attempting to build a model to predict CO emission or molecular gas content (e.g., Gao et al. 2019, 2022; Chown et al. 2021). However, this sense is opposite to that often used to characterize star formation scaling relations, which posit some variant of ${{\rm{\Sigma }}}_{\mathrm{SFR}}\propto {{\rm{\Sigma }}}_{\mathrm{gas}}^{N}$, where gas is the independent variable. Because we will also discuss our results with reference to these star formation scaling relations, we note that, algebraically,

$$\begin{eqnarray}&&{m}_{\text{IR-CO}}=\displaystyle \frac{1}{{m}_{\mathrm{CO} \mbox{-} \mathrm{IR}}}\,{\rm{and}}\,{b}_{\text{IR-CO}}=-\displaystyle \frac{{b}_{\mathrm{CO \mbox{-} IR}}}{{m}_{\mathrm{CO \mbox{-} IR}}},\end{eqnarray} \tag{ 4 }$$

and we will reference m_IR-CO as a close cognate of the slope of star formation scaling relations, often denoted N.

When regressing the CO-to-mid-IR ratio against global properties, we fit

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\displaystyle \frac{{I}_{\mathrm{CO}}}{{I}_{\mathrm{MIR}}}=m\left({\mathrm{log}}_{10}{M}_{\star }-{x}_{0}\right)+b,\end{eqnarray} \tag{ 5 }$$

with M_⋆ in M_⊙ or use the analogous equation substituting SFR/M_⋆ yr⁻¹ for M_⋆. Here, x₀ is an offset that allows us to quote the intercept near the center of the data range in order to decrease covariance between m and b. We use x₀ = 10 when fitting M_⋆ and x₀ = −10 when fitting SFR/M_⋆.

We note a few points. First, our choice to fit the binned data reflects the aim to obtain a general scaling rather than the most precise possible predictor of the existing distribution. In practice, our approach deemphasizes the large amount of high-quality data near the distribution center (see Figure 1 and Table 2), and weights more the data at low and high mid-IR emission. Indeed, inspection shows that our fits to bins do an excellent job at capturing the ridgeline of the data over the intensity range spanned by our bins, which is our goal. The slopes of most of these relations can still easily change by ±0.1 by adopting other methodologies, and our choice to bin by mid-IR and then fit the binned data represents the key choice. We emphasize, however, that we apply exactly the same method to all variable pairs. This allows us to robustly compare results across different mid-IR bands and lines, even when the absolute values have methodological uncertainty. We make an approximate estimate of uncertainties associated with our fits after binning by adding in quadrature the result of (a) bootstrap resampling the bins, and (b) varying the fitting approach among orthogonal distance, ordinary least-squares y versus x, and ordinary least-squares x versus y. We also report the range of x values over which we make binned measurements. This represents the range over which we expect the fit to apply.

In addition to fitting, for each data pair we also measure the Spearman's rank correlation relating x and y for detected (S/N > 5) data, the median ratio ${\mathrm{log}}_{10}y/x$, and the scatter in both the ratio and the residuals of individual data about the fit. The main caveat is that the individual data at 17'' resolution are quite noisy, especially for CO (1–0), resulting in a weaker correlation at these scales. Therefore the radial profile data offer our best high-S/N view of moderately resolved galaxies.

Finally, we note that relations with internally correlated axes should be interpreted carefully, because one risks ascribing physical meaning to something that may be purely mathematical. Our only formally correlated variable pair is CO/MIR versus SFR/M_⋆, where the mid-IR has been used as part of the SFR estimates. The relationship is of interest, and so we plot it, but we caution some of the observed anticorrelation may reflect the built-in correlation due to the variable choice. The rest of our axes are formally independent, though obviously closely linked physically.

Glancing through Table 3 and Figures 2–6 shows our most basic result: for both integrated and resolved galaxies there is a tight relation between I_CO and I_MIR for all scales and choices of band. For both integrated galaxies and radial profiles, the rank correlation coefficients for the relation of CO-to-mid-IR are all ≥0.80 and often >0.9. The lower correlation coefficients for individual regions simply reflect the lower S/Ns in those data. The binned relations are well fit by power laws with slopes in the range ∼0.7–1.2, i.e., close to linear, with the exact value depending on the mid-IR band and CO line choice. Across our integrated and resolved data sets, the correlations span at least three orders of magnitude in intensity from I_MIR ∼ 0.1 MJy sr⁻¹ to I_MIR ≳ 100 MJy sr⁻¹ and I_CO ∼ 0.1 K km s⁻¹ to I_CO ≳ 100 K km s⁻¹. Across this range the scatter in the individual data points about a power-law fit is 0.2–0.3 dex (colored regions in Figure 2).

This CO–mid-IR relationship is among the tightest observational correlations between distinct spectral bands in extragalactic astronomy. Note that these are relations between "distance-independent" surface brightnesses and not luminosities. For luminosities the presence of a correlated distance squared on both axes leads to artificially tighter correlations. Instead, we are presenting these relationships in an "intensive" way, using I_CO and I_MIR. Had we adopted an extensive quantity like luminosity, their tightness and span would rival the far-IR–HCN relation (e.g., Gao & Solomon 2004) or far-IR–radio correlation (e.g., Yun et al. 2001; Bell 2003).

The standard interpretation for the CO–mid-IR correlation is that the mid-IR reflects star formation while CO emission captures the fuel for star formation. Though our focus is observational rather than diagnosing the physical emission mechanisms, we note that optically thick CO emission should be brighter in the presence of significant heating and turbulence, while the optical depth of the mid-IR emission will be higher for larger column densities of gas. In that sense, both bands have some element of tracing both star formation and the gas distribution. The tightness and near-linearity of the correlation lends circumstantial support to the idea that mid-IR and CO often trace closely related parameters. Over the next few years, the ability of JWST and ALMA to resolve galaxies into discrete regions of recent star formation and abundant molecular gas should yield significant insight on this topic (see, e.g., Kim et al. 2023, Leroy et al. 2022b, and Sandstrom et al. 2023 in the "First Results from PHANGS–JWST" Focus Issue).

Table 3 reports our best-fit CO-to-IR slope for each CO line-IR band pair, and we visualize these together in Figure 7. For the most part the results for integrated galaxies, azimuthal averages, and individual regions agree well for a given line-band combination, though there appears to be a mild tendency for the individual region data to show a slightly shallower slope. Given that the integrated-galaxy data sets and the resolved data sets are quite distinct, the good overall agreement gives us confidence that we indeed access the slope of an underlying scaling relation that describes local galaxies fairly well.

There are clear patterns within the slopes for different band-line combinations (Figure 7). For a given IR band, CO (2–1) shows a consistently steeper slope than CO (1–0), which is evident in the offset between red and blue points in Figure 7. Averaging over all bands yields a median offset in slope of ≈ +0.20 from the CO (1–0) relation to the CO (2–1) relation with the same mid-IR band. Several recent studies have found that for local star-forming galaxies the CO (2–1) to CO (1–0) line ratio varies with the local SFR surface density or IR luminosity density as approximately ${R}_{21}\equiv {I}_{\mathrm{CO}}^{2-1}/{I}_{\mathrm{CO}}^{1-0}\propto {{\rm{\Sigma }}}_{\mathrm{SFR}}^{0.15}$ (den Brok et al. 2021; Yajima et al. 2021; Egusa et al. 2022; Leroy et al. 2022). Since the mid-IR is often used to trace Σ_SFR linearly (see Equation (2)), our finding agrees well with these results. In essence, our statistical measurement shows that ${R}_{21}\propto {I}_{\mathrm{MIR}}^{0.20}\sim {{\rm{\Sigma }}}_{\mathrm{SFR}}^{0.20}$ for a very large sample of local galaxies, complementing the more direct studies of small samples or individual galaxies listed above.

We also see offsets among the fitting results for different mid-IR bands. For a fixed choice of line, 12 μm shows a steeper m_CO-IR slope than 22 and 8 μm shows a steeper m_CO-IR slope than 24 μm. Across all lines and samples, the median difference between m_CO-IR for 12 μm and that for 22 μm is +0.18 and the slope difference between 8 μm and 24 μm is also +0.18. The 12 and 8 μm bands both contain a significant (often dominant) contribution from PAH band features, ⁴⁴ while the 22 and 24 μm intensities represent mostly continuum emission. This suggests that the PAH features contribute to systematically steeper CO–mid-IR slopes than the mid-IR continuum (a finding in very good agreement with spectroscopy-based results from Whitcomb et al. 2022).

A possible (albeit nonunique) explanation for this effect would be that the intensity of the PAH emission traces the interstellar medium (ISM) more directly, including perhaps non-CO-emitting gas, while the slightly longer wavelength continuum emission tracks star formation better. Indeed, high-spatial-resolution imaging in the Milky Way shows that the hot dust 24 μm emission appears directly related to star formation, while the diffuse PAH emission usually surrounds the star-forming region (e.g., Watson et al. 2008). A similar explanation has been proposed by Chown et al. (2021) and Whitcomb et al. (2022), and we might detect statistical evidence for this scenario here. Note, however, that differences in how the different mid-IR bands trace the ISM and star formation is likely to be one of degree, with all bands tracing heating by recent star formation or gas to one degree or another. At intermediate scales, PAH-dominated bands still correlate with tracers of star formation and diffuse, stochastically heated dust can still emit mid-IR continuum emission (see Leroy et al. 2022b).

The sense of a steeper CO–mid-IR slope is that there will be less CO relative to mid-IR emission in the faint, diffuse parts of our data set. This could reflect that lower-surface-brightness regions within galaxies are expected to hold a larger fraction of their gas in atomic form and perhaps to have a larger fraction of their molecular gas in a CO-dark phase (e.g., see reviews in Bolatto et al. 2013 and Saintonge & Catinella 2022). These regions are very faint in CO but will still host PAH emission (e.g., Chastenet et al. 2019; Sandstrom et al. 2023).

As discussed in Sections 1 and 2.4, because of the widespread use of the mid-IR in SFR estimates, the slopes that we measure relate closely to the literature on star formation scaling relations (e.g., Kennicutt et al. 2007, Bigiel et al. 2008, Kennicutt & Evans 2012, and Leroy et al. 2013, among many others, all rely heavily on 24 μm data). Our measurements only relate CO and mid-IR emission. They do not include any term accounting for unobscured UV or Hα emission or any correction for "cirrus" contributions to the mid-IR (i.e., IR emission due to heating from older stellar populations) nor do we consider variations in the CO-to-H₂ conversion factor, α_CO. They are thus not physical estimates of the Σ_gas–Σ_SFR relation but direct observables.

What about using CO as a predictor of mid-IR? Using Equation (4), our measured slopes correspond to m_IR-CO ∼ 0.9–1.4. For a given IR band, m_IR-CO is on average 0.27 lower for CO (2–1) compared to CO (1–0). As discussed above, the PAH-bearing bands will yield shallower mid-IR versus CO slopes, with 12 μm slopes for m_IR-CO on average 0.25 lower than those for 22 μm.

Finally, we remark that the relationship between 12 μm and CO (2–1) appears nearly linear, with slope m_CO-IR almost exactly unity and a rank correlation coefficient >0.9 in both our integrated-galaxy sample and our radial profile sample. In analyzing PHANGS–ALMA, Leroy et al. (2021b) identified the 12 μm as the band best correlated with the CO (2–1) out of several tracers of stellar mass and star formation. The very tight observational relationship appears to hold both broadly and for local star-forming galaxies. Therefore 12 μm emission appears to offer the best estimator of the CO (2–1) brightness of a galaxy or part of a galaxy.

We determine the normalization of the I_CO–I_MIR relation, which we report as b_CO-IR in Table 3, and compare across bands and transitions in Figure 8. We see again overall similarity, but also real differences as we vary the choice of line and band.

Before examining the differences, we note that Figure 8 implies a useful zeroth-order approximation that can be used for quick estimates:

$$\begin{eqnarray}&&\displaystyle \frac{{I}_{\mathrm{CO}}}{{I}_{\mathrm{MIR}}}\sim \displaystyle \frac{1\,{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}}{1\,\mathrm{MJy}\,{\mathrm{sr}}^{-1}}.\end{eqnarray} \tag{ 6 }$$

Given that the slopes we derive are all roughly near unity, this is a useful way to estimate CO from mid-IR intensity or vice versa. This rule of thumb is illustrated as a vertical black line in Figure 8, which shows that it can be expected to work at the factor of ∼2 level. For detailed work one should use the more accurate relations for specific line-band combinations in Table 3.

In more detail, for a given IR band, we find the normalization, b_CO-IR, to be higher for CO (1–0) than CO (2–1). This is even more noticeable for the galaxies resolved in CO (1–0) (e.g., black stars in Figure 2). As we discuss in Section 3.4, it seems likely that the root cause of that feature is that the local galaxy CO (1–0) mapping surveys are biased toward CO-bright galaxies. After accounting for this bias, as we describe in Section 3.4, we find b_CO-IR for CO (1–0) to be a median 0.18 dex higher than b_CO-IR for CO (2–1) and the same band. The inverse of this factor, 0.65, closely resembles the typical R₂₁ line ratio measured for nearby star-forming galaxies: R₂₁ ≈ 0.6–0.7 in recent works, spanning range from ≈0.5 to 0.85 (e.g., den Brok et al. 2021; Yajima et al. 2021; Leroy et al.2022). Thus the offset between lines (either after correcting the resolved CO (1–0) sample or from only using the integrated galaxies) appears fully consistent with being simply driven by the average CO line ratio (excitation).

The individual IR bands also differ, reflecting their different brightnesses and consistent with the differences in slopes that we see in Section 3.2. Assuming that these normalizations reflect only the underlying band ratios, this is consistent with 8 μm on average being the brightest band studied, showing intensity ≈1.3 times larger than 22 μm or 24 μm, and ≈1.6 times larger than 12 μm.

Just as the slope links to common star formation scaling relations, the CO-to-mid-IR ratio reflected by b_CO-IR has a close link to the molecular gas depletion time, ${\tau }_{\mathrm{dep}}^{\mathrm{mol}}\equiv {M}_{\mathrm{mol}}/\mathrm{SFR}$. This quantity captures the normalized rate of star formation per unit molecular gas mass, is of broad interest to understanding star formation in galaxies, and is often estimated using combinations of data that include CO line and mid-IR emission (e.g., Leroy et al. 2008; Janowiecki et al. 2017). We note that following the approximate translations to mass and SFR in Equations (1) and (2), this ratio translates to a molecular gas depletion time of τ_dep ≡ M_mol/SFR ∼ 0.8–1.7 Gyr depending on the exact choice of mid-IR band and CO line. This agrees very well with a wide range of recent estimates of these quantities (see, e.g., Leroy et al. 2013; Saintonge et al. 2017; Saintonge & Catinella 2022).

We caution that both these normalizations and our measured slopes will show some scale dependence. A number of studies have shown that at high resolution tracers of star formation and molecular gas in galaxies resolve into distinct distributions. This leads relationships measured on large scales to show increased scatter on small scales (e.g., Schruba et al. 2010; Leroy et al. 2013; Chevance et al. 2020; Kim et al. 2022; Pan et al. 2022). For the case of the mid-IR this "breakdown" in the low-resolution scaling relation appears to be somewhat tempered by the fact that mid-IR emission can act as both a star formation and a gas tracer. In first results comparing PHANGS–JWST (Lee et al. 2023) to PHANGS–ALMA and PHANGS–MUSE (Emsellem et al. 2022), Leroy et al. (2022b) compare CO and Hα to mid-IR emission on ∼70–160 pc scales. Both CO and Hα exhibit strong relationships with the mid-IR at these scales, even though the relationship between CO and Hα shows clear signs of breaking down. As a result, separate relationships appear to relate CO to mid-IR and Hα to mid-IR at high resolution. Analyzing spectroscopic mid-IR data from SINGS, Whitcomb et al. (2022) come to a similar conclusion, that the mid-IR acts as both a gas and a SFR tracer to varying degrees. The relationships that we find here still change at high resolution, but rather than breaking down entirely it appears to resolve into these two separate relationships. CO does still correlate with the mid-IR, simply with a distinct functional relationship and a higher CO-to-mid-IR ratio (see Leroy et al. 2022b).

We have so far focused on the intensity–intensity scaling, but the star formation and molecular gas in galaxies also relates to integrated mass, metallicity, and specific SFR of the system (e.g., see reviews in Tacconi et al. 2020 and Saintonge & Catinella 2022). In Figures 9 and 10 we show the CO-to-mid-IR ratio as a function of M_⋆ and SFR/M_⋆. Table 3 reports the results of these comparisons along with fits to the intensity-scaling relations.

The plots and table show a weak but clear dependence of the CO-to-mid-IR ratio on stellar mass. The sense of the correlation is that there is more CO emission relative to mid-IR emission in more massive galaxies. Lower-mass galaxies also have lower metallicities and fainter CO emission (see Bolatto et al. 2013; Sánchez 2020; Saintonge & Catinella 2022). In principle, both the dust producing mid-IR emission and the molecular gas producing CO become less abundant at low stellar mass (see Bolatto et al. 2013; Galliano et al. 2018). The results here suggest that CO is more quickly depressed in low-mass galaxies than the mid-IR emission, with a bes- fit slope against stellar mass of m ∼ 0.2.

There are a number of plausible explanations for this trend. The star formation per unit gas mass may be higher in low-mass galaxies (and SFR/M_⋆ certainly anticorrelates with M_⋆). More massive galaxies may also be more likely to have a significant fraction of their CO emission from bar-fed nuclear starburst regions. The gas in these regions tends to be more emissive due to increased temperature, dynamically driven line widths, and low opacity, which may cause an increase in CO luminosity for a given molecular mass (e.g., Bolatto et al. 2013; Teng et al. 2022) and lead to a lower CO-to-mid-IR ratio. Alternatively, a larger fraction of mid-IR emission may emerge from dust associated with CO-dark gas and atomic gas in low-mass galaxies compared to high-mass galaxies. There is no a priori physical reason for mid-IR emission to only trace CO-bright molecular gas (e.g., Walterbos & Schwering 1987), so we consider this explanation a likely one. Sandstrom et al. (2023) address this issue in the first results from PHANGS–JWST, and more future highly resolved analysis using JWST and ALMA should help shed light on this topic.

Figure 10 and Table 3 also report the scaling of I_CO/I_MIR with specific star formation rate, SFR/M_⋆. We do caution that because our SFR estimates involve 22 μm emission the axes in this plot are anticorrelated (this is the only such case in the Letter). The plot suggests an overall anticorrelation, especially at 22 μm, with the sense that CO drops relative to mid-IR emission as SFR/M_⋆ increases. Taken at face value, this agrees well with a series of recent results that show that the depletion time of molecular gas anticorrelates with SFR/M_⋆ and distance to the star-forming main sequence (Saintonge et al. 2016; Colombo et al. 2018, 2020). In simple terms, molecular gas as traced by CO appears less efficient at forming stars in galaxies with low specific SFR and more efficient in starburst galaxies. The anticorrelation of CO-to-mid-IR with SFR/M_⋆ in Table 3 and Figure 10 is consistent with this finding, though a more detailed statistical analysis controlling for the correlated axes would help validate our simple correlation measurement.

Figure 10 shows steeper slopes for the CO-to-22 μm ratio correlation with SFR/M_⋆ than for the CO-to-12 μm ratio. Assuming that this is not merely a chance byproduct of internal correlations, it may reflect that the 12–22 μm ratio is varying in response to destruction of PAHs in regions of high specific SFR and intense star formation. This phenomenon is observed within local galaxies (e.g., Chastenet et al. 2019), and the 12 μm-to-22 μm does drop with increasing SFR/M_⋆ (e.g., see Figure 21 in Leroy et al. 2019).

Finally, Figures 9 and 10, along with Figure 2, allow us to examine the selection function for the local mapping surveys used for the resolved analysis. In these three figures we mark the local CO mapping targets as black stars. The surveys that we use for the integrated CO sample reflect the massive, star-forming galaxy population well by design (Bolatto et al. 2017; Saintonge et al. 2017; Colombo et al. 2020; Lin et al. 2020; and see Section 2 and Figure 1), and WISE 12 and 22 μm imaging covers the entire sky, so we expect the correlations arising from those data to suffer the least from sample bias. The figures show that, overall, the CO (2–1) mapping surveys mostly sample the same parameter space as the integrated-galaxy measurements. Indeed, the CO (2–1) sample is heavily dominated by PHANGS–ALMA, which attempted a volume-limited selection of relatively massive galaxies on or near the star-forming main sequence and in that sense represents well the population of star-forming galaxies with ${\mathrm{log}}_{10}{M}_{\star }\,\approx$ 9.5–11 M_⊙. By contrast, the CO (1–0) mapping surveys appear systematically offset toward high CO brightness in all three plots. The local CO (1–0) mapping sample is dominated by COMING, and Sorai et al. (2019) do explicitly state that CO brightness is a key criterion of the selection, so this bias might reasonably be expected.

We attempt a first-order correction for this bias by solving for the median offset between the mapped points (i.e., the black stars) and the best-fit line from the full integrated data set in each fit relation in Figures 2, 9, and 10. On average, the median offset required to bring the CO (1–0) mapping sample in line with the full integrated-galaxy sample is ${\rm{\Delta }}{\mathrm{log}}_{10}{I}_{\mathrm{CO}}\,=$ −0.19 dex (a factor of 1.55). We have applied this shift as a "statistical correction" in Table 3 and Figures 3–8. If we solve for the analogous quantity for CO (2–1) we find a median offset of only 0.04 dex (a factor of 1.10), which we do not apply.

In Sections 3.2 and 3.3 we discussed that the differences in slope and normalization observed when using CO (1–0) and CO (2–1) appear to agree well with recent work on the CO (2–1) to CO (1–0) line ratio, R₂₁. In fact, if we assume that the relations measured for each line accurately reflect the same underlying galaxy population, we can make a much more powerful inference. By dividing the scaling relation for I_{CO 2−1} by that for I_{CO 1−0}, we can predict the line ratio as a function of mid-IR intensity. As mentioned above, several recent direct studies of R₂₁ have found a clear dependence of the line ratio on tracers of the star formation surface density, resulting in ${R}_{21}\sim {{\rm{\Sigma }}}_{\mathrm{SFR}}^{0.15}$ (den Brok et al. 2021; Yajima et al. 2021; Leroy et al. 2022). In contrast to these studies, the results here are statistical in nature. On the other hand, given our large sample these relations should describe the general population of massive, star-forming galaxies near the star-forming main sequence, rather than the smaller samples where this ratio has been studied directly. When relevant, we apply the statistical correction described in Section 3.4 to the CO (1–0) mapping results.

Figure 11 visualizes the resulting predictions, which are as follows:

$$\begin{eqnarray}&&{R}_{21}\approx \left\{\begin{array}{l}0.73\,{I}_{12\,\mu {\rm{m}}}^{0.19}\\ 0.62\,{I}_{22\,\mu {\rm{m}}}^{0.22}\\ 0.62\,{I}_{8\,\mu {\rm{m}}}^{0.26}\\ 0.65\,{I}_{24\,\mu {\rm{m}}}^{0.24}\end{array}\right.,\end{eqnarray} \tag{ 7 }$$

with all intensities in MJy sr⁻¹, and the last two relationships being less accurate because they rely only on the much smaller set of targeted Spitzer mapping data. As illustrated in Figure 11, we recommend capping R₂₁ at a value of R₂₁ = 1 when implementing these relationships, which is the expectation for optically thick emission from gas at uniform temperature. As discussed in Leroy et al. (2022), there is also a likely plausible lower bound, but the exact value of this is uncertain. ⁴⁵

Zoom In Zoom Out Reset image size

Although this is a statistical measurement, Equation (7) should be as robust as any currently available general prescription for R₂₁ in the literature.

During Cycle 1, the PHANGS–JWST Treasury program (Lee et al. 2023) will observe mid-IR emission at 7.7 μm (F770W), 10 μm (F1000W), 11.3 μm (F1130W), and 21 μm (F2100W) from 19 local star-forming galaxies. All of these targets have corresponding high-resolution ALMA data (Leroy et al. 2021b). In Figure 12 we compare the 19 targets of the PHANGS–JWST survey to our full low-resolution data set. The PHANGS–JWST targets all have ALMA CO (2–1) and WISE 12 and 22 μm data, though only about half have IRAC 8 μm and MIPS 24 μm data. In the figure, the red points showing data for PHANGS–JWST targets do a reasonable job of spanning the full range of CO and mid-IR intensities seen in the literature sample and also appear broadly consistent with the CO versus mid-IR relations seen for the broader sample.

Zoom In Zoom Out Reset image size

At the time of submission of this Letter, four of the PHANGS–JWST targets have already been observed by the MIRI instrument on JWST (NGC 628, NGC 1365, NGC 7496, and IC 5332). Lee et al. (2023) in the PHANGS–JWST "First Results" Issue report details of the survey design, data reduction, and quality assurance of these MIRI data. Briefly, each pointing was observed for 89 (F770W), 122 (F1000W), 311 (F1130W), and 322 (F2100W) seconds. Reduction mostly followed the standard JWST procedures, with the sky and instrumental background removed by subtracting a nearby "off" field. After processing by the JWST pipeline, the individual pointings were mosaicked together and the different MIRI bands were compared to set a common background level. Finally, the background level of all bands was anchored to previous wide-field observations by either Spitzer or WISE following the procedure described in Leroy et al. (2022b). Overall, the flux calibration of the JWST data appear consistent with these previous observations and the background level is uncertain at the ±0.1 MJy sr⁻¹ level.

In Figure 13, we show results from JWST for the first four PHANGS–JWST targets. The purple points show PHANGS–ALMA CO (2–1) as a function of mid-IR intensity in four JWST filters. ⁴⁶ The blue points show binned versions of the PHANGS–JWST data and the purple line shows a fit to these bins. For comparison, we also plot the data for individual regions from our larger data set. For the F770W, F1000W, and F1130W bands, we show CO (2–1) as a function of 8 μm emission measured by Spitzer.

Zoom In Zoom Out Reset image size

We compare F770W and F1130W, which are expected to be PAH dominated, to results for the PAH-dominated 8 μm band. Empirically, F1000W tracks these other two MIRI bands closely, though whether F1000W is PAH dominated is less clear, so we also compare it to 8 μm. For F2100W we compare to 24 μm emission, because both bands are expected to be continuum dominated at similar wavelengths. The black solid lines show the fits to this larger sample. Table 4 reports the results of our fits to the JWST data. It also reproduces the 8 and 24 μm individual region fitting results from Table 3 for ease of comparison.

Table 4 and Figure 13 show good overall agreement between the slopes in the JWST data and the comparison data. The PAH-tracing bands all show slopes ∼1.3–1.4, similar to the slope of ∼1.2 measured for the 8 μm. Meanwhile the CO as a function of 21 μm shows a slope of ∼0.9, similar to that relating CO (2–1) to 24 μm. At the most basic level, we confirm that the scalings captured at ∼15'' ∼ 1 kpc by PHANGS–JWST quantitatively match those observed in much larger samples at similar bands using Spitzer. At a practical level, this provides a quantitative link between the more focused, high-resolution JWST observations and the broader, low-resolution mid-IR literature. We do caution that, so far, we have only analyzed four galaxies, and because of this the data from JWST span a limited intensity range. Moreover, essentially all of the bright CO and mid-IR emission within this small sample come from the inner region of a single galaxy, NGC 1365 (see, e.g., Schinnerer et al. 2023). With these caveats in mind, the good agreement between such a small sample of JWST results and extensive previous mid-IR mapping in Figure 13 may be even more impressive.

Physically, the slope of ∼1.2–1.4 for the PAH-tracing F770W, F1000W, and F1130W bands indicates that the mid-IR emission appears brighter relative to CO (2–1) emission at low mid-IR intensity. As discussed in Section 3.2, one plausible interpretation of these slopes would be that mid-IR emission emerges from dust mixed with atomic gas or CO-dark molecular gas or low-excitation gas. This suggests the mid-IR may even act as an effective tracer of gas poorly traced by CO.

The resolution and sensitivity of JWST make this a potentially exciting prospect because we have relatively few ways to trace non-CO-emitting neutral gas at high angular resolution. More validation and calibration are needed to unlock the potential of the mid-IR as a quantitative ISM tracer at high resolution (at low resolution, see Gao et al. 2019, 2022 and Chown et al. 2021), but this work is proceeding quickly. Here, at low resolution we have verified the quantitative link between the PHANGS–JWST results and a larger sample. In this Issue, Sandstrom et al. (2023) carefully examine low-intensity regions in the first PHANGS–JWST data and conclude that the level of mid-IR emission appears consistent with arising from atomic, not molecular gas. Also in this Issue, Leroy et al. (2022b) compare PHANGS–ALMA CO (2–1), PHANGS–MUSE Hα (Emsellem et al. 2022), and JWST mid-IR at ∼70–160 pc scales and make a first attempt at isolating the mid-IR associated with diffuse ISM emission, calibrating the gas-to-mid-IR conversion, and separating ISM-tracing emission from emission driven mainly by heating. In a similar effort using a large set of SINGS spectroscopy, Whitcomb et al. (2022) examine how different continuum bands and mid-IR emission features correlate with CO and star formation tracers, and use this analysis to map out the degree to which parts of the mid-IR spectrum trace either gas or star formation. Taking all of this together, we expect that in the near future mid-IR emission may offer a new, powerful complement to traditional ISM probes like H i, CO, or long-wavelength dust emission.