PHANGS–ALMA: Arcsecond CO(2–1) Imaging of Nearby Star-forming Galaxies

Much of our knowledge about the behavior of the molecular interstellar medium (ISM) in z = 0 galaxies has been established by CO surveys that either integrate over whole galaxies (e.g., the FCRAO survey, Young et al. 1995; AMIGA, Lisenfeld et al. 2011; COLD GASS, Saintonge et al. 2011; ALLSMOG, Bothwell et al. 2014; xCOLD GASS Saintonge et al. 2017; and JINGLE, Saintonge et al. 2018) or resolve the large-scale structure of galaxy disks but do not distinguish individual molecular clouds (e.g., BIMA SONG, Helfer et al. 2003; the Nobeyama CO Atlas, Kuno et al. 2007; HERACLES, Leroy et al. 2009; the JCMT NGLS, Wilson et al. 2012; CARMA STING, Rahman et al. 2012; ATLAS–3D CO, Alatalo et al. 2013; Davis et al. 2014; CARMA EDGE, Bolatto et al. 2017; NRO COMING, Sorai et al. 2019; and ALMAQUEST, Lin et al. 2019).

These surveys have demonstrated a close link between molecular gas and star formation, showing that the location and rate of recent star formation in a galaxy tracks the distribution of molecular gas (e.g., Wong & Blitz 2002; Kennicutt et al. 2007; Bigiel et al. 2008; Leroy et al. 2008; Schruba et al. 2011). Yet despite this good overall correspondence, observations reveal important variations in the amount of star formation per unit molecular gas both among different types of galaxies (e.g., Saintonge et al. 2011; Leroy et al. 2013b; Davis et al. 2014; Huang & Kauffmann 2015, and see Figure 1) and within different regions of the same galaxy (e.g., Leroy et al. 2013b, 2017; Longmore et al. 2013; Meidt et al. 2013; Momose et al. 2013; Utomo et al. 2017; Brownson et al. 2020). The normalized CO emission of galaxies also varies, with CO emission appearing fainter relative to starlight or tracers of recent star formation in low-mass and early-type galaxies (e.g., Young & Scoville 1991; Young et al. 1996; Schruba et al. 2012; Hunt et al. 2015; Saintonge et al. 2016, 2017). This change in brightness arises from both changes in molecular gas content and changes in CO emissivity per unit of molecular gas mass, but the relative magnitude of these effects is uncertain (for a review see Bolatto et al. 2013a).

Zoom In Zoom Out Reset image size

Figure 1 illustrates some of these global trends using data from PHANGS–ALMA (red, see Section 4 for details), xCOLD GASS (Saintonge et al. 2017) as the largest homogeneous unresolved survey, and a compilation of local CO mapping surveys (COMING, the Nobeyama CO Atlas, HERACLES, and a HERA follow-up survey; Kuno et al. 2007; Leroy et al. 2009; Sorai et al. 2019; and A. Schruba et al. 2022, in preparation). The top panels show how the ratios of molecular gas mass to star formation rate and molecular gas mass to stellar mass change across the local galaxy population. The lower panels show the relationship between the average surface densities of molecular gas, stars, and star formation. Together, the four panels of Figure 1 demonstrate the good overall correspondence between star formation, stellar mass, and molecular gas in galaxies, but also illustrate important variations in the molecular content of galaxies normalized by size, star formation rate, or stellar mass. The abundance, structure, and ability of molecular gas to form stars varies across the z = 0 galaxy population.

Unfortunately, low-resolution observations offer limited insight into the physical state of molecular gas. In the Milky Way and its Local Group neighbors, most of the molecular gas resides in GMCs with masses ∼10⁴–10⁷ M_⊙. These clouds have sizes of tens of parsecs and appear dominated by supersonic turbulence (e.g., Solomon et al. 1987; Blitz et al. 2007; Fukui & Kawamura 2010; Roman-Duval et al. 2010; Gratier et al. 2012; Heyer & Dame 2015; Rice et al. 2016; Miville-Deschênes et al. 2017; Schruba et al. 2019).

GMCs do not fill the galaxy disk. In low-resolution extragalactic observations like those mentioned in the previous section, the CO emission from GMCs is diluted with nearby non-CO emitting regions. That is, the intrinsic distribution of CO emission in galaxies is strongly clumped on scales much smaller than the ≳kiloparsec resolution of the previous generation of large CO mapping surveys (e.g., Leroy et al. 2013a). Figure 2 illustrates this phenomenon by showing CO emission from two PHANGS–ALMA targets at two resolutions: 1 kpc and 120 pc resolution. The sharp, clumpy structure that is striking in maps at high resolution (right panels), blurs into faint, low-contrast structures when observed at kiloparsec resolution (left panels).

Zoom In Zoom Out Reset image size

Current models of star formation predict a link between star formation, feedback, and gas properties on the scale of individual GMCs, which can be inferred using high-resolution observations. For example, the mean density and density distribution within a cloud may set the characteristic timescale for star formation (e.g., Padoan & Nordlund 2002; Krumholz & McKee 2005; Hennebelle & Chabrier 2011; Federrath & Klessen 2012; Krumholz & Dekel 2012). The strength of self-gravity in the cloud relative to turbulence and magnetic fields may affect the efficiency with which gas is converted into stars (e.g., Padoan et al. 2012, 2017; Burkhart 2018; Kim et al. 2021b). The density, turbulence, and self-gravity may also determine how the new-born stars cluster (e.g., Kruijssen 2012; Hopkins 2013; Krumholz & McKee 2020; Grudić et al. 2021). The local gas (column) density distribution may also interact with sources of stellar feedback to determine whether a cloud is disrupted or not and how much gas and radiation leaves the system (e.g., Thompson et al. 2005; Walch et al. 2015; Geen et al. 2016; Raskutti et al. 2016, 2017; Thompson & Krumholz 2016; Kim et al. 2018, 2019; Reissl et al. 2018; Geen et al. 2021). Because the timescale, efficiency, and spatial clustering of star formation and feedback have a qualitative impact on the GMC-scale gas properties (e.g., Hopkins et al. 2013; Gentry et al. 2017; Keller & Kruijssen 2020), star formation, stellar feedback, and GMC properties form a complex, multiscale system with many types of physics at play.

Dynamical processes acting on ∼10–1000 pc scales also play a key role in setting the abundance, structure, and ability of molecular gas to form stars (for recent reviews see Dobbs et al. 2014; Krumholz 2014; Chevance et al. 2020c). As a concrete example, the high-resolution images in Figure 2 show the unmistakable imprint of a stellar bar and spiral arms in both galaxies, which is far less obvious in the low-resolution maps of these targets. Spiral arm passage may collect individual quiescent molecular clouds into large star-forming associations (e.g., Koda et al. 2009; Meidt et al. 2015), or trigger phase changes from the atomic to molecular medium (e.g., Dobbs et al. 2014). This can trigger star formation (e.g., Egusa et al. 2017) and/or organize the star-forming structures (e.g., Schinnerer et al. 2017; Elmegreen et al. 2018; Kim et al. 2020; Tress et al. 2020). Meanwhile, gas flows along stellar bars and arms may prevent collapse of the streaming gas, suppressing star formation (e.g., Meidt et al. 2013). These same streaming motions can also redistribute the gas, fuel star formation in the inner parts of the galaxy, or even trigger nuclear starbursts (e.g., Kenney et al. 1992; Sakamoto et al. 1999; Sheth et al. 2005; Schmidt et al. 2016). Collisions between gas clouds may also trigger star formation (e.g., Tan 2000; Inoue & Fukui 2013; Fukui et al. 2021).

A further advantage of high-resolution observations is that they give access to the temporal domain of interstellar processes. Star formation is a dynamic process, with clouds evolving rapidly under the influence of both gravity (e.g., Elmegreen 2000) and "stellar feedback" (e.g., Lee et al. 2016; Kruijssen et al. 2019; Chevance et al. 2020a), a term used to describe the combined influence of ionizing photons, direct and indirect radiation pressure, gas heating, stellar winds, and supernova explosions (e.g., Lopez et al. 2011, 2014; Dale 2015; Rahner et al. 2017, 2019). Over the last decade it has been recognized that the details of the various stellar feedback mechanisms have a large impact on molecular gas properties, star formation rates, and galaxy evolution (e.g., Hopkins et al. 2012; Agertz et al. 2013; Walch et al. 2015; Klessen & Glover 2016; Kim & Ostriker 2017). But many such details remain poorly constrained by observations. The interplay between turbulence and gravity also remains imperfectly understood, with models variously positing short-lived clouds in near freefall (e.g., Vazquez-Semadeni 1994; Elmegreen 2000; Hartmann et al. 2001; Ibáñez-Mejía et al. 2016), star formation proceeding at a steady pace (e.g., Krumholz & McKee 2005), or steadily accelerating star formation within clouds (e.g., Murphy et al. 2011; Lee et al. 2016).

