Cloud-scale Molecular Gas Properties in 15 Nearby Galaxies

From the high-resolution CO imaging described in Section 2, we estimate the molecular gas surface density, Σ, and velocity dispersion, σ, along each sightline at a range of fixed spatial scales. This approach has been advocated by Leroy et al. (2016) based on earlier works by Ossenkopf & Mac Low (2002), Sawada et al. (2012), Hughes et al. (2013b), and Leroy et al. (2013). Recent work analyzing high-resolution and high-sensitivity ALMA data has also adopted similar approaches (e.g., Egusa et al. 2018).

This approach gives us access to all essential physical properties that we would like to measure (e.g., gas surface density, velocity dispersion, dynamical state, and turbulent energy content). Such a "fixed-spatial-scale" (or simply "fixed-scale") approach is nonparametric, minimal in assumptions, and easy to apply to many data sets in a uniform way. Our measurements are easy to replicate in synthetic observations and thus offer a straightforward path for direct comparison between observations and simulations. Finally, this approach characterizes all detected CO emission and allows us to rigorously treat the selection function.

Our fixed-scale approach differs from the cloud identification approaches commonly used in previous studies (e.g., Clumpfind, see Williams et al. 1994; cprops; see Rosolowsky & Leroy 2006). These methods segment CO emission by associating CO emission with local maxima. While this is a useful strategy for identifying isolated structures, there are three overlapping drawbacks that make us prefer the fixed-scale approach. First, the criteria for segmentation are usually not physically motivated. Instead, several recent studies (Pineda et al. 2009; Hughes et al. 2013a; Leroy et al. 2016) have shown that cloud identification methods tend to find beam-sized objects when applied to data sets with moderate resolution and sensitivity. Second, at 45–120 pc resolution, we observe many crowded regions since we target the molecular-gas-rich parts of star-forming galaxies (e.g., galaxy centers, spiral arms, and bar ends). The commonly used segmentation algorithms are not optimized for identifying marginally resolved clouds in this regime. Third, some segmentation algorithms do not characterize all emissions, which then makes the selection function complex. Based on all of these, we believe that the fixed-scale approach is at least as appropriate in this context as cloud identification.

Our fixed-scale approach does not provide information on the spatial extent of the CO-emitting structures. For the rest of this study, we adopt the spatial resolution of the data to be the relevant size scale when, e.g., estimating the mean volume density and virial parameter. We expect this to be a reasonable assumption as long as two conditions hold: (1) the beam size is within the scale-free range of the hierarchical molecular ISM (often taken to be either the disk scale height or the turbulent driving scale, and typically about a hundred or a few hundreds of parsec; see, e.g., Brunt 2003) and (2) bright CO emission fills a large fraction of the beam. When condition (2) is violated, the beam-averaged surface density is diluted by the "dark area" (i.e., beam dilution; see Section 4), and we no longer expect the beam size to represent the relevant size scale. Nevertheless, this concern is not specific to the fixed-scale approach, but generally applicable to most analysis using data sets with moderate spatial resolution.

The sampling scale is an adjustable variable in the fixed-scale approach. We explore the impact of using different sampling scales by changing the linear resolution of the data (see Section 3.2) and repeat our measurements at each resolution. This lets us investigate molecular gas properties as a function of averaging scale.

Before performing any measurements, we preprocess all data sets by convolving them to a set of round Gaussian-shaped beams with fixed linear sizes: 45, 60, 80, 100, and 120 pc (at FWHM). When the native angular resolution of a data set is coarser than the target beam size, we exclude that galaxy from the analysis at that spatial scale.¹⁷ Then, at each resolution, we resample the data sets onto regular square grids so that the pixels Nyquist-sample the new beam, resulting in an (areal) oversampling factor of π/ln2 ≈ 4.53.

In each data set and at each resolution, we identify all sightlines¹⁸ with significant CO emission. To do this, we first identify all 3D regions in the data cube where ≥2 consecutive channels show emission with signal-to-noise ratio (S/N) ≥ 5. Then, we expand this mask in both spatial and spectral directions to include all adjacent pixels that contain CO emission in ≥2 consecutive channels with S/N ≥ 2. This signal identification scheme is similar to those adopted in other works (e.g., cprops; see Rosolowsky & Leroy 2006).

Note that for most of our targets, we detect CO emission at relatively high S/N and much of the emission is spatially connected. In these cases, our sensitivity limit mostly reduces to the two-consecutive 2σ-channel criterion. However, our two low stellar mass PHANGS-ALMA galaxies, NGC 2835 and NGC 5068, have lower overall S/N. For these targets, the two-consecutive 5σ-channel criterion becomes more relevant. We project the selection effects due to these limits into the key parameter space that we study (Section 5.2).

Faint line wings extending to velocities far from the centroid can have an important effect on the line width. With this in mind, we expand the mask along the velocity axis to cover the entire probable velocity range of the CO line. The amount by which we grow the mask depends on the ratio between the intensity at the line peak and at the edge of the mask, and hence varies between sightlines. Based on the peak-to-edge ratio, we expand the mask so that it would extend to ±3σ if the line profile has a Gaussian shape, with our measured peak intensity centered at the peak velocity.

For each sightline in the mask, we calculate the line-integrated intensity I_CO and peak intensity T_peak of the emission within the expanded mask. We further convert these CO line measurements into the molecular gas surface density Σ and velocity dispersion σ following the procedures described in the next two sections (Sections 3.3 and 3.4).

For the line-integrated intensity, peak intensity, and all other measurements, we estimate their statistical errors by performing a Monte Carlo simulation using the data cube itself as a model. We generate 1000 realizations of the data cube by artificially adding random noise to the original cube and repeat all our measurements for each realization. Note that the mask is only generated once for each galaxy using the original cube, and then it is applied to all 1000 mock cubes. We record the rms scatter of the repeated measurements and quote these numbers as the statistical uncertainties.

We estimate a surface density, Σ, from the line-integrated CO intensity, I_CO, along each detected sightline. Throughout this paper, we assume that the low-J CO lines trace mass in a simple way, so that we can estimate Σ from I_CO and an adopted CO-to-H₂ conversion factor (Bolatto et al. 2013). We adopt the following conversion factors α_CO ≡ Σ/I_CO:

$$\begin{eqnarray}&&{\alpha }_{\mathrm{CO}(1-0)}=4.35\,{M}_{\odot }\,{\mathrm{pc}}^{-2}{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1})}^{-1},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\alpha }_{\mathrm{CO}(2-1)}=6.25\,{M}_{\odot }\,{\mathrm{pc}}^{-2}{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1})}^{-1},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\alpha }_{\mathrm{CO}(3-2)}=17.4\,{M}_{\odot }\,{\mathrm{pc}}^{-2}{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1})}^{-1}.\end{eqnarray} \tag{ 3 }$$

These correspond to the Galactic value for CO(1–0) recommended by Bolatto et al. (2013) along with transition-line ratios of CO(2–1)/CO(1–0) = 0.7 (Leroy et al. 2013; Saintonge et al. 2017) and CO(3–2)/CO(1–0) = 0.25 for the Antennae (Ueda et al. 2012; Bigiel et al. 2015). These take the contribution from helium and other heavy elements into account.

A single conversion factor has the merit of showing our measurements directly, but we have good reason to believe that α_CO varies across our sample (Blanc et al. 2013; Sandstrom et al. 2013). Despite a general understanding of the likely variations in α_CO (Bolatto et al. 2013), we still lack a quantitative, observationally verified prescription for α_CO that applies at cloud scales. Moreover, cloud-scale models of the CO-to-H₂ conversion factor often identify the line width of a cloud and/or its density as an important quantity (Maloney & Black 1988; Narayanan et al. 2012). Thus, our measurements may have a complicated, nonlinear interaction with α_CO. We adopt a single α_CO and then discuss possible variations when they become relevant. A more general exploration of α_CO for the CO(2–1) line is a broad goal of the PHANGS-ALMA science program.

