Molecular Cloud Populations in the Context of Their Host Galaxy Environments: A Multiwavelength Perspective

We present a rich, multiwavelength, multiscale database built around the PHANGS-ALMA CO$\,$(2-1) survey and ancillary data. We use this database to present the distributions of molecular cloud populations and sub-galactic environments in 80 PHANGS galaxies, to characterize the relationship between population-averaged cloud properties and host galaxy properties, and to assess key timescales relevant to molecular cloud evolution and star formation. We show that PHANGS probes a wide range of kpc-scale gas, stellar, and star formation rate (SFR) surface densities, as well as orbital velocities and shear. The population-averaged cloud properties in each aperture correlate strongly with both local environmental properties and host galaxy global properties. Leveraging a variable selection analysis, we find that the kpc-scale surface densities of molecular gas and SFR tend to possess the most predictive power for the population-averaged cloud properties. Once their variations are controlled for, galaxy global properties contain little additional information, which implies that the apparent galaxy-to-galaxy variations in cloud populations are likely mediated by kpc-scale environmental conditions. We further estimate a suite of important timescales from our multiwavelength measurements. The cloud-scale free-fall time and turbulence crossing time are ${\sim}5{-}20$ Myr, comparable to previous cloud lifetime estimates. The timescales for orbital motion, shearing, and cloud-cloud collisions are longer, ${\sim}100$ Myr. The molecular gas depletion time is $1{-}3$ Gyr and shows weak to no correlations with the other timescales in our data. We publish our measurements online and expect them to have broad utility to future studies of molecular clouds and star formation.

