Dynamical Equilibrium in the Molecular ISM in 28 Nearby Star-forming Galaxies

We use PHANGS-ALMA CO (2−1) data to trace molecular gas distribution and kinematics across the star-forming disks in all our sample galaxies (A. K. Leroy et al. 2020, in preparation; also see Sun et al. 2018). These CO observations target the actively star-forming area in each of these galaxies, typically covering out to 5–15 kpc in galactic radius. They include data from both the 12 m array and the Morita Atacama Compact Array (ACA; consisting of the 7 m array and four 12 m total power antennas), and therefore should recover emission on all spatial scales. The 1σ noise level in the data is 0.2–0.3 K in a $2.5\,\mathrm{km}\,{{\rm{s}}}^{-1}$ channel. The angular resolution (i.e., beam FWHM) ranges between 10–18, corresponding to 25–120 pc linear resolution at the targets' distances. This allows us to resolve the molecular gas distribution and kinematics at physical scales comparable to the typical size of GMCs (i.e., "cloud scales").

To homogenize the CO data, we convolve all the data cubes to two different linear resolutions, 60 and 120 pc, whenever possible. This allows us to control for resolution-related systematics by comparing results derived at two different linear scales. All CO data in our sample can be matched to 120 pc resolution, whereas only a subset of them (six galaxies) reach 60 pc.

We create CO line intensity (${I}_{\mathrm{CO}};$ or moment-0) maps and line effective width²⁶ (${\sigma }_{\mathrm{CO}}$) maps from the matched resolution CO data cubes. We build these maps by analyzing only significant CO detections in the cube (selected by the "strict" signal masks described in A. K. Leroy et al. 2020, in preparation). This strategy ensures high signal to noise in the derived CO line intensity and line width maps, but discards a fraction of CO flux existing at low signal to noise (see, e.g., Table 2 in Sun et al. 2018). Below, in Section 3.1, we quantify this effect by measuring a CO flux recovery fraction.

Table 2.  List of Key Physical Properties