As most of the galaxies in our sample have relatively high metallicity and PHANGS-ALMA targets their molecular-gas-rich inner regions, we expect most of the molecular gas to produce bright CO emission. That is, "CO-dark" molecular gas should not represent a dominant fraction of the molecular gas mass, and the CO line width should be representative of the true H₂ velocity dispersion. We do include some low-mass spiral galaxies in our analysis (M33, NGC 2835, NGC 5068) and here the contributions of a CO-dark phase may be more significant (though α_CO likely varies by less than a factor of 2; see Gratier et al. 2017).

We assume that our spatial resolution is sufficient to reach the scale of individual GMCs or small collections of GMCs, and we take such structures to be spherically symmetric. Therefore, we do not apply any correction to Σ to account for beam-filling factors or the effect of galaxy disk inclination.

We use the width of the CO emission line to trace the velocity dispersion, σ, of molecular gas along each sightline. For the most part, we expect that the CO line width is driven by turbulent broadening, and thus it directly traces the turbulent velocity dispersion at the spatial scale set by the data resolution.

Several methods exist to estimate the width of an emission line. Common approaches include fitting the line profiles with a Gaussian function or Hermite polynomials, calculating the second moment (i.e., the rms dispersion of the spectrum), or measuring an "effective width" (see below). These methods have varying levels of robustness against noise and make different assumptions about the shape of the line.

For our main results, we use the "effective width"¹⁹ as a proxy for the line width. Following Heyer et al. (2001), we define the effective width as

$$\begin{eqnarray}&&{\sigma }_{\mathrm{measured}}=\displaystyle \frac{{I}_{\mathrm{CO}}}{\sqrt{2\pi }\,{T}_{\mathrm{peak}}}\,,\end{eqnarray} \tag{ 4 }$$

where T_peak is the specific intensity at the line peak (in K). This proxy is less sensitive to noise in line wings than the second moment, and unlike direct profile fitting, it does not a priori assume any particular line shape. However, the effective width can be sensitive to the channel width.

To correct for the broadening caused by a finite channel width and spectral response curve width, we subtract the effective width of the spectral response from the measured effective width (Rosolowsky & Leroy 2006):

$$\begin{eqnarray}&&\sigma \,=\,\sqrt{{\sigma }_{\mathrm{measured}}^{2}-{\sigma }_{\mathrm{response}}^{2}},\end{eqnarray} \tag{ 5 }$$

where σ_response is estimated from the channel width and the channel-to-channel correlation coefficient, following Leroy et al. (2016). Appendix B presents the detailed procedures for estimating σ_response, as well as discussions on the applicability of this "deconvolution-in-quadrature" approach. This broadening correction should be accurate enough in most cases, except when the measured line width becomes close to the channel width. In the following sections, we indicate this regime in our plots and discuss resolution effects, if relevant, when presenting our results.

One important caveat related to the line width near the centers of galaxies is that we do not correct for unresolved rotational motions or contributions from multiple clouds along the line of sight. This "beam-smearing" effect may be important where the rotation curve rises quickly, even at our high resolution. The ongoing efforts of measuring rotation curves from PHANGS-ALMA data (P. Lang et al., in preparation) will allow a more careful treatment of this effect in the future.

Table 2 reports the number of independent measurements and the total emission within the mask, expressed as molecular gas mass, in each target at 45, 80, and 120 pc resolution. In massive disk galaxies, thousands of independent sightlines show CO emission above our sensitivity limit. This number drops to the hundreds in our lower-mass targets. For our adopted CO-to-H₂ conversion factor, the molecular mass implied by the CO flux along each sightline is ∼10⁴–10⁶ M_⊙, and the implied surface density is ∼10–10³ M_⊙ pc⁻².

In total, our detected sightlines include CO emission equivalent to 10⁷–10⁹ M_⊙ of molecular gas per galaxy. Some emission remains outside the mask, however, as it is too faint to be detected at good significance using our CO emission identification method. We estimate the fraction of the CO flux included in our analysis, f_M, by dividing the total flux inside the mask by the sum of the unmasked cube. Because most of our data cubes (all but the Antennae galaxies) incorporate total power data, and most targets have strong enough CO emission, we expect a direct sum of the cube to yield a robust estimate of the total flux.

We report f_M for each data set at each resolution in Table 2. At 80 pc resolution, our completeness is 50%–100% in most galaxies, and slightly less than this in the molecule-poor low-mass galaxies NGC 2835 and NGC 5068. The completeness improves with increasing beam size, reflecting improved surface brightness sensitivity at coarser resolution. It also varies from source to source, depending on the typical brightness of CO in the galaxy (correlating with the molecular gas surface density) and the distance to the galaxy (relating to the surface brightness sensitivity at fixed resolution).

In early studies of the Milky Way cloud population (e.g., Solomon et al. 1987), GMCs are usually assumed to be long lived, roughly virialized structures. Departures from virial equilibrium can be expressed via the virial parameter, α_vir ≡ 2K/U_g, where K and U_g denote the kinetic energy and self-gravitational potential energy. Virialized clouds without surface pressure or magnetic support have α_vir = 1, while marginally bound clouds have α_vir ≈ 2, as do molecular clouds in free-fall collapse (e.g., Ballesteros-Paredes et al. 2011; Camacho et al. 2016; Ibáñez-Mejía et al. 2016).

For nearly spherical clouds, α_vir can be expressed as (Bertoldi & McKee 1992)²⁰

$$\begin{eqnarray}&&{\alpha }_{\mathrm{vir}}\equiv \displaystyle \frac{2K}{{U}_{g}}=\displaystyle \frac{5\,{\sigma }^{2}R}{{fGM}}\,.\end{eqnarray} \tag{ 6 }$$

Here M, R, and σ refer to the cloud mass, radius, and the one-dimensional velocity dispersion. f is a geometrical factor that quantifies the density structure inside the cloud. For spherical clouds with a radial density profile of ρ(r) ∝ r^−γ, f = (1 − γ/3)/(1 – 2γ/5) (Bertoldi & McKee 1992).

Equation (6) implies a relationship between the line width σ, size R, virial parameter α_vir, and surface density Σ of a cloud:

$$\begin{eqnarray}\sigma & = & {\left(\displaystyle \frac{f{\alpha }_{\mathrm{vir}}G}{5}\right)}^{0.5}{R}^{-0.5}{M}^{0.5}\\ & = & {\left(\displaystyle \frac{f{\alpha }_{\mathrm{vir}}G\pi }{5}\right)}^{0.5}{R}^{0.5}{{\rm{\Sigma }}}^{0.5}.\end{eqnarray} \tag{ 7 }$$

Here, Σ is the surface density averaged over the projected area of the cloud on the sky. Equation (7) can be restated as σ²/R ∝ α_vir Σ, so that for a fiducial size–line-width relation σ = v₀ R^0.5, the coefficient v₀ depends on the cloud surface density and virial parameter (Solomon et al. 1987; Heyer et al. 2009). Following Heyer et al. (2009) and the recent extragalactic work discussed in Section 1, the surface densities of molecular clouds are observed to vary in different galactic environments, and Equation (7) has become a key diagnostic for the dynamical state of gas in galaxies.

The derivation above represents a highly idealized view. In reality, the molecular ISM has a complex structure, and we do not expect spherically symmetric clouds with simple density profiles. Nevertheless, if cloud substructure—here parameterized by f—does not vary significantly, we can assume a constant value of f and obtain meaningful relative measurements of the dynamical state of the molecular gas, even if the absolute value of α_vir remains uncertain. This comparative approach has often been used in the extragalactic literature (e.g., Rosolowsky & Leroy 2006; Bolatto et al. 2008), and this is the view we adopt in this work.

When contributions to the gravitational potential other than self-gravity become significant, simply comparing K and U_g does not provide a full description of a cloud's dynamical state. However, insight can still be gained by examining the deviation of cloud line widths and comparing the observed line widths to the expectation for an isolated, self-gravitating cloud.

In particular, the role of external pressure (P_ext) on cloud line widths has been emphasized in recent studies of Galactic and extragalactic clouds (e.g., Heyer et al. 2001; Field et al. 2011; A. Schruba et al., in preparation). For a fixed dynamical state (e.g., virialized, marginally bound) and density profile, we expect a cloud subjected to a high surface pressure to show a larger line width σ than the same cloud without surface pressure. At fixed α_vir, R, and P_ext, the detailed shape of the σ–Σ relation depends on the subcloud density profile (e.g., Field et al. 2011; Meidt 2016). We expect σ to be nearly flat as a function of Σ near the surface density, where P_ext ∼ 0.5πGΣ², that is, where the confinements due to self-gravity and external pressure are comparable in strength. At higher surface densities, the external pressure plays only a modest role, and we expect Equation (7) to still hold, perhaps with a slightly shallower slope.