When observed at sufficient resolution, molecular gas, H ii regions, stellar clusters, and other tracers of star formation and feedback visibly separate (e.g., Kawamura et al. 2009; Onodera et al. 2010; Schruba et al. 2010; Gratier et al. 2012; Schinnerer et al. 2013; Corbelli et al. 2017). To illustrate this effect, we overplot the distribution of Hα and CO emission for four PHANGS–ALMA targets in Figure 3. The distributions of CO, tracing GMCs, and of Hα, tracing recent star formation, mostly track each other at large scales, but are clearly distinct at a spatial resolution of a few times 10 pc. Comparing GMCs to stellar tracers with well-understood ages or life spans allows one to infer the timescales for star formation and feedback from high-resolution imaging (e.g., Kawamura et al. 2009; Kruijssen & Longmore 2014; Corbelli et al. 2017; Kruijssen et al. 2018). This offers the prospect to build a picture of the evolutionary sequence of star formation (e.g., Murray 2011; Lee et al. 2016; Kruijssen et al. 2019), to constrain the feedback mechanisms responsible for cloud destruction (e.g., Chevance et al. 2020b), and to assess the fraction of non-star-forming molecular material (e.g., Schinnerer et al. 2019; Chevance et al. 2020a; Kim et al. 2021a).

Zoom In Zoom Out Reset image size

High-resolution CO observations measure the physical state of the gas, probe crucial dynamical processes, and constrain the timescales for star formation and feedback. So far, most cloud scale studies of normal, nonstarburst galaxies have targeted members of the Local Group (see the review by Fukui & Kawamura 2010). During the past two decades, there have been high spatial resolution, wide-field mapping surveys of the CO emission in the Magellanic Clouds (Fukui et al. 1999; Mizuno et al. 2001; Wong et al. 2011), M31 (Nieten et al. 2006; Rosolowsky 2007; Schruba et al. 2021), M33 (Engargiola et al. 2003; Rosolowsky et al. 2007; Onodera et al. 2012; Druard et al. 2014), and various Local Group dwarf galaxies (e.g., Leroy et al. 2006), as well as small-area mapping of Local Group targets with ALMA (e.g., Rubio et al. 2015; Schruba et al. 2017; Wong et al. 2019).

Unfortunately, the range of galaxy types and dynamical environments in the Local Group is limited. With one massive early-type spiral and a modest number of dwarf galaxies, the Local Group is not representative of the galactic environments where most star formation occurs at z = 0. Local Group galaxies also do not harbor the environmental extremes found in more distant galaxies. While their proximity offers significant advantages in terms of resolution and surface brightness sensitivity, observations in the LMC, M33, and M31 cannot capture the full range of behavior seen in Figure 1.

The lack of diversity in the Local Group is problematic because we know from the low-resolution surveys discussed above that the amount and behavior of molecular gas is closely linked to properties of the host galaxy. The balance between atomic gas and molecular gas depends sensitively on the interstellar pressure and local dust content (e.g., Wong & Blitz 2002; Blitz & Rosolowsky 2006; Leroy et al. 2008; Wong et al. 2013; Schruba et al. 2018; Sun et al. 2020a). The distribution of molecular gas strongly reflects the structure of the stellar disk (e.g., Young et al. 1995; Regan et al. 2001; Leroy et al. 2008; Schruba et al. 2011). These trends also hold across the whole galaxy population (e.g., Young & Scoville 1991; Young et al. 1995; Saintonge et al. 2011, 2017), since the molecular gas content, or at least the CO emission, of a galaxy depends strongly on its mass and metallicity (e.g., Schruba et al. 2012; Bothwell et al. 2014; Hunt et al. 2015; Saintonge et al. 2017).

Mapping CO emission at the scale of individual GMCs across a diverse, representative sample of star-forming galaxies is thus the logical step forward in this field. Before ALMA, however, mapping GMC-scale CO emission from a single normal star-forming galaxy required a major time investment. As a result, mapping studies—especially with the PdBI and OVRO interferometers—typically focused on bright, compact starburst galaxies (e.g., Downes & Solomon 1998) and nuclear regions hosting starbursts and active galactic nuclei (e.g., NUGA and MAIN; García-Burillo et al. 2003; Jogee et al. 2005). After a number of studies targeting individual galaxies or dwarf galaxies (e.g., Rosolowsky & Blitz 2005; Bolatto et al. 2008; Rahman et al. 2011), the CANON survey took an important first step toward synthetic cloud scale imaging of a sample of normal galaxies beyond the Local Group. CANON mapped CO(1−0) emission at ∼2'' resolution over the inner regions of a sample of spiral galaxies (Koda et al. 2009; Donovan Meyer et al. 2012, 2013; Momose et al. 2013). This survey provided important evidence for variations in molecular gas properties as a function of galactic environment, and for the role of spiral arms in GMC formation and evolution.

Subsequently, the PdBI Arcsecond Whirlpool Survey (PAWS; Pety et al. 2013; Schinnerer et al. 2013) mapped M51 at 1'' ≈ 40 pc resolution (improving on a similar effort at 3'' resolution using CARMA by Koda et al. 2009). At this high resolution, the contrast between the molecular gas in M51 and that of Local Group galaxies proved striking (e.g., Hughes et al. 2013a, 2013b). Figure 4 illustrates a similar contrast. It shows that at fixed 150 pc resolution, the surface density and line width of molecular gas vary significantly and systematically as a function of location in the galaxy and among host galaxies. Analysis of the PAWS data helped establish that the cloud scale structure of the molecular ISM depends on dynamical environment and host galaxy properties (Hughes et al. 2013a, 2013b; Colombo et al. 2014b; Leroy et al. 2016) and showed that the local star formation activity in M51 depends on cloud scale ISM structure (Meidt et al. 2013; Leroy et al. 2017). PAWS studies also demonstrated how high-resolution imaging yields insight into the evolution and timescales of individual star-forming regions (Schinnerer et al. 2013; Meidt et al. 2015; Schinnerer et al. 2017).

Zoom In Zoom Out Reset image size

In spite of these important first efforts, the sample size, field of view, and sensitivity of high physical resolution CO observations targeting normal, star-forming galaxies remained limited before ALMA, preventing a synthetic view of molecular gas properties across the full local galaxy population.

We convert from IRAC 3.6 μm or WISE 3.4 μm intensity to stellar mass surface density, Σ_⋆, by multiplying the measured intensity with a local estimate of the mass-to-light ratio, ${{\rm{\Upsilon }}}_{\star }^{3.4}$, so that

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Sigma }}}_{\star }}{1\,{{M}_{\odot }\,\mathrm{pc}}^{-2}}\approx 350\left(\displaystyle \frac{{{\rm{\Upsilon }}}_{\star }^{3.4}}{0.5}\right)\left(\displaystyle \frac{{I}_{3.6\,\mu {\rm{m}}}}{1\,\mathrm{MJy}\,{\mathrm{sr}}^{-1}}\right)\,\cos i\end{eqnarray} \tag{ 1 }$$

for IRAC1 data or

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Sigma }}}_{\star }}{1\,{{M}_{\odot }\,\mathrm{pc}}^{-2}}\approx 330\left(\displaystyle \frac{{{\rm{\Upsilon }}}_{\star }^{3.4}}{0.5}\right)\left(\displaystyle \frac{{I}_{3.4\,\mu {\rm{m}}}}{1\,\mathrm{MJy}\,{\mathrm{sr}}^{-1}}\right)\,\cos i\end{eqnarray} \tag{ 2 }$$

for WISE data. The $\cos i$ factor accounts for the inclination of the galaxy and ${{\rm{\Upsilon }}}_{\star }^{3.4}$ is the near-infrared mass-to-light ratio in units of ${M}_{\odot }\,{L}_{\odot }^{-1}$. The slightly different prefactors between Equations (1) and (2) reflect differences in the bandpass of the two instruments. We assume that modulo the different prefactors, the same mass-to-light ratio ${{\rm{\Upsilon }}}_{\star }^{3.4}$ describes both 3.4 μm and 3.6 μm. We predict ${{\rm{\Upsilon }}}_{\star }^{3.4}$ from an empirical fit relating ${{\rm{\Upsilon }}}_{\star }^{3.4}$ to the local SFR-to-WISE1 color from Leroy et al. (2019),