Quantity Definition	Symbol	Unit	Data Source
Ensemble average of cloud-scale molecular gas properties (see Section 3.3.1)
(measured at $\theta =60$, 120 pc scale; averaged over each kpc-sized aperture; see Section 3)
Average molecular gas surface density (Equations (1) and (2))	${\left\langle {{\rm{\Sigma }}}_{\mathrm{mol},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$	${M}_{\odot }\,{\mathrm{pc}}^{-2}$	PHANGS-ALMA CO (2−1)
Average molecular gas velocity dispersion (Equations (1) and (3))	${\left\langle {\sigma }_{\mathrm{mol},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$	km s−1	PHANGS-ALMA CO (2−1)
CO flux recovery fraction (Section 3.1)	${f}_{\mathrm{CO},\theta \mathrm{pc}}$	⋯	PHANGS-ALMA CO (2−1)
Environmental characteristics (see Section 3.3.2)
(measured at 1 kpc scale)
Kiloparsec-scale molecular gas surface density (Equation (5))	${{\rm{\Sigma }}}_{\mathrm{mol},1\mathrm{kpc}}$	${M}_{\odot }\,{\mathrm{pc}}^{-2}$	PHANGS-ALMA CO (2−1)
Kiloparsec-scale atomic gas surface density (Equation (6))	${{\rm{\Sigma }}}_{\mathrm{atom},1\mathrm{kpc}}$	${M}_{\odot }\,{\mathrm{pc}}^{-2}$	PHANGS-VLA H i 21 cm, etc.
Kiloparsec-scale stellar mass surface density (Equation (7))	${{\rm{\Sigma }}}_{\star ,1\mathrm{kpc}}$	${M}_{\odot }\,{\mathrm{pc}}^{-2}$	${{\rm{S}}}^{4}{\rm{G}}$ IRAC 3.6 μm
Kiloparsec-scale star formation rate surface density (Equation (8))	${{\rm{\Sigma }}}_{\mathrm{SFR},1\mathrm{kpc}}$	M⊙ yr−1 kpc−2	z0MGS NUV+MIR
Pressure estimates (see Section 4)
Average turbulent pressure in molecular gas (Equations (10) and (11))	${\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$	${k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$	PHANGS-ALMA CO (2−1)
Kiloparsec-scale ISM equilibrium pressure (Equations (12)–(14))	${P}_{\mathrm{DE},1\mathrm{kpc}}$	${k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$	All combined
Average cloud-scale equilibrium pressure (Equations (15)–(18))	${\left\langle {P}_{\mathrm{DE},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$	${k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$	All combined
Average weight of molecular clouds (Equations (16) and (17))	${\left\langle {{ \mathcal W }}_{\mathrm{cloud},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$	${k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$	All combined
Average weight of clouds due to self-gravity (Equations (16) and (17))	$\langle {{ \mathcal W }}_{\mathrm{cloud},\,\theta \mathrm{kpc}}^{\mathrm{self}}{\rangle }_{1\mathrm{kpc}}$	${k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$	PHANGS-ALMA CO (2−1)

Download table as:  ASCIITypeset image

For more details regarding the CO data reduction, readers are referred to A. K. Leroy et al. (2020, in preparation), as well as Sun et al. (2018), which adopts a similar data reduction method.

To measure the distribution of gas, stars, and star formation activity in each of our target galaxies, we assemble a multiwavelength supporting data set from a variety of sources. These supporting data typically have much coarser angular resolution than the PHANGS-ALMA CO data, corresponding to linear scales of hundreds to a thousand parsecs.

We use H i 21 cm emission data to trace the atomic gas distribution. We include new data acquired as part of the PHANGS-VLA project (covering 11 targets in our sample; PI: D. Utomo) and the EveryTHINGS project (NGC 4496A; PI: K. Sandstrom). In addition, we use existing data from VIVA (covering seven targets; Chung et al. 2009), THINGS (six targets; Walter et al. 2008), VLA observations associated with HERACLES (NGC 4321 and NGC 4536; Leroy et al. 2013), and an ATCA observation in the literature (NGC 1792; Murugeshan et al. 2019). Most of these H i data have native angular resolution of 15''–25'', which corresponds to linear scales of 0.5–2 kpc at the distances of the targets. Typical column density sensitivities are $(0.5\mbox{--}2)\times {10}^{20}\ {\mathrm{cm}}^{-2}$ at 3σ after integrating over the 10–20 km s⁻¹ line width.

We use Spitzer IRAC 3.6 μm data to trace the stellar mass distribution. These data come from the Spitzer Survey of Stellar Structure in Galaxies (${{\rm{S}}}^{4}{\rm{G}}$; Sheth et al. 2010). For most targets in our sample, we use the ICA 3.6 μm maps (Querejeta et al. 2015), for which the emission from dust has been subtracted. For the remaining three targets (NGC 4571, NGC 4689, and NGC 5134), we use the raw 3.6 μm maps, as these galaxies have global [3.6]–[4.5] colors compatible with an old stellar population (see Section 4.2 in Querejeta et al. 2015). All these maps are masked to remove foreground stars. These data have a native angular resolution of $\sim 2^{\prime\prime}$, which corresponds to linear scales of hundreds of parsecs for our targets.

We combine GALEX near-UV (NUV) data and WISE mid-IR (MIR) data to derive kpc-scale SFR estimates. We use the post-processed GALEX NUV images and WISE 12 μm images from the z = 0 Multiwavelength Galaxy Synthesis (z0MGS; Leroy et al. 2019). These images are corrected for background, aligned for astrometry, and masked to remove foreground stars. Almost all targets in our sample have both GALEX and WISE coverage (except for NGC 4689, for which only WISE is available). We use the data at fixed angular resolution of 75, which corresponds to linear scales $\lesssim 1\,\mathrm{kpc}$.

We homogenize the resolution of these supporting data to a common 1 kpc linear scale. As described above, most data have native resolutions better than 1 kpc, and thus can be directly convolved to this coarser resolution. For some of the H i data, however, the native resolution is coarser. As the resolution of these H i data is not crucial for our analysis (see discussion in Section 4.3), we keep these data at their native resolution in the following analysis.

To complement these supporting data, we also convolve the PHANGS-ALMA CO data to 1 kpc resolution. This provides a tracer of the large-scale molecular gas distribution. The convolution is performed on the CO data cubes, and from these convolved cubes we derive a set of low-resolution CO line intensity maps. This strategy takes advantage of the better surface brightness sensitivity at coarser resolution. In the resulting kpc-scale CO intensity maps, we detect emission at high signal to noise along almost every sightline, and thus expect to recover essentially all emission within the observation footprint.

Altogether, these data provide us with spatially resolved information about the (molecular and atomic) gas distribution, the stellar disk structure, and the local SFR, all of which are measured on a matched 1 kpc spatial scale.

As discussed above, for our CO intensity-weighted averaging scheme to work optimally, the CO observations should be sensitive enough to detect a significant fraction of CO flux at cloud-scale resolution in each aperture. This is indeed the case for the majority of the apertures in our sample, as shown below. When this is not the case, however, the averaged values will suffer from larger uncertainties and reflect only properties of the brightest clouds.

To control for these effects due to the finite sensitivity of CO observations, we quantify a CO flux recovery fraction, ${f}_{\mathrm{CO},\theta \mathrm{pc}}$, within each 1 kpc aperture. We do this by comparing, within the footprint of each aperture, the total flux included in the 60–120 pc resolution CO line intensity map (i.e., within the "strict" signal masks, see Section 2.1) to the total flux in the corresponding CO data cube. We estimate the latter quantity by summing up the cube within a wide, high-completeness signal mask (referred to as the "broad" mask in A. K. Leroy et al. 2020, in preparation).

Our calculation shows that, for the majority (∼70%) of the kpc-scale apertures in our sample, the CO flux recovery fraction on 120 pc scales, ${f}_{\mathrm{CO},120\mathrm{pc}}$, is higher than 50%. Since the CO intensity-weighted averages are calculated from only the detected emission, a recovery fraction of ${f}_{\mathrm{CO}}\gt 50 \%$ means that the majority of the CO emission in that kpc-scale aperture is included in the averaging. First, this assures that the derived averages in these apertures have reasonably small statistical errors ($\lt \sqrt{2}$ times larger than the case of infinite sensitivity given homoscedastic individual measurements). Second, and most importantly, in the cases where the undetected molecular gas has systematically different properties than the detected, the intensity-weighted averages in these high ${f}_{\mathrm{CO}}$ apertures are also much less susceptible to systematic effects due to sampling biases.

Hereafter, when presenting results derived from this CO intensity-weighted averaging approach, we represent the data in darker/lighter colors to denote higher/lower ${f}_{\mathrm{CO}}$. Measurements shown in darker colors are thus more representative of the overall cloud population, whereas those shown in lighter colors only characterize the brightest clouds in the kpc-sized aperture.

Our weighted averaging scheme is a way to quantify the mean molecular gas properties in each kpc-sized aperture. This kpc aperture size is large enough that the CO flux within each aperture might come from various morphological regions of a galaxy (e.g., bulges, bars). In this work, for each of the averaging apertures, we also keep track of the fractional CO flux contribution from different morphological regions.

In detail, we first identify areas covered by morphological structures like bulges and bars (when applicable) in every galaxy. To identify the bulge regions, we use the ${{\rm{S}}}^{4}{\rm{G}}$ structural decomposition results presented by Salo et al. (2015), there referred to as ${{\rm{S}}}^{4}{\rm{G}}$ pipeline 4 results. These results are based on two-dimensional structural decomposition of Spitzer IRAC 3.6 μm images with GALFIT3.0 (Peng et al. 2010). To identify bar regions, we instead use results from the visual identification in Herrera-Endoqui et al. (2015), which have higher quality than the ${{\rm{S}}}^{4}{\rm{G}}$ pipeline 4 results.

Within each kpc-sized aperture, we calculate the fractional contribution in the total CO flux from the bulge and bar regions. Whenever an aperture includes nonzero CO flux coming from bulge or bar regions, we classify it as a "bulge/bar" aperture.²⁷ Otherwise, we classify it as a "disk" aperture. When presenting our results, we show measurements in bulge/bar apertures in orange colors, and those in disk apertures in blue colors.

We use the PHANGS-ALMA CO data and multiwavelength supporting data to estimate physical properties of the molecular gas and its ambient galactic environment. Here, we detail our methods to convert direct observables (e.g., CO intensity) into physical quantities (e.g., molecular gas surface density) in each kpc-sized aperture in our sample. Table 2 lists the key physical quantities we derive, including both average molecular cloud properties and environment characteristics. These physical quantities are the basis of all the following calculations detailed in Section 4.

3.3.1. Cloud-scale Molecular Gas Properties

Following Sun et al. (2018), we estimate cloud-scale molecular gas surface density (${{\rm{\Sigma }}}_{\mathrm{mol},\theta \mathrm{pc}}$) and velocity dispersion (${\sigma }_{\mathrm{mol},\theta \mathrm{pc}}$) from the observed CO (2−1) line intensity (${I}_{\mathrm{CO},\theta \mathrm{pc}}$) and effective width (${\sigma }_{\mathrm{CO},\theta \mathrm{pc}}$), via

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{mol},\theta \mathrm{pc}}=\displaystyle \frac{{\alpha }_{\mathrm{CO}}}{{R}_{21}}\,{I}_{\mathrm{CO},\theta \mathrm{pc}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\sigma }_{\mathrm{mol},\theta \mathrm{pc}}={\sigma }_{\mathrm{CO},\theta \mathrm{pc}}.\end{eqnarray} \tag{ 3 }$$

In Equation (2), ${R}_{21}=0.7$ is the adopted CO (2−1) to CO (1−0) line ratio (Leroy et al. 2013; Saintonge et al. 2017), whereas ${\alpha }_{\mathrm{CO}}$ is a CO-to-H₂ conversion factor,²⁸ whose value varies aperture-by-aperture.

We adopt the following prescription to predict the values of ${\alpha }_{\mathrm{CO}}$ for each aperture in our sample. Similar to the calibration suggested by Accurso et al. (2017) and adopted in the xCOLD GASS survey (Saintonge et al. 2017), we predict ${\alpha }_{\mathrm{CO}}$ via

$$\begin{eqnarray}&&{\alpha }_{\mathrm{CO}}=4.35\,{Z}^{{\prime} -1.6}\ {M}_{\odot }\,{\mathrm{pc}}^{-2}{\left({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}\right)}^{-1}.\end{eqnarray} \tag{ 4 }$$

Here, Z' is the local gas phase abundance normalized to the solar value appropriate for the Pettini & Pagel (2004) metallicity calibration $[12+\mathrm{log}({\rm{O}}/{\rm{H}})=8.69$].

Using Equation (4) to predict ${\alpha }_{\mathrm{CO}}$ requires knowing Z' for every kpc-sized aperture across our sample. However, metallicity measurements existing in the literature only cover a subset of our targets, and they are derived using heterogeneous calibration methods. To ensure a homogeneous coverage of the entire sample, we instead predict Z' in a uniform, empirical way. Using a mass–metallicity relation reported by Sánchez et al. (2019, see Table 1 therein), we first predict Z' at one effective radius (${R}_{{\rm{e}}}$) in each galaxy based on the galaxy global stellar mass (see Table 1). We then extend our prediction to cover the entire galaxy assuming a universal radial metallicity gradient of $-0.1\ \mathrm{dex}/{R}_{{\rm{e}}}$ (Sánchez et al. 2014). Combining this locally predicted Z' with Equation (4), we have a predicted ${\alpha }_{\mathrm{CO}}$ value for every aperture in our sample.

Our choices on the prescriptions for predicting Z' and ${\alpha }_{\mathrm{CO}}$ could affect many of the measured molecular gas properties in this work. To quantify the systematic effects associated with these choices, in Section 6.1 we consider three alternative ${\alpha }_{\mathrm{CO}}$ prescriptions, and compare the quantitative results with those derived based on our "fiducial" prescriptions.

3.3.2. Kpc-scale Environment Characteristics

Internal pressure in molecular gas includes the contributions from thermal and turbulent motion, as well as magnetic fields. Observational evidence, including superthermal CO line widths and the size–line width relation observed within GMCs, suggest that turbulent motion dominates over thermal motion on physical scales comparable to cloud sizes (Larson 1981; Solomon et al. 1987; Heyer & Brunt 2004; also see Heyer & Dame 2015). Numerical simulations of the star-forming ISM on galactic scales also find that the magnetic term is subdominant, typically reaching only ≲50% of the kinetic term in the effective gas pressure (Kim & Ostriker 2017; Pakmor et al. 2017; Su et al. 2017; Hopkins et al. 2020; also see observational evidence presented by Crutcher 1999; Falgarone et al. 2008; Troland & Crutcher 2008; Thompson et al. 2019). Motivated by these findings, we assume in this work that turbulent motion represents the primary source of internal pressure in molecular gas, and treat all other contributions as subdominant. We do not differentiate between turbulent pressure and total internal pressure in molecular gas hereafter.

Turbulent pressure in molecular gas is commonly estimated from volume density and the observed (one-dimensional) velocity dispersion, under the assumption that turbulence is isotropic:

$$\begin{eqnarray}&&{P}_{\mathrm{turb}}={\rho }_{\mathrm{mol}}{\sigma }_{\mathrm{turb},1{\rm{D}}}^{2}.\end{eqnarray} \tag{ 9 }$$

While one can use the observed velocity dispersion ${\sigma }_{\mathrm{mol}}$ along the line of sight as a proxy of ${\sigma }_{\mathrm{turb},1{\rm{D}}}$, ${\rho }_{\mathrm{mol}}$ is not usually directly observed. Here, we convert cloud-scale molecular gas surface density, ${{\rm{\Sigma }}}_{\mathrm{mol}}$, into volume density, ${\rho }_{\mathrm{mol}}$. To do this, we assume a constant-density spherical cloud filling each beam, with the cloud diameter ${D}_{\mathrm{cloud}}$ equal to the beam FWHM³¹ (i.e., 60 or 120 pc, but see Section 6.3 for an alternative approach). The inferred turbulent pressure in molecular gas can then be expressed as

$$\begin{eqnarray}&&{P}_{\mathrm{turb}}={\rho }_{\mathrm{mol}}{\sigma }_{\mathrm{mol}}^{2}=\displaystyle \frac{6{M}_{\mathrm{mol}}}{\pi {D}_{\mathrm{cloud}}^{3}}{\sigma }_{\mathrm{mol}}^{2}=\displaystyle \frac{3{{\rm{\Sigma }}}_{\mathrm{mol}}{\sigma }_{\mathrm{mol}}^{2}}{2{D}_{\mathrm{cloud}}}.\end{eqnarray} \tag{ 10 }$$

The quantity ${M}_{\mathrm{mol}}={{\rm{\Sigma }}}_{\mathrm{mol}}\cdot (\pi {D}_{\mathrm{cloud}}^{2}/4)$ is the total molecular gas mass of the spherical cloud.

We average the estimated cloud-scale turbulent pressure across each kpc-sized region following the same CO flux weighting scheme described in Section 3:

$$\begin{eqnarray}&&{\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}=\displaystyle \frac{{\int }_{A}{P}_{\mathrm{turb},\theta \mathrm{pc}}\,{I}_{\mathrm{CO},\theta \mathrm{pc}}\,{dS}}{{\int }_{A}{I}_{\mathrm{CO},\theta \mathrm{pc}}\,{dS}}.\end{eqnarray} \tag{ 11 }$$

This can be interpreted as the mass-weighted average turbulent pressure within each kpc-sized aperture A.

To compare with previous works, we first estimate the classic, kpc-scale dynamical equilibrium pressure, ${P}_{\mathrm{DE},1\mathrm{kpc}}$. We follow the same basic formalism that has been adopted, with some variations, in many previous works (e.g., Spitzer 1942; Elmegreen 1989; Elmegreen & Parravano 1994; Wong & Blitz 2002; Blitz & Rosolowsky 2004, 2006; Leroy et al. 2008; Koyama & Ostriker 2009; Ostriker et al. 2010; Kim et al. 2011, 2013; Ostriker & Shetty 2011; Shetty & Ostriker 2012; Hughes et al. 2013a; Kim & Ostriker 2015a; Benincasa et al. 2016; Herrera-Camus et al. 2017; Gallagher et al. 2018b; Fisher et al. 2019; Schruba et al. 2019).

This approach models the distribution of gas and stars in a galaxy disk as isothermal fluids in a plane-parallel geometry. For the calculation here, we assume that the (single component) gas disk scale height is much smaller than the stellar disk scale height. We also neglect gravity due to dark matter, as it represents only a minor component in the inner disks of relatively massive galaxies. In this case, we can express ${P}_{\mathrm{DE}}$ as:

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\mathrm{DE},1\mathrm{kpc}} & = & \ \displaystyle \frac{\pi G}{2}\,{{\rm{\Sigma }}}_{\mathrm{gas},1\,\mathrm{kpc}}^{2}\\ & + & {{\rm{\Sigma }}}_{\mathrm{gas},1\mathrm{kpc}}\sqrt{2G{\rho }_{\star ,1\mathrm{kpc}}}\,{\sigma }_{\mathrm{gas},{\rm{z}}}.\end{array}\end{eqnarray} \tag{ 12 }$$

The first term is the weight of the ISM due to the self-gravity of the ISM disk (see, e.g., Spitzer 1942; Elmegreen 1989). The second term is the weight of the ISM due to stellar gravity (see, e.g., Spitzer 1942; Blitz & Rosolowsky 2004). Here, ${{\rm{\Sigma }}}_{\mathrm{gas},1\mathrm{kpc}}={{\rm{\Sigma }}}_{\mathrm{mol},1\mathrm{kpc}}+{{\rm{\Sigma }}}_{\mathrm{atom},1\mathrm{kpc}}$ is the total gas surface density, ${\rho }_{\star ,1\mathrm{kpc}}$ is stellar mass volume density near disk midplane, and ${\sigma }_{\mathrm{gas},{\rm{z}}}$ is the vertical gas velocity dispersion (a combination of turbulent, thermal, and magnetic terms).

Following Blitz & Rosolowsky (2006), Leroy et al. (2008), and Ostriker et al. (2010), we estimate midplane stellar volume densities from the observed surface densities in each kpc-sized aperture:

$$\begin{eqnarray}&&{\rho }_{\star ,1\mathrm{kpc}}=\displaystyle \frac{{{\rm{\Sigma }}}_{\star ,1\mathrm{kpc}}}{4{H}_{\star }}=\displaystyle \frac{{{\rm{\Sigma }}}_{\star ,1\mathrm{kpc}}}{0.54{R}_{\star }}.\end{eqnarray} \tag{ 13 }$$

The first step assumes an isothermal density profile along the vertical direction (i.e., ${\rho }_{* }(z)\propto {{\rm{sech}} }^{2}[z/(2{H}_{* })]$) with ${H}_{\star }$ being the stellar disk scale height (van der Kruit 1988). The second step assumes a fixed stellar disk flattening ratio ${R}_{\star }/{H}_{\star }=7.3$ (Kregel et al. 2002, also see Appendix B). ${R}_{\star }$ is the radial scale length of the stellar disk, for which we adopt the value from the ${{\rm{S}}}^{4}{\rm{G}}$ photometric decompositions of 3.6 μm images (Salo et al. 2015; see column (6) in Table 1).

For ${\sigma }_{\mathrm{gas},{\rm{z}}}$, we calculate the mass-weighted average velocity dispersion of molecular and atomic phases

$$\begin{eqnarray}&&{\sigma }_{\mathrm{gas},{\rm{z}}}={f}_{\mathrm{mol}}{\left\langle {\sigma }_{\mathrm{mol},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}+(1-{f}_{\mathrm{mol}}){\sigma }_{\mathrm{atom}},\end{eqnarray} \tag{ 14 }$$

where ${f}_{\mathrm{mol}}={{\rm{\Sigma }}}_{\mathrm{mol},1\mathrm{kpc}}\,/({{\rm{\Sigma }}}_{\mathrm{mol},1\mathrm{kpc}}+{{\rm{\Sigma }}}_{\mathrm{atom},1\mathrm{kpc}})$ is the fraction of gas mass in the molecular phase. We adopt a fixed atomic gas velocity dispersion ${\sigma }_{\mathrm{atom}}=10\ \mathrm{km}\,{{\rm{s}}}^{-1}$ (see Leroy et al. 2008; Tamburro et al. 2009; Wilson et al. 2011; Caldú-Primo et al. 2013; Mogotsi et al. 2016).

Our adopted assumptions for the ${\rho }_{\star ,1\mathrm{kpc}}$ and ${\sigma }_{\mathrm{gas},{\rm{z}}}$ estimation might introduce systematic biases in the derived ${P}_{\mathrm{DE},1\mathrm{kpc}}$. In Section 6.2, we estimate ${P}_{\mathrm{DE},1\mathrm{kpc}}$ by adopting two alternatives for estimating ${\sigma }_{\mathrm{gas},{\rm{z}}}$ and ${\rho }_{\star ,1\mathrm{kpc}}$, and compare the results with our fiducial ${P}_{\mathrm{DE},1\mathrm{kpc}}$ estimates.

The classic, kpc-scale equilibrium pressure defined in Section 4.2 does not account for gas substructure within each kpc-sized aperture. For atomic gas, surface density fluctuations on the sub-kpc scale are usually moderate (Leroy et al. 2013; also see Bolatto et al. 2011; E. Koch et al. 2020, in preparation), so this issue is likely minor. For molecular gas, however, we expect strong clumping (Leroy et al. 2013). Therefore, gas self-gravity should be significantly enhanced in overdense regions (e.g., in molecular clouds), and the required pressure in molecular gas to balance this enhanced gravity should exceed the classic, kpc-scale pressure estimates.

To account for this, we introduce a modified, cloud-scale dynamical equilibrium pressure, ${\left\langle {P}_{\mathrm{DE},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$. Using the classic formulation as a starting point, we treat the clumpy molecular ISM and diffuse atomic ISM separately, allowing them to have different geometry (also see Ostriker et al. 2010 and Schruba et al. 2019 for similar calculations). We offer a brief summary of this alternative formalism here, but leave a more detailed description of the derivation and adopted assumptions to Appendix A.

In this alternative formalism, we split the total cloud-scale equilibrium pressure into two parts:

$$\begin{eqnarray}\begin{array}{rcl}{\left\langle {P}_{\mathrm{DE},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}} & = & {\left\langle {{ \mathcal W }}_{\mathrm{total},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}\\ & = & {\left\langle {{ \mathcal W }}_{\mathrm{cloud},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}+{{ \mathcal W }}_{\mathrm{atom},1\mathrm{kpc}}.\end{array}\end{eqnarray} \tag{ 15 }$$

The first part, ${\left\langle {{ \mathcal W }}_{\mathrm{cloud},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$, corresponds to the weight of molecular gas, most of which resides in clumpy structures on cloud scales. For simplicity, we assume that all molecular gas is organized into spherical, cloud-like structures. The weight within each individual structure is due to (1) its own self-gravity, (2) the gravity of other molecular structures, and (3) the gravity of stars:

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal W }}_{\mathrm{cloud},\theta \mathrm{pc}} & = & \,{{ \mathcal W }}_{\mathrm{cloud},\theta \mathrm{pc}}^{\mathrm{self}}+{{ \mathcal W }}_{\mathrm{cloud},\theta \mathrm{pc}}^{\mathrm{ext} \mbox{-} \mathrm{mol}}+{{ \mathcal W }}_{\mathrm{cloud},\theta \mathrm{pc}}^{\mathrm{star}}\\ & = & \,\displaystyle \frac{3\pi }{8}G{{\rm{\Sigma }}}_{\mathrm{mol},\theta \mathrm{pc}}^{2}+\displaystyle \frac{\pi }{2}G{{\rm{\Sigma }}}_{\mathrm{mol},\theta \mathrm{pc}}{{\rm{\Sigma }}}_{\mathrm{mol},1\mathrm{kpc}}\\ & & +\,\displaystyle \frac{3\pi }{4}G{\rho }_{\star ,1\mathrm{kpc}}{{\rm{\Sigma }}}_{\mathrm{mol},\theta \mathrm{pc}}{D}_{\mathrm{cloud}}.\end{array}\end{eqnarray} \tag{ 16 }$$

Consistent with our estimation of turbulent pressure in Equation (10), we also assume the cloud diameter ${D}_{\mathrm{cloud}}$ equals the beam FWHM here.

We then adopt the same averaging scheme used in Equation (11) to estimate the (CO flux-weighted) average ${{ \mathcal W }}_{\mathrm{cloud},\theta \mathrm{pc}}$ across each kpc-sized aperture:

$$\begin{eqnarray}&&{\left\langle {{ \mathcal W }}_{\mathrm{cloud},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}=\displaystyle \frac{{\int }_{A}{{ \mathcal W }}_{\mathrm{cloud},\theta \mathrm{pc}}\,{I}_{\mathrm{CO},\theta \mathrm{pc}}\,{dS}}{{\int }_{A}{I}_{\mathrm{CO},\theta \mathrm{pc}}\,{dS}}.\end{eqnarray} \tag{ 17 }$$

The second term in Equation (15), ${{ \mathcal W }}_{\mathrm{atom},1\mathrm{kpc}}$, corresponds to the weight of the smooth extended layer of atomic gas. This weight is due to the gravity of all gas (both atomic and molecular phases) plus the stars, as felt by the atomic layer. Motivated by the relative smoothness of the atomic gas distribution and the coarser resolution of the H i data, we estimate this weight using only kpc-scale measurements, assuming uniform atomic gas surface density within each kpc-sized aperture:

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal W }}_{\mathrm{atom},1\mathrm{kpc}} & = & \,{{ \mathcal W }}_{\mathrm{atom},1\,\mathrm{kpc}}^{\mathrm{self}}+{{ \mathcal W }}_{\mathrm{atom},1\,\mathrm{kpc}}^{\mathrm{mol}}+{{ \mathcal W }}_{\mathrm{atom},1\,\mathrm{kpc}}^{\mathrm{star}}\\ & = & \,\displaystyle \frac{\pi G}{2}{{\rm{\Sigma }}}_{\mathrm{atom},1\,\mathrm{kpc}}^{2}+\pi G{{\rm{\Sigma }}}_{\mathrm{atom},1\mathrm{kpc}}{{\rm{\Sigma }}}_{\mathrm{mol},1\mathrm{kpc}}\\ & & +\ {{\rm{\Sigma }}}_{\mathrm{atom},1\mathrm{kpc}}\sqrt{2G{\rho }_{\star ,1\mathrm{kpc}}}\,{\sigma }_{\mathrm{atom}}.\end{array}\end{eqnarray} \tag{ 18 }$$

Here, we assume that the molecular gas disk is "sandwiched" by the atomic gas disk, and thus the second term above has a prefactor two times larger than that of the first term (see Appendix A for detailed derivation). We adopt ${\sigma }_{\mathrm{atom}}\,=10\ \mathrm{km}\,{{\rm{s}}}^{-1}$, consistent with Section 4.2.

If vertical dynamical equilibrium holds across multiple spatial scales, we would expect molecular gas internal pressure on cloud scales to match the dynamical equilibrium pressure on the same scale. By comparing our measured ${\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$ and best-estimate ${\left\langle {P}_{\mathrm{DE},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$, we can test whether, in a statistical average sense, the molecular ISM in nearby, star-forming disk galaxies can be described by this model.

The top panels in Figure 1 show the average molecular gas turbulent pressure $\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle$ ³² (defined in Section 4.1), measured on $\theta =60$ and 120 pc scales, as a function of the kpc-scale dynamical equilibrium pressure, ${P}_{\mathrm{DE},1\mathrm{kpc}}$ (defined in Section 4.2). Each data point corresponds to one kpc-sized aperture. Darker colors denote higher CO flux recovery fraction (see Section 3), so the $\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle$ measurements are more representative of the bulk molecular gas population within the aperture.

Zoom In Zoom Out Reset image size

At a physical scale of 120 pc, we see $\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle$ values spanning the range 10⁴–${10}^{7}\,{k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$. The corresponding range in ${P}_{\mathrm{DE},1\mathrm{kpc}}$ is 10³–${10}^{6}\,{k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$. For reference, the typical GMC internal pressure (i.e., ${P}_{\mathrm{turb}}$) in the Solar Neighborhood is $\sim {10}^{5}\,{k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$ (see, e.g., Blitz 1993), whereas the estimated local dynamical equilibrium pressure is $\sim {10}^{4}\,{k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$ (see, e.g., Elmegreen 1989). Typical GMC internal pressure in the Galactic Center or nearby galaxy centers is ∼10⁵–${10}^{8}\,{k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$ (Oka et al. 2001; Donovan Meyer et al. 2013; Colombo et al. 2014; Leroy et al. 2015; Sun et al. 2018; Walker et al. 2018; Schruba et al. 2019). Therefore, one may think of our data as spanning from "outer disk" conditions to galaxy centers.

At both 120 and 60 pc resolution, most data points lie above the equality line (solid black line). This suggests that the average internal pressure in molecular gas is usually higher than what is needed to support the weight of a smooth gas disk with its surface density equal to the observed kpc-scale average value.

5.1.1. Quantifying the Overpressurization of the Molecular Gas

5.1.2. Predicting Turbulent Pressure from Kpc-scale Equilibrium Pressure

Since ${P}_{\mathrm{DE},1\mathrm{kpc}}$ contains no infromation about the small-scale gas distribution, it tends to underestimate the true equilibrium pressure on cloud scales. Nevertheless, calculating ${P}_{\mathrm{DE},1\mathrm{kpc}}$ only requires knowing the kpc-scale gas and stellar mass distribution, plus assumptions on the vertical gas velocity dispersion. This makes it possible to estimate ${P}_{\mathrm{DE},1\mathrm{kpc}}$ from low-resolution observations of more distant galaxies, from low-resolution numerical simulations, or even from analytic and semianalytic models of galaxies. If the ISM in other environments follows the same ${\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$–${P}_{\mathrm{DE},1\mathrm{kpc}}$ relation that we observe in Figure 1, then one could use an estimated ${P}_{\mathrm{DE},1\mathrm{kpc}}$ to predict the turbulent pressure on cloud scales.

To make this prediction possible, we fit an empirical ${\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$–${P}_{\mathrm{DE},1\mathrm{kpc}}$ scaling relation. This can be seen as a benchmark relation for the ISM in local star-forming disk galaxies. We derive this relation by fitting a power law to all the disk measurements (blue circles in Figure 1), using the ordinary least square (OLS) method in logarithmic space, and treating ${P}_{\mathrm{DE},1\mathrm{kpc}}$ as the independent variable. This yields best-fit power-law relations (blue dashed lines in Figure 1, top panels):

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}}{{10}^{5}\ {k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}} & = & 3.2\ {\left(\displaystyle \frac{{P}_{\mathrm{DE},1\mathrm{kpc}}}{{10}^{5}{k}_{{\rm{B}}}{\rm{K}}{\mathrm{cm}}^{-3}}\right)}^{1.07},\\ \displaystyle \frac{{\left\langle {P}_{\mathrm{turb},60\mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}}{{10}^{5}\ {k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}} & = & 9.0\ {\left(\displaystyle \frac{{P}_{\mathrm{DE},1\mathrm{kpc}}}{{10}^{5}{k}_{{\rm{B}}}{\rm{K}}{\mathrm{cm}}^{-3}}\right)}^{1.32}.\end{array}\end{eqnarray} \tag{ 19 }$$

We report the scatter around these best-fit relations, as well as the estimated statistical uncertainties on the fitting parameters, in Table 3.

Table 3.  Summary of the Best-fit Power-law Relations in Section 5

Relation	Figures/Equations	Slope	Offset along y-axis	Scatter along y-axis
			at ${10}^{5}\,{k}_{{\rm{B}}}{\rm{K}}\,{\mathrm{cm}}^{-3}$	around Relation
$\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle$–${P}_{\mathrm{DE},1\mathrm{kpc}}$	Figure 1; Equation (19)	$1.07[\pm 0.03]$ a	$0.50[\pm 0.02]$ a dex	0.36 dex
$\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle$–${P}_{\mathrm{DE},1\mathrm{kpc}}$ (debiasb)	Figure 1; Equation (20)	$1.37[\pm 0.05]$ a	$0.60[\pm 0.02]$ a dex	0.32 dex
$\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle$–$\left\langle {P}_{\mathrm{DE},120\mathrm{pc}}\right\rangle$	Figure 2; Equation (21)	$1.02[\pm 0.01]$ a	$-0.12[\pm 0.01]$ a dex	0.17 dex
$\left\langle {P}_{\mathrm{turb},60\mathrm{pc}}\right\rangle$–${P}_{\mathrm{DE},1\mathrm{kpc}}$	Figure 1; Equation (19)	$1.32[\pm 0.06]$ a	$0.96[\pm 0.04]$ a dex	0.31 dex
$\left\langle {P}_{\mathrm{turb},60\mathrm{pc}}\right\rangle$–${P}_{\mathrm{DE},1\mathrm{kpc}}$ (debiasb)	Figure 1; Equation (20)	$1.47[\pm 0.09]$ a	$1.00[\pm 0.04]$ a dex	0.26 dex
$\left\langle {P}_{\mathrm{turb},60\mathrm{pc}}\right\rangle$–$\left\langle {P}_{\mathrm{DE},60\mathrm{pc}}\right\rangle$	Figure 2; Equation (21)	$0.95[\pm 0.02]$ a	$0.06[\pm 0.01]$ a dex	0.13 dex

Notes.

^aAll the quoted errors here are statistical errors estimated from bootstrapping. However, we expect systematic errors to dominate the total uncertainties on these parameters. See Section 6. ^bThese relations are derived in the range ${P}_{\mathrm{DE}}\gt 2\times {10}^{4}\,{k}_{{\rm{B}}}{\rm{K}}\,{\mathrm{cm}}^{-3}$. Caution should be used when extrapolating outside this range.

Download table as:  ASCIITypeset image

We caution that Equation (19) likely has a shallower slope than the actual ${\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$–${P}_{\mathrm{DE},1\mathrm{kpc}}$ relation. This is largely due to the asymmetric data censoring on ${\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$ and ${P}_{\mathrm{DE},1\mathrm{kpc}}$. The PHANGS-ALMA CO observations have higher surface brightness sensitivity at coarser spatial resolution (Section 2.2; also see discussion in Sun et al. 2018). This means that our cloud-scale measurements cannot probe molecular gas surface density as low as our kpc-scale measurements can, and our pressure estimates suffer from a similar censoring effect. The impact of this is even visible in the top left panel in Figure 1: there are few data points with ${\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}\lesssim {10}^{4}\,{k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$, and those few measurements all suffer from low CO recovery fraction. Though not as easily discernible in the top right panel in Figure 1, a similar censoring effect is also present at 60 pc resolution.

To quantify how this data censoring biases our empirical ${\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$–${P}_{\mathrm{DE},1\mathrm{kpc}}$ fit, we calculate another version of the best-fit relation using the same OLS method, but only fitting data points with ${P}_{\mathrm{DE},1\mathrm{kpc}}\gt 2\times {10}^{4}\,{k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$. The impact of the data censoring is much less prominent above this threshold, and thus we expect these fitting results to be less affected. With this fitting strategy, we have (dashed–dotted lines in Figure 1, top panels)

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}}{{10}^{5}\ {k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}} & = & \,4.0{\left(\displaystyle \frac{{P}_{\mathrm{DE},1\mathrm{kpc}}}{{10}^{5}{k}_{{\rm{B}}}{\rm{K}}{\mathrm{cm}}^{-3}}\right)}^{1.37},\\ \displaystyle \frac{{\left\langle {P}_{\mathrm{turb},60\mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}}{{10}^{5}\ {k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}} & = & \,10{\left(\displaystyle \frac{{P}_{\mathrm{DE},1\mathrm{kpc}}}{{10}^{5}{k}_{{\rm{B}}}{\rm{K}}{\mathrm{cm}}^{-3}}\right)}^{1.47}.\end{array}\end{eqnarray} \tag{ 20 }$$

Again, we report statistical uncertainties and residual scatters in Table 3.

For the purpose of predicting ${\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$ from ${P}_{\mathrm{DE},1\mathrm{kpc}}$ in low-resolution observations or simulations, we recommend using Equation (20) in the regime where ${P}_{\mathrm{DE},1\mathrm{kpc}}\gt 2\,\times {10}^{4}\,{k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$. More sensitive observations are needed to pin down this ${\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$–${P}_{\mathrm{DE},1\mathrm{kpc}}$ relation in lower-pressure regimes.

In Section 5.1, we find that ${\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$ exceeds ${P}_{\mathrm{DE},1\mathrm{kpc}}$ in almost all regions across our sample. Our hypothesis is that this reflects molecular gas clumping on small scales, which is left unaccounted for in the kpc-scale ${P}_{\mathrm{DE}}$ estimate. We directly test this hypothesis in Figure 2, where we show the average molecular gas turbulent pressure ${\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$ as a function of the cloud-scale dynamical equilibrium pressure ${\left\langle {P}_{\mathrm{DE},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$ (as derived in Section 4.3). This is a direct "apples to apples" comparison, in the sense that both ${\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$ and ${\left\langle {P}_{\mathrm{DE},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$ are derived on the same $\theta =60$ pc and 120 pc physical scales.

Zoom In Zoom Out Reset image size

We observe a tight, almost linear relation between $\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle$ and $\left\langle {P}_{\mathrm{DE},120\mathrm{pc}}\right\rangle$ across more than three orders of magnitude. We find a similarly strong correlation on 60 pc scale, though with a slightly different normalization. Compared to Figure 1, the observed distribution in Figure 2 shows a relationship much closer to equality (as expected), and much less scatter around the relation as well.

Figure 2 shows that: (1) the dynamical equilibrium model is able to predict the observed turbulent pressure in molecular gas based on the resolved gas and stellar mass distribution in galaxy disks; and (2) to correctly estimate equilibrium pressure within the more clumpy molecular component, it is crucial to account for small-scale density structures, which are only accessible in high spatial resolution observations.

5.2.1. Differentiating Morphological Regions

5.2.2. Quantifying the Turbulent Pressure–Equilibrium Pressure Relation on Cloud Scales

To get a quantitative description of the observed $\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle$–$\left\langle {P}_{\mathrm{DE},\theta \mathrm{pc}}\right\rangle$ relation, we fit a power law to all disk measurements (blue dots). We weight each measurement by its corresponding CO recovery fraction ${f}_{\mathrm{CO}}$, and perform an OLS bisector fit (Isobe et al. 1990) in logarithmic space. Combining all disk measurements, we find the following best-fit relations (shown as blue dashed lines in Figure 2) on 120 and 60 pc scales:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}}{{10}^{5}\ {k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}} & = & \,0.77{\left(\displaystyle \frac{{\left\langle {P}_{\mathrm{DE},120\mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}}{{10}^{5}{k}_{{\rm{B}}}{\rm{K}}{\mathrm{cm}}^{-3}}\right)}^{1.02},\\ \displaystyle \frac{{\left\langle {P}_{\mathrm{turb},60\mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}}{{10}^{5}\ {k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}} & = & \,1.1{\left(\displaystyle \frac{{\left\langle {P}_{\mathrm{turb},60\mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}}{{10}^{5}{k}_{{\rm{B}}}{\rm{K}}{\mathrm{cm}}^{-3}}\right)}^{0.95}.\end{array}\end{eqnarray} \tag{ 21 }$$

The corresponding statistical uncertainties and residual scatters are reported in Table 3.

Thanks to both our large sample size and the tightness of the $\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle$–$\left\langle {P}_{\mathrm{DE},\theta \mathrm{pc}}\right\rangle$ relation, the statistical errors on the best-fit parameters (as quoted in Equation (21)) are very small. However, our estimates for $\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle$ and $\left\langle {P}_{\mathrm{DE},\theta \mathrm{pc}}\right\rangle$ do depend on several assumptions, including the CO-to-H₂ conversion factor, the geometry of the stellar disk, and the geometry of molecular gas structures on small scales. In Section 6, we investigate the systematic uncertainties associated with each of these assumptions.

The normalization of the best-fit $\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle$–$\left\langle {P}_{\mathrm{DE},\theta \mathrm{pc}}\right\rangle$ relation appears to be different (by a factor of 1.4) when estimated from data at different resolution (also see bottom panels in Figure 2). This can be partly attributed to the overall dependence of CO flux recovery fraction on resolution. As mentioned in Section 5.1, the sensitivity is often poorer at higher resolution, which means that only the brightest CO emission remains above the detection limit. For this reason, our estimated turbulent pressure at higher resolution suffers a stronger bias toward the brightest CO peaks, which trace high-density, high-pressure regions. This can qualitatively explain the mildly higher normalization of the $\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle$–$\left\langle {P}_{\mathrm{DE},\theta \mathrm{pc}}\right\rangle$ relation at 60 pc resolution.

5.2.3. Inferring the Dynamical State of Molecular Gas from the ISM Weight Budget

In our formulation (see Section 4.3 and Appendix A), all the terms contributing to the cloud-scale equilibrium pressure $\left\langle {P}_{\mathrm{DE},\theta \mathrm{pc}}\right\rangle$ can be grouped into three classes: (1) the weight of the cloud-scale molecular structures due to their self-gravity (${{ \mathcal W }}_{\mathrm{cloud}}^{\mathrm{self}};$ referred to as the "self-gravity term" hereafter); (2) the weight of these molecular structures due to the gravity associated with external material (including both ${{ \mathcal W }}_{\mathrm{cloud}}^{\mathrm{ext} \mbox{-} \mathrm{mol}}$ and ${{ \mathcal W }}_{\mathrm{cloud}}^{\mathrm{star}};$ referred to as the "external gravity terms"); and (3) the weight of the ambient atomic ISM in the combined potential created by stars and gas (${{ \mathcal W }}_{\mathrm{atom}};$ the "ambient pressure term"). The relative importance of these terms offers clues on a key question: which factor plays a more prominent role in governing the dynamical state of molecular structures like GMCs—is it self-gravity, external gravity, or ambient pressure?

The top panels in Figure 3 show the fractional contribution of the self-gravity term in the total $\left\langle {P}_{\mathrm{DE}}\right\rangle$ estimate, as a function of dynamical equilibrium pressure estimated on 120 pc scale ($\left\langle {P}_{\mathrm{DE},120\mathrm{pc}}\right\rangle$, top-left panel) and on kpc-scale (${P}_{\mathrm{DE},1\mathrm{kpc}}$, top-right panel). We find that the self-gravity term typically accounts for ∼33–70% of the total $\left\langle {P}_{\mathrm{DE}}\right\rangle$, and its fractional contribution exceeds 50% in about half of our disk sample. In the other half of our disk sample, the combination of the external gravity terms and the ambient pressure term dominate the self-gravity term. In this case, the dynamical states of molecular structures are strongly influenced by pressure in the ambient atomic ISM, and/or the gravitational potential created by stars and gas external to a given molecular structure.

Zoom In Zoom Out Reset image size

The black lines and blue shaded regions in Figure 3 represent the running median and 16–84th percentile trends. According to the trends shown in the top panels, the relative importance of the self-gravity term in the total ISM weight appears to increase with increasing $\left\langle {P}_{\mathrm{DE},120\mathrm{pc}}\right\rangle$ (rank correlation coefficient $\rho \,=0.62$, corresponding p-value ≪ 0.001), while it correlates less well with ${P}_{\mathrm{DE},1\mathrm{kpc}}$ ($\rho =0.12$, $p\ll 0.001$).

As discussed above, $\left\langle {P}_{\mathrm{DE},120\mathrm{pc}}\right\rangle$ reflects the pressure within the molecular gas (as seen by its tight correlation with $\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle$), whereas ${P}_{\mathrm{DE},1\mathrm{kpc}}$ represents the pressure in the kpc-scale environment (when neglecting the substructure in the molecular gas). The observed trends in the top two panels in Figure 3 can thus be interpreted as follows: across our sample, the dynamical state of cloud-scale molecular gas structures (hereafter "molecular structures") is strongly related to their internal pressure. Structures with high internal pressure ($\gtrsim {10}^{5}\,{k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$) are more likely to be dominated by self-gravity, whereas those with low internal pressure are more likely to be dominated by external gravity and/or ambient pressure. The large-scale environment pressure, however, offers less predicting power—the correlation is much less monotonic, and the chance of finding self-gravity–dominated molecular structures is about the same in low-pressure environments as in high-pressure environments within our sample.

The above discussion considers the relative importance of the self-gravity term in the total ISM weight budget. Alternatively, one could focus on the gravity felt by the molecular structures, and ask: "What fraction of the total weight of these structures (${{ \mathcal W }}_{\mathrm{cloud},\theta \mathrm{pc}};$ see Equation (16)) is due to their self-gravity, as opposed to external gravity?" To address this question, we plot the fractional contribution of the self-gravity term to the total weight of the cloud-scale molecular structures in the bottom panels in Figure 3. We find that the self-gravity term dominates the total internal weight of cloud-scale molecular structures in most (83%) of our disk sample. That is, in most cases, the observed molecular structures are dense enough to significantly alter the local gravitational potential. Moreover, in the cases when dynamical equilibrium holds and the ambient pressure is negligible, most of these molecular structures would be self-gravitating.

In the cases when self-gravity fails to outweigh external gravity, however, the molecular structures in question are likely not "significant" overdensities. We do not expect these molecular structures to be decoupled from large-scale dynamics, and if this remains the case, these structures might "dissolve" over roughly a galactic dynamical timescale. This picture is in line with recent findings by Chevance et al. (2020), namely that the lifetime of molecular clouds in some cases is driven by the timescales of galactic dynamical processes.

In the bottom panels of Figure 3, the contribution of self-gravity to the total weight of molecular structures shows a positive correlation with $\left\langle {P}_{\mathrm{DE},120\mathrm{pc}}\right\rangle$ (bottom-left; rank correlation coefficient $\rho =0.34$, corresponding p-value ≪ 0.001), and a very mild negative correlation with ${P}_{\mathrm{DE},1\mathrm{kpc}}$ (bottom-right; $\rho =-0.07$, p = 0.006). These trends appear to indicate that molecular structures with high internal pressure and/or in low-pressure environments are less likely to be dominated by external gravity.

We note that our estimates of the relative importance of molecular gas self-gravity may be biased by systematic effects related to sensitivity and (spatial) resolution. The finite sensitivity of the CO data might introduce a selection bias against molecular gas with low surface density in low-density, low-pressure environments (see discussion in Section 5.1). This selection bias offers a likely alternative explanation for the apparently high $\left\langle {{ \mathcal W }}_{\mathrm{cloud},120\,\mathrm{pc}}^{\mathrm{self}}\right\rangle /\left\langle {P}_{\mathrm{DE},120\mathrm{pc}}\right\rangle$ and $\left\langle {{ \mathcal W }}_{\mathrm{cloud},120\,\mathrm{pc}}^{\mathrm{self}}\right\rangle /\left\langle {{ \mathcal W }}_{\mathrm{cloud},120\mathrm{pc}}\right\rangle$ ratios at the low ${P}_{\mathrm{DE},1\mathrm{kpc}}$ end (Figure 3, right column).

On the other hand, the finite spatial resolution of the CO data means that we do not have access to the sub-beam density distribution of molecular gas. Our assumption of a uniform density sphere filling each beam could lead to underestimations of ${{\rm{\Sigma }}}_{\mathrm{mol}}$ and $\left\langle {{ \mathcal W }}_{\mathrm{cloud},120\,\mathrm{pc}}^{\mathrm{self}}\right\rangle$ if the actual sub-beam density distribution is strongly clumped. It is not trivial to predict how this bias would affect the trends we observe in Figure 3, because it remains unclear how different cloud/environment properties affect the clumping of molecular gas below these scales. In the future, CO observations with higher sensitivity and higher spatial resolution targeting low-pressure environments will help to eliminate these systematic effects.

5.2.4. A Physical Picture of Molecular Gas Dynamics on Cloud Scales

In this work, we adopt a metallicity-dependent ${\alpha }_{\mathrm{CO}}$ prescription (see Section 3), which is similar to the prescription suggested by Accurso et al. (2017, hereafter A17). Quite a few alternative prescriptions exist in the literature (e.g., Wolfire et al. 2010; Glover & Mac Low 2011; Feldmann et al. 2012; Narayanan et al. 2012; Bolatto et al. 2013). However, none of these prescriptions provides a concrete, observationally tested estimate of ${\alpha }_{\mathrm{CO}}$ that simultaneously captures the effects of metallicity, radiation field, and gas dynamics.

Our choice of ${\alpha }_{\mathrm{CO}}$ prescription could affect our results. Both ${{\rm{\Sigma }}}_{\mathrm{mol}}$ and ${P}_{\mathrm{turb}}$ are proportional to ${\alpha }_{\mathrm{CO}}$, and the molecular gas self-gravity term in ${P}_{\mathrm{DE}}$ is proportional to ${\alpha }_{\mathrm{CO}}^{2}$. To estimate the amount of systematic uncertainty associated with the choice of ${\alpha }_{\mathrm{CO}}$, we rederive our key measurements using three alternative ${\alpha }_{\mathrm{CO}}$ prescriptions, and compare them with our fiducial prescription. The three alternative prescriptions are: (1) a constant, Galactic conversion factor, (2) a simulation-based ${\alpha }_{\mathrm{CO}}$ calibration suggested by Narayanan et al. (2012, hereafter N12), and (3) an empirical ${\alpha }_{\mathrm{CO}}$ prescription suggested by Bolatto et al. (2013, hereafter B13).

For the constant ${\alpha }_{\mathrm{CO}}$ prescription, we use the Galactic value suggested by B13:

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

For the N12 prescription, we predict ${\alpha }_{\mathrm{CO}}$ in each kpc-sized aperture from the metallicity (Z') and the flux-weighted CO intensity ($\left\langle {I}_{\mathrm{CO}(1-0)}\right\rangle$), following their Equation (11):

$$\begin{eqnarray}\begin{array}{rcl}{\alpha }_{\mathrm{CO},{\rm{N}}12} & = & 8.5\ {M}_{\odot }\,{\mathrm{pc}}^{-2}{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1})}^{-1}\times {Z}^{{\prime} -0.65}\\ & & \times \min \left[1,\ 1.7\times {\left(\displaystyle \frac{\left\langle {I}_{\mathrm{CO}(1-0)}\right\rangle }{{\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{-0.32}\right]\,.\end{array}\end{eqnarray} \tag{ 23 }$$

We note that our implementation includes a factor of 1.36 correction for the mass of heavy elements, which was not included in the original N12 prescription. The quantity $\left\langle {I}_{\mathrm{CO}(1-0)}\right\rangle$ here is estimated from its CO (2−1) counterpart, $\left\langle {I}_{\mathrm{CO}(2-1),120\mathrm{pc}}\right\rangle$, assuming ${R}_{21}=0.7$.

For the B13 prescription, we predict ${\alpha }_{\mathrm{CO}}$ from the metallicity (Z'), typical GMC surface density (${{\rm{\Sigma }}}_{\mathrm{GMC}}$), and kpc-scale total surface density of both gas and stars (${{\rm{\Sigma }}}_{\mathrm{total}}$), following their Equation (31):

$$\begin{eqnarray}\begin{array}{rcl}{\alpha }_{\mathrm{CO},{\rm{B}}13} & = & \ 2.9\ {M}_{\odot }\,{\mathrm{pc}}^{-2}{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1})}^{-1}\\ & & \times \exp \left(\displaystyle \frac{0.4}{Z^{\prime} \,{{\rm{\Sigma }}}_{\mathrm{GMC}}^{100}}\right)\times {\left({{\rm{\Sigma }}}_{\mathrm{total}}^{100}\right)}^{-\gamma },\\ & & \mathrm{with}\ \gamma =\left\{\begin{array}{ll}0.5, & \mathrm{if}\ {{\rm{\Sigma }}}_{\mathrm{total}}^{100}\gt 1\\ 0. & \mathrm{otherwise}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 24 }$$

Here, ${{\rm{\Sigma }}}_{\mathrm{GMC}}^{100}$ and ${{\rm{\Sigma }}}_{\mathrm{total}}^{100}$ are the corresponding mass surface densities normalized to $100\ {M}_{\odot }\,{\mathrm{pc}}^{-2}$. Because neither of them could be derived from our observations without knowing ${\alpha }_{\mathrm{CO}}$ a priori, we set

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{GMC}}=\displaystyle \frac{{\alpha }_{\mathrm{CO}}}{{R}_{21}}\left\langle {I}_{\mathrm{CO}(2-1),120\mathrm{pc}}\right\rangle ,\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{total}}=\displaystyle \frac{{\alpha }_{\mathrm{CO}}}{{R}_{21}}\,{I}_{\mathrm{CO}(2-1),1\mathrm{kpc}}+{{\rm{\Sigma }}}_{{\rm{H}}{\rm{I}},1\mathrm{kpc}}+{{\rm{\Sigma }}}_{\mathrm{star},1\mathrm{kpc}},\end{eqnarray} \tag{ 26 }$$

and then iteratively solve for ${\alpha }_{\mathrm{CO}}$ in each kpc-scale aperture. We note that this iterative approach does not guarantee convergence, and 1.5% of the apertures in our sample do not yield a good solution. We discard the measurements for these apertures from this part of the analysis.

Figure 4 shows the probability density functions (PDFs) of the predicted ${\alpha }_{\mathrm{CO}}$ across our sample for each of the four prescriptions. Our fiducial prescription leads to an ${\alpha }_{\mathrm{CO}}$ distribution that peaks around the Galactic value. This agreement is largely by construction, as the fiducial prescription itself is normalized to the Galactic conversion factor at solar metallicity. The N12 prescription predicts comparatively higher ${\alpha }_{\mathrm{CO}}$ values, with the distribution peaking at around $8.5\ {M}_{\odot }\,{\mathrm{pc}}^{-2}{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1})}^{-1}$. This coincides with the "turning point" where the dependency on $\left\langle {I}_{\mathrm{CO}(1-0)}\right\rangle$ switches off (see Equation (23)). Therefore, the location of this peak is likely determined by the prescription itself rather than the input data. The B13 prescription produces a much wider ${\alpha }_{\mathrm{CO}}$ distribution compared to the other two distributions. This is attributable to the exponential term in Equation (24), which is a stronger dependence on metallicity than any of the other prescriptions. The B13 prescription also tends to predict higher-than-Galactic ${\alpha }_{\mathrm{CO}}$ values. This is likely driven by the relatively low ${{\rm{\Sigma }}}_{\mathrm{GMC}}$ values implied by Equation (25) (the median value in our sample is ${{\rm{\Sigma }}}_{\mathrm{GMC}}\approx 30\ {M}_{\odot }\,{\mathrm{pc}}^{-2}$).

Zoom In Zoom Out Reset image size

We demonstrate how our adopted ${\alpha }_{\mathrm{CO}}$ prescription affects our main conclusions in Figure 5. We show four versions of the $\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle$–$\left\langle {P}_{\mathrm{DE},120\mathrm{pc}}\right\rangle$ relation, each of which corresponds to a different ${\alpha }_{\mathrm{CO}}$ prescription. In all four panels, the data cluster around the line of equality. This is because both $\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle$ and $\left\langle {P}_{\mathrm{DE},120\mathrm{pc}}\right\rangle$ correlate positively with ${\alpha }_{\mathrm{CO}}$, and thus the choice of ${\alpha }_{\mathrm{CO}}$ prescription has less impact on their ratio. However, the choice of ${\alpha }_{\mathrm{CO}}$ prescription does have a more apparent impact on the absolute pressure values. As visible in Figure 5, the N12 and B13 prescriptions both push the whole distribution toward higher values of both pressures (also see Table 4 for quantitative results showing this trend). This is exactly what we would expect from the ${\alpha }_{\mathrm{CO}}$ PDFs: N12 predicts higher ${\alpha }_{\mathrm{CO}}$ on average, and thus higher pressure. B13 tends to predict higher ${\alpha }_{\mathrm{CO}}$ in disk regions, which pushes points near the low-pressure end up to higher pressures.

Zoom In Zoom Out Reset image size

For each prescription, we fit a power-law relation to all measurements from disk regions. We report the results in Table 4. All ${\alpha }_{\mathrm{CO}}$ prescriptions except B13 yield almost linear $\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle$–$\left\langle {P}_{\mathrm{DE},120\mathrm{pc}}\right\rangle$ relations. The B13 prescription instead leads to a superlinear slope of 1.34, with an estimated statistical uncertainty of ∼0.01. This is driven by the higher predicted ${\alpha }_{\mathrm{CO}}$ in low surface density and low metallicity environments.

In summary, adopting a different ${\alpha }_{\mathrm{CO}}$ prescription does not change the conclusion that the observed $\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle$–$\left\langle {P}_{\mathrm{DE},120\mathrm{pc}}\right\rangle$ distribution lies near the equality line. Adopting the B13 prescription makes its slope significantly steeper than linear, while adopting any of the other three prescriptions leads to nearly linear slopes. Adopting different prescriptions does significantly change the observed range of both pressures. These results illustrate the importance of quantifying ${\alpha }_{\mathrm{CO}}$ variations, and motivate future works to provide better constraints on the potential dependence of ${\alpha }_{\mathrm{CO}}$ on key physical properties, including metallicity, radiation field, gas (column) density, and dynamics.

Calculating ${P}_{\mathrm{DE},1\mathrm{kpc}}$ and $\left\langle {P}_{\mathrm{DE},\theta \mathrm{pc}}\right\rangle$ requires knowing the three-dimensional distribution of stars and gas in the galaxy disk. In Sections 4.2 and 4.3, we made a few assumptions to help us infer ${P}_{\mathrm{DE},1\mathrm{kpc}}$ and $\left\langle {P}_{\mathrm{DE},\theta \mathrm{pc}}\right\rangle$ from the observed two-dimensional projected quantities. In our fiducial approach, we assumed that (a) the stellar disk scale height ${H}_{\star }$ is proportional to the radial scale length ${R}_{\star }$, and is not a function of galactocentric radius (i.e., flat stellar disk); and (b) the mass-weighted ISM velocity dispersion is the relevant quantity that sets ISM disk scale height, which in turn sets ${P}_{\mathrm{DE}}$ near the disk midplane. The first assumption is partially motivated by the work by Kregel et al. (2002), and we provide more supporting evidence for this assumption in Appendix B. In this section, we explore the impact of modifying these two key assumptions.

Adopted Stellar Disk Geometry: The assumption of flat stellar disk geometry is widely used in previous works on stellar disk structure (e.g., van der Kruit & Searle 1981; Yoachim & Dalcanton 2006; Comerón et al. 2012), and supported by a recent observational study of edge-on disk galaxies (e.g., see Figure 12 in Comerón et al. 2011). We adopt this assumption as the fiducial choice in this paper (see Section 4.2).

An alternative, commonly considered possibility is a flared disk geometry (Yang et al. 2007; Ostriker et al. 2010). Here, we consider this alternative scenario, and explore whether adopting this alternative affects our conclusion. For this purpose, we re-evaluate ${\rho }_{\star }$ via

$$\begin{eqnarray}&&{\rho }_{\star ,1\,\mathrm{kpc}}^{\mathrm{flared}}=\displaystyle \frac{{{\rm{\Sigma }}}_{\star ,1\mathrm{kpc}}}{0.54{R}_{\star }}\ \exp \left(1-\displaystyle \frac{{r}_{\mathrm{gal}}}{{R}_{\star }}\right).\end{eqnarray} \tag{ 27 }$$

This assumes that the disk scale height flares exponentially at larger ${r}_{\mathrm{gal}}$, which is equivalent to assuming ${H}_{\star }\propto {{\rm{\Sigma }}}_{\star }^{-1}$ (corresponding to a constant stellar velocity dispersion; see Ostriker et al. 2010), and that ${{\rm{\Sigma }}}_{\star }$ drops exponentially as a function of ${r}_{\mathrm{gal}}$.

The top panel in Figure 6 shows the fractional deviation in the ${P}_{\mathrm{DE},1\mathrm{kpc}}$ estimates, when assuming a flared disk shape instead of a flat shape, as a function of ${P}_{\mathrm{DE},1\mathrm{kpc}}$. We find that assuming a flared disk geometry mainly leads to lower ${P}_{\mathrm{DE},1\mathrm{kpc}}$ at the low-pressure end. This trend makes sense, given the structure of galaxy disks. Low ${P}_{\mathrm{DE},1\mathrm{kpc}}$ generally corresponds to large ${r}_{\mathrm{gal}}$, and the flared disk shape will imply lower ${\rho }_{\star }$ in this regime. However, the amplitude of deviation in ${P}_{\mathrm{DE},1\mathrm{kpc}}$ is $\lesssim 0.2$ dex in most cases. This suggests that the deviation from a flat disk shape may lead to a factor of $\lesssim 1.6$ uncertainty on our ${P}_{\mathrm{DE},1\mathrm{kpc}}$ estimates.

Zoom In Zoom Out Reset image size

We also re-evaluate $\left\langle {P}_{\mathrm{DE},120\mathrm{pc}}\right\rangle$ assuming a flared stellar disk geometry. The corresponding $\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle$–$\left\langle {P}_{\mathrm{DE},120\mathrm{pc}}\right\rangle$ relation is quoted in Table 4. The changes in the 16th, 50th, and 84th percentiles of $\left\langle {P}_{\mathrm{DE}}\right\rangle$ show that the flared disk scenario gives lower $\left\langle {P}_{\mathrm{DE}}\right\rangle$ estimates compared to the fiducial scenario, and that this deviation is more significant at the low-pressure end. This leads to a slightly shallower slope (0.99) and a higher normalization (−0.08 dex at ${10}^{5}\,{k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$) for the $\left\langle {P}_{\mathrm{turb}}\right\rangle$–$\left\langle {P}_{\mathrm{DE}}\right\rangle$ relation. Nonetheless, the overall impact on the best-fit parameters is not large.

Therefore, the range of ${P}_{\mathrm{DE}}$ depends on our assumed stellar disk geometry, but the $\left\langle {P}_{\mathrm{turb}}\right\rangle$–$\left\langle {P}_{\mathrm{DE}}\right\rangle$ relation appears reasonably robust. In the near future, work using data from the Multi Unit Spectroscopic Explorer (MUSE; PI: E. Schinnerer) will provide direct measurements of the stellar velocity dispersion in a subset of our targets. This should help improve our knowledge of the three-dimensional stellar disk structure in these targets.

Adopted Gas Velocity Dispersion: When calculating the kpc-scale equilibrium pressure in Section 4.2, we treat the entire ISM as a single component fluid. This motivates us to use the mass-weighted velocity dispersion combining atomic and molecular gas for estimating ${P}_{\mathrm{DE},1\mathrm{kpc}}$ (Equation (14)). Many previous studies have instead adopted a fixed velocity dispersion of ${\sigma }_{\mathrm{gas},{\rm{z}}}\approx 8\mbox{--}11\ \mathrm{km}\,{{\rm{s}}}^{-1}$ (e.g., Blitz & Rosolowsky 2004, 2006; Leroy et al. 2008; Ostriker et al. 2010; Hughes et al. 2013a), which is about the mean observed value for atomic gas at moderate galactocentric radii in nearby galaxy disks (Leroy et al. 2008; Tamburro et al. 2009; Caldú-Primo et al. 2013; Mogotsi et al. 2016).

We compare our approach with the fixed ${\sigma }_{\mathrm{gas},{\rm{z}}}$ approach, again by comparing their corresponding ${P}_{\mathrm{DE},1\mathrm{kpc}}$ estimates. As shown in the lower panel in Figure 6, assuming a fixed ${\sigma }_{\mathrm{gas},{\rm{z}}}=10\ \mathrm{km}\,{{\rm{s}}}^{-1}$ generally leads to higher ${P}_{\mathrm{DE},1\mathrm{kpc}}$ in disk regions relative to our fiducial estimates. This is because the observed molecular gas velocity dispersion at 60–120 pc scales is usually less than $10\ \mathrm{km}\,{{\rm{s}}}^{-1}$ in the disk regions in our sample, and generally smaller at larger ${r}_{\mathrm{gal}}$ (or at low pressure). Nevertheless, the resulting deviation in ${P}_{\mathrm{DE},1\mathrm{kpc}}$ values is again within 0.2 dex in most cases.

We note that our $\left\langle {P}_{\mathrm{DE},\theta \mathrm{pc}}\right\rangle$ estimates treat the molecular and atomic gas separately, and thus do not rely directly on the mass-weighted velocity dispersion (see Section 4.3). Therefore, the discussion above does not apply to these cloud-scale estimates.

In Sections 3 and 4, we adopt an approach that treats the gas in each pixel separately. This approach, which we refer to as the "pixel statistics approach," preserves information from the smallest recoverable scale. We use this approach to derive mean cloud-scale gas properties (e.g., $\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle$) across our sample.

Another popular approach is to segment the observed gas distribution into regions that likely correspond to coherent physical objects. For example, many cloud segmentation algorithms such as CLUMPFIND (Williams et al. 1994) and CPROPS (Rosolowsky & Leroy 2006) group individual voxels into cloud-like objects associated with local maxima. Using the voxels associated with each cloud, one can derive cloud size, velocity dispersion, and luminosity, as well as other higher-order properties.

From such a cloud catalog, one can derive the mass-weighted mean pressure and similar quantities for clouds in each region of each galaxy. We refer to this segmentation-based calculation as the "cloud statistics" approach. The cloud statistics approach accesses the same physics as our pixel statistics measurements, but differs in two important ways. First, the cloud-based approach treats identified objects, rather than resolution elements, as the fundamental structural unit. Second, the cloud-based approach yields size measurements for each of these objects. Most of the cloud literature assumes spherical symmetry, i.e., the projected size of the objects on the sky is assumed to reflect the depth of the objects along the line of sight.

Our data allow both approaches. In this section, we rederive our key measurements using a cloud statistics approach, and then compare the results to those from our pixel statistics measurements.

We use the PHANGS-ALMA CPROPS cloud catalogs, which are derived from the same CO data set that we use (A. Hughes et al. 2020, in preparation; E. Rosolowsky et al. 2020, in preparation). Similar to our calculations, these begin with a set of fixed, 120 pc physical resolution data. The algorithm identifies significant, independent local maxima and then associates emission with each maximum. Next, it measures the size, CO luminosity, and line width associated with each cloud using moment methods, and corrects for biases due to finite sensitivity and resolution. We note that the PHANGS-ALMA CPROPS application uses the "seeded CLUMPFIND" assignment option in CPROPS. This assigns all significant emission to nearby local maxima, and so represents a hybrid between the default CPROPS assignment and the CLUMPFIND assignment schemes. Otherwise, the calculations, as detailed in E. Rosolowsky et al. (2020, in preparation), follow the original CPROPS approach.

Within each kpc-sized aperture, we derive the CO-flux-weighted average turbulent pressure and dynamical equilibrium pressure for all objects identified by CPROPS. To do this, we first write down the expressions for these two quantities using cloud mass, size, and velocity dispersion (i.e., the CPROPS counterparts of Equations (10) and (15)):

$$\begin{eqnarray}&&\left\langle {P}_{\mathrm{turb},\mathrm{CPROPS}}\right\rangle =\left\langle \displaystyle \frac{3{M}_{\mathrm{cloud}}{\sigma }_{\mathrm{cloud}}^{2}}{4\pi {R}_{\mathrm{cloud}}^{3}}\right\rangle ,\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}\begin{array}{rcl}\left\langle {P}_{\mathrm{DE},\mathrm{CPROPS}}\right\rangle & = & \,\left\langle \displaystyle \frac{3\,{{GM}}_{\mathrm{cloud}}^{2}}{8\pi {R}_{\mathrm{cloud}}^{4}}\right\rangle \\ & & \ +\left\langle \displaystyle \frac{{{GM}}_{\mathrm{cloud}}}{2{R}_{\mathrm{cloud}}^{2}}\right\rangle {{\rm{\Sigma }}}_{\mathrm{mol},1\mathrm{kpc}}\\ & & \ +\left\langle \displaystyle \frac{3\,{{GM}}_{\mathrm{cloud}}}{2{R}_{\mathrm{cloud}}}\right\rangle {\rho }_{\star ,1\mathrm{kpc}}\\ & & \ +{{ \mathcal W }}_{\mathrm{atom},1\mathrm{kpc}}.\end{array}\end{eqnarray} \tag{ 29 }$$

The "$\left\langle \right\rangle$" symbol here denotes a CO-flux-weighted average over all CPROPS clouds that have their central coordinate inside the kpc-sized aperture in question.

When substituting the measured cloud parameters from CPROPS into Equations (28) and (29), we pay special attention to two caveats. First, we adopt the same metallicity-dependent conversion factor for the cloud-based analysis as we do for the pixel statistics (see Section 3.3.1). Second, the cloud radius quoted in the CPROPS catalogs is defined as 1.91 times the one-dimensional rms size calculated based on the object's projected intensity distribution on the sky (i.e., following the Solomon et al. 1987 convention to account for clouds being centrally condensed). To enforce better consistency between the cloud and pixel measurements, we convert the radius quoted by CPROPS (${R}_{\mathrm{CPROPS}}$) to the radius of a hypothesized, constant-density spherical cloud via

$$\begin{eqnarray}&&{R}_{\mathrm{cloud}}=\sqrt{5}\ \displaystyle \frac{{R}_{\mathrm{CPROPS}}}{1.91}=1.17\,{R}_{\mathrm{CPROPS}}.\end{eqnarray} \tag{ 30 }$$

Here, the factor of $\sqrt{5}$ is the ratio between the radius of a spherical, constant-density cloud and its projected rms size on the sky (see Equations (11)–(13) in Rosolowsky & Leroy 2006).

Using the ${\alpha }_{\mathrm{CO}}$-corrected cloud mass, the adjusted radius, and the measured velocity dispersion, we derive estimates of $\left\langle {P}_{\mathrm{turb},\mathrm{CPROPS}}\right\rangle$ and $\left\langle {P}_{\mathrm{DE},\mathrm{CPROPS}}\right\rangle$ via Equations (28) and (29). The left panel in Figure 7 shows the relation between these two quantities across our sample. We find that almost all data points lie below the equality line. That is, using the cloud statistics approach and assuming spherical symmetry for the objects, we find ubiquitously lower $\left\langle {P}_{\mathrm{turb},\mathrm{CPROPS}}\right\rangle$ than $\left\langle {P}_{\mathrm{DE},\mathrm{CPROPS}}\right\rangle$. On average, this offset is about 0.66 dex.

Zoom In Zoom Out Reset image size

To understand this apparent discrepancy between the results from cloud statistics and from pixel statistics, we look into the actual measured sizes of the objects in the CPROPS catalogs. With Equation (30) applied, the median value of estimated object diameters across PHANGS-ALMA is $\left\langle 2{R}_{\mathrm{cloud}}\right\rangle \approx 400$ pc, or about 3 times (∼0.5 dex) larger than the beam size (which is 120 pc in this case).

These apparently large cloud sizes are not completely unexpected. Just like many other segmentation algorithms designed to find "clumps," CPROPS tends to recover structures with sizes comparable to or larger than the beam size. This effect has long been noticed and discussed by many previous works (e.g., Verschuur 1993; Hughes et al. 2013a; Leroy et al. 2016). These objects may be real physical structures (e.g., giant molecular associations or filaments). However, the assumption of spherical symmetric is unlikely to hold, because in this case the ∼400 pc diameters are much larger than the ∼100 pc vertical FWHM of the Milky Way molecular gas disk (see Heyer & Dame 2015).

As an ad hoc correction, we rederive the values of $\left\langle {P}_{\mathrm{turb},\mathrm{CPROPS}}\right\rangle$ and $\left\langle {P}_{\mathrm{DE},\mathrm{CPROPS}}\right\rangle$, assuming a modified object geometry. We still use ${R}_{\mathrm{cloud}}$ as the projected size of the object on the sky, but now we assume the line-of-sight depth of the object to be 120 pc. This matches the assumption used for the pixel statistics estimates. This effectively assumes a cylindrical geometry for the identified objects, with their projected shapes on the sky kept the same, but their depth fixed to a constant value. In practice, this means that we derive the cloud surface density via ${{\rm{\Sigma }}}_{\mathrm{cloud}}={M}_{\mathrm{cloud}}/(\pi {R}_{\mathrm{cloud}}^{2})$, and use the cloud surface density, velocity dispersion, and a fixed ${D}_{\mathrm{cloud}}=120$ pc in Equations (10), (15), and (16) to estimate $\left\langle {P}_{\mathrm{turb},\mathrm{CPROPS}}\right\rangle$ and $\left\langle {P}_{\mathrm{DE},\mathrm{CPROPS}}\right\rangle$.

The right panel in Figure 7 shows the relation between the "corrected" $\left\langle {P}_{\mathrm{turb},\mathrm{CPROPS}}\right\rangle$ and $\left\langle {P}_{\mathrm{DE},\mathrm{CPROPS}}\right\rangle$ estimates. In contrast to the results shown in the left panel, we find much better agreement between these "corrected" pressure estimates. The best-fit power-law relation (see Table 4) also becomes much more consistent with the results derived from the pixel statistics approach. These findings suggest that, compared to the spherical symmetry assumption, the assumption of a fixed 120 pc line-of-sight depth is a much better description of the actual geometry of the CPROPS identified objects in the PHANGS-ALMA CO maps.

In summary, the $\left\langle {P}_{\mathrm{turb}}\right\rangle$–$\left\langle {P}_{\mathrm{DE}}\right\rangle$ relation derived from the cloud statistics approach shows consistency with the pixel statistics results, provided that one adopts an appropriate assumption for the geometry of the identified objects. We also emphasize that, when analyzing data with marginal spatial resolution, extra caution should be used when interpreting results of cloud identification algorithms like CPROPS.

In the self-regulated star formation model of Ostriker et al. (2010) and Ostriker & Shetty (2011), an actively star-forming disk is viewed as a (quasi-)steady-state system. Stellar feedback in the form of radiation, winds, and supernovae offsets (in a time-averaged sense) losses of energy and pressure due to cooling and turbulent dissipation. At the same time, pressure maintains (again, in a time-averaged sense) support for the gas against collapse in the gravitational field, balancing both gas self-gravity and external disk gravity. However, a small fraction of the interstellar gas can collapse and form stars locally in regions where support against gravity is insufficient. The result of this localized collapse is the star formation feedback that pressurizes the rest of the ISM, providing internal support. Although this is likely a violent process with alternating episodes of collapse and expansion (see, e.g., Kruijssen et al. 2019; Rahner et al. 2019; Schinnerer et al. 2019; Chevance et al. 2020), dynamical equilibrium is expected when considering the entire ISM across large spatial scales (≫ typical size of GMCs and star-forming regions) and long timescales (≫ typical lifetime of SF cycle). Numerical simulations indeed show that a well-defined quasi-steady-state exists on spatial scales of order ∼1 kpc and timescales of a few hundred Myr, even though the SFR and the ISM properties strongly vary on short space- and timescales (see, e.g., Kim & Ostriker 2017, 2018; Semenov et al. 2017, 2018; Orr et al. 2018).

In this framework, turbulent pressure (${P}_{\mathrm{turb}}$) and dynamical equilibrium pressure (${P}_{\mathrm{DE}}$) are both closely related to a third variable—the local SFR surface density, ${{\rm{\Sigma }}}_{\mathrm{SFR}}$. Here, ${P}_{\mathrm{turb}}$ should be directly proportional to ${{\rm{\Sigma }}}_{\mathrm{SFR}}$, given a fixed momentum injection per star formed (e.g., Ostriker & Shetty 2011). Simultaneously, the equilibrium-state ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ is determined by requiring the sum of turbulent, thermal, and magnetic pressures (individually proportional to ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ due to feedback) to balance ${P}_{\mathrm{DE}}$ (e.g., Ostriker et al. 2010; Kim et al. 2011). In the following subsections, we explore this scenario by showing the relationships between ${P}_{\mathrm{turb}}$ and ${{\rm{\Sigma }}}_{\mathrm{SFR}}$, and between ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ and ${P}_{\mathrm{DE}}$.

7.1.1. SFR Surface Density versus Turbulent Pressure

The left panel in Figure 8 shows the measured cloud-scale average turbulent pressure in molecular gas, ${\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}$, as a function of the kpc-scale average SFR surface density, ${{\rm{\Sigma }}}_{\mathrm{SFR},1\mathrm{kpc}}$. We observe a strong correlation between these two quantities. For all measurements in disk regions (blue dots), we find a rank correlation coefficient of ρ = 0.71 (corresponding p-value ≪ 0.001).

Zoom In Zoom Out Reset image size

In the case that the ISM pressure is feedback-driven, the ratio between ISM turbulent pressure and SFR surface density—where these are averages computed over the same area—reflects to the momentum injection per unit mass of stars formed, ${p}_{\star }/{m}_{\star }$, via

$$\begin{eqnarray}&&{P}_{\mathrm{turb}}=\displaystyle \frac{1}{4}\displaystyle \frac{{p}_{\star }}{{m}_{\star }}\,{{\rm{\Sigma }}}_{\mathrm{SFR}}.\end{eqnarray} \tag{ 31 }$$

The prefactor of 1/4 assumes spherical expansion sites centered on the disk midplane and that the momentum flux to upper and lower halves of the ISM disk translates directly to ISM turbulent pressure (Ostriker & Shetty 2011). While the above picture is idealized and should be modified by details of turbulent injection and dissipation, its prediction agrees with the measured relationship between ${P}_{\mathrm{turb}}$ and ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ in disk simulations of the star-forming multiphase ISM (Kim et al. 2013; Kim & Ostriker 2015a, 2017).

The ratio ${p}_{\star }/{m}_{\star }$ is predicted to range between 10³–${10}^{4}\ \mathrm{km}\,{{\rm{s}}}^{-1}$ for supernova feedback, depending on the ISM properties, spatial and temporal clustering of supernovae, and energy losses due to interface mixing (e.g., Iffrig & Hennebelle 2015; Kim & Ostriker 2015b; Martizzi et al. 2015; Walch & Naab 2015; Kim et al. 2017; El-Badry et al. 2019; Gentry et al. 2019). In Figure 8, we show the predicted ${P}_{\mathrm{turb}}$–${{\rm{\Sigma }}}_{\mathrm{SFR}}$ relation for a range of ${p}_{\star }/{m}_{\star }$ values (dotted lines). Our observed $\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle$–${{\rm{\Sigma }}}_{\mathrm{SFR},1\mathrm{kpc}}$ relation has a normalization that would correspond to a high momentum injection.

One possible explanation for this apparently high momentum injection rate is the clumping of the star formation distribution. Unlike $\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle$, which is estimated on 120 pc scales and then averaged over the kpc-scale aperture, our ${{\rm{\Sigma }}}_{\mathrm{SFR},1\mathrm{kpc}}$ measurements are derived directly on kpc scale. Just like the molecular gas, we expect star formation to cluster on sub-kpc scales (e.g., Grasha et al. 2018, 2019; Schinnerer et al. 2019; Chevance et al. 2020). This will cause the ${{\rm{\Sigma }}}_{\mathrm{SFR},1\mathrm{kpc}}$ values in Figure 8 to appear lower due to the inclusion of area without star formation, and thus it underestimates the actual SFR surface density relevant to feedback momentum injection.

To account for this issue, we introduce a dimensionless prefactor ${C}_{\mathrm{fb}}\gtrsim 1$, which corrects for the artificial dilution of ${{\rm{\Sigma }}}_{\mathrm{SFR},1\mathrm{kpc}}$ compared to ${{\rm{\Sigma }}}_{\mathrm{SFR},\theta \mathrm{pc}}$. Then we have

$$\begin{eqnarray}&&{\left\langle {P}_{\mathrm{turb},120\mathrm{pc}}\right\rangle }_{1\mathrm{kpc}}={C}_{\mathrm{fb}}\,\displaystyle \frac{{p}_{\star }}{4{m}_{\star }}\,{{\rm{\Sigma }}}_{\mathrm{SFR},1\mathrm{kpc}}.\end{eqnarray} \tag{ 32 }$$

The median value of ${C}_{\mathrm{fb}}\,{p}_{\star }/{m}_{\star }$ implied by measurements in our disk sample is $6.3\times {10}^{3}\ \mathrm{km}\,{{\rm{s}}}^{-1}$, and its 16th–84th percentile range is 0.73 dex.

A direct estimation of ${C}_{\mathrm{fb}}$ requires high-resolution data tracing SFR on matched ∼10–100 pc scales. This is currently not available for our entire sample (c.f., Kreckel et al. 2018; Schinnerer et al. 2019; Chevance et al. 2020). However, we can obtain a rough estimate of ${C}_{\mathrm{fb}}$ based on our high-resolution CO data, if we assume that the clustering of star formation matches that of the molecular gas. Given a median area-filling fraction of $\approx 1/3$ in our CO data, we expect ${C}_{\mathrm{fb}}\approx 3$. In this case, most of the $\left\langle {P}_{\mathrm{turb},\theta \mathrm{pc}}\right\rangle$ and ${{\rm{\Sigma }}}_{\mathrm{SFR},1\mathrm{kpc}}$ measurements in the disk sample are consistent with ${p}_{\star }/{m}_{\star }={10}^{3}\mbox{--}{10}^{4}\ \mathrm{km}\,{{\rm{s}}}^{-1}$.

There are many caveats related to this approximate estimation of ${C}_{\mathrm{fb}}$. The spatial and temporal distribution of gas and star formation could change the amount of momentum injected into the gas, and the spatial scales on which most momentum is deposited also matters. Given these complications, a logical next step on this topic will be a detailed, multiscale comparison between observations of gas and star formation at high resolution, as well as numerical simulations with realistic feedback prescriptions.

An alternative explanation for the high ${P}_{\mathrm{turb}}/{{\rm{\Sigma }}}_{\mathrm{SFR}}$ ratio measurements in Figure 8 (left panel) is that other turbulence driving mechanisms might also play a role, at least within a subset of our sample. For example, many measurements in bulge or bar regions show the highest ${P}_{\mathrm{turb}}$ at fixed ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ (orange symbols in Figure 8, left panel). While part of this can be explained by beam-smearing effects, we also expect bar-induced radial inflow and the corresponding conversion of gravitational potential energy to be an additional source of turbulence in the ISM (e.g., through shocks near the center of the stellar bar; Binney et al. 1991; Sormani et al. 2015). More quantitative comparisons between our results and other turbulent driving mechanisms will be carried out in future works.

7.1.2. Equilibrium Pressure versus SFR Surface Density

The right panel in Figure 8 shows kpc-scale SFR surface density, ${{\rm{\Sigma }}}_{\mathrm{SFR},1\mathrm{kpc}}$, as a function of kpc-scale dynamical equilibrium pressure, ${P}_{\mathrm{DE},1\mathrm{kpc}}$. We find a strong correlation between the two quantities for all measurements in disk regions. This supports the idea that, in steadily star-forming disks, the three-dimensional star and gas distribution on kpc scales largely determines the average SFR surface density on the same scale. In any part of the disk where the ISM weight is higher, the steady-state SFR must also increase for feedback to maintain the pressure that matches the dynamical equilibrium pressure.

To quantify this observed correlation, we derive a best-fit power-law relation for all disk measurements, using the OLS bisector method in logarithmic space:

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{SFR}}}{{10}^{-3}\ {M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-2}}=3.2{\left(\displaystyle \frac{{P}_{\mathrm{DE}}}{{10}^{4}{k}_{{\rm{B}}}{\rm{K}}{\mathrm{cm}}^{-3}}\right)}^{0.84}.\end{eqnarray} \tag{ 33 }$$

The rms scatter in ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ around this relation is 0.20 dex. The statistical errors estimated from bootstrapping are 0.01 for the slope and 0.01 dex for the normalization, but we expect systematic errors (on both ${{\rm{\Sigma }}}_{\mathrm{SFR},1\mathrm{kpc}}$ and ${P}_{\mathrm{DE},1\mathrm{kpc}}$) to be a larger source of uncertainty on the fitting results.

The same correlation has been observed in various types of galaxies by many previous studies (e.g., Leroy et al. 2008; Genzel et al. 2010; Ostriker et al. 2010; Ostriker & Shetty 2011; Herrera-Camus et al. 2017; Fisher et al. 2019). To put our measurements into context, we also include in Figure 8 two other data sets derived also from observations of nearby galaxies. Leroy et al. (2008; brown crosses) measure this relation in a sample of 23 nearby galaxies, with each independent measurement representing one kpc-wide radial bin in a galaxy. Herrera-Camus et al. (2017; green plus symbols) investigate the same relation in the H i-dominated regions in 31 KINGFISH galaxies (Kennicutt et al. 2011), treating each kpc-sized region independently.

These previous works adopted slightly different assumptions when deriving ${P}_{\mathrm{DE}}$ from observables. To make a fair comparison between our measurements and theirs, our ${P}_{\mathrm{DE},1\mathrm{kpc}}$ measurements shown in Figure 8 are instead derived by assuming ${\sigma }_{\mathrm{gas},{\rm{z}}}=11\ \mathrm{km}\,{{\rm{s}}}^{-1}$, following exactly the same prescription in Leroy et al. (2008) and Herrera-Camus et al. (2017). We find overall consistency, but our measurements concentrate more toward the high-${P}_{\mathrm{DE}}$, high-${{\rm{\Sigma }}}_{\mathrm{SFR}}$ end. This is because the PHANGS-ALMA CO observations primarily target the star-forming inner part of high-mass disk galaxies, while the Herrera-Camus et al. (2017) sample focuses on the H i–dominated outer disks on purpose, and the Leroy et al. (2008) sample includes many dwarf galaxies.

In addition to previous observational results, we also show in Figure 8 the predicted ${P}_{\mathrm{DE}}$–${{\rm{\Sigma }}}_{\mathrm{SFR}}$ relation from a hydrodynamic simulation (Kim et al. 2013; black dashed line in the right panel):

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{SFR}}}{{10}^{-3}\ {M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-2}}=1.8{\left(\displaystyle \frac{{P}_{\mathrm{DE}}}{{10}^{4}{k}_{{\rm{B}}}{\rm{K}}{\mathrm{cm}}^{-3}}\right)}^{1.13}.\end{eqnarray} \tag{ 34 }$$

This relation has a steeper slope than our best-fit power-law relation (Equation (33)), while it agrees with previous observations in the low ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ regime. We point out that the simulations in Kim et al. (2013) cover a ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ range from 10⁻⁴ to ${10}^{-2}\ {M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-2}$, whereas our sample covers from 10⁻³ to ${10}^{-1}\ {M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-2}$. The shallower slope found in our sample might then reflect some systematic change in the properties of the ISM in massive galaxies and inner disk environments, relative to the ISM in dwarf galaxies and/or outer disk environments.

A shallower ${P}_{\mathrm{DE}}$–${{\rm{\Sigma }}}_{\mathrm{SFR}}$ relation in high-${{\rm{\Sigma }}}_{\mathrm{SFR}}$ environments is seen in several previous works. In the Leroy et al. (2008) sample, some hint of a shallower ${{\rm{\Sigma }}}_{\mathrm{SFR}}$–${P}_{\mathrm{DE}}$ relation is visible from the high-${{\rm{\Sigma }}}_{\mathrm{SFR}}$ measurements. For a sample of local ultra luminous infrared galaxies (ULIRGs) and high-redshift star-forming galaxies, Ostriker & Shetty (2011) found a ${{\rm{\Sigma }}}_{\mathrm{SFR}}$–${P}_{\mathrm{DE}}$ relation with a slope of 0.95 at ${P}_{\mathrm{DE}}\gt {10}^{5}\,{k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$, which is closer to our result.³³ More recently, Fisher et al. (2019) report a much shallower slope of 0.77 at ${P}_{\mathrm{DE}}\gtrsim {10}^{5}\,{k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$ for a sample of local turbulent disk galaxies,³⁴ which are believed to resemble typical star-forming galaxies at z ∼ 1–2. Future studies on the ISM in more local ULIRGs and galaxy centers can provide better constraints on the slope of the ${{\rm{\Sigma }}}_{\mathrm{SFR}}$–${P}_{\mathrm{DE}}$ relation at the high ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ end, and potentially help unveil the physics regulating star formation in the "starburst" regime (see, e.g., Thompson et al. 2005; Ostriker & Shetty 2011; Shetty & Ostriker 2012; Crocker et al. 2018; Krumholz et al. 2018).

Following early suggestions by Elmegreen (1989), the ISM dynamical equilibrium pressure has also been viewed as a determinant of the molecular/atomic phase balance in the ISM. In this scenario, ${P}_{\mathrm{DE}}$ relates closely to the molecular-to-atomic gas ratio, ${R}_{\mathrm{mol}}$, and thus influences the fraction of the ISM in the dense, star-forming phase. This idea has been tested by many subsequent works (e.g., Wong & Blitz 2002; Blitz & Rosolowsky 2006; Leroy et al. 2008), and is commonly adopted as a prescription for determining the molecular gas fraction in semianalytic models of galaxy evolution (e.g., Lagos et al. 2018). Here, we report the observed scaling relation between molecular-to-atomic ratio and the ISM dynamical equilibrium pressure across our sample.

Following previous studies, we use the kpc-scale dynamical equilibrium pressure ${P}_{\mathrm{DE},1\mathrm{kpc}}$ to trace the average ambient pressure in the ISM. To allow a quantitative comparison with previous results, we again use the ${P}_{\mathrm{DE},1\mathrm{kpc}}$ values estimated by assuming a fixed ${\sigma }_{\mathrm{gas},{\rm{z}}}=11\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (similar to Section 7.1.2). We determine the molecular-to-atomic gas ratio from the ratio of our measured molecular and atomic gas surface densities on kpc scales:

$$\begin{eqnarray}&&{R}_{\mathrm{mol},1\mathrm{kpc}}\equiv {{\rm{\Sigma }}}_{\mathrm{mol},1\mathrm{kpc}}/{{\rm{\Sigma }}}_{\mathrm{atom},1\mathrm{kpc}}.\end{eqnarray} \tag{ 35 }$$

We show the relation between ${P}_{\mathrm{DE},1\mathrm{kpc}}$ and ${R}_{\mathrm{mol},1\mathrm{kpc}}$ in Figure 9. Our sample spans nearly two orders of magnitude in ${R}_{\mathrm{mol},1\mathrm{kpc}}$, with most measurements clustering around or above the atomic-to-molecular transition threshold (i.e., ${R}_{\mathrm{mol},1\mathrm{kpc}}\,=1$). We find a positive and statistically significant correlation between ${R}_{\mathrm{mol},1\mathrm{kpc}}$ and ${P}_{\mathrm{DE},1\mathrm{kpc}}$ across our whole sample (Spearman's rank correlation coefficient $\rho =0.58;$ corresponding p-value ≪ 0.001). This strong, positive correlation is qualitatively consistent with previous observations (Wong & Blitz 2002; Blitz & Rosolowsky 2006; Leroy et al. 2008), even though the adopted CO-to-H₂ conversion factor, stellar mass-to-light ratio, and stellar disk geometry vary among studies.

Zoom In Zoom Out Reset image size

We perform an OLS bisector fit on all our disk measurements over the range $0.1\lt {R}_{\mathrm{mol},1\mathrm{kpc}}\lt 10$, following Leroy et al. (2008). This yields a best-fit power-law relation (blue solid line in Figure 9) of

$$\begin{eqnarray}&&{R}_{\mathrm{mol},1\mathrm{kpc}}={\left(\displaystyle \frac{{P}_{\mathrm{DE},1\mathrm{kpc}}}{2.1\times {10}^{4}{k}_{{\rm{B}}}{\rm{K}}{\mathrm{cm}}^{-3}}\right)}^{1.02}.\end{eqnarray} \tag{ 36 }$$

The scatter in ${R}_{\mathrm{mol},1\mathrm{kpc}}$ around this relation is 0.36 dex. The formal statistical errors in the fit are small: 0.02 for the slope α, and 0.01 dex for the threshold pressure ${P}_{0}=2.1\,\times {10}^{4}\,{k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3}$ (at which the ISM transitions from being predominantly atomic to molecular, or vice versa). However, when varying the choice of ${\alpha }_{\mathrm{CO}}$ prescriptions and other assumptions, the best-fit α and P₀ appear systematically uncertain by ∼0.20 and ∼0.15 dex, respectively. Given this level of systematic uncertainty, Equation (36) quantitatively agrees with the ${R}_{\mathrm{mol}}$–${P}_{\mathrm{DE},1\mathrm{kpc}}$ relations reported in previous works (α = 0.73–1.05, ${P}_{0}\,=(1.5\mbox{--}4.5)\times {10}^{4}\,{k}_{{\rm{B}}}\,{\rm{K}}\,{\mathrm{cm}}^{-3};$ see Wong & Blitz 2002; Blitz & Rosolowsky 2006; Leroy et al. 2008, also see black lines in Figure 9).

The weight of a cloud, ${{ \mathcal W }}_{\mathrm{cloud}}$, includes three constituent parts:

$$\begin{eqnarray}&&{{ \mathcal W }}_{\mathrm{cloud}}={{ \mathcal W }}_{\mathrm{cloud}}^{\mathrm{self}}+{{ \mathcal W }}_{\mathrm{cloud}}^{\mathrm{ext} \mbox{-} \mathrm{mol}}+{{ \mathcal W }}_{\mathrm{cloud}}^{\mathrm{star}}.\end{eqnarray} \tag{ A1 }$$

These three terms represent the weights due to the cloud's own self-gravity (${{ \mathcal W }}_{\mathrm{cloud}}^{\mathrm{self}}$), to the gravity of external molecular gas outside the cloud (${{ \mathcal W }}_{\mathrm{cloud}}^{\mathrm{ext} \mbox{-} \mathrm{mol}}$), and to the gravity of stars (${{ \mathcal W }}_{\mathrm{cloud}}^{\mathrm{star}}$). By symmetry, the outer atomic gas layers exert no gravity on any structure in the molecular layer.

For the cloud self-gravity term, we integrate the weight of a sphere with radius ${R}_{\mathrm{cloud}}$ and constant density³⁹ ${\rho }_{\mathrm{cloud}}$:

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal W }}_{\mathrm{cloud}}^{\mathrm{self}} & = & {\int }_{0}^{{R}_{\mathrm{cloud}}}{\rho }_{\mathrm{cloud}}\displaystyle \frac{G{\rho }_{\mathrm{cloud}}(4/3)\pi {r}^{3}}{{r}^{2}}\,{dr}\\ & = & \displaystyle \frac{2\pi }{3}G{\rho }_{\mathrm{cloud}}^{2}{R}_{\mathrm{cloud}}^{2}\\ & = & \displaystyle \frac{3\pi }{8}G{{\rm{\Sigma }}}_{\mathrm{cloud}}^{2}.\end{array}\end{eqnarray} \tag{ A2 }$$

The last step re-expresses ${{ \mathcal W }}_{\mathrm{cloud}}^{\mathrm{self}}$ in terms of the cloud surface density ${{\rm{\Sigma }}}_{\mathrm{cloud}}\equiv {M}_{\mathrm{cloud}}/(\pi {R}_{\mathrm{cloud}}^{2})=(4/3){\rho }_{\mathrm{cloud}}{R}_{\mathrm{cloud}}$.

For the external molecular gas gravity term, we approximate the distribution of all the molecular gas outside the cloud as a slab centered at the disk midplane. This slab has constant volume density ${\bar{\rho }}_{\mathrm{mol}}$ anywhere outside the cloud, and has zero density within the extent of the cloud (i.e., it has a spherical "hole"). For simplicity, we further assume that the vertical half width of this slab is equal to the cloud radius, ${H}_{\mathrm{mol}}={R}_{\mathrm{cloud}}$. The weight of the cloud due to the gravity of this component is thus:

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal W }}_{\mathrm{cloud}}^{\mathrm{ext} \mbox{-} \mathrm{mol}} & = & {{ \mathcal W }}_{\mathrm{cloud},\mathrm{slab}}-{{ \mathcal W }}_{\mathrm{cloud},\mathrm{hole}}\\ & = & {\int }_{0}^{{R}_{\mathrm{cloud}}}{\rho }_{\mathrm{cloud}}\,(4\pi G{\bar{\rho }}_{\mathrm{mol}}z){dz}\\ & & -\,{\int }_{0}^{{R}_{\mathrm{cloud}}}{\rho }_{\mathrm{cloud}}\displaystyle \frac{G{\bar{\rho }}_{\mathrm{mol}}(4/3)\pi {r}^{3}}{{r}^{2}}\,{dr}\\ & = & 2\pi G{\bar{\rho }}_{\mathrm{mol}}{\rho }_{\mathrm{cloud}}{R}_{\mathrm{cloud}}^{2}-\displaystyle \frac{2\pi }{3}G{\bar{\rho }}_{\mathrm{mol}}{\rho }_{\mathrm{cloud}}{R}_{\mathrm{cloud}}^{2}\\ & = & \displaystyle \frac{4\pi }{3}G{\bar{\rho }}_{\mathrm{mol}}{\rho }_{\mathrm{cloud}}{R}_{\mathrm{cloud}}^{2}\\ & = & \displaystyle \frac{\pi }{2}G{\bar{{\rm{\Sigma }}}}_{\mathrm{mol}}{{\rm{\Sigma }}}_{\mathrm{cloud}}.\end{array}\end{eqnarray} \tag{ A3 }$$

The last step re-expresses ${{ \mathcal W }}_{\mathrm{cloud}}^{\mathrm{ext} \mbox{-} \mathrm{mol}}$ in terms of the cloud surface density ${{\rm{\Sigma }}}_{\mathrm{cloud}}$ and the average surface density of the slab ${\bar{{\rm{\Sigma }}}}_{\mathrm{mol}}={\bar{\rho }}_{\mathrm{mol}}(2{H}_{\mathrm{mol}})={\bar{\rho }}_{\mathrm{mol}}(2{R}_{\mathrm{cloud}})$.

For the stellar gravity term, given that ${R}_{\mathrm{cloud}}\ll {H}_{\star }$, we treat the stellar disk as having a uniform density near the midplane, ${\rho }_{\star }$. We thus have:

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal W }}_{\mathrm{cloud}}^{\mathrm{star}} & = & {\int }_{0}^{{R}_{\mathrm{cloud}}}{\rho }_{\mathrm{cloud}}\,(4\pi G{\rho }_{\star }z){dz}\\ & = & 2\pi G{\rho }_{\star }{\rho }_{\mathrm{cloud}}{R}_{\mathrm{cloud}}^{2}\\ & = & \displaystyle \frac{3\pi }{2}G{\rho }_{\star }{{\rm{\Sigma }}}_{\mathrm{cloud}}{R}_{\mathrm{cloud}}.\end{array}\end{eqnarray} \tag{ A4 }$$

The weight of the atomic gas layer, ${{ \mathcal W }}_{\mathrm{atom}}$, also consists of three parts:

$$\begin{eqnarray}&&{{ \mathcal W }}_{\mathrm{atom}}={{ \mathcal W }}_{\mathrm{atom}}^{\mathrm{self}}+{{ \mathcal W }}_{\mathrm{atom}}^{\mathrm{mol}}+{{ \mathcal W }}_{\mathrm{atom}}^{\mathrm{star}}.\end{eqnarray} \tag{ A5 }$$

These terms represent the weights due to the atomic layer's self-gravity (${{ \mathcal W }}_{\mathrm{atom}}^{\mathrm{self}}$), to the gravity of the inner molecular gas layer (${{ \mathcal W }}_{\mathrm{atom}}^{\mathrm{mol}}$), and to the gravity of stars (${{ \mathcal W }}_{\mathrm{atom}}^{\mathrm{star}}$). The calculation of the first and last terms is very similar to that of ${P}_{\mathrm{DE},1\mathrm{kpc}}$ in Section 4.2. We simply quote the results here:

$$\begin{eqnarray}&&{{ \mathcal W }}_{\mathrm{atom}}^{\mathrm{self}}=\displaystyle \frac{\pi }{2}G{{\rm{\Sigma }}}_{\mathrm{atom}}^{2}\end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray}&&{{ \mathcal W }}_{\mathrm{atom}}^{\mathrm{star}}={{\rm{\Sigma }}}_{\mathrm{atom}}\sqrt{2G{\rho }_{\star }}\,{\sigma }_{\mathrm{atom},{\rm{z}}}.\end{eqnarray} \tag{ A7 }$$

Note that the second expression adopts the assumption that ${H}_{\mathrm{mol}}\ll {H}_{\mathrm{atom}}$.

The molecular gravity term, however, is slightly different from the version in the ${P}_{\mathrm{DE},1\mathrm{kpc}}$ calculation. As the entire molecular gas inner layer is assumed to be "sandwiched" by the atomic gas outer layer in this calculation, the latter feels the full gravity of the former, such that:

$$\begin{eqnarray}&&{{ \mathcal W }}_{\mathrm{atom}}^{\mathrm{mol}}=\pi G{\bar{{\rm{\Sigma }}}}_{\mathrm{mol}}{{\rm{\Sigma }}}_{\mathrm{atom}}.\end{eqnarray} \tag{ A8 }$$

In reality, the molecular and atomic medium (or at least the cold atomic phase) are often well-mixed, so this estimate represents an upper limit for the true ${{ \mathcal W }}_{\mathrm{atom}}^{\mathrm{mol}}$ value.

Combining all the above derivations together, we have:

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal W }}_{\mathrm{total}} & = & \,{{ \mathcal W }}_{\mathrm{cloud}}+{{ \mathcal W }}_{\mathrm{atom}}\\ & = & \,{{\rm{\Sigma }}}_{\mathrm{cloud}}\left(\displaystyle \frac{3\pi }{8}G{{\rm{\Sigma }}}_{\mathrm{cloud}}+\displaystyle \frac{\pi }{2}G{\bar{{\rm{\Sigma }}}}_{\mathrm{mol}}+\displaystyle \frac{3\pi }{2}G{\rho }_{\star }{R}_{\mathrm{cloud}}\right)\\ & & \,+{{\rm{\Sigma }}}_{\mathrm{atom}}\left(\displaystyle \frac{\pi }{2}G{{\rm{\Sigma }}}_{\mathrm{atom}}+\pi G{\bar{{\rm{\Sigma }}}}_{\mathrm{mol}}+\sqrt{2G{\rho }_{\star }}\,{\sigma }_{\mathrm{atom},{\rm{z}}}\right).\end{array}\end{eqnarray} \tag{ A9 }$$

This gives the estimated total weight of a molecular cloud and the atomic gas above it in a galaxy disk, or equivalently, the required pressure in this molecular cloud to keep it under dynamical equilibrium.