At low cloud surface densities (P_ext > 0.5πGΣ²), the external pressure exceeds the cloud's self-gravitational pressure. If we assume such a cloud to be in pressure equilibrium, then its internal kinetic energy density (or equivalently, internal turbulent pressure P_turb) scales with the external pressure P_ext (Hughes et al. 2013a). Therefore, we expect the σ–Σ relation to asymptote to an isobaric relation with P_turb ≈ ρσ² ≈ P_ext. For gas clouds with a line-of-sight depth ∼2R, the expected scaling is

$$\begin{eqnarray}\sigma & \approx & {P}_{\mathrm{ext}}^{0.5}{\rho }^{-0.5}\\ & \approx & {(2{P}_{\mathrm{ext}})}^{0.5}{R}^{0.5}{{\rm{\Sigma }}}^{-0.5}.\end{eqnarray} \tag{ 8 }$$

In this case, σ and Σ are inversely correlated because for the same kinetic energy density (thus gas pressure), denser gas should have lower velocity dispersion than less dense gas.

Though we motivate Equation (8) by considering clouds confined by external pressure, this isobaric relation should be a general limit. Whenever the ambient pressure in a medium significantly exceeds self-gravity, we can expect our "cloud" to move toward pressure equilibrium with the surrounding gas. In a realistic cold ISM, this should be the case when the ambient pressure in the disk becomes high relative to the self-gravity of the cold molecular clouds. Then, we expect the clouds to follow an isobaric relation (defined by the local ambient pressure as a "pressure floor"; Keto & Myers 1986; Elmegreen 1989; Field et al. 2011; A. Schruba et al., in preparation), specifying the relationship between the cloud surface density and velocity dispersion. Moreover, the ambient pressure is expected to vary as a function of location in the galaxy in response to the distribution of gas and the potential of the galaxy (e.g., Elmegreen 1989; Wolfire et al. 2003; Ostriker et al. 2010; Herrera-Camus et al. 2017; Meidt et al. 2018). Given this, the low surface density, "pressure-dominated" limit should not be a single σ–Σ relation with a −0.5 slope across the entire galaxy, but rather a group of curves, each defined by the local ambient pressure value.

We measure Σ and σ from data cubes convolved to a common spatial resolution corresponding to the size of a typical Galactic GMC (2R = 45–120 pc). We view these measurements as characterizing the molecular ISM at a scale comparable to these observing beam sizes, r_beam, and thus do not engage in any further structure finding. Following this logic, we will mostly discuss our results by equating each individual beam to a molecular cloud. For a fixed R = r_beam, the expected σ–Σ correlation from Equation (7) should then be

$$\begin{eqnarray}&&\sigma \propto {\alpha }_{\mathrm{vir}}^{0.5}\,{{\rm{\Sigma }}}^{0.5},\end{eqnarray} \tag{ 9 }$$

with R = r_beam now part of the coefficient.

There are a few caveats that should be kept in mind when interpreting our measurement: the unknown line-of-sight depth of the emission, the possible coincidence of physically unassociated structures along the line of sight, and the effect of any mismatch between the beam size and the size of physical structures. Several of these concerns are not exclusively associated with the fixed-scale approach, but are relevant to most studies using position–position–velocity data cubes with moderate resolution.

Line-of-sight depth: For our fixed-scale measurements, the sampled size scale on the sky is known by construction—it is the spatial scale to which we convolve the data. However, we do not have an independent constraint on the line-of-sight depth l. In this paper, we take l ∼ 2r_beam; that is, we assume that the spatial scale sampled along the line of sight is comparable to the scale that we study on the sky.

Because the size R in Equation (7) represents the geometric mean of the size in each of three dimensions, R ∼ (${{r}_{\mathrm{beam}}}^{2}l$)^1/3, we have σ ∝ l^1/6Σ^0.5, which depends weakly on l. Variations in l would add noise or a mild systematic to the slope of our measured σ–Σ relation. However, the effect is expected to be small: a factor of 2 variation in l across a decade in Σ would imply a change in slope of ∼0.05.

Capturing unassociated structures. Line-of-sight depth variations alone have a mild impact, but in the case where the beam samples an ensemble of unassociated objects aligned along the line of sight, we expect larger effects. In this case, we would measure a higher surface density, since the emission from multiple structures is added together. If the line width responds to a larger scale potential well, then we would expect the measured line width to be much larger as well. For a sightline through a heavily populated extended gas disk, we would at a minimum capture any "intercloud" velocity dispersion (e.g., Solomon & de Zafra 1975; Stark 1984; Wilson et al. 2011; Caldú-Primo & Schruba 2016), reflecting either the broader turbulent cascade or the motion of clouds in the larger galactic potential.

In regions with strong velocity gradients, we further expect increased velocity dispersion due to "beam smearing" or velocity fields unresolved by the beam (e.g., Colombo et al. 2014b; Meidt et al. 2018). Such a situation could lead to complex line profiles showing multiple components. We expect this situation to arise most often in the dense inner parts of galaxies, and likely in bars and spiral arms as well, where noncircular motions and shocks are strong and the chance of capturing multiple unassociated structures along one sightline is highest. This situation is also more likely to occur for observations of highly inclined galaxies.

Contamination of bright PSF wings. A related concern arises when the beam samples the edge of a large cloud or the extended "wings" of a beam dominated by a nearby bright object. In both cases, we might expect the measured σ to remain larger, still indicative of the gravitational potential of the whole cloud or the σ value found in the bright, nearby source. However, we do not expect to see molecular clouds with sizes vastly larger than our beam size (45–120 pc), so the main sense of this bias in our data will be related to the wings of the PSF. We might expect to find a mild inflation in σ along faint sightlines near isolated bright sightlines. Given the molecular-gas-rich environment for most of our targets and the observed tight correlation between Σ and σ, we do not expect this bias to play a major role.

Individual unresolved clouds. At the other extreme, we can imagine isolated gas structures on scales much smaller than the beam size (R ≪ r_beam). In the case of a single small cloud within the beam, the beam-averaged surface density no longer traces the cloud's surface density. However, the line width and total gas mass, M = Σ A_beam, are still faithful measurements of the cloud's properties. Thus, from the first half of Equation (7), we can derive the expected relation for individual unresolved, virialized clouds:

$$\begin{eqnarray}&&\sigma \propto {R}^{-0.5}{M}^{0.5}\propto {R}^{-0.5}{{\rm{\Sigma }}}^{0.5}.\end{eqnarray} \tag{ 10 }$$

Here, R is the radius of the cloud, not the beam. We expect R to be positively correlated with M and thus Σ = M/A_beam. We therefore expect the slope of the σ–Σ relation to be shallower than 0.5 when R < r_beam. Moreover, when inferring α_vir from Equation (6), if we are still substituting R with r_beam in this case, we will overestimate α_vir by ∼r_beam/R. This situation is most relevant in low gas density regions where molecular clouds are small and sparse.

Synthesis. We do not expect a significant impact on the measured scaling relation due to variations in the line-of-sight depth. It is also unlikely that the contamination from bright sightlines through PSF wings is significant in our sample. We do expect beam smearing to be important in high-density regions and unresolved structures to be prevalent in low-density regions. Both Σ and σ could be overestimated in the former case, and Σ could be underestimated in the latter case.

Figures 1 and 2 show the distribution of molecular gas mass (as inferred from CO flux) for each galaxy as a function of molecular gas surface density, Σ, measured at 80 and 120 pc resolution, respectively. Figure 3 shows the corresponding distributions as a function of velocity dispersion, σ, measured at 120 pc resolution. For each galaxy, we include all sightlines with detected CO emission across the whole field of view (FoV).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For most of the molecular gas properties considered in this work, their median values often show systematic variation across the spatial scales that we consider, while the shape of their distribution functions and the rank order of galaxies barely vary. Therefore, in the subsequent parts of this section, most of the tables report the median values and widths of all the relevant measurements at 45, 80, and 120 pc scales, while most of the figures only illustrate results at 120 pc scale (where we have available data for all targets).