1. INTRODUCTION Molecular clouds are deeply integrated with their host galaxies by a number of intertwined physical processes. The gas distribution, gravitational potential, radiation field, and feedback-driven flows in the host galaxy regulate molecular cloud formation and evolution (Dobbs et al. 2014;Ballesteros-Paredes et al. 2020). The internal structure and dynamical properties of these clouds in turn set the initial conditions for star formation, which over time reshapes the matter and radiation distribution in the galaxy (McKee & Ostriker 2007;Padoan et al. 2014;Klessen & Glover 2016;Girichidis et al. 2020). These complex interactions lead to strong, observable correlations between the properties of molecular clouds and the local and global properties of the host galaxy. Characterizing these cloud-environment correlations is thus a promising avenue for understanding the physics governing molecular cloud evolution and, consequently, star formation and galaxy evolution.
Observations of molecular clouds in our Galaxy and a number of nearby galaxies have identified various empirical trends manifesting such cloud-environment correlations. Within a galaxy, molecular clouds located closer to the galaxy center appear denser, more massive, and more turbulent (e.g., Oka et al. 2001;Colombo et al. 2014;Freeman et al. 2017;Hirota et al. 2018;Miura et al. 2018;Brunetti et al. 2021, also see Heyer & Dame 2015). Similar trends have been found in galaxy-scale numerical simulations (e.g., Pan et al. 2015;Jeffreson et al. 2020;Treß et al. 2021). Recent observational works also report that more massive and actively star-forming galaxies tend to host clouds with typically larger sizes, masses, surface densities, and velocity dispersions (Hughes et al. 2013a;Leroy et al. 2015Leroy et al. , 2016Schruba et al. 2019;Sun et al. 2020a, but see Bolatto et al. 2008;Fukui & Kawamura 2010;Donovan Meyer et al. 2013).
To proceed from the existing empirical knowledge to more concrete understandings of the cloud-environment correlations, major advances on two issues are necessary. First, the characterization of environmental dependence often stops at a qualitative level in the molecular cloud literature: the "environments" are commonly defined in crude, categorical ways (e.g., galaxy centers, stellar bars, spiral arms), and they are merely considered as a secondary, moderating factor on the scaling relations followed by molecular clouds. To better understand the underlying physics, a more direct approach would be to quantify the dependence of molecular cloud properties on a set of quantitative "environmental metrics," such as the local gas and stellar mass surface density, star formation rate (SFR) surface density, and orbital shear (e.g., Hughes et al. 2013a;Schruba et al. 2019). Second, many previous works (especially earlier ones) had to rely on observations in a small number of galaxies or sub-galactic regions, and thus only probed a limited range of host galaxy properties. While such case studies could yield unique insights for specific targets, only systematic surveys across large galaxy samples can cover a wide, continuous range of host galaxy properties, produce representative population statistics, and make meaningful connections to galaxy evolution models.
The PHANGS-ALMA survey (Leroy et al. 2021a) was designed to address both of these issues. This survey provides sensitive, high resolution, wide fieldof-view CO (2-1) imaging data for ∼90 nearby, highmass, star-forming galaxies. With these galaxies sampled uniformly along the star-forming main sequence, PHANGS-ALMA enables systematic studies of giant molecular clouds (GMCs; M 10 5 M ) across an array of environments where most stars form in the local universe. Furthermore, a rich set of multiwavelength ancillary data furnishes a multifaceted depiction of these host galaxies, making it possible to study their molecular cloud populations in full environmental context. Indeed, one of the core science goals that motivated the PHANGS-ALMA survey was to characterize the dependence of molecular cloud populations on global and local galaxy properties. Studies on this data set have presented population statistics for key molecular cloud properties such as mass, size, surface density, velocity dispersion, and virial parameter (Sun et al. , 2020aRosolowsky et al. 2021, A. Hughes et al. in preparation). They also noted significant variations among galaxies and across morphological regions within galaxies (e.g., centers, bars). In the direct predecessor of this paper, Sun et al. (2020b) conducted a joint analysis on the PHANGS-ALMA CO data and multiwavelength ancillary data. They showed that the variations in molecular gas turbulent pressure can be attributed to the dynamical balance between gravity and internal/external pressure in the gas, as previously argued in Galactic and extragalactic molecular cloud studies (e.g., Field et al. 2011;Hughes et al. 2013a;Schruba et al. 2019).
In this paper, we directly address this core science goal of the PHANGS-ALMA survey. We build on a cross-spatial-scale analysis framework used by Sun et al. (2020b) and calculate population statistics for the molecular cloud properties measured in Sun et al. (2018Sun et al. ( , 2020a; Rosolowsky et al. (2021); and A. Hughes et al. (in preparation). We further cross-match them with a large suite of environmental metrics depicting the local gas and stellar mass distribution, orbital kinematics, morphological structures, and star formation activities in the host galaxy. This allows us to (1) present the full range of cloud populations and host galaxy environments captured by PHANGS-ALMA, (2) delineate the quantitative relationships between cloud characteristics (e.g., mass, surface density, velocity dispersion) and environmental metrics (e.g., gas, star, and SFR surface densities), and (3) identify a subset of relationships that carry unique explanatory/predictive power among all the observed cloud-environment correlations.
Another goal of this paper is to present a set of machine-readable data tables that consolidate all the aforementioned measurements from PHANGS-ALMA and ancillary surveys. These high-level data products have already been used in a number of studies. Herrera et al. (2020) and Barnes et al. (2021) have utilized previous versions of these tables to quantitatively assess several physical mechanisms relevant to molecular cloud and H ii region evolution, including stellar feedback and pressure balance. Several ongoing works also rely on these tables to calculate the star formation efficiency and its link to small-scale turbulence, orbital shear, and disk instabilities (E. Rosolowsky et al. in preparation; J. Sun et al. in preparation; T. Williams et al. in preparation). In this paper, we also utilize the measurements in these tables to calculate a suite of characteristic timescales related to the gravity, turbulent motions, orbital motions, and star formation rate of molecular clouds. The ratios among these timescales provide unique constraints on the viable mechanisms regulating molecular cloud evolution and star formation (also see e.g., Wong 2009;Jeffreson & Kruijssen 2018;Kruijssen et al. 2019a;Chevance et al. 2020a,b;Kim et al. 2021a).
The structure of this paper is as follows. Section 2 describes our galaxy sample and the sources of all data we use. Section 3 elaborates the cross-spatial-scale analysis framework we use to assemble the multiwavelength measurements into a coherent data structure. Section 4 presents the distribution of various molecular cloud population properties and sub-galactic environmental properties measured in this work, and Section 5 characterizes the correlations between these two types of measurements. To demonstrate an application scenario of our rich multiwavelength measurements, Section 6 presents a set of characteristic timescales relevant to molecular cloud evolution and star formation. Finally, Section 7 summarizes all our findings.
2. DATA In this paper, we focus on a sample of 80 galaxies (see Table A1) selected from the full PHANGS-ALMA survey sample (90 galaxies; see Leroy et al. 2021a). We select these galaxies according to two criteria: (1) their PHANGS-ALMA CO (2-1) observations have beam full-width-half-maximum (FWHM) sizes corresponding to physical scales of 150 pc or smaller, so that each beam roughly probes a GMC-sized molecular gas structure; and (2) they are not too heavily inclined (i 75 • ), so that we can unambiguously determine the locations of molecular clouds in the host galaxy. A subset of 28 galaxies in this sample already appeared in Sun et al. (2020b), where they utilized earlier versions of the same observational data sets and data analysis infrastructure.
In addition to the PHANGS-ALMA CO data, our target galaxies have abundant multiwavelength coverage, including radio, mid-/near-infrared (MIR/NIR), opti-cal, and near-/far-ultraviolet (NUV/FUV) data (Sections 2.2, 2.3, 2.4, and 2.7). High-level measurements such as CPROPS object catalogs created from CO data cubes (Rosolowsky et al. 2021), rotation curves derived from CO line kinematics (Lang et al. 2020), and morphological feature masks constructed from near-IR images (Querejeta et al. 2021) are also available for most targets (Sections 2. 1.1, 2.5, and 2.6). These rich ancillary data provide us with comprehensive information about the multiphase ISM, stellar disk structures, star formation, and galactic dynamical properties on kpc scales for our target galaxies.
In the following subsections, we detail the sources and characteristics of all raw data and high-level data products used in this study. We provide a schematic summary of these input data at the top of Figure 1.

PHANGS-ALMA CO Data
We use the PHANGS-ALMA CO (2-1) imaging data 1 (Leroy et al. 2021a) to probe molecular gas properties on 150 pc scales. These data cover the actively starforming area in each galaxy (∼100 kpc 2 on average) and have sufficient depth and resolving power to detect and isolate the CO emission from individual GMCs (with a typical mass of 10 5 M ). They include both interferometric and single dish observations and thus provide sensitivity to emission on all spatial scales. We refer interested readers to Leroy et al. (2021a) for more details regarding sample selection, observational setup, and data product characteristics, and Leroy et al. (2021b) for an in-depth description of data calibration, imaging, and product creation procedures.
In this study, we measure properties of molecular cloud populations from the PHANGS-ALMA CO data using two different approaches. The first approach measures molecular gas properties "object-by-object." In this case, the objects of interest are identified by applying the cloud segmentation algorithm CPROPS to the PHANGS-ALMA CO data cubes (Rosolowsky & Leroy 2006;Rosolowsky et al. 2021, A. Hughes et al. in preparation). The second approach treats the molecular gas as a spatially continuous medium and extracts measurement in a "pixel-by-pixel" fashion directly from the PHANGS-ALMA CO line moment maps, where the beam size corresponds roughly to the typical size of an individual GMC or giant molecular association (Leroy et al. 2016;Sun et al. 2018Sun et al. , 2020a.
While the two approaches access similar physical properties and often lead to similar results (see Sun et al. 2020a;Rosolowsky et al. 2021), they complement each other in important ways. The object-by-object approach treats each identified object as a fundamental structural unit, and by providing size estimates for these H⍺ Line ( Figure 1. Schematics of the data sources, aggregation methods, and the derived physical quantities. For easier navigation within this paper, we also note the section number relevant to each input data product and output quantities.
objects it probes the spatial organization of molecular gas. The pixel-by-pixel method instead treats each resolution element as a fundamental unit, which preserves information from the smallest recoverable scale while remaining agnostic about the organization of the gas on larger scales. When presenting our measurements in Section 4, we compare the two approaches to illustrate how methodological choices could influence the main results.

Object-by-object Measurements
We extract a set of molecular gas measurements for each object tabulated in the PHANGS-ALMA CPROPS catalogs 2 (Rosolowsky et al. 2021, A. Hughes et al. in preparation). These catalogs are constructed by running CPROPS (Rosolowsky & Leroy 2006) on data cubes at a set of common spatial resolutions (60, 90, 120, and 150 pc whenever available). Note that A. Hughes et al. (in preparation) present two versions of CPROPS catalogs for the PHANGS-ALMA sample: one is constructed from data cubes whose noise levels are homogenized among all galaxies, and the other from data cubes with the native noise. We use the latter version in this paper, because (a) we would like to compare the object-by-object measurements to the corresponding pixel-bypixel ones, which were derived from the original data cubes without noise homogenization; and (b) we would like to recover as much CO emission above the noise floor as possible.
For each object identified by the CPROPS algorithm, the catalog records its integrated CO line luminosity, L CO, obj , CO line width, σ CO, obj , and the two-dimensional projected radius on the sky 3 , R 2D, obj . These numbers are calculated after radially extrapolating each object to a hypothetical boundary at 0 K brightness temperature and then deconvolving the beam size and channel width. From these basic observables, we estimate the following physical properties for each object: • Molecular gas mass, M obj . This is derived from the integrated CO (2-1) line luminosity L CO, obj (in units of K km s −1 pc 2 ) via Here R 21 = 0.65 is the adopted CO (2-1) to CO (1-0) line ratio (den Brok et al. 2021;Leroy et al. 2021c), and α CO(1-0) is a varying CO-to-H 2 conversion factor for the CO (1-0) line. By default, we use a metallicitydependent α CO prescription as described in Sun et al. (2020b): α CO(1-0) = 4.35 Z −1.6 M (K km s −1 pc 2 ) −1 . (2) Here, Z is the inferred local gas phase abundance normalized to the solar value (see Section 3.3). While we use Equation 2 as our fiducial prescription, we also calculate α CO using several alternative prescriptions and include them in the data products (see Appendix B).
• Molecular gas surface density, Σ obj . This is derived from M obj and R 2D, obj via Following Rosolowsky et al. (2021), this estimate assumes a two-dimensional Gaussian profile for the projected gas mass distribution and focuses on the area within a FWHM, which includes half of the total gas mass. The cos i term accounts for galaxy inclination (see Table A1) by correcting the derived surface densities to face-on projection. This correction was not present in the formulae used by Rosolowsky et al. (2021) and previous similar work. We motivate this correction in Appendix C.
• One-dimensional velocity dispersion, σ obj . This is derived from the measured CO line width σ CO, obj via Here σ CO, obj comes from the second moment (moment-2) of the CO line profile and is corrected for broadening due to the line spread function (LSF; see Rosolowsky et al. 2021). The extra (cos i) 0.5 term is an empirically determined correction that accounts for the dependence of the observed CO line width on galaxy inclination. Appendix C details the origin of this term and the rationale for including it here. In the following discussions, we will assume that σ obj is dominated by turbulent motions, though this measurement will include additional contributions from thermal or ordered, streaming motions. Fundamentally, it reflects the velocity dispersion along the line of sight direction within each object.
• Three-dimensional mean radius, R obj . This quantity is inferred from R 2D, obj via Here H = 100 pc is an assumed molecular gas disk thickness perpendicular to the galaxy plane (Heyer & Dame 2015), and H cos i would be the expected lineof-sight depth given the disk inclination. Equation 5 assumes a spheroidal geometry when the object diameter on the sky exceeds this line-of-sight depth, and a spherical geometry otherwise. This is similar to the treatments in Rosolowsky et al. (2021), except that here we also correct for galaxy inclination (also see Appendix C). Our adopted value for H is likely uncertain by a factor of ∼ 2 due to variations within a galaxy and among galaxies (e.g., Yim et al. 2014;Bacchini et al. 2019). Systematic trends with galactocentric radius and global galaxy mass are also expected. A fixed value of H = 100 pc cannot capture these variations, which means that our inferred R obj values (and any measurements that rely on them, see bullet points below) are affected accordingly. Nevertheless, the functional form of Equation 5 suggests that at most 1/3 of the fractional uncertainty on H will propagate to R obj , which would only be marginally significant in comparison to other sources of systematic uncertainties (see discussions in Section 4.2 and Appendix E).
• Turbulent pressure, P turb, obj . This is derived from M obj , σ obj , and R obj via Here the mean density ρ obj is derived from R obj and the mass within the FWHM of a two-dimensional Gaussian profile (Equation 16 in Rosolowsky et al. 2021).
• Virial parameter, α vir, obj . This is derived from M obj , σ obj , and R obj via This formula is derived by calculating the kinetic energy (E kin ) and gravitational potential energy (E grav ) for the gas within the two-dimensional FWHM size, assuming a uniform density distribution (consistent with Rosolowsky et al. 2021). With this definition, a virialized object would have α vir, obj = 1, whereas an object in energy equipartition would have α vir, obj = 2. But we note that the virial parameter estimated in this way might not be a complete description of cloud dynamical states if there are strong magnetic field, surface pressure, or external tidal forces (see discussions in, e.g., Ballesteros-Paredes 2006;Sun et al. 2020b;Kim et al. 2021b;Liu et al. 2021).

Pixel-by-pixel Measurements
As an alternative to the object-by-object approach, we also derive molecular gas properties at fixed spatial resolutions (60, 90, 120, and 150 pc whenever available) measured pixel-by-pixel from the PHANGS-ALMA CO moment maps (see Sun et al. 2018Sun et al. , 2020a. The PHANGS-ALMA data reduction pipeline (Leroy et al. 2021b) produces two versions of moment maps: a "broad" version that prioritizes high CO flux completeness through highly inclusive signal masking, and a "strict" version that features high signal-to-noise (S/N) CO line moment measurements thanks to more restrictive masking. For our pixel-by-pixel analysis, we primarily use the "strict" moment maps so that only pixels with reliable measurements are included in our calculation. To account for the lower CO flux completeness of the "strict" maps, we later estimate their flux completeness by comparing the "strict" and "broad" maps, and further correct for sensitivity-induced biases (see Section 3.2.2).
The "strict," beam-matched moment maps provide CO line integrated intensity (moment-0), I CO, pix , CO line effective width (see Heyer et al. 2001), σ CO, pix , and their associated uncertainties for every pixel with detected CO emission. From these basic observables, we derive a set of molecular gas physical properties mirroring those from the object-by-object approach, and estimate their statistical uncertainties through Gaussian error propagation: • Molecular gas surface density, Σ pix . This quantity is derived from the integrated CO (2-1) line intensity I CO, pix (in units of K km s −1 ) for each pixel, via Here R 21 and α CO(1-0) represent the adopted CO line ratio and CO (1-0)-to-H 2 conversion factor as in Equation 1. The same cos i inclination correction from Equation 3 also applies here.
• Molecular gas mass, M pix . We also record the total molecular gas mass captured in each beam 4 via Here A beam = (π/4 ln 2)D 2 beam is the effective area of the beam with a FWHM of D beam (i.e., 60, 90, 120, or 150 pc). No inclination correction is required here since both I CO, pix and A beam are measured/defined in the projected plane of the sky.
• One-dimensional velocity dispersion, σ pix . This quantity is derived from the LSF-corrected CO line width σ CO, pix in each pixel and uses the same inclination correction as Equation 4: Here σ CO, pix represents the CO line effective width, which is a different line width metric than the one based on the second moment used in Equation 4. The effective width is a more robust line width metric than moment-2 at low S/N, but it could give biased results when there are multiple velocity components along the line of sight (Henshaw et al. 2020).
• Three-dimensional mean radius, R pix . We adopt the following three-dimensional size for the gas structure captured in each beam, mirroring Equation 5: Again, the cos i term accounts for galaxy inclination by converting the molecular gas disk thickness (perpendicular to the galaxy plane) to the depth along the line of sight. Note that Equation 11 yields a single R pix value for each given beam size.
• Turbulent pressure, P turb, pix . This is derived from Σ pix and σ pix via This assumes that the gas mass captured in each beam is uniformly distributed within a radius of R pix . This is consistent with the geometrical assumptions adopted in previous studies (e.g., Sun et al. 2020b), yet it leads to an inconsistency with the object-based approach (Equation 6). We comment on this issue in Section 2.1.3.
• Virial parameter, α vir, pix . This is derived from Σ pix and σ pix via This also assumes a spherical geometry and a uniform density distribution within R pix . Similar to the situation with our turbulent pressure estimates, the geometrical assumptions here are not fully consistent with those adopted for the object-based analysis (see Section 2.1.3 for further comments).

Notes on the Common Grounds and Differences between the Object-based and Pixel-base Approaches
The object-based and pixel-based approaches show an apparent symmetry, in the sense that they have many measured quantities in common, such as molecular gas surface density, velocity dispersion, turbulent pressure, and virial parameter. This allows us to make direct comparisons between the two approaches and assess how our methodological choices might influence the quantitative results. However, it is worth emphasizing that, for several reasons, we do not necessarily expect the two approaches to yield exactly the same quantitative results.
First and foremost, the two approaches are motivated by two slightly different views for the structure and geometry of the molecular ISM in galaxies. The objectbased approach views the molecular ISM as a collection of dense, centrally concentrated structures, and the central goal of the CPROPS algorithm is to segment the observed CO emission distribution such that each identified CO-emitting object corresponds to a coherent structure like a GMC or a giant molecular association. The pixel-based approach instead views the molecular ISM as a continuous distribution of gas while being agnostic about its spatial clustering, and the measurement process simply characterizes the gas captured in each beam. In a sense, the two approaches see the same observational data through different lenses, and each attempts to extract measurable properties in a way that is most consistent with its adopted view.
Reflecting these different views, there are also important, practical differences in the methodologies between these two approaches, which make it non-trivial to draw direct comparisons between them. In particular, the object-based approach aims to measure the true size and mass of each identified object by deconvolving the beam and extrapolating the detected part of each object to a hypothetical boundary at 0 K brightness temperature. Such operations could in principle account for biases due to the finite resolution and sensitivity of the observations, but they are implicitly model-dependent and not easily adaptable to fit the pixel-based approach.
The distinct physical models underlying these two approaches are also reflected in the different auxiliary assumptions they adopt when calculating physical quantities. The object-based approach assumes compact, Gaussian-shaped gas distributions and calculate gas surface density, turbulent pressure, and virial parameter for only the half of the gas located within the Gaussian FWHM (Equations 3-7, consistent with Rosolowsky et al. 2021). In contrary, the pixel-based approach considers all the gas mass detected in each beam and assumes it is uniformly distributed within the beam area (Equations 12-13, in line with Sun et al. 2018Sun et al. , 2020a.
Considering these complications, we do not necessarily expect the two approaches to agree in their quantitative results, even though we start from the same CO data cubes and attempt to define measurable properties in a symmetric way.

H i Data
We use interferometric H i 21 cm line data to trace the distribution of neutral atomic gas in each galaxy. These include both new and archival observations taken by the Karl G. Jansky Very Large Array (VLA) and the Australia Telescope Compact Array (ATCA).
Among the 53 galaxies with H i data (see Table A1), 20 have been observed as part of the PHANGS-VLA survey (A. Sardone et al. in preparation). The other galaxies have archival data from either large nearby galaxy surveys such as THINGS (nine galaxies; Walter et al. 2008), VIVA (six galaxies; Chung et al. 2009), HERA-CLES (four galaxies; Leroy et al. 2009), LVHIS (three galaxies; Koribalski et al. 2018), EveryTHINGS (two galaxies; I. Chiang et al. in preparation), or individual case studies with the VLA (seven galaxies) and the ATCA (two galaxies; Murugeshan et al. 2019). These H i data sets have typical angular resolution of 15 −35 (16−84 percentile), which corresponds to linear scales of 0.7−2.8 kpc (see Section 3.3 for further discussions about H i data resolution). The 3σ sensitivity limit ranges 10−100 K km s −1 for the H i line intensity.
Assuming optically thin 21 cm emission, we convert 21 cm line intensity I 21cm to atomic gas surface density Σ atom via Here Σ atom includes the (extra 35%) mass of helium and heavier elements. The cos i term accounts for galaxy inclination.

Near-IR Data
We use near-IR imaging data from the Spitzer Space Telescope and the Wide-field Infrared Survey Explorer (WISE) to trace the old stellar mass distribution (see Table A1, column 8). For 61 galaxies in our sample, we use Spitzer IRAC 3.6 µm images from the S 4 G survey (Sheth et al. 2010). For those without S 4 G data, we instead use WISE W1 band (3.4 µm) images compiled by the z0MGS project . All these data are postprocessed by subtracting background emission, masking foreground stars in the field of view, and convolving the non-Gaussian point spread function (PSF) to a 7.5 Gaussian PSF using appropriate convolution kernels (Aniano et al. 2011).
We convert the stellar continuum intensity at 3.4 µm and 3.6 µm to stellar mass surface density, Σ , locationby-location via Here, Υ 3.4 µm is the stellar mass-to-light (M/L) ratio at 3.4 µm, which should be nearly identical to that at 3.6 µm. We adopt a spatially varying M/L ratio, which was estimated by Leroy et al. (2021a) for all PHANGS-ALMA targets based on an empirical relation between Υ 3.4 µm and the local SFR surface density to 3.4 µm surface brightness ratio.

Mid-IR and UV Data
We use mid-IR images from WISE and far-/near-UV images from the Galaxy Evolution Explorer (GALEX) to trace the distribution of obscured and unobscured star formation. These data are also compiled by the z0MGS project ) and have been postprocessed by subtracting background emission, masking foreground stars, reprojecting to a shared astrometry, and then convolving to a 15 Gaussian PSF.
We combine the mid-IR and UV data and calculate the local star formation rate (SFR) surface density following the prescriptions described in Leroy et al. (2021a). By default we use the combination of GALEX FUV (154 nm) and WISE 22 µm data to calculate the local SFR surface density Σ SFR M yr −1 kpc −2 = 8.9 × 10 −2 I 154 nm MJy sr −1 + 3.0 × 10 −3 I 22 µm MJy sr −1 cos i .
These prescriptions assume a Chabrier initial mass function (IMF; Chabrier 2003) via their calibration against Salim et al. (2016), which is also consistent within ≈ 5% with calibrations using a Kroupa IMF (Kroupa 2001). The quantitative results agree with extinctioncorrected Hα-based SFR estimates from the PHANGS-MUSE survey (Emsellem et al. 2021, F. Belfiore et al. in preparation) at a ∼20−30% level overall, but there is divergence in low SFR regions due to contributions from IR cirrus and/or old stellar populations (Boquien et al. 2016). We refer the reader to Leroy et al. (2021a) for more details on the calibration of these SFR prescriptions.

Rotation Curves
We use rotation curves derived from CO line kinematics by Lang et al. (2020) to characterize galactic orbital kinematics (e.g., orbital period and shear) locally within each galaxy. These rotation curves are measured from the same PHANGS-ALMA CO data set, and therefore cover roughly the same galactocentric radius range as the CO maps themselves. They are available for 62 out of the 80 galaxies.
The rotation curves in Lang et al. (2020) are measured and recorded with finite radial bin sizes (∼150 pc). Due to the sparse distribution of CO detections across the field of view and the likely presence of unaccounted local streaming motions in the gas, the measured circular velocity sometimes fluctuates considerably between adjacent radial bins. These bin-to-bin fluctuations make it challenging to reliably estimate any parameter that depends on the derivative of the rotation curves.
To address this issue, we use a set of functional fitting models constructed from the measured rotation curves (J. Nofech et al. in preparation) rather than the raw measurements themselves. These fitting models adopt the "universal rotation curve" functional form suggested by Persic et al. (1996). The fitting process effectively forces the rotation curve models to be smooth and have physically sensible slopes (i.e., with its logarithmic derivative between -0.5 and 1), while still matching the actual measurements as closely as possible. We visually inspect all fitting results and conclude that the models represent the raw measurements reasonably well.
Based on these best-fit analytical models of the CO rotation curves and the estimated uncertainties on the model parameters, we extract at each radius the circular velocity, V circ , the corresponding angular velocity, Ω circ , the logarithmic derivative of the rotation curve and Oort's A parameter These parameters (and their associated uncertainties) will be used to describe the local galactic dynamical properties at various locations within the target galaxies.

Morphological Environment Masks
We use the environment masks presented in Querejeta et al. (2021) to distinguish different morphological regions in each galaxy. These masks are constructed based on structural decomposition analysis and visual inspection of the IRAC 3.6 µm data (also see Herrera-Endoqui et al. 2015;Salo et al. 2015). The full set of environment masks mark the area covered by morphological features such as galaxy centers, stellar bars, spiral arms, rings, and lenses (see the last panel in Figure 2). The typical width of these environmental masks are set by the physical extent of the corresponding morphological features, which are often 1 kpc wide for stellar bars and spiral arms but can be much narrower for the other features.
In this work, we primarily use these masks to divide each galaxy into two types of environment: the area that falls into galaxy centers and stellar bars (referred to as "center/bar" hereafter), and the remaining outer disk area ("disk" hereafter). We make this distinction because we expect the physical conditions influencing GMCs to be different between these two regimes: the "center/bar" environment often sees galactic dynamics (i.e., gravitational torque and shear) playing a more prominent role, and in some galaxies AGN feedback can significantly impact the molecular gas in its central region.

Other Data
In addition to what has been described above, we also include measurements derived from other data sets in the analysis. These measurements are not presented among the main scientific results in this paper, but they are part of our final data products and they have appeared in publications that used our data products (e.g., Querejeta et al. 2021).
We use continuum-subtracted, narrow-band Hα imaging data to provide alternative estimates of star formation rate in 60 out of our 80 targets. These observations were obtained as part of the PHANGS-Hα survey 5 by either the Wide Field Imager (WFI) on the ESO/MPG 2.2-m Telescope or the Direct CCD on the CIS 2.5-m Irénée du Pont Telescope (A. Razza et al. in preparation). The narrow-band Hα data have been calibrated astrometrically and photometrically, corrected for sky emission, and masked for foreground stars; the continuum contribution was removed based on the associated R-band observations; and the continuum-subtracted data were further corrected for filter transmission and [N ii] contamination.
We combine Hα data with WISE 22 µm data to derive an attenuation-corrected SFR surface density following Calzetti et al. (2007) and Murphy et al. (2011): Σ SFR M yr −1 kpc −2 = 2.7 × 10 13 I Hα erg s −1 cm −2 arcsec −2 + 2.7 × 10 −3 I 22 µm MJy sr −1 cos i , This prescription assumes constant star formation over 100 Myr and a Kroupa IMF (Kroupa 2001 3. CROSS-SPATIAL-SCALE ANALYSIS We adopt a "cross-spatial-scale" analysis framework to connect molecular cloud properties (measured on 60−150 pc scales) to galactic environmental properties (mostly measured on ∼kpc scales). This analysis framework is inspired by a number of previous works (e.g., Sandstrom et al. 2013;Leroy et al. 2016). Briefly, we divide the sky footprint of each galaxy into a set of averaging apertures, within which we aggregate high resolution molecular gas measurements to characterize the underlying cloud population. We also attempt to build a full inventory of ancillary measurements to characterize various aspects of the host galaxy itself. In this way, we assemble the diverse set of observational data described in Section 2 into a coherent, multiwavelength database. An early version of this database was constructed by Sun et al. (2020b), with its subsequently improved versions used in several publications (e.g., Herrera et al. 2020;Jeffreson et al. 2020;Barnes et al. 2021;Querejeta et al. 2021;Stuber et al. 2021). The source code for database construction, including generic tools for aggregating measurements from maps and catalogs into the existing database, is available on GitHub 6 , and a copy of the version used in this article is published on Zenodo (Sun 2022).

Defining Averaging Apertures
We divide the sky footprint of each galaxy into a set of hexagonal apertures, as illustrated in Figure 2. These apertures form a regular tiling in the plane of the sky, with a "central" aperture positioned right at the galaxy center. Adjacent apertures are separated by a linear distance of 1.5 kpc, which implies that each aperture has a projected area of 1.95 kpc 2 on the sky. The configuration of the hexagonal apertures here is analogous to the "solution pixels" used in Sandstrom et al. (2013), except that the apertures in the current work do not overlap with each other.
For a complete coverage of the galaxy footprint, we include all apertures covering out to r gal = 1.5r 25 in each galaxy, where r 25 is the galaxy radius defined by its 25 mag/arcsec 2 isophote (in B band; see Table A1). This way, the constructed database for each galaxy includes almost all valid measurements from all data sets described in Section 2. However, this work focuses on the correlation of molecular clouds and their galactic environments, and thus we will only present results from a subset of apertures that enclose non-zero signals in the PHANGS-ALMA CO moment maps.
In addition to the hexagonal tiling method described above, we also run a parallel line of analysis with a different binning scheme. Specifically, we define a series of radial bins that are 500 pc in width and again cover out Figure 2. This figure showcases a subset of the multiwavelength data that we assemble for the galaxy NGC 628. The top row displays the PHANGS-ALMA CO (2-1) line intensity map (tracing molecular gas), the THINGS VLA H i 21-cm line intensity map (tracing neutral atomic gas), and the S 4 G Spitzer IRAC 3.6 µm image (tracing stellar mass). The bottom row displays the GALEX 154 nm image (tracing unobscured star formation), the WISE 22 µm image (tracing obscured star formation), and the PHANGS environment mask. In each panel, the scale bar at the lower-right corner shows the spatial extent of 1 kpc, whereas the white ellipse at the lower-left gives the beam size (except for the last panel). The white grids demarcate the hexagonal apertures (1.5 kpc in size) in which we extract molecular cloud population statistics and build a comprehensive inventory of host galaxy structural, kinematic, and star formation properties.
to r gal = 1.5 r 25 in each galaxy. Assembling measurements in these radial bins allows us to rigorously calculate their radial profiles, but at the expense of losing all non-axisymmetric information. We publish the data products from these radial profile calculations together with those from the hexagonal aperture analysis (see Appendix F). We do not present the results of this parallel line of analysis in this paper, but we expect the distributions of most measurements to be consistent with the hexagonal aperture averages once we use consistent weighting schemes (e.g., by the enclosed area or molecular gas mass; see Section 4) for each aperture/ring.

Aggregating Molecular Cloud Measurements
Within each aperture, we calculate the ensemble average of molecular cloud measurements using a molecular gas mass-weighted averaging scheme. This is equiva-lent to a CO intensity-weighted averaging, because the α CO value is calculated per aperture rather than per object/pixel in this study (see Appendix B). We use a " " symbol to denote this averaging operation: Here, X i, θ pc represents a molecular gas property measured for the i-th object or pixel at θ pc resolution (θ = 60, 90, 120, 150); it can be any of the objector pixel-based measurements defined in Sections 2.1.1 and 2.1.2. M i is the molecular gas mass associated with the object or pixel for which X i, θ pc is measured. The summation in Equation 23 includes all detected objects/ pixels with their center coordinates located inside the sharp boundary of the averaging aperture. In this case, each object/pixel belongs to a unique averaging aperture, and thus the averaging results in adjacent apertures are independent by construction. Based on Equation 23, we can also estimate statistical uncertainties for the population-averaged cloud properties through Gaussian error propagation. We take into account the uncertainties on both the quantity to be averaged, X i, θ pc , and the weight, M i . When aggregating the pixel-by-pixel measurements, we further consider the built-in correlation between adjacent pixels and scale the estimated uncertainty of the population average according to the oversampling factor.
Our aperture averaging scheme resembles the one adopted by Leroy et al. (2016) but differs from that approach in important ways. In that work, the averaging is performed via a Gaussian kernel convolution, in which case the averaging result at any given location has a nonzero response to molecular clouds far away from that location. This response pattern is designed to replicate the Gaussian beam of low-resolution data sets, and thus it may be preferable for rigorous calculations combining cloud-scale and kpc-scale measurements (e.g., Leroy et al. 2017;Utomo et al. 2018, see also L. Neumann in preparation). However, such an extended response pattern can lead to built-in correlations between averaging results at adjacent locations. More importantly, it can yield biased population statistics when, for example, studying a region with little molecular gas next to a very gas-rich region (such as a galaxy center). Since one of the main goals of this work is to derive reliable molecular cloud population statistics, we deem the "sharp boundary" scheme more appropriate here and will use it consistently for calculating both cloud population statistics and host galaxy properties (also see Section 3.3).

Molecular Gas Clumping Factor
The averaging operation described above essentially extracts the (mass-weighted) expectation value of a molecular gas property from its probability distribution within each averaging aperture. But one can also extract other types of statistics from the same distribution, such as the standard deviation of the gas surface density distribution (which quantifies the inhomogeneity of the medium), or the slope of the GMC mass function. These other types of statistics can also provide unique observational constraints on the physical processes driving molecular cloud formation and evolution.
As part of the analysis done for this work, we calculate the molecular gas "clumping factor," which is a dimensionless characterization of the surface density inhomogeneity in each aperture : Here Σ i, θ pc is the molecular gas surface density measured in the i-th pixel at θ pc resolution. Similar to Equation 23, the summation includes all pixels with CO detections within the averaging aperture, and N pix is the total number of such pixels. The right hand side of Equation 24 can be interpreted as the ratio between the mass-weighted mean and the area-weighted mean of molecular gas surface density in the limit of infinite sensitivity (see Leroy et al. 2013). We note that c pix is a measure of the width (i.e., second moment) of the surface density distribution among many similar parametrizations in the literature (e.g., the smoothness index and the Gini coefficient; see Davis et al. 2022, and references therein). To measure this type of parameter reliably, a careful treatment of nondetections is particularly important. We describe our strategy to handle non-detections in Section 3.2.2, and we illustrate the amplitude of the necessary corrections in Section 4.2 and Appendix D.

CO Flux Completeness and Corrections
The ensemble-average molecular cloud properties (Equation 23) and the molecular gas clumping factor (Equation 24) are both calculated based on pixels/ objects that are detected in the PHANGS-ALMA CO data. For these calculations to reflect the true statistics of the entire molecular cloud population in each region, the CO detections need to be reasonably complete, such that they represent a significant portion of the underlying cloud population.
Our object-and pixel-based measurements come from the PHANGS-ALMA CPROPS catalogs and the "strict" moment maps (Sections 2.1.1 and 2.1.2), thus the completeness of our analysis is determined by the completeness of these data products. Both data products adopt similar signal identification criteria to extract high-confidence CO detections in the original data cubes (Leroy et al. 2021b;Rosolowsky et al. 2021), which ensure reliable CO line measurements for the detected pixels/objects. However, this comes at the price of excluding faint CO emission, which renders these data products incomplete in terms of both flux coverage and area coverage.
The extent of this effect can be quantified by the CO flux completeness, f flux , and area coverage fraction, f area , of the CPROPS catalog or the strict moment maps (see tables 15 and 16 in Leroy et al. 2021a). Here, we calculate f flux and f area for each averaging aperture and report these values along with the ensemble-average molecular cloud properties. Specifically, within the footprint of each aperture, we calculate f area by comparing the total area covered by CO detections in the "strict" moment-0 map to the total area of the aperture. We calculate f flux by comparing the total CO flux included in the "strict" moment-0 map to that in the corresponding "broad" moment-0 map. The latter map is constructed with much more inclusive signal identification criteria than the strict map and has nearly 100% flux completeness (for more details, see Leroy et al. 2021b).
The incomplete CO flux and area coverage of the "strict" moment maps and CPROPS catalogs introduces a selection bias in our analysis. The sense of this bias is that we miss places where CO emission is too faint to meet the masking criteria (e.g., areas occupied by small, low-mass molecular clouds or a diffuse gas component). This selection bias affects many of the ensemble-average cloud properties calculated in this study, and is particularly severe for the clumping factor (see Section 4.2 and Appendix D).
To account for this systematic bias, we introduce a correction factor for our measurements in each aperture based on the f flux and f area values in that aperture (see Appendix D for detailed derivations). We assume that the CO intensity distribution (or equivalently, molecular gas surface density distribution) has a lognormal shape within each averaging aperture, and that the aforementioned selection bias prevents us from detecting CO emission below an intensity threshold. Under these assumptions, we can solve for the width of the lognormal intensity distribution as well as its centroid (relative to the intensity threshold) from f flux and f area . This in turn allows us to calculate the appropriate correction factors to apply to the ensemble-average molecular cloud surface density, Σ obj and Σ pix , and the clumping factor, c pix . We then assume that the correction factor calculated for Σ obj and Σ pix also applies to the ensemble-average cloud mass and turbulent pressure, but we leave molecular cloud size, velocity dispersion, and virial parameter uncorrected. Though these latter quantities likely do suffer sensitivity-induced selection biases (e.g., see the illustration of selection functions in Sun et al. 2018), the appropriate functional forms of their completeness corrections remain uncertain at present. Finally, we scale the corresponding (statistical) uncertainty for each ensemble-average value by the same correction factor.
As illustrated in Section 4.2 and Appendix D, thanks to the relatively high flux completeness of the PHANGS-ALMA CO data, the correction factors on the average cloud surface densities and the clumping factor are both moderate (<0.3 dex for 90% of the apertures with CO detections). Nevertheless, we do expect our completeness correction scheme to be less reliable for apertures with low f flux and/or f area , in which case the extrapolation is done based on very few measurements. For this reason, we will exclude apertures with low f flux or f area when performing analyses that requires accurate cloud population statistics in individual apertures (see Section 5 and Appendix D).

Aggregating Local Environmental Metrics
In addition to the compilation of ensemble-average molecular cloud properties described above, we assemble an inventory of "environmental metrics" that delineate various host galaxy local properties within each averaging aperture. This inventory covers orbital kinematic properties (derived from rotation curves), gas-phase metallicity (predicted from scaling relations), surface densities of molecular gas, atomic gas, stellar mass, and SFR (estimated from multiwavelength imaging data), and morphological environment information (inherited from environmental masks).
We generally use two schemes to integrate these environmental metrics into the databases of aperture-wide statistics. For those metrics that are calculated analytically (e.g., galactocentric coordinates, metallicity) or interpolated from analytical models (e.g., rotation curverelated properties), we directly record their values at the location of the aperture center. For those metrics that rely on two-dimensional images, we use the native resolution images and calculate the unweighted average among all pixels inside the sharp boundary of each aperture. This latter scheme is consistent with the averaging scheme we used for aggregating molecular cloud properties (modulo the different weighting), and thus allows for direct comparisons between the averaging results.
We elaborate the specific treatment for each type of environmental metric below: • Coordinates. For each hexagonal aperture, we record its central R.A. and Dec. coordinates. Then based on the center coordinates, inclination angle, position angle, and the distance of the galaxy (see Table A1), we calculate the deprojected galactocentric radius, r gal , (in kpc units) at the aperture center and the deprojected azimuthal angle, φ gal , in the galaxy plane with respect to the major axis direction. These coordinates uniquely determine the location of each aperture both on the sky and in the deprojected galaxy plane.
• Orbital kinematics. We report local orbital kinematic properties for apertures in the galaxy sample and galactocentric radius range covered by the rotation curve measurements from Lang et al. (2020). As detailed in Section 2.5, these orbital properties include the circular velocity, V circ , angular velocity, Ω circ , logarithmic derivative of the rotation curve, β, and Oort's A parameter, A Oort . They are calculated by interpolating the functional fitting model of the rotation curves at the location of the aperture center.
• Metallicity. We report the predicted gas-phase metallicity in each aperture using a prescription similar to the one described in Sun et al. (2020b), but with a few methodological improvements. In short, we first infer the metallicity at r gal = 1.0 r e in each galaxy based on a galaxy global mass-metallicity relationship (Sánchez et al. 2019), and then extrapolate to all r gal assuming a fixed radial metallicity gradient of −0.1 dex/r e within each galaxy (Sánchez et al. 2014).
For better methodological consistency with the original references, here we approximate the galaxy ef-fective radius as r e ≈ 1.68 r disk , where r disk is the stellar disk scale length. We also elevate the global stellar masses in Table A1 by 0.1 dex before substituting their values into the mass-metallicity relationship. We refer interested readers to Appendix B for more details about these adjustments.
• Molecular gas surface density (kpc-scale). We report the area-weighted mean molecular gas surface density, Σ mol , in each kpc-scale aperture. We emphasize the distinction between this measurement and the massweighted average of molecular cloud surface density, Σ [pix|obj] , defined in Section 3.2. The area-weighted mean Σ mol here is calculated from the total CO flux inside the hexagonal boundary of each kpc-sized aperture divided by its total deprojected area. For this particular calculation, we use the native resolution "broad" moment-0 map to ensure a high flux completeness (see Section 3.2.2). We then use the same metallicity-dependent CO-to-H 2 conversion factor to convert CO line intensity into mass surface density unit, as we do in Equation 8. We note that our methodology for calculating this kpc-scale aperture averaged Σ mol is different from the one used in Sun et al. (2020b). There, the kpc-scale Σ mol was derived via convolving the CO moment-0 maps to a fixed 1 kpc resolution and then sampling the convolved maps at the aperture centers. As discussed above, the new averaging scheme in this paper leads to better methodological consistency with our calculation of molecular cloud population statistics.
• Atomic gas surface density. We report the areaweighted mean atomic gas surface density, Σ atom , in all apertures for which we have H i data (see Table A1). This is calculated in the same way as the area-weighted mean Σ mol : we divide the total H i 21 cm line flux inside the hexagonal aperture by the aperture area, and then convert it to mass surface density unit via Equation 14.
Since the H i data resolution is typically comparable to or coarser than our adopted aperture size, our calculated Σ atom might not reflect the true atomic gas surface density inside the sharp aperture boundaries, but rather a slightly "smoothed" version of it. However, the atomic gas distribution is usually much smoother than the molecular gas (e.g., see Leroy et al. 2013), and Σ atom only plays a minor role throughout this paper. The resolution degradation is thus not a serious concern for the following analysis.
• Stellar mass surface density. We report the areaweighted mean stellar mass surface density, Σ , in each aperture. We calculate Σ via Equations 15 or 16 based on the mean WISE 3.4 µm or IRAC 3.6 µm surface brightness at 7.5 resolution within sharp aperture boundaries. We determine the stellar M/L ratio, Υ 3.4 µm , for each aperture by sampling the M/L ratio maps from Leroy et al. (2021a) at the location of the aperture center.
• SFR surface density. We report the area-weighted mean SFR surface density, Σ SFR , in each aperture. This is primarily calculated via Equations 17-19 based on the best available UV/IR data combination (see Table A1) at 15 resolution. We note that this resolution could approach the averaging aperture size in the more distant targets in our sample, in which case concerns about correlated measurements could again arise. To evaluate these concerns, we compare the UV/IR-based Σ SFR measurements with Hα-based measurements (the latter includes data at much higher angular resolution; see Section 2.7). We find quantitatively consistent results at Σ SFR 10 −3 M yr −1 kpc −2 , which is the range of interest in this paper (see Figure 3 below).
• Morphological environment. We keep track of the morphological regions each averaging aperture inhabits in the host galaxy (see Section 2.6). Because of the kpc-scale sizes of these apertures, some of them could stretch across multiple morphological regions.
To deal with this ambiguity, we calculate the fraction of CO flux originating from each morphological region (especially galaxy centers and stellar bars) relative to the sum over the entire aperture. We then classify all apertures that have a non-zero 7 CO flux contribution from galaxy centers or stellar bars as "center/bar" apertures, and all the remainder as "disk" apertures.

Outcome of the Cross-spatial-scale Analysis
Our analysis yields a rich value-added database for each of the 80 galaxies listed in Table A1. These databases present the molecular cloud populations residing in each galaxy, along with the large-scale gas and stellar mass distribution, kinematic information, morphological structures, and star formation activities of the galaxy disk itself. Together, these high-level measurements have a broad range of applications (see Section 7). They are published in the form of machinereadable tables online (see Appendix F).
Our databases include 46,628 apertures in total. These apertures collectively cover the footprint of every target galaxy out to a galactocentric radius limit of 1.5 r 25 . The majority of these apertures have local environmental measurements derived from multiwavelength data (such as UV, and IR), yet only a smaller subset of them have valid molecular gas measurements from PHANGS-ALMA CO data. This is because the footprint of the PHANGS-ALMA survey is often more confined and covers only the inner, molecular gas-rich part of the galaxy disk (see Leroy et al. 2021a). Since this paper focuses on linking the molecular cloud population to their local environment, in the following sections, we restrict ourselves to a subset of 3,383 apertures that are inside the PHANGS-ALMA survey footprint and show detectable CO emission in the 60-150 pc scale "strict" moment maps. Nonetheless, the full set of 46,628 apertures will be included in the public data release given the rich information provided by the multiwavelength ancillary data alone.

DISTRIBUTIONS OF AVERAGE MOLECULAR CLOUD PROPERTIES AND SUB-GALACTIC ENVIRONMENTS IN PHANGS-ALMA
In this section, we characterize the distributions of region-averaged molecular cloud properties and host galaxy local properties across the full PHANGS-ALMA data set. To do this, we use the databases constructed in Sections 3 and focus on 3,383 apertures with CO measurements from PHANGS-ALMA (including 2,724 apertures classified as "disk" and the remainder as "center/bar"). In the main text, we will only present the statistics of molecular cloud measurements at 150 pc scales, which is the best common resolution achievable for all galaxies. Quantitative comparisons across different resolutions are shown in Appendix E.

Sub-galactic Environments Probed by PHANGS-ALMA
Our multiwavelength measurements provide a multifaceted depiction of the range of local galactic environments probed by the PHANGS-ALMA survey. To this end, Figure 3 shows the histograms of 12 local environmental metrics across 3,383 apertures. We also calculate statistics such as the median value and 16−84 percentile range for each environmental metric and tabulate them in Table 1. These statistics are calculated from the histogram using two different weighting schemes: simple counting of the number of apertures or weighting each aperture by the molecular gas mass it encloses. The first scheme treats all apertures equally, and the calculated statistics reflect a typical kpc-sized area covered by PHANGS-ALMA; the latter scheme instead treats each unit of gas mass equally, and the calculated statistics reflect the local environment in which most molecular gas resides.
Below we split the 12 environmental metrics into four topical groups and comment on the corresponding histograms and statistics.
• Galactocentric radii. The PHANGS-ALMA CO measurements cover a wide radial range in terms of both absolute and normalized r gal (panels a and b). When weighting by the number of apertures, we find median values and ±1σ ranges of r gal = 5.4 +3.1 −2.7 kpc and r gal /r disk = 1.8 +1.0 −0.8 across all apertures. We find smaller values when weighting each aperture by its encircled molecular gas mass. This reflects that the molecular gas distribution typically peaks toward the galaxy center, and thus apertures at smaller radii often enclose more molecular gas mass.
We note that the r gal histogram appears "quantized" simply due to the fixed 1.5 kpc linear size of the hexagonal apertures and their tiling pattern on the sky. This behavior is not obvious in the r gal /r disk histogram because the normalization factor r disk varies among galaxies, which effectively "smooth" the histogram.
• Kinematic properties. For the subset of apertures located in the 62 galaxies with CO kinematic measurements, we report the distributions of orbital angular velocity (panel c) and Oort's A parameter (panel d ). Weighting all apertures equally, we find typical Ω circ = 32 +16

−8
Gyr −1 and A Oort = 14 +6 −4 km s −1 kpc −1 across galaxy disks, which translate to an orbital period of ∼ 200 Myr and a local shearing timescale of ∼ 70 Myr (i.e., the reciprocal of A Oort ; also see Section 6). These values suggest that the kinematic properties of a typical kpc-sized area probed by PHANGS-ALMA are very similar to those of the Solar Neighborhood (Ω circ = 27.8 Gyr −1 and A Oort = 15.3 km s −1 kpc −1 ; Bovy 2017a).
We also find that apertures located in galaxy centers and stellar bars show systematically higher Ω circ and A Oort values. This is expected from their locations at smaller r gal and the stronger shear often observed in these environments.
• Galaxy disk mass components. Weighting all apertures in galaxy disks equally, we find typical surface densities of Σ = 65 +80 −34 M pc −2 and Σ mol = 5.4 +10. 8 −3.7 M pc −2 (panels e and f ). Among the 53 galaxies with H i 21 cm line data (see Table A1), we find a typical total gas surface density of Σ gas = Σ mol + Σ atom = 13 +14 −7 M pc −2 (panel g). This gives a typical gas fraction of f gas = Σ gas /(Σ + Σ gas ) = 0.16 +0.14 −0.07 (panel i ) and a molecular fraction of f mol = Σ mol /Σ gas = 0.53 +0.22 Examining the corresponding molecular gas massweighted statistics for galaxy disks, we find that most molecular gas mass resides in environments with even higher surface densities (Σ = 110 +140 −60 M pc −2 , Σ mol = 17 +29 −11 M pc −2 , Σ gas = 28 +30 −15 M pc −2 ) and molecular fraction (f mol = 0.71 +0.13 −0.20 ). For comparison, these gas surface densities are likely higher than the averaged value across any kpc-sized neighborhood in our Galaxy (e.g., Nakanishi & Sofue 2006;Spilker et al. 2021, though the central 1 kpc might be an ex- Figure 3. The full range of host galaxy local properties sampled by PHANGS-ALMA, outlined here by stacked histograms of 12 "local environmental metrics" across 3,383 apertures. The panels show: (a) galactocentric radius, (b) galactocentric radius normalized by the disk scale length, (c) orbital angular velocity, (d) Oort's A parameter, (e) stellar surface density, (f) molecular gas surface density, (g) total gas surface density, (h) SFR surface density, (i) gas fraction, (j) molecular fraction of the gas, (k) molecular gas depletion time, and (l) total gas depletion time. We use blue and orange colors to distinguish the contributions from the "disk" and the "center/bar" subsamples in the summed histogram (black solid outline). The symbols and horizontal error bars at the top show the median value and ±1σ range (i.e., 16−84 percentile range) within each subsample. The two symbol types correspond to two different weighting schemes for calculating the median values and percentiles: weighting by number of apertures (open squares) versus weighting by molecular gas mass (solid circles). ception given uncertainties in the conversion factor there).
• Star formation activity. The typical range of SFR surface density of "disk" apertures, when weighted by simple number counts, is Σ SFR = 3.5 +6.8 This is again comparable to the estimated Solar Neighborhood SFR surface density at the present day (Σ SFR = 1.7 × 10 −3 M yr −1 kpc −2 ; Bovy 2017b). Combined with the measured Σ mol and Σ gas in our sample, this implies typical depletion times of t dep, mol = Σ mol /Σ SFR = 1.5 +1.0 −0.7 Gyr for the molecular gas (panel k ) and t dep, gas = Σ gas /Σ SFR = 3.3 +2.2 −1.3 Gyr for the total gas (panel l ).
In comparison, the molecular gas mass-weighted statistics reveal that most of the molecular gas mass resides in more actively star-forming environments with Σ SFR = 9.6 +20.9 −6.4 × 10 −3 M yr −1 kpc −2 . Yet associated gas depletion times appear similar to the aperture number-weighted values, with t dep, mol = 1.8 +1.0 −0.7 Gyr and t dep, gas = 2.8 +1.7 −1.1 Gyr. In other words, the SFR surface density is proportionally higher in these environments as their gas surface densities are.
In summary, the PHANGS-ALMA survey covers an wide variety of host galaxy local environments. When weighting all apertures equally, the most representative local environment in our sample closely resembles the Solar Neighborhood in many aspects. Comparatively, most molecular gas mass is hosted in regions that are closer to the galaxy center, have higher surface densities of stars, gas, and SFR, and are possibly not matched by any kpc-scale regions in our Galaxy.

Molecular Cloud Populations in PHANGS-ALMA
Our calculations aggregate individual molecular cloud measurements to yield mass-weighted average properties for each aperture (Section 3.2). The distributions of these population-averaged measurements offer a comprehensive portrait of how cloud populations vary across PHANGS-ALMA. This is demonstrated by Figure 4, which shows histograms of population-averaged cloud properties measured from object/pixel-based approaches (originally measured at 150 pc). These histograms include all 3,383 apertures with pixel-based data and 2,784 apertures with object-based data. The latter number is smaller because CPROPS uses a slightly more stringent criterion for identifying objects in CO data cubes (objects made of too few cube pixels are rejected even when they satisfy the S/N criteria; Rosolowsky et al. 2021).
Similar to Section 4.1, here we calculate the median values and 16−84 percentile ranges with two different weighting schemes (Table 2). We note that because most variables depicted in Figure 4 (except c pix ) are already mass-weighted averages within individual apertures, an additional mass-weighted averaging step across all apertures can be interpreted as the mass-weighted average cloud properties combining all clouds in all galaxies in PHANGS-ALMA. In comparison, the other weighting scheme (i.e., same weight for all apertures) gives us a view of the typical molecular cloud population likely to be found at a random location in a PHANGS-ALMA galaxy. As results from the latter weighting method is less straightforward to interpret, below we focus mostly on the mass-weighted statistics in our discussion.
• Molecular cloud mass and size. The mass-weighted average molecular cloud mass at 150 pc resolution is 6.9 +10.5 −4.5 × 10 6 M for all clouds in galaxy disks (panel a). This value is high compared to the typical mass of molecular clouds in the Milky Way (e.g., Rice et al. 2016;Miville-Deschênes et al. 2017;Colombo et al. 2019), but is consistent with numbers measured in nearby galaxy studies (e.g., Hughes et al. 2013a). As pointed out by Rosolowsky et al. (2021), the finite resolution and sensitivity of the PHANGS-ALMA CO data limit our ability to identify molecular clouds with mass 10 5 M . Specifically, the 60-150 pc resolution of the PHANGS-ALMA data would lead to individual, moderate size clouds being blended into a single object by CPROPS (also see discussions on the resolution-dependence of average cloud mass in Appendix E). That said, our completeness correction can partly remedy sensitivity-related biases by compensating for isolated, less massive clouds undetected in the high resolution CO observations (the histograms for the corrected measurements extend to lower values than those for the uncorrected measurements).
In line with this consideration, the mass-weighted cloud effective radius measured at the same resolution spans 90 +19 −15 pc, which slightly exceeds half the beam FWHM size (panel b). This is consistent with a series of previous studies, all of which found that cloud segmentation algorithms tend to identify objects comparable to or slightly larger than the beam size (e.g., Verschuur 1993;Pineda et al. 2009;Hughes et al. 2013a;Leroy et al. 2016).
• Molecular cloud surface density. At 150 pc resolution, the mass-weighted molecular cloud surface density in galaxy disks is 78 +124 −50 M pc −2 from the object-based approach and 33 +60 −22 M pc −2 from the pixel-based approach (both weighted by gas mass; see panels d and h). These values are on the low end of the surface density distribution of Galactic molecular clouds (e.g., Colombo et al. 2019). This likely reflects the coarser spatial resolution of our data compared to most Galactic studies, which means our measurements would be "diluted" by low column density sightlines within each beam.
The quantitative differences between the objectand pixel-based approaches reflect that the they attempt to measure fundamentally different quantities (see Section 2.1.3). In particular, the CPROPS algorithm attempts to measure the true surface density of an identified object. Therefore, it includes additional de-convolution and extrapolation procedures, which lead to smaller cloud sizes and larger masses, and thus larger surface densities. The pixel-by-pixel analysis instead measures the surface density pointby-point from a contiguous molecular gas distribution. Without a priori expectation for the gas spatial distribution, it does not perform any de-convolution but simply extracts measurements at the resolution of the observations. Given these differences, for marginally resolved clouds the pixel-based approach would simply yield the native, beam-averaged value at the data resolution, whereas the object-based approach would yield higher surface densities as a result of the deconvolution.
Compared to the results for disk apertures, the cloud populations in galaxy centers and stellar bars have much higher mass-weighted mean surface densities of 210 +560 −160 M pc −2 (object-based) or 200 +840 −160 M pc −2 (pixel-based). Such a trend has been highlighted in previous works on the same galaxies (Sun et al. , 2020aRosolowsky et al. 2021) Figure 4. The "demographic profile" of molecular cloud populations captured in PHANGS-ALMA, illustrated here by stacked histograms of 11 population-averaged molecular cloud measurements derived from either object-or pixel-based approaches. As described in Sections 2.1 and 3, these quantities represent aperture-wise mass-weighted averages and their derivation accounts for the effect of galaxy inclination and finite data sensitivity. The panels show: (a) object molecular gas mass, (b) object radius, (c) pixel-wise molecular gas clumping factor, (e-g) object-based molecular gas surface density, velocity dispersion, turbulent pressure, and virial parameter, (h-k) pixel-based molecular gas surface density, velocity dispersion, turbulent pressure, and virial parameter. The derivation of these properties accounts for the effect of galaxy inclination (Section 2.1.1-2.1.2). We have also applied completeness corrections on a subset of these quantities to offset sensitivity-related biases (Section 3.2.2). The histograms of uncorrected measurements (gray, unfilled) are shown in contrast to those of the corrected measurements (black). and is also consistent with observations in our Galaxy (Oka et al. 2001). We do caution that these results are more sensitive to the choice of α CO prescriptions (see Appendix B). Several lines of evidence suggest lower α CO in galaxy centers (e.g., see Bolatto et al. 2013;Sandstrom et al. 2013;Israel 2020;Teng et al. 2021) and our fiducial prescription only mildly depresses α CO near galaxy centers.
The quantitative discrepancies between the two approaches here can also be explained by methodological differences. As described above, the objects identified by CPROPS are often slightly larger than the beam size. We would then expect larger CO line width measurements from the object-based approach due to both the size-line width relation in molecular clouds (e.g., Larson 1981) and additional contributions from galaxy rotation and large-scale gas streaming motions. This explains the sense of deviation for the measurements in galaxy disks. However, in places where molecular clouds are spatially crowded (such as galaxy centers), the ppv space segmentation in CPROPS helps to demarcate clouds that fall along the same line of sight but are separable in velocity space, whereas the pixel-based analysis simply measures the line effective width and thus cannot tell them apart (see e.g., Henshaw et al. 2016). This explains why the pixel-based approach yields higher velocity dispersions with a wider spread in these environments.
The sense of deviation between the two approaches here is similar to that of the surface density and velocity dispersion measurements. It is most apparent near the low pressure end, where the distribution almost always exceeds 10 4 K cm −3 in the object-based statistics but extends to below 10 3 K cm −3 in the pixel-based statistics. This aligns with the intuition that the object-based calculations focus on the denser inner portion of molecular clouds (i.e., the half of gas within the FWHM of each object), whereas the pixelbased analysis treats every chunk of molecular gas equally at fixed resolution, and thus can reflect the behavior of the lower pressure, more diffuse gas when it dominates the mass budget. (object-based; panel g) or 1.6 +0.9 −0.6 (pixel-based; panel k ). If taken at face value, these values would suggest that the dynamical state of molecular clouds are somewhere between virial equilibrium (α vir = 1) and energy equipartition (α vir = 2; see also Sun et al. 2018Sun et al. , 2020aRosolowsky et al. 2021). However, systematic uncertainties related to the sub-resolution gas distribution and the CO-to-H 2 conversion factor are especially concerning for the α vir measurements given the narrow dynamic range. More, it can be challenging to determine the true dynamical state of the ob-served gas structures from a measured α vir alone (e.g., Ballesteros-Paredes et al. 2011;Ibáñez-Mejía et al. 2016;Lu et al. 2020). The most conservative conclusions from these data are that molecular cloud populations show a relatively narrow range of dynamical states and appear near energy equipartition across a wide range of environments.
Comparing the mass-weighted average α vir between cloud populations in galaxy disks versus those in galaxy centers and stellar bars, the pixel-based statistics indicate higher α vir values for the latter, while the object-based statistics show essentially no difference. This can be explained by the same line-of-sight blending effect that drive the pixel-based velocity dispersion measurements higher in centers and bars (also see Henshaw et al. 2016;Kruijssen et al. 2019b). We suggest to prefer the object-based results in this case, but again caution that we may be overestimating α CO in these regions. If this is the case, the correct objectbased values would also suggest higher α vir in bars and galaxy centers.
We find a completeness-corrected clumping factor of 1.9 +1.1 −0.4 (uniform weight per aperture) or 1.9 +1.0 −0.4 (weighted by molecular gas mass in each aperture) across our sample (panel c). These values are markedly smaller than those reported in Leroy et al. (2013, median value of ∼7 at 20-300 pc resolution), which were calculated from early CO observations targeting a handful of very nearby galaxies (including several CO-poor Local Group members). This is partly due to the much higher sensitivity of the PHANGS-ALMA data set, and partly due to the improvement in the treatment of non-detections. The small clumping factors we derive suggest that the molecular gas is much less clumpy than reported in previous studies on a smaller set of galaxies.
In summary, we observe substantial variations in the molecular cloud population-averaged properties across all apertures in our sample. The mass-weighted average of all clouds in PHANGS-ALMA yields high masses, large sizes, and low surface densities compared to the cloud population in our Galaxy. The contrast of cloud populations in galaxy disks and those in center/bar environments are qualitatively consistent with findings in previous extragalactic and Galactic studies. While the object-based and pixel-based results display qualitatively similar trends, there exist important quantitative discrepancies, which reflect their different measurement approaches in expected ways. Finally, after correcting for the completeness of CO detections, we find that the molecular gas in PHANGS-ALMA galaxies appears significantly less clumpy than previously determined from low sensitivity CO observations of a few very nearby galaxies. Note-The median value and ±1σ range (i.e., 16−84 percentile range) of the galactic environmental metrics shown in Figure 3. Results in columns (2) to (4) are calculated by weighting each aperture equally, whereas those in columns (5) to (7) are calculated by weighting each aperture by the molecular gas mass it contains. Note-The median value and ±1σ range (i.e., 16−84 percentile range) of the population average molecular cloud properties shown in Figure 4. Results in columns (2) to (5) are calculated by weighting each aperture equally, whereas those in columns (6) to (9) are calculated by weighting each aperture by the molecular gas mass it contains. Columns (2) and (6) correspond to measurements without applying completeness corrections (see Section 3.2.2).

ENVIRONMENTAL DEPENDENCE OF MOLECULAR CLOUD POPULATIONS IN PHANGS-ALMA
In the previous section, we see strong variations in both the molecular cloud populations and the host galaxy local environments from aperture to aperture. These variations are also known to correlate with each other -numerous studies have shown that the physi-cal properties of molecular clouds depend on their host galaxy environment in various ways (see discussions in Section 1). The rich set of measurements derived in this study allow for a systematic characterization of such cloud-environment connections across an unprecedented range of physical conditions. In this section, we first summarize the basic, pairwise correlations between the population-averaged cloud properties and the lo-cal/global host galaxy properties (Section 5.1). We then perform a variable selection method to identify a subset of host galaxy properties carrying the most predictive power, and to construct empirical predictive models for the molecular cloud properties (Section 5.2).
Both the pairwise correlation analysis and the variable selection procedure requires high quality measurements for the cloud population statistics. The latter also needs to be applied on a consistent sample of apertures that all have the relevant variables measured to good quality. To meet these requirements, this part of the analysis works with a subsample of 871 apertures from 42 galaxies. These apertures are selected because: (1) they have the most complete multiwavelength data coverage, such that none of the cloud population statistics or host galaxy environmental metrics are missing; and (2) the PHANGS-ALMA CO data have reasonably high flux completeness and area coverage fraction inside these apertures (f flux > 0.5 and f area > 0.2, see Appendix D), such that our measured cloud population statistics represent a significant portion of the molecular gas residing in that area. This down-selection primarily restricts our sample to regions with higher Σ mol (with all apertures weighted equally, the median value and 16-84 percentile range is 14 +17 −8 M pc −2 for the sub-sample, in comparison to 6.0 +13.3 −4.1 M pc −2 for the parent sample). Consequently, it also tends to select apertures in more massive, molecular gas-rich galaxies and at smaller r gal within each galaxy (4.4 +2.4 −1.8 kpc for the sub-sample versus 5.4 +3.1 −2.7 kpc for the parent sample). Nevertheless, the selected apertures still cover a considerable range of the relevant parameter space to allow for the following analyses.

Pairwise Correlation
To provide an empirical characterization of the observed cloud-environment connections, we calculate Spearman's rank correlation coefficients for each pair of population-averaged molecular cloud property and host galaxy environmental metric within the sub-selected 871 apertures. For this analysis, we consider seven local environmental metrics (r gal , Σ mol , Σ atom , Σ , Σ SFR , Ω circ , and A Oort ), each of which holds a unique piece of information unavailable to all other variables. We complement these local properties with five global galaxy properties (r disk , M mol , M atom , M , and SFR), so that our correlation analysis can also capture galaxy-to-galaxy trends. These global properties are measured by (Leroy et al. 2021a) for all PHANGS-ALMA galaxies.
The top left part of Figure 5 summarize the outcome of this pairwise analysis. Broadly speaking, we find significant correlations between most pairs of populationaveraged cloud properties and environmental metrics. The signs of the correlation coefficients indicate that an average molecular cloud tends to be denser, more massive, more turbulent, and more strongly self-gravitating at places that are closer to the galaxy center, have higher gas and stellar content, show more active star formation, and feature shorter orbital period and stronger local shear. The correlations of population-averaged cloud properties versus global galaxy properties are weaker, but most of them are statistically significant. These findings agree well with previous observations targeting individual galaxy or smaller galaxy samples (see references in Section 1).
Beyond these general trends, we highlight three interesting patterns in Figure 5. First, most molecular cloud properties show the strongest correlation with the (kpcscale) aperture-averaged molecular gas surface density, Σ mol . While some systematic effects (e.g., uncertainties in the CO-to-H 2 conversion factor or calibrations of the raw data) could influence both the independent and dependent variables here, these correlations, including the ones regarding cloud surface densities, still carry real information about multi-scale structures in the molecular gas. Specifically, the correlation strength between the cloud-scale and kpc-scale surface densities partly reflects the inhomogeneity of the molecular gas on spatial scales between 1.5 kpc (aperture size) and 60−150 pc (data resolution). The limiting case of a perfect correlation appears only if the gas distribution is completely homogeneous, or if it is structured in a way that the clumping factor is the same in all apertures (which is not far from the reality given the narrow range of clumping factors observed across our sample; see Section 4.2). On the contrary, we would expect no correlation if all molecular gas is concentrated into small, isolated clouds with fixed surface densities.
Second, for the cloud radius and virial parameter, the correlations with environmental metrics are weaker compared to the other cloud properties. This is mainly due to the narrow dynamic range of these quantities in our data. For the cloud radius, the narrow range is somewhat imposed by the limited data resolution and our adopted object identification algorithm (see Section 2.1.1 and 4.2). For the virial parameter, the narrow range across our sample is more intrinsic and reflects the relative uniformity of the cloud dynamical state (see Sun et al. 2018Sun et al. , 2020bRosolowsky et al. 2021).
Third, when comparing some of the local environmental metrics to their corresponding "integrated", galaxy global measurements (i.e., Σ mol -M mol , Σ atom -M atom , Σ -M , and Σ SFR -SFR), the correlation coefficients for the latter are always smaller. One possible explanation is that the correlations between cloud populations and their local environment are more fundamental, to the extent that all galaxy-to-galaxy trends arise as their consequences. In other words, the apparent relationships between cloud populations and global galaxy properties might be completely mediated by the local properties. This hypothesis is challenging to test based on the pairwise correlation coefficients alone, as the strengths of mutual correlations between the independent variables (a.k.a., multicollinearity) are not explicitly modeled. Top left: Spearman's rank correlation coefficients between the population-average molecular cloud properties (dependent variables) and the host galactic properties (independent variables). Darker red/blue colors indicate stronger positive/ negative correlations. The number in each entry is the corresponding correlation coefficient, with black/gray font colors indicating p-values smaller/larger than 0.001. Bottom left: Variance inflation factors (VIFs) calculated for the independent variables. Larger VIFs indicates higher multicollinearity, i.e., stronger mutual correlations among the independent variables. Right: Example scatter plots illustrating the correlations between three molecular cloud property-host galaxy property pairs.
The issue of multicollinearity is in fact a general concern that impacts more than the local-global quantity pairs identified above. Many local environmental metrics considered here are known to follow scaling relations (e.g., see reviews by Kennicutt & Evans 2012;Sánchez et al. 2021), and the same is true for the global galaxy properties (e.g., Saintonge & Catinella 2022). This issue poses challenges to determining whether there are any secondary trends on top of the cloud-environment relationships with the largest correlation coefficients. To quantify the severity of this issue, the bottom left part of Figure 5 shows the variance inflation factors (VIFs) of the same independent variable set. A larger VIF means that a larger fraction of variations in that particular independent variable can be explained by the other independent variables. The VIFs for many variables exceed commonly adopted cutoff values of 5−10 (James et al. 2013), signaling strong multicollinearity.
In the following section, we address this issue of multicollinearity with an information criterion-based variable selection method.

Variable Selection
The goal of this section is to identify a subset of environmental metrics that are most directly relevant for setting each population-averaged cloud property. We attempt to distinguish the most fundamental cloudenvironment correlations from the ones that likely arise as indirect consequences of covariance among environment metrics. Here we distinguish these underlying relations through variable selection. For each cloud property (as a target variable), we compose an empirical predictive model using a minimal set of environmental metrics (feature variables) that carry most predictive power. This approach has the advantage of removing irrelevant feature variables while optimizing prediction accuracy in the face of multicollinearity. While this approach is still limited by the precision at which we can estimate each quantity of interest, it is an effective way to collapse a high-dimensional data set into concise and highly interpretable predictive models.

Variable Selection Methodology
The basis of this analysis is a multivariable linear regression in the logarithmic space. That is, we restrict the model functional forms to simple linear combinations of logarithmic variables (including an intercept term), which are equivalent to products of power-laws of the original variables (with a normalization constant). The regression is done independently for each populationaveraged cloud property as a target variable, using all the environmental metrics in Figure 5 as available features. Although we have applied inclination corrections to the measured cloud properties, we still include cos i as an extra feature to capture residual trends with inclination when they are present. After converting our feature and target variables to their logarithms, we mediansubtract all features to further reduce correlations between the fitting variables (i.e., power-law slopes versus the normalization constant).
With this regression setup, we perform a lasso model fit (Tibshirani 1996) and use the Bayesian Information Criterion (BIC; Schwarz 1978) for model selection. This is implemented with the LassoLarsIC function in the scikit-learn Python library. In detail, for a linear predictive model with the formŷ i = β 0 + m j=1 β j x ij , the lasso regression minimizes the following objective function: Here i = 1, 2, ..., n is the index of data (index of averaging aperture in our case) and j = 1, 2, ..., m is the index of features in the model. The α parameter is a nonnegative hyper-parameter, so that the second term in Equation 25 adds a penalty for the use of any non-zero slope in the fitted model. This particular "regularization" term is the reason that the lasso as a regression method can also be used for variable selection.
The lasso regression yields a best-fit model that minimizes Equation 25 for each choice of the α parameter. To guide subsequent model selection, we calculate the BIC value for each such model via where d is the number of features with non-zero slope in that particular model, and σ 2 is the noise variance for the target variable. This expression of BIC is derived from its formal definition of d ln(n) − 2 ln(L) assuming Gaussian error, withL being the maximum of the likelihood function. It is consistent with other definitions in the literature (e.g., James et al. 2013) up to irrelevant constants.
For the noise variance σ 2 , the contribution from statistical uncertainties is generally small for our measurements (typical fractional error ∼1-10%). We thus expect several sources of systematic uncertainties on the level of at least 0.1-0.3 dex to dominate the total noise variance. These sources include (but are not limited to) the estimated α CO and the adopted R 21 values, the unknown sub-resolution gas spatial and kinematic structures, and calibrations of the ALMA data. Considering these uncertainties, we conservatively use a noise variance of σ 2 = 0.1 dex 2 (i.e., about a factor of 2) for the cloud masses, sizes, surface densities, and velocity dispersions. We use a larger variance of σ 2 = 0.25 dex 2 (i.e., about a factor of 3) for the turbulent pressures and virial parameters, as they are particularly sensitive to the assumed geometry of the gas.
To complete the model selection procedure, we compare the BIC of all candidate models to their minimum (BIC min ) and identify a subset of candidate models that satisfy ∆BIC ≡ BIC − BIC min ≤ 10 (see e.g., Kass & Raftery 1995, for justifications of this ∆BIC = 10 threshold). Among this subset, we select the one model that corresponds to the largest α value, which typically includes the fewest features. This last step allows us to further suppress any less relevant features, as there is less strong evidence for their inclusion in the model.

Variable Selection Outcomes
With the lasso regression and the BIC-based model selection, we find a set of "preferred" power-law predictive models, whose analytical forms are tabulated in Table 3. We also report in Table 3 the model residual scatter, the model coefficient of determination (R 2 , which quantifies the model explanatory power), and the BIC difference between the preferred model and a null model with only the normalization term. Figure 6 illustrates the full path of the lasso regression up to the "preferred" model. The model residual scatter reduces as each new feature variable is added into the model.
We find that only a small number (0−4) of environmental metrics are included in the preferred model for each molecular cloud property, which means that most of the correlations in Figure 5 can be attributed to a more concise set of fundamental correlations. This result is in sharp contrast to the impression one would get from the correlations in Figure 5, which indicates ubiquitous, significant trends for nearly all local and global galaxy properties. Evidently, for most galactic properties considered in this work, their apparent correlations with molecular cloud properties are potentially explicable via their covariance with other galaxy properties, such that the predictive models need only a few variables. Once the modulating effect of those few variables are accounted for, we see no evidence that the remainder play a significant role at the current precision level of our measurements.  Note-The first column lists the power-law predictive models given by the lasso regression and BIC-based model selection (Section 5.2). The second column shows the residual scatters around these models. The third column quotes the coefficients of determination, i.e., the fraction of variation in the dependent variable that is explained by the model. The last column records the BIC difference between the selected model and a null model with only the normalization term.
Our variable selection exercise allows us to draw some interesting conclusions based on the functional form of the power-law predictive models in Table 3. First of all, the absence of global galaxy properties (except inclination angle, see below) in these models implies that their correlations with cloud properties are not fundamental: these correlations probably originate from the tighter connections between molecular clouds and their local (sub-galactic) environment. As star-forming galaxies follow various scaling relations, galaxies with larger size, mass, and SFR would include more sub-galactic regions with higher mass and SFR surface densities, which subsequently entails cloud populations with higher average masses, surface densities, velocity dispersions, and turbulent pressures. This is likely the driver of the observed systematic variations in molecular cloud properties from galaxy to galaxy (e.g., Hughes et al. 2013a;Sun et al. 2020a;Rosolowsky et al. 2021).
Furthermore, the galactocentric radius r gal does not appear in any of the predictive models either. Given that most physical properties of the host galaxy (including many not considered in this work) are strong functions of galactocentric radius, its general absence in the predictive models is rather encouraging. It suggests that those environmental metrics included in the models are doing a decent job in capturing most systematic trends: they likely make better proxies than r gal for many relevant physical quantities not considered in this work (e.g., radiation field, magnetic field, and cosmic ray strength).
The predictive models for specific molecular cloud properties also provide insights into various aspects of molecular cloud formation and evolution. Below we comment on these models individually: • Molecular cloud mass. The average molecular cloud mass shows primary dependence on the kpc-scale molecular gas surface density and galaxy inclination: together, these two quantities can explain 70% of all variations. On the one hand, we can make sense of the former dependence in light of the theoretical expectation that gravitationally unstable gas disks tend to fragment into objects at a specific mass scale (i.e., the Toomre mass). This mass scale is often linked to the local gas surface density and the disk vertical scale height H via M T ≈ πH 2 Σ mol (e.g., Murray et al. 2010). While H is not available for variable selection, our derived predictive models for molecular cloud mass do exhibit slightly sub-linear dependencies on Σ mol . This is consistent with the expected anti-correlation between Σ mol and H, as Σ mol rapidly declines with galactocentric radius, while H mildly increases with it in the inner part of nearby disk galaxies (e.g., Yim et al. 2011Yim et al. , 2014. On the other hand, a non-trivial portion of the trends with Σ mol and cos i could originate from observational and methodological limitations. The nearly constant cloud radii given by CPROPS, in combination with small clumping factors (see Section 4.2), would naturally produce strong correlations between masses of the identified clouds and the large-scale Σ mol . The additional cos i dependence also signifies stronger source blending in more inclined galaxies, despite CPROPS's attempt to deblend based on information in the velocity space (see Section 2.1.1). Figure 6. The full path of the lasso regression up to the "preferred" model for each population-averaged cloud property as a target variable. The model residual scatter decreases as each new feature variable is added into the model. The "preferred" model (large black dot) is selected by comparing the BIC of all models along the full regression path, as explained in Section 5.2.1. Its coefficient of determination (R 2 , see text label) can be calculated from the model residual scatter and the total scatter in the target variable (horizontal dotted line).
Therefore, any attempt to interpret the predictive model should also take these non-physical factors into consideration.
• Molecular cloud radius. The BIC-based variable selection favors the "null" model for the average molecular cloud radius, i.e., the one that includes only a normalization term and nothing else. In other words, there is no strong evidence that any of the environmental metrics considered here can effectively predict the (small) variations in the cloud radius. This agrees with the notion that the sizes of the CPROPS-identified objects are more influenced by algorithm-related factors (e.g., deblending criteria) and data characteristics (e.g., beam size; see Figure E1) than physical properties of the gas or the host galaxy environment (also see Hughes et al. 2013a;Rosolowsky et al. 2021).
• Molecular cloud surface density. For this quantity (measured with either object-based or pixel-based ap-proach), we find significant, secondary trends with Σ SFR on top of the prominent correlations with the large-scale Σ mol . A possible explanation is that regions with more clumpy molecular gas (i.e., larger clumping factor, higher Σ 150 pc /Σ mol ratio) are possibly more subject to gravitational instabilities and the gas there has more chance to form stars. This points at a possibility that knowing the molecular gas clumping factor on 60-150 pc scales could allow for better prediction of Σ SFR at a given Σ mol on kpc scales (i.e., improving upon the Schmidt-Kennicutt relation;Schmidt 1959;Kennicutt 1998a). It is also worth noting that the preferred models for the average cloud surface densities do not include a cos i dependent term. Though partly by construction, this still affirms the effectiveness of the inclination correction on the cloud surface densities (see Section 2.1.1 and 2.1.2 as well as Appendix C).
• Molecular gas clumping factor. We do not separately construct a power-law predictive model for the clumping factor because this quantity can be well approximated by Σ pix /Σ mol when the data sensitivity is sufficiently high (see Section 3.2.1). Instead, one can easily derive a predictive model for this quantity by dividing the model for Σ pix by Σ mol .
• Molecular cloud velocity dispersion. The preferred models for this quantity include the same environmental metrics, Σ mol and Σ SFR , as the models for the cloud surface densities. Since molecular gas surface density and velocity dispersion correlate strongly even on a cloud-to-cloud level (see discussions in Section 1), it is not surprising that the same environmental metrics turn out to be most relevant for both molecular cloud properties after population averaging. However, the preferred predictive models for the cloud velocity dispersion explain a much less fraction of its total observed scatter (32-52%) compared to those for the cloud surface density (74-85%). At least part of this is attributable to the latter fraction being exaggerated, because the cloud-scale and kpc-scale molecular gas surface densities tend to co-vary for many non-physical reasons (e.g., relying on the same conversion factor). It is also possible that the physical drivers of gas velocity dispersion variations are less well-captured by the set of environmental metrics included in this work (e.g., gas inflow rate is another possible driver of velocity dispersion variations; see Krumholz & Burkert 2010).
• Molecular cloud turbulent pressure. The functional forms of the preferred models for this cloud property are broadly consistent with the expectation from the models for the average cloud surface density and velocity dispersion. Interestingly, the preferred model for the pixel-based turbulent pressure also includes an extra term depending on the large-scale stellar mass surface density, Σ , albeit with a small power-law index (0.08). Adding this term improves the model R 2 by only a small amount (from 0.70 to 0.71; see Figure 6), but it lowers the model BIC by more than 10, which means that our data clearly favor the model with an extra dependence on Σ . This extra dependence is in line with theoretical models proposing that molecular clouds can be influenced by the external gravitational potential of the host galaxy stellar disk (e.g., Meidt et al. 2018Meidt et al. , 2020Sun et al. 2020b;Liu et al. 2021).
• Molecular cloud virial parameter. For this quantity, the preferred models for the object-based and pixelbased results show the largest deviation. The preferred model for the object-based measurement includes four galaxy properties (Σ mol , Ω circ , A Oort , and cos i), whereas that for the pixel-based only includes A Oort . This difference is probably related to the narrower dynamic range in α vir, pix, 150pc than in α vir, obj, 150pc (see Section 4.2). Moreover, the appearance of Ω circ , A Oort , and cos i in the models points at potential influence of galactic rotation on the inferred cloud dynamical state. Particularly, if the measured velocity dispersion includes contributions from differential galactic rotation (i.e., beam smearing), it would lead to positive correlations with Ω circ and A Oort because they reflect the strength of the differential motion, and with cos i because the beam smearing effect is more prominent in more inclined galaxies.
In summary, through the lasso regression and BICbased model selection we compose power-law predictive models for all population averaged molecular cloud properties. These models capture the primary cloudenvironment correlations with at most four environmental metrics as independent variables. The most commonly involved environmental metrics in these models are the large-scale molecular gas surface density, Σ mol , and SFR surface density, Σ SFR . Furthermore, the general absence of global galaxy properties in these models suggests that galaxy-to-galaxy variations in molecular cloud populations might be the mere consequences of their tighter connections with sub-galactic environments.

CHARACTERISTIC TIMESCALES OF MOLECULAR CLOUD EVOLUTION IN PHANGS-ALMA
Molecular cloud formation and evolution are influenced by a number of physical processes including turbulence driving and cascade, gravitational collapse, galactic rotation and shearing motions, cloud-cloud collisions, and gas depletion due to star formation. These processes not only operate over a vast span of spatial scales, but also feature different characteristic timescales. Estimating their timescales across diverse galactic environments allows us to demonstrate the balance (or not) between these processes and infer how star formation is regulated under distinct physical conditions (Wong 2009;Jeffreson & Kruijssen 2018;Kruijssen et al. 2019a;Chevance et al. 2020a,b;Kim et al. 2021a).
In this section, we estimate six different characteristic timescales as a use case demonstration for our rich multiwavelength measurements. We detail the definition and derivation of each timescale in Section 6.1, and compare the quantitative results in Section 6.2.

Timescale Definitions
• Free-fall time, t ff . This is the timescale for a molecular cloud to collapse in free-fall due to self-gravity, provided no pressure support to counterbalance it. We estimate this timescale from the mean volume density of molecular clouds under the assumption of spherical symmetry. The population-average free fall time,t ff , is subsequently calculated via (27) Here M mol and R cloud are the cloud mass and radius, estimated from either object-or pixel-based approaches 8 . The " " symbol denotes the same massweighted averaging scheme as defined in Section 3.1. We define this population-averaged free fall time as a mass-weighted harmonic mean so that it appropriately reflects the overall timescale for the whole cloud population 9 (also see Jeffreson & Kruijssen 2018;Utomo et al. 2018). We apply a similar completeness correction to this measurement as we did for the other population-averaged cloud properties (see Appendix D for more detail).
• Turbulence crossing time, t cr . This is the timescale for the turbulent flow to cross the span of a molecular cloud. We drive it from the cloud radius and the (onedimensional) turbulent velocity dispersion, and then calculate the mass-weighted harmonic mean as Here R cloud and σ mol are the radius and onedimensional velocity dispersion of individual molecular clouds, again derived from either object-based or pixel-based analyses (see Sections 2.1.1 to 2.1.2). Under this definition, the crossing time is related to the free-fall time and the virial parameter via t ff /t cr ≈ 0.50 α 0.5 vir . Therefore, the crossing time of a virialized molecular cloud would be roughly two times longer than its free-fall time.
• Orbital time, t orb . This is the period of the orbital revolution around the galaxy center. We derive it from the orbital angular velocity measured from the CO rotation curves (see Section 2.5): • Shearing time, t shear . This is the timescale for two objects to move closer/farther by a unit length azimuthally, given that they are on two circular orbits separated radially by the same unit length. It equals the reciprocal of Oort's A parameter (see Section 2.5): 8 We use M mol = M obj /2 for the object-based measurements to be consistent with our calculations in Section 2. 1.1. 9 This mass-weighted harmonic mean can be very convenient in the following scenario: if all clouds form stars on their corresponding free-fall timescale with the same efficiency per free-fall time ( ff ), one can easily derive the total SFR of a cloud population via ff Mtot/t ff , where Mtot is the total gas mass held by the cloud population, andt ff the population-averaged free fall time defined by the mass-weighted harmonic mean. . (30) • Cloud-cloud collision time, t coll . Most generally, this is the timescale for any particular molecular cloud to collide with another cloud (i.e., it is not the timescale for such collisions to happen within a given area). We estimate this timescale following a simplified model of shear-induced collision (Tan 2000). The key assumptions are that molecular clouds are randomly distributed in each aperture, and that cloud-cloud collision happens only when clouds catch up with other clouds on adjacent circular orbits due to orbital shear. In this scenario, we can estimate a populationaveraged collision time,t coll , via: Here v shear R cloud /t shear is the shear velocity of two orbits separated radially by R cloud , the average impact parameter among all collisions. λ mfp (2R cloud N cloud ) −1 is the mean free path of cloudcloud collisions given a linear cross section of 2R cloud and an area number density of N cloud . The extra factor of 2 on the numerator accounts for the fact that the other cloud can be located on either an inner orbit (smaller r gal ) or an outer orbit (larger r gal ) relative to the cloud in question (Tan 2000). Equation 31 makes it straightforward to estimatē t coll from measurable quantities in the object-based approach (or specifically, N cloud as area density of identified objects and R 2 cloud as mass-weighted average of object radius squared). Yet it is not trivial to measure them with the pixel-based approach in a totally symmetric way. Alternatively, we follow a similar line of reasoning as Tan (2000) to approximatē t coll from other pixel-based measurements: The second step assumes all molecular gas is concentrated into clouds with characteristic surface densities Σ cloud and radii R cloud , such that Σ mol ≈ N cloud (Σ cloud πR 2 cloud ). The third step follows from the definition of the molecular gas clumping factor, c pix , as the contrast between molecular gas surface densities on cloud scales and kpc scales (see discussion in Section 3.2. 1 and Leroy et al. 2013).
We caution that the simplifying assumptions involved in Equations (31) and (32) likely bias ourt coll estimates high. In reality, molecular clouds are not evenly distributed and have random motions in addition to circular rotation (also see Dobbs et al. 2015). Possible blending of multiple clouds in a single beam or a single identified object can also lead to longer estimated collision timescales than reality.
• Molecular gas depletion time, t dep, mol . This is the timescale to convert all molecular gas into stars at the current SFR, provided no other sources or sinks for the gas: Figure 7 shows the statistical distributions of all six timescales, including their variants derived from objector pixel-based measurements. For reference, we also mark the range of molecular cloud lifetimes, t life , measured from the spatial (de-)correlation of molecular gas and young star tracers in a subset of PHANGS-ALMA targets (Chevance et al. 2020a,b;Kim et al. 2021a). Table 4 summarizes the median values and 1σ ranges of all our estimated timescales.

Timescale Comparisons
Based on Figure 7 and Table 4, we can identify three distinct groups of timescales separated by roughly an order of magnitude apart from each other. The first group consists oft ff andt cr , both around 5−20 Myr. These timescales correspond to physical processes taking place inside molecular clouds. The median values of t ff and t cr differs by less than a factor of two, which is no more than a restatement that most molecular clouds have virial parameters of order unity. Furthermore, they appear comparable to (or only slightly shorter than) t life , as seen in previous observations and simulations 10 (Fukui et al. 1999;Elmegreen 2000;Kawamura et al. 2009;Murray 2011;Grudić et al. 2018;Kim et al. 2018;Kruijssen et al. 2019a;Benincasa et al. 2020;Chevance et al. 2020a). This is not inconsistent with our estimated α vir = 1−2, as even clouds in free-fall collapse can yield apparent virial parameters of ∼2 (e.g., Ballesteros-Paredes et al. 2011;Camacho et al. 2016).
The second group of timescales consists of t shear , t orb , andt coll , all of which are ∼100 Myr. These timescales characterize dynamical processes taking place on kpc scales or even over entire galaxies. The order of magnitude contrast between them and the cloud "internal" timescales discussed above implies that the effects of galactic-scale dynamics on individual molecular clouds are likely modest. More specifically, t shear t cr ort ff indicates that shearing motions are generally small on 10 We note that some studies in the literature argue for a longer molecular cloud lifetime on the order of 100 Myr (e.g., Scoville & Hersh 1979;Koda et al. 2009). These studies typically use the timescales of galactic dynamical processes (such as orbital time or spiral arm crossing time) as anchoring points to derive molecular cloud lifetimes, which might partly explain why their estimated cloud lifetimes are comparatively longer. cloud scales relative to motions generated by turbulence or gravitational collapse (at least in most regions targeted by PHANGS-ALMA; also see Utreras et al. 2020); t orb t life means that molecular clouds can only last a small fraction of a complete orbital revolution around the galaxy center (see Chevance et al. 2020a); andt coll t life suggests that cloud-cloud collisions do  Figure 7. The x-axes of all panels are scaled to show the same logarithmic range. Measurements in the "disk" and "center/bar" samples are denoted by data points in blue and orange, respectively. We display Spearman's rank correlation coefficient for all measurements at the lower right corner in each panel, with black (gray) font color indicating a p-value smaller (larger) than 0.001.
The last group consists of t dep, mol by itself, which is about 1−3 Gyr across our sample. This range is consistent with many previous studies on the molecular dominated regions in nearby, massive, star-forming galaxies (e.g., Bigiel et al. 2008;Leroy et al. 2008;Utomo et al. 2017;Muraoka et al. 2019;Ellison et al. 2021). The large ratios between t dep, mol and all other characteristic timescales reaffirm the notion that star formation is inefficient: the implied star formation efficiency is 0.5−1% per free-fall time or turbulence crossing time (see Evans et al. 2014;Lee et al. 2016;Vutisalchavakul et al. 2016;Utomo et al. 2018, J. Sun et al. in preparation), ∼1% per cloud lifetime (Kruijssen et al. 2019b;Chevance et al. 2020aChevance et al. , 2021Kim et al. 2021a), and ∼10% per orbital revolution or cloud-cloud collision (Silk 1997;Kennicutt 1998b).
Beyond the typical values of the estimated timescales and their ratios, we also examine which timescales correlate the best with t dep, mol . Figure 8 shows that all the other timescales we considered show weak to no correlations with t dep, mol (judging from the small correlation coefficients). The only statistically significant trends are with t shear and t orb . They exhibit mild positive correlations with t dep, mol with coefficients of ρ = 0.18 and 0.14, quantitatively consistent with the results in Wong (2009). Since these two timescales only differ by a factor of π/(1 − β), and the measured 1 − β has a narrow dynamic range across our sample and a relatively large uncertainty, it is expected that t shear and t orb contain virtually the same amount of information and have similar predictive power for t dep, mol . 7. SUMMARY This work examines the fundamental correlations between molecular clouds and their host galaxy environments in 80 nearby, massive, star-forming galaxies targeted by the PHANGS-ALMA survey. It directly addresses one of the core science questions that motivated the PHANGS-ALMA survey: how do molecular cloud populations depend on local and global properties of the host galaxy? Taking advantage of the large, representative galaxy sample and the homogeneous, high-quality data provided by PHANGS-ALMA, we provide a first systematic description of the environmental dependence of the cloud populations residing in typical star-forming environments across the local universe.
To achieve this overarching goal, we use PHANGS-ALMA CO (2-1) imaging data products (Sun et al. , 2020aLeroy et al. 2021a,b) and CPROPS-based object catalogs (Rosolowsky et al. 2021, A. Hughes et al. in preparation) to determine a rich set of molecular gas properties on 60−150 pc scales. We further complement these molecular cloud scale measurements with  Note-The median value and 1σ range (i.e., 16−84 percentile range) of the characteristic timescales shown in Figure 7. Results in columns (2) to (4) are calculated by weighting each aperture equally, whereas those in columns (5) to (7) are calculated by weighting each aperture by the molecular gas mass it contains.
multiwavelength observations covering UV, optical, IR, and radio bands (e.g., Leroy et al. 2019;Querejeta et al. 2021, A. Razza et al. in preparation, A. Sardone et al. in preparation), as well as high-level data products including rotation curves (Lang et al. 2020) and morphological feature catalogs (Querejeta et al. 2021). Together, these ancillary data present the kpc-scale gas and stellar mass distributions, star formation activities, kinematic properties, and morphological structures in the host galaxy (see Figure 1 for a schematic summary). Following the cross-spatial-scale analysis framework developed by Sun et al. (2020b), we divide the sky footprint of each target galaxy into a series of hexagonal apertures, each 1.5 kpc in size (Figure 2). We aggregate cloud-scale molecular gas measurements within each aperture, and then calculate the mass-weighted, population-averaged properties. We also compile measurements of host galaxy properties as area-weighted averages across the ∼kpc-scale apertures. Our analysis covers 46,628 apertures in total, and 3,383 apertures with both cloud population measurements and host galaxy measurements. We publish these rich multiwavelength measurements online in machine-readable formats (see Appendix F).
Utilizing these databases, we construct basic statistical profiles for both the molecular cloud populations probed by the PHANGS-ALMA survey and the kpcscale sub-galactic environments they inhabit. We quantify empirical correlations between cloud population properties and host galaxy environmental metrics. We further perform a data-driven variable selection technique and identify a small subset of environmental metrics as primary predictors of the cloud population statistics. Our main findings are as follows: 1. The PHANGS-ALMA survey samples a wide range of host galaxy local properties. This is illustrated by the broad distributions of galactocentric radii, orbital kinematic properties, as well as kpc-scale gas, stellar, and star formation rate surface densities among all the kpc-scale averaging apertures ( Figure 3 and Table 1). Judged purely by aperture number counts, the most typical sub-galactic environment in our sample closely resembles the Solar Neighborhood. In comparison, most molecular gas mass is hosted in apertures with higher gas, stellar, and SFR surface densities, which are likely not matched by any kpc-scale region in the Milky Way.
2. Molecular cloud populations vary substantially across kpc-scale regions. This is seen in population-averaged cloud properties such as mass, surface density, velocity dispersion, and turbulent pressure from aperture to aperture ( Figure 4 and Table 2). These populationaveraged measurements have been corrected for the effect of galaxy inclination and finite data sensitivity with novel methods. We conclude that variations of cloud properties within and among galaxies are not merely random scatter from cloud to cloud, but reflect systematic change across sub-galactic environments.
3. Cloud population average properties appear significantly correlated with many local and global host galaxy properties ( Figure 5). The sense of these correlations indicate that cloud populations with higher average mass, surface density, and turbulence strength prefer galactic environments at smaller galactocentric radii, higher gas, star, and SFR surface densities, shorter orbital period, and stronger shear. Similar trends are also present with global galaxy properties.
4. Our BIC-based variable selection analysis yields a set of power-law predictive models that capture the most prominent trends for each cloud population-averaged property ( Table 3). The small number of indepen-dent variables appeared in these models suggests that most cloud-environment correlations can be reduced to the primary dependencies on a few local environmental metrics, especially on the kpc-scale molecular gas and SFR surface densities. The absence of global galaxy properties in these predictive models suggests that the correlations between molecular clouds and their local kpc-scale environment are more fundamental, and that galaxy-to-galaxy variations might arise merely as their consequences.
The rich multiwavelength measurements derived in this work have broad applications. We demonstrate one application scenario by deriving and comparing a set of characteristic timescales relevant to molecular cloud evolution and star formation (Figure 7 and Table 4). This further inquiry leads to the following findings: 5. The molecular cloud population average free-fall time and turbulent crossing time are around 5−20 Myr, comparable to typical molecular cloud lifetimes estimated in a subset of our target galaxies (Kruijssen et al. 2019b;Chevance et al. 2020aChevance et al. , 2021Kim et al. 2021a). These results support the notion that, when averaged across co-spatial populations, typical molecular clouds have virial parameters of order unity and only live for a few "internal" dynamical timescales.
6. The characteristic timescales of galactic-scale dynamical processes (including orbital revolution, shearing, and cloud-cloud collision) are around 100 Myr, or about an order of magnitude longer than the cloud "internal" timescales or their estimated lifetimes. This contrast seems to suggest that galactic dynamical processes would have to be highly efficient to have a pronounced impact on molecular clouds throughout their short lifetime.
7. The molecular gas depletion time ranges 1−3 Gyr across our sample, implying star formation efficiencies of 0.5−1% per cloud free-fall time or crossing time, ∼1% per cloud lifetime, and ∼10% per orbital time or cloud-cloud collision time.
8. Among all molecular cloud internal timescales and galactic dynamical timescales we considered, only orbital time and shearing time show statistically significant (yet weak) correlations with the molecular gas depletion time (Figure 8).
Our rich multiwavelength measurements have already supported multiple observational studies on PHANGS galaxies. These studies cover a broad range of topics, including the dynamical equilibrium of the ISM (Sun et al. 2020b), pressures in H ii regions (Barnes et al. 2021), morphological features in the stellar disks (Querejeta et al. 2021), nuclear gas outflows (Stuber et al. 2021), and the molecular gas-star formation cycle (Pan et al. 2022). We also expect similar applications in future studies on PHANGS targets examining molecular cloud lifecycle (J. Kim et al. in preparation), molecular cloud star formation efficiency per free-fall time (J. Sun et al. in preparation), and galaxy disk global instabilities (T. Williams et al. in preparation).
Beyond the projects mentioned above, we expect these databases to be useful for many purposes, and we highlight as few of them here. (1) Our cloud population measurements can be directly compared to similar measurements in other types of galaxies, such as dwarf galaxies (Mizuno et al. 2001;Leroy et al. 2006;Schruba et al. 2017;Imara & Faesi 2019), starburst galaxies (Ueda et al. 2012;Brunetti et al. 2021;Krieger et al. 2021), bulge-dominated early-type galaxies (Utomo et al. 2015;Espada et al. 2019;Liu et al. 2021), or even lensed, highz galaxies (Dessauges-Zavadsky et al. 2019). Such comparisons could highlight commonalities and differences among star-forming environments across the universe.
We plan to keep maintaining and improving the databases, thereby making them a long-term reference for the community. In particular, crucial next steps will come from incorporating measurements of ionised gas and stellar populations from the PHANGS-MUSE survey (Emsellem et al. 2021), as well as star clusters from the PHANGS-HST survey . Future versions of these databases will be released at the same location online as indicated in Appendix F.
This work was carried out as part of the PHANGS collaboration. The work of JS is partially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) through the Canadian Institute for Theoretical Astrophysics (CITA) National Fellowship. The work of JS, AKL, and DU is partially supported by the National Science Foundation (NSF) under Grants No. 1615105, 1615109, and 1653300. The work of JS and AKL is partially supported by the National Aeronautics and Space Administration (NASA) under ADAP grants NNX16AF48G and NNX17AF39G. EWR acknowledges the support of the NSERC, funding reference number RGPIN-2017-03987. EWK acknowledges support from the Smithsonian Institution as a Submillimeter Array Fellow. GAB gratefully acknowledges support by the ANID BASAL project FB210003. HAP acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 694343)  This work is based in part on observations made with NSF's Karl G. Jansky Very Large Array (VLA; project code: 14A-468, 14B-396, 16A-275, 17A-073, 18B-184). VLA is also operated by NRAO.
This work is based in part on observations made with the Australia Telescope Compact Array (ATCA). ATCA is part of the Australia Telescope National Facility, which is funded by the Australian Government for operation as a National Facility managed by CSIRO.
This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA.
This work makes use of data products from the Widefield Infrared Survey Explorer (WISE), which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by NASA.
This work is based in part on observations made with the Galaxy Evolution Explorer (GALEX). GALEX is a NASA Small Explorer, whose mission was developed in cooperation with the Centre National d'Etudes Spatiales (CNES) of France and the Korean Ministry of Science and Technology. GALEX is operated for NASA by the California Institute of Technology under NASA contract NAS5-98034.
This work is based in part on data gathered with the CIS 2.5m Irénée du Pont Telescope and the ESO/MPG 2.2m Telescope at Las Campanas Observatory, Chile.
This work has made use of the NASA/IPAC Infrared Science Archive (IRSA) and the NASA/IPAC Extragalactic Database (NED), which are funded by NASA and operated by the California Institute of Technology. Relevant datasets to this work include S4G team (2010) and Sandstrom (2019).

APPENDIX
A. GALAXY SAMPLE We list our galaxy sample in Table A1.      B. PRESCRIPTIONS FOR THE METALLICITY AND THE CO-TO-H 2 CONVERSION FACTOR In this work, we adopt empirical relation-based prescriptions to infer a local gas-phase metallicity (Section 3.3) and its associated CO-to-H 2 conversion factor (Sections 2.1.1 to 2.1.2) in each kpc-sized aperture. Here we detail these prescriptions and the rationale behind them.

B.1. Metallicity
To account for possible variations of CO-to-H 2 conversion factor across our galaxy sample, a key first step is to get homogeneous and reliable metallicity estimates. Although extensive compilations of (global and resolved) metallicity measurements for nearby galaxies exist in the literature (e.g., Pilyugin et al. 2014;De Vis et al. 2019), we do not yet have a uniform sample of resolved metallicity measurements with the same calibration scheme for all PHANGS-ALMA targets. In this work, we instead rely on two well-calibrated scaling relations to capture the general trends of metallicity variation across our sample.
Here Z (r e ) is the local gas-phase abundance at r gal = r e normalized by the solar value [12+log (O/H) = 8.69], Z (r gal ) is the normalized abundance at arbitrary r gal , and M is the galaxy global stellar mass derived by assuming a Chabrier IMF (Chabrier 2003). Note that these scaling relations are appropriate for abundance measurements adopting the O3N2 calibration in Pettini & Pagel (2004). While Equations B1 and B2 are identical to the formulae used in Sun et al. (2020b), we make two methodological improvements when applying them in this work. First, we elevate the M values in Table A1 by 0.1 dex before inserting them into Equation B1. According to Sánchez et al. (2019), this 0.1 dex offset can largely correct for systematic effects due to (a) differences between Salpeter and Chabrier IMFs, and (b) the finite aperture size of their IFU data (see their Appendix A and Figure A1). Second, we estimate r e by multiplying the stellar disk scale length, r disk , by a factor of 1.68. This step mirrors the procedure for deriving r e in Sánchez et al. (2019). Overall, these two methodological changes improve the self-consistency of our metallicity prescription.
As a sanity check, we compare our predictions to the observed two-dimensional metallicity distributions in 18 galaxies in the PHANGS-MUSE survey (Williams et al. 2021, also see Emsellem et al. 2021). Modulo the uncertain translation between different metallicity calibration schemes (i.e., O3N2 versus S-cal; see Pettini & Pagel 2004;Pilyugin & Grebel 2016), the predictions and the actual measurements show similar median values (within 0.05 dex) across this subsample, although the dynamic range of our predictions appears narrower than the observed range (0.12 dex versus 0.21 dex).

B.2. CO-to-H 2 Conversion Factor
The CO-to-H 2 conversion factor, α CO , is expected to depend strongly on metallicity (e.g., Wilson 1995;Arimoto et al. 1996;Israel 1997;Wolfire et al. 2010;Glover & Mac Low 2011;Feldmann et al. 2012;Narayanan et al. 2012;Schruba et al. 2012;Amorín et al. 2016;Accurso et al. 2017;Gong et al. 2020). Within a star-forming galaxy, α CO can vary by more than a factor of 2 (Leroy et al. 2011;Blanc et al. 2013;Sandstrom et al. 2013), with part of it attributable to metallicity variations (at least in the low temperature "outer disk" regime). These considerations motivate us to use a metallicity-dependent α CO prescription as a fiducial choice in this work. Nevertheless, we also consider a few other prescriptions and provide these alternative estimates in the published datasets.
Our fiducial estimate follows the same metallicity-dependent prescription as described in Sun et al. (2020b): where Z is the predicted local metallicity from Equations B1 and B2. The adopted power-law slope in Equation B3 is motivated primarily by the metallicity-dependent part of the xCOLD GASS calibration (Accurso et al. 2017), whereas the normalization is anchored to the Galactic value at solar metallicity (including the gas mass contribution by helium; see Bolatto et al. 2013). This prescription gives similar predictions to many other prescriptions in the literature (e.g., Genzel et al. 2012;Schruba et al. 2012;Amorín et al. 2016) within the metallicity range probed in this work (e.g., see Figure 6 in Accurso et al. 2017).
Beyond this fiducial α CO prescription, we calculate four alternative prescriptions, following and expanding on Sun et al. (2020b). The first is simply a constant value matching the Galactic average: The second prescription follows Narayanan et al. (2012) and infers α CO from both metallicity (Z ) and the fluxweighted CO (2-1) line intensity ( I CO(2−1) ): The above formula is adapted from equation (11) in Narayanan et al. (2012) with two notable distinctions. First, the original formula depends on I CO(1−0) , whereas our formula converts that dependence to I CO(2−1) assuming a line ratio of R 21 = 0.65 (den Brok et al. 2021;Leroy et al. 2021c). Second, we increase the normalization by a factor of 1.36 to correct for helium contribution. The third prescription follows Bolatto et al. (2013) and infers α CO from metallicity (Z ), molecular cloud surface density (proxied by Σ mol, pix ), and kpc-scale total surface density including both gas and stellar mass (Σ total ): Since calculating Σ mol, pix and Σ total relies on knowing α CO in the first place, we solve for α CO iteratively until the output of Equation B6 converges to the assumed value for calculating Σ mol, pix and Σ total (also see Equations 24-26 in Sun et al. 2020b). We combine the above prescriptions for α CO(1-0) with the adopted CO (2-1)-to-CO (1-0) line ratio R 21 = 0.65 (den Brok et al. 2021;Leroy et al. 2021c) to get the appropriate conversion factor for the CO (2-1) line, α CO(2-1) . We then apply these α CO(2-1) values on a per aperture basis -that is, we assume a constant conversion factor within the ∼kpc scale extent of each averaging aperture. These treatments largely follow Sun et al. (2020b).
As another improvement over Sun et al. (2020b), here we add a fourth alternative prescription following Gong et al. (2020). This simulation-motivated prescription considers the dependence on metallicity (Z ), CO line integrated intensity (I CO ), and the physical beam size (D beam ). It directly predicts the conversion factor for the CO (2-1) line without relying on a separately assumed R 21 value. The original formula is expressed in number column density convention (i.e., Equation 4b in Table 3 We convert it into mass surface density units via α CO(2-1) M pc −2 (K km s −1 ) −1 = 2.18 X CO(2-1) 10 20 cm −2 (K km s −1 ) −1 . (B8) Since this prescription is calibrated at 100 pc scales, we derive the α CO(2-1) values pixel-by-pixel at 60-150 pc scales and then calculate CO line intensity weighted mean values within the ∼kpc-scale averaging apertures. This is done for all four resolution levels considered in this work (60, 90, 120, and 150 pc). We include our estimated α CO(2-1) values from all aforementioned prescriptions in the published databases. This allows for easy conversions if the reader wish to adopt an alternative prescription instead of our fiducial choice. Concretely, our molecular cloud measurements scale with the adopted conversion factor as M ∝ Σ ∝ P turb ∝ α CO(2-1) and α vir ∝ α −1 CO(2-1) ; the kpc-scale molecular gas surface density scales as Σ mol ∝ α CO(2-1) ; the timescale measurements scale ast ff ∝ α −0.5 CO(2-1) and t dep, mol ∝ α CO(2-1) . The impact of different prescriptions on some of the molecular gas measurements is examined in detailed by Sun et al. (2020b).
We expect future works to improve the handling of the CO-to-H 2 conversion factor even further. Particularly, combining a varying R 21 (either observed directly or predicted based on similar observations; Leroy et al. 2021c) with an R 21 -dependent α CO prescription (Gong et al. 2020) would allow us to better capture the gas excitation temperature variations, especially in galaxy centers where this effect becomes very pronounced. We also defer a more thorough comparisons between the different α CO prescriptions to a subsequent paper.
C. INCLINATION CORRECTIONS In Section 2.1.1 and 2.1.2, we introduce inclination corrections on several molecular cloud properties measured at 60-150 pc scales. These corrections represent an important methodological change, as many previous works (including Sun et al. 2018Sun et al. , 2020aRosolowsky et al. 2021, using the same data set) have assumed spherical geometry for observations on these spatial scales and, consequently, have not applied such corrections. This spherical approximation assumes that at the resolution scale, the structure of the molecular ISM is isotropic. While this assumption is common throughout the literature, it is not well tested. Furthermore, the 60-150 pc resolutions are also comparable to the thickness of the molecular gas disk, which implies anisotropy. In this appendix, we show that there is an inclination dependence in our measurements and motivate specific functional forms for empirical corrections.
The primary motivation for introducing these corrections is that the observed molecular gas properties (such as surface density and velocity dispersion) at ∼ 100 pc scales show apparent correlations with the host galaxy inclination angle. For example, A. Hughes et al. (in preparation) found such trends in the measured properties of CPROPS-identified objects among the PHANGS-ALMA galaxy sample. A parallel examination of the pixel-by-pixel measurements derived from the same dataset (Sun et al. 2020a) also reveals a similar trend (see Figure C1 left panel). Since inclination angle is not an intrinsic property of galaxies, the presence of these inclination-related trends signifies systematic biases in the molecular gas measurements derived with both CPROPS-based and pixel-by-pixel approaches.
To quantify these biases, we carry out a modified version of the variable selection analysis in Section 5.2. We suppress the inclination corrections for all molecular cloud properties (i.e., the cos i terms in Equations 1-13) while keeping cos i in the list of feature variables. This re-analysis yields a new set of power-law predictive models similar to those in Table 3, but many of these models now carry an extra (cos i) β term. Specifically, in the predictive models for Σ obj and Σ pix , the power of the cos i terms are close to β = −1.0; while in the models for σ obj and σ pix they are close to β = −0.5. In other words, to eliminate the apparent cos i dependence in these models (which should not be present were these models reflecting purely physical trends), one would have to multiply cos i to the measured surface densities and (cos i) 0.5 to the velocity dispersions. This motivates the correction terms in Equations 3, 4, 8, and 10.
We also inspect the effects of these inclination corrections on the original pixel-by-pixel measurements in Sun et al. (2020a) without doing any cloud population averaging. Figure C1 compares the surface density-velocity dispersion relation before and after applying the inclination corrections. Without these corrections, the median velocity dispersion at given surface densities appears tends to be higher in galaxies with high inclination angles (smaller cos i). This trend largely disappears when we apply the inclination corrections to both axes. This result suggests that, despite being motivated by the observed cloud-environment correlations, our adopted inclination corrections can also remove the unphysical inclination dependence in the relationships among molecular cloud internal properties.
Our adopted inclination corrections are empirical and data-driven, but their functional forms have physical implications. For the surface densities, a multiplicative term of cos i is the exact correction one would use for disk-like structures with their orientations aligned with the whole galaxy disk. Our proposed interpretation is that the interstellar gas forms filamentary networks, which preferentially align with the large-scale galaxy disk even at ∼ 100 pc scales (Zucker et al. 2018). This preference in orientation should eventually disappear at smaller scales, but evidently, the PHANGS-ALMA observations do not yet reach the transitional spatial scale. Figure C1. The relationship between molecular gas surface density (Σ mol ) and velocity dispersion (σ mol ) measured pixel-bypixel at 150 pc scales, without and with inclination corrections (left and right, respectively). Both plots are made using the published data tables in Sun et al. (2020a). The colors of the data points represent the cos i values of the host galaxies (see colorbar), but have been median-filtered to bring out the overall trend across the parameter space. The four colored lines in each panel show the running median of σ mol at fixed Σ mol for galaxies in four different inclination bins. The left panel (without inclination corrections) reveals a mild but statistically significant trend of elevated σ mol at fixed Σ mol in galaxies with higher inclination (i.e., smaller cos i). This trend largely disappears in the right panel (with the inclination corrections applied).
For the velocity dispersions, there are at least two effects that can produce some inclination dependence: (1) contributions from ordered, in-plane motions of the gas (e.g., beam smearing), and (2) anisotropy of the gas velocity dispersion (usually with the in-plane components larger than the vertical component; see e.g., Jeffreson et al. submitted). Both effects would predict higher velocity dispersions at higher inclination angles, which is consistent with the direction of the adopted correction factor, but neither would call for a specific functional form of (cos i) 0.5 for this correction. Alternatively, if one assumes that the velocity structure of the turbulent gas can be described by a line width-size relation of σ mol (l) ∝ l 0.5 (e.g., Solomon et al. 1987), and that this relationship still holds at 60-150 pc scales, then the varying line-of-sight depth with inclination (i.e., l ∝ 1/ cos i) could imply σ mol ∝ (cos i) −0.5 , which matches our empirical result.
We stress that all the inclination-dependent trends we identify above are real measurements that only emerge in statistical analysis of many galaxies at similarly high resolution. They imply that the molecular gas structures are clearly anisotropic at 60-150 pc scales. Existing and future surveys at even higher resolution (e.g., in very nearby galaxies or in CO-bright sub-regions of PHANGS galaxies) can help us extend this analysis to smaller spatial scales, where the transition from anisotropic to isotropic structures presumably occurs. This transitional spatial scale can be further compared with estimates of the gas disk scale height from independent methods (e.g., Koch et al. 2020, Jeffreson et al. submitted).
D. COMPLETENESS CORRECTIONS In Section 3.2.2, we identify a systematic bias affecting the population-averaged molecular cloud properties. This systematic bias originates from the incomplete CO flux recovery in the PHANGS-ALMA "strict" moment maps and the associated CPROPS catalogs. We introduce a completeness correction to account for this bias, making use of the measured CO flux completeness, f flux , and area coverage fraction, f area , in each averaging aperture. In this appendix, we present the mathematical derivation of this completeness correction method.
We assume that the intrinsic CO intensity probability distribution function (PDF) within each averaging aperture follows a lognormal distribution: This assumption is motivated by observational constraints on the CO intensity distribution in nearby galaxies (e.g., Hughes et al. 2013b;Leroy et al. 2016;Egusa et al. 2018;Sun et al. 2018Sun et al. , 2020a. We further assume that the CO emission included in the "strict" moment maps constitutes everything above a threshold intensity I CO, th (i.e., the maps have a sharp, well-defined sensitivity limit), and that the "strict" maps capture all emission along the sightlines with CO detections. Under these assumptions, we can express f flux and f area in terms of σ int ,Ī CO, int , and I CO, th : Notice that these two relations allow us to inversely solve for σ int and ln(I CO, th /Ī CO, int ) and express their combinations in terms of the measurable quantities f area and f flux : Here "erf −1 " stands for the inverse error function. We then calculate the appropriate correction factors for the population-averaged molecular cloud properties and express them in terms of f area and f flux . The three population-averaged properties we consider here are the the fluxweighted average cloud surface density, Σ (Section 2.1.1-2.1.2), the molecular gas clumping factor, c pix (Section 3.2.1), and the flux-weighted average of the reciprocal of free-fall time, t −1 ff (Section 6.1).
where f compl (β) = +∞ I CO, th I β CO P (I CO ) dI CO +∞ 0 I β CO P (I CO ) dI CO = 1 2 1 − erf ln(I CO, th /Ī CO, int ) − βσ 2 int √ 2σ int Note that the second steps in Equations D6-D8 is valid because we adopt a constant α CO within each aperture. Figure D1 shows the joint distribution of f area and f flux for all 3,383 hexagonal apertures, with color-codes reflecting the amplitude of the derived correction factors F correct, Σ and F correct, c . The f area and f flux values for ∼68% of apertures are consistent with intrinsic lognormal I CO -PDFs with σ int ranging from 0 to 1.0 dex. This broadly matches the observed I CO -PDF widths of 0.2 to 0.6 dex in sub-regions of nearby galaxies (e.g., Hughes et al. 2013b;Sun et al. 2018Sun et al. , 2020a The implied I CO PDF widths for another ∼12% of apertures exceed 1.0 dex, which seem too high to be physical. These apertures tend to have low f area , possibly suggesting that these apertures include true "empty" areas devoid of CO emission, in which case the lognormal PDF assumption is no longer appropriate. Finally, the remaining ∼20% of apertures have f flux < f area . Given the high f area in most of these regions, it is likely that there is missing CO flux along sight lines with CO detections. We calculate the completeness correction for these apertures by assuming an ad-hoc f flux value equal to f area .
The color trends in Figure D1 indicate that apertures with low f flux or f area would require very significant completeness corrections, which means that the CO detections in these apertures are too "unrepresentative" of the underlying cloud population for us to extract reliable statistics. This motivates us to select a subsample of apertures with high Figure D1.
Area coverage fraction, farea, versus CO flux recovery fraction, f flux , on 150 pc scales within each of the 3,383 hexagonal apertures. Black lines show the expected farea vs. f flux relations for sensitivity-limited CO observations given lognormal-shaped intrinsic ICO-PDFs (different line styles correspond to lognormal PDFs with 1σ widths ranging from 0 to 1.0 dex). The colors of the data points reflect the amplitude of the appropriate correction factors on Σ150pc (left) and cpix, 150pc (right) according to Equations D6 and D7. A red box highlights the parameter space with f flux > 0.5 and farea > 0.2, which are the criteria for selecting the subsample of 871 apertures with high fidelity measurements in Section 5.
f flux and high f area for the more careful correlation analyses in Section 5. The selection criteria, f flux > 50% and f area > 20%, are also illustrated in Figure D1. Among all apertures that meet these criteria, the correction factors on Σ have a median value of 0.03 dex and a maximum of 0.2 dex, which means that the uncorrected Σ values are already close to the inferred true population average. For the same set of apertures, the correction factors on c pix have a smaller median of 0.01 dex but a much larger maximum of 0.6 dex. The few apertures with very large correction factors are those with high f flux but low f area , where a simple lognormal PDF is likely not a good description of the underlying CO intensity distribution.
E. RESOLUTION DEPENDENCE In this work, we derive molecular cloud measurements from the PHANGS-ALMA CO data at four common resolution levels: 60, 90, 120, and 150 pc. The main text focuses on results derived at 150 pc scale so as to cover the full sample of apertures while keeping the presentation succinct. In this appendix, we draw comparisons across all four spatial scales to illustrate the scale-dependence of the measured molecular cloud properties.
The molecular cloud measurements are available for fewer and fewer apertures/galaxies as the resolution goes from 150 pc to 60 pc (see Leroy et al. 2021a). To control for changes in the aperture/galaxy sample and any associated selection effects, here we focus on a subset of 328 apertures in 15 galaxies with data at all four spatial scales. Besides, the data sensitivity also drops as the beam size decreases, which could also leads to systematic biases in our molecular cloud population statistics. We address this issue by applying completeness corrections to all population-averaged measurements according to the flux completeness and area coverage fraction of the CO moment maps at each resolution (see Section 3.2.2). Figure E1 illustrates that most molecular cloud properties presented in Figure 4 show some level of resolution dependence. In detail, the average molecular cloud mass and radius both increase strongly as the data resolution gets coarser (the former is most likely driven by the latter). This reinforces the conclusion that the sizes of the CPROPS-identified objects are primarily set by data resolution rather than physical properties of the gas distribution (see Section 5.2). The slope of the R ojb trend appears sub-linear ( R obj ∝ D 0.6 beam ), which seems to suggest that the molecular gas structure is not completely scale-free between 60-150 pc. The average cloud surface density mildly decreases with beam size, as expected from more beam dilution. The average cloud velocity dispersion increases with beam size, but with a power-law slope of ∼0.20-0.23, shallower than the line width-size relation for molecular clouds in the Milky Way disk (e.g., Solomon et al. 1987). The average cloud turbulent pressure mildly declines toward coarser resolution, whereas the average virial parameter largely remains roughly constant. Both trends are predictable from Figure E1. Resolution dependence of the molecular cloud population-averaged properties shown in Figure 4. In each panel, the four "violins" represent histograms of the corresponding cloud properties at 60, 90, 120, and 150 pc resolutions. These histograms are made from a common set of apertures for which measurements at all four resolutions are available. The black open square and vertical bar indicate the median value and 16-84 percentile range (assigning equal weight per aperture). their functional relations with cloud mass, size, and velocity dispersion. Finally, as the molecular gas surface density distribution becomes more homogeneous at coarser resolutions, the estimated clumping factor diminishes accordingly.
F. MACHINE-READABLE TABLES We publish our high-level measurements in the form of machine-readable data tables via the PHANGS CANFAR storage 13 . A permanent copy of the version used in this article will also be available via the CANFAR Data Publication Service 14 . This version includes two types of tables. The first records the hexagonal aperture measurements, and the second presents the radial bin measurements (see Section 3.1). Here we provide column-by-column descriptions of the hexagonal aperture tables in Table F1. For clarity, we also add links to the relevant sections and equations for each derived quantity. We do not separately provide column-by-column descriptions for the radial bin tables, but simply note that the set of columns therein are almost identical, except that the radial bin tables lack the RA, DEC, and phi_gal columns.
We plan to keep improving these data tables and release subsequent versions via the same CANFAR storage 13 . Future versions will cover a larger sample of galaxies and include more measurements derived from other data sets, such as H ii regions and stellar populations from the PHANGS-MUSE survey (Emsellem et al. 2021) and star clusters from the PHANGS-HST survey ).  Statistical error on logarithmic derivative of CO rotation curve Zprime Gas-phase metallicity relative to solar (Eq. B1-B2) alpha_CO21_S20 s M K −1 km −1 pc −2 CO (2-1)-to-H2 conversion factor (fiducial; Sun+20; Eq. B3) alpha_CO21_N12 s M K −1 km −1 pc −2 CO (2-1)-to-H2 conversion factor (Narayanan+12; Eq. B5) Statistical error on flux-weighted mean reciprocal of free-fall time @ X * pc