$$\begin{eqnarray}{{\rm{\Upsilon }}}_{\star }^{3.4}\,\left[{M}_{\odot }\,{L}_{\odot }^{-1}\right]=\left\{\begin{array}{ll}0.5 & \mathrm{if}\,Q\lt a\\ 0.5+b\left(Q-a\right) & \mathrm{if}\,a\lt Q\lt c\\ 0.2 & \mathrm{if}\,Q\gt c\end{array}\right.\end{eqnarray} \tag{ 3 }$$

where $Q={\mathrm{log}}_{10}\mathrm{SFR}/({\nu }_{{\rm{W}}1}{L}_{\nu ,{\rm{}}\,{\rm{W1}}})$ and ν_W1 refers to the frequency corresponding to the 3.4 μm central wavelength of the WISE1 band. Here a = −11.0, b = −0.375, and c = −10.2.

Placing PHANGS–ALMA in a larger context: Leroy et al. (2019) derived Equation (3) from an empirical fit relating the measurements and fitting results in the GSWLC. Salim et al. (2016) fit multiband near-ultraviolet, optical, and IR photometry using the CIGALE population synthesis code (Boquien et al. 2019) and the Bruzual & Charlot (2003) models to estimate SFR and M_⋆ for ≳650,000 SDSS galaxies. They showed excellent agreement with the stellar masses in the previous standard JHU-MPA Value Added Catalog (see Kauffmann et al. 2003; Tremonti et al. 2004). ⁴⁶ Using a similar framework, Leroy et al. (2019) also estimated M_⋆ and SFR for ∼15,000 galaxies likely to lie within 50 Mpc. Thus, adopting this approach to estimate stellar mass (and SFR below) maximizes our ability to place PHANGS–ALMA objects in a larger context. In addition to matching the SDSS main galaxy sample via the GSWLC, this specifically places PHANGS–ALMA on a nearly matched mass scale to the GASS and xCOLD GASS surveys (Saintonge et al. 2017; Catinella et al. 2018). These also built on SDSS stellar mass estimates from the JHU-MPA Value Added Catalog and leveraged combinations of GALEX and WISE for SFR estimates. Because we also construct these estimates for all local galaxies with a CO map, we are in a position to compare these galaxies, PHANGS–ALMA, and a larger set of literature galaxy property measurements.

Using SFR/M_⋆ to predict the mass-to-light ratio: Equation (3) specifically uses the ratio of SFR-to-WISE1 luminosity, i.e., a "specific star formation rate−like" quantity, to predict ${{\rm{\Upsilon }}}_{\star }^{3.4}$. In the GSWLC as well as other theoretical and empirical studies (e.g., Bell et al. 2003; Kannappan et al. 2013; Telford et al. 2020), some cognate of the specific star formation rate is a strong predictor of the age of the stellar population and thus variations in the near-infrared mass-to-light ratio. Because calculating the specific star formation rate requires knowing M_⋆, this creates a degenerate situation. Therefore we predict ${{\rm{\Upsilon }}}_{\star }^{3.4}$ from the ratio of SFR-to-WISE1 luminosity. Based on comparison with the GSWLC, this predicts ${{\rm{\Upsilon }}}_{\star }^{3.4}$ with ∼0.1 dex scatter. See Appendix A of Leroy et al. (2019) for more details.

Radial Profile-based ${{\rm{\Upsilon }}}_{\star }^{3.4}$ : We use Equation (3) to predict a radial profile of ${{\rm{\Upsilon }}}_{\star }^{3.4}$. We first match the resolution of all these data sets at 15''. Then, we calculate radial profiles of WISE1 emission and SFR, either from FUV+WISE4, NUV+WISE4, or only WISE4 emission as available. Then we apply Equation (3) to all rings. We smooth the resulting profile with a 2 kpc kernel to avoid noise in the estimate. In the outer parts of galaxies where we do not recover SFR at sufficient signal-to-noise ratio, we adopt the median ${{\rm{\Upsilon }}}_{\star }^{3.4}$ calculated outside 0.5 r₂₅.

This radial profile−based approach reflects the reality that stellar populations change as a function of galactocentric radius (e.g., Watkins et al. 2016; Dale et al. 2020, among many other examples). Because of the use of low-resolution azimuthal averages, it will be less sensitive to mass-to-light ratio variations associated with young stars or dust, e.g., in spiral arms or bars (see Querejeta et al. 2015). Overall, we found that the azimuthal averages offered a good compromise between isolating distinct regions of a galaxy and achieving a good signal-to-noise ratio. This refines the approach in Leroy et al. (2019). They used global colors that trace SFR/M_⋆ to estimate ${{\rm{\Upsilon }}}_{\star }^{3.4}$. The largest difference arises for targets with powerful nuclear sources, either starburst or AGN, which can affect the global colors. Any such bright, contaminating source can bias the global colors to suggest a high SFR/M_⋆ and effectively obscure an underlying older stellar population.

Validation against PHANGS–MUSE: PHANGS–MUSE (Emsellem et al. 2021; P.I. E. Schinnerer) mapped 19 PHANGS–ALMA targets with the MUSE integral field unit on the Very Large Telescope. PHANGS–MUSE delivers estimates of the stellar mass based on fitting the 4850–7000 Å spectrum of individual regions using stellar population synthesis models. Specifically, the fitting uses the E-MILES libraries (Vazdekis et al. 2016) and adopts the same Chabrier (2003) initial mass function (IMF) as the Salim et al. (2016, 2018) fits, so that the IMF is matched across both treatments. These fits offer stellar mass estimates independent of the near-infrared based approaches above. The PHANGS–MUSE data cover only a limited number of targets and do not span the entire stellar disk of their targets. Therefore, we make use of these to explore the uncertainty in our ${{\rm{\Upsilon }}}_{\star }^{3.4}$ and Σ_⋆ values by comparing to an independent estimate with a distinct set of systematic uncertainties. After convolving the PHANGS–MUSE internal release v2 to 15'' resolution, we match the astrometric grid of our IRAC- and WISE-based measurements to compare Σ_⋆ from MUSE to our estimate. We also calculate the ${{\rm{\Upsilon }}}_{\star }^{3.4}$ implied by MUSE by dividing the MUSE-based stellar mass by the WISE1 intensity.

We compare the MUSE results to our fiducial measurements in Figure 8. The left panel shows overall excellent agreement between our near-infrared-based Σ_⋆ and the MUSE-based values. Over two orders of magnitude in Σ_⋆, the ratio between the two methods remains roughly fixed, with only a modest offset and moderate scatter. Specifically, we find a median offset of ∼0.08 dex or a factor of ∼1.2, and ∼0.11 dex or 30% point-to-point scatter.

Zoom In Zoom Out Reset image size

The right panel of Figure 8 shows that the MUSE fitting results in a similar trend of ${{\rm{\Upsilon }}}_{\star }^{3.4}$ versus SFR/M_⋆ as what we adopt for the near-infrared data. The PHANGS–MUSE fitting implies a steady decline in ${{\rm{\Upsilon }}}_{\star }^{3.4}$ as a function of increasing SFR/M_⋆. The slope of the decline resembles that applied to the near-infrared maps (blue). The higher Σ_⋆ found in the right panel means that the actual ${{\rm{\Upsilon }}}_{\star }^{3.4}$ values implied are shifted to ∼20% higher values than those we predict based on Leroy et al. (2019). Mostly, the offset appears to be a simple multiplicative shift. The "reddest," lowest SFR/M_⋆ regions in the MUSE data also show higher ${{\rm{\Upsilon }}}_{\star }^{3.4}$, approaching 1 ${M}_{\odot }\,{L}_{\odot }^{-1}$, than implied by our adopted formula.

Overall, Figure 8 suggests excellent agreement between two independent ways to estimate the mass. The PHANGS–MUSE results provide a key validation of our adopted ${{\rm{\Upsilon }}}_{\star }^{3.4}$ treatment. To transfer from the "PHANGS near-infrared" approach used to compute sample properties to the PHANGS–MUSE system one would need to scale the near-infrared-based stellar mass or Σ_⋆ up by a factor of ∼1.2. Alternatively, one could scale the PHANGS–MUSE mass down by a similar factor.

Figure 9 and Table 5 provide further checks on our M_⋆ estimates. Our fiducial masses use IRAC or WISE near-infrared maps with the radially varying estimate of the mass-to-light ratio described above. In Figure 9 we adjust this approach and compare alternative estimates to our values. From left to right, we compare to (1) masses calculated only from integrated colors and luminosities by Leroy et al. (2019); (2) masses calculated using only IRAC or WISE data, i.e., checking whether adopting one data set or the other makes a difference; (3) masses calculated without any variable mass-to-light ratio, i.e., with a single fixed ${{\rm{\Upsilon }}}_{\star }^{3.4}$ of 0.35 ${M}_{\odot }\,{L}_{\odot }^{-1}$ applied to all data; and (4) masses calculated applying a fixed mass-to-light ratio of 0.5 ${M}_{\odot }\,{L}_{\odot }^{-1}$ to the contaminant-corrected "old stellar" ⁴⁷ maps of Querejeta et al. (2015). The final column, (5), compares the surface density, Σ_⋆, estimated in our adopted approach to that derived from PHANGS–MUSE spectral modeling. Unlike the other cases, which consider whole galaxies, this comparison treats each 15'' line of sight independently.

Zoom In Zoom Out Reset image size

Table 5. Comparison among Stellar Mass Estimates

Other Estimate	log10 Median a	log10 Scatter
	(dex)	(dex)
Integrated photometry	−0.07	0.06
WISE1 versus IRAC1	0.00	0.03
WISE1 + fixed ${{\rm{\Upsilon }}}_{\star }^{3.4}$	0.06	0.13
S4G ICA + fixed ${{\rm{\Upsilon }}}_{\star }^{3.4}$	0.00	0.11
PHANGS–MUSE b	0.08	0.11

Notes. "Integrated photometry": our fiducial estimate divided by those from Leroy et al. (2019) based only on integrated photometry. "WISE1 versus IRAC1": comparing results using WISE1 versus IRAC1 and otherwise treating the data identically. "WISE1 + fixed ${{\rm{\Upsilon }}}_{\star }^{3.4}$ ": using WISE1 but with only a fixed mass-to-light ratio of 0.35 ${M}_{\odot }\,{L}_{\odot }^{-1}$, not the radially varying or galaxy-dependent prescription used in our fiducial estimates. "S⁴G ICA + fixed ${{\rm{\Upsilon }}}_{\star }^{3.4}$ ": "old stellar emission" maps from Querejeta et al. (2015) using a fixed ${{\rm{\Upsilon }}}_{\star }^{3.4}$ of 0.5 ${M}_{\odot }\,{L}_{\odot }^{-1}$. The table reports comparisons for a sample of local galaxies that have CO mapping including all PHANGS–ALMA targets. The results restricting to just PHANGS–ALMA are almost identical. See also Figure 9.

^a Reported as log₁₀ of the median of the other estimate over our estimate for galaxies where both are available. ^b For PHANGS–MUSE only the statistics refer to a line-of-sight by line-of-sight comparison at 15'' resolution. For all other cases we refer to comparison of integrated M_⋆.

Download table as:  ASCIITypeset image

Table 5 summarizes the numerical results. Qualitatively, we see that shifting from the estimates based on integrated photometry from Leroy et al. (2019) to resolved estimates of ${{\rm{\Upsilon }}}_{\star }^{3.4}$ here increases the overall stellar mass by only ∼20% on average but can matter much more in a few cases. These tend to be cases with strong nuclear concentrations, but also a few cases where our more careful treatment of the data here has improved the masking. Besides a few edge cases, it makes a negligible difference to the overall M_⋆ whether we use IRAC or WISE as long as the treatment of the data remains the same. The adoption of variable ${{\rm{\Upsilon }}}_{\star }^{3.4}$ does matter. Without the variable ${{\rm{\Upsilon }}}_{\star }^{3.4}$, masses scatter by ±0.13 dex or 35% relative to our derived stellar masses. Similarly, if we integrate the contaminant-corrected "old stellar emission" component maps from Querejeta et al. (2015) and apply a fixed ${{\rm{\Upsilon }}}_{\star }^{3.4}$, we find ±0.11 dex or 30% scatter relative to our adopted values.

In Figure 9 and Table 5 we have already chosen our typical mass-to-light ratios to remove systematic offsets between our adopted approach and the comparison data. Thus two key points from the comparison are:

• 1.  
  We find ${{\rm{\Upsilon }}}_{\star }^{3.4}=0.35$ ${M}_{\odot }\,{L}_{\odot }^{-1}$ on average.
• 2.  
  When using the "old stellar emission" maps derived by Querejeta et al. (2015) from combining IRAC 3.6 and 4.5 μm emission, we suggest to apply ${{\rm{\Upsilon }}}_{\star }^{3.4}=0.5$ ${M}_{\odot }\,{L}_{\odot }^{-1}$ for consistency with our adopted system. The difference between this value and the average ${{\rm{\Upsilon }}}_{\star }^{3.4}=0.35$ ${M}_{\odot }\,{L}_{\odot }^{-1}$ reflects the removal of contaminants from the Querejeta et al. (2015) maps.

Although we set the median to agree, the scatter between our adopted system and either the Querejeta et al. (2015) or fixed ${{\rm{\Upsilon }}}_{\star }^{3.4}$ is still significant. In general, our checks show that any substantial shift in methodology induces ∼0.1 dex, or ≈25%, scatter in the measurements. This uncertainty dramatically exceeds any statistical uncertainty associated with the photometry and we suggest adopting ±0.1 dex as a realistic systematic uncertainty on stellar masses of individual galaxies (in good agreement with the level of uncertainty due to age variations found in the IRAC-based ICA maps by Meidt et al. 2014). It also seems reasonable to adopt a comparable uncertainty to describe the overall calibration of our stellar mass scale.

Despite these uncertainties, we emphasize again that (1) our calibration scheme anchors our mass scale to the Salim et al. (2016, 2018) and Leroy et al. (2019) estimates, and through these to work on the SDSS main galaxy sample, and (2) that our adopted ${{\rm{\Upsilon }}}_{\star }^{3.4}$ treatment seems to qualitatively agree with results from PHANGS–MUSE, while showing a mild 0.08 dex systematic offset.

We measure R_e and l_⋆ from the 75 resolution estimates of Σ_⋆ described in Section 4.3. Using these data, we construct a radial profile using the orientations in Table 3. In order to lower our sensitivity to the adopted inclination and resolution, we consider only points within ±75° of the major axis, except for the central 30''. Within the central 30'', we include all position angles. ⁴⁸ Then, we measure the azimuthally averaged intensity using the mean within 30'' and the median intensity in each ring at larger radii. The use of the median suppresses any unmasked foreground stars at the expense of also losing some sensitivity to galactic structure. We use the mean in the central aperture to retain full sensitivity to any bright nuclear emission when calculating R_e. Next, we subtract a background value from the radial profile, which is just the median of all radial bins outside 2 r₂₅. Finally, we integrate these profiles assuming azimuthal symmetry and use this integral to identify the half-mass radius, R_e, the radius encompassing 90% of the mass, R₉₀, and the stellar scale length, l_⋆. The integral for calculating R_e extends out to 2 r₂₅, which we verified by eye in each case to represent an apparent convergence of the radial profile. We report the estimated R_e and l_⋆ in Table 3.

To fit l_⋆, we identify a fitting range over which the radial profile appears exponential. By default, this range spans from 0.4 to 1.0 r₂₅, i.e., excluding roughly the central two scale lengths. We examined each profile by eye and adjusted this range to reflect where the disk appears exponential over a large range of radii outside the galaxy center. This approach does a good job of matching sophisticated disk/bulge decompositions (e.g., Salo et al. 2015). This also closely resembles the approach adopted by Sánchez et al. (2014), a study that we use to estimate metallicity gradients.

We measure effective radii and stellar scale lengths. The 25th magnitude B-band isophotal radius also remains widely used, particularly for studies of local galaxies. In Table 6 and Figure 10 we report translations between different size measurements for our galaxies: half-mass radius, R_e, profile-based scale length, l_⋆, and isophotal radius, r₂₅ (de Vaucouleurs et al. 1991; Makarov et al. 2014, we adopt from RC3 via HyperLEDA). For reference, for a pure exponential disk, l_⋆ should relate to the effective radius via

$$\begin{eqnarray}&&{R}_{{\rm{e}}}^{{\rm{d}}{\rm{i}}{\rm{s}}{\rm{k}}}=1.68\,{l}_{\star }\,{\rm{o}}{\rm{r}}\,{l}_{\star }=0.60\,{R}_{{\rm{e}}}^{{\rm{d}}{\rm{i}}{\rm{s}}{\rm{k}}}.\end{eqnarray} \tag{ 4 }$$

In practice, we find that our measured R_e and l_⋆ deviate from this ideal case. We show this in Figure 11 and report the typical ratio and scatter for our measurements in Table 6. We find that in practice, the relationship in our sample can be better described via

$$\begin{eqnarray}&&{R}_{{\rm{e}}}=1.41\,{l}_{\star }\,{\rm{o}}{\rm{r}}\,{l}_{\star }=0.71\,{R}_{{\rm{e}}}.\end{eqnarray} \tag{ 5 }$$

The difference between these two formulae simply reflects that our targets are not ideal exponential disks, and that many galaxies harbor compact inner structures. These drive R_e to smaller values but do not affect l_⋆. As Figure 11 shows, even Equation (5) still misses a significant minority of outliers that have high l_⋆ compared to R_e, reflecting strong central mass concentrations relative to an exponential disk.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The 25th magnitude B-band isophotal radius, r₂₅, remains in common use when studying nearby galaxies. Table 6 and Figure 10 show that r₂₅ exhibits a well-defined scaling with our two size measures:

$$\begin{eqnarray}&&{R}_{{\rm{e}}}\approx {r}_{25}/3.08\,{\rm{a}}{\rm{n}}{\rm{d}}\,{l}_{\star }\approx {r}_{25}/4.3\end{eqnarray} \tag{ 6 }$$

on average, both with ∼0.14 dex scatter.

We compare our estimated sizes to literature size estimates for the same galaxies and show the results in Table 6 and Figure 10. We compare our measured sizes to:

• 1.  
  The 25th magnitude B-band isophotal radius from RC3 via HyperLEDA (de Vaucouleurs et al. 1991; Makarov et al. 2014). This remains the most widely used measure of galaxy size for studies of local galaxies, and the results of this comparison have already been used to derive Equation (6).
• 2.  
  Effective radii estimated by the 2MASS Large Galaxy Atlas (LGA; Jarrett et al. 2003). They measured R_e at J, H, and K bands for many of our galaxies and a large sample of additional galaxies. These are profile-based sizes using near-infrared data and we expect our measurements to match them well.
• 3.  
  Effective radii estimated from the S⁴G near-infrared data by Muñoz-Mateos et al. (2015). These are also profile-based measurements, in many cases using the same data with slightly distinct processing, assumptions, and orientations. We expect to match these values well.
• 4.  
  The results of single- and multi-component fits to the S⁴G data by Salo et al. (2015). They present one-component Sérsic profiles, which yield an effective radius. They also fit exponential disk components to most targets. When they fit an expdisk component, their h _r captures the exponential scale length and represents the most direct comparison for our l_⋆.

Overall, Figure 10 and Table 6 show that our measured sizes match measurements from other sources well. As expected, we match previous profile-based estimates of the size using near-infrared data well. In 2MASS, we find best agreement with the J-band size, with a median ratio of 2MASS-to-PHANGS R_e of 0.99. We find slightly smaller sizes in 2MASS H and K bands, likely reflecting the poorer sensitivity in those bands, but the agreement is also excellent for those bands. Our match to the Muñoz-Mateos et al. (2015) sizes is also good; we find an 0.95 ratio on average.

On average, the effective radii from the Sérsic model fits by Salo et al. (2015) are ∼21% larger than our R_e. This likely reflects the fact that fitting a single Sérsic profile to the entire galaxy de-emphasizes the contribution of any distinct, massive central component. By-eye checks of galaxies with strongly discrepant size estimates between us and Salo et al. (2015) reveal that such galaxies often have central mass concentrations. These same discrepancies arise comparing the Salo et al. (2015) values to Muñoz-Mateos et al. (2015) ones. Those two measurements use the same data, but the single-component fitting approach of Salo et al. (2015) tends to yield larger R_e than the profile fits of Muñoz-Mateos et al. (2015), a 26% higher R_e for matched galaxies, on average, with ∼0.1 dex scatter.

Despite the good match on average, we find just under ∼0.1 dex or about 25% scatter between our measurements and those of Muñoz-Mateos et al. (2015) or Jarrett et al. (2013), and we adopt ±0.1 dex as a realistic 1σ uncertainty on R_e. The scatter between estimates reflects differences in the adopted orientations and sensitivity to many aspects of the exact method used. The method of defining the total flux, choice to focus on the major axis or whole galaxy, choice to apply a radially varying mass-to-light ratio, and a host of other minor decisions all affect the derived R_e. To some extent, this reflects that galaxies are not axisymmetric, thin, tilted disks with a single orientation and a simple radial profile. It also highlights how R_e , which depends on an estimate of the total flux and focuses on the inner, often complex, part of the galaxy is a less stable size measure than an isophotal radius for local galaxies.

Figure 10 and Table 6 show that our stellar scale lengths agree reasonably well with those from Salo et al. (2015). Their sophisticated two-dimensional decomposition tends to yield ∼11% larger scale lengths than our fits to the radial profiles. Our approach focuses on fitting the part of the profile roughly in the range 0.4−1.0 r₂₅ that looks approximately exponential. This might lead us to emphasize more cleanly declining parts of the disk and avoid regions with flatter profiles. Regardless, the overall agreement with Salo et al. (2015) is also good, and l_⋆ is likely somewhat more uncertain than R_e.

Tables 7 and 8 and Figure 12 illustrate the uncertainty in our estimated SFR. Table 7 reports the scatter and bias implied by adopting a single linear calibration and attempting to match the full Salim et al. (2018) data set. Here "bias" refers to the rms scatter in the implied calibration factor, C, as a function of M_⋆ and SFR. To calculate this, Leroy et al. (2019) binned the Salim et al. (2018) data into cells of SFR and M_⋆, solved for C in each cell, and then calculated the rms scatter across all populated cells in the SFR–M_⋆ plane. Thus, this number captures the degree of systematic uncertainty implied by applying a single value of C to the whole galaxy population. This is discussed in detail by Leroy et al. (2019), but here we emphasize two key points. First, using WISE4 only shows notably higher bias than WISE4 combined with a UV band. This reflects the fact that low-mass galaxies, in particular, often show less obscuration compared to high-mass galaxies. The importance of a hybrid approach to star formation rate tracers has been reviewed by Kennicutt & Evans (2012) and Calzetti (2013).

Zoom In Zoom Out Reset image size

Table 7. Adopted SFR Calibrations and Uncertainties

Band	${\mathrm{log}}_{10}C$ a	Rms Scatter b	Rms Bias c
	$\left(\displaystyle \frac{{M}_{\odot }\,{\mathrm{yr}}^{-1}}{\mathrm{erg}\,{{\rm{s}}}^{-1}}\right)$	(dex)	(dex)
FUV d	−43.42	0.1	0.09 e
NUV d	−43.24	0.1	... e
WISE4 with FUV	−42.73	0.17	0.06
WISE4 with NUV	−42.79	0.18	0.08
WISE4 alone	−42.63	0.17	0.15
WISE3 with FUV	−42.79	0.21	0.44
WISE3 with NUV	−42.86	0.22	0.47
WISE3 alone	−42.70	0.2	0.37

Notes. From Leroy et al. (2019) based on the GSWLC (Salim et al. 2016, 2018). These coefficients represent a best effort to anchor our SFR estimates to both a large population of local galaxies (Leroy et al. 2019) and the very large SDSS-based catalog of Salim et al. (2016, 2018).

^a C is the factor to convert ν L_ν to SFR (notation follows Kennicutt & Evans 2012). ^b Galaxy-to-galaxy scatter of C, in dex, from comparison to Salim et al. (2018). See Leroy et al. (2019) for more details. ^c "Bias" here refers to the systematic scatter, in dex, among the empirical calibration coefficients to be applied to the WISE data across the SFR–M_⋆ plane. This is calculated by binning the Salim et al. (2018) data into cells of SFR and M_⋆, calculating C for each cell, and then measuring the rms scatter across the binned data. This indicates the amount of systematic uncertainty in the empirical calibration of the tracer. See Leroy et al. (2019) for more details. ^d The coefficients applied to FUV and NUV do not depend on which IR band is used. ^e The bias in the calibration coefficient to convert extinction-corrected FUV emission from the Salim et al. (2018) fits to SFR for cells in SFR–M_⋆ space with at least 50 galaxies. The NUV calibration factor is extrapolated from FUV, so no independent bias is calculated.

Download table as:  ASCIITypeset image

The other important point from Table 7 relates to the use of WISE3. WISE3 data have better sensitivity and resolution than the WISE4 data, and below we will use them to aperture correct our CO data. However, as Table 7 shows, the appropriate calibration to translate WISE3 to SFR shows strong trends as a function of stellar mass and specific star formation rate. This is consistent with the significant contribution of polycyclic aromatic hydrocarbons (PAHs) to the WISE3 band and the known dependence of PAH emission on metallicity (e.g., Engelbracht et al. 2006) and radiation field (e.g., Chastenet et al. 2019). This is a strong effect, greater than a factor of 2. We do note that this will mainly affect extremes in the galaxy population, i.e., low-mass, low-metallicity galaxies and very low and high SFR/M_⋆ galaxies.

For some PHANGS–ALMA science applications, the higher resolution and improved signal-to-noise ratios of the WISE3 data will be important. Indeed, we used WISE3 to define the field of view for exactly these reasons. When we need such a ratio to be quantitative, we will generally estimate a WISE3-to-WISE4 ratio for that galaxy or location and use this as an empirical correction to apply the better-calibrated WISE4-based SFR prescriptions.

Figure 12 and Table 8 mostly bear out these trends within our own data. In this figure we vary the adopted band combination and compare to our fiducial SFR estimate. We also compare the value from Leroy et al. (2019) based only on integrated photometry. Both the table and the figure show that all band combinations have consistent calibrations on average. However, if we do not hybridize WISE with some UV band, then we find a higher galaxy-to-galaxy scatter, ∼0.1 dex using only WISE4. If we substitute WISE3 for WISE4, we also see ∼0.1 dex scatter relative to our fiducial estimates. Finally, while using WISE3 alone yields a matched result on average, though, we find almost 0.2 dex galaxy-to-galaxy scatter relative to our fiducial estimates.

For 19 of our targets, PHANGS–MUSE (E. Emsellem et al. 2021, in preparation) offers an independent, high-quality estimate of Σ_SFR from Hα recombination line emission corrected for extinction using the Balmer decrement. An in-depth exploration of the differences between these two methods is presented by F. Belfiore et al. (2021, in preparation). Here, we present a brief comparison of Σ_SFR estimated from the PHANGS–MUSE Hα and Hβ maps at 15'' resolution and the UV+IR Σ_SFR estimates available for the whole sample.

For this comparison, we compare our best Σ_SFR estimates, which are from FUV+WISE4 for all 19 MUSE targets, to Σ_SFR estimated from Hα and Hβ measured from PHANGS–MUSE at 15'' resolution. We use the two recombination lines to estimate A_Hα , the extinction affecting the Hα line, by assuming a screen geometry, a Cardelli et al. (1989) extinction curve, and Case B recombination (see F. Belfiore et al. 2021, in preparation for more details). Then we place the two 15'' resolution Σ_SFR maps on the same astrometric grids and compare results.

Figure 13 shows the results of this comparison. As the left panel shows, the two estimates correlate extremely well, showing a Spearman rank correlation coefficient of ≈0.95 despite using completely different approaches to estimate Σ_SFR. However, the estimators do yield numerically different estimates. On average the FUV+WISE4 estimate is 20%–30% higher than the Hα-based estimate. Specifically, treating all points equally yields a median Hα-to-FUV+WISE4 ratio of ≈0.73, or −0.14 dex, with ≈0.2 dex, or 60% rms scatter about the median. Considering the binned measurements, shown in blue, we find a median ratio of 0.84 or −0.07 dex with ≈0.1 dex rms scatter in the median value from bin to bin. Given the independent approaches to estimate Σ_SFR, this offset still represents reasonable agreement, and any average offset could be accounted for by adjusting the empirically derived coefficient on the WISE4 term in the FUV+WISE4 indicator.

Zoom In Zoom Out Reset image size

The Figure shows that the offset and scatter between the two tracers do reflect an average offset and a systematic trend, such that the two tracers agree well at Σ_SFR ≳ 10⁻² M_⊙ yr⁻¹ kpc² but show increasing divergence below this value. Specifically, at low Σ_SFR, the FUV+WISE4 SFR indicator predicts a higher Σ_SFR than the Hα-based estimate. Figure 13 illustrates a systematic trend with Σ_SFR, and F. Belfiore et al. (2021, in preparation) show that a similar trend holds as a function of the ratio Σ_SFR/Σ_⋆, a quantity that traces the local specific star formation rate.

The sense of the trend agrees with the expectation that the 22 μm emission captured in the WISE4 band may contain significant emission not directly associated with star formation, and that this contamination will be stronger in regions with low star formation activity (e.g., see Groves et al. 2012; Leroy et al. 2012; Boquien et al. 2016; Simonian & Martini 2017). Given the multiwavelength coverage in PHANGS, it may be possible to correct these low-intensity sightlines using physical or empirical prescriptions for the IR cirrus (e.g., Leroy et al. 2012; Davis et al. 2014). Alternatively, high-resolution mid-infrared observations from the James Webb Space Telescope offer a path forward to morphologically and structurally isolate emission associated with star-forming regions (e.g., Calzetti et al. 2005; Boquien et al. 2016).

We refer the reader to F. Belfiore et al. (2021, in preparation) for more details but emphasize three main points from the comparison. First, the overall scaling between the two tracers appears very good, and the median offset of 20%–30% represents good agreement given the independent approaches. Second, the clear systematic trends suggest that the FUV+WISE4 maps be treated with caution at low Σ_SFR ∼ 10⁻³ M_⊙ yr⁻¹ kpc². A cirrus correction to the WISE4 data might improve the situation. However, third, the coefficients on the FUV+WISE4 estimator have been set to ensure average consistency with the GSWLC and the larger SDSS sample (Salim et al. 2016, 2018). For this work, we use the coefficients "as is" and note the comparison to PHANGS–MUSE as an important indicator of the systematic uncertainty.

Our maps only cover the area of active star formation, defined by an IR intensity contour (Section 5). Though we expect most CO emission to lie within our target area, faint CO emission can continue well into the outer disks of galaxies (e.g., Young et al. 1995; Braine et al. 2007; Schruba et al. 2011). To account for this, we derive an aperture correction that can be applied to our measured integrated CO(2−1) to yield the global L_CO for the galaxy.

After considering several options, we chose to calculate this aperture correction by using the WISE3 map as a template for CO emission. We calculate the sum of WISE3 emission over the whole galaxy and the sum of WISE3 emission over just the region with CO coverage. In Table 4 we report this WISE3-based aperture correction, but note that we leave it to the user to actually apply the correction. The L_CO that we report reflects only the direct integration of our maps.

Table 9 and Figure 14 show the motivation for choosing WISE3. Inside the PHANGS–ALMA footprint, we construct several possible "templates" that could be used for the aperture correction. Here, a template refers to a quantity that might linearly trace CO intensity. We consider WISE3 intensity, WISE4 intensity, Σ_SFR from FUV+WISE4, and Σ_⋆ estimated from WISE1+${{\rm{\Upsilon }}}_{\star }^{3.4}$. We also consider three plausible versions of an exponential disk model: one adopting a scale length of 0.2 r₂₅ (Young et al. 1995; Leroy et al. 2008; Bolatto et al. 2017), one adopting our fit exponential scale length, and one adopting 0.6 R_e (see Equation (4)). These size-based corrections resemble those used, e.g., by Lisenfeld et al. (2011) and Saintonge et al. (2011).

Zoom In Zoom Out Reset image size

Table 9. PHANGS–ALMA CO(2−1) Aperture Corrections

Template	Correction a	Correlation b
WISE3 (adopted)	1.25 (1.10–1.57)	0.82
WISE4	1.16 (1.04–1.59)	0.77
SFR (FUV+W4) c	1.19 (1.07–1.84)	0.79
WISE1+${{\rm{\Upsilon }}}_{\star }^{3.4}$	1.41 (1.19–1.86)	0.53
Exponential disk (0.2 r25)	1.16 (1.07–1.55)	0.43
Exponential disk (l⋆)	1.36 (1.10–1.90)	0.41
Exponential disk (0.6 Re)	1.18 (1.05–1.47)	0.18

Notes. This table reports results only for PHANGS–ALMA. The correction gives the implied aperture corrections. Correlation indicates how well each band linearly predicts CO(2−1) emission for each individual galaxy inside the coverage.

^a Median and 16th to 84th percentile range for the aperture correction to be applied to the PHANGS–ALMA data. ^b Linear correlation coefficient relating the template band to CO intensity in the region with coverage. Each galaxy is first normalized by the median template-to-CO ratio. ^c Best SFR estimate. This uses FUV+WISE4 when available. Otherwise we adopt NUV+WISE4 if available and WISE4 otherwise.

Download table as:  ASCIITypeset image

Table 9 shows that WISE3 exhibits the strongest linear correlation coefficient with CO intensity of any considered quantity. After scaling the data for each galaxy by the median WISE3-to-CO ratio for that galaxy, the WISE3 emission exhibits a linear correlation coefficient of 0.82 with CO(2−1) intensity across the PHANGS sample. This stunningly strong correlation has also been seen and investigated in detail by Chown et al. (2021) in the EDGE–CALIFA CO(1−0) data. We confirm their finding that I_CO shows a strong, linear correlation with WISE3 within galaxies. After removing a galaxy-to-galaxy offset, the correlation between WISE3 emission and CO(2−1) intensity is even stronger than that relating CO to Σ_SFR estimated from FUV and WISE4 emission. Figure 14 shows this correlation. After galaxy-by-galaxy normalization, WISE3 tracks CO(2−1) emission almost linearly with <0.1 dex scatter across three decades in intensity. We do caution that, as discussed above, the ratio of WISE3 emission to SFR shows strong galaxy-to-galaxy variations. These variations correlate with metallicity and specific star formation rate (Leroy et al. 2019). Our application avoids much of this issue by removing an overall galaxy-to-galaxy scaling. The PHANGS–ALMA data offer the possibility to explore the empirical and physical nature of the WISE3–CO correlation in much more detail in future works.

WISE3 offers an outstanding option to derive aperture corrections and we derive a median correction of ∼1.25 for PHANGS–ALMA, with most targets yielding corrections between 1.1 and 1.5. In other words, the WISE3 emission implies that our maps usually cover ∼70%–90% of the total CO emission. Table 9 and the right panel of Figure 14 show that we would have arrived at similar aperture corrections using only WISE4 emission or our FUV+WISE4-based Σ_SFR estimates. Our stellar mass maps imply slightly higher corrections, on average, but Σ_⋆ also correlates notably less well with I_CO than the SFR and mid-infrared tracers. The right panel in Figure 14 shows that had we adopted one of these other plausible templates, we might expect the aperture correction to scatter by ±10% (the gray band) relative to our preferred WISE3 values. We consider that the true uncertainty on the aperture correction is slightly less than this because of the strength of the WISE3–CO correlation.

Aperture corrections that assume an exponential disk structure have the large advantage that they can be applied with almost no multiwavelength data. However, simple exponential disk models show a notably worse correlation with I_CO than any of the other templates. Among these, the template using our fitted l_⋆ and the one using a scale length of 0.2 r₂₅ offer approximately equally good corrections. On average, these disk models yield about the same magnitude of correction factor as our other templates. However, Table 9 implies that they offer less precision, and Figure 14 shows that the derived correction using an l_⋆ scale length exponential disk scatters by 10%–20% relative to the WISE3 correction.

Figure 15 compares our derived CO fluxes and intensities to previous CO(2−1) mapping using the IRAM 30 m telescope. ⁵⁰ In total, PHANGS–ALMA covers 18 galaxies that have also been mapped in the CO(2−1) line using the HERA array receiver on the IRAM 30 m telescope (Schuster et al. 2007). These maps are part of the HERACLES survey (Leroy et al. 2009, 2013b) and a follow-up survey (A. Schruba et al. 2022, in preparation). They cover a wider area than the ALMA maps, but with higher noise levels. The calibration of the IRAM 30 m maps depends on "chopper wheel" observations and a main beam efficiency correction, while both the ALMA interferometric and total power data are pinned to the monitoring of ALMA calibrators with monitored fluxes.

Zoom In Zoom Out Reset image size

The IRAM 30 m maps have native resolution of 133, so we must smooth the PHANGS–ALMA data make a direct comparison. For this exercise, we used a 17'' version of both cubes. We aligned both cubes for each galaxy onto the same astrometric and spectroscopic grid. We then apply an identical mask to each galaxy in each survey to create matched resolution, matched velocity-interval integrated intensity maps, along with associated uncertainty maps.

The top left panel in Figure 15 compares the integrated fluxes between the IRAM and ALMA maps. For the calculation in the top left panel, we apply the WISE3-based aperture corrections discussed above. We find a median PHANGS-to-HERA ratio of 0.97, or −0.01 dex, with a robustly estimated galaxy-to-galaxy scatter of 0.06 dex, or about 15%. The right panel shows galaxy-integrated ratios. We calculate these only integrating over the area where both maps have coverage. These also show good agreement. The statistical uncertainty in these ratios is low, <1% in almost all cases, so the observed scatter reflects mostly calibration or other systematic uncertainties. The 16th to 84th percentile range, illustrated in gray, spans from 0.87 (−0.06 dex) to 1.06 (0.025 dex) with the median 0.97 (−0.01 dex). These suggest an overall ≲5% difference in the calibration scale with about ±10%–15% calibration uncertainty in any given galaxy. We mark four galaxies with strongly discrepant fluxes in Figure 15: NGC 3351 (PHANGS-to-HERA ratio of 0.74), NGC 3627 (1.24), NGC 4579 (0.80), and NGC 4689 (1.32). The discrepancy in NGC 3627, which was observed by the IRAM 30 m under nonideal conditions shortly after HERA commissioning, has been noted and explored by den Brok et al. (2021). They also noted a similar discrepancy between new and old IRAM 30 m mapping projects for NGC 5194 (M51), which is not in our sample. den Brok et al. suggest that issues with the HERA data most likely account for these strong outliers.

The bottom rows in Figure 15 explore the pixel-to-pixel correspondence between the PHANGS and HERA data at 17'' resolution. Here, we only consider where both maps have coverage and I_CO > 1 K km s^{−1 51} in both maps, a total of 1161 independent 17'' lines of sight. Again, we find good overall agreement, here spanning two decades in CO intensity. The bottom right panel shows the distribution of ratios for individual pixels. Here, we find a median ratio of 0.93 (−0.03 dex) and a robustly estimated scatter of ±0.11 dex or about ±28%. Thus, the pixel data indicate a slight offset in overall flux scale, now with PHANGS ∼7% fainter than HERA.

Based only on the statistical errors associated with individual data, we expect a scatter of about ±0.05 dex or about ±12%, a range that we indicate by a blue bar in Figure 15. This is significantly smaller than the measured scatter, again indicating the presence of systematic uncertainties. The blue histogram shows a model that treats these systematic uncertainties as a ±20% point-to-point multiplicative uncertainty and also incorporates the measured statistical errors. Overall this matches the measured scatter well.

This ±20%, or ±0.11 dex, point-to-point scatter is somewhat higher than the galaxy-to-galaxy scatter seen in the top row. Indeed, if we subtract the galaxy-to-galaxy scatter and the statistical scatter from the measured scatter in quadrature we are left with about ±0.06–0.07 dex (about ±15%) point-to-point unaccounted scatter. Our best estimate is that this reflects uncertainties in the image reconstruction in both data sets. The HERA array receiver data suffer from variable pixel gains (see Leroy et al. 2009; den Brok et al. 2021), while the ALMA data have been partially reconstructed from interferometric imaging. Based on a detailed analysis and visual inspection of the five galaxies in den Brok et al. (2021), the HERA pixel gain uncertainties appear to be responsible for the most serious issues, and we primarily attribute this unaccounted scatter to the HERA data. This will appear, e.g., as mapping artifacts like striping or variations in the gain across the map. Such issues are common in single-dish maps, especially those made with array receivers, and the identified magnitude of uncertainty, ±15%, agrees with that inferred from jackknife tests carried out on the HERACLES data (Leroy et al. 2009).

Overall, this comparison paints the following picture:

• 1.  
  The two data sets agree very well to first order.
• 2.  
  There appears to be a ∼3%–7% systematic offset in the flux scale between PHANGS and the HERA data, with the PHANGS data showing slightly lower flux. Our best estimate is that the ALMA data are better calibrated.
• 3.  
  We observe ∼15% galaxy-to-galaxy scatter in the ratio of ALMA-to-HERA fluxes, which reflects the combined calibration uncertainty in the ALMA and IRAM 30 m data. This is largely consistent with the 10%–15% calibration uncertainty expected for the HERA data (Leroy et al. 2009; den Brok et al. 2021) and the 5%–10% uncertainty associated with ALMA band 6 observations.
• 4.  
  At 17'' resolution, the point-to-point scatter between the two data sets is ∼0.11 dex. This is higher than expected from only the statistical noise or the galaxy-to-galaxy flux calibration uncertainties. We consider that this 0.11 dex scatter reflects roughly equal contributions from statistical uncertainty, flux calibration uncertainty, and point-to-point gain variations in the single-dish map. The latter arises from pixel-to-pixel gain variations in the array receiver (see den Brok et al. 2021) and uncertainties in interferometric image reconstruction.

We target the J = 2 → 1 rotational transition of CO, CO(2−1). This line is readily excited at temperatures and densities that characterize molecular clouds, and its spatial distribution in nearby galaxies correlates well with the fundamental CO(1−0) transition (e.g., see Braine & Combes 1992; Leroy et al. 2009, 2013b; Koda et al. 2020; den Brok et al. 2021). This choice represents a key practical element of our observing strategy. Given typical observed CO(2−1)-to-CO(1−0) line ratios of R₂₁ ≈ 0.5−1 (Leroy et al. 2013b; den Brok et al. 2021) and considering the array configurations and system temperatures at ALMA, mapping CO(2−1) is a factor of ∼2–4 times faster than mapping CO(1−0) to the same mass surface density sensitivity and angular resolution. Moreover, the ν = 230.538 GHz rest frequency of ¹²CO(2−1) lies in a favorable part of the atmosphere for ALMA, and can be observed effectively almost any time that ALMA is on sky.

We visualize the spectral setup of our observations in Figure 19. We observed CO(2−1) with a bandwidth of 937.5 MHz and a channel width of Δν = 244 kHz. After taking into account the online Hanning smoothing, this corresponds to a native channel width of ∼0.6 km s⁻¹ and a total velocity coverage of ∼1200 km s⁻¹. This channel width easily resolves emission lines from individual molecular clouds, which have a typical full line width of ∼10 km s⁻¹. The velocity coverage easily encompasses the full emission from each of our targets. In practice, we image the cubes at a velocity resolution of ∼2.5 km s⁻¹ to reduce data volume and improve the signal-to-noise ratio and quality of the deconvolution.

Zoom In Zoom Out Reset image size

We allocated the other correlator resources in several different ways over the course of the project. During the Large Program and the two main pilot projects (see Table 2), we configured the spectral setup to observe the millimeter continuum with the other three windows. We configured each to have the maximum possible 1.875 GHz bandwidth. We placed one of these windows to cover the C¹⁸O(2–1) line at a spectral resolution of 2 MHz, corresponding to ∼5–6 km s⁻¹ after taking the Hanning smoothing into account. The C¹⁸O isotopologue is several hundred times less abundant than ¹²C¹⁶O (e.g., Wilson & Rood 1994) and the C¹⁸O(1–0) line is often 30–100 times fainter than CO(1−0) in local galaxies (e.g., Jiménez-Donaire et al. 2019). Therefore the C¹⁸O line represents a target of opportunity that may be recovered in a few bright regions or via stacking. Though we imaged these data, this line is not included in the first PHANGS–ALMA data release.

In one early set of observations targeting NGC 0628, also illustrated in Figure 19, we targeted CS(5–4) instead of C¹⁸O(2–1). The main practical difference is that this places the main CO(2−1) line in the lower sideband instead of the upper sideband, which affects the spectral variation of the noise (see Leroy et al. 2021b).

ALMA began allowing more flexible Band 6 spectral setups in Cycle 6 after the execution of most PHANGS–ALMA observations. Therefore, in the cases of some of the PHANGS–ALMA extensions using the ACA 7 m array, we also targeted the ¹³CO(2−1) line. The ¹³CO isotopologue is only ∼30–120 times less abundant than ¹²CO (e.g., Wilson & Rood 1994; Milam et al. 2005), and the ¹³CO line is only ∼10 times fainter than ¹²CO (e.g., Roman-Duval et al. 2016; Cormier et al. 2018). This makes the ¹²CO/¹³CO ratio a useful probe of optical depth and excitation (e.g., Cormier et al. 2018). As shown in Figure 19, observing ¹³CO and C¹⁸O together requires using two somewhat narrower 937.5 MHz windows.

Finally, we note that ALMA employs "Doppler setting" rather than online Doppler tracking and records data in the topocentric frame. Our observed bandwidth is large enough that there is never any problem keeping the line emission from the galaxy in the full spectral window. However, this does imply that we must employ spectral regridding to construct the final velocity cubes.

We target a FWHM resolution of 1'', which corresponds to ∼100 pc resolution at 20 Mpc or ∼80 pc at the outer edge of our selection function (see Section 3). We picked this resolution so that an individual beam matches the size of an individual massive GMC or association of massive molecular clouds. During Cycle 5, when the PHANGS–ALMA Large Program was approved, this resolution at Band 6 and ν ≈ 230 GHz corresponded to the C43-2 array configuration. In practice, because ALMA executes observations within some tolerance, many of our observations also occurred in configurations similar to C43-1 or C43-3. Several factors related to data processing mean that our final data products have typical angular and physical resolutions slightly coarser than 1'' (see Section 5.2).

Our galaxies typically have full optical sizes of several arcminutes (Section 4). CO emission often shows extended, complex structure across the full area of a galaxy. The C43-2 main configuration recovers emission only from scales of ≲10''. Thus short- and zero-spacing information is important to accurately reconstruct the intensity distribution (e.g., see the case of M51; Koda et al. 2009; Pety et al. 2013; and see details for PHANGS–ALMA data processing in Leroy et al. 2021b). Therefore, we observed all targets with both components of the Morita Atacama Compact Array (ACA), the compact array of 7 m dishes and the 12 m single-dish telescopes used for "total power."

The ratio of observing time among the main array, 7 m array, and total power antennas followed the standard observatory recommendations. For our fiducial configuration C43-2, the ALMA Cycle 5 technical handbook recommended an observing time ratio of 1:5:8.5 among the 12 m main array, the 7 m array, and the total power antennas.

In extension projects, we observed a subset of very nearby targets using only the ACA 7 m and total power arrays. ⁵² These projects targeted galaxies within 5 Mpc, and typically at ∼3–4 Mpc, so that the ∼7'' resolution of the 7 m array at ν = 230 GHz already corresponds to 100–140 pc, reasonably matched to the main array resolutions for more distant targets.

PHANGS–ALMA integrates for 20–30 s per pointing, quickly covering each field in each large mosaic with a modest integration time. To arrive at this strategy, we targeted a fiducial rms noise level of 7.5 mJy beam⁻¹ or better at the frequency of the redshifted CO(2−1) line and a 5 km s⁻¹ channel width. At our nominal 1'' target resolution, this target sensitivity corresponds to an rms brightness temperature sensitivity of ∼0.17 K.

We can also express our target sensitivity as a mass surface density sensitivity via:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{mol}}={\alpha }_{\mathrm{CO}}^{1-0}\,{R}_{21}^{-1}\,{I}_{\mathrm{CO}\ (2-1)}\,\cos i.\end{eqnarray} \tag{ 10 }$$