5.1.1. Central and Disk Distributions

5.1.2. Width and Median of the Distributions

5.1.3. Scale Dependence and Clumping

In Section 4, we argued that the scaling relation between surface density, Σ, and line width, σ, reflects the physical state of the gas structure in each beam. If molecular gas in galaxies tends toward a universal dynamical configuration, then we expect to observe a correlation (e.g., fixed α_vir) or anticorrelation (fixed P_turb) between these two quantities. The gas maps (Figure 12) and distribution functions that we present already reveal some similarities between the distributions of Σ and σ, indicating the existence of a positive correlation.

In Figure 4, we show the line width, σ, as a function of surface density, Σ, at 120 pc resolution across our main sample. In each panel, we plot all measurements in one galaxy as colored or white filled circles. For reference to the larger population, we also plot all measurements across this main sample as gray contours (showing data density levels including 60%, 90%, and 99% of points) in the background.

Zoom In Zoom Out Reset image size

To highlight the difference between the central regions and the disks, we separate sightlines into disk and center populations. We calculate the median and rms scatter in σ in 0.2 dex wide bins of fixed Σ for both populations and show them as black filled squares or unfilled triangles with error bars.

Our selection criteria introduce a bias into the analysis (Section 3.2). Since detection depends on the S/N in a channel-by-channel sense, we preferentially pick out sightlines with narrow line profiles at a fixed line-integrated intensity. This can affect our measurement at low Σ, near the sensitivity limit. In Figure 4, we plot the sensitivity limits imposed by our selection criteria as yellow and red shaded regions.²¹ These represent the two-consecutive 5σ-channel and two-consecutive 2σ-channel criteria, respectively (see Section 3.2).

Our completeness will be less than 100% throughout the yellow hatched region, and it rapidly drops to zero inside the red region. As discussed above, in most targets, the two-at-2σ criterion represents the relevant case over most of the area in most of our targets, and completeness is still reasonably high in the yellow region. Thus, there is only a sharp edge at the boundary of the red region. But for NGC 2835 and NGC 5068, the two-at-5σ criterion is more restrictive and completeness through the yellow region is lower.

We also label the σ_measured = Δv threshold of each data set with a horizontal dotted line. Although we account for the broadening of the line due to the finite channel width and spectral response curve via Equation (5), we do not expect σ to be reliable much below this value.

Overall, Figure 4 shows the equivalent of the σ²/r ∝ Σ scaling relation (Heyer et al. 2009, Section 4) for around 30,000 independent beams spanning 12 nearby star-forming galaxies. Although this relationship is usually studied for individual clouds or a subcloud structure, Figure 4 shows that a version of this σ–Σ relation holds sightline by sightline across our main sample. The relation spans three to four orders of magnitude in surface density and two decades in line width, and appears to be a fundamental property of the molecular ISM at cloud scales.

5.2.1. Star-forming Galaxy Disks

In Figure 4, the diagonal lines show the σ ∝ Σ^0.5 expectation for gas structures with α_vir = 1 (dashed) and 2 (dashed–dotted), assuming our fiducial sub-beam geometries, no surface pressure term, and taking the beam size as the relevant size scale. By eye, the disk population of these star-forming galaxies all show scaling relations that are approximately parallel to these lines, which implies a roughly constant α_vir within each galaxy disk region.

To provide a more quantitative description of the observed scaling relations, we model the data for each galaxy disk at each resolution with a power law of the form

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\left(\displaystyle \frac{\sigma }{\mathrm{km}\,{{\rm{s}}}^{-1}}\right)=\beta {\mathrm{log}}_{10}\left(\displaystyle \frac{{\rm{\Sigma }}}{{10}^{2}\,{M}_{\odot }\,{\mathrm{pc}}^{-2}}\right)+A,\end{eqnarray} \tag{ 11 }$$

where β is the power-law index and A is the normalization of the fit at Σ = 10² M_⊙ pc⁻². We also include an intrinsic scatter along the σ direction in our model, for which we denote the rms scatter in the logarithmic space as Δ_intr.

To take both the selection effect and the measurement uncertainties into account, we find the best-fit model parameters and their associated uncertainties using a Markov Chain Monte Carlo (MCMC) method (as implemented in emcee by Foreman-Mackey et al. 2013). We include only the data points above our two-consecutive 5σ selection criterion (i.e., the blue points in Figure 4 or the rows flagged as "complete" in Table 3) as the input for MCMC, and also take into account this truncation in our Bayesian model. We assume that for each galaxy at each resolution, the statistical uncertainties in σ and Σ, as well as the correlation in their uncertainties, can be described by the estimated values from the Monte Carlo simulation that we describe in the last paragraph of Section 3.2. We report all derived model parameters in Table 6. A more detailed description of the MCMC setup, definition of the priors and likelihood functions, as well as the distribution–correlation plots are presented in Appendix C.

Table 6.  Best-fit Parameters of the σ–Σ Scaling Relation in Galaxy Disks

Galaxy	At 45 pc Resolution			At 80 pc Resolution			At 120 pc Resolution
	β	A	Δintr(Δtot)	β	A	Δintr(Δtot)	β	A	Δintr(Δtot)
NGC 628	0.50a	0.63a	0.00(0.10)	0.38a	0.74a	0.03(0.11)	0.34a	0.79a	0.02(0.11)
NGC 1672	⋯	⋯	⋯	⋯	⋯	⋯	0.63a	0.91a	0.08(0.12)
NGC 2835	0.51a	0.60a	0.00(0.09)	0.49a	0.74a	0.00(0.09)	${0.45}_{-0.09}^{+0.09}$	0.83a	0.00(0.08)
NGC 3351	⋯	⋯	⋯	0.56a	0.76a	0.00(0.09)	${0.41}_{-0.06}^{+0.06}$	0.82a	0.00(0.09)
NGC 3627	⋯	⋯	⋯	0.44a	0.84a	0.10(0.14)	0.42a	0.92a	0.10(0.15)
NGC 4254	⋯	⋯	⋯	⋯	⋯	⋯	0.34a	0.81a	0.08(0.11)
NGC 4303	⋯	⋯	⋯	⋯	⋯	⋯	0.41a	0.79a	0.08(0.10)
NGC 4321	⋯	⋯	⋯	⋯	⋯	⋯	0.50a	0.82a	0.07(0.10)
NGC 4535	⋯	⋯	⋯	⋯	⋯	⋯	0.53a	0.83a	0.07(0.10)
NGC 5068	N/Ab	N/Ab	N/Ab	N/Ab	N/Ab	N/Ab	0.50a	0.74a	0.04(0.11)
NGC 6744	⋯	⋯	⋯	0.50a	0.69a	0.05(0.09)	0.47a	0.76a	0.04(0.10)
M51	0.49a	0.69a	0.09(0.11)	0.40a	0.83a	0.10(0.11)	0.36a	0.90a	0.10(0.11)
M31	${0.4}_{-0.1}^{+0.1}$	0.66a	0.00(0.23)	${0.4}_{-0.1}^{+0.1}$	${0.77}_{-0.05}^{+0.05}$	0.00(0.25)	$-{0.0}_{-0.3}^{+0.2}$	${0.7}_{-0.1}^{+0.1}$	0.00(0.18)
M33	⋯	⋯	⋯	${0.5}_{-0.2}^{+0.2}$	0.79a	0.00(0.31)	${0.4}_{-0.4}^{+0.3}$	${0.9}_{-0.1}^{+0.1}$	0.00(0.24)
Antennae	⋯	⋯	⋯	0.50a	0.71a	0.12(0.14)	0.47a	0.81a	0.13(0.15)
PHANGS+M51	0.48a	0.66a	0.07(0.11)	0.47a	0.85a	0.10(0.14)	0.37a	0.85a	0.12(0.13)