Here, ${\alpha }_{\mathrm{CO}}^{1-0}$ is the CO(1−0) conversion factor, R₂₁ is the CO(2−1)-to-CO(1−0) line ratio in kelvin, i is the inclination of the galaxy, and I_{CO (2−1)} is the line-integrated CO(2−1) intensity in K km s⁻¹. For a standard Milky Way ${\alpha }_{\mathrm{CO}}^{1-0}=4.35$ M_⊙ pc⁻² (K km s⁻¹)⁻¹ (Bolatto et al. 2013a), a CO(2−1)-to-CO(1−0) ratio of R₂₁ = 0.65 (Leroy et al. 2013b; den Brok et al. 2021), and a galaxy with inclination i,

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{mol}}\,[{M}_{\odot }\,{\mathrm{pc}}^{-2}]=6.7\,{I}_{\mathrm{CO}\ (2-1)}\,[{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}]\,\cos i.\end{eqnarray} \tag{ 11 }$$

Then, for our target per-channel sensitivity of 0.17 K and if we assume a total line width of Δv = 5 km s⁻¹ for the spectral line, the 1σ noise in Σ_mol is Σ_mol ≈ 5.7 M_⊙ pc⁻² before any inclination correction.

For comparison, typical GMC surface densities range from a few times 10 M_⊙ pc⁻² to a few times 100 M_⊙ pc⁻². Therefore, this sensitivity level allows us to detect every location where a GMC fills a reasonable fraction of the synthesized beam at good significance.