Notes. For each galaxy at each resolution level, we report (1) the best-fit power-law slope β (see Equation (11)), (2) the best-fit power-law normalization A, and (3) the scatter in σ at fixed Σ around the best-fit scaling relation (in units of dex; we report both the estimated intrinsic scatter, Δ_intr, and the total observed rms scatter, Δ_tot). The last row shows the best-fit values when combining the main sample together (see Section 5.2.1 and Appendix C).

^aThe MCMC sampling suggests that the statistical uncertainties on these β and A estimates are smaller than 0.05 (i.e., insignificant compared to the systematic uncertainties). We suggest adopting a total uncertainty of 0.05 for these β and A estimates (see Section 5.2.1). ^b"N/A" means the MCMC sampling does not converge.

Download table as:  ASCIITypeset image

At 80 pc resolution, we find best-fit power-law slopes of β = 0.34–0.63 in individual galaxy disks. The best-fit normalization at Σ = 100 M_⊙/pc² is 4–8 km s⁻¹. If we combine the data at 80 pc resolution for all galaxies in the main sample (i.e., the row named "PHANGS+M51" in Table 6), we find a best-fit σ–Σ relation of

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\left(\displaystyle \frac{\sigma }{\mathrm{km}\,{{\rm{s}}}^{-1}}\right)=0.47\,{\mathrm{log}}_{10}\left(\displaystyle \frac{{\rm{\Sigma }}}{{10}^{2}\,{M}_{\odot }\,{\mathrm{pc}}^{-2}}\right)+0.85\end{eqnarray} \tag{ 12 }$$

for this whole sample of 12 galaxies. This relation holds for sightlines with molecular gas surface density over the range Σ ∼ 20–2000 M_⊙ pc⁻².

Due to the huge number of data points involved, the formal statistical errors on the best-fit power-law slope and zero point are often quite small. However, by comparing these MCMC results with those obtained from other fitting strategies (e.g., binning by Σ then fitting, changing weighting schemes), we find that different methods produce results that differ by ∼0.05 for both β and A. Since the best-fit results depend mostly on which fitting scheme we adopt, for comparisons with other studies or applications in theoretical models/numerical simulations as empirical relations, we adopt a typical uncertainty of 0.05 in β and A. For the few cases where we find >0.05 statistical errors from the MCMC method, we instead explicitly report these values in Table 6.

We find the intrinsic scatter in σ at fixed Σ to be less than 0.10 dex for all galaxies at all resolutions. In a few cases (NGC 628 at 45 pc; NGC 2835, NGC 3351, M31, and M33 at all resolutions), the MCMC modeling returns a best-fit value of zero for the intrinsic scatter. Given that the measurement uncertainties in σ and Σ are comparable to the total observed rms scatter around the best-fit relation, these results most likely imply that we are overestimating the measurement uncertainties in σ and/or Σ for these galaxies. For clarity, in Table 6 we also report the total rms scatter around the best-fit relation, which includes both the intrinsic scatter and the measurement uncertainties.

Our power-law fits suggest that molecular gas has similar dynamical properties in star-forming galaxy disks. Considering the uncertainty in β, about half of the best-fit values are consistent with 0.5, as expected for resolved, self-gravitating structures (see Section 4). For the remaining targets, most show shallower slopes and high σ in the low-Σ regime. We discuss the possible origin of such behavior in Section 5.2.5.

5.2.2. Star-forming Galaxy Centers

5.2.3. High Surface Density Regime: The Antennae Galaxies

Our main sample, consisting of 11 PHANGS-ALMA targets and M51, emphasizes relatively massive galaxies (M_⋆ = 10¹⁰–10¹¹ M_⊙) on the star-forming main sequence. These galaxies represent the typical environment for star formation in the local universe. Nevertheless, our ultimate goal is to achieve a quantitative, homogeneously analyzed picture that covers the full range of conditions found in galaxies, from molecule-poor outer disks and dwarf galaxies to gas-rich turbulent disks at high redshift. With the aim of extending our sample toward these extremes, we included the Antennae galaxies, the nearest major merger; M31, a quiescent massive spiral; and M33, a star-forming dwarf spiral. These targets offer clues about the behavior of the scaling relations that we measure outside the active regions of disk galaxies.

In the right panel of Figure 5, we show the σ–Σ relation measured in the interacting region of the Antennae galaxies (orange and white filled circles), on top of the main sample population (gray contours in the background). Molecular gas in this regime mostly populates the upper-right corner in the σ–Σ parameter space. Such extraordinarily high surface density and velocity dispersion are expected and have been noted before, as the galaxy merger brings a huge amount of gas into the interacting region and concentrates it into a small area (Wilson et al. 2003). The complex kinematics of the collision can create large velocity dispersions (Wei et al. 2012).

Zoom In Zoom Out Reset image size

The measured σ–Σ scaling relation for the Antennae lies close to the extrapolation of the average relation in the disks of star-forming galaxies, albeit with a slightly larger scatter (0.12–0.13 dex; see the last row in Table 6). In other words, the measured high line widths match the expectation given the high surface densities and approximately fixed α_vir to first order (see also Leroy et al. 2016). In this case, the main effect of the major merger may be to drive the internal pressure of the gas to high values (see Section 5.4 below) while not substantially altering its observed dynamical state.

We caution that additional caveats, including the CO-to-H₂ conversion factor and line ratios (Wilson et al. 2003; Zhu et al. 2003; Schulz et al. 2007), might have stronger impact on the Antennae galaxies than our other targets. Zhu et al. (2003) suggested that α_CO could be a factor of 2–4 smaller than the Galactic value in this interacting region (but see, e.g., Wilson et al. 2003). If α_CO is indeed smaller than the Galactic value that we adopt, then it means that the true Σ values could be lower than our estimates, and the points could shift leftward in the σ–Σ space. Moreover, in such interacting systems, many of the brightest regions show complex, multicomponent line profiles (Herrera et al. 2012; Johnson et al. 2015), which makes the CO line width measurements trickier. Our effective width approach de-emphasizes the component-to-component width of the line, but some spectral decomposition is likely necessary in future works.

5.2.4. Low Surface Density Regime: M31 and M33

M31 and M33 represent low gas surface density, H i-dominated environments. As illustrated in the left two panels of Figure 5, both of these galaxies have sightlines showing much lower Σ than targets in the main sample (also see Table 4), as well as higher σ compared to an extrapolation of the σ–Σ trend for the main sample. Such behavior suggests that the fixed-scale σ–Σ relation deviates at low Σ from the fixed-α_vir lines.

The distinction between these two targets and other targets in our sample is largely due to the difference in their CO map depth. Because of the proximity of the Local Group targets, CO emission has been surveyed at relatively high sensitivity. Comparing the position of the yellow/red hatched regions in the M31 or M33 panels in Figure 5 to any of the panels in Figure 4, one can see that the M31 and M33 sightlines probe a part of the σ–Σ parameter space that is inaccessible for our other targets at typically 10–20 times the distance.

As a more quantitative way to illustrate this sensitivity effect, we degraded the M33 data to the sensitivity of our more distant targets. We add random noise into the M33 data cube to artificially degrade it, so that its noise level matches that in the data cube for NGC 2835, a galaxy with similar global properties. We then repeat our whole analysis (including regenerating the new mask for signal identification) on these degraded data cubes and measure σ and Σ for all identified sightlines.

We show the result for these mock measurements in Figure 6. After we elevate the noise level, the recovered molecular gas mass faction in M33 drastically drops from 63% to 8%. This low detection rate agrees with the results from previous observations, in which case low-sensitivity interferometer observations only recovered a small fraction (20%) of the CO flux in the galaxy (see Rosolowsky et al. 2007 and references therein).

Zoom In Zoom Out Reset image size

The sightlines that we do detect come from a handful of dense, bright regions. The sensitivity limit imposes a selection bias against low-Σ and high-σ sightlines, and thus, the only remaining detections are those with low σ at the highest Σ. Consequently, at this lower sensitivity, only the sightlines including the most apparently bound gas structures enter the analysis, and the σ–Σ relation becomes close to the α_vir = 2 line.