Under the same assumptions for ${\alpha }_{\mathrm{CO}}^{1-0}=4.35$ M_⊙ pc⁻² (K km s⁻¹)⁻¹, R₂₁ = 0.65, and fiducial line width, our target sensitivity corresponds to 1σ point-mass sensitivity of about 4 × 10³ M_⊙, 1.5 × 10⁴ M_⊙, and 4.4 × 10⁴ M_⊙ at distances of 5, 10, and 17 Mpc, respectively. Our chosen target sensitivity is motivated by the canonical value for the characteristic mass of GMCs in the inner Milky Way, M_mol ≳ 10⁵ M_⊙ (e.g., Blitz 1993), and thus, even in our most distant targets, we expect to detect individual moderately massive GMCs (M_mol ≳ 10⁵ M_⊙) at 2−3σ significance. Note that this Δv = 5 km s⁻¹ is the fiducial value used for sensitivity calculations. As described above, we observed at finer spectral resolution than this and our nominal data products are constructed using ∼2.5 km s⁻¹ channels.

The achieved surface brightness sensitivity depends sensitively on the final beam size of the image, and our final resolutions tend to be slightly coarser than 1''. Therefore, in practice our images tend to have somewhat better sensitivity than this 0.17 K. We discuss the achieved properties of the PHANGS–ALMA data more below.