As the low sensitivity and low detection rate could bias the inferred average α_vir value, we expect that similar issues also affect our measurements of NGC 2835 and NGC 5068, the two PHANGS-ALMA objects with the faintest CO emission and lowest CO flux recovery fractions (30%–40% at 80 pc). Though we do not know for certain where the undetected emission in these targets lies in the σ–Σ space, the most natural expectation would be that it occupies a similar part of the parameter space as it does in M31 and M33. We believe that the systematic deviation of M31 and M33 in the σ–Σ parameter space, as compared to the PHANGS-ALMA targets, implies a flattening of the σ ∝ Σ^0.5 scaling that is common in low molecular gas density environments (Σ ≲ 30 M_⊙ pc⁻²). Limited sensitivity prevents such detections outside the Local Group.

Existing observations in the Milky Way do support this prediction. To show this, we take the molecular cloud catalog published by Miville-Deschênes et al. (2017) as a reference sample that is observed at a much better sensitivity. To replicate the effect of beam dilution in our sample, for all Milky Way clouds with size R smaller than our beam radius r_beam = 80 pc, we calculate their "beam-averaged" surface density as Σ' = Σ R²/r_beam² and plot their σ–Σ' relation in Figure 6 as gray dots. These beam-diluted Milky Way data suggest that we would have also seen a large population of high-σ, low-Σ measurements if we were to observe our own Galaxy at a high sensitivity but at fixed 80 pc spatial resolution. Note that this exercise does not take into account the possibility of catching multiple clouds in a beam, and thus the beam-dilution effect might be slightly stronger than reality.

5.2.5. Interpretations of $\sigma$ Excess at Low Σ

Our fixed-scale Σ and σ measurements provide access to the dynamical state of molecular gas in each beam, commonly expressed via the virial parameter α_vir. Figure 4 and Table 6 show that, to first order, α_vir ∼ σ²/Σ varies modestly across our targets. Meanwhile, the complementary quantity, P_turb ∼ Σ σ² shows an enormous range (see Section 5.4). Both quantities have physical significance, and α_vir−P_turb represents a useful parameter space for diagnosing the state of the molecular gas. Specifically, most modern star formation theories predict a strong dependence of the star formation efficiency per free-fall time on the virial parameter and the turbulent pressure (Kruijssen 2012; Krumholz et al. 2012; Hennebelle & Chabrier 2013; Padoan et al. 2017, among many others). In this section and the next, we examine the distributions of α_vir and P_turb in our sample.

Following the discussion in Section 4, we take the beam size as the relevant size scale (i.e., R = r_beam, thus R = 40 pc at 80 pc resolution), adopt f = 10/9 (appropriate for a density profile of ρ(r) ∝ r⁻¹), and infer α_vir following Equation (6):