This target sensitivity drove our main array integration times to low values, ≈12–24 s per pointing. Integrating on a single pointing under typical Band 6 conditions, achieving 7.5 mJy beam⁻¹ per 5 km s⁻¹ channel takes ALMA ∼30 s using 43 antennas. After accounting for the effect of overlapping mosaic pointings, the actual nominal integration time per field dropped to ∼12 s. Quantization of ALMA integration times and scheduling block creation tended to increase this per field integration time to ≈18–24 s, which set our final sensitivity (see Section 7.1.2).

Our targets all have a large extent compared to the 30'' primary beam of the 12 m ALMA antennas at ν = 230 GHz. We cover them using mosaics consisting of many individual pointings, typically 100–450 per galaxy. To select the exact target area, we use mid-infrared emission, which we take to indicate the presence of recent star formation and thus the likely presence of molecular gas (see Section 4.6.1 for more discussion on this topic).

Specifically, when designing the field of view for the Large Program, we aimed to cover the full area of each galaxy that shows a WISE3 intensity above 0.5 MJy sr⁻¹ at 75 resolution. The pilot programs targeted the area with WISE4 intensity above 1 MJy sr⁻¹ at 15'' resolution, which yielded qualitatively similar results. Based on the tight scaling between mid-infrared intensity and CO intensity discussed in Section 4.6.1, this design implies that we will cover all parts of the galaxy where ${I}_{\mathrm{CO}}^{2-1}\gtrsim 1$ K km s⁻¹ at low resolution. Roughly, this means we cover all regions where we expect molecular gas to dominate the cold gas budget in the galaxy (e.g., see Leroy et al. 2008) and also all regions where the average intensity over a large area exceeds our 1σ sensitivity at high resolution.