$$\begin{eqnarray}{\alpha }_{\mathrm{vir}} & = & \displaystyle \frac{5\,{\sigma }^{2}{r}_{\mathrm{beam}}}{{fG}\,M}=\displaystyle \frac{5}{{fG}}\displaystyle \frac{{\sigma }^{2}{r}_{\mathrm{beam}}}{{\rm{\Sigma }}\,{A}_{\mathrm{beam}}}=\displaystyle \frac{5\mathrm{ln}2}{\pi {fG}}\displaystyle \frac{{\sigma }^{2}}{{\rm{\Sigma }}\,{r}_{\mathrm{beam}}}\\ & \approx & 5.77{\left(\displaystyle \frac{\sigma }{\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{2}{\left(\displaystyle \frac{{\rm{\Sigma }}}{{M}_{\odot }{\mathrm{pc}}^{-2}}\right)}^{-1}{\left(\displaystyle \frac{{r}_{\mathrm{beam}}}{40\mathrm{pc}}\right)}^{-1}.\end{eqnarray} \tag{ 13 }$$

While our adopted prefactors might be subject to systematic errors, this mainly renders the absolute value of α_vir uncertain. We expect the relative sense of our inferred α_vir values to be relatively robust.

5.3.1. Disk and Center Distributions

Figure 7 shows the distribution of molecular gas mass as a function of α_vir, measured at 120 pc resolution. For each galaxy, we plot the distribution for gas in the whole field of view (solid curves) as well as that in only the disk region (dashed curves). Table 7 lists the median value and 16%–84% range of α_vir in each galaxy.

Zoom In Zoom Out Reset image size

Table 7.  Properties of the α_vir Distribution Function

Galaxy	At 45 pc Resolution			At 80 pc Resolution			At 120 pc Resolution
	Disk	Disk	Center	Disk	Disk	Center	Disk	Disk	Center
	Median	16%–84%	Median	Median	16%–84%	Median	Median	16%–84%	Median
	αvir	Width	αvir	αvir	Width	αvir	αvir	Width	αvir
NGC 2835	2.1	0.40	3.5	2.6	0.40	4.3	2.8	0.35	4.3
NGC 5068	1.3	0.47	1.0	1.6	0.52	1.2	1.7	0.48	1.4
NGC 628	2.7	0.41	3.0	2.9	0.48	2.8	2.8	0.51	2.3
NGC 1672	⋯	⋯	⋯	⋯	⋯	⋯	3.2	0.47	2.1
NGC 3351	⋯	⋯	⋯	2.5	0.37	2.2	3.0	0.38	2.5
NGC 3627	⋯	⋯	⋯	3.1	0.62	6.3	2.9	0.63	6.3
NGC 4535	⋯	⋯	⋯	⋯	⋯	⋯	1.9	0.45	2.2
NGC 4254	⋯	⋯	⋯	⋯	⋯	⋯	1.6	0.55	1.2
NGC 4303	⋯	⋯	⋯	⋯	⋯	⋯	1.5	0.47	1.7
NGC 4321	⋯	⋯	⋯	⋯	⋯	⋯	1.9	0.47	1.6
M51	2.4	0.47	1.9	2.1	0.51	1.6	1.9	0.53	1.4
NGC 6744	⋯	⋯	⋯	1.7	0.42	⋯	1.8	0.44	⋯
M31	7.6	0.62	⋯	10.2	0.67	⋯	10.7	0.68	⋯
M33	⋯	⋯	⋯	8.3	0.64	7.1	9.3	0.62	7.8
Antennae	⋯	⋯	⋯	1.5	0.59	⋯	1.3	0.61	⋯

Note. For each galaxy at each resolution, we report (1) the median α_vir value by gas mass for the "disk" population, (2) the full width of the 16%–84% gas mass range of α_vir distribution for the "disk" population (in units of dex), and (3) the median α_vir value by gas mass for the "center" population.

Download table as:  ASCIITypeset image

Most galaxies (excluding M31 and M33) have distributions centered near α_vir ≈ 1.5–3.0, with a 16%–84% width of 0.4–0.65 dex. Because a fixed α_vir corresponds to a slope β = 0.5 in Figure 4, the width of the distributions in Figure 7 captures a mixture of the scatter about the best-fit relation and the deviation from β = 0.5. As we might expect from the results above, the distribution of α_vir appears more uniform among the high-mass star-forming galaxies in the main sample, when focusing only on their disks. In this case, the distribution resembles a log-normal one but usually with a mild skew toward higher α_vir values.

Compared to the large range in Σ, α_vir shows a narrow range of values across our sample. Specifically, while the median value of Σ varies substantially from galaxy to galaxy, the median α_vir varies much less, suggesting that most molecular gas share a common dynamical state across massive, star-forming disk galaxies. As shown above and in the literature, the variations in α_vir that do exist can be linked to the environment. For example, studies by the PAWS survey have shown a link between α_vir and dynamical environment in M51 (Meidt et al. 2013; Colombo et al. 2014a). But in the disks of most of our target galaxies, the strong correlation between Σ and σ with small scatter implies that the variation in α_vir is much smaller than the variations in the surface density, line width, and their combination—the internal turbulent pressure (see Section 5.4).

We find more variation in α_vir when including the galaxy centers. Figure 7 shows a wide distribution of α_vir associated with the central regions of NGC 1672 and NGC 3351. Another two barred galaxies, NGC 2835 and NGC 3627, show relatively high α_vir in their central regions compared to disks (by 25% to 100%). This captures the displacement of the σ–Σ relation to higher line width in galaxy centers. As emphasized above, the line width in central regions may include contributions from bulk gas motions like rotation. Therefore, the interpretation of these high α_vir should be made with some caution. Nevertheless, these line widths are consistent with the idea that the presence of a bar leads to efficient turbulence driving within the inner Lindblad resonance (e.g., Krumholz & Kruijssen 2015).

The peak positions of the α_vir distributions lie between 1.5–3.0, given f = 10/9 and R = r_beam. As discussed in Section 4, marginally bound or free-falling gas should have K ≈ U_g (i.e., "energy equipartition"), which implies α_vir ≈ 2. Thus, our inferred α_vir values indicate K ∼ U_g but with a modest excess of turbulent kinetic energy compared to self-gravitational potential energy. However, considering the uncertainties in the prefactors on α_vir (Equation (13)), which reflects our limited ability to infer an absolute molecular mass and calculate the true gravitational potential, our derived absolute values of α_vir are correspondingly uncertain. We suggest reading Figure 7 as showing a remarkable degree of uniformity in the dynamical state of molecular gas across the disks of galaxies and among galaxies, given matched assumptions and treatment. Given the associated uncertainties, our measurement does not distinguish between the various physical states dominated by self-gravity (e.g., marginally bound, spherical free-fall, hierarchical collapse). In Section 6, we discuss the path toward accurate α_vir estimates.

5.3.2. Outliers

5.3.3. Scale Dependence

Comparing α_vir at different resolutions can reveal the clumpiness and kinematic structure of molecular gas across physical scales. For example, in a homogeneous turbulent medium with a fiducial velocity dispersion–size relation σ ∝ r^0.5, we expect α_vir to show little scale dependence (as the σ and r_beam dependences cancel each other out in Equation (13)). In the other extreme, i.e., in the case of a single isolated cloud unresolved at any accessible resolution, beam dilution will cause Σ to scale with the beam size as ${{r}_{\mathrm{beam}}}^{-2}$, and thus α_vir ∝ Σ⁻¹ r_beam⁻¹ ∝ r_beam.

In Figure 8, we compare the α_vir distribution in the disks of our targets at different resolutions. We consider each galaxy at 45, 60, 80, 100, and 120 pc resolution (when accessible) and take the beam size as the relevant length scale. As in Figure 7, we show the median by mass and 16%–84% mass range of α_vir. We highlight α_vir = 1 and α_vir = 2 with dashed and dashed–dotted vertical lines.

Zoom In Zoom Out Reset image size

Though both σ and Σ depend on spatial scale, the inferred gas dynamical state varies much more weakly. A few of our targets (NGC 2835, NGC 5068, NGC 3351, NGC 6744, M31, and M33) show higher median α_vir at coarser resolution, while the others show no trend or sometimes an inverse trend. These trends are statistically significant due to the large number of independent measurements. However, the changes in the median α_vir never exceed 0.2 dex across available resolutions, which is mild compared to the 0.4 dex change in beam size (from 45 to 120 pc). Thus, the observed scale dependence in α_vir is much weaker than the expected α_vir ∝ r_beam for unresolved isolated clouds, and in at least half of our targets, α_vir shows no or inverse correlation with the sampling scale.

As mentioned above, this analysis is subject to observational selection effects because the fraction of molecular gas mass that enters our calculation decreases as we move from coarser to finer resolution (Table 2). Because the total gas mass being analyzed is not conserved across resolutions, an artificial scale dependence can be introduced by systematically excluding more gas with higher α_vir at finer resolution. In Section 5.2, we measure the slope of the σ–Σ relation to be β < 0.5, which means that gas in low-Σ regions tends to have higher α_vir. Consequently, losing these low-Σ sightlines in finer resolution maps can bias the median α_vir to lower values and thus artificially produce a positive correlation between α_vir and r_beam.

Taking these observational biases into account, our best estimate is that our observations do not provide strong evidence for a dependence of α_vir on spatial scale. Our measurements appear to capture gas with kinetic energy similar to but slightly greater than the self-gravitational potential energy across the scales that we consider (with the exception of M31 and M33).

Besides the virial parameter, our measurements also allow us to infer the internal turbulent pressure, P_turb, or equivalently, the kinetic energy density, in molecular gas. Comparing P_turb to self-gravity, the external pressure in the ambient medium and disk structure offers more insights into the dynamical state of the molecular ISM.

The internal pressure in molecular gas with line-of-sight depth ∼2R can be expressed as

$$\begin{eqnarray}&&{P}_{\mathrm{turb}}\approx \rho {\sigma }^{2}\approx \displaystyle \frac{1}{2R}{\rm{\Sigma }}{\sigma }^{2}.\end{eqnarray} \tag{ 14 }$$

Under the assumption of R = r_beam, we estimate P_turb from²²

$$\begin{eqnarray}&&{P}_{\mathrm{turb}}/{k}_{B}\approx 61.3\,{\rm{K}}\,{\mathrm{cm}}^{-3}\left(\displaystyle \frac{{\rm{\Sigma }}}{{M}_{\odot }\,{\mathrm{pc}}^{-2}}\right){\left(\displaystyle \frac{\sigma }{\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{2}{\left(\displaystyle \frac{{r}_{\mathrm{beam}}}{40\mathrm{pc}}\right)}^{-1}.\end{eqnarray} \tag{ 15 }$$

Figure 9 shows the distribution of P_turb measured at 120 pc scale in all 15 galaxies. In each panel, the distribution for the entire field of view is shown by the solid curve, while that for the disk and center regions are shown by the dashed and dotted curves. The measured median values and widths of the P_turb distribution are reported in Table 8.

Zoom In Zoom Out Reset image size

Based on the functional form of Equation (14) and the observed σ–Σ correlation, one would expect the properties of the P_turb distributions to resemble those of Σ and σ. Indeed, in Figure 9, we see a qualitatively similar behavior to that in Figures 1–3, including the multimodal shape of the distribution functions and the contrast between disk and central populations.

Lines of fixed turbulent pressure P_turb, a.k.a. isobars, appear as diagonal lines (σ ∝ Σ^−0.5) running from top left to bottom right in Figure 4. These lines lie almost perpendicular to the actual distribution of data. In other words, P_turb represents the reprojection of the observed σ–Σ distribution along its principal axis. Sightlines located up and to the right in the plot (e.g., galaxy centers or the densest regions in disks) correspondingly have larger Σσ², and thus higher turbulent pressure.

As expected, we measure the 16%–84% range for P_turb to be wide even in individual galaxy disks (typical widths are ∼1–2 dex). The distributions for the centers of strongly barred galaxies sometimes show two to three orders of magnitude excess in median P_turb, which makes them appear as a distinct peak in the P_turb distribution plots (NGC 1672, NGC 3351, NGC 3627, NGC 4535, NGC 4303, NGC 4321). Such a large dynamical range of turbulent pressure in the molecular ISM has been claimed for both the Milky Way and other nearby galaxies (see Hughes et al. 2013a; Leroy et al. 2015; A. Schruba et al., in preparation, for compilations).

The dynamical range in P_turb across our sample is further boosted by strong intergalaxy variations. We find P_turb ∼ 10³–10⁵ K cm⁻³ in the disks of low-mass galaxies, while P_turb ∼ 10⁵–10⁷ K cm⁻³ for high-mass ones. Such a correlation between P_turb and galaxy stellar mass can be either a consequence of the Σ–M_⋆ correlation we mentioned in Section 5.1 or a more fundamental relation that links the local gas properties to the ambient galactic environment. This topic will be investigated in a forthcoming paper. Here we emphasize that P_turb varies dramatically within individual galaxies, and strongly and systematically among galaxies. The internal pressure of a patch of molecular gas is clearly a strong function of environment.

This work reveals that the molecular gas surface density Σ shows a wide distribution across the local star-forming galaxy population at fixed spatial scales. Recent molecular cloud surveys also find variations in the cloud average surface density across galactic and extragalactic systems (Heyer et al. 2009; Donovan Meyer et al. 2013; Colombo et al. 2014a; Leroy et al. 2015; Egusa et al. 2018, among many others).

The σ–Σ scaling relation we find under the fixed-scale framework (see Figure 10) is the analog of the σ²/R ∝ Σ relation frequently quoted in recent cloud studies. This relation is often discussed in the context of molecular clouds in virial equilibrium. However, subsequent works have suggested alternative dynamical states, including marginally bound, free-falling (Ballesteros-Paredes et al. 2011; Camacho et al. 2016; Ibáñez-Mejía et al. 2016), or pressure-confined clouds (Heyer et al. 2001; Oka et al. 2001; Field et al. 2011; Hughes et al. 2013a; A. Schruba et al., in preparation). Each of these scenarios predicts a particular value of α_vir or some dependence of α_vir on environment.

We found typical values of 1.5–3.0 for the virial parameter α_vir ∝ σ²/Σ, which implies that the molecular ISM often possesses slightly more kinetic energy than gravitational potential energy at the scales we study. However, uncertainties in the prefactors and sub-beam distribution make the normalization of α_vir a poor discriminant among these ideas. One of the main goals of this paper is to lay the groundwork for a direct correlation of molecular gas structure with disk structure, kinematics, and host galaxy properties. We expect such analysis to help distinguish among these scenarios.

We do not see evidence in our sample for the strikingly low line widths recently observed in the inner regions of two early-type galaxies by Davis et al. (2017, 2018). As they discuss, their measurements would fall far below our measured σ–Σ scaling relation, which implies α_vir < 1. This region appears almost totally unpopulated by data in our sample. Their measurements carefully accounted for the rotation of the gas disk and had higher spatial resolution (29 and 13 pc, respectively) than we achieve here. The easiest explanation to reconcile the two sets of observations seems to be invoking differences in sub-beam structure or line-of-sight depth in one or both data sets, but several other effects (e.g., missing short-spacing information, choice of CO excitation line) may also play a role. This remains a topic for future research.

We observe a wide range of turbulent pressure P_turb ∝ Σσ² across our sample. This wide spread in P_turb agrees with previous observations, which reveal vastly different internal pressures in different parts of the Milky Way (e.g., Bertoldi & McKee 1992; Heyer et al. 2001; Oka et al. 2001; Field et al. 2011) and in selected extragalactic systems across the cosmic distance scale (e.g., Swinbank et al. 2012; Kruijssen & Longmore 2013; Livermore et al. 2015). There appears to be no doubt that different galaxies drive their molecular gas to vastly different pressures at cloud scales.

Though not a unique explanation, a correlation of P_turb with environment might be expected if molecular clouds represent overdensities in a self-regulated disk. In self-regulated disk models, the gas maintains radial (e.g., Silk 1997) or vertical (e.g., Elmegreen 1989) force balance in the galaxy potential. The ISM regulates itself around this equilibrium state, with feedback from star formation coupling the large-scale gas velocity dispersion and density. In recent years, this idea has been explored in a series of works using increasingly sophisticated simulations and analytic theories (Koyama & Ostriker 2009; Ostriker et al. 2010; Kim et al. 2011, 2013; Ostriker & Shetty 2011; Kim & Ostriker 2015).

Focusing on equilibrium in the vertical direction (i.e., perpendicular to the disk), these theories predict that over large spatial scales and long timescales, internal pressure in the ISM will balance its weight inside the disk's gravitational potential. In this case, the large-scale dynamic equilibrium pressure, P_DE, may represent some local baseline value that is also relevant to the internal pressure of molecular clouds. This equilibrium pressure correlates positively with the stellar and gas surface densities and thus is higher in massive galaxies and the inner parts of galaxies (Wong & Blitz 2002; Blitz & Rosolowsky 2004, 2006; Leroy et al. 2008).

If P_DE represents a long-term and large-scale average, we expect that the molecular gas that we observe represents an over-pressurization relative to this environment-dependent baseline value. This is qualitatively consistent with our observed trend of higher P_turb in high stellar mass galaxies and in galaxy centers, but it also can be directly tested by combining our data with measurements of local disk structure at various spatial scales. We will present a quantitative investigation of this scenario in a forthcoming paper.

In the previous sections, we present our results in terms of physical quantities, such as the molecular gas surface density, Σ, and velocity dispersion, σ. These quantities are inferred from two direct observables, the surface brightness and effective width of a low-J CO line, mostly CO(2–1). The simple translations between CO line measurements and physical quantities also means that our results can be easily converted back to their observable form whenever necessary. This, together with our fixed-spatial-scale approach, allows for a straightforward comparison between our observations and predictions of the structure of CO emission from simulations of galaxies or analytic theory. Given the recent progress in modeling CO emission from simulations (e.g., Smith et al. 2014; Pan et al. 2015; Duarte-Cabral & Dobbs 2016; Gong et al. 2018), we expect this to be a fruitful application of our measurements in the next few years. In this sense, our results represent the current best available measurements of the distribution of CO emission from star-forming disk galaxies at 45–120 pc resolution and should offer a valuable benchmark to identify the physics needed to produce a realistic molecular medium.

In terms of observables, the correlation shown in Figure 4 is fundamentally a correlation between the CO line-integrated intensity I_CO and the effective width. Since the effective width is directly proportional to I_CO by construction (Equation (4); ignoring the insignificant broadening correction for this discussion), one may wonder how much of the observed σ–Σ correlation results from this by-construction correlation. To address this question, we show in Figure 11 the correlation between peak temperature T_peak and line width σ, which are statistically independent quantities for resolved line profiles. Here we only include the 12 targets with CO(2–1) data to make a fair comparison.

Zoom In Zoom Out Reset image size

Given that Σ ∝ I_CO ∝ σT_peak, an intrinsic σ ∝ Σ^0.5 relation should manifest itself as a slope = 1 correlation between the effective width σ and line peak intensity T_peak (a.k.a. "peak temperature"). Conversely, if the σ–Σ relation is purely artificial, we expect no correlation between σ and T_peak. The observed correlation across the PHANGS sample is indeed consistent with the expectation of σ ∝ T_peak, although the relationship has a larger scatter and, as shown in Appendix D, suffers from more severe selection effects than the σ–Σ correlation. This result confirms that our observed line-width–intensity relation is not artificial in nature, but rather conveys physical information about cloud-scale molecular gas kinematics. We show the galaxy-by-galaxy σ–T_peak relation in Appendix D.

As a side note, the low T_peak found in M33 (widely below 0.1 K) further illustrates the existence of severe beam-dilution effect in this galaxy. As we have mentioned in Section 4 and Section 5.2.4, along with the influence of missing gas and external pressure (A. Schruba et al., in preparation), this beam dilution likely contributes to producing the shallow σ–Σ relation at low Σ.

Much of our interpretation relies on the adopted CO-to-H₂ conversion factor, α_CO, for the appropriate transition. We do expect α_CO to vary between the central regions and disks of our targets (Israel 2009; Sandstrom et al. 2013) and CO-dark gas may play an important role at low Σ (see, e.g., Smith et al. 2014). We adopt a fixed α_CO for this study in order to link our measurements closely to observables. We expect improved measurements of CO excitation and comparison to high-resolution dust maps to help refine the interpretation of our measurements in the coming years. As mentioned above, comparisons to the predicted CO emission from recent simulations (e.g., Smith et al. 2014; Pan et al. 2015; Duarte-Cabral & Dobbs 2016; Gong et al. 2018) should be increasingly feasible in the future, and our approach makes this comparison relatively straightforward.

This analysis uses the initial PHANGS-ALMA sample. This sample is currently being expanded to cover almost all southern, nearby, low-inclination, star-forming galaxies. With the final sample of 74 galaxies, it will be possible to make definitive measurements of distribution functions and scaling relations across the whole local star-forming galaxy population.

In addition, many aspects of our analysis would benefit from better knowledge of the sub-beam gas distribution and the line-of-sight depth through the gas layer. This would improve our translation from surface to volume density (and so pressure), allow for better absolute estimates of the virial parameter, and help distinguish low beam filling from variations in physical properties. We consider this to be a strong argument for even more highly resolved case studies of the nearby galaxy population.

Last but not least, we expect to continue to improve our methodology. The ongoing parallel effort following a cloud identification scheme (E. Rosolowsky et al., in preparation) will provide a complementary description of the molecular gas structures, whereas more sophisticated characterization of CO intensity (e.g., forward-modeling at the cube level) will help improve both the statistical rigor and completeness of our measurements.