As illustrated in Figure 20, we defined the actual coverage region for each galaxy as a rotated rectangle, with the position angle and extent matched to the WISE contour by eye. In a few cases, a rectangle was not a good match to the exact WISE contour, e.g., because a small amount of mid-infrared emission would force us to dramatically expand the mosaic and thus the integration time. In these cases, we made a judgment call to adjust the exact field of view to a more practical value.

Zoom In Zoom Out Reset image size

In practice, as Section 4.6.1, Table 9, and Figure 20 show, our adopted field of view typically encloses 80% of the WISE3 emission and almost 90% of the WISE4 emission from each target. The aperture corrections in Tables 4 give a more detailed estimate of how well our field of view captures the area of interest for each target. These values should indicate the fraction of total CO emission captured by our field of view. Qualitatively, our maps are more compact than those made by array receivers on single-dish telescopes (e.g., HERACLES; Leroy et al. 2009), but they do cover almost all of the active star formation and molecular gas in each target.

Within each rectangle, the ALMA observing tool places hexagonally packed fields with pointings spaced by $\lambda /\sqrt{3}D$ or ∼13''. During each of our observing cycles, ALMA allowed no more than 150 pointings per "Science Goal." This restriction placed a maximum size on individual mosaics, and also defined the field of view for the 7 m and TP observations. Many of our targets required two or even three of these maximum-sized 150 pointing mosaics to cover the mid-infrared defined regions. For these targets, we observed the galaxy in several "parts." All parts shared the same spectral setup, but we defined each part as a separate Science Goal with a distinct field of view. The fields of view of different Science Goals targeting the same galaxy overlapped one another by ∼10'' along the shared edge to ensure uniform sensitivity. The separate Science Goals were scheduled, observed, and mostly reduced separately. The left panel in Figure 20 illustrates a three-part case, NGC 2997. In this galaxy, ALMA observed each of the three rectangular fields separately. The figure demonstrates the overlap between fields that we use to achieve a good final image of the whole galaxy. The right panel shows NGC 3059, where a single <150 pointing mosaic covers the whole galaxy. Both panels show the rotated rectangular fields that we use and the final coverage of the different arrays, with the total power data slightly extended relative to the ACA 7 m data and both covering the footprint of the 12 m data.

As we discuss below, observing the galaxies in separate parts led to uneven resolution in some cases. Our processing matches the angular resolution of the parts at the coarsest common beam, which led to a modest overall loss in resolution in our data products. In practice, a modest difference in surface brightness sensitivity between the individual parts of a mosaic is likely to be the most noticeable consequence of these multipart observations.

Following standard ALMA practice, the total power observations used an "on-the-fly mapping" approach to cover the area of interest. Individual observations used multiple antennas to cover the same area. The total power observations were descended from the same Science Goals as the interferometric observations. As a result, when a galaxy was observed in multiple parts spread across several Science Goals, the total power observations were similarly divided to cover individual parts of the galaxy.

We periodically reviewed the ALMA archive to find all observations that match our PHANGS–ALMA setup. We searched for observations of nearby, d ≲ 20 Mpc, galaxies targeting CO(2−1) and covering all or most of the galaxy's star-forming area at ∼100 pc resolution. As of our latest review, in fall 2020, we identified four galaxies that closely matched the PHANGS–ALMA setup and included these in our sample. The targets are NGC 1365, NGC 5128 (Centaurus A), NGC 5236 (M83), and NGC 7793; and we list the project codes and PIs in Table 2. We include NGC 1365 "as is" from the archive. For Centaurus A, the CO(2−1) part of the total power data in the archive was not usable due to a problem in the observations. We obtained new ACA mapping in ALMA Cycle 7 (P.I. C. Faesi, see Table 2). NGC 5236 was observed over two cycles using nine separate Science Goals. The data have a large extent, but otherwise match the PHANGS–ALMA setup exactly. The archival NGC 7793 observations obtained only 12 m main array data. We obtained 7 m and total power data in ALMA Cycle 7 (P.I. C. Faesi, see Table 2).

For this version of a PHANGS–ALMA data release, we considered only the CO(2−1) line and only processed wide-area maps of star-forming galaxies. By including the CO(1−0) and CO(3−2) lines and relaxing the areal coverage requirements, many more archival observations of molecular line emission in nearby galaxies could be included. The PHANGS–ALMA data processing pipeline (Leroy et al. 2021b) can straightforwardly handle these data, and we intend to return to this in future work.